samer@0: /***********************************************************************\ samer@0: * samer@0: * File: RngStream.c for multiple streams of Random Numbers samer@0: * Language: ANSI C samer@0: * Copyright: Pierre L'Ecuyer, University of Montreal samer@0: * Date: 14 August 2001 samer@0: * License: GPL version 2 or later samer@0: * samer@0: * Notice: Please contact P. L'Ecuyer at samer@0: * for commercial purposes. samer@0: * samer@0: \***********************************************************************/ samer@0: samer@0: samer@0: #include "RngStream.h" samer@0: #include samer@0: #include samer@0: #include samer@0: #include samer@0: samer@0: /*---------------------------------------------------------------------*/ samer@0: /* Private part. */ samer@0: /*---------------------------------------------------------------------*/ samer@0: samer@0: typedef struct rng_state *RngStream; samer@0: samer@0: #define norm 2.328306549295727688e-10 samer@0: #define m1 4294967087.0 samer@0: #define m2 4294944443.0 samer@0: #define a12 1403580.0 samer@0: #define a13n 810728.0 samer@0: #define a21 527612.0 samer@0: #define a23n 1370589.0 samer@0: samer@0: #define two17 131072.0 samer@0: #define two53 9007199254740992.0 samer@0: #define fact 5.9604644775390625e-8 /* 1 / 2^24 */ samer@0: samer@0: samer@0: samer@0: /* Default initial seed */ samer@0: static double nextSeed[6] = { 12345, 12345, 12345, 12345, 12345, 12345 }; samer@0: samer@0: samer@0: /* The following are the transition matrices of the two MRG components */ samer@0: /* (in matrix form), raised to the powers 2^0, 2^76, and 2^127, resp.*/ samer@0: samer@0: static double A1p0[3][3] = { samer@0: { 0.0, 1.0, 0.0 }, samer@0: { 0.0, 0.0, 1.0 }, samer@0: { -810728.0, 1403580.0, 0.0 } samer@0: }; samer@0: samer@0: static double A2p0[3][3] = { samer@0: { 0.0, 1.0, 0.0 }, samer@0: { 0.0, 0.0, 1.0 }, samer@0: { -1370589.0, 0.0, 527612.0 } samer@0: }; samer@0: samer@0: samer@0: static double A1p76[3][3] = { samer@0: { 82758667.0, 1871391091.0, 4127413238.0 }, samer@0: { 3672831523.0, 69195019.0, 1871391091.0 }, samer@0: { 3672091415.0, 3528743235.0, 69195019.0 } samer@0: }; samer@0: samer@0: static double A2p76[3][3] = { samer@0: { 1511326704.0, 3759209742.0, 1610795712.0 }, samer@0: { 4292754251.0, 1511326704.0, 3889917532.0 }, samer@0: { 3859662829.0, 4292754251.0, 3708466080.0 } samer@0: }; samer@0: samer@0: static double A1p127[3][3] = { samer@0: { 2427906178.0, 3580155704.0, 949770784.0 }, samer@0: { 226153695.0, 1230515664.0, 3580155704.0 }, samer@0: { 1988835001.0, 986791581.0, 1230515664.0 } samer@0: }; samer@0: samer@0: static double A2p127[3][3] = { samer@0: { 1464411153.0, 277697599.0, 1610723613.0 }, samer@0: { 32183930.0, 1464411153.0, 1022607788.0 }, samer@0: { 2824425944.0, 32183930.0, 2093834863.0 } samer@0: }; samer@0: samer@0: samer@0: samer@0: samer@0: samer@0: /*-------------------------------------------------------------------------*/ samer@0: samer@0: static double MultModM (double a, double s, double c, double m) samer@0: /* Compute (a*s + c) % m. m must be < 2^35. Works also for s, c < 0 */ samer@0: { samer@0: double v; samer@0: long a1; samer@0: v = a * s + c; samer@0: if ((v >= two53) || (v <= -two53)) { samer@0: a1 = (long) (a / two17); samer@0: a -= a1 * two17; samer@0: v = a1 * s; samer@0: a1 = (long) (v / m); samer@0: v -= a1 * m; samer@0: v = v * two17 + a * s + c; samer@0: } samer@0: a1 = (long) (v / m); samer@0: if ((v -= a1 * m) < 0.0) samer@0: return v += m; samer@0: else samer@0: return v; samer@0: } samer@0: samer@0: samer@0: /*-------------------------------------------------------------------------*/ samer@0: samer@0: samer@0: static void print_vec(double *p) samer@0: { samer@0: int i; samer@0: for (i=0; i<6; i++) printf("%.0lf ",p[i]); samer@0: printf("\n"); samer@0: } samer@0: samer@0: samer@0: static void MatVecModM (double A[3][3], double s[3], double v[3], double m) samer@0: /* Returns v = A*s % m. Assumes that -m < s[i] < m. */ samer@0: /* Works even if v = s. */ samer@0: { samer@0: int i; samer@0: double x[3]; samer@0: for (i = 0; i < 3; ++i) { samer@0: x[i] = MultModM (A[i][0], s[0], 0.0, m); samer@0: x[i] = MultModM (A[i][1], s[1], x[i], m); samer@0: x[i] = MultModM (A[i][2], s[2], x[i], m); samer@0: } samer@0: memcpy(v,x,3*sizeof(double)); samer@0: } samer@0: samer@0: samer@0: static void MatMatModM (double A[3][3], double B[3][3], double C[3][3], samer@0: double m) samer@0: /* Returns C = A*B % m. Work even if A = C or B = C or A = B = C. */ samer@0: { samer@0: int i, j; samer@0: double V[3], W[3][3]; samer@0: for (i = 0; i < 3; ++i) { samer@0: for (j = 0; j < 3; ++j) V[j] = B[j][i]; samer@0: MatVecModM (A, V, V, m); samer@0: for (j = 0; j < 3; ++j) W[j][i] = V[j]; samer@0: } samer@0: memcpy(C,W,9*sizeof(double)); samer@0: } samer@0: samer@0: static void MatTwoPowModM (double A[3][3], double B[3][3], double m, long e) samer@0: /* Compute matrix B = (A^(2^e) % m); works even if A = B */ samer@0: { samer@0: int i, j; samer@0: samer@0: /* initialize: B = A */ samer@0: if (A != B) { samer@0: for (i = 0; i < 3; i++) { samer@0: for (j = 0; j < 3; ++j) samer@0: B[i][j] = A[i][j]; samer@0: } samer@0: } samer@0: /* Compute B = A^{2^e} */ samer@0: for (i = 0; i < e; i++) samer@0: MatMatModM (B, B, B, m); samer@0: } samer@0: samer@0: samer@0: /* returns the next integer in the sequence samer@0: * new state goes in g1 samer@0: * g1 can be the same as g, resulting in samer@0: * an update in place. samer@0: */ samer@0: double next(RngStream g, RngStream g1) samer@0: { samer@0: double p1, p2; samer@0: samer@0: // printf("before: "); print_vec(g->Cg); samer@0: /* Component 1 */ samer@0: p1 = a12 * g->Cg[1] - a13n * g->Cg[0]; samer@0: p1 -= m1*floor(p1/m1); samer@0: g1->Cg[0] = g->Cg[1]; samer@0: g1->Cg[1] = g->Cg[2]; samer@0: g1->Cg[2] = p1; samer@0: samer@0: /* Component 2 */ samer@0: p2 = a21 * g->Cg[5] - a23n * g->Cg[3]; samer@0: p2 -= m2*floor(p2/m2); samer@0: g1->Cg[3] = g->Cg[4]; samer@0: g1->Cg[4] = g->Cg[5]; samer@0: g1->Cg[5] = p2; samer@0: samer@0: // printf("after: "); print_vec(g1->Cg); samer@0: /* Combination */ samer@0: return (p1 > p2) ? (p1 - p2) : (p1 - p2 + m1); samer@0: } samer@0: samer@0: samer@0: /*-------------------------------------------------------------------------*/ samer@0: static int CheckSeed (unsigned seed[6]) samer@0: { samer@0: /* Check that the seeds are legitimate values. Returns 0 if legal seeds, samer@0: * i:1..6 if seed[i] is bad samer@0: * 8 if seeds 1..3 are zero samer@0: * 9 if seeds 4..6 are zero samer@0: */ samer@0: int i; samer@0: samer@0: for (i = 0; i < 3; ++i) { samer@0: if (seed[i] >= m1) return i+1; samer@0: } samer@0: for (i = 3; i < 6; ++i) { samer@0: if (seed[i] >= m2) return i+1; samer@0: } samer@0: if (seed[0] == 0 && seed[1] == 0 && seed[2] == 0) return 8; samer@0: if (seed[3] == 0 && seed[4] == 0 && seed[5] == 0) return 9; samer@0: samer@0: return 0; samer@0: } samer@0: samer@0: samer@0: /*---------------------------------------------------------------------*/ samer@0: /* Public part. */ samer@0: /*---------------------------------------------------------------------*/ samer@0: samer@0: samer@0: void RngStream_InitState(RngStream g) { samer@0: memcpy(g->Cg,nextSeed,6*sizeof(double)); samer@0: //printf("init: "); print_vec(g->Cg); samer@0: } samer@0: samer@0: void RngStream_InitJump(struct rng_jump *j, int e) samer@0: { samer@0: void *a1, *a2; samer@0: samer@0: if (e>=127) { a1=A1p127; a2=A2p127; e-=127; } samer@0: else if (e==76) { a1=A1p76; a2=A2p76; e-=76; } samer@0: else { a1=A1p0; a2=A2p0; } samer@0: samer@0: memcpy(j->T1,a1,9*sizeof(double)); samer@0: memcpy(j->T2,a2,9*sizeof(double)); samer@0: MatTwoPowModM(j->T1,j->T1,m1,e); samer@0: MatTwoPowModM(j->T2,j->T2,m2,e); samer@0: } samer@0: samer@0: void RngStream_DoubleJump(struct rng_jump *j1, struct rng_jump *j2) samer@0: { samer@0: MatMatModM(j1->T1,j1->T1,j2->T1,m1); samer@0: MatMatModM(j1->T2,j1->T2,j2->T2,m2); samer@0: } samer@0: samer@0: void RngStream_Advance(struct rng_jump *j, RngStream g1, RngStream g2) samer@0: { samer@0: MatVecModM (j->T1, &g1->Cg[0], &g2->Cg[0], m1); samer@0: MatVecModM (j->T2, &g1->Cg[3], &g2->Cg[3], m2); samer@0: } samer@0: samer@0: double RngStream_Float(RngStream g, RngStream g1) { samer@0: return norm*next(g,g1); samer@0: } samer@0: samer@0: double RngStream_Double(RngStream g, RngStream g1) { samer@0: double u = RngStream_Float(g,g1); samer@0: u += fact*RngStream_Float(g1,g1); samer@0: return (u < 1.0) ? u : (u - 1.0); samer@0: } samer@0: samer@0: double RngStream_Raw(RngStream g, RngStream g1) { samer@0: return next(g,g1); samer@0: } samer@0: samer@0: /*-------------------------------------------------------------------------*/ samer@0: samer@0: int RngStream_SetState (RngStream g, unsigned seed[6]) samer@0: { samer@0: int i; samer@0: if (CheckSeed (seed)) samer@0: return -1; /* FAILURE */ samer@0: for (i = 0; i < 6; ++i) samer@0: g->Cg[i] = (double)seed[i]; samer@0: samer@0: return 0; /* SUCCESS */ samer@0: } samer@0: samer@0: samer@0: void RngStream_GetState (RngStream g, unsigned seed[6]) samer@0: { samer@0: int i; samer@0: for (i = 0; i < 6; ++i) seed[i] = (unsigned)g->Cg[i]; samer@0: } samer@0: