連載
» 2014年10月21日 10時00分 UPDATE

無償ソフトで技術計算しよう【プログラミング応用編】(1):FreeMatで魔方陣を作ってみよう (1/2)

「プログラミング応用編」では、プログラミングコンテストの問題などをFreeMatでプログラム作成していく。今回は、魔方陣を作ってみよう。

[伊藤孝宏,MONOist]

 将棋のルールを習っても、将棋ができるとは限りません。同じように、プログラムを作成できるようになるには、コマンドの使い方の他に、実際にプログラムを作成してみることが必要です。今回から応用編として3回に渡り、プログラミングコンテストの問題などを例に実際にプログラムを作成してみます。

 筆者もプログラミングは得意ではありませんが、一緒にいくつかのプログラムを作成していきましょう。手始めに、MATLABにはあるけど、FreeMatにはない関数を作成してみます。

筆者注:FreeMatはコマンド入力後にエンターキーを押すとコマンドを実行します。本連載ではエンターキーの記述を省略しますが、操作の際にはコマンド入力後にエンターキーが必要です。



FreeMatでエンジニアリング計算:
入門編
グラフィックス編
プログラミング基礎編
プログラミング処理編

最大公約数、最小公倍数を求める

 MATLABは2つの数値の最大公約数を求めるコマンド「gcd」を持っています。これを関数mファイルで作成してみましょう。

 最大公約数は複数の公約数のうち最大のもので、例えば、30と24の公約数は1,2,3,6ですから、最大公約数は6です。一方、最小公倍数は公倍数のうち最小のものです。整数aと整数bの最大公約数gcdと最小公倍数lcmとの間には、a×b=gcd×lcmの関係があるので、最大公約数が求まれば、lcm=a×b/gcdにより、最小公倍数が求められます。例えば、30と24の最小公倍数は、30*24/6=120です。

 コンピュータで最大公約数を求めるには、ユークリッドの互除法というアルゴリズムを用います。ユークリッドの互除法は、2つの整数m、n(m>n)において、mとnの最大公約数は、m-nとnとの最大公約数に等しいという性質を用いて、お互いに引き算を行い、最終的に両者が等しくなったときの値が最大公約数となります。具体的には、mとnとが等しくない間は「m>nの場合、m=m-n。それ以外は、n=n-m」を繰り返し、m=nとなったら「mおよびnは最大公約数に等しい」という動作にします。

 以上をFreeMatでプログラミングすると、ex320.mとなります。

function gl=ex320(a,b)
    m=a;n=b;
    while(m~=n)
        if(m>n)
            m=m-n;
        else
            n=n-m;
        end
    end
    gl=[m,a*b/m];
ex320.m

>>「ex320.m」ダウンロード

 ex320.mでは、a*b/mとして、最小公倍数を求めて、最大公約数、最小公倍数の2列の配列として返します。

 例えば、128と72の最大公約数を求めるには、コマンドウィンドウで下記のように入力します。すると、最大公約数は8、最小公倍数は1152という結果が得られます。

--> ex320(128,72)
ans = 8  1152

 「ユークリッドの互除法」は余りを用いても同じような結果が得られます。つまり、mをnで割った余りをkとして、kが0でない間は以下を繰り返します。m=n、n=kとして、kが0となったら、mおよびnは最大公約数に等しくなります。

 以上をプログラムにすると、ex321.mとなります。ここで、mod(m,n)は、mをnで割った余りを返します。

function gl=ex321(a,b)
    m=a;n=b;k=1;
    while(k~=0)
        k=mod(m,n);
        m=n;n=k;
    end
    gl=[m,a*b/m];
ex321.m

>>「ex321.m」ダウンロード

 差を用いる方法よりも、剰余を用いる方法の方が計算効率が良くなります。tic、tocコマンドを用いて、ex320.mとex321.mの計算に要する時間を計測してみると、下記のように、ex321.mの方が計算効率が良いことが分かります。

--> tic;ex320(128,72);time=toc
time = 0.1480
--> tic;ex321(128,72);time=toc
time = 0.0040

素数を求める

 素数は、「1と自分自身以外に約数を持たない数」で、2、3、5、7、11などです。MATLABは素数を求めるコマンド「primes」を持っていますが、FreeMatにはありません。そこで、任意の整数nまでの素数を求める関数mファイルを作成してみます。

 素数を求めるアルゴリズムとして、「エラトステネスの篩(ふるい)」を用います。

 エラトステネスの篩では、以下の手順で素数を探索します。

 「リストに2からnまでの整数を用意する。リストの先頭の数を素数リストに移動し、その倍数をリストから除外する。以上の操作をリストの先頭の値がnの平方根に達するまで行う」

 上の動作をプログラムにしたのが、ex322.mです。

function p=ex322(n)
    pr=ones(1,n);pr(1)=0;
    for k=2:round(n^0.5)
        if(pr(k)==1)
            for m=2*k:n
                if(mod(m,k)==0)
                    pr(m)=0;
                end
            end
        end
    end
    p=find(pr);
ex322.m

>>「ex322.m」ダウンロード

 始めに、pr=ones(1,n)として1行n列の要素1の配列を作成します。prは要素番号が篩に残っているかどうかを表し、要素が1の番号が篩に残り、0であれば除外されます。番号1は素数ではないので、pr(1)=0として除外します。pr(k)が1であれば、mod(m,k)が0となるkは、pr(k)=0として、篩から除外します。以上をn0.5の整数分まで繰り返します。すると、番号が素数では要素が1となる配列prが得られるので、find(pr)として、要素が1の部分の番号だけを取り出すと、素数が得られます。

 例えば、10までの素数を求めるには、コマンドウィンドウで下記のように入力すると、2、3、5、7が得られます。

--> ex322(10)
ans = 2 3 5 7

素因数分解を行う

 正の整数を素数の積に分解することを素因数分解と呼びます。例えば、42=2*3*7と素因数分解できます。MATLABはfactorという素因数分解コマンドを持っています。

 素因数分解は以下のアルゴリズムを用います。nを2で割り切れなくなるまで割る。次に、3で同様のことを行います。以降、4、5、6……と同様のことを行います。ここで、素数以外の数、例えば、4では、既に2で割り切れているので、4で割り切れることはなく、4は素因数にはなりません。こうして、nが1になるまで以上を繰り返すと、割っていった数が素因数となります。

 以上をプログラムにしたのが、ex323.mです。nがaで割り切れた場合(剰余が0)、素因数の配列fにaを追加し、新たなnとして、nをn/aで置き換えます。割り切れない場合は、aをa+1で置き換えます。

function f=ex323(n)
    a=2;f=[];
    while(n~=1)
        if(mod(n,a)==0)
            f=[f,a];
            n=n/a;
        else
            a=a+1;
        end
    end
ex323.m

>>「ex323.m」ダウンロード

 例えば、126の素因数分解は、下記のようにコマンドウィンドウに入力すると、2、3、3、7に分解されることが分かります。

--> ex323(126)
ans = 2 3 3 7
       1|2 次のページへ

Copyright© 2017 ITmedia, Inc. All Rights Reserved.