1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
|
/*
* PURPOSE
* uniform random number generator developed by Pierre
* Lecuyer based on a clever and tested combination of
* two linear congruential sequences
*
* s1 <- a1*s1 mod m1 , a1 = 40014, m1 = 2147483563
* s2 <- a2*s2 mod m2 , a2 = 40692, m2 = 2147483399
*
* output <- s1-s2 mod (m1 - 1)
*
* so output is in [0, 2147483561], period about 2.3 10^18
*
* The state is given by (s1, s2). In case of a user
* modification of the state we must have :
*
* s1 in [1, m1-1]
* s2 in [1, m2-1]
*
* ORIGIN
* The basic code is provided at the Luc Devroye 's home page.
* Modifications by Bruno Pincon (in particular added routines
* to set and get the state, and modify the generator to get
* exactly s1-s2 mod (m1 - 1) for "coherence" with the others
* generators : provides numbers in [0, MaxRngInt(generator)]
* (see NOTE some lines after)
*
*/
#include "../graphics/Math.h" /* to use sciprint */
#include <math.h>
/* initial default state (seeds) : */
static long s1 = 1234567890 ;
static long s2 = 123456789 ;
unsigned long clcg2()
{
register long k,z;
/* s1 = a1*s1 mod m1 (Schrage 's method) */
k= s1 /53668;
s1 =40014*(s1%53668)-k*12211;
if (s1 < 0) s1 += 2147483563;
/* s2 = a2*s2 mod m2 (Schrage 's method) */
k=s2/52774;
s2=40692*(s2%52774)-k*3791;
if (s2 < 0) s2 += 2147483399;
/* final step : z = s1-s2 mod m1-1 */
z = s1 - s2; /* here z is in [2-m2,m1-2] */
if (z < 0) z += 2147483562;
/* NOTE : in the original implementation the final test is :
* if (z < 1) z += 2147483562;
*
* which is not exactly z = s1-s2 mod (m1 - 1)
*
* This is also why it is different from the version used by
* randlib.
*/
return( (unsigned long) z );
}
int set_state_clcg2(double g1, double g2)
{
if ( g1 == floor(g1) && g2 == floor(g2) &&
1 <= g1 && g1 <= 2147483562 &&
1 <= g2 && g2 <= 2147483398 )
{
s1 = (long) g1;
s2 = (long) g2;
return ( 1 );
}
else
{
sciprint("\n\r bad seeds for clcg2, must be integers with s1 in [1, 2147483562]");
sciprint("\n\r and s2 in [1, 2147483398]\n\r");
return ( 0 );
}
}
void get_state_clcg2(double g[])
{
g[0] = (double) s1;
g[1] = (double) s2;
}
|
