#include "TrueBasic.h"

#ifdef _MSC_VER
float __stdcall RND(int *seed)
{
	return rnd(*seed);
}
#endif

float rnd(int seed)
{
// 混合法を併用した線形合同法の乱数発生プログラム
// "Numerical Recipes" の第二版の ran1 を基礎にしている
#define a 16807
#define m 2147483647  // 2^31 - 1
#define rm 1.0/m
#define q 127773      // m = a*q + p
#define p 2836
#define n 32
#define ndiv (1 + (m-1)/n)
#define rmax (1.0 - 1.2e-7)
   static int r[n+1],r0,r1;
   int j,k;
   if (seed < 0) seed = 584287;
   if (seed != 0)
// 乱数表の初期化
      {
       r1 = abs(seed);
    for (j = n+8;j>=1;j--)
      {
      k = r1/q;
      r1 = a*(r1-k*q) - p*k;
      if (r1 < 0) r1 = r1 + m;
      if (j < n) r[j] = r1;
        } 
      r0 = r[1];
      } 
// 初期化を行わないときの開始点
// オーバーフローすることなく r1 = mod(a*r1,m) を計算する
    k = r1/q;
    r1 = a*(r1 - k*q) - p*k;
    if (r1 < 0) r1 = r1 + m;
    j = 1 + r0/ndiv;
    r0 = r[j];
    r[j] = r1;
    return (float)min(rm*r0,rmax);
}