/*
* PURPOSE
* clcg4 generator stuff
*
* AUTHORS
* The following code is from L'Ecuyer and Andres "A Randow Number based
* on the combinaison of Four LCG" (distributed at the Pierre L'Ecuyer
* home page with a corresponding paper).
*
* NOTES
* The original code was slightly modified by Bruno Pincon for inclusion
* in Scilab.
*
* list of main modifs :
*
* - lot of routine 's names have changed to have some kind of
* uniformity with the others generators
*
* - add a var is_init so that initialisation is performed inside
* this module (to simplify the interface). And bring modif in
* the different routines :
* if (!is_init) then proceed to initialisation ...
*
* - add a routine advance_state_clcg4 (for compatibility with the
* old package (Scilab used this feature))
*
* - I have change the generator (clcg4 routine) so as it has the
* form (1) in place of (2) (see the joined paper of L'Ecuyer &
* Andres) :
*
* From the 4 LCG :
*
* x_{j,n} = a_j * x_{j,n-1} mod m_j 0 <= j <= 3
*
* The output with form (2) (original form in this code) :
*
* z_n = ( sum_j delta_j * x_{j,n} / m_j ) mod 1
*
* have been changed in the form (1) :
*
* z_n = ( sum_j delta_j * x_{j,n} ) mod m_1 (then u_n = z_n / m_1)
*
* to have some "uniformity" with all the others generators (which
* gives integers). Also it is better for the uin(a,b) generation
* to start from integers.
*/
/*---------------------------------------------------------------------*/
/* clcg4.c Implementation module */
/*---------------------------------------------------------------------*/
#include "../graphics/Math.h" /* to use sciprint */
#include "clcg4.h"
#include <math.h> /* for floor */
/***********************************************************************/
/* Private part. */
/***********************************************************************/
#define H 32768 /* = 2^15 : use in MultModM. */
static long aw[4], avw[4], /* a[j]^{2^w} et a[j]^{2^{v+w}}. */
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];
/* Initial seed, previous seed, and current seed. */
static int is_init = 0;
static long v_default = 31;
static long w_default = 41;
static long MultModM (long s, long t, long M)
/* Returns (s*t) MOD M. Assumes that -M < s < M and -M < t < M. */
/* See L'Ecuyer and Cote (1991). */
{
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;
}
return R;
}
void comp_aw_and_avw(long v, long w)
{
int i, j;
for (j = 0; j < 4; j++)
{
aw [j] = a [j];
for (i = 1; i <= w; i++)
aw [j] = MultModM (aw [j], aw [j], m[j]);
avw [j] = aw [j];
for (i = 1; i <= v; i++)
avw [j] = MultModM (avw [j], avw [j], m[j]);
}
}
void init_clcg4(long v, long w)
{
/* currently the scilab interface don't let the user chooses
* v and w (always v_default and w_default) so this routine
* is in the "private" part (also because initialisation is
* always perform inside this module, depending of the var
* is_init)
*/
double sd[4] = {11111111., 22222222., 33333333., 44444444.};
comp_aw_and_avw(v, w);
set_initial_seed_clcg4(sd[0], sd[1], sd[2], sd[3]);
}
int verif_seeds_clcg4(double s0, double s1, double s2, double s3)
{
/* verify that the seeds are "integers" and are in the good range */
if ( s0 == floor(s0) && s1 == floor(s1) &&
s2 == floor(s2) && s3 == floor(s3) &&
1 <= s0 && s0 <= 2147483646 &&
1 <= s1 && s1 <= 2147483542 &&
1 <= s2 && s2 <= 2147483422 &&
1 <= s3 && s3 <= 2147483322 )
return ( 1 );
else
return ( 0 );
}
void display_info_clcg4()
{
/* display the seeds range (in case of error) */
sciprint("\n\r bad seeds for clcg4, must be integers with s1 in [1, 2147483646]");
sciprint("\n\r s2 in [1, 2147483542]");
sciprint("\n\r s3 in [1, 2147483422]");
sciprint("\n\r s4 in [1, 2147483322]");
}
/*---------------------------------------------------------------------*/
/* Public part. */
/*---------------------------------------------------------------------*/
int set_seed_clcg4(int g, double s0, double s1, double s2, double s3)
{
if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
if ( verif_seeds_clcg4(s0, s1, s2, s3) )
{
Ig [0][g] = (long) s0; Ig [1][g] = (long) s1;
Ig [2][g] = (long) s2; Ig [3][g] = (long) s3;
init_generator_clcg4(g, InitialSeed);
sciprint("\n\r => be aware that you have may lost synchronization");
sciprint("\n\r between the virtual gen %d and the others !", g);
sciprint("\n\r use grand(\"setall\", s1, s2, s3, s4) if you want recover it.");
return ( 1 );
}
else
{
display_info_clcg4();
return ( 0 );
}
}
void get_state_clcg4(int g, double s[4])
{
int j;
if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
for (j = 0; j < 4; j++) s [j] = (double) Cg [j][g];
}
void init_generator_clcg4(int g, SeedType Where)
{
int j;
if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
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];
}
}
void advance_state_clcg4(int g, int k)
{
long int b[4];
int i, j;
if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
for ( j = 0 ; j < 4 ; j++ )
{
b[j] = a[j];
for ( i = 1 ; i <= k ; i++ )
b[j] = MultModM( b[j], b[j], m[j]);
Ig[j][g] = MultModM ( b[j], Cg[j][g], m[j] );
}
init_generator_clcg4(g, InitialSeed);
}
int set_initial_seed_clcg4(double s0, double s1, double s2, double s3)
{
int g, j;
if (! is_init ) comp_aw_and_avw(v_default,w_default);
if ( ! verif_seeds_clcg4(s0, s1, s2, s3) )
{
display_info_clcg4();
return ( 0 );
};
is_init = 1;
Ig [0][0] = (long) s0;
Ig [1][0] = (long) s1;
Ig [2][0] = (long) s2;
Ig [3][0] = (long) s3;
init_generator_clcg4(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]);
init_generator_clcg4(g, InitialSeed);
}
return ( 1 );
}
unsigned long clcg4(int g)
{
/* Modif Bruno : the generator have now the form (1) in place of (2) */
long k,s;
double u;
if (! is_init ) {init_clcg4(v_default,w_default); is_init = 1; };
/* advance the 4 LCG */
s = Cg [0][g]; k = s / 46693;
s = 45991 * (s - k * 46693) - k * 25884;
if (s < 0) s = s + 2147483647; Cg [0][g] = 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;
s = Cg [2][g]; k = s / 15499;
s = 138556 * (s - k * 15499) - k * 3979;
if (s < 0) s = s + 2147483423; Cg [2][g] = s;
s = Cg [3][g]; k = s / 43218;
s = 49689 * (s - k * 43218) - k * 24121;
if (s < 0) s = s + 2147483323; Cg [3][g] = s;
/* final step */
u = (double)(Cg[0][g] - Cg[1][g]) + (double)(Cg[2][g] - Cg[3][g]);
/* we must do u mod 2147483647 with u in [- 4294966863 ; 4294967066 ] : */
if (u < 0) u += 2147483647;
if (u < 0) u += 2147483647;
if (u >= 2147483647) u -= 2147483647;
if (u >= 2147483647) u -= 2147483647;
return ((unsigned long) u );
}
syntax highlighted by Code2HTML, v. 0.9.1