※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

プロジェクトオイラーという数学問題の載ったサイトの問題を堀江伸一こと私がProlog言語で解いていくページ。
Prologは配列と優先順位付キューがなく、std::setの自前実相も使い出が悪いのでそれらがないと効率的に解けない問題をPrologで解くことは難しいです。
そういう場合C++に逃げたりしますが基本Prologでときます。

Problem 131 「素数と立法数の関係」 †

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20131
いくつかの素数pでは, ある正の整数nが存在して, n^3+pn^2が立方数になる.
例えば, p = 19のときには, 8^3+19×8^2=12^3である.
このような性質を持つ各素数について, nの値は一意に定まる. また, 100未満の素数では4つしかこの性質を満たさない.
この性質を持つ100万未満の素数は何個あるだろうか?

解法
この問題は効率的な解法を思いつかなかったので数学掲示板で回答を教えてもらいました。
最大公約数をとることに気付けば自力でもとけたかもしれません。


らすかるさん というかたの回答。
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、、、がある桁数でpで割り切れたとき1111、、、を9倍したものは999、、、
これに1を足してpで割ると、1余る。
最小の桁数はかならずp-1の約数桁となり10^pの約数桁の余りは計算を早くできる。
かつ探す約数は2と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、、、を素数pで割った時割り切れる最短の長さが2と5以外の素因数を持っていたらその数は条件を満たしません。
なのですが私のコードはとても遅いのでたぶん、間違ってはないが頭の悪い解法なのでしょう。
検索して正しい解法を調べるべきですね?
割り切ったレピュニット数を9倍すると999、、、99
それに1足した数を素数pで割ると余りが1になる長さがありそれはフェルマーの小定理よりp-1桁。
すると答えは素数p-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 「素数ペアの結合」 †

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20134
連続する素数 p1 = 19, p2 = 23 について考える. 1219 は末尾の桁が p1 からなり p2 で割り切られる最小の数であることが確かめられる.
実際, p1 = 3, p2 = 5 を除けば, 全ての p2 > p1 なる連続する素数のペアについて, 末尾の桁が p1 からなり p2 で割り切られる数 n が存在する. S を n の最小のものであるとする.
5 ≤ p1 ≤ 1000000 を満たす連続する素数のペア全てに対し ∑ S を求めよ.

解法
19と23なら
100*x+19と分解でき。
19 mod 23=-4
ですから100x mod 23=4となればよい。
オイラーの定理より
100^22 mod 23=1
ですから。
xの一つの答えとして
4*100^21
となる。
これをmod演算の中で解くと最小のxが見つかる
この問題はPrologで速度が出なかったのでC++でもコードを書いてみた。
Prologのほうはコードのほとんどの計算時間が100万以下の素数を求める時間で消費されています。
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 「ピタゴラスタイル」 †

http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%20139
ピタゴラス数を題材にした問題。
詳細はリンク先を参照のこと。
問200くらいまでは結構普通に解ける問題が多いと思うのでそこまではコードを掲載。
それ以上の問題は今後解き方や考え方だけ掲載しようと思ってる。


解法
取り合えず答えが見たかったので、最初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の組から求まる原始ピタゴラス数が答えの元となります。
そして原始ピタゴラス数が求まればそれを自然数倍に相似拡大した三角形は全部この問題の条件を満たす。
かつ三角形が原始ピタゴラス数のとき直角の2辺が1差のものしかこの問題の条件を満たさない。
出てきた答えは以上のような不思議で単純な関係があったのでなぜこれが成り立つか考えてみたが自力ではちょっと考え付きませんでした。

以下はYahoo知恵袋でこの問題についてaerile_reさんというかたに教えていただいた内容を要約したものです。

aerile_reさんによる解説

a^2+b^2=c^2
b-a=kとしここでcがkの倍数であると仮定します。
kは既約なピタゴラス数の性質より奇数となります。
すると
a^2+(a+k)^2=c^2
展開して整理すると
2a^2=c^2-2ka-k^2となりcはkの倍数であると仮定したので
aはkの倍数となります。
bはa+kだったので必然的にbはkの倍数であるとなり,a,b,cがすべてkの倍数となり、既約であるという条件と矛盾します。
よってkは1しかありえません。
  • 解説要約終わり
ここから先MiとNiがペル数になるという条件もあるのですがこれはよくわかりませんでした。

解説の部分まででも十分計算量が落ちているので今のところはここで満足している状態です。

解法
辺の差が1差ですので
1 か -1=m^2-n^2-2mn
としてnを任意の定数としてnを1から計算しmの2次方程式としてとくとm=n+sqrt(2n^2 (+か-) 1)
あとはこれが整数かつピタゴラス数の数であり周長が10^8以下であると確認し、直角三角形の自然数倍の相似拡大の個数を数えて集計すれば答えとなります。

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]).