じばるどーね!

趣味、思想、妄想、学問、ごった煮ブログ

【プロジェクトオイラー】数学の難問をプログラムの力で解く!!Part3(個人的に面白かった問題)

  こんにちは、ここ数日インフルに倒れていたじばるどーねのとぼです。プロジェクトオイラーの記事です。またある程度解き進めるまでは、このシリーズはとりあえず今回で最後にしようと思います。

 

f:id:tb9810w124816:20180325171613p:plain

オイラープロジェクトについての詳しい説明は省きます。

 

今回は、個人的に序盤で解いて楽しかった問題、為になった問題を4つ紹介します。あ、ちなみに書き溜め記事ってのもあるし、なんだかんだこの記事ではC++を使っています。 

 

 

1. Problem12 Highly divisible trianglar number

f:id:tb9810w124816:20180409093206p:plain

  問題文は上記のようになっています。今回の問題は英語が比較的読みやすい気もしますが、分からなくても例のwikiが使えるので安心です。問題は「三角数の中で500を超える正の約数を持つ最小の数はいくつか」というものです。三角数とは、自然数nを用いて
n*(n+1)/2の形で表される自然数のことを言います。
 
<この問題の面白いところ>
 個人的にこの問題が面白いと 思った点は、計算量を工夫次第で減らせる点です。プロジェクトオイラーとしての問題というより、この類の問題が競技プログラミングで出題されたら少し工夫できていいな、と思ったまでです。
 まず、通常ある自然数Nの約数の個数を数えたい場合、
1から√N以下の整数全てでNが割れるか(約数であるかどうか)試せば十分です。なぜなら、Nがある数iを約数に持つならば、必ず同時にN/iも約数に持つはずであり、i≠√Nのとき明らかにiとN/iのうちの
一方は√N未満の数、他方は√Nより大きい数となるからです。したがって、√N未満で割り切れた整数の数を2倍して、最後に√Nで(for文の条件式をi*i<=Nとすればよい)割り切れれば1を足す、という作業をすればよいのです。
 さて、以上から三角数N=n*(n+1)/2の約数の個数を数えるのにかかる計算量はO(√N)=O(√n²)=O(n)と分かりました。しかし、Nが三角数であることに着目すれば、もう少しオーダーを減らせるのではないでしょうか?
 まず三角数はある自然数nによってn*(n+1)/2と書き表される数のことを言いました。整数論でよく用いられる性質でnとn+1は互いに素(1以外に公約数を持たない)という性質があるのですが、よくみると三角数は互いに素である2数を掛け算していることが分かります。さらにn,n+1のいずれか一方は偶数であることも分かるので、nが偶数の時、Nは互いに素な2数n/2,n+1の積であり、nが奇数のとき、Nは互いに素な2数n,(n+1)/2の積であると結論づけられるのです!互いに素ということは約数にダブりがないので、Nの約数の総数はこの互いに素な2数の約数の個数の積で表されることになります!
 以上の考察から、計算量はO(√n+√n)=O(√n)=O(⁴√N)にまで減らすことができました。この計算量の減少は問題が三角数であったからできたことであり、これで相当大きい三角数(n=10¹⁶)も高速で計算できるようになりました。
 
int main() {
     ll n1=1, n2=1 ;  //n1,n2は互いに素な2数の各約数の個数
     int t = 2;  //t*(t-1)/2=N(答え)
     while (n1*n2<=500) { //約数が500を超えない間
        n1 = n2;  //n1<n2とすると一個前のn2がそのままn1になる
        n2 = 0;  //n2のカウントをリセット
        int po = (t + 1);  //上限
        if (po % 2 == 0) po /= 2; //n/2 OR (n+1)/2の判定
        for (int i = 1; i*i <= po; ++i) if (po%i == 0) n2+=2;
        if(i*i==po) n2--;   //n2が平方数だった時の考慮
        t++;
      }
     cout << (t*(t-1))/2 << endl;   //答え出力
     return 0;
}
 なお、ll はlong long 型のことです!!
 細かいですがソースコードではO(2√n)からO(√n)くらいに減らす工夫をしています。
 

2.Problem33 Digit cancelling fraction

f:id:tb9810w124816:20180409100657p:plain

 お茶目な問題です。49/98=4/8となる理由を昔のアホな数学者が「49と98の共通の文字9を取ったからだ!」と勘違いしてしまったそうですが、40/50=4/5などの自明な(当たり前な)例を除いてこのような2桁/2桁の分数は4つしかないそうで、この4つを約分して積をとったときの分母の値を求めよ、というものです。

 

<この問題の面白いところ>

   「 9を取り除くから49/98=4/8」っていう問題の発想が単純に面白いですよね。気持ちとしては「40-32÷2=4!」並みの面白さを感じます。ただし、この問題、「自明な例(non-trivial examples)」の定義があいまいなんですが、おそらく44/33,20/10など、「10の倍数分の10の倍数」、「11の倍数分の11の倍数」の2パターンのことでしょう。これ以外にある解法を探していきます。全探索で十分でしょう。条件判定が少し難しいですが、手元で少し計算すると、

i) (10a+b)/(10a+c)=b/c

ii)(10a+b)/(10c+b)=a/c

iii)(10a+b)/(10c+a)=b/c

の3パターンは自明な例を除いてあり得ないことが分かります(証明略。)

よって残りのパターン

iV)(10a+b)/(10b+c)=a/c

のみを考えればよいことになり、分数のままでコンピュータが計算するのは難しいので、通分などして整数の条件に直すと

 10a(c-b)=c(a-b)

が成り立つことが条件であると分かります。これでほとんど問題が解けました。

 

<ソースコード>

int main() {
    int nu = 1, den = 1;
    for (int a = 1; a <= 9; ++a) {
        for (int b = 1; b <= 9; ++b) {
             for (int c = 1; c <= 9; ++c) {
                   if(a!=b&&b!=c) //自明な例11の倍数を除く
                   if (10 * a*(c - b) == c * (a - b)) {  //条件判定
                       nu *= a;  //答えを分子に掛ける
                      den *= c; //答えを分母に掛ける
                      cout << 10 * a + b << "/" << 10 * b + c << endl;
                    }
               }
          }
      }
     int t = gcd(nu, den);  //nu(分子)とden(分母)の約分
     cout << den / t << endl;
     return 0;
}
 ちなみにこの問題は答えも美しいので是非やってみてください。

 

3.Problem41 Pandigital prime

f:id:tb9810w124816:20180409103349p:plain

  「1からn(nは9以下の自然数)が1回ずつ各桁に現れる数をn桁パンデジタル数と呼ぶとき、最大のパンデジタル素数を求めよ」という問題です。英語圏の言葉遊び「パングラム」(a~zまでのすべてのアルファベットを1回ずつ使って文章を作る遊び)から由来しているのでしょう。

 

<この問題の面白いところ>

 この「パングラム数」という言葉、私はプロジェクトオイラーで初めて聞いたのですが、プロジェクトオイラーでは頻出のテーマで、「難しい整数問題を解いている感じ」がしてとても好きです。まあ、「難しい」のは人力で解いた場合なんですけども!また、以前ちょっと触れたnext_permutationが大活躍する問題でもあります。オイラープロジェクトやってたからnext_permutationを知れて、その結果AGC22のAが(実質Bもオイラープロジェクトのおかげで解けている)無事解けたと考えると感慨深いです。

 

<ソースコード>

まず、軽く考察を。パンデジタル素数なので当然3や9で割れてはいけません。しかし9桁のパンデジタル数、8桁のパンデジタル数は各桁の総和がそれぞれ45,36となり、9の倍数の性質から、確実に9で割り切れてしまいます。そこで、最大のパンデジタル素数の候補となるのは7桁のパンデジタル数であるということが分かります。

あとは1~7の組み合わせを全部試せばよいのですが、ここでnext_permutationというC++の機能を使うと簡単に順列の全通りの組み合わせが列挙できます。

 なお、コード中に登場するiota(v.begin(),v.end(),val)ですが、iotaは引数valが加算可能な型であれば、v[0]にvalを入れ、その後v.end()までv[i]にval+iを格納していきます。整数や文字型などに適用可能で、順列の初期値生成に適しています。

int main() {
     ans = 0;
     vector<int> v(7);
     iota(v.begin(), v.end(), 1); //最初は1234567
     do {  //next_permutationで順列ループ
           int f = 0; //素数判定
           int a = 0;
           int dig=1;
           for (int i = 0; i < 7; ++i) {
                a += dig * v[i]; //答えを組み込む
                dig *= 10;
           }
           for (int i = 2; i*i <= a; ++i) {
                 if (a%i == 0) {
                    f = 1; break;
                 }
           }
          if (!f) ans = max(ans, a);  //最大値の更新
      } while (next_permutation(v.begin(),v.end()));

      cout << ans << endl;
      return 0;
}
 

4. Problem59 XOR decryption

f:id:tb9810w124816:20180409105149p:plain

 この問題は滅茶苦茶面白い問題なのですが、英語よわの私には文面が辛すぎるので日本語wikiをガン見しました。

f:id:tb9810w124816:20180409105308p:plain

 すると、やはり滅茶苦茶面白そうな問題であることが分かります。

 

<この問題の面白いところ>

 「この問題を通して一つの実在の暗号をプログラムを使って解いている」っていう体験がなんかハッカーっぽくてかっこいい(くそこなみかん)。でも確かにワクワクする問題だと思います。また、XORの特殊な性質「(A XOR B) XOR B =A」が利用されているのも面白いです(この性質は真理値表を書いてみると明らかになります)。XORのこの性質から、Bを鍵とすることでAを暗号化し、再び複合することもできることが明らかになります。まずA XOR B=Xで暗号化して送信者が受信者に送り、受信者は鍵Bを用いて

X XOR B =Aとすれば元の文章(平文)を解読できるというわけです。

 プロジェクトオイラーにはこの他にも暗号技術にかかわる問題が出題されており、好奇心をくすぐります。

 

<ソースコード>

 でもこの問題は自分が鍵を知らない状態で平文を解読するものです。まあ、この方がハッカー感あっていいですよね(くそこなみ)。ただし、文章は一般的な英単語から構成された文章であること、そして暗号は3桁の英語小文字であることが分かっています。3桁とはザルですね!セキュリティがなっていません()。26*26*26で全探索されてしまいます。とはいえ、26*26*26のパターンをすべて人間の目で確認して正しく解読できているか判断するのは不可能ではありませんが、大変面倒です。そこで、英文におけるアルファベットの統計的な出現率を利用します。"e"が最頻出であることは知っていましたが、それだけでは解読にやや不十分であったので、以下のurlを参考に出現率が6%を超える上位8文字の個数を文中からカウントしていき、その合計が最大となった鍵を正しい鍵とすることにしました。

文字頻度表

上位8文字はe,a,t,i,o,s,n,rでした。eat Rosin. (松やにを食え。)

さて、ソースコードは以下になります。

int main() {
        FILE *fp;
        char gomi[20];
        int pass[200000];
        fp = fopen("in.txt", "r");
        ll sc;
        int pp = 0;
        while (fscanf(fp,"%d,",&pass[pp]) != EOF) {
           pp++;
           cout << (char)pass[pp-1];
        }

       cout << endl << endl;
       int ecounter = 0;
       char id[3] = { "" };
      for (int i = 0; i < 26; ++i) {
           for (int j = 0; j < 26; ++j) {
                 for (int k = 0; k < 26; ++k) {
                        int pip = 0;
                        for (int h = 0; h < pp; ++h) {
                            

                             if(h%3==0)pass[h] ^= (i + 97);
                             else if (h % 3 == 1)pass[h] ^= (j + 97);
                             else pass[h] ^= (k + 97);
                            

                            if (pass[h] == 'e'|| pass[h] == 'a'|| pass[h] == 't'|| pass[h] == 'i'|| pass[h] == 'o'|| pass[h] == 's'|| pass[h] == 'n'|| pass[h] == 'r') pip++;
                            }
                           if (ecounter < pip) {
                                ecounter = pip;
                                id[0] = i + 97, id[1] = j + 97, id[2] = k + 97;
                           }
                           for (int h = 0; h < pp; ++h) {
                                    if (h % 3 == 0)pass[h] ^= (i + 97);
                                   else if (h % 3 == 1)pass[h] ^= (j + 97);
                                   else pass[h] ^= (k + 97);
                            }
                     }
              }
       }
      ans = 0;
      for (int h = 0; h < pp; ++h) {
             pass[h] ^= id[h%3];
             ans += pass[h];
             cout << (char)pass[h];
      }
     cout << endl << endl;
     std::fclose(fp);
     cout << ans << endl;
     return 0;
}

 一部凄まじく汚いコードですが、許してください。

 また、Disucussionでは英文にほぼ確実に含まれている"a"や"the"といった冠詞をカウントしていく手法を取っている方もいました。様々な解読法があり、非常に面白い問題だと思います。

 ちなみに解読画面はこんな感じになりました。

f:id:tb9810w124816:20180409114734p:plain

 

 5.まとめ

  今回で、一旦プロジェクトオイラーの記事は区切りをつけようと思います。この記事で興味を持ってくれた方が一人でもいれば幸いです。私はこれからも少しづつProject Eulerを進めていくつもりです。

 

yadrfu.hatenablog.com

yadrfu.hatenablog.com