1 /* This program by D E Knuth is in the public domain and freely copyable.
2 * It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
3 * (or in the errata to the 2nd edition --- see
4 * http://www-cs-faculty.stanford.edu/~knuth/taocp.html
5 * in the changes to Volume 2 on pages 171 and following). */
6
7 /* N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
8 included here; there's no backwards compatibility with the original. */
9
10 /* This version also adopts Brendan McKay's suggestion to
11 accommodate naive users who forget to call ranf_start(seed). */
12
13 /* If you find any bugs, please report them immediately to
14 * taocp@cs.stanford.edu
15 * (and you will be rewarded if the bug is genuine). Thanks! */
16
17 /************ see the book for explanations and caveats! *******************/
18 /************ in particular, you need two's complement arithmetic **********/
19
20 /* This version has been changed by Jon Nordby to allow to create multiple
21 * independent generator objects. All changes made to this file are considered
22 * to be in the public domain. */
23
24 #include "config.h"
25
26 #include "rng-double.h"
27
28 #include <stdlib.h>
29
30 /* the following routines are adapted from exercise 3.6--15 */
31 /* after calling ranf_start, get new randoms by, e.g., "x=ranf_arr_next()" */
32
33 #if 0
34 /* original settings */
35 #define QUALITY 1009 /* recommended quality level for high-res use */
36 #define TT 70 /* guaranteed separation between streams */
37 #define KK 100 /* the long lag */
38 #define LL 37 /* the short lag */
39 #else
40 /* low quality settings, seems to work for MyPaint */
41 /* (Disclaimer: I don't understand what those numbers do, I just reduced them. --maxy) */
42 #define QUALITY 19
43 #define TT 7
44 #define KK 10
45 #define LL 7
46 #endif
47
48 #define is_odd(s) ((s)&1)
49 #define mod_sum(x,y) (((x)+(y))-(int)((x)+(y))) /* (x+y) mod 1.0 */
50
51 const double ranf_arr_dummy=-1.0;
52 const double ranf_arr_started=-1.0;
53
54 struct RngDouble {
55 double ran_u[KK]; /* the generator state */
56 double ranf_arr_buf[QUALITY];
57 double *ranf_arr_ptr; /* the next random fraction, or -1 */
58 };
59
60 void
rng_double_get_array(RngDouble * self,double aa[],int n)61 rng_double_get_array(RngDouble *self, double aa[], int n)
62 {
63 register int i,j;
64 for (j=0;j<KK;j++) aa[j]=self->ran_u[j];
65 for (;j<n;j++) aa[j]=mod_sum(aa[j-KK],aa[j-LL]);
66 for (i=0;i<LL;i++,j++) self->ran_u[i]=mod_sum(aa[j-KK],aa[j-LL]);
67 for (;i<KK;i++,j++) self->ran_u[i]=mod_sum(aa[j-KK],self->ran_u[i-LL]);
68 }
69
70
71 RngDouble *
rng_double_new(long seed)72 rng_double_new(long seed)
73 {
74 RngDouble *self = (RngDouble *)malloc(sizeof(RngDouble));
75
76 self->ranf_arr_ptr=(double *)&ranf_arr_dummy;
77
78 rng_double_set_seed(self, seed);
79
80 return self;
81 }
82
83 void
rng_double_free(RngDouble * self)84 rng_double_free(RngDouble *self)
85 {
86 free(self);
87 }
88
89 void
rng_double_set_seed(RngDouble * self,long seed)90 rng_double_set_seed(RngDouble *self, long seed)
91 {
92 register int t,s,j;
93 double u[KK+KK-1];
94 double ulp=(1.0/(1L<<30))/(1L<<22); /* 2 to the -52 */
95 double ss=2.0*ulp*((seed&0x3fffffff)+2);
96
97 for (j=0;j<KK;j++) {
98 u[j]=ss; /* bootstrap the buffer */
99 ss+=ss; if (ss>=1.0) ss-=1.0-2*ulp; /* cyclic shift of 51 bits */
100 }
101 u[1]+=ulp; /* make u[1] (and only u[1]) "odd" */
102 for (s=seed&0x3fffffff,t=TT-1; t; ) {
103 for (j=KK-1;j>0;j--)
104 u[j+j]=u[j],u[j+j-1]=0.0; /* "square" */
105 for (j=KK+KK-2;j>=KK;j--) {
106 u[j-(KK-LL)]=mod_sum(u[j-(KK-LL)],u[j]);
107 u[j-KK]=mod_sum(u[j-KK],u[j]);
108 }
109 if (is_odd(s)) { /* "multiply by z" */
110 for (j=KK;j>0;j--) u[j]=u[j-1];
111 u[0]=u[KK]; /* shift the buffer cyclically */
112 u[LL]=mod_sum(u[LL],u[KK]);
113 }
114 if (s) s>>=1; else t--;
115 }
116 for (j=0;j<LL;j++) self->ran_u[j+KK-LL]=u[j];
117 for (;j<KK;j++) self->ran_u[j-LL]=u[j];
118 for (j=0;j<10;j++) rng_double_get_array(self, u,KK+KK-1); /* warm things up */
119 self->ranf_arr_ptr=(double *)&ranf_arr_started;
120 }
121
122
123 double
rng_double_cycle(RngDouble * self)124 rng_double_cycle(RngDouble *self)
125 {
126 rng_double_get_array(self, self->ranf_arr_buf, QUALITY);
127 self->ranf_arr_buf[KK]=-1;
128 self->ranf_arr_ptr=self->ranf_arr_buf+1;
129 return self->ranf_arr_buf[0];
130 }
131
132 double
rng_double_next(RngDouble * self)133 rng_double_next(RngDouble *self)
134 {
135 return ((*self->ranf_arr_ptr>=0) ? *(self->ranf_arr_ptr)++ : rng_double_cycle(self));
136 }
137