bashで遊ぶ6

はじめに

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

Colaboratoryとは

Colaboratoryについて(まとめ)

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

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

tea63 – Qiita

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

これまで

Project Euler Q1 ー Project Euler Q4

bashで遊ぶ

Project Euler Q5 ー Project Euler Q11

bashで遊ぶ2

Project Euler Q12 ー Project Euler Q18

bashで遊ぶ3

Project Euler Q19 ー Project Euler Q25

bashで遊ぶ4

Project Euler Q26 ー Project Euler Q30

bashで遊ぶ5

今回

Project Euler Q31 ー Project Euler Q36

 
 
 

Project Euler 31 硬貨の和

Project Euler Q31 【硬貨の和】 – Qiita

これ凄いよなあ.

%%bash
time {
seq 0 200 |
awk '{print "1 "$1}' |
awk '{for(i=0;i<=200/2;i++){v=$1*$2+2*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/5;i++){v=$1+5*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/10;i++){v=$1+10*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/20;i++){v=$1+20*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/50;i++){v=$1+50*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/100;i++){v=$1+100*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/200;i++){v=$1+200*i;if(v==200){print v}}}' |
wc -l
}

73682

real 0m4.990s
user 0m9.651s
sys 0m0.107s

%%bash
time {
seq 0 200 |
awk '{for(i=0;i<=200/2;i++){v=$1+2*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/5;i++){v=$1+5*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/10;i++){v=$1+10*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/20;i++){v=$1+20*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/50;i++){v=$1+50*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/100;i++){v=$1+100*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/200;i++){v=$1+200*i;if(v==200){print v}}}' |
awk 'END{print NR}'
}

73682

real 0m5.193s
user 0m10.118s
sys 0m0.091s

Pythonで解いてから見直していて気付いたけど,
ループ構造を逆にするだけで速くなる筈.
また,2ポンドの場合は,1通りしか無いので省略すると(結果+1にすると),

%%bash
time {
seq 0 200 |
awk '{for(i=0;i<=200/100;i++){v=$1+100*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/50;i++){v=$1+50*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/20;i++){v=$1+20*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/10;i++){v=$1+10*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/5;i++){v=$1+5*i;if(v<=200){print v}}}' |
awk '{for(i=0;i<=200/2;i++){v=$1+2*i;if(v==200){print v}}}' |
awk 'END{print NR+1}'
}

73682

real 0m2.045s
user 0m2.151s
sys 0m0.010s

 

Pythonの場合:

これは部分和問題でNP困難かつNP完全な問題.
解く方法としてはDPな話.
真面目に解く(一般化して考える)場合,メモ化DFSとか
pulpとかで解ける形に定式化して組合せ最適化を考えるだろうか.

まずはNumpyでクルードに解いてみたい.
(大体10億万通りを超えると使い物にならないけど)
(今回のケースだとそのままやると10億通り超える)
(何気に1億通り以下ならそこそこ使える)
(1000万通り以下ならmsオーダー)

2ポンドの場合は1通りしか無いので計算は不要.
1ペンス以外の合計が2ポンド(200ペンス)以内の場合,
そこに1を足していけば必ず2ポンド(200ペンス)が作れる.
つまり,1ペンス以外の合計が2ポンド以下になる組み合わせだけ求めれば良い.

import numpy as np

def p31_numpy():
  a = np.arange(0, 200, 2, dtype=np.uint16)[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis, np.newaxis]
  b = np.arange(0, 200, 5, dtype=np.uint16)[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis]
  c = np.arange(0, 200, 10, dtype=np.uint16)[:, np.newaxis, np.newaxis, np.newaxis]
  d = np.arange(0, 200, 20, dtype=np.uint16)[:, np.newaxis, np.newaxis]
  e = np.arange(0, 200, 50, dtype=np.uint16)[:, np.newaxis]
  f = np.arange(0, 200, 100, dtype=np.uint16)
  return (a + b + c + d + e + f <= 200).sum() + 7

p31_numpy()

73682

%timeit p31_numpy()

10 loops, best of 3: 51.6 ms per loop

こんなに雑な解き方でも割と高速に求まる.
たかだか数千万オーダーなら簡便な方法として結構お薦め.

ちなみに,1まで含めて求めると10億オーダーになるので,

1 loop, best of 3: 13.4 s per loop

という感じになる.一応,解けるけど,非常に遅い(後メモリが流石にヤバい).

いくらなんでもこれで終わらせるのは宜しくないので,DFSで考える.

def count(lst, t):
  if len(lst) == 0 or t == 0:
    return 1
  else:
    res, c = 0, lst[0]
    for i in range(t//c + 1):
      res += count(lst[1:], t-c*i)
    return res


def p31():
  clst = [100, 50, 20, 10, 5, 2]
  target = 200
  return count(clst, target) + 1


p31()

73682

%timeit p31()

10 loops, best of 3: 23.4 ms per loop

与えるリストの順番を逆にすると処理時間が8倍位になる

同じ処理を繰り返しているのが分かるのでメモ化したいが,
「functools.lrucache」は引数にリストとかあると使えない.
(辞書のキーにリストを指定できないから)

整数分割数の問題としてみた場合,もっとスマートに解けそう.

 
 

Project Euler 32 パンデジタル積

Project Euler Q32 【パンデジタル積】 – Qiita

%%bash
time {
seq 9876 |
awk -v FS= '{s=1;for(i=1;i<=NF;i++){s*=gsub($i,$i)};if(s==1){print}}' |
factor |
awk 'NF!=2{print $1*1}' |
awk '{for(i=1;i<=$1;i++){if($1%i==0){print i"@"$1/i,$1}}}' |
grep -v "0" |
awk -v FS= 'length==11{s=1;for(i=1;i<=NF;i++){s*=gsub($i,$i)};if(s==1){print}}' |
uniq -f1 |
awk '{s+=$2}END{print s}'
}

45228

real 0m3.862s
user 0m3.972s
sys 0m0.024s

Pythonの場合:

適当な数字を掛けてざっくり考えてみると,1-9の順列を出した後,
考えるのは1×4桁か2×3桁位で良くて,探索範囲は4,000-9,000位.
かなり強い仮定を置くけど,後は普通に計算すれば良い.
見た目は汚くても桁の計算はこれ位なら書いた方が速い.

import itertools


def p32():
  res, mem = 0, {}
  for x in itertools.permutations(range(1, 10)):
    if 3 < x[0] < 9:
      a = 1000*x[0] + 100*x[1] + 10*x[2] + x[3]
      b = 10*x[4] + x[5]
      c = 100*x[6] + 10*x[7] + x[8]
      d = x[4]
      e = 1000*x[5] + 100*x[6] + 10*x[7] + x[8]
      if a == b*c:
        if a not in mem:
          res += a
          mem[a] = (b, c)
      elif a == d*e:
        if a not in mem:
          res += a
          mem[a] = (d, e)
  return res


p32()

45228

%timeit p32()

10 loops, best of 3: 195 ms per loop

 
 

Project Euler 33 桁消去分数

Project Euler Q33 【桁消去分数】 – Qiita

%%bash
time {
seq 12 98 |
awk '{for(i=10;i<$1;i++){if($1%11!=0&&i%11!=0){print i,$1}}}' |
awk '{for(i=1;i<=10;i++){printf $0" "};print ""}' |
awk '{for(i=1;i<=9;i++){gsub(i,"",$(2*i+1));gsub(i,"",$(2*(i+1)));print $1,$2,$(2*i+1),$(2*(i+1))}}' |
awk '$1$2!=$3$4&&$4!=0&&int($3/$4)==0' |
awk '$1==$3*$2/$4' |
awk 'BEGIN{a=1;b=1}{a*=$3;b*=$4}END{print a,b}' |
awk '{a=$1;b=$2;for(i=$1;2<=i;i--){if(a%i==0&&b%i==0){a/=i;b/=i}};print a,b}'
}

1 100

real 0m0.051s
user 0m0.058s
sys 0m0.009s

Pythonの場合:

まあ,たかだか2桁なので総当りで最大公約数で割りながら
一致を探しても計算量は知れているけど,
手計算でも解けそうというか,むしろ手計算で解く問題.

答えは4つと分かっている訳だし,少し書き出してみると,
同じ側の桁を消去すると,「1より小さい」を満たさない.
分子の左側,分母の右側を消去すると,明らかに取り得る値がない.
残るは分子の右側,分母の左側を消去した時.

\frac{10 a + b}{10 b + c} = \frac{a}{c}

これを計算すれば,「c=5」か「b-a=5」の時だけ.

c=5の時,a=1,2(a<5で3,4が不適なのは自明)とみていくと,
a = 1, 2
b = 9, 6
c = 5, 5

b-a=5の時,a<5なのでa=1, 2, 3, 4を計算すれば,
a = 1, 2, 3, 4
b = 6, 7, 8, 9
c = 4, 5.6, 6.857142857142857, 8
となる.

従って,

\frac{19}{95},\frac{26}{65},\frac{16}{64},\frac{49}{98}

「1 / (16/64 * 26/65 * 19/95 * 49/98)」は「100」
 
 

Project Euler 34 桁の階乗

Project Euler Q34 【桁の階乗】 – Qiita

%%bash
time {
seq 3 2540160 |
awk -v FS= '{printf $0" ";s=0;for(i=1;i<=NF;i++){x=1;for(j=1;j<=$i;j++){x*=j};s+=x};print s}' |
awk '$1==$2{s+=$1}END{print s}'
}

40730

real 0m12.620s
user 0m14.090s
sys 0m0.276s

Pythonの場合:

Q30と基本的には(というかほとんど)一緒.
辞書で変換テーブルを作る(再計算しない).

「math.factorial(2) + math.factorial(9) * 6 = 2177282」

で階乗の和が追いついていないので上限はそれ以下だと分かるが,
ここでは,取り敢えず上限を「100000」としてまずは計算してみる.
(必ずしも全数探索する必要はないよねの精神)

import math

def p34(max_iter=100000, disp=False):
  d = {i: math.factorial(i) for i in range(10)}
  res_s = set()

  res = 0
  for x in range(10, max_iter):
    q, rem = divmod(x, 10)
    if q in d:
      temp = d[q] + d[rem]
    else:
      temp = d[rem]
      while q > 0:
        q, rem = divmod(q, 10)
        temp += d[rem]
    d[x] = temp
    if x == temp:
      res += temp
      res_s.add(temp)
  if disp:    return res, res_s
  return res


p34()

40730

%timeit p34()

10 loops, best of 3: 46.6 ms per loop

ちなみに,

p34(max_iter=3000000, disp=True)

(40730, {145, 40585})

7桁までループすると処理の9割以上は無駄という悲しい事になる.

 
 

Project Euler 35 巡回素数

Project Euler Q35 【巡回素数】 – Qiita

%%bash
time {
seq 2 1000000 |
factor |
awk 'NF==2{print $2}' |
grep -v "[0|4|6|8]" |
grep -v "[2|5]." |
awk '{for(i=1;i<=length;i++){printf substr($1$1,i,length)" "};print ""}' |
while read line;do factor $line;echo "@";done |
tr '\n' ' ' |
tr '@' '\n' |
awk 'NF==(length($1)-1)*2' |
wc -l
}

55

real 0m1.113s
user 0m0.509s
sys 0m0.364s

Pythonの場合:

循環素数の問題.本気で高速化したい場合,Q27みたいにやるべきだろうけど,
100万未満の素数ってたかだか78498個しかないので,
愚直に回しても割と速いかもしれない.

100未満は既に求まっているので求めない.
偶数や5が含まれると,下1桁に来た時に割り切れる(合成数になる;約数がある)ので除く.
(0,2,4,5,6,8以外 → 1,3,7,9の順列という方向で攻めても良いかも)

from sympy import isprime, sieve

def p35():
  def circular(p):
    dt = []
    while p > 0:
      a = p % 10
      if a in [0, 2, 4, 5, 6, 8]:
        return 0
      else:
        dt.append(a)
      p //= 10
    res, n = [], len(dt)
    for i in range(n):
      temp = 0
      for j, x in enumerate(dt[i:] + dt[:i]):
        temp += x * 10**j
      if isprime(temp):
        res.append(temp)
      else:
        return 0
    return res

  mem = dict.fromkeys([2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, 97])
  for p in sieve.primerange(100, 1000000):
    if p not in mem:
      res = circular(p)
      if res:
        mem.update(dict.fromkeys(res))
  return len(mem)

p35()

55

%timeit p35()

10 loops, best of 3: 171 ms per loop

 
 

Project Euler 36 二種類の基数による回文数

Project Euler Q36 【二種類の基数による回文数】 – Qiita

%%bash
time {
seq 999999 |
paste -d" " - <(seq 999999 | rev) |
awk '$1==$2{print $1,int(log($1)/log(2))}' |
awk '{printf $1" ";v=$1;for(i=$2;0<=i;i--){printf int(v/2^i);v-=int(v/2^i)*2^i};print ""}' |
awk '{printf $0" ";s=1;for(i=1;i<=int(length($2)/2);i++){s*=(substr($2,i,1)==substr($2,length($2)+1-i,1))};print s}' |
awk '$3==1{s+=$1}END{print s}'
}

872187

real 0m1.074s
user 0m0.620s
sys 0m0.371s

Pythonの場合:

Pythonって文字列の操作が結構得意な言語だと思う.
簡単に処理できるって意味で,文字列の処理はかなり遅いけど.
例えば,Q35の循環素数を作るのも文字列に変換すれば簡単だし.
(遅くなると思うからやらなかったけど)

def p36():
  res = 0
  for i in range(1000000):
    s, bs = str(i), bin(i)[2:]
    if s == s[::-1] and bs == bs[::-1]:
      res += i
  return res

p36()

872187

%timeit p36()

1 loop, best of 3: 606 ms per loop

「bin(i)[2:]」って一見ダサいけど,
多分2進数に変換する方法としては一番効率が良い筈.
(format関数より20%位速くて,fプリフィックスより10%位速い)

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

bashで遊ぶ6 への4件のフィードバック

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

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

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

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

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