計算数学演習おまけ

このおまけは、C言語でプログラムを書くことが多い研究室において、授業で扱いきれなかった有益になりそうな情報を羅列しているページです。
別に、某F研において、「研究で使うのは基本的に授業で習った範囲で大概カバーできるよ」と言ったにも関わらず、学生から「いやいや、乱数とか複数のファイル出力とか配列を関数の引数にするとか習ってねぇよ」と怒られたから、というわけではないです。たぶん。
作成中だったり、追記されていくと思いますので、必要な項目があれば催促してください。こんな項目欲しいとかあれば、メール等でご相談ください。できる限り対応します。

擬似乱数発生器の使い方（Mersenne Twister）

完全にランダムな試行によって発生した数を乱数といいます。例えば、サイコロを振って出た目は1〜6の整数値を返す乱数です。（既存の）コンピュータでは完全にランダムな試行ということは行えないので、乱数を模した疑似乱数を扱います。

C言語には標準で擬似乱数を発生させる関数randや関数srandが用意されていますが、その発生方法は線形合同法を基にしています。ちょこっとお試しで使う分には問題ないのですが、出てくる値が周期的であり、しかもその周期はそんなに長くありません。そのため、乱数を大量に使用することが多い科学計算では適切ではありません。

そこでここでは比較的周期が長く、かつ高速なMersenne Twisterを利用します。Mersenne Twisterは様々なソフトウェアでも採用されています。開発者は松本眞氏と西村拓士氏で、2人の名前の頭文字のMとTとなるようにしたとかしてないとか。

Mersenne Twosterの使い方は以下の通りです（簡易版ですので、詳しい人は自分に合ったやり方でやってください）。

  1. Mersenne Twisterのページからmt19937ar.sep.tgzをダウンロードし、解凍します（tgzの解凍が必要なので、Windowsの場合は適宜ソフトを使用）。
  2. フォルダごと作業フォルダ（今から書くであろうプログラムがある場所）に持ってきます。
  3. プログラム（ここではtest.c）を書きます。このとき注意点は以下の点
    • プロトタイプ宣言の参照：プログラムの最初の方に
      #include "mt19937ar.sep/mt19937ar.h"
      を加える（必要な関数のプロトタイプ宣言が詰まっている。< >ではなく” “でくくっていることに注意）。
    • 擬似乱数の初期化：main関数の最初の方で
      unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4; init_by_array(init, length);を追加（0x123, 0x234, 0x345, 0x456の部分を変えると乱数値も変わる。lengthは624未満にすること。）
    • 疑似乱数を使う際は
      double a = genrand_real1();

      とすれば区間[0,1]の一様乱数が生成される。
      同様に
      double a = genrand_real2();
      で区間[0,1)の一様乱数が、
      double a = genrand_real3();
      で区間(0,1)の一様乱数が生成される。
      また、
      double a = genrand_res53();
      でより長周期の区間[0,1)の一様乱数が生成される。

  4. コンパイルの際には
    gcc test.c mt19937ar.sep/mt19937ar.c -o test
    のようにmt19937ar.sep/mt19937ar.cを加えるとコンパイルできます。その他詳細はmt19937ar.sep内にあるreadme-mt.txtやmtTest.cを参照のこと。

さらにMeresenne Twisterを改良したdSFMT（Double precision SIMD-oriented Fast Mersenne Twister）もあるので、研究に使うようであればそちらの方の使用を推奨。

ファイル出力・入力

ファイル出力

授業では標準出力をする関数printfを使いました。これで問題がない場合は良いのですが、例えば授業内であった例でいうと、入力促進文もファイルに出力されてしまう、といったように、ファイル出力もしたいし標準出力もしたい、という場合があると思います。このようなときに、ファイル出力関数であるfprintfを使います。使い方は以下の通りです。

  1. 出力ファイル名を決めます（ここではoutput.datとします）。
  2. main関数の最初の方で以下のように追記します。
    FILE *fout_data = fopen("output.dat", "w");
    fopen関数の1個目引数はファイル名、2個目の引数は書き込みを表す”w”または追加読み書きを表す”a”になります。「fout_data」という変数名（正確には変数ではない）は自由に決められますので、変更した場合は、以下の部分も同じように変更してください。
    定義とfopenを分けたい場合は次のようにします。

    FILE *fout_data;
fout_data = fopen("output.dat", "w");
  3. ファイル出力部分で
    fprintf(fout_data, "%lf %lf\n", a, b)のように記載します。fprintf関数であることと1個目の引数が先ほど定義したfout_dataであること以外は、printf関数と同じような使い方です。
  4. main関数の最後の方で
    fclose(fout_data);
    を追記します。

一連の流れは2.でファイルを開き、3.で出力し、4.でファイルを閉じる、となります。4.を忘れがちなので注意が必要です（正常に出力されなかったり、メモリが解放されないなどの不具合が起きます）。

複数のファイルに出力したい場合は2.のfout_dataの部分をfout_data1やfout_data2のように複数定義し、出力先に合わせて3.のfprintf内のfout_dataを変更し、定義した分だけfcloseをすればOKです。

ファイル入力

ファイル出力と同様に、ファイル入力も可能です。この場合もscanfではなくfscanfを使います。

  1. 入力ファイル名を決めます（ここではinput.datとします）。
  2. main関数の最初の方で以下のように追記します。
    FILE *fin_data = fopen("input.dat", "r");
    fopen関数の1個目引数はファイル名、2個目の引数は読み込みを表す”r”になります。「fin_data」という変数名（正確には変数ではない）は自由に決められますので、変更した場合は、以下の部分も同じように変更してください。
  3. ファイル入力部分で
    fscanf(fin_data, "%lf %lf\n", &a, &b)
    のように記載します。fscanf関数であることと1個目の引数が先ほど定義したfin_dataであること以外は、scanf関数と同じような使い方です。
  4. main関数の最後の方で
    fclose(fin_data);
    を追記します。

補足もファイル出力と同様です。ただしファイル入力の場合は、入力するファイルが存在しない場合もありますので、そのときにエラーを検知する仕組みがあると尚良いです。詳細は下の空の型voidの後半に記載していますので、そちらも参照してください。

配列の初期化のTips

実は配列では、定義と同時のときに限り下記のような初期化ができます。

double a[3] = {0.0, 1.0, 3.0};
char b[4] = "test";

1個目は、実数配列aの要素を個別に定義するのではなく、まとめて定義しています。また2個目では、文字列の定義をしています。

これを分けて次のように定義することはできません。

double a[3];
a = {0.0, 1.0, 3.0};
char b[4];
b = "test";

文字列のTips

文字列に変数の値を参照したい場合は、sprintf関数を使い、以下のようにします。

char b[64];
sprintf(b,"a0=%lf_a1=%lf_a2=%lf", a[0], a[1], a[3]);

これで文字列bにa[0]~a[2]を反映した文字列が代入されます。sprintf関数の使い方は、1個目の引数に格納したいchar型配列を指定する以外はprintfと同じですが、結果を出力するのではなくchar型配列に格納します。

これとファイル出力を組み合わせると、次のようなこともできます。

int main() {
    double a[3] = {0.0, 1.0, 3.0};
    char fout_name[64];
    sprintf(fout_name, "a0=%lf_a1=%lf_a2=%lf.dat", a[0], a[1], a[2]);
    FILE *fout_data = fopen(fout_name, "w");
    fprintf(fout_data, "test\n");
    fclose(fout_data);
    return 0;
}

これにより、出力ファイル名に変数の値を反映させることが可能です。

空の型void

授業では関数型として変数型と同様にintやdoubleを使いました（例えばint beki(int a, int b)やdouble myabs(double a)など）。
実用上、これらの型が使えればプログラミングには何の問題もないのですが、何も返さない関数型というのもあります。例えば次のような関数は、標準出力のみ行い、戻り値は必要ありません。

void myprint(int a) {
    printf("a=%d\n", a);
}

このような場合ではvoid型で定義することができます。

また、引数についてもvoid型を使うことができます。すでに皆さんが散々main関数の定義をint main()と書いてきたと思いますは、void型であることを明確にしてint main(void)としても良いです。つまり、引数としてvoid型を使う場合は、省略可能になります。ただし「全くの同義」ではありません。関数を呼び出す際、引数にvoidが明示的に書かれている関数で引数を指定するとエラーになりますが、引数のvoidが省略されている関数で引数を指定すると、エラーにはなりません（コンパイルでwarningは出る）。このような誤用が許されるプログラムは多くのプログラミング言語では推奨されていません。voidを明示的に書く癖をつけておいても良いかもしれません。

一方で、（紹介しておいてなんですが）プログラミング言語によっては戻り型としてvoid型を推奨していないこともあります。その理由として、その関数が正常に動作したか検証が必要な場合があるからです。例えば次のようなプログラムを考えてみましょう。

#include <stdio.h>

void myscanprint();

int main(void) {
    myscanprint();
    return 0;
}

void myscanprint(void) {
    int a;
    FILE *fin_data;
    if ((fin_data = fopen("input_sample.dat", "r")) == NULL){
        /*何もしない*/
    }
    else {
        fscanf(fin_data, "%d", &a);
        fclose(fin_data);
    }
    printf("%d\n", a);
}

myscanprintの中で、ファイルが正しく開けた場合は17,18行目に進み、ファイル入力が行われるのですが、ファイルが正しく開けなかった場合、fin_dataにはNULLが入ります（NULL型については、授業ではポインタを教えていないので、ざっくり言うと、fin_dataが指し示す先には何もない（何も示していない）状態）。その結果、aは初期化されておらず、適当な数値が入った状態で出力されます。
しかし、実行した側から見ると、このプログラムが正常にファイル入力を行ったかどうかはわかりません。こういったときに、次の「プログラムが正しく動いているかどうかを判断するための戻り値」があると便利です。

#include <stdio.h>

int myscanprint();

int main(void) {
    int errorcheck;
    errorcheck = myscanprint();
    if (errorcheck != 0) {
        printf("warning: ファイル入力は正しく行われていないかもしれません。\n");
    }
    return 0;
}

int myscanprint(void) {
    int a, errorcode;
    FILE *fin_data;
    if ((fin_data = fopen("input_sample.dat", "r")) == NULL){
        errorcode = 1;
    }
    else {
        fscanf(fin_data, "%d", &a);
        fclose(fin_data);
        errorcode = 0;
    }
    printf("%d\n", a);

    return errorcode;
}

すでにお気づきの人もいるかもしれませんが、main関数の最後に何故return 0;を書くのか… それは下のmain関数の引数と戻り値で説明します。

main関数の引数と戻り値

色々な関数を定義してきたと思いますが、関数といえば戻り値と引数です。ところでmain関数も関数です。ということはmain関数の引数と戻り値についても考えてみましょう。

そもそもとして、関数は他の関数から呼び出されることが前提で作ります。では、main関数を他の関数から呼び出すことはあるのでしょうか？また、他の関数に値を返すことがあるのでしょうか？main関数とはプログラム実行時に最初に実行される関数ですので、他の関数に呼び出されることはありません。main関数が呼び出されるのは、そのプログラムの外、すなわちターミナル上や別のブログラムからです。

main関数の引数

利点

まずはmain関数で引数が指定できると、どのような利点があるのか説明します。
1つ目は、コンパイルの手間が省けることです。main関数の引数は実行時に指定することができます。つまり、「変数の中身を変えるたびにコンパイルする」という作業がなくなるわけです。これによってエディタとターミナルを行ったり来たりすることなく、変数の値を変更することができます。
2つ目に（かなりアドバンスドな内容に踏み込みますが）、シェルスクリプトと組み合わせるとパラメータのグリッドサーチ（バカパラとも呼ばれる）が簡単にできるようになる、ということです。
これらについては後ほどの具体例を見てもらえればと思います。

使い方

main関数の引数については、int main(int argc, char **argv)の形式がよく使われます。argcは引数の数で、argvに2次元配列の形式で文字列が格納されます（ちなみにargはargument（引数）に由来）。コンパイラによってはさらに進化していたりしますが、基本的にはこの形だけ覚えておけば十分です（コンパイラ独自の形とかにすると環境が変わった際に不便なので…）。
次のサンプルプログラムarg_main_sample.cを考えましょう。

arg_main_sample.c

#include <stdio.h>

int main(int argc, char **argv) {
    int i;
    for (i = 0; i < argc; i++) {
        printf("%s\n", argv[i]);
    }
    return 0;
}

（注：%sは文字”列”（つまり配列形式）の出力時のフォーマット）
このプログラムを

gcc arg_main_sample.c -o arg_main_sample
./arg_main_sample

のように、ただただコンパイルして実行すると

./arg_main_sample

と表示されます。次にそのまま

./arg_main_sample this is a sample.

として実行してみてください。結果は、

./arg_main_sample
this
is
a
sample.

のようになったと思います。この結果から、「ターミナルで入力した内容がプログラムに反映された」という事実は少なくともわかると思います。

今回のプログラムの実行時の引数と、argcとargvの関係について図にまとめます。

実行時には実行プログラム名を含む空白で区切られた文字列がそれぞれ配列に格納され、それらをさらに配列として組み合わせているため、2次元配列の形になっています。arg_main_sample.cでは、それらの配列を文字列として表示しています。

引数が文字列の場合は良いですが、引数が数値の場合は一手間が必要です。例えば先ほどのarg_main_sample.cを数値の計算ができるように、次のように変えてみましょう。

#include <stdio.h>

int main(int argc, char **argv) {
    double x[argc], sum = 0.0;
    int i; for (i = 0; i < argc; i++) {
        x[i] = argv[i];
        sum += x[i];
    }
    printf("%lf\n", sum);
    return 0;
}

入力した数値を全て足し合わせるプログラムになります。ところが、このままコンパイルすると、次のようなエラーが出ます（コンパイラによって多少異なります）。

arg_main_sample2.c:6:8: error: assigning to 'double' from incompatible type 'char *'
    6 |                 x[i] = argv[i];
      |                      ^ ~~~~~~~

要約すると、「x[i]はdoubleなのに、argv[i]は文字列（文字の配列）なのはけしからん」と怒られているわけです。つまり、数値を入れたつもりが、プログラムは文字列として受け取っているわけです。そこで関数atofを使って次のように書き換えます。

arg_main_sample2.c

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv) {
    double x[argc], sum = 0.0;
    int i; for (i = 0; i < argc; i++) {
        x[i] = atof(argv[i]);
        sum += x[i];
    }
    printf("%lf\n", sum);
    return 0;
}

注意点として、関数atofは通常は定義されていませんので、2行目にあるようにstdlib.hをインクルードする必要があります。また、今回はx[i]がdouble型ですので、関数atofを使っていますが、int型の場合は関数atoiを使ってください。

試しにarg_main_sample2.cをコンパイルし、

./arg_main_sample2 1 15 5e-2

と打ってみてください。16.050000と表示されます。
（”5e-2″は5×10-2のこと。便利なので覚えておきましょう。）

main関数の戻り値

上の空の型voidで説明した通り、戻り値が必要ない関数でも、その関数が正しく終了したかわかるようにint型にした方が便利です。このことはmain関数についても当てはまります。
本来main関数は必ず実行される関数で、他の関数への戻り値が必要ない関数です。しかし、プログラムが正常に終了したのかどうか、どこか異常なことが起きていないか、という判断ができる方が便利でしょう。そのため、（少なくともこの授業では）main関数はvoid型ではなくint型を使ってします（もちろん、int型だけでなく、何らかの計算した数値を返すようにdouble型にすることも可能です）。基本的には、正常終了の場合は0を、異常終了の場合はエラーの種類に応じて整数値を返すようにプログラムを構築すると良いでしょう。

ターミナル上（sh, bash, csh, zshなど）でmain関数の戻り値を確認するには、プログラムを実行後に

echo $?

と打つと確認できます（ちなみにfishの場合はecho $status（備忘録））。詳しくは述べませんが、echoはターミナル上での”printf”に相当するようなものと考えて下さい。また、$?は直前の実行結果が代入されている変数です。詳しいことは「シェルスクリプト」「bash」「zsh」などで検索すると、プログラミングでの手間を楽にしてくれる方法が見つかるかもしれません。

関数の引数に配列を入れたい・返したい

関数の引数に配列を入れたいときもあると思います。例えば、配列の最大値を返す関数とか、並べ替えて配列を返して欲しいとか。関数と配列が絡むと、結構面倒なことが多いです。本当はメモリアクセスとかポインタとアドレスとか色々習った上でやるべきなのですが、そうも言っていられないので、やり方だけ簡潔に書いておきます。

引数にする場合

まず配列を引数とする場合には、以下のようになります。

#include <stdio.h>

double f(double *, int); // プロトタイプ宣言

int main(void) {
    double x[] = {1,2,3,4,5,6,7,8,9,10};
    printf("%lf\n", f(x, 10));
    return 0;
}

// 実際の関数の定義
double f(double *x, int N) {
    int i;
    double mean = 0.0;
    for (i = 0; i < N; i++) {
        mean += x[i];
    }
    mean /= N;
    return mean;
}

通常のdouble型の変数であればプロトタイプ宣言時の仮引数はdoubleだけですが、*もついています。また、配列の要素数も引数に入れておくことをお勧めします。実際の関数の定義部分でも、同様です。
使う際には、7行目にあるように、[ ]などは付けずに引数にします。配列の要素がint型やchar型の場合も適宜変数型を変更すれば同様に使えます。

注意点：仮引数とした配列を変更すると、もともとの関数の配列も変更されるので注意！

先程のプログラムを改変し、次のようなプログラムを見てみましょう。

double f(double *, int); // プロトタイプ宣言

int main(void) {
    int i;
    double x[] = {1,2,3,4,5,6,7,8,9,10};

    for (i = 0; i < 10; i++) {
        printf("%lf\n", x[i]);
    }
    
    printf("average = %lf\n", f(x, 10));
    
    for (i = 0; i < 10; i++) {
        printf("%lf\n", x[i]);
    }
    
    return 0;
}

// 実際の関数の定義
double f(double *x, int N) {
    int i;
    double mean = 0.0;
    for (i = 0; i < N; i++) {
        x[i] *= 2.0;
        mean += x[i];
    }
    mean /= N;
    return mean;
}

大きな違いは、関数fの中で和を計算する前にx[i] *= 2.0;としているところです。また、結果の出力の前後で配列の値も表示しています。結果を見てみると、

1.000000
2.000000
3.000000
4.000000
5.000000
6.000000
7.000000
8.000000
9.000000
10.000000
average = 11.000000
2.000000
4.000000
6.000000
8.000000
10.000000
12.000000
14.000000
16.000000
18.000000
20.000000

となっています。不思議なことに、関数fを実行したあとでmain関数内での配列xまでそれぞれ2倍されてしまっています。このあたりのことを正しく理解しようとすると、ポインタの概念が必要になりますが、ざっくりとイラストを使って説明します。

引数が変数の場合（図上側）は、関数を超えて変数を持ち込むのではなく、関数に入る前に値をコピーしていると考えてください。したがって、呼び出し元の関数の変数（この場合main関数の変数x）の値が呼び出し先の関数の仮引数（この場合関数fの変数a）にコピーされますが、呼び出し先の関数の仮引数を変更しても呼び出し元の関数の変数（この場合main関数の変数x）には影響を与えません。
一方、引数が配列の場合（図下側）は、変数の場合とは異なり、配列の（先頭の）データが置いてある「メモリ上の場所」が受け渡されます。そのため、呼び出し先の関数の仮引数にも「メモリ上の場所」が入ることになるのですが、これは関数をまたがって共通のものになります。そのため、呼び出し先の関数内で「メモリ上の場所」を使って配列の要素の値を変更すると、呼び出し元の関数の配列も変更されてしまいます。

このような事態は、呼び出し先の関数で、仮引数の配列を別の配列に格納することで避けられます（下のサンプル参照）。

#include <stdio.h>

double f(double *, int); // プロトタイプ宣言

int main(void) {
    int i;
    double x[] = {1,2,3,4,5,6,7,8,9,10};

    for (i = 0; i < 10; i++) {
        printf("%lf\n", x[i]);
    }
    
    printf("average = %lf\n", f(x, 10));
    
    for (i = 0; i < 10; i++) {
        printf("%lf\n", x[i]);
    }
    
    return 0;
}

// 実際の関数の定義
double f(double *x, int N) {
    int i;
    double mean = 0.0;
    double y[N];
    
    for (i = 0; i < N; i++) {
        y[i] = x[i];
    }
    for (i = 0; i < N; i++) {
        y[i] *= 2.0;
        mean += y[i];
    }
    mean /= N;
    return mean;
}

分岐が多い場合にはswitch-case文

分岐といえばif文ですが、計算結果などの値に応じて分岐が複数ある場合にはswitch文が見やすくなります。例えば、次の2つの構文は同じになります。

switch (式) {
    case 0:
        式==0のときの処理;
        break;

    case 1:
        式==1のときの処理;
        break;
    .
    .
    .
    default:
        いずれも該当しなかった場合の処理;
    break;
}
if (式 == 0) {
    式==0のときの処理;
}
else if (式 == 1) {
    式==1のときの処理;
}
.
.
.
else {
    いずれも該当しなかった場合の処理;
}

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University