bashで遊ぶ5

はじめに

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

 
 

Project Euler 26 逆数の循環節 その1

Project Euler Q26 【逆数の循環節 その1】 – Qiita

-Mオプションが使えないと計算に失敗する

%%bash
time {
seq 3 999 |
grep -v "^5$" |
factor |
awk 'NF==2{printf $2" ";for(i=1;i<=999;i++){v=(10^i)%$2;if(v==1){printf i" ";break}};print ""}' |
sort -k2,2n |
tail -1
}

997 309

real 0m0.019s
user 0m0.015s
sys 0m0.002s

%%bash
time {
seq 3 999 | factor | egrep -v '2|5' |
awk 'NF==2{printf $2" ";for(i=1;i<=999;i++){v=(10^i)%$2;if(v==1){printf i" ";break}};print ""}' |
sort -k2,2n |
tail -1

}

997 309

real 0m0.016s
user 0m0.009s
sys 0m0.005s

Pythonの場合:

フェルマーの小定理とかさっぱり覚えていないけど
要は分母dに素数と2,5の倍数が含まれていなければ循環節で
その長さは(10進数の場合)10**nをdでひたすら割れば良いと

import sys
from functools import lru_cache

sys.setrecursionlimit(10000)


@lru_cache(maxsize=None)
def f(d, i=1, m=10):
  if m % d != 1:
    return f(d, i+1, m*10)
  else:
    return i


def p26():
  res = 0
  for d in range(2, 1000):
    if d % 2 == 0 or d % 5 == 0:
      pass
    else:
      temp = f(d)
      if temp > res:
        res = temp
  return res


p26()

982

%timeit p26()

10000 loops, best of 3: 172 µs per loop

 
 

Project Euler 27 二次式素数

Project Euler Q27 【二次式素数】 – Qiita

%%bash
time {
seq 3 1000 |
factor |
awk 'NF==2' |
awk '{for(i=-499;i<=500;i++){v=2*i-1;if(2<=1+v+$2){print v,$2}}}' |
awk '{for(i=0;i<=100;i++){v=i^2+$1*i+$2;if(v<2){break};for(j=2;j<=int(sqrt(v))+1;j++){if(v%j==0){next}};print i,v,$1,$2}}' |
uniq -f2 -c |
sort -n |
tail -1 |
awk '{print $4*$5}'
}

-59231

real 0m0.940s
user 0m1.039s
sys 0m0.056s

Pythonの場合:

1.n=0が素数 → bは素数

bは[2, 1000)の素数

2.n=2が素数 → 2a+b+4が素数 → bは奇数

bは[3, 1000)の素数

3.n=1が素数 → a+b+1が素数 → aは奇数かつa>-b

a>-bでどちらも奇数なのでaは[2-b, 1000)の奇数

そして

{n}^{2} + {a}^{n} + b

という式の形をみると「n = b」の時

n \left(a + n + 1\right)

になって式の結果がbの倍数になる
合成数になって素数じゃなくなるので
それ以上ループを回す必要性がない事が分かる

from sympy import isprime, sieve


def f(n, a, b):
  return n**2 + a*n + b


def p27():
  max_a, max_b = 0, 0
  res = 0
  for b in sieve.primerange(3, 1000):
    for a in range(2-b, 1000, 2):
      temp = 0
      for n in range(b):
        if isprime(f(n, a, b)):
          temp += 1
        else:
          break
      if temp > res:
        res = temp
        max_a, max_b = a, b
  return max_a * max_b

p27()

-59231

%timeit p27()

1 loop, best of 3: 496 ms per loop

まあここまでやっても結構遅いんだけどね(´・ω・`)
(「sieve.primerange」や「isprime」がボトルネック)

Numpyでクルードに解いたら10msオーダー出せそう
(コードを書き換えるのが面倒なのでやらない)

せっかくシンプルなループ構造なので
Cythonでどこまで高速化できるか試したい

実行環境がオンボロPCなのでColaboratoryの結果と比較はしたら駄目

上記コードをオンボロPCでそのまま実行したら

%timeit p27()

1.07 s ± 7.76 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

これを基準に何%高速化できるかをみたい

%%cython
from sympy import sieve
import array
from cpython cimport array

cdef int f(int n, int a, int b):
    return n**2 + a*n + b


cdef int _p27_cython(int[:] prime_arr, int p_num):
    cdef int max_a = 0
    cdef int max_b = 0
    cdef int temp = 0 
    cdef int res = 0
    cdef int a, b, bi, start, temp_val
    for bi in range(1, p_num-1):
        b = prime_arr[bi]
        start = 2 - b
        for a in range(start, 1000, 2):
            temp = 0
            for n in range(b):
                temp_val = f(n, a, b)
                for bi in range(1, p_num-1):
                    if temp_val == prime_arr[bi]:
                        temp += 1
                        break
                else:
                    break
            if temp > res:
                res = temp
                max_a, max_b = a, b
    return max_a * max_b


def p27_cython():
    sieve.extend(1000)
    cdef array.array p_arr = sieve._list
    cdef int[:] cp_arr = p_arr
    num = len(p_arr)
    return _p27_cython(cp_arr, num)
p27_cython()

-59231

%timeit p27_cython()

65.1 ms ± 877 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

これを参考に

Working with Python arrays – Cython 0.28a0 documentation

簡単にちょちょっとやってこんなもん
(というか僕にはこれが限界……(´・ω・`))

Sympy.sieveがボトルネックになっているのと
array.arrayの部分をvectorに置き換えたらもっと速くなるかなと
手間を掛けずにもうちょっと高速化する場合
SympyのソースコードをコピペしてCythonで最適化すれば
Cで書く並みの速さにならないかなって(思うだけでやらない!)

何にせよ,非常に簡単に

1.07 s ± 7.76 ms per loop

65.1 ms ± 877 µs per loop

になった

ColaboratoryでCythonが使えたらもっと面白いのになあ……

いやホントColaboratory流行れ
(そして調べたら何でも出てくる様になって欲しい……)
 
 

 
 

Project Euler 28 螺旋状に並んだ数の対角線

Project Euler Q28 【螺旋状に並んだ数の対角線】 – Qiita

%%bash
time {
seq -s" " 2 $((1001**2)) |
awk '{for(i=1;i<=NF;i++){printf $i" ";if($i==int(sqrt($i))^2 && ($i-1)%8==0){print ""}}}' |
awk 'BEGIN{s=1}{for(i=1;i<=4;i++){s+=$(2*i*NR)}}END{print s}'
}

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

real 0m0.300s
user 0m0.273s
sys 0m0.049s

むう

mawkって制約多い?というか単純に古いのだろうか……

「MSYS2/MinGW-w64」に入っているオークは

$ awk –version
GNU Awk 4.1.4, API: 1.1 (GNU MPFR 3.1.5-p2, GNU MP 6.1.2)
Copyright (C) 1989, 1991-2016 Free Software Foundation.

だったので試してみる


$ time {
> seq -s" " 2 $((1001**2)) |
> awk '{for(i=1;i awk 'BEGIN{s=1}{for(i=1;i }

669171001

real 0m3.170s
user 0m2.388s
sys 0m0.342s

Pythonの場合:

あーなんだっけあれを思い出す.
ドラキュラ伯爵みたいな名前のやつ.
ヴラドの螺旋?

これだ

ウラムの螺旋

今の場合は対角線の合計を出すので右回りか左回りか関係ないのでコードを流用してみる.

(以下,過去のコードを流用するため左回りで考えている)

def cell(n, x, y):
  d = 0
  x = x - (n - 1) // 2
  y = y - n // 2
  l = 2 * max(abs(x), abs(y))
  d = (l * 3 + x + y) if y >= x else (l - x - y)
  return (l - 1)**2 + d


def ulam(n=1001):
  return [[v for v in [cell(n, x, y) for x in range(n)]] for y in range(n)]


def p28():
  ulam_lst = ulam()
  cl, cr = 0, len(ulam_lst[0]) - 1
  res = 0
  for rlst in ulam_lst:
    a, b = rlst[cl], rlst[cr]
    if a == b:
      res += a
    else:
      res += a + b
    cl += 1
    cr -= 1
  return res

p28()

669171001

%timeit p28()

1 loop, best of 3: 1.25 s per loop

おっそ……明らかに螺旋リストがボトルネックになっている.

ちょっとだけちゃんと考えると,11 x 11を眺めると

ulam(11)

[[101, 100, 99, 98, 97, 96, 95, 94, 93, 92, 91],
[102, 65, 64, 63, 62, 61, 60, 59, 58, 57, 90],
[103, 66, 37, 36, 35, 34, 33, 32, 31, 56, 89],
[104, 67, 38, 17, 16, 15, 14, 13, 30, 55, 88],
[105, 68, 39, 18, 5, 4, 3, 12, 29, 54, 87],
[106, 69, 40, 19, 6, 1, 2, 11, 28, 53, 86],
[107, 70, 41, 20, 7, 8, 9, 10, 27, 52, 85],
[108, 71, 42, 21, 22, 23, 24, 25, 26, 51, 84],
[109, 72, 43, 44, 45, 46, 47, 48, 49, 50, 83],
[110, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82],
[111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121]]

対角線に着目すると

a = [[lst[i], lst[10-i]] for i, lst in enumerate(ulam(11))]
b = [[a[i+1][0]-lst[0], a[i+1][1]-lst[1]] for i, lst in enumerate(a[:-1])]
print(a)
print(b)

[[101, 91], [65, 57], [37, 31], [17, 13], [5, 3], [1, 1], [9, 7], [25, 21], [49, 43], [81, 73], [121, 111]]
[[-36, -34], [-28, -26], [-20, -18], [-12, -10], [-4, -2], [8, 6], [16, 14], [24, 22], [32, 30], [40, 38]]

それぞれ階差数列になっている事が分かる.

これをこのまま書き下すと,

def p28_2():
  ludiff, rudiff = 4, 2
  lldiff, rldiff = 8, 6
  lu = 1 + ludiff
  ru = 1 + rudiff
  ll =  1 + lldiff
  rl = 1 + rldiff
  diffdiff = 8
  res = 1 + lu + ru + ll + rl
  for _ in range(1, 500):
    ludiff += diffdiff
    rudiff += diffdiff
    lldiff += diffdiff
    rldiff += diffdiff
    lu += ludiff
    ru += rudiff
    ll += lldiff
    rl += rldiff
    res += lu + ru + ll + rl
  return res

p28_2()

669171001

%timeit p28_2()

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

ちょっと汚いコードなので,少しプリティにすると,

def p28_3():
  ludiff, rudiff, lldiff, rldiff, diffdiff = 4, 2, 8, 6, 8
  adn = ludiff + rudiff + lldiff + rldiff
  temp = 4 + adn
  res = 1 + temp
  for _ in range(1, 500):
    adn += 4 * diffdiff
    temp +=  adn
    res += temp
  return res

p28_3()

669171001

%timeit p28_3()

10000 loops, best of 3: 42.6 µs per loop

 
 

Project Euler 29 個別のべき乗

Project Euler Q29 【個別のべき乗】 – Qiita

Colaboratoryのmawkでは-Mオプションが使用できないので答えが合わない.

%%bash
time {
seq 2 100 |
awk '{for(i=2;i<=100;i++){print $1^i}}' |
sort -nu |
wc -l
}

8275

real 0m0.021s
user 0m0.020s
sys 0m0.003s

「MSYS2/MinGW-w64」のgawkで試す.


$ time {
> seq 2 100 |
> awk '{for(i=2;i sort -nu |
> wc -l
> }
9220

real    0m0.587s
user    0m0.374s
sys     0m0.030s

 
 
「sort -nu」だと結果がおかしかったけど,
「sort | uniq」だとそれらしい答えが求まる.

%%bash
time {
seq 2 100 |
awk '{for(i=2;i<=100;i++){print $1^i}}' | 
sort -n | 
uniq |
wc -l
}

9181

real 0m0.021s
user 0m0.016s
sys 0m0.003s

%%bash
time {
seq 2 100 |
awk '{for(i=2;i<=100;i++){print $1^i}}' | 
sort -g | 
uniq |
wc -l
}

9181

real 0m0.050s
user 0m0.045s
sys 0m0.004s

%%bash
time {
seq 2 100 |
awk '{for(i=2;i<=100;i++){print $1^i}}' | 
sort -n | 
uniq |
awk 'END{print NR}'
}

9181

real 0m0.024s
user 0m0.016s
sys 0m0.004s

多分,指数表記で丸め込まれてダブっている部分があるのかなと.

%%bash
time {
seq 2 100 |
awk '{for(i=2;i<=100;i++){printf("%.f\n", $1^i)}}' |
sort -n | 
uniq |
wc -l
}

9183

real 0m0.084s
user 0m0.073s
sys 0m0.008s

%%bash
time {
seq 2 100 |
awk '{for(i=2;i<=100;i++){printf("%.f\n", $1^i)}}' |
sort -n | 
uniq |
awk 'END{print NR}'
}

9183

real 0m0.087s
user 0m0.078s
sys 0m0.006s

ビンゴ!

mawkは計算できないできないと思っていたけど,
意外に大きな数字まで計算ができるんだなあと.
浮動小数の問題は常に付きまとうけども,今までの問題であった
「gawkではできるけど,mawkではできない」の多くは,
「gawkは整数は整数だけど,mawkは全部浮動小数」からくる問題で,
特に浮動小数点数の表現と誤差の問題というよりは,
指数表記になっていて意図した比較或いは処理ができていないという事にありそう.

Pythonの場合:

import itertools

def p29():
  return len(set(a**b for a, b in itertools.product(range(2, 101), repeat=2)))


p29()

9183

%timeit p29()

100 loops, best of 3: 7.8 ms per loop

 
 

Project Euler 30 各桁の5乗

Project Euler Q30 【各桁の5乗】 – Qiita

%%bash
time {
seq 999999 |
awk -v FS= '{s=0;for(i=1;i<=NF;i++){s+=$i^5;printf $i};printf " "s"\n"}' |
awk 'BEGIN{s=-1}$1==$2{s+=$1}END{print s}'
}

443839

real 0m1.538s
user 0m1.841s
sys 0m0.090s

Pythonの場合:

全然良い方法が思い浮かばないので愚直に考える.

まず,探索範囲は10からmax_iter = 9**5 * 6だろう.

また,何度も累乗の計算をするのも無駄なので,
辞書で変換テーブルを作り,途中計算も辞書に追加しながら
ながら処理を進めれば多少はマシだろう.

def p30():
  max_iter = 9**5 * 6
  d = {i: i**5 for i in range(10)}

  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
  return res


p30()
%timeit p30()

10 loops, best of 3: 177 ms per loop

辞書を使うとCythonで最適化するのが凄く難しいけど,
そもそもCで書けば途中計算保存しなくても10倍位は速くなるんだろうな.

 
 
 

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