プロジェクトオイラー 問７３

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

問７３「真既約分数をソートした集合内の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の分数を真既約分数として
カウントアップする方法です。
　しかしこの方法では私のＰＣでは１分半かかり１分ルールを守れませんでした。
が、自力ではここまでだったのでそのまま書いておきます。

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

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

２．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の場合、その分数は真既約分数です。


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

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

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

上記を利用した解答例は以下です。畳んでいます。

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

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

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

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

・最初の１つの分数から隣接する大きい値の分数を求める方法は以下です。
　以下の連立方程式を解きます。
　┌ 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から１つずつ大きな分数を計算していき、
　範囲の最大値y1/y2になるまで個数をカウントアップします

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