プロジェクトオイラー 問１１０

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

問１１０「ディオファントス逆数 その2」
次の等式で x, y, n は正の整数である.

1/x + 1/y = 1/n

n = 1260 では 113 の異なる解があり, この n が解の個数が 100 を超える最小の値である.

解の数が 4,000,000 を超える最小の n を求めよ.

注: この問題は Problem 108 を非常に難しくしたケースである. 総当り法で解ける範囲を超えているので, 賢い解き方が求められる.






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






私の回答例は以下の通りです。

def P(n):
	L, i = [], 1
	if n>0:
		L.append(2)
		while len(L)<n:
			i += 2
			for j in L:
				if not i%j: break
			else:
				L.append(i)
	return L

def mul(L): return reduce(lambda x,y:x*y, L)

def f(m):
	import math
	tmax = int(math.log(m*2, 3))
	p = P(tmax+1)
	nb, kb, Lb, L = 0, 0, [], []
	while len(L)<=tmax:
		L = [1]*(len(L)+1)
		M = [i*2+1 for i in L]
		k = (mul(M)+1)//2
		while k<=m:
			for i in xrange(len(L)):
				if tmax<=1 or (i<len(L)-1 and L[i]==L[i+1]):
					L[i] += 1
					for j in xrange(i):
						if L[j]>L[i]: L[j] = L[i]
					break
			else:
				L = [1]*(len(L)+1)

			M = [i*2+1 for i in L]
			k = (mul(M)+1)//2

		N = [p[i]**j for i,j in enumerate(L)]
		n = mul(N)
		if nb==0 or n<=nb: nb, kb, Lb = n, k, L
	return nb, kb, Lb

m=4000000
n, k, L = f(m)
print "n :", n
print "k :", k
print "L :", L

問１０８ではnを１つずつ増やしながら解の個数kを求めてkが閾値を超えたら終了というロジックで解きましたが、解の個数kが巨大な数になると、素因数分解の処理が重くて処理が終わりません。

まず、問１０８からわかっていることは以下のとおりです。
・解の個数kは、n²を２つの約数の積で表したパターンの個数に等しい。

・n²の約数の個数iについて、n²を素因数分解した結果に含まれる各素因数に対して、採用0個の場合も含めたパターン数は各要素数+1であり、これを全素因数の分だけ掛け合わせた積が、約数の個数iということになります。
　例：72の約数の個数は、72=2³＋3², べき乗の値から、(3+1)x(2+1)=12

・n²の約数の個数iとして、２つの約数の積でn²になる組合せの個数kとして、
　k = (i+1)//2  （小数以下切捨て）です。

上記からnと解の個数kの関係を、例にあるn=1260, k=113の場合を手で計算してみると以下の通りです。

n = 1260
↑積　　↓因数分解 
素数とべきリスト値の累乗リストN = [2², 3², 5¹, 7¹] = [4, 9, 5, 7]
↑小さい順の素数にべきリストLの値で累乗
べきリストL = [2, 2, 1, 1]
↓各要素の２乗
べき２乗リスト = [4, 4, 2, 2]
↓各要素に+1する
約数の個数リストM = [5, 5, 3, 3]
↓各要素の積
約数の個数 = 5*5*3*3 = 225
↓
解の個数k = (225+1)//2 = 113

問108ではnから出発して１つずつ因数分解をしていましたが、
一見して明らかに一番取り扱い易そうなのは、上記のべきリストLです。
べきリストLを調整しつつ、kを算出してしきい値判定して、nを特定する方法でいきます。

しきい値mとして、m<kである最小のnを求めるので、
m<kを満たすkが特定されたときにnが最小であるためには、
べきリストLは、左の要素≧右の要素 ...(a)

べきリストLの要素数をtに固定すると、
nが最小になるのは、全要素が1のとき。例　[1,1, ... ,1]　1がt個ある。
べきリストLの要素数tの最大値tmaxは、全要素が1でm=kのときのtの値です。
このとき約数の個数リストMは1*2+1=3だけがtmax個ある状態で、
解の個数k = (3^tmax)//2 = m（しきい値）となります。
3^tmax = m*2
両辺の対数をとって、
tmax * log3 =log(m*2)
tmax = log(m*2)/log3 = log₃(m*2) （小数以下切捨て）...(b)

１．関数P(n)
・小さい方からn個の素数リストを返します。
・Problem108で作成した関数を流用します。詳しくはProblem108を参照。

２．関数mul(L)
・リストLの全要素の積を返します。
・reduce文で、第１引数の関数f、第２引数にリストLを設定しています。
・関数fは「lambda x,y:x*y」です。
　無名関数lambdaを使って、引数１と引数２の積を計算します。
・reduce文では関数fの第２引数にリストLの要素を順に渡し、
　関数fの第１引数に関数処理結果を設定します。
　このため、関数fでfの第１引数と第２引数を元に計算することで、累積的な処理が記述できます。
　reduce文は、累積的は処理を短く記述できますが、少しややこしい構文です。

３．関数f(m)
・1/x + 1/y = 1/n（x,y,nは正の整数）となる解の個数kが引数mを超える最小値を求め、n、k、べきリストLを返します。
・べきリストLの要素数の最大値tmaxは、上記(b)をそのまま実装します。
・素数リストpは、関数Pでtmax個の素数を準備します。

・外側のwhileループでは、べきリストの要素数がtmax以下の場合にべきリストを調整して解の候補を探します。
・べきリストLの初期値は全要素1で要素数をループの都度１つずつ増やします。

・内側のwhileループでは、解の個数kがそのしきい値m以下の場合にべきリストを調整して解の候補を探します。

・現在のべきリストに基づいて次のべきリストを調整します。
・iのforループで、べきリストの要素を右から順に見ていきます。
　まず、べきリストの最大要素数が1以下、または「右隣と同じ値」である最初の要素を+1します。
　例：[1] → [2]
　例：[4, 3, 2, 2, 1, 1] → [4, 3, 3, 2, 1, 1]
　さらに、+1した要素の左側は、L[i]にクリアします。
　例：[4, 3, 3, 2, 1, 1] → [3, 3, 3, 2, 1, 1]
　-1する前のリストは当ループが進めばまた出てきます。
　
・forループが最後まで終了したら、べきリストLの要素数を１つ増やして全要素1にクリアします。
　for文のelse句は、for文のループがbreakで途中で抜けずに最後まで終了した直後に実行する処理を記述します。
　例：[3, 2, 1] → [1, 1, 1, 1]

・このようにしてべきリストLを作成したら、チェックします。
　約数の個数リストM、解の個数kを計算します。

・解の個数kがしきい値を超えて内側のwhileループが終了したら、素数とべきリスト値の累乗リストN、そして求める値nを計算し、バックアップしてある値と比較して、最小のnの場合にバックアップするようにします。




解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
n : 9350130049860600
k : 4018613
L : [3, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1]