日本語翻訳サイトは、プロジェクトオイラー日本語 でネット検索です。
問１０８「ディオファントス逆数 その1」
次の等式で x, y, n は正の整数である.
1/x + 1/y = 1/n
n = 4 では 3 つの異なる解がある.
1/5 + 1/20 = 1/4
1/6 + 1/12 = 1/4
1/8 + 1/8 = 1/4
解の数が 1000 を超える最小の n を求めよ.
注: この問題は Problem 110 の易しいケースである. こちらを先に解く事を強く勧める.
注意！！！
自力で解きたい人へ
以降の記述には解法に関するネタバレが含まれます。
私の回答例は以下の通りです。
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 g(n, P):
s, d = n, {}
for i in P:
while s>=i:
if s%i: break
else:
while not s%i:
d[i], s = d.get(i, 0)+1, s//i
if s<i: break
return d
def f(m):
import math
i = int(math.log(2*m-1, 2))
n, k, p = 1, 1, P(i)
while k<=m:
n += 1
d = g(n*n, p)
i = 1
for j in d.values(): i *= j+1
k = (i+1)//2
return n, k
m=1000
n, k = f(m)
print "n :", n
print "k :", k
問題の式を変形してみます。
1/x + 1/y = 1/n
両辺にnxyを掛けて
ny + nx = xy
xy - nx - ny = 0 ...a
式aを２つの数の積で表して差分調整することにします。
xy, nx, nyを展開式の項として持つ式で最も簡単な式は以下の通り。
(x-n)(y-n) = xy - nx - ny + n2 ...b
式bに式aを代入して、
(x-n)(y-n) = n2 ...c
n2を２つの正の整数の積で表し、それぞれ+nすればxとyになります。
つまり、当問題でチェックすべき「解の数」は、
n2を２つの約数の積で表したパターンの個数です。
また、２つの約数の積が元の数になる組合せは、
約数を順に並べて小さい方からと大きい方からそれぞれ１つずつ採用し、
中央値がある場合は中央値を２つ採用した組合せです。
なので、n2の約数の個数i、該当する約数の組合せの個数kとして、以下の通りです。
k = (i+1)/2 （小数点以下切捨て）...d
また、約数の個数iは素因数分解した素因数ごとの採用パターン数の積です。
n2を素因数分解した結果の素因数ごとの要素数jとして、
その採用0個の場合も含めたパターン数j+1について、
全ての素因数の分の積を計算します。
この約数の個数iから、約数ペアの個数kを求め、指定値mを超えるまで処理を繰り返します。
また、素因数分解にあたり素数をいくつまで用意したらいいか検討します。
素因数分解で使用する素数の個数の上限値tについて、
この分だけ最小素数2を累乗しても約数の個数i以上となるので、
2t ≧ i
両辺の対数を取って、
tlog2 ≧ log(i)
t ≧ log(i)/log2 = log2i
式dを代入して、
t ≧ log2(2k-1)
t = log2(2k-1) （小数点以下切り捨て）...e
素因数分解の準備として、素因数分解で使用する素数の個数の上限値tの分だけ素数を用意しておきます。
１．関数P(n)
・小さい方からn個の素数リストを返します。
・リストLは素数リスト、変数iは素数候補です。
・リストLに最初の素数2を設定します。
・3以上の奇数を素数候補として、リストLにある素数で順に割り算します。
割り算の余りが無い時点で素数候補iは素数ではないので、
forループをbreakして次の素数候補へ処理がすすみ、
リストLのすべての素数で割り算の余りがあれば、forループが最後まで回るのでelse句へ進んで素数リストLに素数候補iを追加します。
・上記を指定した素数の個数になるまで続けます。
２．関数g(n, P)
・Problem47で作成した関数を流用します。詳しくはProblem47を参照。
・引数nを素数リストPを使って素因数分解し、
素数をキー、キーごとの分解した個数を値とする辞書（連想配列）を返します。
３．関数f(m)
・1/x + 1/y = 1/n（x,y,nは正の整数）となる解の個数kがmを超える最小値を求め、nとkを返します。
・対数関数を使う準備としてmathモジュールをインポートし式eを実装します。
math.log関数の引数は、真数、底の順に設定します。
変数nは問題文の右辺1/nのn、変数kはチェック対象の解の個数です。
・whileループで、解の個数kが指定値mを超えるまで処理を繰り返します。
・変数dは関数gによるn2に対して素数リストpを使用して因数分解した結果の辞書（連想配列）です。
・ループ関数jで因数分解辞書dから、素因数の値に関係なく、素因数ごとの要素数を設定します。
変数iに、該当する素因数の個数j個に0個採用のケースも含めたj+1を掛けてためます。
・変数kとして、式dを実装します。//関数は割り算の商（整数値）です。
解答はこのすぐ下の行です。文字の色を白にしてます。選択状態にすると見えます。
n : 180180
k : 1013
