bashで遊ぶ9

はじめに

Project Eulerを通してbashを写経しながら
実行結果を眺めたりググったり
mawkとgawkの違いに戸惑いつつ
Colaboratoryを主要な実行環境とする

Colaboratoryとは

Colaboratoryについて(まとめ)

ついでにPythonでも遊ぶ
(使えるものは使って効率重視,合理性を目的にしたい)
(実装効率、時間効率、空間効率の順に重視)
(目標はミリセコンドオーダー以下;ターゲットであってコミットメントでは無い)
(予想される処理時間より効率が良ければそれでOKとする)
(ColaboratoryではCython, Numbaを使うのが難しいのでそういうのは無し)
(ただCython欲が出た時にはオンボロPCで別途試してみる事もある)
(Project Eulerを解く事自体は目的ではない)

参考にさせて貰っているもの

tea63 – Qiita

さんのBashで「Project Euler」を解く記事

これまで

Project Euler 1 ー Project Euler 4

bashで遊ぶ

Project Euler 5 ー Project Euler 11

bashで遊ぶ2

Project Euler 12 ー Project Euler 18

bashで遊ぶ3

Project Euler 19 ー Project Euler 25

bashで遊ぶ4

Project Euler 26 ー Project Euler 30

bashで遊ぶ5

Project Euler 31 ー Project Euler 36

bashで遊ぶ6

Project Euler 37 ー Project Euler 43

bashで遊ぶ7

Project Euler 44 ー Project Euler 49

bashで遊ぶ8

 

今回

Project Euler 50 ー Project Euler 53

 
 
Project Euler Q50 【連続する素数の和】 – Qiita

mawk(Colaboratoryに入っているやつ)では

%%bash
time {
seq 1000000 |
factor |
awk 'NF==2{printf $2" "}' |
awk '{for(i=1;i<=NF-1;i++){s=0;c=0;for(j=i+1;j<=NF;j++){if(1000000<=s+$j){break}else{s+=$j;c+=1}};print s,c}}' |
awk '21<$2' |
sort -k2,2nr |
awk '$0=$1' |
factor |
awk 'NF==2{print $2}' |
head -1
}

awk: program limit exceeded: maximum number of fields size=32767
FILENAME=”-” FNR=1 NR=1

real 0m0.322s
user 0m0.414s
sys 0m0.131s

gawkで試す。

$ ▶ time {
> seq 1000000 |
> factor |
> awk 'NF==2{printf $2" "}' |
> awk '{for(i=1;i<=NF-1;i++){s=0;c=0;for(j=i+1;j<=NF;j++){if(1000000 awk '21 sort -k2,2nr |
> awk '$0=$1' |
> factor |
> awk 'NF==2{print $2}' |
> head -1
> }
997651

real    0m2.500s
user    0m2.848s
sys 0m0.128s

Pythonの場合:

100万未満の素数列を出す
前1,000個位の累積和を出す
累積和の後ろから後ろー前して素数になったらそこが最大

だろう。前1,000個でいいかどうかはまずcumsumして確かめる。

import numpy as np
from sympy import sieve

sieve.extend(1000000)

p_arr = np.array(sieve._list)
p_arr[:550].cumsum()[-10:]

array([ 978042, 981959, 985878, 989801, 993730, 997661, 1001604,
1005551, 1009518, 1013507])

1,000個どころか547個目で100万を超える様だ。

なので546個の累積和を出しておいて
累積和[::-1] – 累積和のブーリアンを求めれば終わり

import numpy as np
from sympy import sieve


def p50():
  sieve.extend(1000000)
  p_arr = np.array(sieve._list)
  cum_arr = p_arr[:546].cumsum()
  ind1, ind2 = np.where(np.isin((cum_arr[-200:] - cum_arr[:200][:, np.newaxis]), p_arr[-200:]))
  return cum_arr[-200:][ind2.max()] - cum_arr[:200][ind1.min()]


p50()

997651

%timeit p50()

100 loops, best of 3: 4.7 ms per loop

Numpyで解ける問題は本当に楽

 
 

Project Euler Q51 【素数の桁置換】 – Qiita

gensubはgawkの拡張機能っぽい

 $ ▶ time {
> seq 100000 999999 |
> factor |
> awk 'NF==2{print $2}' |
> awk -v FS= '{delete s;printf $0" ";for(i=1;i awk '{print $1,$3*1,$4*1,$5*1}' |
> awk '{for(i=2;i<=NF;i++){if($i<2){continue};printf $1" ";for(j=0;j factor |
> tr ' ' '_' |
> xargs -n11 |
> tr '_' ' ' |
> awk '{for(i=1;i<=NF-1;i++){if(0 awk 'NF==9&&length==(6*9+9){print $1}'
> }
121313

real    0m38.546s
user    0m6.016s

Pythonの場合:

どう考えても時間が掛かりそうでやる前から萎える……

・素数は3の倍数ではない(p mod 3 = 1 or 2)
置換桁の合計mod 3 = 0 → 合計に影響しない
置換桁の合計mod 3 = 1 or 2 → いずれかが3の倍数になる

従って置換桁が3の倍数以外では最大7通りしか得られないので、
置換桁は3の倍数(3桁か6桁)だと考えられる。

・1の位は置換しない
10以上の素数の1の位は常に1,3,7,9のいずれか。
8通りはあり得ないので除外。

・置換桁が3桁で1の位は置換できないので素数は5桁か6桁

でも5桁が例示されているって事は5桁じゃなさそうだけどね

・8通りなので0, 1, 2のいずれかで置換できるパターンがある

上記の条件から実装していくだけ。



……でまあ実装すべく色々考えてみたんだけど、
正直この問題はもう考えたくない……疲れた。

前問5分位で書けて「楽勝じゃん」と思った後のこれ。

しかも、結局は効率の良いコードを思い浮かばなくて、
ネットで調べて人様のコードを参考にさせて頂いた。

Project Euler – Problem 51 – 漫ろ草 SozoroGusa

元のコードがそもそも凄く効率的で、素数生成部だけsympy.sieveで置き換えたけど、
Python2で実行して500 ms、Python3(xrange→rangeに)で実行して790 msだった。
(自分でどんだけ頭ひねっても1.2秒の後のこれだったので凄く感動した\(^o^)/)

んで、そこから最適化してみたら、

from collections import defaultdict
from functools import lru_cache
from sympy import sieve


def rplnum(num, pos):
  for i in range(10):
    res = ''
    for j, k in enumerate(str(num)):
      res += str(i) if j in pos else k
    yield int(res)


@lru_cache(maxsize=None)
def getmember(p, d_rp):
  dd = defaultdict(list)
  for v, k in enumerate(str(p)):
    dd[k].append(v)
  for k, v in sorted(dd.items()):
    if len(v) == d_rp:
      yield int(k), v


def p51():
  d_rp, dt = 3, (5, 6)
  for i in dt:
    ps = set(sieve.primerange(10**(i-1), 10**i))
    for p in ps:
      for val, pos in getmember(p, d_rp):
        if val > 2 or pos[-1] >= i - 1:
          break
        cnt = 0
        for u in rplnum(p, pos):
          if u not in ps:
            cnt += 1
          if cnt > 2:
            break
        if cnt == 2:
          return p


p51()

121313

%timeit p51()

10 loops, best of 3: 38.9 ms per loop

素数列はリストよりも集合の方がだいぶ速くなる。
(変換処理とin演算子での処理時間が大きく変わってくる)
このケースでは”.join()よりも+オペレータで回す方が多少速い。
defaultdictは呼び出しコストが大きいので、
毎回呼び出す場合はメモ化関数では「dict.get(k, []) + [v]」の方が速いけど、
メモ化関数ではdefaultdictを使う方が速い。
細々した違いだけでも約500 ms位変わった(約790 → 330 ms)けど、
とにかくメモ化で10倍位変わる。
文字列と整数の変換は一見高コストにみえるけど、
モジュロも高コストなので、計算コストに大差なく、
それだったら文字列を用いる方が実装コストが低い。

それにしても、最後放り投げたとは云え、前問との落差で心折られた。
ここまで頑張っても、多分あれだろう。素数リスト作った後、愚直にループ回して
Cythonでがっつり最適化したらそっちの方が速いんだろうな……(´・ω・`)
(ColaboratoryではCython/Numbaが使えな(使えるようにできない))

 
 
 
 
 

Project Euler Q52 【置換倍数】 – Qiita

mawkで実行できない問題が増えてきた。。

Pythonの場合:

心が折れた後の息抜き問題。

sum(x) = 3*sum(x)なのでxは3の倍数だと考えられる
x, 2xで同じ数を含むという事は3の倍数でかつ2倍して桁上がりする様な数
桁上がりする事を考慮すると9の倍数(と云えるか確信はないが)
それぞれ1桁目が違う値を持つのは6桁以上
特に6桁だけを考えると探索範囲は[100000, 166667)
5の倍数があるから0か5が含まれる

def p52():
  for i in range(100008, 170000, 9):
    s = set(str(i))
    if( ('0' in s or '5' in s) and
        s == set(str(2*i)) and
        s == set(str(3*i)) and
        s == set(str(4*i)) and
        s == set(str(5*i)) and
        s == set(str(6*i))
      ):
      return i
  return 'not found'


p52()

142857

%timeit p52()

100 loops, best of 3: 5.21 ms per loop

ifの中の条件の順番を入れ替えるだけでちょっとだけ速くなるけど。

 
 
Project Euler Q53 【組み合わせ選択】 – Qiita

%%bash
time seq 23 100 |
awk '{for(i=1;i<=$1;i++){print $1,i,$1-i}}' |
awk '{s=1;t=1;u=1;for(i=1;i<=$1;i++){s*=i};for(i=1;i<=$2;i++){t*=i};for(i=1;i<=$3;i++){u*=i};print s/t/u}' |
awk '1000000<$1' |
wc -l

4075

real 0m0.073s
user 0m0.070s
sys 0m0.001s

Pythonの場合:

Numpy/SciPyの頼もしさ
100 x 100通りなのでブロードキャストでガサッと処理できる

import numpy as np
from scipy.special import comb

def p53():
  n = np.arange(1, 101)
  return (comb(n[:, np.newaxis], n) > 1000000).sum()

p53()

4075

%timeit p53()

1000 loops, best of 3: 826 µs per loop

カテゴリー: 未分類 パーマリンク

bashで遊ぶ9 への1件のフィードバック

  1. ピンバック: bashで遊ぶ10 | 粉末@それは風のように (日記)

コメントは受け付けていません。