Problem 29 「個別のべき乗」 †
2 ≤ a ≤ 5 と 2 ≤ b ≤ 5について, ab を全て考えてみよう:
- 22=4, 23=8, 24=16, 25=32
- 32=9, 33=27, 34=81, 35=243
- 42=16, 43=64, 44=256, 45=1024
- 52=25, 53=125, 54=625, 55=3125
これらを小さい順に並べ, 同じ数を除いたとすると, 15個の項を得る:
4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125
2 ≤ a ≤ 100, 2 ≤ b ≤ 100 で同じことをしたときいくつの異なる項が存在するか?
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2029
Prolog言語
計算量
BigO(99*6*２)
重複する可能性のある数字は
2,4,8,16,32,64
3,9,27,81
5,25
6,36
7,49
10,100
で
2と4は指数だけに注目すると２倍の関係だし
２と8も指数だけに注目すると３倍の関係。
、、、
他の行も同様。
これを考えて全組み合わせから重複分を引いていくと答えが出ます。
code
Prolog
a,bの上限を変えて実行した結果。
2<=a<=100,2<=b<=100
17 ?- time(main29).
9183
% 6,352 inferences, 0.000 CPU in 0.001 seconds (0% CPU, Infinite Lips)
true.
2<=a<=1000,2<=b<=1000
18 ?- time(main29).
977358
% 118,272 inferences, 0.031 CPU in 0.021 seconds (149% CPU, 3784704 Lips)
2<=a<=10000,2<=b<=10000
19 ?- time(main29).
99347607
% 2,040,518 inferences, 0.391 CPU in 0.403 seconds (97% CPU, 5223726 Lips)
2<=a<=100000,2<=b<=100000
21 ?- time(main29).
9981236306
% 29,087,619 inferences, 5.656 CPU in 6.108 seconds (93% CPU, 5142562 Lips)
詳細な解法
aを素因数分解してpiを素数,riを指数とします。
a=p1^r1*,,,*pn^rn
と考えます。
f(a)をaをベクトルに変える関数だと考えます。
a=2^3*5^2*7^1
f(a)=[3,0,2,1]
a=2^3*3^1*5^2
f(a)=[3,1,2]
a=3^2
f(a)=[0,2]
素数を小さいほうから一つずつ次元を対応させて指数がその次元のベクトル成分だとします。
すると、a1^b1=a2^b^2となるのは自然数k1、k2を考えて
k1f(a1)=k2f(a2)
となる場合だけであることがわかります。
これを考えると
ｇｃｄ（ｒ１,r2,,,,rn）=1
となるものの定数倍で表現されるものだけが互いに重複する可能性があります。
よって
2<=a<=100
2<=b<=100
なら
2,4,8,16,32,64
3,9,27,81
5,25
6,36
7,49
10,100
で上のどの行も一番小さい数字をベクトルにするとその行の残りは一番小さい数字の定数倍ベクトルとなります。
a=11以上は２乗すると１００を超えるので考える必要がありません。
定数倍だけ考えればいいので
5,25
6,36
7,49
10,100
は定数倍の世界では同じ計算になります。
2<=bであることに注意して、一番短いベクトルの長さを1としてこれをｂ倍したもの
2,3,,,,99,100
4,6,,,,198,200
で数字の出現する回数を数えて 各数字の出現回数-1の総計が
5,25
6,36
7,49
10,100
で作れた数字が重複した回数となります。
3なら3,9,2,7,81と4つですから４つめまで考えて。
2,3,4、、、99,100
4,6,8、、、198,200
6,9、、、、、297,300
8,12、、、、396,400
で数字の出現する回数を数えれば
3,9,27,81の重複回数がわかります。
2の場合も6の倍数で２～100倍までを計算するだけで同じです。
この手法で
2=<a<=10万
2=<b<=10万
という小さな数程度まではそれなりの計算時間で計算できます。
私は素朴に実装したので２＝＜a<=100万,2=<b=<100万はメモリが足りませんでした。
賢い人なら計算の高速化やメモリ使用量の大幅低減ができるでしょう。
prolog
%堀江伸一
limit_a(100001):-!.
limit_b(100001):-!.
add_next(B,_,List,List):-
limit_b(B),
!.
add_next(B,M,[[N1,C]|Rest],[[N1,C]|Result]):-
N1<B*M,
!,
add_next(B,M,Rest,Result).
add_next(B,M,[[N1,C]|Rest],[[N1,C1]|Result]):-
N1=:=B*M,
!,
C1 is C+1,
B1 is B+1,
add_next(B1,M,Rest,Result).
add_next(B,M,List,[[BM,1]|Result]):-
!,
B1 is B+1,
BM is B*M,
add_next(B1,M,List,Result).
count([],0):-!.
count([[_,C]|Rest],Result):-
!,
count(Rest,Re),
Result is Re+C-1.
dp(_,[],_,Ans):-!,
write(Ans).
dp(M,[[M,Mult]|Rest],List,Ans):-
!,
add_next(2,M,List,List1),
count(List1,Count),
Ans1 is Ans-Mult*Count,
M1 is M+1,
dp(M1,Rest,List1,Ans1).
dp(M,Mults,List,Ans):-
!,
add_next(2,M,List,List1),
M1 is M+1,
dp(M1,Mults,List1,Ans).
divs(N,Div,0,N):-N mod Div>0,!.
divs(N,Div,ResultC,ResultN):-
!,
N1 is N//Div,
divs(N1,Div,ReC,ResultN),
ResultC is ReC+1.
gcd(-1,B,B):-!.
gcd(0, B, G) :- G is abs(B).
gcd(A, B, G) :- A =\= 0, R is B mod A, gcd(R, A, G).
gcds(1,_,GCD,GCD):-!.
gcds(N,Div,_,1):-
N<Div^2,
!.
gcds(N,Div,GCD,Result):-
N mod Div=:=0,
!,
divs(N,Div,GCD1,N1),
gcd(GCD,GCD1,GCD2),
Div1 is Div+1,
gcds(N1,Div1,GCD2,Result).
gcds(N,Div,GCD,Result):-
!,
Div1 is Div+1,
gcds(N,Div1,GCD,Result).
seedA(N):-
limit_a(Limit),
between(2,Limit,N),
(Limit=<N^2->!,fail;true),
gcds(N,2,-1,GCDs),
GCDs=:=1.
max_r(AR,_,0):-
limit_a(Limit),
Limit=<AR,
!.
max_r(AR,A,Result):-
!,
AR1 is AR*A,
max_r(AR1,A,Re),
Result is Re+1.
a_max_r(Nums,R):-
member(E,Nums),
max_r(E,E,R).
countR([],[R,Count],[[R,Count]]):-!.
countR([R|Rest],[R,Count],Result):-
!,
Count1 is Count+1,
countR(Rest,[R,Count1],Result).
countR([R|Rest],[R1,Count1],[[R1,Count1]|Result]):-
!,
countR(Rest,[R,1],Result).
main29:-
X is 512*1024*1024,set_prolog_stack(global,limit(X)),
findall(N,seedA(N),Nums),
findall(R,a_max_r(Nums,R),Rs),
msort(Rs,[Top|Rs1]),
countR(Rs1,[Top,1],RList),
limit_a(A),
limit_b(B),
S is (A-2)*(B-2),
dp(1,RList,[],S).
0 件のコメント:
コメントを投稿