2012年11月19日月曜日

プロジェクトオイラー 問65

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。


問65「e についての連分数である近似分数の100項目の分子の桁の合計を求めよ」
2の平方根は無限連分数として書くことができる。
            1
√2 = 1 + -------------------
                1
          2 + ---------------
                    1
              2 + -----------
                        1
                  2 + -------
                      2 + ...

無限連分数である√2 = [1;(2)]と書くことができるが、
(2)は2が無限に繰り返されることを示す。

同様に、√23 = [4;(1,3,1,8)]。
平方根の部分的な連分数の数列から良い有理近似が得られることが分かる。

√2の近似分数について考えよう。
    1   2
1 + - = -
    2   3

      1     7
1 + ----- = -
        1   5
    2 + -
        2

      1         17
1 + --------- = --
          1     12
    2 + -----
            1
        2 + -
            2

      1             41
1 + ------------- = --
          1         29
    2 + ---------
             1
        2 + -----
                1
            2 + -
                2

従って、√2の近似分数からなる数列の最初の10項は:
1, 3/2, 7/5, 17/12, 41/29, 99/70, 239/169, 577/408, 1393/985, 3363/2378, ...


もっとも驚くべきことに、数学的に重要な定数、
e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]。

eの近似分数からなる数列の最初の10項は:
2, 3, 8/3, 11/4, 19/7, 87/32, 106/39, 193/71, 1264/465, 1457/536, ...

10項目の近似分数の分子の桁を合計すると1+4+5+7=17である。
eについての連分数である近似分数の100項目の分子の桁の合計を求めよ。
-----



注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。





私の解答例は以下です。畳んでいます。
def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a

def h(n):
	if n<0: return 0
	L = [2]
	for i in xrange(n-1):
		L.append({True:(i//3+1)*2, False:1}[i%3==1])
	return L

def g(L):
	b, c = L.pop(), 1
	for a in L[::-1]:
		x, y = a*b+c, b
		s = gcd(x, y)
		b, c = x//s, y//s
	return b, c

def f(n):
	L = h(n)
	x, y = g(L)
	s = str(x)
	M = [int(s[i]) for i in xrange(len(s))]
	return sum(M), x, y

n = 100
print f(n)


4つの関数を使っています。

1.関数gcd(a, b)
・「ユークリッドの互除法」を使用して、引数a,bの最大公約数を返します。
 説明は問5を参照してください。


2.関数h(n)
・n項までのネイピア数の連分数リストを返します。
 但し問題文ではセミコロンで示している区切り文字も全てカンマ区切りにします。
 e = [2; 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]
  →[2, 1,2,1, 1,4,1, 1,6,1 , ... , 1,2k,1, ...]


・if n<0: return 0
 nが0未満は未定義なので0を返します。


・L = [2]
 初期値は連分数の固定部分なので2です。


・for i in xrange(n-1):
 ループ変数iには0以上n-1未満のn個の整数が順に設定されます。


・L.append({True:(i//3+1)*2, False:1}[i%3==1])
 「1,2k,1, ...」を発生させ、リストLに設定します。
 %演算子は割り算の余りです。
 [i%3==1]で3つごとに余り1になればTrue、他はFalseです。
 これを前の{}内の辞書要素のうち当てはまる値を設定します。
 //演算子は割り算の商です。True、つまり、3で割り余り1になる場合、
 「(i//3+1)*2」で上記式の2kにあたる部分の値を設定し、その他は1を設定します。


3.関数g(L):
・リストLとして連分数の固定部分aと循環部分b1,b2,...を受け取り、
 その近似分数x/yの分子x、分母yを返します。式で書くと以下の通りです。
      1             x
a + ------------- = -
           1        y
    b1 + --------
               1
         b2 + ---
              ...


・b, c = L.pop(), 1
 リストのpop関数はリストの要素を取り出して、それをリストから削除します。
 引数未指定の場合は最後の要素を取り出して、それをリストから削除します。
 初期値として、bに連分数の循環部分の最後、cに1を設定します。


・for a in L[::-1]:
 リストLを逆順に1つずつ取り出しループ変数aに設定します。


・x, y = a*b+c, b
    c   a*b+c   x
a + - = ----- = -
    b     b     y
上記式を実装しました。


・s = gcd(x, y)
 計算したx,yが約分できることがあるので、関数gcdで最大公約数を求めておきます。


・b, c = x//s, y//s
 先ほど求めたx,yをその最大公約数で割り、できる限り簡単な分数に直します。
 これを次のb,cとしてループ変数aのループを最後まで回します。


・return b, c
 連分数リストLから算出したその近似分数の分子と分母を返します。

4.関数f(n):
・ネイピア数のn項目までの連分数の近似分数の分子の各桁の合計、及びその分子、分母を返します。


・L = h(n)
 関数h(n)でn項までのネイピア数の連分数リストを取得しLとします。


・x, y = g(L)
 関数g(L)で上記リストLの近似分数の分子と分母を取得し、x,yとします。


・s = str(x)
 L = [int(s[i]) for i in xrange(len(s))]
 上記で求めた分子xを文字型にして、1桁ずつ数値にしたリストをMとします。


・return sum(M), x, y
 分子の各桁の合計値、分子、分母を返します。


5.関数の外
・n = 100
 print f(n)
 問題に合うように100を関数fに渡し戻り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
272,
6963524437876961749120273824619538346438023188214475670667,
2561737478789858711161539537921323010415623148113041714756

(追記)
・20130127 ネタバレ注意の追記など

2012年11月11日日曜日

プロジェクトオイラー 問64

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。

問64「N ≤ 10000 について奇数の周期をもつ連分数はいくつあるか?」
平方根は連分数の形で表したときに周期的であり, 以下の形で書ける:
           1
N = a0 + --------------------
                1
         a1 + ---------------
                     1
              a2 + ----------
                          1
                   a3 + -----
                         ...

例えば, √23を考えよう.
                             1              1
√23 = 4 + 23 - 4 = 4 + ------------ = 4 + ---------------
                             1                 √23 - 3
                          --------         1 + --------
                          √23 - 4                  7
となる.
この操作を続けていくと,
              1
√23 =  4 + --------------------------
                   1
             1 + ---------------------
                        1
                  3 + ----------------
                             1
                       1 + -----------
                                   1
                             8 + -----
                                  ...
を得る.
操作を纏めると以下になる:
a0 = 4, 1/(√23-4) =  (√23+4)/7  = 1 + (√23-3)/7
a1 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
a2 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
a3 = 1, 7/(√23-4) = 7(√23+4)/7  = 8 + (√23-4)
a4 = 8, 1/(√23-4) =  (√23+4)/7  = 1 + (√23-3)/7
a5 = 1, 7/(√23-3) = 7(√23+3)/14 = 3 + (√23-3)/2
a6 = 3, 2/(√23-3) = 2(√23+3)/14 = 1 + (√23-4)/7
a7 = 1, 7/(√23-4) = 7(√23+4)/7  = 8 + (√23-4)

よって, この操作は繰り返しになることが分かる.
表記を簡潔にするために,
√23 = [4;(1,3,1,8)] と表す.
(1,3,1,8)のブロックは無限に繰り返される項を表している.

最初の10個の無理数である平方根を連分数で表すと以下になる.
√2=[1;(2)], period=1
√3=[1;(1,2)], period=2
√5=[2;(4)], period=1
√6=[2;(2,4)], period=2
√7=[2;(1,1,1,4)], period=4
√8=[2;(1,4)], period=2
√10=[3;(6)], period=1
√11=[3;(3,6)], period=2
√12= [3;(2,6)], period=2
√13=[3;(1,1,1,1,6)], period=5
N ≦ 13で奇数の周期をもつ平方根は丁度4つある.

N ≦ 10000 について奇数の周期をもつ平方根が何個あるか答えよ.
-----



注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。






私の解答例は以下です。畳んでいます。
def g(n, b, c, L, M):
	if (b, c) in M: return L
	if n==b*b: return L
	M.append((b, c))
	C, t = (n-b*b)/c, int(n**(1.0/2))
	u = (b+t)%C
	L.append((b+t)//C)
	return g(n, t-u, C, L, M)

def f(n):
	t = 0
	for i in xrange(n+1):
		s = g(i, int(i**(1.0/2)), 1, [], [])
		if len(s)%2: t += 1
	return t

n = 10000
print f(n)


以下のように考えました。
  c                     √n-B
------ ・・・(A),   A + ----- ・・・(B)
√n-b                     C
として、(A)を連分数(B)にすることを考えます。
なお、(B)で計算したAの値を循環部分として保存し、
B, Cの値を次世代のb, cとして再帰的に計算します。
まず(A)を有理化します。
  c          c*(√n+b)     c*(√n+b)
------ = --------------- = ----------- ・・・(A)'
√n-b     (√n-b)*(√n+b)    (n-b*b)
連分数展開後の(B)の形に合わせるようにて分母Cを決めます。
    (n-b*b)
C = ------- ・・・(C)
       c
(A)'に(C)を当てはめ、
(√n+b)
------- ・・・(D)
   C
(D)から、連分数(B)のAに相当する部分をくくりだします。
そのために、√nの整数部分tとして、
t = int(√n) ・・・(E)
(b+t)   √n-t
----- + ------ ・・・(F)
  C       C
ですが、連分数(B)のAは整数のため、
Aは(F)の第1項の整数部分で、余りuは第2項に含めます。
A = (b+t)//C ・・・(G)
u = (b+t)%C  ・・・(H)
(G)(H)を(F)に反映して
 (b+t)        √n-t +u
{----- - u} + ---------
   C               C
             √n - (t-u)
= (b+t)//c + -----------
                 C
以上より、A,B,Cをまとめておくと以下の通りです。
(E)(H)のt,uを使って、
A = (b+t)//C ・・・(G)
B = t-u      ・・・(I)
C = (n-b*b)/c・・・(C)

2つの関数を使っています。

1.関数g(n, b, c, L, M)
  c                        √n-B
------ を連分数展開して、A + ------ (分数部分の絶対値は1未満)
√n-b                         C
に式変形したときの、循環部分のAのリストLを取得します。
リストMにはbとcの組、タプル(b,c)を貯めます。
例)√23 = [4;(1,3,1,8)]の場合、[1,3,1,8]を返します。
・関数内に自分自身を呼び出す記述のある再帰関数です。

・if (b, c) in M: return L
 まず再帰の終了条件です。
 引数b,cの組が今までにでてきていれば、循環部分がひと回りしたことになるので処理を終了し循環部分のAのリストLを呼び出し元に戻します。

・if n==b*b: return L
 この条件の場合、入力値の分母が0でゼロ割り算エラーのため連分数にならないのでこの時点でのリストLを返却します。

・M.append((b, c))
 ここまでくれば、引数b,cは連分数の循環しない部分なので、タプル(b,c)をリストMに貯めます。

・C, t = (n-b*b)/c, int(n**(1.0/2))
 上記説明の(C)(E)を実装しました。
 **演算子は累乗で、平方根を求める部分を1/2乗として記述しました。
 整数で1/2とするとint型のため値が0になるので、1.0とfloat型にすることで
 1/2がfloat型の0.5になります。分母の2を2.0にしても同じことです。
 共通モジュールのmathのsqrt関数を使っても同じことです。
 import math
 t = math.sqrt(n)

・u = (b+t)%C
 上記説明の(H)を実装しました。

・L.append((b+t)//C)
 循環部分リストLには、上記説明の(G)を実装しました。

・return g(n, t-u, C, L, M)
 連分数展開を続けます。
 現世代のB,Cが次世代のb,cなので、上記説明の(C)(I)を実装します。


2.関数f(n)
・問題に合う、n以下で奇数の周期を持つ平方根の個数を返します。

・t = 0
 個数カウンタtを初期クリアしておきます

・for i in xrange(n+1):
 ループ変数iとして、0以上n+1未満(n以下)の整数を発生させます。

・s = g(i, int(i**(1.0/2)), 1, [], [])
 関数gで連分数展開したときの循環部分リストを取得しsとします。
 第1引数は連分数展開する平方根の元の値なので、ここではループ変数iです。
 第2引数は√iからくくりだせる最大整数なので、ここでは√iの整数部分です。
 第3引数には1、第4第5引数には空のリストを設定します。

・if len(s)%2: t += 1
 %演算子は割り算の余りです。
 len(s)は関数gで取得した連分数循環部分リストsの要素数です。
 この要素数が奇数ということは2で割った余りが1で、1は論理的にTrueなので、上記if文が成り立ちます。このとき、個数カウンタtを1つカウントアップします。

・return t
 変数iのループが終了したら、個数カウンタを呼び出し元に返します。

3.関数の外
・n = 10000
 print f(n)
 問題に合うように10000を関数fに渡し戻り値を表示します。
 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
1322


(追記)
・20121208 式(C)の分子分母が逆だったので修正。説明文の中に2箇所あり。
 誤)C = c/(n-b*b)・・・(C)
 ↓
 正)C = (n-b*b)/c・・・(C)
・20130127 ネタバレ注意の追記など

2012年11月2日金曜日

プロジェクトオイラー 問63

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。

問63 「n 乗して得られる n 桁の正整数は何個あるか?」
5桁の数 16807 = 75は自然数を5乗した数である.
同様に9桁の数 134217728 = 89も自然数を9乗した数である.

自然数をn乗して得られるn桁の正整数は何個あるか?

-----

注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。






私の解答例は以下です。畳んでいます。
def f():
	import math
	return [(a,b,a**b) for a in xrange(1,10) for b in xrange(1, int(1/(1-math.log(a, 10)))+1)]

L = f()
print len(L),L


自然数aをb乗してb桁の正整数ということは、桁数の区切りは10の累乗値なので以下の式が成り立ちます。
0 ≦ 10(b-1) ≦ ab < 10b ・・・(A)
a>1、b>1 ・・・(B)

まず、aの変域を求めます。(A)の右部分の2式より、
1 ≦ a < 10 ・・・(C)

次にbの変域を求めます。
(A)の左部分をbについて解きます。
まず、左部分の2式の対数をとって、
(b-1)*log10 ≦ b*(log a)

次にbを左辺に寄せて、
(b-1)/b ≦ (log a)/log10

対数の底を10にして、
1-(1/b) ≦ log10a

1-log10a ≦ 1/b

b ≦ 1/(1-log10a)

(B)よりb>0なので、bの変域は
1 ≦ b ≦ 1/(1-log10a) ・・・(D)

上記(C)(D)の範囲で二重ループさせます。

2つの関数を使っています。

1.関数f()
・import math
 数学用関数モジュールmathをインポートして使えるようにしておきます。

・return [(a,b,a**b) for a in xrange(1,10) for b in xrange(1, int(1/(1-math.log(a, 10)))+1)]
 内包表記内に二重ループを入れています。
 まず前の「for a in xrange(1,10)」で、aが1以上10未満の整数でループします。
 そしてループ変数aを後ろのfor文へ渡します。
 後ろの「for b in xrange(1, int(1/(1-math.log(a, 10)))+1)」で、
 bが1以上1/(1-math.log10a)未満の整数でループします。
 math.log()関数は対数関数です。第2引数は対数の底です。
 int関数で切り捨てて整数にします。
 このa,bの値を先頭の式、(a,b,a**b)に渡します。
 これは問題に合うa, b, aのb乗の値のタプルです。
 このタプルを内包表記でリストにため、そのまま呼び出し元に返します。

2.関数の外
・L = f()
 print len(L),L
 関数f()で問題に合うa, b, aのb乗のタプルを格納したリストを取得しLとします。
 求める個数、とそのリストを表示します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
49

(追記)
・20130127 ネタバレ注意の追記など

2012年10月31日水曜日

プロジェクトオイラー 問62

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。

問62 「立方数になるような桁の置換をちょうど5つもつ最小の立方数を求めよ」
立方数 41063625 (3453) は, 桁の順番を入れ替えると2つの立方数になる:
56623104 (3843) と 66430125 (4053) である.
41063625は, 立方数になるような桁の置換をちょうど3つもつ最小の立方数である.
立方数になるような桁の置換をちょうど5つもつ最小の立方数を求めよ.


-----





注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。





私の解答例は以下です。畳んでいます。

def f(n):
	j = 0
	while True:
		j += 1
		L = [i*i*i for i in xrange(10**j)]
		M = [[str(s).count(str(t)) for t in xrange(10)] for s in L]
		for s in M:
			t = M.count(s)
			if t>=n:
				P = [i for i, u in enumerate(M) if s==u]
				N = [L[i] for i in P]
				return N[0], N, P

n = 5
print f(n)

最初は立方数の1つひとつの順列を作り、それが立方数かチェックしようとしましたが、1分では処理が終わりません。
そこで、立方数の1つひとつを構成する数字別個数リストを作って、同じリストになるものの数が指定個数になるかチェックするようにしました。

関数を1つ使っています。

1.関数f(n)
・j = 0
 変数jは10の累乗ごとに区切ってチェックするのでその累乗値カウンタです。初期クリアします。

・while True:
 無限ループです。

・j += 1
 累乗カウントをカウントアップします。

・L = [i*i*i for i in xrange(10**j)]
 Lは立方数リストです。10**jは10のj乗で、0からこの値までを1区切りにして処理を続けます。

・M = [[str(s).count(str(t)) for t in xrange(10)] for s in L]
 Mは数字別個数リストのリストです。二重の内包表記です。
 まず外側の内包表記の「for s in L」で先ほどの立方数リストLから1つずつ取り出しfor文の前の内側の内包表記に渡します。
 次に内側の内包表記の「for t in xrange(10)」で、ループ変数tに0以上10未満の整数が1つずつ設定されfor文の前に渡されます。
 「str(s).count(str(t))」で、立方数sに含まれる0以上10未満のいずれかのtの個数を数えます。
 例えば、s=1006012008の場合、内側の内包表記は[5,2,1,0,0,0,1,0,1,0]です。
 このような内側の内包表記が数字別個数リストで、外側の内包表記でそれをさらにリストにためています。

・for s in M:
 Mは数字別個数リストのリストから1つずつ取り出した数字別個数リストをループ変数sに設定します。

・t = M.count(s)
 リストMにある、ループ変数sの数字別個数リストと同じものの件数を数えtとします。

・if t>=n:
 上記tが指定個数n以上ならば、問題に合います。

・P = [i for i, u in enumerate(M) if s==u]
 Pは問題にあう立方数の元の数のリストです。
 enumerate関数で0からの連番を発生させます。このインデックス値が立方数の元の数です。

・N = [L[i] for i in P]
 リストPの値をインデックス値に立方数リストLから値を集め、リストNとします。

・return N[0], N, P
 問題に合う最小の立方数、桁を置換した立方数のリスト、その元の数のリストの3つを呼び出し元に返します。
 
2.関数の外
・n = 5
 print f(n)
 問題に合うように指定個数を5個として関数fで関連情報リストも含めて取得し、表示します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。(127035954683L,
[127035954683L, 352045367981L, 373559126408L, 569310543872L, 589323567104L],
[5027, 7061, 7202, 8288, 8384])

(追記)
・20130127 ネタバレ注意の追記など

2012年10月30日火曜日

プロジェクトオイラー 問61

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。

問61 「6つの巡回する4桁の図形数となる唯一の集合の和を求めよ」
三角数, 四角数, 五角数, 六角数, 七角数, 八角数は多角数であり, それぞれ以下の式で生成される.

三角数P3,n=n(n+1)/2 1, 3, 6, 10, 15, ...
四角数P4,n=n2 1, 4, 9, 16, 25, ...
五角数P5,n=n(3n-1)/2 1, 5, 12, 22, 35, ...
六角数P6,n=n(2n-1) 1, 6, 15, 28, 45, ...
七角数P7,n=n(5n-3)/2 1, 7, 18, 34, 55, ...
八角数P8,n=n(3n-2) 1, 8, 21, 40, 65, ...

3つの4桁の数の順番付きの集合 (8128, 2882, 8281) は以下の面白い性質を持つ.
  1. この集合は巡回的である. 最後の数も含めて, 各数の後半2桁は次の数の前半2桁と一致する
  2. それぞれ多角数である: 三角数 (P3,127=8128), 四角数 (P4,91=8281), 五角数 (P5,44=2882) がそれぞれ別の数字で集合に含まれている
  3. 4桁の数の組で上の2つの性質をもつはこの組だけである.
同じように, 6つの4桁の数からなる順番付きの集合で,
  1. 集合は巡回的である
  2. 三角数, 四角数, 五角数, 六角数, 七角数, 八角数が全て表れる唯一の集合を求め, その和を求めよ.

-----

注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。






私の解答例は以下です。たたんでいます。
def p3(n): return n*(n+1)//2
def p4(n): return n*n
def p5(n): return n*(3*n-1)//2
def p6(n): return n*(2*n-1)
def p7(n): return n*(5*n-3)//2
def p8(n): return n*(3*n-2)

def h(f, a, b):
	m = __import__(__name__)
	L = []
	for i in xrange(b):
		s = getattr(m, f)(i)
		if a<=s<b: L.append(s)
		elif b<=s: break
	return L

def g(L, M=[], n=2):
	if len(L)<=len(M): return M
	for s in L[len(M)]:
		if M and str(M[-1])[-n:]!=str(s)[:n]: continue
		if len(L)-len(M)>1:
			M = g(L, M[:]+[s], n)
			if len(L)<=len(M): break
			if s==L[len(M)][-1]: break
		elif str(s)[-n:]==str(M[0])[:n]: return M+[s]	
	else:
		M.pop()
	return M

def f(a, b, n=4, m=2):
	import itertools
	P = [h("p"+str(i), 10**(n-1), 10**n) for i in xrange(a, b+1)]
	for t in itertools.permutations(xrange(b+1-a)):
		L = [P[int(i)] for i in t]
		M = g(L, n=m)
		if M: break
	return sum(M), M, [str(i+a)+"kakusu" for i in t]

if __name__=="__main__":
	s, M, t = f(a=3, b=8, n=4, m=2)
	print s, M
	print t

角数の関数6個と他3個の関数を使っています。

1.関数p3(n)~p8(n)
・三角数から八角数の定義式のとおりです。

2.関数h(f, a, b)
・当pythonソースファイル内の関数fにb未満の値を渡し、a以上b未満の戻り値のリストを返します。

・m = __import__(__name__)
 __import__関数は、指定したモジュールを実行時に動的にインポートする関数です。
 _name__変数には、このプログラムのpythonソースのモジュール名(ファイル名の拡張子よりも前の部分)が設定されています。
 通常のimport関数はインポートするモジュール名を事前に固定文字列で指定しておきますが、自分のファイル自身をインポートする場合はファイル名変更しても影響を受けないように、動的にインポートします。
 後で関数を動的にインポートするための準備です。
 変数mには動的にインポートしたモジュールオブジェクトが設定されます。
 
・for i in xrange(b):
 ループ変数iは0以上b未満の範囲で1ずつ増加します。

・s = getattr(m, f)(i)
 getattr(m, f)関数は、モジュールmの関数fを呼び出します。
 変数iは関数fの引数です。
 当pythonソースファイル中の関数fに先ほどのループ変数iを渡し、その戻り値を変数sとします。

・if a<=s<b: L.append(s)
 elif b<=s: break
 上記の値sが、a以上b未満ならリストLに追加し、b以上なら終了します。

3.関数g(L, M=[], n=2)
・リストLは三角数などの角数リストが格納されたリストです。
 リストMはリストLの要素である角数リストを順にn桁で巡回的連結した値を格納したのリストです。
 リストLの全部の要素リストをn桁で巡回的連結した場合、その連結した値のリストを返します。

・if len(L)<=len(M): return M
 自分自身を呼び出す再帰関数なので、まず再帰呼び出し終了条件です。
 リストLの要素リストの個数以上に、順に巡回的連結した値の個数があれば順に巡回的に連結できた値のリストMを返し処理終了です。

・for s in L[len(M)]:
 len(M)は巡回的連結済みの個数であり、リストLのインデックス値としては次の巡回的連結要素を探すリストのインデックス値です。
 リストLから次に巡回的に連結する値を探す対象の角数リストを選び、ループ変数sとします。

・if M and str(M[-1])[-n:]!=str(s)[:n]: continue
 リストMに要素があり、その最後の要素がループ変数sにn桁で巡回的連結しなければsは対象外なので、次のループ変数sにループを進めます。
 リストMの最後の要素とループ変数sがn桁で巡回的連結している場合は以下の処理を続けます。

・if len(L)-len(M)>1:
  M = g(L, M[:]+[s], n)
  if len(L)<=len(M): break
 巡回的連結値の個数とリストLの要素数の差が1より大きい、つまり、リストLの最後までいっていない場合、さらに自分自身である関数g()を再帰呼び出しして、巡回的連結可能な値を探します。
 この再帰呼び出しのときの引数はリストLと巡回判定桁数nはそのままの値で、巡回的連結値リストMにはsを追加した状態で次を探します。
 再帰呼び出しから戻ってきたら、最初の行と同じ終了判定をします。

・elif str(s)[-n:]==str(M[0])[:n]: return M+[s]
 巡回的連結値の個数よりも、リストLの要素数が多くないということは、これらの個数が等しく、リストLの要素である角数リストすべてから巡回的連結値を探し出した状態なので、最後のチェックとして先頭要素との巡回的連結で連結の輪が閉じれるかを判定します。まだリストMに追加前のループ変数sからリストMの最後の要素にn桁で巡回的連結できればリストMが完成しこれを呼び出し元に返します。
 M+[s]は、リストMと要素sだけのリストの連結です。appendして返すことと同じです。

・else:
 M.pop()
 このelseはループ変数sのfor文のelseです。
 for文が途中でbreakせず、ループ変数の最後の要素まで処理し終った後に行う処理です。
 ここでは呼び出し元に返すこともなくfor文が終わったということなので、現在にリストMの最後の要素が巡回的連結先が無いことが確定しているので、このリストMの最後の要素をpop()関数で取り除きます。
 
・return M
 最後まで処理が来た場合は、リストMをそのまま返します。

4.f(a, b, n=4, m=2)
・n桁のa角数からb角数を1つずつ選びm桁で巡回的につながる組を返します。

・import itertools
 標準モジュールitertoolsを使えるようにしておきます。

・P = [h("p"+str(i), 10**(n-1), 10**n) for i in xrange(a, b+1)]
 リストPは、関数h()による角数リストをためたリストです。内包表記で書いています。
 for文でループ変数iにa以上b以下の整数を1つずつ設定します。
 ループ変数iは前に送られ、関数h()を呼び出しその戻り値を設定します。
 関数hの第1引数で、"p"+str(i)、つまり三角数などの関数、「p3」などを指定します。第2引数と第3引数で戻り値がn桁になるようします。例えばn=3の場合、100、1000を引き渡します。
 
・for t in itertools.permutations(xrange(b+1-a)):
 ループ変数tは、後ででてくるリストLのindex番号の組合せに使います。
 permutations()は順列を返す関数です。
 例えばrange(3)ならば、012, 021, 102, 120, 201, 210と順に返します。
 xrange(b+1-a)で、角数リストの個数分だけ0からの連番を発生させます。
 
・L = [P[int(i)] for i in t]
 上記tで先ほどのリストPの中を組み替えて、6つの角数リストを任意の順に組み替えます。
・M = g(L, n=m)
 if M: break
 関数g()を呼び出し、空のリスト(論理的にFalse)ならば、ループを続け次の
 順列の角数リストで巡回的連結可能な値を探しにいきます。

・return sum(M), M, [str(i+a)+"kakusu" for i in t]
 巡回的連結がした値の合計値、その値リスト、およびどの角数リストからとってきたかを呼び出し元に返します。

5.関数の外
・s, M, t = f(a=3, b=8, n=4, m=2)
 print s, M
 print t
 三角数から八角数まで各角数のうち4桁の値を1つずつ2桁で循環的連結できる値の和、そのときの値のリスト、どの角数リストからとってきたかのリストの3つを受け取り表示します。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
28684 [8256, 5625, 2512, 1281, 8128, 2882]
['3kakusu', '4kakusu', '7kakusu', '8kakusu', '6kakusu', '5kakusu']


(追記)
・20130127 ネタバレ注意の追記など

2012年10月9日火曜日

プロジェクトオイラー 問60

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。

問60 「任意の2つの素数を繋げたときに別の素数が生成される, 5つの素数の組を求めよ」
素数3, 7, 109, 673は非凡な性質を持っている.
任意の2つの素数を任意の順で繋げると, また素数になっている.
例えば, 7と109を用いると, 7109と1097の両方が素数である.
これら4つの素数の和は792である.
れは, このような性質をもつ4つの素数の組の和の中で最小である.

任意の2つの素数を繋げたときに別の素数が生成される, 5つの素数の組の和の中で最小のものを求めよ.


-----

注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。





最初にコーディングしたときにfor文の五重ループになってしまったので、どうしても再帰関数に改良したかったことと、処理速度の向上とで、1ヶ月近くかかってしまいました。
上限値を設けてやっと1分ルールを満たしているので、まだまだロジックが甘いのですが、行き詰まってしまったこともあり、現在のプログラムのまま提示しておきます。

答えは合っていますが、問題では「最小のもの」を求めるところ、題意を満たす「最初のもの」を求めるプログラムです。

私の解答例は以下です。畳んでいます。
def h(n):
	if n==2: return True
	if n<2 or (not n%2): return False
	for i in xrange(3, n, 2):
		if i*i>n: break
		if not n%i: return False
	return True

def e(n):
	i = n+1+int(n%2)
	while not h(i): i += 2
	return i

def g(L, t):
	t = str(t)
	for s in L:
		s = str(s)
		if not h(int(s+t)): return False
		if not h(int(t+s)): return False
	return True

def f(n, L=[], p=0):
	if len(L)>=n: return L
	while p<10**(n-1):
		p = e(p)
		if g(L, p):
			L = f(n, L+[p], p)
			if len(L)>=n: break
	else:
		L.pop()
	return L

n = 5
s = f(n)
print sum(s), s


4つの関数を使っています。

1.関数h(n)
・素数判定です。問58で作成したものを流用しました。
 素数ならTrue、素数でなければFalseを返します。

2.関数e(n)
・引数を超える最小の素数を返します。n>1のとき有効です。

・i = n+1+int(n%2)
 変数iが素数候補です。
 %演算子は割り算の余りです。割り切れれば0、割り切れないと1なので、
 式全体では偶数なら1足し、奇数なら2足すことになります。
 後続処理の準備として、素数候補iを奇数にしておきます。

・while not h(i): i += 2
 関数h()で素数判定をして、素数でなければ2ずつ増やしてチェックします。
 前準備でこの処理にくるときには素数候補iは奇数にしているので、
 iが偶数で無限ループになることはありません。

・return i
 ここまでの処理でnを超える最初の素数としてiを返却します。

3.関数g(L, t)
・引数リストLの各要素と引数tを任意の順に2つ連結した値が、すべて素数かどうかを判定します。

・t = str(t)
 後で連結するので、引数tを文字型にします。

・for s in L:
 引数リストLから順にループ変数sに設定します。

・s = str(s)
 後で連結するので、ループ変数sを文字型にします。

・if not h(int(s+t)): return False
 if not h(int(t+s)): return False
 上記の文字列sと文字列tを前後とその入れ替えの2パターンで連結し、
 int関数で数値化して、関数hで素数判定します。
 素数でなければ即終了でFalseを返します。

・return True
 引数リストLと引数tでここまでで処理を抜けていなければ、そのすべての組合せが素数なので、Trueを返します。

4.関数f(n, L=[], p=0)
・問題に合う値のリストを返します。
 自分で自分を呼び出す記述のある再帰関数です。
 引数nは求める素数の個数です。
 引数Lは素数を格納するリストで初期値は空リストです。
 引数pは最後にチェックした素数で、初期値は0とします。

・if len(L)>=n: return L
 再帰関数の終了条件です。
 リストLに求める個数n個以上の要素がたまったら処理終了で、リストLを返します。

・while p<10**(n-1):
 対象とする素数が青天井では処理が終わらないので、とりあえず求める素数の個数分の桁数、10**(n-1)、未満の素数で、処理することにしました。
 求める素数が2個なら10以内の素数の組で最小組[3,7]が見つかり、
 同じく3個なら100以内の素数の組で最小組[3, 37, 67]が見つかり、
 同じく4個なら問題文にあるように1000以内の素数の組が見つかります。

・p = e(p)
 得られている最後の素数pを関数eで処理し、次に大きい素数を改めてpとします。

・if g(L, p):
 関数gで引数リストLの各要素と引数pを任意の順に2つ連結した値が、すべて素数かどうかを判定します。Falseの場合はwhileループの次の値へ進みます。

・L = f(n, L+[p], p)
 1つ前の関数g()ですべて素数の場合、次に大きい素数pを候補として、さらに題意にある素数を探して追加するため、当関数f()を再帰呼び出しします。
 なお、第2引数の「L+[p]」は、再帰呼び出しで引数Lで受けるので、お試しでリストLの末尾に追加することになります。

・if len(L)>=n: break
 再帰呼び出しから戻り、リストLの要素数がn以上ならば探索終了でwhileループから抜けます。

・else:
 L.pop()
 このelse節はインデント(行頭のずらし位置)から明らかなように、for文のelseです。直前のif文のelse節ではありません。
 for文のelse節の処理タイミングは、forループの最後の値まで処理が行われた場合のループ終了直後です。
 ここまでの処理でbreakで抜けない場合、再帰呼び出しの際にお試しで追加したpが求める値ではないということになります。
 このため、関数pop()で末尾の要素を削除しておきます。

・return L
 処理の最後にリストLを呼び出し元に戻します。

5.関数の外
・n = 5
 s = f(n)
 print sum(s), s
 問題に合うように引数の値を5として関数f()を呼び出します。
 返されたリストをsとして、そのリスト要素合計、及び参考としてリスト要素を表示します。


 
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
26033 [13, 5197, 5701, 6733, 8389]


(追記)
・20130127 ネタバレ注意の追記など

2012年9月12日水曜日

プロジェクトオイラー 問59

プロジェクトオイラーの問題をpythonでやってみます。
日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索してください。

問59「XOR暗号化された暗号文を解読できるか?」
各文字はそれぞれ一意のコードに割り当てられている.よく使われる標準としてASCII (American Standard Code for Information Interchange) がある.ASCIIでは, 大文字A = 65, アスタリスク (*) = 42, 小文字k = 107というふうに割り当てられている.

モダンな暗号化の方法として, テキストファイルの各バイトをASCIIに変換し, 秘密鍵から計算された値とXORを取るという手法がある.
XOR関数の良い点は, 暗号化に用いたのと同じ暗号化鍵でXORを取ると平文を復号できる点である.
65 XOR 42 = 107であり, 107 XOR 42 = 65である.

破られない暗号化のためには, 鍵は平文と同じ長さのランダムなバイト列でなければならない.ユーザーは暗号文と暗号化鍵を別々の場所に保存する必要がある.また, もし一方が失われると, 暗号文を復号することは不可能になる.
悲しいかな, この手法はほとんどのユーザーにとって非現実的である.

そこで, 鍵の変わりにパスワードを用いる手法が用いられる.パスワードが平文より短ければ (よくあることだが), パスワードは鍵として繰り返し用いられる.
この手法では, 安全性を保つために十分長いパスワードを用いる必要があるが, 記憶するためにはある程度短くないといけない.

この問題での課題は簡単になっている.
暗号化鍵は3文字の小文字である.
cipher1.txtは暗号化されたASCIIのコードを含んでいる.
また, 平文はよく用いられる英単語を含んでいる.

この暗号文を復号し, 平文のASCIIでの値の和を求めよ.

-----


注意!!!
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。






私の解答例は以下です。畳んでいます。
def g(s, n):
	import itertools
	for t in itertools.product(s, repeat=n): yield "".join(t)

def f(sIn, a, n):
	L = [int(t) for t in open(sIn).read().split(",")]
	for s in g(a, n):
		M = ([ord(t) for t in s]*len(L))[:len(L)]
		N = [i^j for i,j in zip(L, M)]
		R = [chr(i) for i in N]
		t = "".join(R)
		if " the " in t: return sum(N), s, t

import os
sIn = os.path.join(os.path.dirname(__file__), "cipher1.txt")
a = "abcdefghijklmnopqrstuvwxyz"
n = 3
print f(sIn, a, n)


2つの関数を使っています。

1.関数g(s, n)
・文字列sからn個選んだ直積を連結した文字列を返します。
 直積とは、重複を許し順番が違えば別のものとして扱う組み合わせです。
 例えば、abcから2つ選んだ直積は、aa,ab,ac,ba,bb,bc,ca,cb,ccです。

・import itertools
 標準モジュールitertoolsを使えるようにしておきます。

・for t in itertools.product(s, repeat=n): yield "".join(t)
 標準モジュールitertoolsにあるproduct関数で、文字列sからをn個選んだ直積の構成要素のタプルを生成し、for文のループ変数tに設定します。
 forループの中で、タプルtを区切り文字無しで連結し、yield文で返します。
 yield文で返すと、返した時点で処理内容が一旦凍結され、再度呼び出されたときに処理を再開して続きの処理による値を返します。
 このような動きをするyield文で返す関数をジェネレータ関数といいます。
 なお、単発呼び出しで終了するreturn文で返す関数はイテレータ関数といいます。
 また、product関数はpython2.6以上です。
 python2.5にもitertoolsモジュールはありますが、product関数が含まれていません。

2.関数f(sIn, a, n)
・主処理です。
・入力ファイルsInの中の文字列を、指定した文字列aからn個選んだ秘密鍵sによって順にxor演算で復号した数値列N、及びこれをasciiコード扱いして得られる復号文字列tがあるとき、数値列Nの合計値、キー文字列s、復号文字列tを返します。

・引数sInは入力ファイルのフルパス、文字列aは秘密鍵の構成文字候補の文字列、nは秘密鍵の桁数です。

・L = [int(t) for t in open(sIn).read().split(",")]
 暗号文の各数値をもつリストLを作成します。
 入力ファイルsInを開き、その内容をすべて読み込んで、カンマで区切った文字列を
 ループ変数tに設定します。
 このtをfor文の前に送り、int関数で数値化してリストLとします。
 なお、開いたファイルオブジェクトはpythonがタイミングをみて閉じくれます。

・for s in g(a, n):
 関数gで文字列aからn個選んだ秘密鍵候補をループ変数sとします。

・M = ([ord(t) for t in s]*len(L))[:len(L)]
 Mは、秘密鍵候補sの構成要素の数値を暗号文の長さまで繰り返した復号用鍵のタプルです。
 内包表記のループ変数tとして秘密鍵候補sの文字1つずつを、for文の前に送り、ord関数でasciiコードとして対応する数値にしてリスト化します。
 このリストをlen(L)倍、つまり暗号文の長さ倍して、秘密鍵を順に繰り返し連結した状態にしてから、[:len(L)]で暗号文の長さの分だけ取り出しておきます。
 
・N = [i^j for i,j in zip(L, M)]
 Nは復号化した数値のリストです。
 zip関数で上記のリストLとリストMから1つずつ組み合わせて取り出し、ループ変数i,jとし、内包表記のfor文の前に送り、^演算子で数値j,jのxorを計算しリストにため、リストNとします。

・R = [chr(i) for i in N]
 復号数値リストNから1つずと取り出しchr関数でascii文字化してリストにため、リストRとします。
 
・t = "".join(R)
 上記の文字列リストRの要素を区切り文字無しで連結し、復号文tとします。

・if " the " in t: return sum(N), s, t
 復号が正しいか判定します。
 判定は「よく用いられる英単語」があるという問題文に従い、半角空白ではさんだ「 the 」としてみました。結果的にこれで判定できました。
「 the 」があれば、復号化した数値リストNの要素の合計、秘密鍵、復号文tを返します。

3.関数の外
・import os
 ファイルを扱うので、標準モジュールosを使えるように取り込んでおきます。

・sIn = os.path.join(os.path.dirname(__file__), "cipher1.txt")
 暗号文のフルパスです。
 「__file__」は実行するpythonファイルのフルパスです。
 os.path.dirname関数で親フォルダまで取得し、ファイル名とともに、
 os.path.joinでパス表示を連結しました。
 つまり、暗号文は実行するpythonファイルと同じフォルダに置いて実行するように記述しています。
 
・a = "abcdefghijklmnopqrstuvwxyz"
 秘密鍵の候補です。問題文から小文字のアルファベットです。

・n = 3
 秘密鍵の桁数です。問題文から3としました。

 

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
107359
秘密鍵はgod。平文の内容は、聖書ヨハネによる福音書の第1章第1節から第14節。

(追記)
・20130129 ネタバレ注意の追記など