C言語でrand関数と自作関数で乱数の生成方法を教えて！
こういった悩みにお答えします．
本記事の信頼性
- リアルタイムシステムの研究歴12年．
- 東大教員の時に，英語でOS（Linuxカーネル）の授業．
- 2012年9月～2013年8月にアメリカのノースカロライナ大学チャペルヒル校（UNC）コンピュータサイエンス学部で客員研究員として勤務．C言語でリアルタイムLinuxの研究開発．
- プログラミング歴15年以上，習得している言語: C/C++，Python，Solidity/Vyper，Java，Ruby，Go，Rust，D，HTML/CSS/JS/PHP，MATLAB，Assembler (x64，ARM)．
- 東大教員の時に，C++言語で開発した「LLVMコンパイラの拡張」，C言語で開発した独自のリアルタイムOS「Mcube Kernel」をGitHubにオープンソースとして公開．
- 2020年1月～現在はアメリカのノースカロライナ州チャペルヒルにあるGuarantee Happiness LLCのCTOとしてECサイト開発やWeb/SNSマーケティングの業務．2022年6月～現在はアメリカのノースカロライナ州チャペルヒルにあるJapanese Tar Heel, Inc.のCEO兼CTO．
- 最近は自然言語処理AIとイーサリアムに関する有益な情報発信に従事．
- （AI全般を含む）自然言語処理AIの論文の日本語訳や，AIチャットボット（ChatGPT，Auto-GPT，Gemini（旧Bard）など）の記事を50本以上執筆．アメリカのサンフランシスコ（広義のシリコンバレー）の会社でプロンプトエンジニア・マネージャー・Quality Assurance（QA）の業務委託の経験あり．
- （スマートコントラクトのプログラミングを含む）イーサリアムや仮想通貨全般の記事を200本以上執筆．イギリスのロンドンの会社で仮想通貨の英語の記事を日本語に翻訳する業務委託の経験あり．
こういった私から学べます．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
乱数
乱数はランダムな数を意味します．
乱数の特徴として，偏りがないことや再現性がないことが要求されますが，コンピュータの計算による乱数生成では次の計算が予想できてしまいます．
IoT等でよく使われる暗号計算でセキュリティを高めるためには，予測性が低い（安全性が高い）乱数生成が必要とされています．
一様乱数と一様分布
本記事で紹介する乱数は，一様乱数で生成されるとします．
一様乱数とは，確率変数 の値の如何に関わらず確率密度関数が一定の値をとるような分布（一様分布）で生成される乱数のことです．
一様分布は，以下の数式になります．
$$p(x) = 1 \ (0 < x < 1)$$
rand関数による乱数の生成
|
1 2 3 |
int rand(void); int rand_r(unsigned int *seedp); void srand(unsigned int seed); |
rand関数は，[0, RAND_MAX]の範囲の擬似乱数整数を返します．
RAND_MAXは私の環境のstdlib.hで以下のように定義されています．
|
1 2 |
/* The largest number rand will return (same as INT_MAX). */ #define RAND_MAX 2147483647 |
rand_r関数は，rand関数をリエントラント可能（再入可能）に改良したものです．
リエントラント可能とは，複数のスレッドが同時にアクセスしても正常に動作することをいいます．
rand_r関数の引数seedpは，rand_r関数の呼び出し間で状態を保持するために使用されるunsigned int型のポインタです．
seedpが指す整数と同じ初期値でrand_r関数を呼び出し，呼び出し間でその値が変更されなければ，同じ疑似乱数系列を生成します．
srand関数は，疑似乱数整数の新しい種を引数に設定します．
この種を基に疑似乱数整数を生成します．
同じ種を設定すると同じ疑似乱数整数を生成します．
rand関数による乱数の生成コード
rand関数による乱数の生成コードは以下になります．
rand関数の返り値\(x\)は\(0 \leq x \leq RAND\_MAX\)なので，12行目の「(rand() + 1.0) / (RAND_MAX + 2.0)」で規格化（正規化）します．
規格化すると，get_uniform関数の返り値\(x\)が\(0 < x < 1\)になります．
ここで，get_uniform関数の返り値\(x\)に0と1は含まれないことに注意して下さい．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #define NR_LOOPS 5 double get_uniform(void) { return (rand() + 1.0) / (RAND_MAX + 2.0); } int main(void) { int i; for (i = 0; i < NR_LOOPS; i++) { printf("get_uniform() = %lf\n", get_uniform()); } return 0; } |
実行結果は以下になります．
|
1 2 3 4 5 6 7 |
$ gcc rand.c $ a.out get_uniform() = 0.840188 get_uniform() = 0.394383 get_uniform() = 0.783099 get_uniform() = 0.798440 get_uniform() = 0.911647 |
rand/srand関数による乱数の生成コード
rand関数では，上記のように同じ疑似乱数整数しか生成できません．
そこで，srand関数で疑似乱数整数の新しい種を引数に設定することで，異なる擬似乱数整数を生成できます．
rand/srand関数による乱数の生成コードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #define NR_LOOPS 5 double get_uniform(void) { return (rand() + 1.0) / (RAND_MAX + 2.0); } int main(void) { int i; srand(100); for (i = 0; i < NR_LOOPS; i++) { printf("get_uniform() = %lf\n", get_uniform()); } return 0; } |
実行結果は以下になります．
rand関数のみを利用した場合と異なる擬似乱数整数を生成していることがわかります．
|
1 2 3 4 5 6 7 |
$ gcc srand.c $ a.out get_uniform() = 0.315598 get_uniform() = 0.284943 get_uniform() = 0.240601 get_uniform() = 0.484127 get_uniform() = 0.375793 |
rand/srand関数による時刻で変動する乱数の生成コード
rand/srand関数を利用したコードにより生成する擬似乱数整数は変わりましたが，種の値が100と一定なので同じ乱数しか生成できません．
そこで，種の値を時刻に設定することで，実行毎に生成する乱数を変更します．
rand/srand関数による時刻で変動する乱数の生成コードは以下になります．
15行目でtime関数を呼び出して，その返り値をsrand関数の引数に設定します．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <time.h> #define NR_LOOPS 5 double get_uniform(void) { return (rand() + 1.0) / (RAND_MAX + 2.0); } int main(void) { int i; srand((unsigned int) time(NULL)); for (i = 0; i < NR_LOOPS; i++) { printf("get_uniform() = %lf\n", get_uniform()); } return 0; } |
実行結果の例は以下になります．
実行毎に生成する乱数が変動していることがわかります．
※time関数の返り値は秒なので，実行毎に1秒以上の間隔を空けましょう．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
$ gcc srand_time.c $ a.out get_uniform() = 0.873856 get_uniform() = 0.951136 get_uniform() = 0.521113 get_uniform() = 0.404735 get_uniform() = 0.028890 $ a.out get_uniform() = 0.590964 get_uniform() = 0.278478 get_uniform() = 0.128135 get_uniform() = 0.045421 get_uniform() = 0.898250 |
drand48_r関数による乱数の生成コード
|
1 |
int drand48_r(struct drand48_data *buffer, double *result); |
rand_r関数のseedp引き数が指す値により提供される状態は非常に小さな空間なので，弱い疑似乱数生成器（予測性が高い擬似乱数生成器）になってしまいます．
drand48_r関数は，rand_r関数より大きな空間で擬似乱数を生成します．
rand_r関数の空間サイズは32ビット（または16ビット）なのに対して，drand48_r関数の空間サイズは48ビットになります．
ただし，drand48_r関数はGCC拡張なので，Visual Studioでは使えないことに注意して下さい．
drand48_r関数による乱数の生成コードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <time.h> #define NR_LOOPS 5 struct drand48_data rand_buf; double get_uniform(void) { double d; drand48_r(&rand_buf, &d); return d; } int main(void) { int i; srand48_r(time(NULL), &rand_buf); for (i = 0; i < NR_LOOPS; i++) { printf("get_uniform() = %lf\n", get_uniform()); } return 0; } |
実行結果は以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
$ gcc drand48_r.c $ a.out get_uniform() = 0.057465 get_uniform() = 0.333993 get_uniform() = 0.884453 get_uniform() = 0.010115 get_uniform() = 0.751871 $ a.out get_uniform() = 0.928268 get_uniform() = 0.038584 get_uniform() = 0.622898 get_uniform() = 0.475636 get_uniform() = 0.740057 |
自作関数による乱数の生成
自作関数による乱数の生成方法として以下の2つを紹介します．
- 線形合同法による乱数生成
- メルセンヌ・ツイスタによる乱数生成
線形合同法による乱数生成
線形合同法は，以下の漸化式で乱数を生成する擬似乱数生成器です．
$$X_{n+1} = (A * X_n + B)\ \% \ M$$
ここで，\(A\)，\(B\)，\(M\)は定数で，\(M>A\)，\(M>B\)，\(A>0\)，\(B \geq 0\)です．
線形合同法による乱数の生成コードは以下になります．
ここで，定数\(A=48271\)，\(B=0\)，\(M=INT\_MAX (=2147483647)\)と設定します（ステファン・パークとキース・ミラーの設定値）．
\((x + 1) / (M + 1)\)と規格化することで一様乱数を生成できます．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <limits.h> #include <time.h> #define A 48271 #define B 0 #define M INT_MAX #define NR_LOOPS 5 static unsigned long x; double get_uniform(void) { x = (A * x + B) % M; return (double)(x + 1.0) / (M + 1.0); } int main(void) { int i; x = (unsigned long) time(NULL); for (i = 0; i < NR_LOOPS; i++) { printf("get_uniform() = %lf\n", get_uniform()); } return 0; } |
実行結果は以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
$ gcc lcg.c $ a.out get_uniform() = 0.310967 get_uniform() = 0.700302 get_uniform() = 0.263659 get_uniform() = 0.097595 get_uniform() = 0.025294 $ a.out get_uniform() = 0.310990 get_uniform() = 0.785334 get_uniform() = 0.865012 get_uniform() = 0.989207 get_uniform() = 0.993250 |
メルセンヌ・ツイスタによる乱数生成
メルセンヌ・ツイスタとは，線形合同法の多次元で規則的に分布（不均等分布）するという欠点を解決した疑似乱数生成器です．
不均等分布すると連続する値間の相関性が高くなり，統計的にランダムになるため予測性が低くなります．
メルセンヌ・ツイスタによる乱数の生成コードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <time.h> #include "SFMT.h" #define NR_LOOPS 5 sfmt_t sfmt; double get_uniform(void) { return sfmt_genrand_real3(&sfmt); } int main(void) { int i; sfmt_init_gen_rand(&sfmt, (unsigned int) time(NULL)); for (i = 0; i < NR_LOOPS; i++) { printf("get_uniform() = %lf\n", get_uniform()); } return 0; } |
上記のコードをビルドするためには，メルセンヌ・ツイスタのライブラリから必要なコードのみを抽出したものをこちらからダウンロードして下さい．
mt.zipファイルはunzipコマンドで解凍します．
|
1 |
$ unzip mt.zip |
mt.zipファイルの中身は以下になります．
- SFMT-common.h, SFMT-params.h, SFMT-params19937.h, SFMT.c, SFMT.h：メルセンヌ・ツイスタのライブラリ
- mt.c：メルセンヌ・ツイスタのテストプログラム
- Makefile：上記のコードのビルド用
※メルセンヌ・ツイスタのライブラリはとても長いコードなので説明は省略します．
上記のコードをビルドするMakefileは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 |
# # Mersenne Twister # # Author: Hiroyuki Chishiro # SRCS = SFMT.c mt.c OBJS = $(SRCS:.c=.o) TARGET = mt CC = gcc RM = rm -f CFLAGS += -Wall CFLAGS += -DSFMT_MEXP=19937 .PHONY: all clean all: $(TARGET) $(TARGET): $(OBJS) $(CC) -o $@ $(OBJS) $(CFLAGS) $(LDFLAGS) clean: $(RM) -r $(TARGET) $(OBJS) *~ |
mt.zipファイルを解凍した時に作成したmtディレクトリに移動し，makeでビルドして実行すると，以下のようになります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
$ cd mt $ make gcc -Wall -DSFMT_MEXP=19937 -c -o SFMT.o SFMT.c gcc -Wall -DSFMT_MEXP=19937 -c -o mt.o mt.c gcc -o mt SFMT.o mt.o -Wall -DSFMT_MEXP=19937 $ mt get_uniform() = 0.921015 get_uniform() = 0.259150 get_uniform() = 0.199722 get_uniform() = 0.140602 get_uniform() = 0.239370 $ mt get_uniform() = 0.618190 get_uniform() = 0.838736 get_uniform() = 0.827866 get_uniform() = 0.426207 get_uniform() = 0.526618 |
モンテカルロ法で円周率の計算
モンテカルロ法とは，乱数を利用した試行を繰り返すことにより近似解を求める手法です．
メルセンヌ・ツイスタを利用してモンテカルロ法で円周率を求めるコード一式をこちらからダウンロードして同様に展開して下さい．
|
1 |
$ unzip mt_monte_carlo.zip |
mt.zipとの主な違いはmt.cがmt_monte_carlo.cに変更したことで，そのコードは以下になります．
37行目で「(double) count / NR_LOOPS * 4」と計算する理由は，モンテカルロ法で円周率を求めるの記事がわかりやすいので参考にして下さい．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <time.h> #include "SFMT.h" #define NR_LOOPS 100000000 sfmt_t sfmt; double get_uniform(void) { return sfmt_genrand_real3(&sfmt); } int main(void) { int i; int count = 0; double x, y, z; sfmt_init_gen_rand(&sfmt, (unsigned int) time(NULL)); for (i = 0; i < NR_LOOPS; i++) { x = get_uniform(); y = get_uniform(); z = x * x + y * y; if (z <= 1.0) { count++; } } printf("pi = %lf\n", (double) count / NR_LOOPS * 4); return 0; } |
実行結果は以下になります．
3回実行したところ，円周率の値は3.141592の近似解になっていることがわかります．
|
1 2 3 4 5 6 7 8 9 10 11 |
$ cd mt_monte_carlo $ make gcc -Wall -DSFMT_MEXP=19937 -c -o SFMT.o SFMT.c gcc -Wall -DSFMT_MEXP=19937 -c -o mt_monte_carlo.o mt_monte_carlo.c gcc -o mt_monte_carlo SFMT.o mt_monte_carlo.o -Wall -DSFMT_MEXP=19937 $ mt_monte_carlo pi = 3.141621 $ mt_monte_carlo pi = 3.141517 $ mt_monte_carlo pi = 3.141241 |
getrandomシステムコールで乱数の取得
|
1 |
ssize_t getrandom(void *buf, size_t buflen, unsigned int flags); |
getrandomシステムコールは，bufが指すバッファに最大buflenのランダムなバイトを詰めます．
これらのバイトは，ユーザ空間の乱数生成器のシードや暗号演算の目的で利用できます。
getrandomシステムコールは，エントロピープール（Linuxが保持するバッファで複数の環境ノイズ）から乱数を取得します．
エントロピープールが保持する複数の環境ノイズには，入力デバイスやストレージで発生する割り込みを元にした情報が格納されます．
エントロピープールの利用可能なエントロピー数を知りたい場合は，以下のコマンドを実行して下さい．
|
1 2 |
$ cat /proc/sys/kernel/random/entropy_avail 2403 |
エントロピープールの上限値を知りたい場合は，以下のコマンドを実行して下さい．
|
1 2 |
$ cat /proc/sys/kernel/random/poolsize 4096 |
デフォルトでは，getrandomシステムコールは/dev/urandomからエントロピーを引き出します．
flags引数にGRND_RANDOMを設定すると/dev/randomに変更できます．
また，flags引数にGRND_NONBLOCKを設定するとノンブロッキング状態になります．
ノンブロッキング状態では，/dev/randomと/dev/urandomの両方の場合で乱数を取得できないとエラーになり，すぐに-1が返ります．（つまり，乱数が取得できるまで待機しません．）
/dev/randomと/dev/urandomの違いは，エントロピープールで利用可能なエントロピーがなくなった場合の振る舞いです．
エントロピーがなくなった場合，/dev/randomはエントロピーが一定の数になるまでブロックしますが，/dev/urandomは過去に利用した乱数を再利用します．
なので，/dev/randomより/dev/urandomの方が利便性は高いですが，安全性は低いというトレードオフがあります．
また，/dev/randomは暗号論的擬似乱数生成器に利用できますが，/dev/urandomや上記で紹介した擬似乱数生成器は暗号論的擬似乱数生成器に利用できないことに注意して下さい．
getrandomシステムコールを利用するコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <stdlib.h> #include <string.h> #include <unistd.h> #include <sys/random.h> #define BUFSIZE 8 int main(int argc, char *argv[]) { ssize_t size; unsigned char buf[BUFSIZE]; int flags; if (argc != 2) { fprintf(stderr, "Usage: %s [0|n|r|nr]\n", argv[0]); exit(1); } if (strcmp(argv[1], "0") == 0) { /* read from /dev/urandom with blocking * if the entropy pool has not yet been initialized. */ flags = 0; } else if (strcmp(argv[1], "n") == 0) { /* read from /dev/urandom with blocking. * return -1 if the entropy pool has not yet been initialized. */ flags = GRND_NONBLOCK; } else if (strcmp(argv[1], "r") == 0) { /* read from /dev/random with blocking if no random bytes are available. */ flags = GRND_RANDOM; } else if (strcmp(argv[1], "nr") == 0) { /* read from /dev/random without blocking. * return -1 if no random bytes are available. */ flags = GRND_NONBLOCK | GRND_RANDOM; } else { fprintf(stderr, "Error: unknown options: %s\n", argv[1]); exit(2); } if ((size = getrandom(buf, BUFSIZE, flags)) == -1) { perror("getrandom"); exit(3); } write(1, buf, size); return 0; } |
実行結果は以下になります．
getrandomシステムコールで取得するデータはバイナリデータなので，Base64でエンコードしてテキストデータに変換しています．
利用可能なエントロピーがない場合はブロッキングされるか，ノンブロッキングの場合は-1を返すことに注意して下さい．
|
1 2 3 4 5 6 7 8 9 |
$ gcc getrandom.c $ a.out 0 | base64 7y3YhivbUHg= $ a.out n | base64 MVc5ChLnLfA= $ a.out r | base64 S15xGjCf $ a.out nr | base64 HaVy8PKE |
Base64を知りたいあなたはこちらからどうぞ．
まとめ
C言語でrand関数と自作関数で乱数を生成する方法を紹介しました．
rand関数と自作関数はget_uniform関数を呼び出すことで，生成した乱数を取得できるコードにしています（getrandomシステムコール以外）．
乱数を生成する方法をまとめると以下になります．
- rand関数（標準ライブラリ）
- rand関数：お手軽に使いやすいが，低い安全性
- drand48_r関数：rand関数より高い安全性だが，GCC拡張（Visual Studioでは利用不可）
- 自作関数
- 線形合同法：お手軽に使いやすいが，低い安全性
- メルセンヌ・ツイスタ：使うのに手間がかかるが，高い安全性
- getrandomシステムコール
- /dev/random：暗号論的擬似乱数生成器に利用できるが，エントロピーがないとブロッキング
- /dev/urandom：暗号論的擬似乱数生成器に利用できないが，エントロピーがないと乱数を再利用
上記のトレードオフを踏まえて，自分の環境に適切な擬似乱数生成器を利用しましょう．
また，乱数の利用例として，モンテカルロ法による円周率を計算しました．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
