randlib: Fixed bug #8560
[scilab.git] / scilab / modules / randlib / src / c / fsultra.c
1 /* 
2  * 
3  * Copyright (C) 2010 - DIGITEO - Michael Baudin
4  * Copyright (C) 2004 - Bruno Pincon
5  * Copyright (C) 1992 - Arif Zaman
6  * Copyright (C) 1992 - George Marsaglia
7  * 
8 FSU - ULTRA     The greatest random number generator that ever was
9                 or ever will be.  Way beyond Super-Duper.
10                 (Just kidding, but we think its a good one.)
11
12 Authors:        Arif Zaman (arif@stat.fsu.edu) and
13                 George Marsaglia (geo@stat.fsu.edu).
14
15 Date:           27 May 1992
16
17 Version:        1.05
18
19 Copyright:      To obtain permission to incorporate this program into
20                 any commercial product, please contact the authors at
21                 the e-mail address given above or at
22
23                 Department of Statistics and
24                 Supercomputer Computations Research Institute
25                 Florida State University
26                 Tallahassee, FL 32306.
27
28 See Also:       README          for a brief description
29                 ULTRA.DOC       for a detailed description
30
31 -----------------------------------------------------------------------
32 */ 
33 /*
34    File: ULTRA.C
35
36    This is the ULTRA random number generator written entirely in C.
37
38    This may serve as a model for an assembler version of this routine.
39    The programmer should avoid simply duplicating and instead use the
40    usual assembler features to increase the speed of this routine.
41
42    Especially the subroutine SWB should be replaced by the one
43    machine instruction (usually called subtract-with-borrow) that
44    is available in almost every hardware.
45
46    For people not familiar with 8086 assembler, it may help to
47    consult this when reading the assembler code. This program should
48    be a dropin replacement for the assembler versions, but is about
49    half as fast.
50 */
51
52 /* Slight modifications by Bruno Pincon (4 december 2004) for inclusion 
53    in scilab and nsp:
54
55    1/ in scilab we use only i32bit output ( renamed here fsultra )
56       and  I have deleted the others;
57
58    2/ only one array is now used (swbseed which is renamed
59       swb_state) and the xor with the linear congruential generator
60       is done only just before the output.
61
62    3/ add a var is_init (to say if the generator is initialised)
63
64    4/ add routine to set/get the state
65
66 */
67 #include <math.h>             /* to use floor    */
68 #include "others_generators.h"
69 #include "sciprint.h"
70 #include "localization.h"
71
72 #define N  37           /* size of table        */
73 #define N2 24           /* The shorter lag      */
74
75 static int is_init=0;  
76 static int swb_state[N];          /* state of the swb generator */
77 static int swb_index=N;            /* an index on the swb state */
78 static int swb_flag;               /* the carry flag for the SWB generator */
79 static unsigned int cong_state;   /* state of the congruential generator */
80
81 /* for this generator the state seems completly defined by:
82 swb_state, swb_index, swb_flag (which define the state of the swb generator)
83 cong_state (which defines the state of the congruential generator)
84 */
85
86 /* those are the default for the simple initialisation routine */
87 static  double DEFAULT_SEED1= 1234567.0, DEFAULT_SEED2=7654321.0; 
88
89
90
91
92 /* SWB is the subtract-with-borrow operation which should be one line
93 in assembler code. This should be done by using the hardware s-w-b
94 operation in the SWBfill routine (renamed advance_state_swb here).
95
96 What has been done here is to look at the msb of x, y and z=x-y-c.
97 Using these three bits, one can determine if a borrow bit is needed
98 or not according to the following table:
99
100 msbz=0  msby=0  msby=1          msbz=1  msby=0  msby=1
101
102 msbx=0  0       1               msbx=0  1       1
103 msbx=1  0       0               msbx=1  0       1
104
105 PS: note that the definition is very carefully written because the
106 calls to SWB have y and z as the same memory location, so y must
107 be tested before z is assigned a value.
108 */
109 #define SWB(c,x,y,z) c = (y<0) ? (((z=x-y-c) < 0) || (x>=0)) : (((z=x-y-c) < 0) && (x>=0));
110
111 static void advance_state_swb()
112
113     int i;
114     /*
115     *  The following are the heart of the system and should be
116     *  written is assembler to be as fast as possible. It may even make sense
117     *  to unravel the loop and simply wirte 37 consecutive SWB operations!
118     */
119     for (i=0;  i<N2; i++) 
120         SWB(swb_flag,swb_state[i+N-N2],swb_state[i],swb_state[i]);
121     for (i=N2; i<N;  i++) 
122         SWB(swb_flag,swb_state[i  -N2],swb_state[i],swb_state[i]);
123     swb_index = 0;
124 }
125
126 /* set_state_fsultra_simple initializes the state from 2 integers  
127
128 it defines the constants and fills the swb_state array one bit at
129 a time by taking the leading bit of the xor of a shift register
130 and a congruential sequence. The same congruential generator continues
131 to be used as a mixing generator for the Subtract-with-borrow generator
132 to produce the `ultra' random numbers
133
134 Since this is called just once, speed doesn't matter much and it might
135 be fine to leave this subroutine coded just as it is.
136
137 PS:     there are quick and easy ways to fill this, but since random number
138 generators are really "randomness amplifiers", it is important to
139 start off on the right foot. This is why we take such care here.
140 */
141
142 int set_state_fsultra_simple(double s1, double s2)
143
144     unsigned int shrgx, tidbits=0;
145     int i, j;
146
147     if (    (s1 == floor(s1) && 0.0 <= s1 && s1 <= 4294967295.0)
148         && (s2 == floor(s2) && 0.0 <= s2 && s2 <= 4294967295.0) )
149     {
150         cong_state = ((unsigned int) s1)*2 + 1;
151         shrgx = (unsigned int) s2;
152         for ( i=0 ; i<N ; i++)
153         {
154             for ( j=32 ; j>0 ; j--)
155             { 
156                 cong_state = cong_state * 69069;
157                 shrgx = shrgx ^ (shrgx >> 15);
158                 shrgx = shrgx ^ (shrgx << 17);
159                 tidbits = (tidbits>>1) | (0x80000000 & (cong_state^shrgx));
160             }
161             swb_state[i] = tidbits;
162         }
163         swb_index = 0;
164         swb_flag = 0;
165         advance_state_swb();  /* pour retrouver la même séquence que ds scilab V3.0 */
166         is_init = 1;
167         return 1;
168     }
169     else
170     {
171         sciprint(_("\nBad seed for fsultra, must be integers in [0, 2^32-1]\n"));
172         return 0;
173     }
174 }
175
176 int set_state_fsultra(double *s)
177
178     double try;
179     int i;
180
181     try = s[0];
182     if ( floor(try) != try || try < 0.0  ||  try > (double) N)
183     {
184         sciprint(_("\nThe first component of the fsultra state, must be an int in [0, %d]\n"),N);
185         return 0;
186     }
187     swb_index = (int) try;
188
189     try = s[1];
190     if ( try != 0.0  &&  try != 1.0)
191     {
192         sciprint(_("\nThe second component of the fsultra state, must be 0 or 1\n"));
193         return 0;
194     }
195     swb_flag = (int) try;
196
197     try = s[2];
198     if ( floor(try) != try  ||  try <= 0 ||  try > 4294967295.0 )
199     {
200         sciprint(_("\nThe third component of the fsultra state, must be an int in [1, 2^32-1]\n"));
201         return 0;
202     }
203     cong_state = (unsigned int) try;
204
205     /* no verif here ... */
206     for (i = 0 ; i < N ; i++) 
207         swb_state[i] = (int) (((unsigned int) s[i+3]) & 0xffffffff);
208
209     is_init = 1;
210     return 1;
211 }
212
213
214 /*  to return the state at the scilab level  */
215 void get_state_fsultra(double s[])
216 {
217     int i;
218
219     if ( ! is_init )
220         set_state_fsultra_simple(DEFAULT_SEED1, DEFAULT_SEED2);
221
222     s[0] = (double)  swb_index;
223     s[1] = (double)  swb_flag;
224     s[2] = (double)  cong_state;
225     for (i = 0 ; i < N ; i++) 
226         s[i+3] = (double) (unsigned int) swb_state[i];
227 }
228
229 unsigned int fsultra()
230 {
231     if (swb_index >= N)  /* generate N words at one time */
232     { 
233         if ( ! is_init )
234             set_state_fsultra_simple(DEFAULT_SEED1, DEFAULT_SEED2);
235         else
236             advance_state_swb();
237     }
238     return (swb_state[swb_index++] ^ (cong_state = cong_state * 69069));
239 }