randlib: Fixed bug #8560
[scilab.git] / scilab / modules / randlib / src / c / kiss.c
1 /* 
2  * 
3  * Copyright (C) 2010 - DIGITEO - Michael Baudin
4  * Copyright (C) 2004 - Bruno Pincon
5  * Copyright (C) 1999 - G. Marsaglia
6  * 
7  *   PURPOSE
8  *      the kiss generator of G. Marsaglia
9  *      generate random integers (uint) in [0, 2^32 - 1]
10  *      the state is given by 4 integers (z, w, jsr, jcong)
11  *
12  *   NOTES
13  *      The code was given by G. Marsaglia at the end of  a
14  *      thread  concerning  RNG  in C in several newsgroups
15  *      (whom sci.math.num-analysis) "My offer of RNG's for
16  *      C  was an invitation to dance..."
17  *
18  *      Slight modifications by Bruno Pincon for inclusion in
19  *      Scilab (added set/get state routines)
20  *
21  *      kiss is made of combinaison of severals  others but  
22  *      they  are not interfaced at the scilab level.
23  *
24  *      Need that it is assumed that the 
25  *         unsigned int arithmetic is the classic 32 bits 
26  *         unsigned arithmetic modulo 2^32 (ie all is exact
27  *         modulo 2^32) 
28  *
29  *      License: This code is released under public domain.
30  *      (Compatible with CeCILL)
31  */
32 #include <math.h>             /* to use floor    */
33
34 #include "sciprint.h"
35 #include "others_generators.h"
36 #include "localization.h"
37
38 /* The Marsaglia 's macros : */
39 #define znew  (z=36969*(z&65535)+(z>>16))
40 #define wnew  (w=18000*(w&65535)+(w>>16))
41 #define MWC   ((znew<<16)+wnew )
42 #define CONG  (jcong=69069*jcong+1234567)
43 #define SHR3  (jsr^=(jsr<<17), jsr^=(jsr>>13), jsr^=(jsr<<5))
44 #define KISS  ((MWC^CONG)+SHR3)
45
46 /*  the kiss 's state  (any int in [0,2^32-1] are OK ?) */
47 static unsigned int z=362436069, w=521288629, jsr=123456789, jcong=380116160;
48
49 unsigned int kiss()
50 {
51     return ( KISS );
52 }
53
54 int set_state_kiss(double g1, double g2, double g3, double g4)
55 {
56     if (g1 == floor(g1) && g2 == floor(g2) && 
57         g3 == floor(g3) && g4 == floor(g4) &&  
58         0.0 <= g1 && g1 <= 4294967295.0 &&
59         0.0 <= g2 && g2 <= 4294967295.0 &&
60         0.0 <= g3 && g3 <= 4294967295.0 &&
61         0.0 <= g4 && g4 <= 4294967295.0 )
62     {
63         z = (unsigned int) g1;
64         w = (unsigned int) g2;
65         jsr = (unsigned int) g3;
66         jcong = (unsigned int) g4;
67         return ( 1 );
68     }
69     else
70     {
71         sciprint(_("Bad seeds for kiss, must be integers in [0,2^32-1]\n"));
72         return ( 0 );
73     }
74 }
75
76 void get_state_kiss(double g[])
77 {
78     g[0] = (double) z;
79     g[1] = (double) w;
80     g[2] = (double) jsr;
81     g[3] = (double) jcong;
82 }
83