/*----------------------------------------------------------------------------*/ /* ael.h */ /* */ /* Module fournissant un generateur (combine de 4 LCG) des valeurs pseudo - */ /* - aleatoires uniformement distribuees sur (0,1) */ /* */ /*----------------------------------------------------------------------------*/ #include #include #include #include "aecl.h" /**********************************************************************************/ /* Partie privee du module. */ /**********************************************************************************/ #define H 32768 /* = 2^15 : utile pour la fonction MultModM */ static long V,W; static long aw[4],avw[4], /* a[j]^W et a[j]^{VW} */ a[4] = { 45991, 207707, 138556, 49689 }, m[4] = { 2147483647, 2147483543, 2147483423, 2147483323 }; static long Ig[4][Maxgen+1], Lg[4][Maxgen+1], Cg[4][Maxgen+1]; /* Le germe initial,le dernier */ /* germe et le germe courant. */ static double c[4] = { 4.6566128774142008e-10, -4.6566131029282642e-10, 4.6566133631368265e-10, -4.6566135799773169e-10 }, d = 2.2117727679125529e-17; /* Les coefficients pour Genval*/ static short j; static long MultModM ( long s, long t, long M ) ; static long PowerModM ( long s, long e, long M ); /* ****************************************************************************** */ /* long MultModM( long s, long t, long M ) */ /* Retourne (s*t) MOD M. Suppose que -M < s < M, -M < t < M. */ /* Voir L'Ecuyer et Cote (1991). */ /* Adapte de MULTLI.MOD (latmrg) */ /* ****************************************************************************** */ static long MultModM( long s, long t, long M ) { long R, S0, S1, q, qh, rh, k ; if( s < 0 ) s += M ; if( t < 0 ) t += M ; if( s < H ) { S0 = s ; R = 0 ; } else { S1 = s/H ; S0 = s - H*S1 ; qh = M/H ; rh = M - H*qh ; if( S1 >= H ) { S1 -= H ; k = t/qh; R = H*(t - k*qh) - k*rh ; while( R < 0 ) R += M ; } else R = 0 ; if( S1 != 0 ) { q = M/S1 ; k = t/q ; R -= k*(M - S1*q) ; if( R > 0 ) R -= M ; R += S1*(t - k*q) ; while( R < 0 ) R += M ; } k = R/qh ; R = H*(R - k*qh) - k*rh ; while( R < 0 ) R += M ; } if( S0 != 0 ) { q = M/S0 ; k = t/q ; R -= k*(M - S0*q) ; if( R > 0 ) R -= M ; R += S0*(t - k*q) ; while( R < 0 ) R += M ; } /* printf ("Ici \n"); */ return R ; } /* ****************************************************************************** */ /* long PowerModM( long s, long e, long M ) */ /* Retourne (s^e) MOD M. Suppose que -M < s < M. */ /* ****************************************************************************** */ static long PowerModM ( long s, long e, long M ) { long c; if ( e == 0 ) return 1; if ( e == 1 ) return s; c = 1; do { if ( e % 2 ) c = MultModM ( s, c, M ); s = MultModM ( s, s, M ); e = e / 2; } while ( e > 1 ); return MultModM ( s, c, M ); } /*--------------------------------------------------------------------------------*/ /* Partie visible du module. */ /*--------------------------------------------------------------------------------*/ /**********************************************************************************/ /* GetState permet d'examiner l'etat courant du generateur g. Cet etat est retour-*/ /* ne dans dans le tableau s ( 4 valeurs - etats de chacun de LCG utilise. */ /**********************************************************************************/ void GetState ( Gen g, long s[4] ) { for (j = 0; j < 4; j++) Cg [j][g] = s [j]; InitGenerator (g, InitialSeed); } void WriteState ( Gen g ) { for (j = 0; j < 4; j++) { printf(" g = %u",g); printf("\n Ig [%u",j,"] = %lu",Ig[j][g]); printf(" Lg [%u",j,"] = %lu",Lg[j][g]); printf(" Cg [%u",j,"] = %lu",Cg[j][g]); } } /**********************************************************************************/ /* SetSeed initialise le germe initial Ig du generateur g a s. Les valeur du */ /* tableau s doivent satisfaire : 1 <=s[0]<= 2147483646 ; 1 <=s[1]<= 2147483542 */ /* 1 <=s[2]<= 2147483422 ; 1 <=s[3]<= 2147483322 */ /* L'etat courant du generateur g est recalcule. */ /**********************************************************************************/ void SetSeed ( Gen g, long s[4] ) { for (j = 0; j < 4; j++) s [j] = Cg [j][g]; } /**********************************************************************************/ /* Reinitialise le generateur g en fonction de la valeur de Where. */ /**********************************************************************************/ void InitGenerator ( Gen g, SeedType Where ) { for (j = 0; j < 4; j++) { switch (Where) { case InitialSeed : Lg [j][g] = Ig [j][g]; break; case NewSeed : Lg [j][g] = MultModM (aw [j], Lg [j][g], m [j]); break; case LastSeed : break; } Cg [j][g] = Lg [j][g]; } } /**********************************************************************************/ /* SetInitialSeed initialise le germe initial Ig du generateur 0 a s. Les valeur */ /*du tableau s doivent satisfaire : 1 <=s[0]<= 2147483646 ; 1 <=s[1]<= 2147483542 */ /* 1 <=s[2]<= 2147483422 ; 1 <=s[3]<= 2147483322 */ /* Les germes initials d'autres generateurs sont recalcules. */ /**********************************************************************************/ void SetInitialSeed ( long s[4] ) { Gen g; for (j = 0; j < 4; j++) Ig [j][0] = s [j]; InitGenerator (0, InitialSeed); for (g = 1; g <= Maxgen; g++) { for (j = 0; j < 4; j++) Ig [j][g] = MultModM (avw [j], Ig [j][ g-1], m [j]); InitGenerator (g, InitialSeed); } } /**********************************************************************************/ /* Init initialise le package : V =v0 ; W = w0 . */ /**********************************************************************************/ void Init ( long v0, long w0 ) { long sd[4] = { 11111111,22222222, 33333333, 44444444 }; V = v0; W = w0; for (j = 0; j < 4; j++) { aw [j] = PowerModM ( a [j], W, m[j]); avw [j] = PowerModM (aw [j], V, m[j]); } SetInitialSeed ( sd ); } /**********************************************************************************/ /* Fourni une valeur aleatoire suivant la loi uniforme (0,1) en utilisant le */ /* generateur g */ /**********************************************************************************/ double Genval ( Gen g ) { long k,s; double u; u = d; s = Cg [0][g]; k = s / 46693; s = 45991 * (s - k * 46693) - k * 25884; if (s < 0) s = s + 2147483647; Cg [0][g] = s; u = u + c [0] * s; s = Cg [1][g]; k = s / 10339; s = 207707 * (s - k * 10339) - k * 870; if (s < 0) s = s + 2147483543; Cg [1][g] = s; u = u + c [1] * s; if (u < 0) u = u + 1.0; s = Cg [2][g]; k = s / 15499; s = 138556 * (s - k * 15499) - k * 3979; if (s < 0.0) s = s + 2147483423; Cg [2][g] = s; u = u + c [2] * s; if (u > 1.0) u = u - 1.0; s = Cg [3][g]; k = s / 43218; s = 49689 * (s - k * 43218) - k * 24121; if (s < 0) s = s + 2147483323; Cg [3][g] = s; u = u + c [3] * s; if (u < 0) u = u + 1.0; return (u); } /**********************************************************************************/ /* InitDefault initialise le package : V =2^20 ; W = 2^30 */ /**********************************************************************************/ void InitDefault () { Init (1024 * 1024, 1024 * 1024 * 1024); }