bashで遊ぶ10

はじめに

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 44 ー Project Euler 49

bashで遊ぶ8

Project Euler 50 ー Project Euler 53

bashで遊ぶ9

 

今回

Project Euler 54 ー Project Euler 59

そろそろ,これで一旦先に進めるのはやめて,一から見直していこうかな.
特にmawk(Colaboratoryに入っているオーク)ではできないのが多いので,
コードを参考にしつつ,改めて自分で一から考えたい.
(これ等エントリについて,ちょっとずつ加筆,修正していきたい)

 

Project Euler 54 ポーカーハンド

 
Project Euler Q54 【ポーカーハンド】 – Qiita

Colaboratoryにテキストファイルをダウンロードする.

!wget https://projecteuler.net/project/resources/p054_poker.txt

これもmawkでは実行できなかったので,また後日考えたい.

Pythonの場合:

8C TS KC 9H 4S 7D 2S 5D 3S AC
5C AD 5D AC 9C 7C 5H 8D TD KS
3H 7H 6S KC JS QH TD JC 2D 8S
TH 8H 5C QS TC 9H 4D JC KS JS
7C 5H KC QH JD AS KH 4C AD 4S
5H KS 9C 7D 9H 8D 3S 5D 5C AH
・
・
・

という形で非常に綺麗なテキストデータなので,
if文に流し込んであげれば解けるし,
1000行しかないテキストデータなので,
どんな処理をしても数十msオーダーで解けそうだけど,
実にPythonらしくないというか…….

でもまあ,他にやり方思い浮かばない.

from collections import Counter


def p54(strings):
  def calc(val, suits):
    s_val, s_suits, c = set(val), set(suits), Counter(val)
    k, v = list(zip(*sorted(c.items(), key=lambda x: (-x[1], -x[0]))))
    if len(s_suits) == 1: #All cards of the same suit
      if all(val[i-1]+1 == val[i] for i in range(1, 5)):
        return 1000 if 14 in s_val else 800 + k[0] #Royal or Straight flush
      else:
        return 500 + k[0] #Flush
    else:
      if len(s_val) == 2:
        return 700 + k[0] if 4 in v else 600 + k[0] #Four card or Full house
      elif len(set(val)) == 3:
        return 300 + k[0] if 3 in v else 200 + k[0] #Three card or Two pair
      elif len(s_val) == 4:
        return 100 + k[0] #One pair
      else:#Straight or High card
        return 400 + k[0] if all(val[i-1]+1 == val[i] for i in range(1, 5)) else k[0]

  d = {'2': 2, '3': 3, '4': 4, '5': 5, '6': 6, '7': 7, '8': 8, 
          '9': 9, 'T': 10, 'J': 11, 'Q': 12, 'K': 13, 'A': 14}

  res = 0
  for lines in slst:
    x = [[(s[0], s[1]) for s in lines]]
    val, suits = [list(zip(*subli)) for subli in x][0]
    val1, suits1 = sorted(d[k] for k in val[:5]), suits[:5]
    val2, suits2 = sorted(d[k] for k in val[5:]), suits[5:]
    pt1 = calc(val1, suits1)
    pt2 = calc(val2, suits2)
    if pt1 > pt2:
      res += 1
  return res


with open('p054_poker.txt', 'r') as f:
  text = f.read()

strings = text.splitlines()
slst = [s.split() for s in text.splitlines()]
p54(slst)

376

%timeit p54(slst)

10 loops, best of 3: 25.2 ms per loop

 
 

Project Euler 55 Lychrel数

Project Euler Q55 【Lychrel数】 – Qiita

%%bash
time {
seq 9999 |
awk '{n[NR,0]=$1;printf $1" ";for(i=1;i<=49;i++){for(j=length(n[NR,i-1]);0<j;j--){t[NR,i]=(t[NR,i] substr(n[NR,i-1],j,1))};n[NR,i]=n[NR,i-1]+t[NR,i];printf n[NR,i]" "};print ""}' |
awk '{for(i=2;i<=NF;i++){for(j=length($i);0<j;j--){t[NR,i]=(t[NR,i] substr($i,j,1))};if($i==t[NR,i]){next}};print $1}' |
wc -l
}

586

real 0m14.072s
user 0m14.723s
sys 0m0.127s

これもmawkだから解けない系かな.
mawkは計算できる上限が31bit位な上,全てを浮動小数で計算するとか何とか.

Pythonの場合:

Lychrel number – Wikipedia

196までLychrel数が無いという事は既知の事実らしいので,
事前情報として用いて,探索範囲は[197, 10000)とする.
(まあ処理時間に有意な差はみられないが)

def p55():
  def isLychrel(val):
    for _ in range(50):
      val += int(str(val)[::-1])
      if val == int(str(val)[::-1]):
        return 0
    return 1

  return sum(map(isLychrel, range(197, 10000))) + 1


p55()

249

%timeit p55()

10 loops, best of 3: 56.8 ms per loop

 
 

Project Euler 56 もっとべき乗の数字和

Project Euler Q56 【もっとべき乗の数字和】 – Qiita

mawkではry

Pythonの場合:

from itertools import permutations

def p56():
  return max(sum(map(int, str(a**b))) for a, b in permutations(range(2, 100), 2))

p56()

972

%timeit p56()

10 loops, best of 3: 137 ms per loop

非常にシンプルに解けるけど,重複計算していて効率が悪い.
整数で処理する形に変更し,メモ化してあげる.
ただし,桁出しの関数をそのままメモ化すると
メモリがサチって答えおかしくなるので,
入力(整数)と出力(その桁を合計したもの)を
対応させた辞書(メモ)だけになる様にする.

from itertools import permutations
from functools import lru_cache


def gen(val):
  if val < 10:
    yield val
  else:
    yield from gen(val//10)
    yield val % 10


@lru_cache(maxsize=None)
def getintdt(val):
  return sum(gen(val))


def p56_2():
  return max(getintdt(a**b) for a, b in permutations(range(2, 100), 2))


p56_2()

972

%timeit p56_2()

100 loops, best of 3: 8.04 ms per loop

Project Euler 57 平方根の近似分数

Project Euler Q57 【平方根の近似分数】 – Qiita

%%bash
time {
seq 1000 |
awk 'BEGIN{s=b=1} {c=b+s;t=2*b+s;print c,t;b=c;s=t}' |
awk 'length($1)<length($2)' |
wc -l
}

72

real 0m0.004s
user 0m0.004s
sys 0m0.000s

mawkは31bitまでしか計算できないので,答えが出ない(途中でサチる).

awkを使わなかったらどうなんだろう?と思って,
自分でも書いてみたけど(上のを展開しただけだけど),

%%bash
time {
s=1
b=1
cnt=0
for i in `seq 1000`
do
c=$(($b+$s))
t=$((2*$b+$s))
if [ ${#t} -gt ${#c} ]
then
cnt=$((cnt+1))
fi
b=$c
s=$t
done
echo $cnt
}

256

real 0m0.023s
user 0m0.020s
sys 0m0.002s

途中から負値が飛び飛びに出ているので,オーバーフローを起こしている.
算術式展開の場合は計算できる値の上限は「2**63 – 1」っぽいな.
任意精度で処理する方法はさっぱりだけど,
これまでにmawkで解けなかった問題については,
置き換えれば一部解けそうだなという事は分かった.

Pythonの場合:

def p57(nmr=1393, dnm=985, n=9, res=1):
  for _ in range(n, 1001):
    nmr = 2 * dnm + nmr
    dnm = nmr - dnm
    if len(str(nmr)) > len(str(dnm)):
      res += 1
  return res

p57()

153

%timeit p57()

100 loops, best of 3: 2.89 ms per loop

 
 

Project Euler 58 螺旋素数

Project Euler Q58 【螺旋素数】 – Qiita

%%bash
time {
seq 99999 |
awk '{print 4*$1^2-2*$1+1,4*$1^2+4*$1+1,4*$1^2+1,4*$1^2+2*$1+1}' |
factor |
awk '{print NF==2}' |
paste -d" " - - - - |
awk '{for(i=1;i<=4;i++){s+=$i};print s/(1+4*NR)}' |
awk '$1<0.1{print 2*NR+1}' |
head -1
}

26241

real 0m0.126s
user 0m0.162s
sys 0m0.022s

「factor」は本当に速いなあ.
Cプリミティブと同等な速度を叩き出す様に感じる.

Pythonの場合:

Q28と同じ様な問題で,こちらは完全にウラムの螺旋の話.

ウラムの螺旋

ダウンロード (2)

図をみると分かり易いが,素数列の出現には規則性があって,
対角,水平,垂直方向に密集して何かしらの模様にみえる.

Q28で対角線上の値を求めているので,そのコードをそのまま流用して,
それが素数かどうか判定して割合が10%未満になった所で打ち切れば良い.

元のコードが汚いので酷い感じになるけど.

from sympy import isprime

def p58():
  ludiff, rudiff, lldiff, rldiff, diffdiff = 4, 2, 8, 6, 8
  lu = 1 + ludiff # 5 is prime
  ru = 1 + rudiff # 3 is prime
  ll =  1 + lldiff # 9
  rl = 1 + rldiff # 7 is prime
  total, cnt = 4, 3
  for i in range(3, 15000):
    ludiff += diffdiff
    rudiff += diffdiff
    lldiff += diffdiff
    rldiff += diffdiff
    lu += ludiff
    ru += rudiff
    ll += lldiff
    rl += rldiff
    total += 4
    cnt += isprime(lu) + isprime(ru) + isprime(ll) + isprime(rl)
    if cnt / total < 0.1:
      break
  return 2 * (i - 2) + 1
 

p58()

26241

%timeit p58()

1 loop, best of 3: 295 ms per loop

綺麗にするとこうなるけど,

from sympy import isprime

def p58_2():
  cnt = 0
  for i in range(1, 15000):
    cnt += (isprime(4*i**2-2*i+1) 
               + isprime(4*i**2+4*i+1)
               + isprime(4*i**2+1) 
               + isprime(4*i**2+2*i+1))
    if cnt / (4*i+1) < 0.1:
      break
  return 2 * i + 1
 

p58_2()

26241

%timeit p58_2()

1 loop, best of 3: 294 ms per loop

isprimeがボトルネックなので,速さは余り変わらない.
ただ,isprimeは凄くPythonicなコードで,
ピュアにやろうとしたらこれがリーズナブルチョイスだろう.
自分で実装したって劣化版しかできない.絶対に.
素数列が数十万以下であればin sieveの方が速いかも.
今回みたいに何十億って値まで素数判定する場合は,
isprimeの方が遥かに効率的(というかin sieveできないだろう).

最初の方では厳密に素数判定する必要はそもそもなくて,
確率的に素数を判定すればもっと高速化できそう.

フェルマーテストでどうだろう?と思ったけど,

import math

def Fermat_test(val, a=5, b=17):
  if val == 2: 
    return True
  elif val%2 == 0 or val == 1:
    return False
  elif math.gcd(a, val) != 1:
    return False
  elif pow(a, val-1, val) != 1:
    return False
  elif math.gcd(b, val) != 1:
    return False
  elif pow(b, val-1, val) != 1:
    return False
  return True


def p58_3():
  cnt = 0
  for i in range(1, 15000):
    cnt += (Fermat_test(4*i**2-2*i+1) 
               + Fermat_test(4*i**2+4*i+1)
               + Fermat_test(4*i**2+1) 
               + Fermat_test(4*i**2+2*i+1))
    if cnt / (4*i+1) < 0.1:
      break
  return 2 * i + 1
 

p58_3()

26247

%timeit p58_3()

10 loops, best of 3: 211 ms per loop

精度はあまり良くなくてたまたま一致しただけ.
そして,速さがはっきり言ってそんなに変わらない…….

SympyにはMiller-Rabin素数判定法も用意されている様なので試してみる.

参考:
ミラー–ラビン素数判定法 – Wikipedia

Ntheory Class Reference – SymPy 1.1.1 documentation

from sympy.ntheory.primetest import mr

def p58_4():
  cnt = 0
  for i in range(1, 15000):
    cnt += (mr(4*i**2-2*i+1, [2, 3]) 
               + mr(4*i**2+4*i+1, [2, 3])
               + mr(4*i**2+1, [2, 3]) 
               + mr(4*i**2+2*i+1, [2, 3]))
    if cnt / (4*i+1) < 0.1:
      break
  return 2 * i + 1
 

p58_4()

26241

%timeit p58_4()

1 loop, best of 3: 496 ms per loop

遅くなった.

「from sympy.ntheory.primetest.isprime」
の実効性の高さを再確認しただけだった.
手間を掛けずに高速化を図るのは無理そう.

「sympy.ntheory.primetest.isprime」のソースコードをみたら当たり前だった.
まず分かりきっている範囲を決定論的に調べ,
値が非常に大きい範囲では確率論的(Miller-Rabin)に調べている.

以前,「in sieve」の方が速いと云ったけど,数万以下であれば,
isprimeでもsieve.searchによる判別があるけど,
それまでに判別条件を通る過程で幾分か遅くなるが,
一定以上の値だと「in sieve」は使い物にならない位遅い.

 
 

Project Euler 59 XOR暗号解読

Problem 59 「XOR暗号解読」 – PukiWiki

(追記,加筆修正予定)

Pythonの場合:

!wget https://projecteuler.net/project/resources/p059_cipher.txt

してColaboratoryに暗号文をダウンロードする.

%%bash
cat p059_cipher.txt |
wc -m

3,203

3,203文字程度の非常に短い文章.
何かの一節かアブストラクトだろう.

この手の頻出単語は「スペース」か,e,a,i,o,tとかなんかその辺.
事前にスペースが削除されていなければ,
そして多分元の文章を(スペースを削除したり難読化の為の処理をせず)
そのまま暗号化しているとすれば,
十中八九スペースに当たりをつけてやれば良さそう.
(英語の頻出ワードや頻度分析について,Wikipediaが詳しい)

from collections import Counter


def p59(cphlst):
  key_l = 3
  splst = [[cval for i, cval in enumerate(slst) if i % key_l == j] for j in range(key_l)]
  encrypt_keys = [Counter(cval).most_common(1)[0][0] for cval in splst]
  dcrypt_keys = [ord(' ') ^ i for i in encrypt_keys]
  key_word = ''.join(chr(cval) for cval in dcrypt_keys)
  
  res, decrypted_text = 0, ''
  for i, cval in enumerate(slst):
    dcrypt_ord = dcrypt_keys[i%3] ^ cval
    res += dcrypt_ord
    decrypted_text += chr(dcrypt_ord)
  return res, key_word, decrypted_text


with open('p059_cipher.txt', 'r') as f:
  text = f.read()

cphlst = list(map(int, text.split(',')))
p59(cphlst)

(107359,
‘god’,
“(The Gospel of John, chapter 1) 1 In the beginning the Word already existed. He was with God, and he was God. 2 He was in the beginning with God. 3 He created everything there is. Nothing exists that he didn’t make. 4 Life itself was in him, and this life gives light to everyone. 5 The light shines through the darkness, and the darkness can never extinguish it. 6 God sent John the Baptist 7 to tell everyone about the light so that everyone might believe because of his testimony. 8 John himself was not the light; he was only a witness to the light. 9 The one who is the true light, who gives light to everyone, was going to come into the world. 10 But although the world was made through him, the world didn’t recognize him when he came. 11 Even in his own land and among his own people, he was not accepted. 12 But to all who believed him and accepted him, he gave the right to become children of God. 13 They are reborn! This is not a physical birth resulting from human passion or plan, this rebirth comes from God.14 So the Word became human and lived here on earth among us. He was full of unfailing love and faithfulness. And we have seen his glory, the glory of the only Son of the Father.”)

%timeit p59(cphlst)

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

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