C LANGUAGE TECHNOLOGY

【C言語】rand関数と自作関数で乱数の生成【モンテカルロ法で円周率の計算】

2021年6月6日

悩んでいる人

C言語でrand関数と自作関数で乱数の生成方法を教えて!

こういった悩みにお答えします.

本記事の信頼性

  • リアルタイムシステムの研究歴12年.
  • 東大教員の時に,英語でOS(Linuxカーネル)の授業.
  • 2012年9月~2013年8月にアメリカのノースカロライナ大学チャペルヒル校(UNC)コンピュータサイエンス学部で客員研究員として勤務.C言語でリアルタイムLinuxの研究開発.
  • プログラミング歴15年以上,習得している言語: C/C++PythonSolidity/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関数による乱数の生成

rand関数は,[0, RAND_MAX]の範囲の擬似乱数整数を返します.

RAND_MAXは私の環境のstdlib.hで以下のように定義されています.

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は含まれないことに注意して下さい.

実行結果は以下になります.

rand/srand関数による乱数の生成コード

rand関数では,上記のように同じ疑似乱数整数しか生成できません.

そこで,srand関数で疑似乱数整数の新しい種を引数に設定することで,異なる擬似乱数整数を生成できます.

rand/srand関数による乱数の生成コードは以下になります.

実行結果は以下になります.

rand関数のみを利用した場合と異なる擬似乱数整数を生成していることがわかります.

rand/srand関数による時刻で変動する乱数の生成コード

rand/srand関数を利用したコードにより生成する擬似乱数整数は変わりましたが,種の値が100と一定なので同じ乱数しか生成できません.

そこで,種の値を時刻に設定することで,実行毎に生成する乱数を変更します.

rand/srand関数による時刻で変動する乱数の生成コードは以下になります.

15行目でtime関数を呼び出して,その返り値をsrand関数の引数に設定します.

実行結果の例は以下になります.

実行毎に生成する乱数が変動していることがわかります.

※time関数の返り値は秒なので,実行毎に1秒以上の間隔を空けましょう.

drand48_r関数による乱数の生成コード

rand_r関数のseedp引き数が指す値により提供される状態は非常に小さな空間なので,弱い疑似乱数生成器(予測性が高い擬似乱数生成器)になってしまいます.

drand48_r関数は,rand_r関数より大きな空間で擬似乱数を生成します.

rand_r関数の空間サイズは32ビット(または16ビット)なのに対して,drand48_r関数の空間サイズは48ビットになります.

ただし,drand48_r関数はGCC拡張なので,Visual Studioでは使えないことに注意して下さい.

drand48_r関数による乱数の生成コードは以下になります.

実行結果は以下になります.

自作関数による乱数の生成

自作関数による乱数の生成方法として以下の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)\)と規格化することで一様乱数を生成できます.

実行結果は以下になります.

メルセンヌ・ツイスタによる乱数生成

メルセンヌ・ツイスタとは,線形合同法の多次元で規則的に分布(不均等分布)するという欠点を解決した疑似乱数生成器です.

不均等分布すると連続する値間の相関性が高くなり,統計的にランダムになるため予測性が低くなります.

メルセンヌ・ツイスタによる乱数の生成コードは以下になります.

上記のコードをビルドするためには,メルセンヌ・ツイスタのライブラリから必要なコードのみを抽出したものをこちらからダウンロードして下さい.

mt.zipファイルはunzipコマンドで解凍します.

mt.zipファイルの中身は以下になります.

  • SFMT-common.h, SFMT-params.h, SFMT-params19937.h, SFMT.c, SFMT.h:メルセンヌ・ツイスタのライブラリ
  • mt.c:メルセンヌ・ツイスタのテストプログラム
  • Makefile:上記のコードのビルド用

※メルセンヌ・ツイスタのライブラリはとても長いコードなので説明は省略します.

上記のコードをビルドするMakefileは以下になります.

mt.zipファイルを解凍した時に作成したmtディレクトリに移動し,makeでビルドして実行すると,以下のようになります.

モンテカルロ法で円周率の計算

モンテカルロ法とは,乱数を利用した試行を繰り返すことにより近似解を求める手法です.

メルセンヌ・ツイスタを利用してモンテカルロ法で円周率を求めるコード一式をこちらからダウンロードして同様に展開して下さい.

mt.zipとの主な違いはmt.cがmt_monte_carlo.cに変更したことで,そのコードは以下になります.

37行目で「(double) count / NR_LOOPS * 4」と計算する理由は,モンテカルロ法で円周率を求めるの記事がわかりやすいので参考にして下さい.

実行結果は以下になります.

3回実行したところ,円周率の値は3.141592の近似解になっていることがわかります.

getrandomシステムコールで乱数の取得

getrandomシステムコールは,bufが指すバッファに最大buflenのランダムなバイトを詰めます.

これらのバイトは,ユーザ空間の乱数生成器のシードや暗号演算の目的で利用できます。

getrandomシステムコールは,エントロピープール(Linuxが保持するバッファで複数の環境ノイズ)から乱数を取得します.

エントロピープールが保持する複数の環境ノイズには,入力デバイスやストレージで発生する割り込みを元にした情報が格納されます.

エントロピープールの利用可能なエントロピー数を知りたい場合は,以下のコマンドを実行して下さい.

エントロピープールの上限値を知りたい場合は,以下のコマンドを実行して下さい.

デフォルトでは,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システムコールを利用するコードは以下になります.

実行結果は以下になります.

getrandomシステムコールで取得するデータはバイナリデータなので,Base64でエンコードしてテキストデータに変換しています.

利用可能なエントロピーがない場合はブロッキングされるか,ノンブロッキングの場合は-1を返すことに注意して下さい.

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社で自分に合うスクールを見つけましょう.後悔はさせません!

-C LANGUAGE, TECHNOLOGY
-, , , , , , , , ,