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 448111 : _32to64(ulong a, ulong b) { u64 v = a; return (v<<32)|b; }
36 : static void
37 13724897 : _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 193822188 : block(void)
47 : {
48 193822188 : const int r = 64;
49 193822188 : const int a = 33, b = 26, c = 27, d = 29, s = 53;
50 : u64 t, v, w;
51 193822188 : xorgen_i = (xorgen_i+1)&(r-1);
52 193822188 : t = state[xorgen_i];
53 193822188 : v = state[(xorgen_i+(r-s))&(r-1)]; /* index is (i-s) mod r */
54 193822188 : t ^= t<<a; t ^= t>>b; /* (I + L^a)(I + R^b) */
55 193822188 : v ^= v<<c; v ^= v>>d; /* (I + L^c)(I + R^d) */
56 193822188 : w = t^v;
57 193822188 : return state[xorgen_i] = w; /* update circular array */
58 : }
59 :
60 : /* v > 0 */
61 : static void
62 440916 : init_xor4096i(u64 v)
63 : {
64 440916 : const int r = 64;
65 : int k;
66 :
67 28647759 : for (k = r; k > 0; k--) {/* avoid correlations for close seeds */
68 28206843 : v ^= v<<10; v ^= v>>15; /* recurrence has period 2**64-1 */
69 28206843 : v ^= v<<4; v ^= v>>13;
70 : }
71 28632099 : for (xorgen_w = v, k = 0; k < r; k++) { /* initialise circular array */
72 28191183 : v ^= v<<10; v ^= v>>15;
73 28191183 : v ^= v<<4; v ^= v>>13;
74 28191183 : state[k] = v + (xorgen_w+=weyl);
75 : }
76 : /* discard first 4*r results */
77 111443103 : for (xorgen_i = r-1, k = 4*r; k > 0; k--) (void)block();
78 435190 : }
79 :
80 : void
81 378486 : pari_init_rand(void) { init_xor4096i(1UL); }
82 :
83 : static u64
84 82812382 : rand64(void)
85 : {
86 82812382 : u64 v = block();
87 82825067 : xorgen_w += weyl; /* update Weyl generator */
88 82825067 : return v + (xorgen_w ^ (xorgen_w>>27));
89 : }
90 :
91 : /* One random number uniformly distributed in [0..2**BIL) is returned, where
92 : * BIL = 8*sizeof(ulong) = 32 or 64. */
93 : ulong
94 94984 : pari_rand(void) { return rand64(); }
95 :
96 : void
97 110704 : setrand(GEN x)
98 : {
99 110704 : const int r2 = numberof(state);
100 : long i, lx;
101 : u64 v;
102 : GEN xp;
103 110704 : if (typ(x)!=t_INT) pari_err_TYPE("setrand",x);
104 110704 : if (signe(x) <= 0) pari_err_DOMAIN("setrand","n", "<=", gen_0, x);
105 110697 : lx = lgefint(x);
106 110697 : if (lx == 3) { v = x[2]; init_xor4096i(v); return; }
107 : #ifndef LONG_IS_64BIT
108 6896 : if (lx == 4)
109 : {
110 1 : v = _32to64(*int_W(x,1),*int_W(x,0));
111 1 : init_xor4096i(v); return;
112 : }
113 : #endif
114 48271 : xp = int_LSW(x);
115 : #ifdef LONG_IS_64BIT
116 41376 : if (lx != 2 + r2+2)
117 6 : pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
118 2689050 : for (i = 0; i < r2; i++, xp = int_nextW(xp)) state[i] = *xp;
119 41370 : xorgen_w = *xp; xp = int_nextW(xp);
120 : #else
121 6895 : if (lx != 2 + 2*r2+3)
122 1 : pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
123 448110 : for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
124 441216 : state[i] = _32to64(*int_nextW(xp), *xp);
125 6894 : xorgen_w = _32to64(*int_nextW(xp), *xp); xp = int_nextW(int_nextW(xp));
126 : #endif
127 48264 : xorgen_i = (*xp) & 63;
128 : }
129 :
130 : GEN
131 1458690 : getrand(void)
132 : {
133 1458690 : const int r2 = numberof(state);
134 : GEN x;
135 : ulong *xp;
136 : long i;
137 1458690 : if (xorgen_i < 0) init_xor4096i(1UL);
138 :
139 : #ifdef LONG_IS_64BIT
140 1250530 : x = cgetipos(2+r2+2); xp = (ulong *) int_LSW(x);
141 81284330 : for (i = 0; i < r2; i++, xp = int_nextW(xp)) *xp = state[i];
142 1250530 : *xp = xorgen_w; xp = int_nextW(xp);
143 : #else
144 208160 : x = cgetipos(2+2*r2+3); xp = (ulong *) int_LSW(x);
145 13530400 : for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
146 13322240 : _64to32(state[i], int_nextW(xp), xp);
147 208160 : _64to32(xorgen_w, int_nextW(xp), xp); xp = int_nextW(int_nextW(xp));
148 : #endif
149 1458690 : *xp = xorgen_i? xorgen_i: 64; return x;
150 : }
151 :
152 : /* assume 0 <= k <= BITS_IN_LONG. Return uniform random 0 <= x < (1<<k) */
153 : long
154 23997516 : random_bits(long k) { return rand64() >> (64-k); }
155 :
156 : /********************************************************************/
157 : /* */
158 : /* GENERIC ROUTINES */
159 : /* */
160 : /********************************************************************/
161 :
162 : /* assume n > 0 */
163 : ulong
164 40370307 : random_Fl(ulong n)
165 : {
166 : ulong d;
167 : int shift;
168 : #ifdef LONG_IS_64BIT
169 34032323 : int SHIFT = 0;
170 : #else
171 6337984 : int SHIFT = 32;
172 : #endif
173 :
174 40370307 : if (n == 1) return 0;
175 :
176 40108724 : shift = bfffo(n); /* 2^(BIL-shift) > n >= 2^(BIL-shift-1)*/
177 : /* if N a power of 2, increment shift. No reject */
178 40108724 : if ((n << shift) == HIGHBIT) return rand64() >> (SHIFT+shift+1);
179 : for (;;) {
180 55746237 : d = rand64() >> (SHIFT+shift); /* d < 2^(64-shift) uniformly distributed */
181 : /* reject strategy: proba success = n 2^(shift-64), in [1/2, 1[ */
182 55760326 : if (d < n) return d;
183 : }
184 : }
185 :
186 : /* assume N > 0, see random_Fl() for algorithm. Make sure that 32-bit and
187 : * 64-bit architectures produce the same integers (consuming random bits
188 : * by packets of 64) */
189 : GEN
190 2883010 : randomi(GEN N)
191 : {
192 2883010 : long lx = lgefint(N);
193 : GEN x, d;
194 : int shift;
195 :
196 2883010 : if (lx == 3) return utoi( random_Fl(N[2]) );
197 :
198 254341 : shift = bfffo(*int_MSW(N));
199 : /* if N a power of 2, increment shift */
200 254341 : if (Z_ispow2(N) && ++shift == BITS_IN_LONG) { shift = 0; lx--; }
201 254342 : x = cgetipos(lx);
202 207891 : for (;;) {
203 462233 : GEN y, MSW = int_MSW(x), STOP = MSW;
204 : #ifdef LONG_IS_64BIT
205 1016938 : for (d = int_LSW(x); d != STOP; d = int_nextW(d)) *d = rand64();
206 376515 : *d = rand64() >> shift;
207 : #else
208 85718 : if (!odd(lx)) STOP = int_precW(STOP);
209 : /* STOP points to where MSW would in 64-bit */
210 194497 : for (d = int_LSW(x); d != STOP; d = int_nextW(d))
211 : {
212 108779 : ulong a, b; _64to32(rand64(), &a,&b);
213 108779 : *d = b; d = int_nextW(d);
214 108779 : *d = a;
215 : }
216 : {
217 85718 : ulong a, b; _64to32(rand64() >> shift, &a,&b);
218 85718 : if (d == MSW) /* 32 bits needed */
219 56060 : *d = a;
220 : else
221 : { /* 64 bits needed */
222 29658 : *d = b; d = int_nextW(d);
223 29658 : *d = a;
224 : }
225 : }
226 : #endif
227 462233 : y = int_normalize(x, 0);
228 462233 : if (abscmpii(y, N) < 0) return y;
229 : }
230 : }
231 :
232 : GEN
233 546425 : random_F2x(long d, long vs)
234 : {
235 546425 : ulong db, dl = dvmduBIL(d,&db);
236 546425 : long i, l = 2 + dl + !!db;
237 546425 : GEN y = cgetg(l,t_VECSMALL); y[1] = vs;
238 : #ifdef LONG_IS_64BIT
239 937528 : for (i=2; i<l; i++) uel(y,i) = rand64();
240 : #else
241 78100 : for (i=2; i<l-1; i+=2)
242 : {
243 136 : u64 v = rand64();
244 136 : uel(y,i) = (ulong) v;
245 136 : uel(y,i+1) = (ulong) (v>>32);
246 : }
247 77964 : if (i<l) uel(y,i) = (ulong) rand64();
248 : #endif
249 546425 : if (db) uel(y,l-1) &= ((1UL<<db)-1UL);
250 546425 : return F2x_renormalize(y,l);
251 : }
252 :
253 : GEN
254 59 : random_zv(long n)
255 : {
256 59 : GEN y = cgetg(n+1, t_VECSMALL);
257 : long i;
258 94981 : for (i=1; i<=n; i++) uel(y,i) = pari_rand();
259 59 : return y;
260 : }
261 :
262 : GEN
263 4473 : randomr(long prec)
264 : {
265 : pari_sp av;
266 : long b;
267 : GEN x, y;
268 4473 : if (prec <= 2) return real_0_bit(0);
269 4445 : x = cgetr(prec); av = avma;
270 4445 : b = prec2nbits(prec);
271 4445 : y = randomi(int2n(b));
272 4445 : if (!signe(y)) return real_0_bit(b);
273 4445 : affir(y, x); shiftr_inplace(x, - b);
274 4445 : set_avma(av); return x;
275 : }
276 :
277 : static GEN
278 156772 : polrandom(GEN N) /* assume N!=0 */
279 : {
280 156772 : long i, d = lg(N);
281 156772 : GEN z = leading_coeff(N);
282 156772 : GEN y = cgetg(d,t_POL);
283 156772 : y[1] = evalsigne(1) | evalvarn(varn(N));
284 716086 : for (i=2; i<d; i++) gel(y,i) = genrand(z);
285 156772 : return normalizepol_lg(y,d);
286 : }
287 :
288 : GEN
289 5592018 : genrand(GEN N)
290 : {
291 : GEN z;
292 5592018 : if (!N) return utoi( random_bits(31) );
293 5591339 : switch(typ(N))
294 : {
295 747460 : case t_INT:
296 747460 : switch(signe(N))
297 : {
298 : pari_sp av;
299 : GEN d;
300 287000 : case 1:
301 287000 : return randomi(N);
302 460453 : case -1:
303 460453 : av = avma; N = addiu(N, 1); d = subui(1, shifti(N, 1));
304 460453 : return gc_INT(av, addii(N, randomi(d)));
305 7 : default: pari_err_DOMAIN("random","N","=",gen_0,gen_0);
306 : }
307 7 : case t_REAL:
308 7 : return randomr(realprec(N));
309 680421 : case t_INTMOD:
310 680421 : z = cgetg(3, t_INTMOD);
311 680421 : gel(z,1) = icopy(gel(N,1));
312 680421 : gel(z,2) = randomi(gel(N,1)); return z;
313 184493 : case t_FFELT:
314 184493 : return ffrandom(N);
315 156772 : case t_POL:
316 156772 : if (signe(N)==0) return pol_0(varn(N));
317 156772 : return polrandom(N);
318 3822179 : case t_VEC:
319 3822179 : if (lg(N) == 3)
320 : {
321 3576174 : pari_sp av = avma;
322 3576174 : GEN a = gel(N,1), b = gel(N,2), d;
323 : long s;
324 3576174 : if (typ(a) != t_INT) a = gceil(a);
325 3576174 : if (typ(b) != t_INT) b = gfloor(b);
326 3576174 : if (typ(a) != t_INT || typ(b) != t_INT) pari_err_TYPE("random", N);
327 3576174 : d = subii(b,a); s = signe(d);
328 3576174 : if (s <= 0)
329 : {
330 10066 : if (!s) return gc_GEN(av, a);
331 0 : pari_err_TYPE("random([a,b]) (a > b)", N);
332 : }
333 3566108 : if (lgefint(d) == 3)
334 : {
335 3564687 : ulong z = d[2];
336 3564687 : if (z != ~0UL) return gc_INT(av, addiu(a, random_Fl(z + 1)));
337 : }
338 1421 : return gc_INT(av, addii(a, randomi(addiu(d,1))));
339 : }
340 246005 : return ellrandom(N);
341 7 : default:
342 7 : pari_err_TYPE("genrand",N);
343 : return NULL;/*LCOV_EXCL_LINE*/
344 : }
345 : }
|
