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

Problem 26 「逆数の循環節 その1」 †

単位分数とは分子が1の分数である. 分母が2から10の単位分数を10進数で表記すると次のようになる.

1/2 = 0.5
1/3 = 0.(3)
1/4 = 0.25
1/5 = 0.2
1/6 = 0.1(6)
1/7 = 0.(142857)
1/8 = 0.125
1/9 = 0.(1)
1/10 = 0.1

0.1(6)は 0.166666... という数字であり, 1桁の循環節を持つ. 1/7 の循環節は6桁ある.

d < 1000 なる 1/d の中で小数部の循環節が最も長くなるような d を求めよ.

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

解法

正しい答えは出たがあまり自信がありませんが。

こちらのサイトを参考に計算量を抑える方法を試してみました。

http://www-cc.gakushuin.ac.jp/~851051/student06/simana.pdf

また分数を定数倍しても循環節の長さはかわらないとＷｉｋｉに書いてあったので２か5の倍数は計算において無視しています。

上限値10000にした場合
7 ?- time(main26).
[9967,9966]
% 8,569,507 inferences, 1.328 CPU in 1.325 seconds (100% CPU, 6452335 Lips)

上限値１０万にした場合
8 ?- time(main26).
[99989,99988]
% 177,557,842 inferences, 27.125 CPU in 27.117 seconds (100% CPU, 6545911 Lips)

今回のコードでは計算速度を実感したかったので２から上限まですべて検討していますが。
数字を上限から0へと探して見つかった最長循環列の長さが数字より大きくなったら計算を止める
とすれば計算時間は大幅に短縮できます。
ページ下のほうに逆から探索したコードがあります。

divs(N,D,N):-N mod D>0,!.

divs(N,D,Result):-

       !,

       N1 is N//D,

       divs(N1,D,Result).

fai(1,_,Fai,Fai):-!.

fai(N,D,Fai,Result):-

       N<D*D,

       !,

       Result is (Fai*(N-1))/N.

fai(N,D,Fai,Result):-

       N mod D=:=0,

       !,

       divs(N,D,N1),

       Fai1 is (Fai*(D-1))/D,

       D1 is D+(D mod 2),

       fai(N1,D1,Fai1,Result).

fai(N,D,Fai,Result):-

       !,

       D1 is D+1+(D mod 2),

       fai(N,D1,Fai,Result).

mod_pow(N,0,_,_,N):-!.

mod_pow(N,R,Div,Pow10,Result):-

       R mod 2=:=1,

       !,

       N1 is (N*Pow10) mod Div,

       R1 is R//2,

       Pow10_2 is (Pow10*Pow10) mod Div,

       mod_pow(N1,R1,Div,Pow10_2,Result).

mod_pow(N,R,Div,Pow10,Result):-

       !,

       R1 is R//2,

       Pow10_2 is (Pow10*Pow10) mod Div,

       mod_pow(N,R1,Div,Pow10_2,Result).

calc1(N1,R,R):-N1 mod R=:=0.

calc1(N1,R,R1):-N1 mod R=:=0,!,R1 is N1//R.

seed1(R1,Result,Type):-R1 mod 2=:=1,!,member([Result,Type],[[-1,-1],[1,1]]).

seed1(R1,Result,Type):-R1 mod 2=:=0,!,member([Result,Type],[[1,-1],[1,1]]).

list(N,[R1,Type,Mod]):-

       N mod 2>0,

       N mod 5>0,

       fai(N,2,N,N1),

       between(1,N1,R),

       (R*R>N1->!,fail;true),

       calc1(N1,R,R1),

        seed1(R1,Seed,Type),

       mod_pow(Seed,R1,N,10,Mod).

toF(U,U):-U mod 2=:=0,!.

toF(U,U):-U mod 4=:=2.

toF(U,U1):-U mod 4=:=2,!,U1 is U//2.

toF(U,U1):-U1 is U*2.

search_list2(List,U):-

       member([U,1,1],List),

       toF(U,U1),

       member([U1,-1,1],List),

       !.

search_list(N,Len):-

       !,

       findall(E,list(N,E),List),

       msort(List,List1),

       search_list2(List1,Len),

       !

       .

search(1000,Ans):-!,write(Ans).

search(N,[_,AnsLen]):-

       search_list(N,AnsLen1),

       AnsLen<AnsLen1,

       !,

       N1 is N+1,

       search(N1,[N,AnsLen1]).

search(N,Ans):-

       N1 is N+1,

       search(N1,Ans).

main26:-

       search(2,[0,0]).




















逆から探索するコード。
100億を上限としても
଀% 6,267,118 inferences, 0.969 CPU in 0.974 seconds (99% CPU, 6469283 Lips)
と１秒で答えが出る。

divs(N,D,N):-N mod D>0,!.
divs(N,D,Result):-
       !,
       N1 is N//D,
       divs(N1,D,Result).

fai(1,_,Fai,Fai):-!.
fai(N,D,Fai,Result):-
       N<D*D,
       !,
       Result is (Fai*(N-1))/N.
fai(N,D,Fai,Result):-
       N mod D=:=0,
       !,
       divs(N,D,N1),
       Fai1 is (Fai*(D-1))/D,
       D1 is D+(D mod 2),
       fai(N1,D1,Fai1,Result).
fai(N,D,Fai,Result):-
       !,
       D1 is D+1+(D mod 2),
       fai(N,D1,Fai,Result).


mod_pow(N,0,_,_,N):-!.
mod_pow(N,R,Div,Pow10,Result):-
       R mod 2=:=1,
       !,
       N1 is (N*Pow10) mod Div,
       R1 is R//2,
       Pow10_2 is (Pow10*Pow10) mod Div,
       mod_pow(N1,R1,Div,Pow10_2,Result).

mod_pow(N,R,Div,Pow10,Result):-
       !,
       R1 is R//2,
       Pow10_2 is (Pow10*Pow10) mod Div,
       mod_pow(N,R1,Div,Pow10_2,Result).

calc1(N1,R,R):-N1 mod R=:=0.
calc1(N1,R,R1):-N1 mod R=:=0,!,R1 is N1//R.

seed1(R1,Result,Type):-R1 mod 2=:=1,!,member([Result,Type],[[-1,-1],[1,1]]).
seed1(R1,Result,Type):-R1 mod 2=:=0,!,member([Result,Type],[[1,-1],[1,1]]).

list(N,[R1,Type,Mod]):-
       N mod 2>0,
       N mod 5>0,
       fai(N,2,N,N1),
       between(1,N1,R),
       (R*R>N1->!,fail;true),
       calc1(N1,R,R1),
        seed1(R1,Seed,Type),
       mod_pow(Seed,R1,N,10,Mod).

toF(U,U):-U mod 2=:=0,!.
toF(U,U):-U mod 4=:=2.
toF(U,U1):-U mod 4=:=2,!,U1 is U//2.
toF(U,U1):-U1 is U*2.

search_list2(List,U):-
       member([U,1,1],List),
       toF(U,U1),
       member([U1,-1,1],List),
       !.

search_list(N,Len):-
       !,
       findall(E,list(N,E),List),
       msort(List,List1),
       search_list2(List1,Len),
       !
       .


search(N,[AnsN,AnsLen]):-N<AnsLen,!,write([AnsN,AnsLen]).
search(N,[_,AnsLen]):-
       search_list(N,AnsLen1),
       AnsLen<AnsLen1,
       !,
       N1 is N-1,
       search(N1,[N,AnsLen1]).
search(N,Ans):-
       N1 is N-1,
       search(N1,Ans).
main26:-
       search(10000000000,[0,0]).