プロジェクトオイラーという数学問題の載ったサイトの問題を堀江伸一こと私がProlog言語で解いていくページ。
Prologは配列と優先順位付キューがなく、std::setの自前実相も使い出が悪いのでそれらがないと効率的に解けない問題をPrologで解くことは難しいです。
そういう場合C++に逃げたりしますが基本Prologでときます。
Problem 131 「素数と立法数の関係」 †
解法
この問題は効率的な解法を思いつかなかったので数学掲示板で回答を教えてもらいました。
最大公約数をとることに気付けば自力でもとけたかもしれません。
らすかるさん というかたの回答。
u^3-n^3=pn^2
u=bg, n=ag, gはuとnの最大公約数で1＜a＜bとすると
g(b^3-a^3)=pa^2
b^3-a^3はaと互いに素だからgがa^2の倍数。
g=a^2hとおくと
h(b^3-a^3)=p
b^3-a^3＞1だから、h=1,b^3-a^3=(b-a)(b^2+ab+a^2)=p
よって b-a=1,b^2+ab+a^2=p
すなわち (a+1)^2+a(a+1)+a^2=p
整理して 3a^2+3a+1=p
またh=1なのでg=a^2、従ってn=a^3,u=a^2(a+1)
元の式に代入すると
{a^2(a+1)}^3-(a^3)^3=(3a^2+3a+1)a^6
という恒等式になります。
3a^2+3a+1＜1000000 を解くと a＜(√133333-1)/2＜577 なので
a=1～576について3a^2+3a+1が素数になるかどうかを調べればいいですね。
多少工夫するとしたら、
a≡1,5 (mod 7) のとき 3a^2+3a+1 は7の倍数になりますので
候補の数を5/7には減らせます。
また個数が412個以下であることもわかります。
a=576のとき3a^2+3a+1=997057=(素数)なので
1000000以下ではこれが最大。
でした。
これに従ったコードは以下の通り。
これはエレガントな解法であやかりたいものです。
not_prime(N):-N<2,!.
not_prime(N):-
Limit is floor(sqrt(N)),
between(2,Limit,M),
N mod M=:=0,
!.
is_prime(N):-not(not_prime(N)).
test(P):-
between(1,576,A),
P is 3*A^2+3*A+1,
is_prime(P).
main131:-
findall(A,test(A),Ans),
write(Ans),
length(Ans,Len),
write(Len).
Problem 132 「巨大なレピュニットの因数」 †
http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20132
1のみからなる数をレプユニットという. R(k) を長さ k のレプユニットとする.
例えば, R(10) = 1111111111 = 11×41×271×9091 となり, 素因数の和は9414となる.
R(10＾9) の最初の40個の素因数の和を求めよ.
解法
11111、、、がある桁数でｐで割り切れたとき1111、、、を9倍したものは999、、、
これに1を足してｐで割ると、１余る。
最小の桁数はかならずｐ-1の約数桁となり10＾ｐの約数桁の余りは計算を早くできる。
かつ探す約数は２と5の倍数でなくてはいけないのでそれだけ探す。
これで結構高速化と思って実装、何か遅いと思ったらnot_prime述語を書き間違えていた。
notPrime(P):-
between(2,P,N),
(N*N > P -> !,fail;true),
P mod N=:=0,
!.
isPrime(P):-
not(notPrime(P)).
isOKNum(P):-
isPrime(P),
P mod 2>0,
P mod 5>0.
modPow(_,_,0,Sum,Sum):-!.
modPow(P,Pow10,Num,Mult,Result):-
Num mod 2=:=1,
!,
Num1 is Num//2,
Pow10_1 is (Pow10*Pow10) mod P,
Mult1 is (Mult*Pow10) mod P,
modPow(P,Pow10_1,Num1,Mult1,Result).
modPow(P,Pow10,Num,Mult,Result):-
Num1 is Num//2,
!,
Pow10_1 is (Pow10*Pow10) mod P,
modPow(P,Pow10_1,Num1,Mult,Result).
bai2or5(P,P1):-0 < (P mod 5) ,P1 is P*2.
bai2or5(P,P1):-P1 is P*5.
f(Div,P,P,_):-
P =< Div,
!,
fail.
f(_,_,MinLen,MinLen):-
10^9 mod MinLen=:=0,
!.
f(Div,P,MinLen,Result):-
PM is P-1,
PM mod Div=:=0,
!,
modPow(P,10,Div,1,Amari1),
(Amari1=:=1 -> min(MinLen,Div,MinLen1);MinLen1 is MinLen),
bai2or5(Div,Div1),
f(Div1,P,MinLen1,Result).
min(A,B,A):-A<B,!.
min(_,B,B):-!.
search(_,40,AnsList,AnsList):-!.
search(P,Count,AnsList,Result):-
isOKNum(P),
findall(Len,f(1,P,P,Len),MinLens),
length(MinLens,Len),
Len>0,
!,
Count1 is Count+1,
P1 is P+1,
search(P1,Count1,[P|AnsList],Result).
search(P,Count,AnsList,Result):-
P1 is P+1,
!,
search(P1,Count,AnsList,Result).
sum([],0):-!.
sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
main132:-
search(5,0,[],AnsList),
write(AnsList),
sum(AnsList,Ans),
write(Ans).
Problem 133 「レピュニットの非因数」 †
1のみからなる数をレピュニットという. R(k) を長さ k のレピュニットとする. 例えば, R(6) = 111111 となる.
R(10＾n) というレピュニットについて考える.
R(10), R(100), R(1000) は 17 では割り切れないが, R(10000) は 17 で割り切られる. さらに, R(10＾n) が 19 で割り切られるような n は存在しない. 驚くべきことに, R(10＾n) の因数となりうる100未満の素数は 11, 17, 41, 73 の4個のみである.
R(10＾n) の因数となりえない100000未満の素数の和を求めよ.
1111、、、を素数ｐで割った時割り切れる最短の長さが２と5以外の素因数を持っていたらその数は条件を満たしません。
なのですが私のコードはとても遅いのでたぶん、間違ってはないが頭の悪い解法なのでしょう。
検索して正しい解法を調べるべきですね？
割り切ったレピュニット数を9倍すると999、、、99
それに1足した数を素数ｐで割ると余りが1になる長さがありそれはフェルマーの小定理よりp-1桁。
すると答えは素数ｐ-1の約数だけを調べればよい。
これで10秒まあまあの速度です。
notPrime(P):-
between(2,P,N),
(P<N*N -> !,fail;true),
P mod N=:=0,
!.
isPrime(P):-
not(notPrime(P)).
modPow(_,_,0,Sum,Sum):-!.
modPow(P,Pow10,Num,Mult,Result):-
Num mod 2=:=1,
!,
Num1 is Num//2,
Pow10_1 is (Pow10*Pow10) mod P,
Mult1 is (Mult*Pow10) mod P,
modPow(P,Pow10_1,Num1,Mult1,Result).
modPow(P,Pow10,Num,Mult,Result):-
Num1 is Num//2,
!,
Pow10_1 is (Pow10*Pow10) mod P,
modPow(P,Pow10_1,Num1,Mult,Result).
yakusu(P,Div,P1):-P mod Div=:=0,P1 is P // Div.
yakusu(P,Div,Div):-P mod Div=:=0.
f(Div,P,_):-
Div2 is Div*Div,
Div2>P,
!,
fail.
f(Div,P,P1):-
PM is P-1,
yakusu(PM,Div,P1),
modPow(P,10,P1,1,Amari),
Amari=:=1.
f(Div,P,Result):-
!,
Div1 is Div+1,
f(Div1,P,Result).
min(A,B,A):-A<B,!.
min(_,B,B):-!.
div2or5(P,P1):-P mod 5=:=0,!,P1 is P//5.
div2or5(P,P1):-P mod 2=:=0,!,P1 is P//2.
is2or5(1):-!.
is2or5(P):-div2or5(P,P1),
is2or5(P1).
arrayMin([X],X):-!.
arrayMin([X|Xs],Result):-arrayMin(Xs,Re),min(X,Re,Result).
searchMin(P):-
findall(Len,f(1,P,Len),Lens),
arrayMin(Lens,MinLen),
!,
not(is2or5(MinLen)).
all_search(P):-
between(11,100000,P),
isPrime(P),
searchMin(P).
sum([],0):-!.
sum([X|Xs],Result):-sum(Xs,Re),Result is Re+X.
main133_1:-
findall(P,all_search(P),Ps),
sum(Ps,Ans),
Ans1 is Ans+2+3+5+7,
write(Ans1).
Problem 134 「素数ペアの結合」 †
解法
19と23なら
100*x+19と分解でき。
19 mod 23=-4
ですから100ｘ mod 23=4となればよい。
オイラーの定理より
100^22 mod 23=1
ですから。
ｘの一つの答えとして
4*100^21
となる。
これをmod演算の中で解くと最小のｘが見つかる
この問題はＰｒｏｌｏｇで速度が出なかったのでＣ++でもコードを書いてみた。
Ｐｒｏｌｏｇのほうはコードのほとんどの計算時間が１００万以下の素数を求める時間で消費されています。
C++ time 0.532sec
Prolog time 13.397sec
c++版
#include<stdio.h>
#include<iostream>
#include<string.h>
#include<time>
const int Limit =1000090;
bool is_prime[Limit];
void prime_list(){
memset(is_prime,true,sizeof(is_prime));
is_prime[0]=is_prime[1]=false;
int add;
for(int i=2;i<Limit;i++){
if(is_prime[i]==false)continue;
if(i%2==0)add=i;
else add=i*2;
for(int j=i+add;j<Limit;j+=add){
is_prime[j]=false;
}
}
}
__int64 base(int p1){
__int64 Base=1;
while(Base<p1)Base*=10;
return Base;
}
__int64 mod_pow(__int64 p1,__int64 p2){
__int64 Base,Pow;
Base=Pow=base(p1);
__int64 AllPow=1,Sa=p2-p1,R=p2-2;
while(R>0){
if(R % 2==1){
AllPow=(AllPow*Pow) % p2;
}
Pow=(Pow*Pow) % p2;
R/=2;
}
return ((AllPow*Sa) % p2)*Base+p1;
}
int main(){
clock_t start,end;
start = clock();
prime_list();
__int64 ans=0,T;
int p2;
for(int p1=5;p1<1000*1000;p1+=2){
if(is_prime[p1]==false)continue;
for(p2=p1+2;is_prime[p2]==false;p2+=2){
}
ans+=mod_pow(p1,p2);
}
end= clock();
std::cout<<ans<<"\n";
std::cout<<(double)(end-start)/CLOCKS_PER_SEC<<"秒かかりました";
}
prolog版
not_prime(2):-!,fail.
not_prime(3):-!,fail.
not_prime(N):-
Limit is floor(sqrt(N)),
between(2,Limit,D),
(N mod D)=:=0,
!.
is_prime(N):-not(not_prime(N)).
base(P,10):- P<10,!.
base(P,100):- P<100,!.
base(P,1000):- P<1000,!.
base(P,10000):- P<10000,!.
base(P,100000):- P<100000,!.
base(P,1000000):-P<1000000,!.
mod_pow(0,_,_,Result,Result):-!.
mod_pow(R,P2,Pow,PowAll,Result):-
R mod 2=:=1,
!,
R1 is R//2,
Pow1 is (Pow*Pow) mod P2,
PowAll1 is (PowAll*Pow) mod P2,
mod_pow(R1,P2,Pow1,PowAll1,Result).
mod_pow(R,P2,Pow,PowAll,Result):-
!,
Pow1 is (Pow*Pow) mod P2,
R1 is R//2,
mod_pow(R1,P2,Pow1,PowAll,Result).
searchN(P2,Base,P1,Result):-
Sa is P2-P1,
P22 is P2-2,
%X is (Base^(P2-2)*Sa) mod P2,
mod_pow(P22,P2,Base,1,T),
X is (T*Sa) mod P2,
Result is X*Base+P1.
searchP([P1,_],Ans):-1000000=<P1,!,write(Ans).
searchP([P1,P2],Ans):-
is_prime(P2),
!,
(P2 mod 1000<10->write([P2]),nl;true),
base(P1,Base),
searchN(P2,Base,P1,Re),
Ans1 is Ans+Re,
P3 is P2+1,
searchP([P2,P3],Ans1).
searchP([P1,P2],Ans):-
P3 is P2+1,
searchP([P1,P3],Ans).
main134:-
searchP([5,7],0).
Problem 139 「ピタゴラスタイル」 †
解法
取り合えず答えが見たかったので、最初Wikiに書いてある通りの原始ピタゴラス数の求め方で全探索しました。
取り立てて遅いというわけではないがちょっと遅い処理になりました。
出てきた答えを見ると、
Wikiの原始ピタゴラス数を求める関数a=M^2-N^2,b=2MN,c=M^2+N^2として
この問題の条件を満たすMi,NiはM1=2,N1=1として
Mi+1=2Mi+Ni
Ni+1=Mi
MiとNiの組から求まる原始ピタゴラス数が答えの元となります。
そして原始ピタゴラス数が求まればそれを自然数倍に相似拡大した三角形は全部この問題の条件を満たす。
かつ三角形が原始ピタゴラス数のとき直角の２辺が１差のものしかこの問題の条件を満たさない。
出てきた答えは以上のような不思議で単純な関係があったのでなぜこれが成り立つか考えてみたが自力ではちょっと考え付きませんでした。
以下はYahoo知恵袋でこの問題についてaerile_reさんというかたに教えていただいた内容を要約したものです。
aerile_reさんによる解説
a^2+b^2=c^2
b-a=kとしここでcがｋの倍数であると仮定します。
ｋは既約なピタゴラス数の性質より奇数となります。
すると
a^2+(a+k)^2=c^2
展開して整理すると
2a^2=c^2-2ka-k^2となりｃはｋの倍数であると仮定したので
aはkの倍数となります。
ｂはa+kだったので必然的にｂはkの倍数であるとなり,a,b,cがすべてｋの倍数となり、既約であるという条件と矛盾します。
よってkは１しかありえません。
ここから先MiとNiがペル数になるという条件もあるのですがこれはよくわかりませんでした。
解説の部分まででも十分計算量が落ちているので今のところはここで満足している状態です。
解法
辺の差が１差ですので
1 か -1=m^2-n^2-2mn
としてｎを任意の定数としてｎを１から計算しｍの２次方程式としてとくとm=n+sqrt(2n^2 (+か-) 1)
あとはこれが整数かつピタゴラス数の数であり周長が１０＾８以下であると確認し、直角三角形の自然数倍の相似拡大の個数を数えて集計すれば答えとなります。
calc1(N,T):-T is 2*N*N+1.
calc1(N,T):-T is 2*N*N-1.
gcd(0, B, G) :- G is abs(B).
gcd(A, B, G) :- A =\= 0, R is B mod A, gcd(R, A, G).
sum([],0):-!.
sum([[_,Perm]|Rest],Result):-sum(Rest,Re),Result is Re+Perm.
ok(N,[[M,N,A,B,C],Perm]):-
calc1(N,T),
T1 is floor(sqrt(T)),
T=:=T1*T1,
M is N+T1,
M>N,
1=:=(M-N) mod 2,
gcd(M,N,1),
A is M^2-N^2,
B is 2*M*N,
C is M^2+N^2,
All is A+B+C,
10^8>All,
Perm is (10^8-1)//All.
roopN(N,_):-
M is N+1,
10^8=<2*M*(M+N),
!,
fail.
roopN(N,Result):-
ok(N,Result).
roopN(N,Result):-
N1 is N+1,
roopN(N1,Result).
main139:-
findall(Ans,roopN(1,Ans),Answers),sum(Answers,Ans1),
write([ans,Ans1]).
最終更新：2014年02月12日 16:40
