bashで遊ぶ8

はじめに

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

bashで遊ぶ6

Project Euler Q37 ー Project Euler Q43

bashで遊ぶ7

 

今回

Project Euler Q44 ー Project Euler Q49

Project Euler Q44 【五角数】 – Qiita

%%bash
time {
seq 9999 |
awk '{print $1*(3*$1-1)/2}' |
awk '{for(i=NR-1;1<=i;i--){v=i*(3*i-1)/2;p=$1+v;m=$1-v;if((1+sqrt(1+24*p))%6==0&&(1+sqrt(1+24*m))%6==0){print m;break}}}' |
sort -k1,1n |
head -1
}

5482660

real 0m14.218s
user 0m14.183s
sys 0m0.009s

五角数の性質を考えた時,小さい方から求めると
最初にみつかった時点で処理を打ち切る事ができる.
ただseqは指定した回数イテレートをするので,
seqで9999回イテレートしているから処理が凄く遅い.
実際は2,100-2,200位で答えが求まるので,

%%bash
time {
seq 2200 |
awk '{print $1*(3*$1-1)/2}' |
awk '{for(i=NR-1;1<=i;i--){v=i*(3*i-1)/2;p=$1+v;m=$1-v;if((1+sqrt(1+24*p))%6==0&&(1+sqrt(1+24*m))%6==0){print m;break}}}' |
sort -k1,1n |
head -1
}

5482660

real 0m0.670s
user 0m0.663s
sys 0m0.001s

seqはbreakできなさそうだからforに置き換えたらどうかと思った.

%%bash
time {
echo 1 | 
awk '{for(j=$1;j<=9999;j++){n=j*(3*j-1)/2;for(i=j-1;0<i;i--){v=i*(3*i-1)/2;p=n+v;m=n-v;if((1+sqrt(1+24*p))%6==0&&(1+sqrt(1+24*m))%6==0){print m;print i;print j;break2}}}}'
}

5482660
1020
2167

real 0m7.049s
user 0m7.042s
sys 0m0.002s

速くはなったけど,「seq 2200」と比べると全然速さが違う(´・ω・`)

Pythonの場合:

逆関数を出しても「ループを回す」という事に変わりはないので,
Python向きでは無い処理になるだろうなと.
この処理,Cなら数msオーダーなんだろうな.逆に云えば,
非常に単純なループ構造だからこれは絶対Cython案件…….
でも,ColaboratoryではCython使えないので略.

探索範囲は1億ループ以下なので,
Q31の時みたいにNumpyでクルードに解くのも可能.
(多分,500-800 msで解ける;過去の結果から推計)

ここではピュアに.

import math

def ispntg(val):
  return (1 + math.sqrt(1 + 24*val)) % 6 == 0


def p44(max_iter=9000):
  for i in range(1, max_iter):
    k = i * (3 * i - 1) // 2
    for l in range(i-1, 0, -1):
      j = l * (3 * l - 1) // 2
      if ispntg(k-j) and ispntg(k+j):
        return k - j, i, l


p44()

(5482660, 2167, 1020)

%timeit p44()

1 loop, best of 3: 1.24 s per loop

えー……おっそ……(´・ω・`)

line_profilerで確認すると,

!pip install -q line_profiler
%load_ext line_profiler
%lprun -f p44 p44()
Timer unit: 1e-06 s

Total time: 5.36466 s
File: 
Function: p44 at line 5

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     5                                           def p44(max_iter=9000):
     6      2167       1050.0      0.5      0.0    for i in range(1, max_iter):
     7      2167       1387.0      0.6      0.0      k = i * (3 * i - 1) // 2
     8   2348008    1010078.0      0.4     18.8      for l in range(i-1, 0, -1):
     9   2345842    1371464.0      0.6     25.6        j = l * (3 * l - 1) // 2
    10   2345842    2980678.0      1.3     55.6        if ispntg(k-j) and ispntg(k+j):
    11         1          1.0      1.0      0.0          return k - j, i, l

やっぱりって感じだけど.平方根の計算は重すぎる.

少しだけ簡単にできる範囲で最適化を考えると,
例えば内側のループだけnumpyで一括処理して,
trueが出た所で外側のイテレートを終えるとかでも良さそう.

import numpy as np


def ispntg(k, j):
  a = np.mod(1 + np.sqrt(1 + 24*(k-j)), 6) == 0
  b = np.mod(1 + np.sqrt(1 + 24*(k+j)), 6) == 0
  return np.logical_and(a, b)


def p44_2(max_iter=9000):
  for i in range(1, max_iter):
    k = i * (3 * i - 1) // 2
    l = np.arange(i-1, 0, -1)
    j = l * (3 * l - 1) // 2
    cond = ispntg(k, j)
    if cond.any():
      return k-j[cond.argmax()], i, l[cond.argmax()]
  return 'not found'


p44_2()

(5482660, 2167, 1020)

%timeit p44_2()

1 loop, best of 3: 190 ms per loop

ループな処理なので,numbaとかするだけで数十msになりそう.
(Colaboratory以下略)

 
 

Project Euler Q45 【三角数, 五角数, 六角数】 – Qiita

%%bash
time {
seq 99999 |
awk '{print $1*(2*$1-1)}' |
awk '{print $1,(1+sqrt(1+24*$1))%6}' |
awk '$2==0' |
awk '{print $1,(-1+sqrt(1+8*$1))%2}' |
awk '$2==0{print $1}' |
tail -1
}

1533776805

real 0m0.104s
user 0m0.146s
sys 0m0.011s

Pythonの場合:

単純な問題なので,愚直にループを回しても割とそんなに遅くはない.
試しに三角数,五角数,六角数を求めて差分集合をとっても,
遅い遅いと云われるPythonでも

10 loops, best of 3: 67 ms per loop

だった.

でもまあ,1つ前で五角数かチェックする関数を書いているんだから,
せっかくだからそれを利用する.
六角数なら三角数なのだから,六角数を求めて,
それをQ44の関数でチェックするだけ.

def ispntg(val):
  return (1 + math.sqrt(1 + 24*val)) % 6 == 0


def p45(max_iter=100000):
  for i in range(144, max_iter):
    h = i * (2 * i - 1)
    if ispntg(h):
      return h

p45()

1533776805

%timeit p45()

100 loops, best of 3: 14.4 ms per loop

Q44でnumpyを使っているし,これもnumpyで解けば,

import numpy as np

def p45_numpy(max_iter=100000):
  i = np.arange(144, max_iter)
  h = i * (2 * i - 1)
  return h[(np.mod(1 + np.sqrt(1 + 24*(h)), 6) == 0).argmax()]

p45_numpy(), p45_numpy(30000)

(1533776805, 1533776805)

%timeit p45_numpy()
%timeit p45_numpy(30000)

100 loops, best of 3: 2.04 ms per loop
1000 loops, best of 3: 624 µs per loop

単純な問題なのでかなり高速に解ける.

 
 

Project Euler Q46 【もうひとつのゴールドバッハの予想】 – Qiita

%%bash
time {
seq 3 2 9999 |
factor |
awk 'NF!=2{print $1*1}' |
awk '{printf $1" ";for(i=1;i<=(sqrt($1/2));i++){printf $1-2*i^2" "};print ""}' |
awk '{printf $1" "NF-1" ";for(i=2;i<=NF;i++){for(j=2;j<=sqrt($i);j++){if($i%j==0||$i==1){g[NR]+=1;break}}};print g[NR]}' |
awk '$2==$3{print $1}' |
head -1
}

5777

real 0m1.023s
user 0m1.087s
sys 0m0.009s

Pythonの場合:

import math
from sympy import sieve

def p46(max_iter=10000):
  sieve.extend(max_iter)
  def issquare(val):
    a = math.sqrt(val/2)
    return a == int(a)

  for i in range(33, max_iter, 2):
    for p in sieve.primerange(2, i+1):
      if issquare(i-p):
        break
    else:
      return i


p46()

5777

%timeit p46()

1 loop, best of 3: 273 ms per loop

import numpy as np
from sympy import sieve

def p46_numpy(max_iter=10000):
  def issquare(val):
    a = np.sqrt(val / 2)
    return a == a.astype(int)

  sieve.extend(max_iter)
  p_arr = np.array(sieve._list)

  for i in range(33, max_iter, 2):
    cond = i - p_arr[p_arr <= i]
    if issquare(cond).any():
      pass
    else:
      return i


p46_numpy()

5777

%timeit p46_numpy()

10 loops, best of 3: 41.9 ms per loop

 
 

Project Euler Q47 【異なる素因数】 – Qiita

%%bash
time {
seq 999999 |
factor |
awk '{delete a;for(i=2;i<=NF;i++){a[$i]=1};if(length(a)==4){print $0}}' |
awk '{print $1*1,$1+1,$1+2,$1+3}' |
factor |
awk '{delete a;for(i=2;i<=NF;i++){a[$i]=1};if(length(a)==4){print $0}}' |
awk '{print $1*1}' |
sort -k1,1n |
uniq -c |
awk '$1==4{print $2-3}' |
head -1
}

awk: line 1: illegal reference to array a
awk: line 1: illegal reference to array a

real 0m0.007s
user 0m0.002s
sys 0m0.004s

Pythonの場合:

下限は235*7で,そっからは強引にやるしか思い浮かばない.
sympy.factorintを使うと物凄く楽.ただし,物凄く遅い.

from sympy import factorint


def p47(max_iter=1000000):
  cnt = 0
  for i in range(210, max_iter):
    if len(factorint(i)) == 4:
      cnt += 1
      if cnt == 4:
        return i - 3
    else:
      cnt = 0
  return 'not found'


p47()

134043

%timeit p47()

1 loop, best of 3: 1.1 s per loop

そもそも,素因数分解をする必要なんて無くて,
4つの素数で割れるかどうかのチェックだけで良いので,
「sympy.factorint」を使うのはオーバーヘッドが大きい.

numpyを使って素数4つで割れる所だけピックアップして,
シーケンシャルパターンの検出には以前のテクニックを使おう.
(値が[0, 256)なら非常に高速(多分最速)に処理できる方法)

import numpy as np
from sympy import sieve

def isdiv(i, p_arr):
  cond = i[:, np.newaxis] / p_arr
  return cond == cond.astype(int)


def p47_2():
  i = np.arange(210, 200000)
  p_arr = np.array(sieve._list[:200])
  cond = (isdiv(i, p_arr).sum(1) == 4).astype(int).tolist()
  return 210 + bytes(cond).index(bytes((1, 1, 1, 1)))

p47_2()

134043

%timeit p47_2()

1 loop, best of 3: 862 ms per loop

思ったより遅い……(´・ω・`)

%lprun -f p47_2 p47_2()

Timer unit: 1e-06 s

Total time: 0.919054 s
File:
Function: p47_2 at line 9

Line # Hits Time Per Hit % Time Line Contents

 9                                           def p47_2():
10         1        726.0    726.0      0.1    i = np.arange(210, 200000)
11         1         39.0     39.0      0.0    p_arr = np.array(sieve._list[:200])
12         1     916410.0 916410.0     99.7    cond = (isdiv(i, p_arr).sum(1) == 4).astype(int).tolist()
13         1       1879.0   1879.0      0.2    return 210 + bytes(cond).index(bytes((1, 1, 1, 1)))
import numpy as np
from sympy import sieve


def p47_3():
  i = np.arange(210, 200000)
  p_arr = np.array(sieve._list[:200])
  cond1 = i[:, np.newaxis] / p_arr
  cond2 = (cond1 == cond1.astype(int)).sum(1)
  cond3 = np.convolve(cond2 == 4, [1, 1, 1, 1], 'valid')
  return 210 + cond3.argmax()

p47_3()

134043

%timeit p47_3()

1 loop, best of 3: 835 ms per loop

余り変わらない……(´・ω・`)

%lprun -f p47_3 p47_3()
Timer unit: 1e-06 s

Total time: 0.871085 s
File: 
Function: p47_3 at line 5

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
     5                                           def p47_3():
     6         1        636.0    636.0      0.1    i = np.arange(210, 200000)
     7         1         27.0     27.0      0.0    p_arr = np.array(sieve._list[:200])
     8         1     367772.0 367772.0     42.2    cond1 = i[:, np.newaxis] / p_arr
     9         1     500826.0 500826.0     57.5    cond2 = (cond1 == cond1.astype(int)).sum(1)
    10         1       1527.0   1527.0      0.2    cond3 = np.convolve(cond2 == 4, [1, 1, 1, 1], 'valid')
    11         1        297.0    297.0      0.0    return 210 + cond3.argmax()

計算コストの50%以上が「cond1 == cond1.astype(int)」……無駄過ぎる.
「numpy.remainder()」であまりの計算をする様に変更.

import numpy as np
from sympy import sieve


def p47_4():
  i = np.arange(210, 200000)
  p_arr = np.array(sieve._list[:200])
  cond1 = (np.remainder(i[:, np.newaxis], p_arr) == 0).sum(1)
  cond2 = np.convolve(cond1 == 4, [1, 1, 1, 1], 'valid')
  return 210 + cond2.argmax()

p47_4()

134043

%timeit p47_4()

1 loop, best of 3: 794 ms per loop

また,ちょっとだけ速くなった.

まあ,4000万ループに対してモジュロを求めてこの速度ならまあまあか.

strideを決めて途中で打ち切れば倍位は速くできそう.

 
 

Project Euler Q48 【自身のべき乗(self powers)】 – Qiita

%%bash
time {
seq 1000 |
awk '{a=1;for(i=1;i<=$1;i++){a=substr(a*$1,length(a*$1)-10+1,10)};s+=a}END{print substr(s,length(s)-10+1,10)}'
}

inf

real 0m0.730s
user 0m0.724s
sys 0m0.004s

ほげえ

mawkは31ビットが上限らしい.

Pythonの場合:

Pythonの場合,どうせ任意精度で計算される訳だから,
がっさーと全部足して文字列の最後10文字
抜き取る方がリーズナブルかなと.

例えば,Cみたいなループ処理が得意な言語で指数処理を落とし込めるなら,
10**10で割りながら加算すれば10桁が求まるだろうけど,
ループは苦手,任意精度で整数計算は高コストなPythonではかなり遅くなると思う.

取り敢えず,幾つか比較を.

def p48():
  res = 0
  for i in range(1, 1001):
    res += i**i
  s = str(res)
  return int(s[-10:])


def p48_2(mod = 10**10):
  return sum(i**i for i in range(1, 1001)) % mod


def p48_3(mod = 10**10):
  return sum(pow(i, i, mod) for i in range(1, 1001)) % mod


p48(), p48_2(), p48_3()

(9110846700, 9110846700, 9110846700)

%timeit p48()
%timeit p48_2()
%timeit p48_3()

100 loops, best of 3: 9.48 ms per loop
100 loops, best of 3: 9.38 ms per loop
100 loops, best of 3: 2.03 ms per loop

前言撤回.途中計算でもモジュロした方が圧倒的に速くなる.
勿論,ビルトイン関数を用いているからこそであって,
この処理をべた書きしたらその分遅くなると思うけど.

 
 
Project Euler Q49 【素数数列】 – Qiita

Pythonの場合:

問題がよく分からないので,試しに素数列に3330を足してみると,
足された結果も一定の間隔で素数になっている事が分かるので,
刻みは3330で良さそうと強い仮定を置いて解いてみる.

+3330の刻み幅なので,1つ目の一番小さい素数の探索範囲は,
一番小さい素数は[1488, 3000]位で良さそう.

from sympy import sieve, isprime

def p49():
  diff = 3330
  for p1 in sieve.primerange(1488, 3000):
    p2 = p1 + diff
    p3 = p2 + diff
    if isprime(p2) and isprime(p3):
      ps1 = set(str(p1))
      ps2 = set(str(p2))
      ps3 = set(str(p3))
      if (len(ps2.symmetric_difference(ps1)) == 0 and 
          len(ps3.symmetric_difference(ps1)) == 0):
        return str(p1) + str(p2) + str(p3)
  return 'not found'

p49()

‘296962999629’

%timeit p49()

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

 
 

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

bashで遊ぶ8 への2件のフィードバック

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

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

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