Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 : /********************************************************************/
15 : /* */
16 : /* PSEUDO-RANDOM INTEGERS */
17 : /* */
18 : /********************************************************************/
19 : #include "pari.h"
20 : #include "paripriv.h"
21 : /********************************************************************/
22 : /* XORGEN (Richard P. Brent) */
23 : /* http://wwwmaths.anu.edu.au/~brent/random.html */
24 : /* (initial adaptation to PARI/GP by Randall Rathbun) */
25 : /********************************************************************/
26 : /* Adapted from xorgens.c version 3.04, Richard P. Brent, 20060628 (GPL).
27 : * 32-bit or 64-bit integer random number generator with period at
28 : * least 2**4096-1. It is assumed that "ulong" is a 32-bit or 64-bit integer */
29 :
30 : #ifdef LONG_IS_64BIT
31 : typedef ulong u64;
32 : #else
33 : typedef unsigned long long u64;
34 : static u64
35 448241 : _32to64(ulong a, ulong b) { u64 v = a; return (v<<32)|b; }
36 : static void
37 13017651 : _64to32(u64 v, ulong *a, ulong *b) { *a = v>>32; *b = v&0xFFFFFFFF; }
38 : #endif
39 : static THREAD u64 state[64];
40 : static THREAD u64 xorgen_w;
41 : static THREAD int xorgen_i;
42 : /* weyl = odd approximation to 2^64*(sqrt(5)-1)/2. */
43 : static const u64 weyl = (((u64)0x61c88646U)<<32)|((u64)0x80b583ebU);
44 :
45 : static u64
46 168300820 : block(void)
47 : {
48 168300820 : const int r = 64;
49 168300820 : const int a = 33, b = 26, c = 27, d = 29, s = 53;
50 : u64 t, v, w;
51 168300820 : xorgen_i = (xorgen_i+1)&(r-1);
52 168300820 : t = state[xorgen_i];
53 168300820 : v = state[(xorgen_i+(r-s))&(r-1)]; /* index is (i-s) mod r */
54 168300820 : t ^= t<<a; t ^= t>>b; /* (I + L^a)(I + R^b) */
55 168300820 : v ^= v<<c; v ^= v>>d; /* (I + L^c)(I + R^d) */
56 168300820 : w = t^v;
57 168300820 : return state[xorgen_i] = w; /* update circular array */
58 : }
59 :
60 : /* v > 0 */
61 : static void
62 381860 : init_xor4096i(u64 v)
63 : {
64 381860 : const int r = 64;
65 : int k;
66 :
67 24804144 : for (k = r; k > 0; k--) {/* avoid correlations for close seeds */
68 24422284 : v ^= v<<10; v ^= v>>15; /* recurrence has period 2**64-1 */
69 24422284 : v ^= v<<4; v ^= v>>13;
70 : }
71 24790296 : for (xorgen_w = v, k = 0; k < r; k++) { /* initialise circular array */
72 24408436 : v ^= v<<10; v ^= v>>15;
73 24408436 : v ^= v<<4; v ^= v>>13;
74 24408436 : state[k] = v + (xorgen_w+=weyl);
75 : }
76 : /* discard first 4*r results */
77 96205435 : for (xorgen_i = r-1, k = 4*r; k > 0; k--) (void)block();
78 378621 : }
79 :
80 326435 : void pari_init_rand(void) { init_xor4096i(1UL); }
81 :
82 : static u64
83 72472650 : rand64(void)
84 : {
85 72472650 : u64 v = block();
86 72485074 : xorgen_w += weyl; /* update Weyl generator */
87 72485074 : return v + (xorgen_w ^ (xorgen_w>>27));
88 : }
89 :
90 : /* One random number uniformly distributed in [0..2**BIL) is returned, where
91 : * BIL = 8*sizeof(ulong) = 32 or 64. */
92 : ulong
93 71364 : pari_rand(void) { return rand64(); }
94 :
95 : void
96 103718 : setrand(GEN x)
97 : {
98 103718 : const int r2 = numberof(state);
99 : long i, lx;
100 : u64 v;
101 : GEN xp;
102 103718 : if (typ(x)!=t_INT) pari_err_TYPE("setrand",x);
103 103718 : if (signe(x) <= 0) pari_err_DOMAIN("setrand","n", "<=", gen_0, x);
104 103711 : lx = lgefint(x);
105 103711 : if (lx == 3) { v = x[2]; init_xor4096i(v); return; }
106 : #ifndef LONG_IS_64BIT
107 6898 : if (lx == 4)
108 : {
109 1 : v = _32to64(*int_W(x,1),*int_W(x,0));
110 1 : init_xor4096i(v); return;
111 : }
112 : #endif
113 48285 : xp = int_LSW(x);
114 : #ifdef LONG_IS_64BIT
115 41388 : if (lx != 2 + r2+2)
116 6 : pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
117 2689830 : for (i = 0; i < r2; i++, xp = int_nextW(xp)) state[i] = *xp;
118 41382 : xorgen_w = *xp; xp = int_nextW(xp);
119 : #else
120 6897 : if (lx != 2 + 2*r2+3)
121 1 : pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
122 448240 : for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
123 441344 : state[i] = _32to64(*int_nextW(xp), *xp);
124 6896 : xorgen_w = _32to64(*int_nextW(xp), *xp); xp = int_nextW(int_nextW(xp));
125 : #endif
126 48278 : xorgen_i = (*xp) & 63;
127 : }
128 :
129 : GEN
130 1386728 : getrand(void)
131 : {
132 1386728 : const int r2 = numberof(state);
133 : GEN x;
134 : ulong *xp;
135 : long i;
136 1386728 : if (xorgen_i < 0) init_xor4096i(1UL);
137 :
138 : #ifdef LONG_IS_64BIT
139 1188846 : x = cgetipos(2+r2+2); xp = (ulong *) int_LSW(x);
140 77274905 : for (i = 0; i < r2; i++, xp = int_nextW(xp)) *xp = state[i];
141 1188847 : *xp = xorgen_w; xp = int_nextW(xp);
142 : #else
143 197882 : x = cgetipos(2+2*r2+3); xp = (ulong *) int_LSW(x);
144 12862330 : for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
145 12664448 : _64to32(state[i], int_nextW(xp), xp);
146 197882 : _64to32(xorgen_w, int_nextW(xp), xp); xp = int_nextW(int_nextW(xp));
147 : #endif
148 1386729 : *xp = xorgen_i? xorgen_i: 64; return x;
149 : }
150 :
151 : /* assume 0 <= k <= BITS_IN_LONG. Return uniform random 0 <= x < (1<<k) */
152 : long
153 17090866 : random_bits(long k) { return rand64() >> (64-k); }
154 :
155 : /********************************************************************/
156 : /* */
157 : /* GENERIC ROUTINES */
158 : /* */
159 : /********************************************************************/
160 :
161 : /* assume n > 0 */
162 : ulong
163 38131449 : random_Fl(ulong n)
164 : {
165 : ulong d;
166 : int shift;
167 : #ifdef LONG_IS_64BIT
168 32119235 : int SHIFT = 0;
169 : #else
170 6012214 : int SHIFT = 32;
171 : #endif
172 :
173 38131449 : if (n == 1) return 0;
174 :
175 37859842 : shift = bfffo(n); /* 2^(BIL-shift) > n >= 2^(BIL-shift-1)*/
176 : /* if N a power of 2, increment shift. No reject */
177 37859842 : if ((n << shift) == HIGHBIT) return rand64() >> (SHIFT+shift+1);
178 : for (;;) {
179 52590798 : d = rand64() >> (SHIFT+shift); /* d < 2^(64-shift) uniformly distributed */
180 : /* reject strategy: proba success = n 2^(shift-64), in [1/2, 1[ */
181 52601059 : if (d < n) return d;
182 : }
183 : }
184 :
185 : /* assume N > 0, see random_Fl() for algorithm. Make sure that 32-bit and
186 : * 64-bit architectures produce the same integers (consuming random bits
187 : * by packets of 64) */
188 : GEN
189 6389898 : randomi(GEN N)
190 : {
191 6389898 : long lx = lgefint(N);
192 : GEN x, d;
193 : int shift;
194 :
195 6389898 : if (lx == 3) return utoi( random_Fl(N[2]) );
196 :
197 193812 : shift = bfffo(*int_MSW(N));
198 : /* if N a power of 2, increment shift */
199 193812 : if (Z_ispow2(N) && ++shift == BITS_IN_LONG) { shift = 0; lx--; }
200 193812 : x = cgetipos(lx);
201 145660 : for (;;) {
202 339472 : GEN y, MSW = int_MSW(x), STOP = MSW;
203 : #ifdef LONG_IS_64BIT
204 809378 : for (d = int_LSW(x); d != STOP; d = int_nextW(d)) *d = rand64();
205 272282 : *d = rand64() >> shift;
206 : #else
207 67190 : if (!odd(lx)) STOP = int_precW(STOP);
208 : /* STOP points to where MSW would in 64-bit */
209 155321 : for (d = int_LSW(x); d != STOP; d = int_nextW(d))
210 : {
211 88131 : ulong a, b; _64to32(rand64(), &a,&b);
212 88131 : *d = b; d = int_nextW(d);
213 88131 : *d = a;
214 : }
215 : {
216 67190 : ulong a, b; _64to32(rand64() >> shift, &a,&b);
217 67190 : if (d == MSW) /* 32 bits needed */
218 37258 : *d = a;
219 : else
220 : { /* 64 bits needed */
221 29932 : *d = b; d = int_nextW(d);
222 29932 : *d = a;
223 : }
224 : }
225 : #endif
226 339472 : y = int_normalize(x, 0);
227 339472 : if (abscmpii(y, N) < 0) return y;
228 : }
229 : }
230 :
231 : GEN
232 543290 : random_F2x(long d, long vs)
233 : {
234 543290 : ulong db, dl = dvmduBIL(d,&db);
235 543290 : long i, l = 2 + dl + !!db;
236 543290 : GEN y = cgetg(l,t_VECSMALL); y[1] = vs;
237 : #ifdef LONG_IS_64BIT
238 931782 : for (i=2; i<l; i++) uel(y,i) = rand64();
239 : #else
240 77838 : for (i=2; i<l-1; i+=2)
241 : {
242 135 : u64 v = rand64();
243 135 : uel(y,i) = (ulong) v;
244 135 : uel(y,i+1) = (ulong) (v>>32);
245 : }
246 77703 : if (i<l) uel(y,i) = (ulong) rand64();
247 : #endif
248 543291 : if (db) uel(y,l-1) &= ((1UL<<db)-1UL);
249 543291 : return F2x_renormalize(y,l);
250 : }
251 :
252 : GEN
253 49 : random_zv(long n)
254 : {
255 49 : GEN y = cgetg(n+1, t_VECSMALL);
256 : long i;
257 71351 : for (i=1; i<=n; i++) uel(y,i) = pari_rand();
258 49 : return y;
259 : }
260 :
261 : GEN
262 4473 : randomr(long prec)
263 : {
264 : pari_sp av;
265 : long b;
266 : GEN x, y;
267 4473 : if (prec <= 2) return real_0_bit(0);
268 4445 : x = cgetr(prec); av = avma;
269 4445 : b = prec2nbits(prec);
270 4445 : y = randomi(int2n(b));
271 4445 : if (!signe(y)) return real_0_bit(b);
272 4445 : affir(y, x); shiftr_inplace(x, - b);
273 4445 : set_avma(av); return x;
274 : }
275 :
276 : static GEN
277 156478 : polrandom(GEN N) /* assume N!=0 */
278 : {
279 156478 : long i, d = lg(N);
280 156478 : GEN z = leading_coeff(N);
281 156478 : GEN y = cgetg(d,t_POL);
282 156478 : y[1] = evalsigne(1) | evalvarn(varn(N));
283 714532 : for (i=2; i<d; i++) gel(y,i) = genrand(z);
284 156478 : return normalizepol_lg(y,d);
285 : }
286 :
287 : GEN
288 5592020 : genrand(GEN N)
289 : {
290 : GEN z;
291 5592020 : if (!N) return utoi( random_bits(31) );
292 5591341 : switch(typ(N))
293 : {
294 744744 : case t_INT:
295 744744 : switch(signe(N))
296 : {
297 : pari_sp av;
298 : GEN d;
299 285453 : case 1:
300 285453 : return randomi(N);
301 459284 : case -1:
302 459284 : av = avma; N = addiu(N, 1); d = subui(1, shifti(N, 1));
303 459284 : return gerepileuptoint(av, addii(N, randomi(d)));
304 7 : default: pari_err_DOMAIN("random","N","=",gen_0,gen_0);
305 : }
306 7 : case t_REAL:
307 7 : return randomr(realprec(N));
308 680421 : case t_INTMOD:
309 680421 : z = cgetg(3, t_INTMOD);
310 680421 : gel(z,1) = icopy(gel(N,1));
311 680421 : gel(z,2) = randomi(gel(N,1)); return z;
312 187901 : case t_FFELT:
313 187901 : return ffrandom(N);
314 156478 : case t_POL:
315 156478 : if (signe(N)==0) return pol_0(varn(N));
316 156478 : return polrandom(N);
317 3821783 : case t_VEC:
318 3821783 : if (lg(N) == 3)
319 : {
320 3576104 : pari_sp av = avma;
321 3576104 : GEN a = gel(N,1), b = gel(N,2), d;
322 3576104 : if (typ(a) != t_INT) a = gceil(a);
323 3576104 : if (typ(b) != t_INT) b = gfloor(b);
324 3576104 : if (typ(a) != t_INT || typ(b) != t_INT) pari_err_TYPE("random", N);
325 3576104 : d = subii(b,a);
326 3576104 : if (signe(d) < 0) pari_err_TYPE("random([a,b]) (a > b)", N);
327 3576104 : return gerepileuptoint(av, addii(a, randomi(addiu(d,1))));
328 : }
329 245679 : return ellrandom(N);
330 7 : default:
331 7 : pari_err_TYPE("genrand",N);
332 : return NULL;/*LCOV_EXCL_LINE*/
333 : }
334 : }
|
