プロジェクトオイラー Problem 64 「奇数周期の平方根」 †

Problem 64 「奇数周期の平方根」 †

平方根は連分数の形で表したときに周期的であり, 以下の形で書ける:

√N = a0 + 1 / (a1 + 1 / (a2 + 1 / (a3 + ...)))

例えば, √23を考えよう.

√23 = 4 + √23 - 4 = 4 + 1 / (1 / (√23 - 4)) = 4 + 1 / (1 + (√23 - 3) / 7)

となる.

この操作を続けていくと,

√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 について奇数の周期をもつ平方根が何個あるか答えよ.

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2064




解法
連分数について人生で一度も教育を受けたことがないので
数列を眺めて直感で無理矢理漸化式を立てて計算した。
出鱈目だな俺。
とりあえず答えはあった。

gcd(0, B, G) :- G is abs(B).
gcd(A, B, G) :- A =\= 0, R is B mod A, gcd(R, A, G).

calc_next(N,Up,Dell,D2,NextDell,[Add,NextDell,D2]):-
!,
D1 is N-Dell*Dell,
gcd(Up,D1,G1),
NextUp is Up//G1,
D2  is D1//G1,
Add is ((floor(sqrt(N))+Dell)*NextUp)//D2,
NextDell is Add*D2-Dell.

calc(N,Up,Dell,List):-
calc_next(N,Up,Dell,NextUp,NextDell,E),
(member(E,List) ->
!,length(List,Len),
Len mod 2=:=1;
append(List,[E],List1),
calc(N,NextUp,NextDell,List1)).



search(1):-
between(1,10000,N),
N1 is sqrt(N),
not(N1=