コンピュータ「オイラと数学しよう」 ――オイラープロジェクトの問題にチャレンジ（その2）：無償ソフトで技術計算しよう【プログラミング応用編】（3）（2/3 ページ）

「プログラミング応用編」では、プログラミングコンテストの問題などをFreeMatでプログラム作成していく。今回も「オイラープロジェクト」の問題を解こう。

[伊藤孝宏，MONOist]

問題6

Problem 6　 Sum square difference

The sum of the squares of the first ten natural numbers is, 1² + 2² + ... + 10² = 385

The square of the sum of the first ten natural numbers is, (1 + 2 + ... + 10)² = 55² = 3025

Hence the difference between the sum of the squares of the first ten natural numbers and the square of the sum is 3025-385 = 2640.

Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

　意味は、「1〜10までの自然数の二乗和は385になります。一方、和の二乗は3025になります。従って、和の二乗と二乗和との差は2640となります。1〜100までの自然数の和の二乗と二乗和との差を求めなさい」。

　1〜10までの配列nを作り、和の二乗は(sum(n))^2となり、二乗和はドット演算子を使い、sum(n.^2)で求まります。答えは下記のようになります。これもサービス問題ですね。

--> n=1:10;(sum(n))^2-sum(n.^2)
ans = 2640
--> n=1:100;(sum(n))^2-sum(n.^2)
ans = 25164150

問題7

Problem 7　10001st prime

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.

What is the 10 001st prime number？

　意味は、「最初から6番目までの素数は、2、3、5、7、11、13となり、6番目の素数は13です。では、10001番目の素数を求めなさい」。

　エラトステネスの篩を使います。今回の場合、n番目の素数と、あらかじめ素数の数が分かっているので、ps=[2,3,5,7,11,zeros(1,n-5)]として、素数の配列psを用意します。既に分かっている素数2、3、5、7、11の後は0がnまで続く配列で、0の部分に素数が格納されていきます。nnを11から増加させ、psで割り切れるかどうかで素数であるかを調べます。

　ここで、偶数はあきらかに素数にはならないので、nnは11から2ずつ増加させます。素数かどうかの判別は、find(mod(nn,ps)==0)として、割り切れる、つまりmodが0となる位置を調べます。mod(nn,0)はnnとなるので、まだ素数が格納されていない部分でmodが0となることはありません。

　modが0となる要素が見つからない場合、空の配列が返るので、「isempty」コマンドで検出します。isemptyコマンドは、引数が空配列だと1を返し、それ以外では0を返します。modに0がない、すなわち、素数の場合は、psに追加し、cntを1増加させます。cntがnになるまで、以上の操作を繰り返すと、n個の素数が得られ、目的の素数は素数の配列psのn番目の要素になります。

function p7=ex331(n)
    ps=[2,3,5,7,11,zeros(1,n-5)];
    cnt=5;nn=11;
    while(cnt~=n)
        if(isempty(find(mod(nn,ps)==0)))
            cnt=cnt+1;ps(cnt)=nn;
        end
        nn=nn+2;
    end
    p7=ps(n);

ex331.m

＞＞「ex331.m」ダウンロード

　例で試してみると、下記のように、正しいことが分かります。

--> ex331(6)
ans = 13

　本番で試してみると、残念ながら、FreeMat自体の処理速度が遅いせいなのか、答えが得られませんでした。そこで、また掟（おきて）破りのようで申し訳ないのですが、MATLABでex331.mを動作させてみました。答えは以下が得られました。

>> ex331(10001)
ans = 104743