2013年3月10日日曜日

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

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

問74「60個の循環しない項を持つ階乗列となる数を特定せよ」

145は各桁の階乗の和が145と自分自身に一致することで有名である.

1! + 4! + 5! = 1 + 24 + 120 = 145

169の性質はあまり知られていない. これは169に戻る数の中で最長の列を成す.
このように他の数を経て自分自身に戻るループは3つしか存在しない.

169 → 363601 → 1454 → 169
871 → 45361 → 871
872 → 45362 → 872

どのような数からスタートしてもループに入ることが示せる.
例を見てみよう.

69 → 363600 → 1454 → 169 → 363601 (→ 1454)
78 → 45360 → 871 → 45361 (→ 871)
540 → 145 (→ 145)

69から始めた場合, 列は5つの循環しない項を持つ.
また100万未満の数から始めた場合最長の循環しない項は60個であることが知られている.

100万未満の数から開始する列の中で, 60個の循環しない項を持つものはいくつあるか?

-----


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




def h(n):
	c = 1
	for i in xrange(2, n+1): c *= i
	return c

def g(n, d): return sum([d[s] for s in str(n)])

def f(n, m):
	d = {}
	for i in xrange(10): d[str(i)] = h(i)
	N = []
	for i in xrange(n):
		s, M = i, []
		while (s not in M) and (i<=s):
			M.append(s)
			s = g(s, d)
		t = len(M)
		if s<i: t += N[s]
		N.append(t)
	return sum([1 for i in N if i==m])

n = 1000000
print f(n, 60)

1.関数h(n)
・nの階乗値を返します。再帰呼び出しは使わず単純ループで掛け算します。
・再帰呼び出しにした場合、n≧1000で「maximum recursion depth exceede」エラーになるので、再帰呼び出しの深さの限界値をsys.setrecursionlimit()で設定する必要があります。

2.関数g(n, g)
・数字nの各桁の階乗値の和を返します。
・dは1桁数字の階乗値を持つ辞書です。但し、キーは文字型です。
・数字nを文字型にして一文字ずつ階乗値辞書から階乗値を取得し、最後にsum関数で合計します。

3.関数f(n, m)
・n未満の数字のうち、循環しない各桁階乗和をmもつものの個数を返します。
・辞書dとして、関数gで使用する階乗値辞書を作っておきます。
・リストNには、i番目の要素に数字iの循環しない各桁階乗和の個数を入れます。
・ループ変数iはn未満の数を昇順に1つずつとります。
・sは各桁階乗和の計算対象の候補です。初期値はループ変数iです。
・リストMには、各桁階乗和を1つずつ進めながら循環部分になるまで持ちます。
・whileループで、計算対象候補sがすでにリストMにあれば循環部分突入なので終了です。
 また計算対象候補sがループ変数i未満なら、その先の循環部分までの個数はリストNにあるのでやはりループは終了です。
・tで各桁階乗和をためたリストMの個数を求め、さらに計算対象候補sがループ変数i未満だった場合はそれ以降の循環部分に達するまでの個数を足します。
・上記の状態のtが数字iでの、循環しない各桁階乗和の個数です。リストNに追加しておきます。

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

2013年3月5日火曜日

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

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

問73「真既約分数をソートした集合内の1/3と1/2の間に何個の分数があるか?」

nとdを正の整数として, 分数 n/d を考えよう. 
n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≦ 8について既約分数を大きさ順に並べると, 以下を得る:
1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
1/3と1/2の間には3つの分数が存在することが分かる.

では, d ≦ 12,000 について真既約分数をソートした集合では,
1/3 と 1/2 の間に何個の分数があるか?

-----


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




最初に考えた方法は、
分子候補と分母候補の二重ループ内で、
ユークリッドの互除法でその分母子の最大公約数が1の分数を真既約分数として
カウントアップする方法です。
 しかしこの方法では私のPCでは1分半かかり1分ルールを守れませんでした。
が、自力ではここまでだったのでそのまま書いておきます。

def gcd(a, b):
	#if a<b: a, b = b, a
	while b: a, b = b, a%b
	return a
def f(n, x1, x2, y1, y2):
	c = 0
	for i in xrange(2, n+1):
		a, b = i*x1//x2+1, i*y1//y2
		if i*y1%y2>0: b += 1
		for j in xrange(a, b):
			if gcd(i, j)==1: c += 1
	return c

s = f(12000, 1, 3, 1, 2)
print s


1.gcd(a,b)
・「ユークリッドの互除法」を使用して最大公約数を返します。
 説明は問5を参照してください。
 関数fから呼び出されるとき、既約分数の分母、分子の順でに引数にセットすれば、引数は必ずa>bなので、割る数と割られる数を入れ替える必要はなくなるため、最初の行はコメントアウトしました。

2.f(n, x1, x2, y1, y2)
・分母がn以下の真既約分数のうち、x1/x2より大きくy1/y2より小さいものの個数を返します。
・ループ変数iが分母、jが分子です。
・分母iは2からnの整数です。
・分子jの開始値aと終了値bは、それぞれ分母iのx1/x2倍とy1/y2倍になります。
 //演算子は商の整数部を示します。
 開始値はx1/x2より大きいということで、i*x1//x2は含まないので、開始値aを+1しておきます。
 終了値はy1/y2より小さいということで、j=i*y1//y2のときに、
 この真既約分数全体で丁度y1/y2と同じとき、
 つまりはi*y1/y2で余りがあって割り切れない場合は対象に含めます。
 そこで余りの%演算子を使った、i*y1//y2が0より大きければ終了値bを+1します。
・gcd関数で分母子の最大公約数が1の場合、その分数は真既約分数です。


自分のPCでも1分ルールを守りたかったのですがやり方を思いつかなかったので、ネット検索で以下のことがわかりましたので参考にしました。

・分母の最大値nを与えて昇順に並べた真既約分数の数列を「ファレイ数列」という。
・ファレイ数列の隣接する分数では以下の定理が成り立つ。
(定理1)小さい方からa/b, c/dが隣接する場合、bc-ad = 1。
(定理2)小さい方からa/b, c/d, e/fが隣接する場合、c/d = (a+e)/(b+f) 
証明はネット検索してください。

定理(2)はファレイ数列で連続する3つの分数の真ん中の分数は、
左右の分数の、分子の和/分母の和という、隣接3要素間の関係を簡潔に示しています。
これを使えば、連続する2つの分数からその1つ先の隣接する分数が計算でき、
ファレイ数列を指定範囲の分だけ生成して数えることができます。

上記を利用した解答例は以下です。畳んでいます。
def f(n, x1, x2, y1, y2):
	a, b, cnt = x1, x2, 0
	c = (n*a+1)//b
	while (b*c-1)%a: c -= 1
	d = (b*c-1)//a
	while (c, d)!=(y1, y2):
		cnt+=1
		k = (n+b)//d
		a, b, c, d = c, d, k*c-a, k*d-b
	return cnt

s = f(12000, 1, 3, 1, 2)
print s

1.f(n, x1, x2, y1, y2)
・分母がn以下の真既約分数のうち、x1/x2より大きくy1/y2より小さいものの個数を返します。

・5行目までで、範囲の最小値x1/x2の左隣の分数c/dを求めます。
 このまま範囲全体のループにしても解答は求まりますが、cをカウントダウンしているwhileループのループ回数が大きく処理が遅いので、最初の値を求めることだけに留めます。

・6行目のwhileループで、隣接するa/bとc/dからその1つ大きい分数を求め、
 分母子の組ごと1つずつ大きい方に進めて範囲の最大値y1/y2になるまで、
 件数をカウントアップします。

・最初の1つの分数から隣接する大きい値の分数を求める方法は以下です。
 以下の連立方程式を解きます。
 ┌ bc-ad = 1(a,b,c,dは正の整数)(定理1)より。...(1)
 └ d≦n     (d, nは正の整数)...(2)

 (1)(2)より、
 d = (bc-1)/a ≦n、正の整数...(3)

 また、(1)より、bc ≦ an+1
 c ≦ (an+1)/b、正の整数...(4)
 分母子c, dはそれぞれの値が大きいほどその分数c/dはa/bとの差が小さくなります。
 なので、分子cはその最大値候補から1ずつ減らしたループの中で最初の整数条件を満たした値です。

 そこで、cは(4)より、(an+1)/bから1ずつ減らしたループで、
 停止条件は(3)より、(bc-1)/aが整数値かをその余りの有無で判定します。
 こうして分子cが決まれば、これを(3)に代入して分母dも決まります。
 これで、最初の分数、a/bとc/dが決まります。

・後は分母子c,dが、求める範囲の最大値y1/y2になるまで隣接分数を求めていきます。
 小さい方からa/b, c/d, e/fが隣接する場合、
 (定理2)より、c/d = (a+e)/(b+f)
 c/dが真既約分数なので、(a+e)/(b+f)の約分できる値をkとして、
 以下の連立方程式を解きます。
 ┌ kc = a+e
 └ kd = b+f
 
 ┌ e = kc-a ...(5)
 └ f = kd-b ≦ n(分母なので)...(6)
 (6)より、k ≦ (n+b)/d ...(7)
 
 これで、範囲の最小値x1/x2から1つずつ大きな分数を計算していき、
 範囲の最大値y1/y2になるまで個数をカウントアップします

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


2013年2月16日土曜日

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

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

問72「真既約分数の集合は何個の要素を持つか?」

nとdを正の整数として, 分数 n/d を考えよう. 
n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≦ 8について既約分数を大きさ順に並べると, 以下を得る:
1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
この集合は21個の要素をもつことが分かる.

d ≦ 1,000,000について, 真既約分数の集合は何個の要素を持つか.

-----


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





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

def f(n):
	L = range(n+1)
	for i in xrange(2, n):
		if L[i]<i: continue
		for j in xrange(i, n+1, i):
			L[j] /= i
			L[j] *= i-1
	return sum(L[2:])

n = 1000000
print f(n)


分母nの真既約分数の個数は「nと互いに素なn未満の数」の個数と同じです。
これは問69のトーティエント関数φ(n)の定義そのものです。
つまり求める真既約分数の個数は、φ(2)からφ(1000000)までの和です。

φ(n)の計算ですが、
「nと互いに素なn未満の数」は、言い換えれば、1からnまでのn個の数字の中から素因数で割り切れる数を外していき最後まで残った数字の個数です。

例えば
nの素因数に2があれば、1からnのうち、2,4,6,...のように偶数を外すので
個数はn*(1/2)個減り、残りはn*(1 - 1/2)=n*(1/2)になります。

nの素因数に3があれば、1からnのうち、3,6,9,...のよう3の倍数を外すので
個数はn*(1/3)個減り、残りはn*(1 - 1/3)=n*(2/3)になります。

このとき外す数からカウントすると、例えば6なら2の倍数と3の倍数で重複してカウントされてしまいます。そこで、残る数から考えると上記の累積値として求まります。

つまり、nの素因数が2,3,・・・,mの場合、素因数で割っていって残る数の個数φ(n)は、
 n*(1 - 1/2)*(1 - 1/3)*・・・*(1 - 1/m)
= n*(1/2)    *(2/3)    *・・・*(m-1)/m
となります。

そこで1からnまでの一つひとつの数字について、
素因数を求めて上記のようにφ(n)を計算し総和を求めてみましたが、
このロジックのまま実装すると1分ルールを守れませんでした。
素因数を何度も求める部分が重いようでした。

そこでさらにひと工夫です。
φ(n)をnの順に求めるのではなく素因数候補の方を順に反映していくようにしてみます。
つまり、mを素因数とするφ(i)すべてに対して、(m-1)/m倍を累積していきます。
素因数候補が2ならば、φ(2)から2刻みにφ(4), φ(6)・・・を途中計算として1/2倍します。
素因数候補が3ならば、φ(3)から3刻みにφ(6), φ(9)・・・を途中計算として2/3倍します。
素因数候補がmならば、φ(m)からm刻みにφ(m*2), ・・・を途中計算として(m-1)/m倍します。

1.関数f(n)
・φ(0)からφ(n)のリストLを返します。リストLのi番目の要素L[i]=φ(i)です。
・φ(n)の値のリストLとして、最初に、i番目の要素に値iをセットした値nまでのリスト[0, 1,..., n]を準備します。
・ループ変数iが素因数候補です。
 iを小さい値から順に処理したとき、
 L[i]==iならば、上記の(m-1)/mが反映されていないので素因数です。
 L[i]<iならば、すでにiより小さい値で(m-1)/mされてるので素因数ではありません。
・なお、(m-1)/m倍するとき、この分数を先に計算するとfloat型になります。
 int型のまま計算する方が処理時間が短くすむため、(m-1)倍と1/m倍は別々に計算しています。
・ループ変数jがφ(n)の引数nです。
・ループ処理がすべて終わったら、リストLの2番目以降の要素の和を求めます。


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

2013年1月30日水曜日

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

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

問71「真既約分数を昇順に並べ上げよ」
nとdを正の整数として, 分数 n/d を考えよう. 
n<d かつ HCF(n,d)=1 のとき, 真既約分数と呼ぶ.

d ≦ 8について既約分数を大きさ順に並べると, 以下を得る:
1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8
3/7のすぐ左の分数は2/5である.

d ≦ 1,000,000について真既約分数を大きさ順に並べたとき, 
3/7のすぐ左の分数の分子を求めよ.

-----


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





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

def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a

def f(n, a, b):
	L = [0,0,1]
	s = 1.0*a/b
	for i in xrange(1, n+1):
		j = int(i*s)
		if not j: continue
		while gcd(i, j)>1:
			j -= 1
		u = s - 1.0*j/i
		if 0<u<L[2]: L = [j, i, u]
	return L

n, a, b = 1000000, 3, 7
print f(n, a, b)


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

2.g(n,a,b)
 n以下の分母を持つ真規約分数を大きさ順に並べたときに、
 a/bのすぐ左の分数についての、分子、分母、その分数のa/bと差を返します。
 
 まず、基準値となるa/bを計算しsとします。
 1.0を掛けているのは整数型ではなく浮動小数点型で計算するためです。
 float関数で型変換するよりも速いです。
 ループ変数iは求める分数の分母で、jは分子です。
 分子は、分母と基準値を掛け、int関数で整数部をとることで算出します。
 jが0の場合、分数にならなので分母を1つ進めます。
 また、j/iが真既約分数にならない場合、真既約分数になるまでjを1つずつ小さい値にします。
 真既約分数の判定は、分母子の最大公約数が1であることとしました。
 こうして分数j/iが定まったところで、基準値との差分を求めuとします。
 uがここまで処理したうちの最小値ならば、分子j、分母i、差分uをリストLの値と差し替えてためます。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
[428570, 999997, 1.4285757138354782e-07]

2013年1月26日土曜日

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

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


問70「φ(n) が n の置換となる n を調べ上げよ」
オイラーのトーティエント関数 φ(n) (ファイ関数とも呼ばれる) とは, 
n 未満の正の整数で n と互いに素なものの個数を表す. 
例えば, 1, 2, 4, 5, 7, 8 は9未満で9と互いに素であるので, φ(9) = 6 となる. 
1 は全ての正の整数と互いに素であるとみなされる. よって φ(1) = 1 である.

面白いことに, φ(87109)=79180 であり, 87109は79180を置換したものとなっている.

1 < n < 10 7 で φ(n) が n を置換したものになっているもののうち, 
n/φ(n) が最小となる n を求めよ.

-----


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






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

def p(n): L = [0,0]+[1]*(n-1) i = 2 while i*i<=n: while not L[i]: i += 1 for j in xrange(i+i, n+1, i): L[j] = 0 i += 1 return [i for i in xrange(n+1) if L[i]] def f(n): N = [n, n, []] P = p(n//2) w = int(n**(0.5)) for i, u in enumerate(P): if w<u: break z = n//u for j in xrange(len(P)): if z<P[j]: break for v in P[j: i: -1]: s, t = u*v, (u-1)*(v-1) if sorted(str(s))==sorted(str(t)): a = 1.0*s/t if a<N[1]: N = [s, a, [u,v]] break return N n = 10000000 print f(n)


問69の発展問題ということで、同様にはなく、数学的に別の表現ができないかを考えます。
まず、問69とは逆に、n/φ(n)を最小にするためには、nをなるべく小さくしφ(n)はなるべく大きくします。
また、問69の検討を考慮すると、
・nよりもφ(n)が大きく効いてくるので、φ(n)を大きくする。
・φ(n)の最大値は、nが素数のときのφ(n)=n-1なので、なるべく大きな素数を持つ。

ここで、nとφ(n)が置換の関係との条件から、素数では問題に合わないことがわかります。
それは、nとn-1は必ず1文字は異なり、置換の関係にならないからです。

そこで素数2つの積で、φ(n)との置換の関係を満たすnということになります。
それは、nが複数の素数の積ならば、その因数である素数の倍数以外には割りきれないため、φ(n)が大きい値のままですみます。
また、この場合の因数として使用する素数の件数なるべく少ない方がよいということになります。

1.関数p(n)
・n以下の素数リストを返します。問37と同じです。

2.関数f(n)
 n未満の整数で、問題に合うn、n/φ(n)、nを構成する素数2つのリストを返します。
 候補の素数の範囲は、2つの素数の積がn未満ということから、n//2以下で十分です。
 素数候補のリストから、
 ループ変数uで小さい方から1つずつ取り出し、
 ループ変数vで大きい方から1つずつ取り出します。
 素数uと素数vの積とそのφ値をs,tとします。
 置換の関係のチェックはs,tを1文字ずつ取り出してそのソート状態が一致することとします。
 リストNに関係情報をためて、n/φ(n)の小さい値がきたら入替えます。

解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
[8319823, 1.0007090511248113, [2339, 3557]]

2013年1月20日日曜日

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

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



問69「n/φ(n) が最大となる n の値を求めよ」
オイラーのトーティエント関数, φ(n) [時々ファイ関数とも呼ばれる]は,
n と互いに素な n 未満の数の数を定める.
たとえば, 1, 2, 4, 5, 7, そして8はみな9未満で9と互いに素であり, φ(9)=6.


n互いに素な数φ(n)n/φ(n)
2112
31,221.5
41,322
51,2,3,441.25
61,523
71,2,3,4,5,661.1666...
81,3,5,742
91,2,4,5,7,861.5
101,3,7,942.5


n ≦ 10 では n/φ(n) の最大値は n=6 であることがわかる.
n ≦ 1,000,000で n/φ(n) が最大となる値を見つけよ.


-----

まず、問題文のとおりにプログラミングしてみました。
n未満の整数を1つずつnとともに「エラトステネスのふるい」にかけることで
その最大公約数が1ならば互いに素と判定してカウントし、n/φ(n)を計算してみました。
でもこれではn=1,000,000どころか、n=10,000でさえ1分を超え、1分ルールを守れません。


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



#最大公約数を求める(エラトステネスのふるい)
def gcd(a, b):
	if a<b: a, b = b, a
	if b: return gcd(b, a%b)
	else: return a
#n未満でnと互いに素な数の個数
def g(n):
	if n: L = [i for i in xrange(n) if gcd(n, i)<=1]
	else: L = [0]
	return len(L)
def f(n):
	L = [float(i)/g(i) for i in xrange(n+1)]
	m = max(L)
	return L.index(m), m

問題文のまま実装したのではロジック改良でクリアできるレベルではないので、
プログラミング以前に、まずこの問題を数学的に別の表現に言い換えることを考えました。

n/φ(n)を最大値にするためには、nをなるべく大きくしφ(n)はなるべく小さします。

問題文中の表をよくみると、まずnが素数以外であることがわかります。
nが素数の場合、1以上n未満の整数はすべて互いに素なので
φ(n)=n-1になり、他の値と比べてnに対するφ(n)の値が明らかに大きすぎます。

さらに、nが偶数のときにφ(n)の値が小さいことがわかります。
理由は、nが偶数ならば、n未満の偶数すべてnとは互いに素にならないためです。
同様にnが3の倍数ならば、n未満の3の倍数はすべてnと互いに素になりません。
つまり、pを素数とすると、
nがpの倍数ならば、n未満のpの倍数はすべてnと互いに素にならずφ(n)にカウントされません。
また、nの素因数にpを持てば、pの倍数は因数として持たなくてよいです。
つまり、φ(n)の最小値となるnは、
n未満の素数すべてを素因数として1つずつ持つ値となります。

また、n/φ(n) の分母の値の増え方は1ずつですが、
nに素数因子を1つ持つごとのφ(n)の減り方は明らかに1個よりみ大きいので、
分子φ(n)の減少の方が分母nの増加よりも大きく効いてきます。

上記から、φ(n)が最小値となるnは、
素数を小さい順に掛け、その積がn未満の最大値になるような値となります。

ちなみに、問題文のn≦10では、素数は、2, 3, 5, 7,・・・なので、
2×3=6で、n=6のときにn/φ(n) が最大値になります。



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



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

1.h(n)
 引数nが素数ならばその数を、素数以外ならば1を返します。
 問58で作成した素数判定を修正し、素数自身または1を返却することで、
 呼び出し元で簡単に素数だけ掛け算できるようにしました。
 詳細は問58を参照してください。
 
・if n%2 or n==2:
 「n%2」はnを2で割った余りで、論理的にTrueになるのは1の場合、つまり奇数です。
 if文の中は、nが2か奇数の場合の処理です。

・for i in xrange(3, n, 2):
 ループ変数iは3以上n未満の奇数です。2はループ変数の増分です。

2.f(n)
・r = 1
 求める数をクリアします。掛け算していくので初期値は1です。

・for i in xrange(n):
 ループ変数iは0以上n未満の整数です。

・r *= h(i)
 関数h(i)の戻り値を累積して掛け算していきます。

・if r>n: return r//i
 掛け算累積値がnを超えたら、最後に掛けた数で割って
 n未満の掛け算累積値を返します。

3.関数の外
・1,000,000を関数f(n)に渡し、戻り値を表示します。


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

2012年12月23日日曜日

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

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


問68「"magic" 5-gon ring に当てはまる16桁の数字列は何か?」
下に示す図のようなものを"magic" 3-gon ringという.
これは1~6の数字を当てはめて, 各列の数字の和が9となっている.
これを例として説明する.

外側のノードのうち一番小さいものの付いた列(例では4,3,2)から時計回りに回ってそれぞれ列の数字を3つ連ねて説明する.
例えば例のものは4,3,2; 6,2,1; 5,1,3という組で説明することができる.
1~6の数字を当てはめて, 各列の数字の和が等しくなるものは次の8通りある.
合計
94,2,3; 5,3,1; 6,1,2
94,3,2; 6,2,1; 5,1,3
102,3,5; 4,5,1; 6,1,3
102,5,3; 6,3,1; 4,1,5
111,4,6; 3,6,2; 5,2,4
111,6,4; 5,4,2; 3,2,6
121,5,6; 2,6,4; 3,4,5
121,6,5; 3,5,4; 2,4,6
この組の各数字を連結して, 9桁の数字で表すことができる.
例えば, 上の図のものは4,3,2; 6,2,1; 5,1,3であるので432621513である.
さて, 下の図に1~10の数字を当てはめ, 各列の数字の和が等しくなる"magic" 5-gon ringを作って, それを表す16桁または17桁の数字のうち, 16桁のものの最大の数字を答えよ.
(注, 3つの場合の例を見ても分かる通り, 列の始まりの数字を比べた時一番小さい数字で始まる列から時計回りに繋げるという条件のもとで文字列を生成する必要があります. この条件下で最大となる数字を答えてください. )

-----


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






私の解答例は以下です。畳んでいます。
def g(M, Li, Lt):
	t = min(Lt)
	for L in M:
		if t in L: break
	L = [t]+sorted([s for s in L if s!=t], reverse=True)
	R = [L]
	while len(R)<len(M):
		for P in M:
			if sorted(P)!=sorted(L) and (L[-1] in P): break
		L = [s for s in P if s in Lt] + [L[-1]] \
		+ [s for s in P if s in Li and s!=L[-1]]
		R.append(L)
	return "".join([str(s) for L in R for s in L])

def f(n):
	import itertools
	N = n*2
	Q = []
	for Li in itertools.combinations(xrange(1, N+1), n):
		Li, Lt = list(Li), [s for s in xrange(1, N+1) if s not in Li]
		Sa = N*(N+1)//2 + sum([i for i in Li])
		LL = [list(L)+[Sa//n-sum(L)] for L in itertools.combinations(Li, 2)]
		for M in itertools.combinations(LL, n):
			La = [s for L in M for s in L]
			if sorted(La)==sorted(list(Li)+list(Li)+Lt): Q.append(g(M, Li, Lt))
	t = min([len(s) for s in Q])
	return max([int(s) for s in Q if len(s)==t])

n = 5
print f(n)

関数を2つ使用しています。

1.関数g(M, Li, Lt)
・枝リストから、一番外側の数字が最小である枝から時計回りに、
 全ての枝のノードの連結値のうちの最大値を返します。
・引数Mは枝(ノードのリスト)のリスト、Liは内側ノードのリスト、Ltは外側ノードのリストです。

・t = min(Lt)
 tは外側ノードの最小値

・for L in M:
  if t in L: break
 for文で枝リストMからノードリストをループ変数Lとして取り出します。
 if文で外側ノードの最小値がある枝の場合、for文を抜けます。
 ノードリストLは外側ノードの最小値を持つ枝で、最初の枝ということになります。

・L = [t]+sorted([s for s in L if s!=t], reverse=True)
 最初の枝のノードリストLを並び替えます。
 並び順は、外側ノード、その後に内側ノードの大きい順です。
 これで最初の枝のノード並びが確定します。
 sorted文で、内側ノードを大きい順に並び替えています。
 sorted文の第1引数のリストは内包表記しています。
 ノードリストLからループ変数としてノードsを取り出し、後ろのif文に回します。
 if文でsとtが異なる場合、つまり内側ノードの場合、sorted対象のリストの要素としてためます。
 そして第2引数reverse=Trueで、逆順(大きい順)に並べるように指定します。
 外側ノードtを[]でリスト化し、内側ノードリストと結合して、改めてリストLとします。

・R = [L]
 ノード並びが確定したので、最初の枝を結果リストRにためます。

・while len(R)<len(M):
 結果リストRの要素数が元の枝の数になるまで、2本目以降の枝の処理をします。

・for P in M:
  if sorted(P)!=sorted(L) and (L[-1] in P): break
 前の枝の末尾ノードは次の枝にだけ共有されていることを利用して、次の枝Pを求めます。
 「for P in M」で枝リストMから枝をループ変数Pとします。
 if文で枝P枝と1つ前の枝Lとをそれぞれソートして異なること、そして前の枝Lの末尾ノードL[-1]が存在するような枝Pになったところで処理を抜けます。

・L = [s for s in P if s in Lt] + [L[-1]] \
  + [s for s in P if s in Li and s!=L[-1]]
 枝のノードを外側からの順に並び替えます。
 並び替え後ノードリストLは、Pの次のループで前の枝として扱います。

 まず、[s for s in P if s in Lt]は外側ノード1つだけのリストです。
 内包表記です。まず枝Pからノードを1つずつ取り出しループ変数sとして
 後ろのif文に回します。そして外側ノードリストLtにある場合だけ先頭に回します。

 次に、[L[-1]]は1つ前の枝の最後の値だけのリストです。枝Pでは2番目の要素です。

 最後に、[s for s in P if s in Li and s!=L[-1]]は上記2つ以外のノードです。
 内包表記です。まず枝Pからノードを1つずつ取り出しループ変数sとして
 後ろのif文に回します。そして内側ノードリストLiにあり、先ほどの2番目の要素と異なる場合だけ採用し先頭に回します。

 この文全体で上記3要素それぞれだけのリスト3つを+演算子で1つのリストに結合します。

・R.append(L)
 枝のノードが実際の並びになったので、結果リストRに追加してためます。

・return "".join([str(s) for L in R for s in L])
 すべての枝のノードを並び直したので、全ノードを文字扱いして順に連結た値を返します。
 内包表記です。前のfor文で結果リストRから枝Lを取り出し後ろに回します。
 後ろのfor文で枝Lからノードsを取り出し先頭に回します。
 先頭で文字型に変換し、その全要素をリストにします。
 そしてこのリストをjoin関数で区切り文字なしで連結します。

2.関数f(n):
・マジックn角形リングについて、問題で求めるノードの連結値を返します。
 具体的には、マジックn角形リングの外側から内側へ各ノードの値を連結した枝のうち、一番外側の数字が最小の枝から時計回りに、全ての枝数字を連結して最小桁数の最大値を返します。
 引数nは枝の本数でもあります。

・import itertools
 組合せの関数を使用するために、標準モジュールの「itertools」を使えるようにしておきます。

・N = n*2
 Nはノード数です。マジック五角形はノードが10個あります。

・Q = []
 Qは返却する連結値候補リストです。初期クリアしておきます。

・for Li in combinations(xrange(1, N+1), n):
 Liはすべての内側ノードのリストです。
 xrange関数で1からNまでの整数を発生させ、combinations関数でそのうちn個ずつの組合せをLiとします。この時点ではLiはタプルです。

・Li, Lt = list(Li), [s for s in xrange(1, N+1) if s not in Li]
 Liは内側ノードリストで、タプルからリストに変換しておきます。
 Ltは外側ノードリストです。1からNまで整数のうち内側ノードリストに無い値のリストです。
 内包表記です。1からNまで整数をループ変数sとして後ろに回し、if文で内側ノードリストLiに無い値だけ、先頭に回してリストにためます。

・Sa = N*(N+1)//2 + sum([i for i in Li])
 Saはライン総和です。ライン総和には内側ノード2回分と外側ノード1回分が使われるので、「ノード総和+内側ノード総和」としています。
 「N*(N+1)//2」はノード総和です。1からNまでの整数の和の公式そのままです。
 //演算子は割り算で商の整数部を返します。
 「sum([i for i in Li])」は内側ノードLiの総和です。

・LL = [list(L)+[Sa//n-sum(L)] for L in combinations(Li, 2)]
 リストLLは枝候補リストです。内包表記で記述しています。
 まず内側ノードLiからノード2つの組合せタプルをループ変数Lとして先頭に回します。
 list(L)は直後にリストを連結するため、編集不可のタプルから編集可能なリストへ変換しておきます。
 「Sa//n」はライン総和Saを枝数nで割ることで枝ライン1つ分の和のことです。
 [Sa//n-sum(L)]は、この枝1つ分の和から内側ノード2つ分の和sum(L)を差し引くことで、 この枝の外側ノードの値となります。
 list(L)+[Sa//n-sum(L)の全体で、枝候補リストになります。

・for M in combinations(LL, n):
 枝候補リストLLから、枝本数分の組合せを取り出した枝候補組合せリストをループ変数Mとします。

・La = [s for L in M for s in L]
 Laは枝候補組合せリストMの全ノードリストです。内包表記しています。
 枝候補組合せリストMから枝Lを取り出し後ろに回します。
 そして取り出した枝Lからノードsを取り出し、先頭に回しそのままリストにためます。

・if sorted(La)==sorted(Li+Li+Lt): Q.append(g(M, Li, Lt))
 上記の枝候補組合せリストの全ノードについて、内側ノードLiが2回と外側ノードLtが1回使われているかチェックします。
 リストLaと、内側ノードLi2つ分と外側ノードLtの連結リストについて、
 それぞれsorted関数で小さい順に並べて、まったく同じになるかをチェックします。
 チェックOKならば、関数gに枝候補組合せリストM、内側ノードリストLi、外側ノードリストLtを
 渡してます。
 そして、関数gの戻り値である、外側ノード最小値からのノード連結値が最大になるようにを並べ替えた
 連結値を連結値候補リストQにためます。

ここまでで全ての枝候補組合せについて、その連結値候補を取得しました。

・t = min([len(s) for s in Q])
 tは連結候補の最小桁数です。
 連結値候補リストQから連結値をループ変数sとして取り出し先頭に回し、
 len関数で桁数を求めてリストにためます。
 min関数でそのリストの最小値を取得します。

・return max([int(s) for s in Q if len(s)==t])
 連結値候補リストQから連結値をループ変数sとして取り出し、後ろのif文に回します。
 そして最小桁数と一致する場合だけ先頭に回し、int関数で数値化してリストにためます。
 max関数でそのリストの最大値を取得し、呼び出し元に返却します。

上記処理は1秒程度で終わるので1分ルールは守っていますが、
処理速度をさらに向上するアイデアは以下の通りです。
・ライン総和Saがライン数nで割り切れない組は処理を飛ばす。
・枝候補リストLLでの外側ノード値が、外側ノードリストLtにない組は処理を飛ばす。
・枝候補リストLLでの外側ノード値で、外側ノードリストLtの全要素を使ってなければ処理を飛ばす。

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

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