【プロジェクトオイラー】数学の難問をプログラムの力で解く！！Part2(序盤で使える道具・知識)
こんにちは。じばるどーねのとぼです。今回も「Project Euler（プロジェクトオイラー）」について書いていきたいと思います。...そういえば、「プロジェクトオイラー」のこと間違えて「オイラープロジェクト」って言っちゃうの私だけ？？？
さて、Project Eulerがどういうサイトであるかは前回紹介しました。
とりあえず、プログラムで解く数学の問題がたくさんあると。（ざっくり）
ですが、やっぱ解くだけじゃ直に飽きてしまうかもしれません。
Project eulerでは問題を解くモチベーションにつながる、Level やAward というシステムがあります！！
Levelは0からスタートして、問題を25題解くごとに1上昇します。Levelが1以上になると、自分のLevelのコミュニティに所属できます。さらに、現在ですと、Level2以上になるとサイトに名前が載ります。ついでに、今の私のLevelは2です。解いた問題数は51問でかなりひよっこLevel2ですが。まあ、だんだんと解ける問題が少なくなってきて辛いわけですが、頑張ってLevel3にはなりたいですね。
Level10以上位から超強い領域だと感じています...。
さて、そんなProject Eulerですが、今回は序盤（Level1になるのに必要な）で使える基本的なプログラミングや数学の知識について、実際の問題を見ながら説明していきたいと思います。
1.ユークリッドの互除法
名前だけなら聞いたことがあるという方も多いのではないでしょうか。ユークリッドの互除法は「2つの自然数(整数)の最大公約数を求める」のに使うアルゴリズムです。
例えば、「7920と2475の最大公約数を求めよ」のような問題は、ユークリッドの互除法によって簡単にすることができます。
ユークリッドの互除法のアルゴリズムの基本的な流れを以下に示します。
１．自然数xとyの最大公約数gcd(x,y)を求める(x≧yとする)
２．もしxがyで割り切れるならばgcd(x,y)=yである （終了）
３. xがyで割り切れないならばgcd(x,y)はxをyで割ったときの余りとyの最大公約数
に等しい。xをyで割った余りをx%yと書き表すと、gcd(x,y)=gcd(y,x%y)となる。
４. x=y,y=x%yに置き換えて2.に戻る
以上になります。ここで、yで割ったときの余りは0以上y-1以下の自然数であるので、必ずy>x%yとなることから、4.でx=y,y=x%yに置き換えてもx≧yとなります。しかも、毎回の繰り返しで必ずx,yが小さくなるので、有限ステップでアルゴリズムが終了することも分かります。
では、先ほどの例題「7920と2475の最大公約数を求めよ」をユークリッドの互除法を使って解いてみましょう。
求める最大公約数をgcd(7920,2475)とします。
7920÷2475= 3 余り 495だからgcd(7920,2475)=gcd(2475,495)
2475÷495= 5 余り0で割り切れるからgcd(2475,495)=495
よってgcd(7920,2475)=495です。（終了）
...どうでしょうか。一個一個割れるか試すのではなかなか思いつかない495という答えも比較的簡単に出せてしまいます。
更に、ユークリッドの互除法には以下の２つの応用が利きます。
１．自然数x,yの最小公倍数lcm(x,y)も一緒に求められる。
最小公倍数lcm(x,y)は、要はxとyに共通している約数にx,y固有の約数を掛け合わせたものになりますから(この方法より小さくできる公倍数はない)、
lcm(x,y) = gcd(x,y)×(x÷gcd(x,y))×(y÷gcd(x,y)) = x×y÷gcd(x,y)
ということになります。gcdが求まった時点でlcmも計算１つで求まってしまうわけですね！例えば7920と2475の最小公倍数なら、7920×2475÷495=7920×5=39600より答えは39600となるわけです。
２．複数の自然数の最大公約数にも応用できる
自然数がx,y,zと増えて、この３つの最大公約数gcd(x,y,z)を求める場合でも、結局gcd(x,y,z)はgcd(x,y)とzの最大公約数であるとみなせるので、
gcd(x,y,z)=gcd(gcd(x,y),z)となりわざわざ3つの公約数を見つける必要がなくなります。自然数が4つ、５つ、６つと増えても同様です。例えば7920と2475と231の最大公約数は、「7920と2475の最大公約数」である495と231の最大公約数と等しいため、
495÷231=2余り33,231÷33=7余り０ということで、答えが33であると分かるのです。
<実装>
//ユークリッドの互除法の概略コード
gcd(x,y){
if x<y : swap(x,y) //x≧yとなるように交換する
if x%y==0 : return y //xがyで割り切れたら答えはy
else : return gcd(y,x%y) //gcd(x,y)=gcd(y,x%y)である
}
とまあこんな感じです。概略の作法あんま分かってないので勘弁してください。
<Project Eulerにおける例題>
「1から20までのすべての整数で割り切れる最小の数を求めよ」という問題です。
もし自分で考えたい方はここでストップ。いますぐProject Eulerに登録して解こう！
（数行先で説明）
この問題、つまりは1から20の最小公倍数を求めよって話ですよね。
なので、答えは1~20の最小公倍数をかけていけばいいです。
したがって、はじめ答えをans=1として、lcm(x,y)=x×y÷gcd(x,y)の性質よりiを2から20まで順に増やしていってans=ans×i÷gcd(ans,i)とすれば完了です。
//EP5の解法
ll gcd(ll x, ll y) { //ユークリッドの互除法
if (x < y) swap(x, y);
if (y == 0) return x;
return gcd(y, x%y);
}
int main() {
ll ans = 1; //答え(llは約10^18までの整数を扱える型(これはC++。Pyの方がよい)
for (int i = 2; i <= 20; ++i) {
ll t = gcd(ans, i); //今の答え(lcm)とiの最小公倍数を答えに掛ける
ans *= (i / t);
}
cout << ans << endl;
return 0;
}
２.素数判定（エラトステネスの篩）
素数判定問題は有名ですが、プロジェクトオイラーでも頻出。いかに早く素数列を作れるかが重要だと思います。が私はとてもよわいので今のところおっそいよわよわな生成法しかしていません。競プロだとまずいかもしれないです。Project Eulerの場合は時間制約がないので正解できます。
速いエラトステネスの篩については以下のサイトをご参照ください。
まあ、ある数が素数かどうか判定したいとしますね。
「2525123は素数ですか？」という問題に対して「2で割れない、3で割れない、4で割れない、5で割れない、６で割れない、７で割れない、...」と全部試していったら時間かかりすぎる！ちなみに、2525123は素数です。
「nが素数ですか？」という問題に対して1からnまですべてで割ってみることで判定したらオーダーO(n)。まあ、少し考えると、もし自然数xが素数でないならある自然数z,yを用いてx=z×yと書き表され、z,yの少なくとも一方は√ｘ以下であることから1～√ｎまでの数で割って調べればいいとわかりこれでオーダーはO(√n)に減らせるんですが、まだ遅い。そこでエラトステネスの篩（ふるい）が登場します。
よく考えてみてください。ある自然数nが2で割れないと分かった時点で4や6や8などの偶数で割れないことは明らかです。わざわざ試す必要がない。同様に3で割れないと分かった時点で6,9,12,15などの3の倍数で割れないことは明らか。わざわざ試す必要がないので省略します。つまるところ√n以下の素数で割れるかだけ調べればよくなります。このように「2の倍数でない、3の倍数でない、5の倍数でない、...」と徐々に素数かどうかの判定を絞っていく様が自然数という砂を篩に掛けているようなので、このアルゴリズムを「エラトステネスの篩」といいます。
vectorというのは動的な配列です。つまり、配列の要素数を宣言しなくてもいい配列です。このvectorに素数を格納していきます。
vectorのerstが0のときは素数でない、1のときは素数であることを意味します。ですが、このコードぶっちゃけ全然ふるえていません。精々探索を先に話した方法で√nまでにしたことと(forの条件式dx*dx<=nより)、偶数をスキップしたことくらいです。3の倍数も飛ばすべきかもしれません。ただ、5の倍数も,7の倍数も、と条件を増やしすぎると却って条件式が長くなり遅くなる場合もあります。
//エラトステネスの篩 prime(n)[k]で1~nまで探索して(k+1)番目の素数を出力
vector <int> erstlist;
void prime(int n) {
vector<int> erst;
for (int i = 0; i < n; ++i) erst.PB(i);
erst[0] = erst[1] = 0;
for (int dx = 2; dx*dx <= n; dx++) {
if (erst[dx] == NULL) continue;
for (int i = dx * 2; i < n; i += dx) erst[i] = NULL;
}
for (auto i : erst) if (i) erstlist.PB(i);
return;
}
ぶっちゃけエラトステネスの篩難しいです。ふえ～～～。誰か教えて。
でも、素数ネタはProject Eulerにいっぱいあります。もうほんと、１番多いかもしれない。
＜例題：簡単＞
例えばProblem10は単純に2百万未満の素数をすべて足せ、というもの。まあ、素数表を作って投げておしまいです。
＜例題：ちょい難＞
まず英語がきついので、前回紹介した日本語wikiを使いましょう。
（弱い。でも少しは自分で頑張って読んでみてからwikiは見るし、そうした方がいい）
すると、問題の意味が分かります()
要するにnの整二次式の1次項の係数と定数項をうまく変えて、nに0から順に数字を入れていったときになるべく連続で値が素数になるようにしなさい、という問題です。当然、素数表がないと解けません！！逆に素数表さえあればあとは全探索するだけ！！...とおもって解いたら私は答えが出力されるのに100秒くらいかかるクソコードを生成しました。
まあ、bの絶対値は素数になるなどの工夫をして多少は探索量減らしたんですが...Discussionにはもっと良いコードがありますね。
<おまけ>
最強の方法：時間かかってもいいから素数表を生成し、ライブラリと同じ感覚で素数表配列を予め作っておく。これで実際の問題に対してO(1)で素数表が作れるよ()
３.modと冪・階乗
冪とか階乗とかで無茶苦茶大きい数について下n桁を求めよ、とか各桁の数字を総和を求めよとかとかいう類の問題もプロジェクトオイラー序盤では結構あります。無論、単純に計算しているだけでは桁が大きすぎてオーバーフローしてしまいます。
ここでは、２つの問題のパターン
１．大きな数の下１０桁を求めるパターン
２．大きな数の各桁の総和を求めるパターン
について解法を説明します。なお、２．の方が１．の上位互換です。
まず、1についてざっと説明すると、下１０桁は結局10000000000で割ったあまりですので、計算がオーバーフローしないタイミングで逐次%1000000000（10000000000で割ったあまりを求める）してゆけばいいだけです。いちいち10000000000ってかくのが面倒なのではじめにBIG=10000000000とでも宣言しておきます。
例えば2の1000乗の下１０桁を求めたいとすると以下のようなC++プログラムが書けます。
int main() {
ll BIG = 1000000000;
ll ans = 1; //2^0=1
for (int i = 0; i < 1000; ++i) ans = (ans * 2) % BIG;
cout << ans << endl;
return 0;
}
ll BIG = 1000000000;
ll ans = 1; //2^0=1
for (int i = 0; i < 1000; ++i) ans = (ans * 2) % BIG;
cout << ans << endl;
return 0;
}
ちなみにPythonだと2の1000はprint(2**1000)と打つだけで出力されます。オーバーフローしません。それに一瞬で出ます。Python強くない？
☝Pythonでprint(2**1000)とだけ書いて出力した結果(2の1000乗)
実行結果、2の1000の下10桁は「668069376」でした。
なお、問題のパターンとしては以下に説明する2.のパターンですが2の1000乗にまつわる問題はProject Eulerの16問目にあります。また97問目には巨大素数を下10桁を求める問題があり、今回紹介した方法をそのまま用いて突破できます。
次に２についてですが、そもそもにおいてこういった問題はC++ではなくPythonを使えば問題ありません。Project Eulerは時間制約がないですし、Pythonなら実際2の1000乗などの大きい数もオーバーフローしないで計算できますし計算速度も1秒かからない。草。でも、それでもC++でやりたければ方法があります。プログラムに筆算を指せればいいのです。そして各桁を１つの配列要素にして繰り上がりは変数carryで管理します。例えば、C++で80!（1から80までの整数を全てかけたもの）を計算してその各桁の和を求めたいときは以下のようなコードを書けばよいことになります。
//80!の計算
int main() {
ll sum = 0;
int ans[1000] = { 0 }; //念のため1000桁用意
ans[0] = 1; //はじめ最下位は1(0!)
for (int i = 1; i <= 80; ++i) { //1~80までかけていく
int carry = 0; //桁の繰り上がり
for (int j = 0; j < 1000; ++j) {
ans[j]=(ans[j]*i)+carry; //各桁i倍して繰り上がりを足す
carry = ans[j] / 10;
ans[j] %= 10;
}
}
for (int i = 0; i < 1000; ++i) sum += ans[i];
cout << sum << endl;
return 0;
}
なお、Pythonだと80!ごとき、BIGで割らなくても計算できます。以下のような簡単で短いコードを書くだけで全桁計算できます。
ans=1
for i in range(80):
ans*=(i+1)
else:
print(ans)
for i in range(80):
ans*=(i+1)
else:
print(ans)
で、例えば80!の各桁の総和を求めたいとすれば、以下のようなコードを追加します。
sum=0
while ans>=1:
sum+=ans%10
ans//=10
else:
print("The sum is",sum)
while ans>=1:
sum+=ans%10
ans//=10
else:
print("The sum is",sum)
この実行結果が以下です。
ただし、1つ注意点があります。ans//=10の部分をans/=10とするとansが浮動小数点として計算され、小数点以下の切り捨てがされません。Python初心者の私としては驚きました、これ。
でも除算「/」の代わりに整数範囲での除算「//」という演算子があるので無問題ですね！！
。。。(-ω-;)ウーン。Pythonつよいなぁ。
なお、問題設定を100!の桁総和に変えただけの問題がProject Eulerの20問目にあります。
4.まとめ
Project EulerやるにはPythonが有利。機械学習もPython!! なんなら私程度のレベルなら競プロの問題もPythonでいい気がする（時間制約がC++でないと困難であるような定数倍問題はそもそもなかなか解けない）。でもまあ、素数表ならC++の方が早く作れる。でも、極論言えば素数表は一回手元で作ってそれを配列かなんかに初期値として入れとけばいい。あれ？、Python神じゃん！！Python勉強します。皆さんも一緒にPython勉強しましょう。
5.関連記事
Project Euler記事第一弾です。あと、ここまで書いてきて今更気づいたんですが、Project Euler序盤で使える最強の武器「next_permutation」の話し書き忘れました。辞書順問題とかに使えます。これを含め、豊富なライブラリがあるからやっぱC++もいいなぁ...。次回書きます、続く。
Part1
Part3
