動機

C言語で１から乱数を返す関数をつくってみました．

C言語には標準で乱数を生成するrand()関数が用意されています．ただ，この関数は１つの乱数しか返さないため，モンテカルロ法のような乱数を大量に使う場合少々不便です．
そのため，空の配列と要素数を与えるとその配列に乱数を入れてくれるような関数があれば便利だと思い，実際に作ってみました．

せっかくなので，乱数の生成に関するアルゴリズムにも少々触れながら解説していきたいとおもいます．

コード

rand()の拡張版: randint(uint32_t *xp, int nums)

まず第一ステップは，整数の擬似乱数の配列を返す関数の作成です．
rand()関数をfor文で回すのもバカバカしいので，乱数の生成部分も自分で書きました．

ところでみなさんはコンピュータがどのようにして乱数を生成しているか知っていますか？
実はある初期値(seed)に対して定められた変形を行っているだけなのです．
正確には乱数ではないので，「擬似乱数」とよばれます．

今回は数あるアルゴリズムの中から，単純でありながら周期が非常に長い"xorshft32"を採用しました．
詳しい説明はこちらに譲りますが，なんとこのアルゴリズム，名前の通りXORとビットシフトのみで構成されています．
これだけで乱数が生成されるのは不思議ですね．

コードはこちら

randint

int randint(uint32_t *xp, int nums)
{
  int i;
  uint32_t x;
  clock_t clk;
  time_t sec;
  static uint32_t seed = 0;

  if(seed == 0) {
    clk = clock();
    time(&sec);
    seed = clk * sec;
  }

  for(i=0;i<nums;i++){
    x = (i > 0) ? xp[i-1] : seed;
    x = x ^ (x << 13);
    x = x ^ (x >> 17);
    x = x ^ (x << 15);
    xp[i] = x;
  }
  seed = x;
  return i;
}

一様乱数の生成: random(double *xp, int nums)

こちらは0から1までの一様乱数（どの数も同じ確率で出てくる乱数）を生成する関数です．
randintで生成した乱数に対して0x7fffffffでマスクすることで下位31ビットを取り出し，これを0x80000000で割ることで0から1までの数値に変換しています．

コードはこちら

random

int random(double *xp, int nums) 
{
  int i;
  uint32_t mask, maxint;
  uint32_t yp[nums];

  mask = 0x7fffffff;
  maxint = 0x80000000;
  randint(yp, nums);
  #pragma omp parallel
  {
    #pragma omp for
    for (i = 0; i<nums; i++) {
      uint32_t my;
      my = yp[i] & mask;
      xp[i] = (double)my / maxint;
    }
  }
  return i;
}

正規乱数の生成: randn(double *xp, int nums)

最後に，平均0，標準偏差1の正規分布に従う正規乱数を生成する関数です．
正規乱数の生成にはボックス＝ミュラー法を用いました．これは一様乱数から正規乱数を得るアルゴリズムです．2つの一様乱数を用意し，片方の乱数にはlogとsqrtを，もう片方の乱数には三角関数をかませ，両者をかけあわせるとあら不思議，正規乱数が出来上がってしまうのです．

コードはこちら

randn

int randn(double *xp, int nums){
  int i, nums2;
  nums2 = (nums % 2 == 0) ? nums : nums + 1;
  double yp[nums2];
  random(yp, nums2);
  #pragma omp parallel
  {
    #pragma omp for
    for (i = 0; i<nums2; i += 2) {
      double sqrt_log_x, theta, cos_y, sin_y;
      sqrt_log_x = sqrt(-2.0 * log(yp[i]));
      theta = 2.0 * M_PI * yp[i+1];
      cos_y = cos(theta);
      sin_y = sin(theta);
      xp[i] = sqrt_log_x * cos_y;
      if (i < nums) {
	xp[i+1] = sqrt_log_x * sin_y;
      }
    }
  }
  return  nums;
}

最後に

いかがでしたでしょうか．乱数って調べてみると奥が深いんですね．

さて，余談ですが，ところどころ#pragma ompのような呪文が書かれているのに気づかれたでしょうか．
知っている方も多いと思いますが，実はこの呪文，付け加えるだけで勝手に並列処理をしてくれるんです．
ヘッダーに#include <omp.h>を付け加え，コンパイルの時にオプションで-fopenmpを指定すると有効化されます．ぜひお試しください．私の実行環境(Intel Core i7)では約２倍の速度がでました．

あと，余計なお世話かもしれませんが，math.hを使用するときはコンパイル時にオプション-lmを指定しないとエラーを吐かれるので注意してくださいね．ではでは．

myrand.c

# include <stdio.h>
# include <stdint.h>
# include <time.h>
# include <math.h>
# include <omp.h>


int randint(uint32_t *xp, int nums)
{
  int i;
  uint32_t x;
  clock_t clk;
  time_t sec;
  static uint32_t seed = 0;

  if(seed == 0) {
    clk = clock();
    time(&sec);
    seed = clk * sec;
  }

  for(i=0;i<nums;i++){
    x = (i > 0) ? xp[i-1] : seed;
    x = x ^ (x << 13);
    x = x ^ (x >> 17);
    x = x ^ (x << 15);
    xp[i] = x;
  }
  seed = x;
  return i;
}

int random(double *xp, int nums) 
{
  int i;
  uint32_t mask, maxint;
  uint32_t yp[nums];

  mask = 0x7fffffff;
  maxint = 0x80000000;
  randint(yp, nums);
  #pragma omp parallel
  {
    #pragma omp for
    for (i = 0; i<nums; i++) {
      uint32_t my;
      my = yp[i] & mask;
      xp[i] = (double)my / maxint;
    }
  }
  return i;
}
  
int randn(double *xp, int nums){
  int i, nums2;
  nums2 = (nums % 2 == 0) ? nums : nums + 1;
  double yp[nums2];
  random(yp, nums2);
  #pragma omp parallel
  {
    #pragma omp for
    for (i = 0; i<nums2; i += 2) {
      double sqrt_log_x, theta, cos_y, sin_y;
      sqrt_log_x = sqrt(-2.0 * log(yp[i]));
      theta = 2.0 * M_PI * yp[i+1];
      cos_y = cos(theta);
      sin_y = sin(theta);
      xp[i] = sqrt_log_x * cos_y;
      if (i < nums) {
	xp[i+1] = sqrt_log_x * sin_y;
      }
    }
  }
  return  nums;
}


int main(){
  int i;
  int nums = 1000;
  double x[nums];
  randn(x, nums);
  for (i = 0; i < nums; i++) {
    printf("%lf\n", x[i]);
  }
  return 0;
}

追記

randint関数がほぼ同時に呼ばれた時にseedが更新されないバグを直しました．