1 #line 2 "../src/kernel/none/mp_indep.c"
2 /* $Id: mp_indep.c 12097 2010-01-29 11:57:49Z bill $
3
4 Copyright (C) 2000 The PARI group.
5
6 This file is part of the PARI/GP package.
7
8 PARI/GP is free software; you can redistribute it and/or modify it under the
9 terms of the GNU General Public License as published by the Free Software
10 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
11 ANY WARRANTY WHATSOEVER.
12
13 Check the License for details. You should have received a copy of it, along
14 with the package; see the file 'COPYING'. If not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
16
17 /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */
18 ulong
invrev(ulong b)19 invrev(ulong b)
20 {
21 static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
22 ulong x = b + tab[b & 7]; /* b^(-1) mod 2^4 */
23
24 /* Newton applied to 1/x - b = 0 */
25 #ifdef LONG_IS_64BIT
26 x = x*(2-b*x); /* one more pass necessary */
27 #endif
28 x = x*(2-b*x);
29 x = x*(2-b*x); return x*(2-b*x);
30 }
31
32 void
affrr(GEN x,GEN y)33 affrr(GEN x, GEN y)
34 {
35 long lx,ly,i;
36
37 y[1] = x[1]; if (!signe(x)) return;
38
39 lx=lg(x); ly=lg(y);
40 if (lx <= ly)
41 {
42 for (i=2; i<lx; i++) y[i]=x[i];
43 for ( ; i<ly; i++) y[i]=0;
44 return;
45 }
46 for (i=2; i<ly; i++) y[i]=x[i];
47 /* lx > ly: round properly */
48 if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
49 }
50
51 GEN
ishiftr(GEN x,long s)52 ishiftr(GEN x, long s)
53 {
54 long ex,lx,n;
55 if (!signe(x)) return gen_0;
56 ex = expo(x) + s; if (ex < 0) return gen_0;
57 lx = lg(x);
58 n=ex - bit_accuracy(lx) + 1;
59 return ishiftr_lg(x, lx, n);
60 }
61
62 GEN
icopy_spec(GEN x,long nx)63 icopy_spec(GEN x, long nx)
64 {
65 GEN z;
66 long i;
67 if (!nx) return gen_0;
68 z = cgeti(nx+2); z[1] = evalsigne(1)|evallgefint(nx+2);
69 for (i=0; i<nx; i++) z[i+2] = x[i];
70 return z;
71 }
72
73 GEN
mului(ulong x,GEN y)74 mului(ulong x, GEN y)
75 {
76 long s = signe(y);
77 GEN z;
78
79 if (!s || !x) return gen_0;
80 z = muluispec(x, y+2, lgefint(y)-2);
81 setsigne(z,s); return z;
82 }
83
84 GEN
mulsi(long x,GEN y)85 mulsi(long x, GEN y)
86 {
87 long s = signe(y);
88 GEN z;
89
90 if (!s || !x) return gen_0;
91 if (x<0) { s = -s; x = -x; }
92 z = muluispec((ulong)x, y+2, lgefint(y)-2);
93 setsigne(z,s); return z;
94 }
95
96 /* assume x > 1, y != 0. Return u * y with sign s */
97 static GEN
mulur_2(ulong x,GEN y,long s)98 mulur_2(ulong x, GEN y, long s)
99 {
100 long m, sh, i, lx = lg(y), e = expo(y);
101 GEN z = cgetr(lx);
102 ulong garde;
103 LOCAL_HIREMAINDER;
104
105 y--; garde = mulll(x,y[lx]);
106 for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
107 z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
108
109 sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
110 if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m);
111 z[1] = evalsigne(s) | evalexpo(m+e); return z;
112 }
113
114 GEN
mulsr(long x,GEN y)115 mulsr(long x, GEN y)
116 {
117 long s, e;
118
119 if (!x) return gen_0;
120 s = signe(y);
121 if (!s)
122 {
123 if (x < 0) x = -x;
124 e = expo(y) + (BITS_IN_LONG-1)-bfffo(x);
125 return real_0_bit(e);
126 }
127 if (x==1) return rcopy(y);
128 if (x==-1) return negr(y);
129 if (x < 0)
130 return mulur_2((ulong)-x, y, -s);
131 else
132 return mulur_2((ulong)x, y, s);
133 }
134
135 GEN
mulur(ulong x,GEN y)136 mulur(ulong x, GEN y)
137 {
138 long s, e;
139
140 if (!x) return gen_0;
141 s = signe(y);
142 if (!s)
143 {
144 e = expo(y) + (BITS_IN_LONG-1)-bfffo(x);
145 return real_0_bit(e);
146 }
147 if (x==1) return rcopy(y);
148 return mulur_2(x, y, s);
149 }
150
151 #ifdef KARAMULR_VARIANT
152 static GEN addshiftw(GEN x, GEN y, long d);
153 static GEN
karamulrr1(GEN y,GEN x,long ly,long lz)154 karamulrr1(GEN y, GEN x, long ly, long lz)
155 {
156 long i, l, lz2 = (lz+2)>>1, lz3 = lz-lz2;
157 GEN lo1, lo2, hi;
158
159 hi = muliispec_mirror(x,y, lz2,lz2);
160 i = lz2; while (i<lz && !x[i]) i++;
161 lo1 = muliispec_mirror(y,x+i, lz2,lz-i);
162 i = lz2; while (i<ly && !y[i]) i++;
163 lo2 = muliispec_mirror(x,y+i, lz2,ly-i);
164 if (signe(lo1))
165 {
166 if (ly!=lz) { lo2 = addshiftw(lo1,lo2,1); lz3++; }
167 else lo2 = addii(lo1,lo2);
168 }
169 l = lgefint(lo2)-(lz3+2);
170 if (l > 0) hi = addiispec(hi+2,lo2+2, lgefint(hi)-2,l);
171 return hi;
172 }
173 #endif
174
175 /* set z <-- x*y, floating point multiplication.
176 * lz = lg(z) = lg(x) <= ly <= lg(y), sz = signe(z). flag = lg(x) < lg(y) */
177 INLINE void
mulrrz_i(GEN z,GEN x,GEN y,long lz,long flag,long sz)178 mulrrz_i(GEN z, GEN x, GEN y, long lz, long flag, long sz)
179 {
180 long ez = expo(x) + expo(y);
181 long i, j, lzz, p1;
182 ulong garde;
183 GEN y1;
184 LOCAL_HIREMAINDER;
185 LOCAL_OVERFLOW;
186
187 if (lz > KARATSUBA_MULR_LIMIT)
188 {
189 pari_sp av = avma;
190 #ifdef KARAMULR_VARIANT
191 GEN hi = karamulrr1(y+2, x+2, lz+flag-2, lz-2);
192 #else
193 GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
194 #endif
195 garde = hi[lz];
196 if (hi[2] < 0)
197 {
198 ez++;
199 for (i=2; i<lz ; i++) z[i] = hi[i];
200 }
201 else
202 {
203 shift_left(z,hi,2,lz-1, garde, 1);
204 garde <<= 1;
205 }
206 if (garde & HIGHBIT)
207 { /* round to nearest */
208 i = lz; do z[--i]++; while (z[i]==0 && i > 1);
209 if (i == 1) { z[2] = HIGHBIT; ez++; }
210 }
211 z[1] = evalsigne(sz)|evalexpo(ez);
212 #if 0
213 {
214 GEN U;
215 KARATSUBA_MULR_LIMIT = 100000;
216 U = mulrr(x, y);
217 KARATSUBA_MULR_LIMIT = 4;
218 if (!gequal(U, z)) pari_err(talker,"toto");
219 }
220 #endif
221 avma = av; return;
222 }
223 if (lz == 3)
224 {
225 if (flag)
226 {
227 (void)mulll(x[2],y[3]);
228 garde = addmul(x[2],y[2]);
229 }
230 else
231 garde = mulll(x[2],y[2]);
232 if (hiremainder & HIGHBIT)
233 {
234 ez++;
235 /* hiremainder < (2^BIL-1)^2 / 2^BIL, hence hiremainder+1 != 0 */
236 if (garde & HIGHBIT) hiremainder++; /* round properly */
237 }
238 else
239 {
240 hiremainder = (hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
241 if (garde & (HIGHBIT-1))
242 {
243 hiremainder++; /* round properly */
244 if (!hiremainder) { hiremainder = HIGHBIT; ez++; }
245 }
246 }
247 z[1] = evalsigne(sz) | evalexpo(ez);
248 z[2] = hiremainder; return;
249 }
250
251 if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
252 lzz=lz-1; p1=x[lzz];
253 if (p1)
254 {
255 (void)mulll(p1,y[3]);
256 garde = addll(addmul(p1,y[2]), garde);
257 z[lzz] = overflow+hiremainder;
258 }
259 else z[lzz]=0;
260 for (j=lz-2, y1=y-j; j>=3; j--)
261 {
262 p1 = x[j]; y1++;
263 if (p1)
264 {
265 (void)mulll(p1,y1[lz+1]);
266 garde = addll(addmul(p1,y1[lz]), garde);
267 for (i=lzz; i>j; i--)
268 {
269 hiremainder += overflow;
270 z[i] = addll(addmul(p1,y1[i]), z[i]);
271 }
272 z[j] = hiremainder+overflow;
273 }
274 else z[j]=0;
275 }
276 p1 = x[2]; y1++;
277 garde = addll(mulll(p1,y1[lz]), garde);
278 for (i=lzz; i>2; i--)
279 {
280 hiremainder += overflow;
281 z[i] = addll(addmul(p1,y1[i]), z[i]);
282 }
283 z[2] = hiremainder+overflow;
284
285 if (z[2] < 0)
286 ez++;
287 else
288 {
289 shift_left(z,z,2,lzz, garde, 1);
290 garde <<= 1;
291 }
292 if (garde & HIGHBIT)
293 { /* round to nearest */
294 i = lz; do z[--i]++; while (z[i]==0 && i > 2);
295 if (z[i] == 0) { z[2] = (long)HIGHBIT; ez++; }
296 }
297 z[1] = evalsigne(sz) | evalexpo(ez);
298 }
299
300 GEN
mulrr(GEN x,GEN y)301 mulrr(GEN x, GEN y)
302 {
303 long flag, ly, lz, sx = signe(x), sy = signe(y);
304 GEN z;
305
306 if (!sx || !sy) return real_0_bit(expo(x) + expo(y));
307 if (sy < 0) sx = -sx;
308 lz = lg(x);
309 ly = lg(y);
310 if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
311 z = cgetr(lz);
312 mulrrz_i(z, x,y, lz,flag, sx);
313 return z;
314 }
315
316 GEN
mulir(GEN x,GEN y)317 mulir(GEN x, GEN y)
318 {
319 long sx = signe(x), sy, lz;
320 GEN z;
321
322 if (!sx) return gen_0;
323 if (!is_bigint(x)) return mulsr(itos(x), y);
324 sy = signe(y);
325 if (!sy) return real_0_bit(expi(x) + expo(y));
326 if (sy < 0) sx = -sx;
327 lz = lg(y); z = cgetr(lz);
328 mulrrz_i(z, itor(x,lz),y, lz,0, sx);
329 avma = (pari_sp)z; return z;
330 }
331
332 /* written by Bruno Haible following an idea of Robert Harley */
333 long
vals(ulong z)334 vals(ulong z)
335 {
336 static char tab[64]={-1,0,1,12,2,6,-1,13,3,-1,7,-1,-1,-1,-1,14,10,4,-1,-1,8,-1,-1,25,-1,-1,-1,-1,-1,21,27,15,31,11,5,-1,-1,-1,-1,-1,9,-1,-1,24,-1,-1,20,26,30,-1,-1,-1,-1,23,-1,19,29,-1,22,18,28,17,16,-1};
337 #ifdef LONG_IS_64BIT
338 long s;
339 #endif
340
341 if (!z) return -1;
342 #ifdef LONG_IS_64BIT
343 if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
344 #endif
345 z |= ~z + 1;
346 z += z << 4;
347 z += z << 6;
348 z ^= z << 16; /* or z -= z<<16 */
349 #ifdef LONG_IS_64BIT
350 return s + tab[(z&0xffffffff)>>26];
351 #else
352 return tab[z>>26];
353 #endif
354 }
355
356 GEN
divsi(long x,GEN y)357 divsi(long x, GEN y)
358 {
359 long p1, s = signe(y);
360 LOCAL_HIREMAINDER;
361
362 if (!s) pari_err(gdiver);
363 if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
364 hiremainder=0; p1=divll(labs(x),y[2]);
365 if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
366 if (s<0) p1 = -p1;
367 return stoi(p1);
368 }
369
370 GEN
divir(GEN x,GEN y)371 divir(GEN x, GEN y)
372 {
373 GEN z;
374 long ly;
375 pari_sp av;
376
377 if (!signe(y)) pari_err(gdiver);
378 if (!signe(x)) return gen_0;
379 ly = lg(y); z = cgetr(ly); av = avma;
380 affrr(divrr(itor(x, ly+1), y), z);
381 avma = av; return z;
382 }
383
384 void
mpdivz(GEN x,GEN y,GEN z)385 mpdivz(GEN x, GEN y, GEN z)
386 {
387 pari_sp av = avma;
388 long tx = typ(x), ty = typ(y);
389 GEN q = tx==t_INT && ty==t_INT? divii(x,y): mpdiv(x,y);
390
391 if (typ(z) == t_REAL) affrr(q, z);
392 else
393 {
394 if (typ(q) == t_REAL) pari_err(gdiver);
395 affii(q, z);
396 }
397 avma = av;
398 }
399
400 GEN
divsr(long x,GEN y)401 divsr(long x, GEN y)
402 {
403 pari_sp av;
404 long ly;
405 GEN z;
406
407 if (!signe(y)) pari_err(gdiver);
408 if (!x) return gen_0;
409 ly = lg(y); z = cgetr(ly); av = avma;
410 affrr(divrr(stor(x,ly+1), y), z);
411 avma = av; return z;
412 }
413
414 GEN
modii(GEN x,GEN y)415 modii(GEN x, GEN y)
416 {
417 switch(signe(x))
418 {
419 case 0: return gen_0;
420 case 1: return remii(x,y);
421 default:
422 {
423 pari_sp av = avma;
424 (void)new_chunk(lgefint(y));
425 x = remii(x,y); avma=av;
426 if (x==gen_0) return x;
427 return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
428 }
429 }
430 }
431
432 void
modiiz(GEN x,GEN y,GEN z)433 modiiz(GEN x, GEN y, GEN z)
434 {
435 const pari_sp av = avma;
436 affii(modii(x,y),z); avma=av;
437 }
438
439 GEN
divrs(GEN x,long y)440 divrs(GEN x, long y)
441 {
442 long i,lx,garde,sh,s=signe(x);
443 GEN z;
444 LOCAL_HIREMAINDER;
445
446 if (!y) pari_err(gdiver);
447 if (!s) return real_0_bit(expo(x) - (BITS_IN_LONG-1)+bfffo(y));
448 if (y<0) { s = -s; y = -y; }
449 if (y==1) { z=rcopy(x); setsigne(z,s); return z; }
450
451 z=cgetr(lx=lg(x)); hiremainder=0;
452 for (i=2; i<lx; i++) z[i] = divll(x[i],y);
453
454 /* we may have hiremainder != 0 ==> garde */
455 garde=divll(0,y); sh=bfffo(z[2]);
456 if (sh) shift_left(z,z, 2,lx-1, garde,sh);
457 z[1] = evalsigne(s) | evalexpo(expo(x)-sh);
458 if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
459 return z;
460 }
461
462 GEN
truedvmdii(GEN x,GEN y,GEN * z)463 truedvmdii(GEN x, GEN y, GEN *z)
464 {
465 pari_sp av;
466 GEN r, q, *gptr[2];
467 if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
468
469 av = avma;
470 q = dvmdii(x,y,&r); /* assume that r is last on stack */
471
472 if (signe(r)>=0)
473 {
474 if (z == ONLY_REM) return gerepileuptoint(av,r);
475 if (z) *z = r; else cgiv(r);
476 return q;
477 }
478
479 if (z == ONLY_REM)
480 { /* r += sign(y) * y, that is |y| */
481 r = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
482 return gerepileuptoint(av, r);
483 }
484 q = addis(q, -signe(y));
485 if (!z) return gerepileuptoint(av, q);
486
487 *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
488 gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(pari_sp)r,gptr,2);
489 return q;
490 }
491 GEN
truedvmdis(GEN x,long y,GEN * z)492 truedvmdis(GEN x, long y, GEN *z)
493 {
494 pari_sp av = avma;
495 long r;
496 GEN q = divis_rem(x,y,&r);
497
498 if (r >= 0)
499 {
500 if (z == ONLY_REM) { avma = av; return utoi(r); }
501 if (z) *z = utoi(r);
502 return q;
503 }
504 if (z == ONLY_REM) { avma = av; return utoi(r + labs(y)); }
505 q = gerepileuptoint(av, addis(q, (y < 0)? 1: -1));
506 if (z) *z = utoi(r + labs(y));
507 return q;
508 }
509
510 /* 2^n = shifti(gen_1, n) */
511 GEN
int2n(long n)512 int2n(long n) {
513 long i, m, d, l;
514 GEN z;
515 if (n < 0) return gen_0;
516 if (n == 0) return gen_1;
517
518 d = n>>TWOPOTBITS_IN_LONG;
519 m = n & (BITS_IN_LONG-1);
520 l = d + 3; z = cgeti(l);
521 z[1] = evalsigne(1) | evallgefint(l);
522 for (i = 2; i < l; i++) z[i] = 0;
523 *int_MSW(z) = 1L << m; return z;
524 }
525 /* To avoid problems when 2^(BIL-1) < n. Overflow cleanly, where int2n
526 * returns gen_0 */
527 GEN
int2u(ulong n)528 int2u(ulong n) {
529 ulong i, m, d, l;
530 GEN z;
531 if (n == 0) return gen_1;
532
533 d = n>>TWOPOTBITS_IN_LONG;
534 m = n & (BITS_IN_LONG-1);
535 l = d + 3; z = cgeti(l);
536 z[1] = evalsigne(1) | evallgefint(l);
537 for (i = 2; i < l; i++) z[i] = 0;
538 *int_MSW(z) = 1L << m; return z;
539 }
540
541 /* actual operations will take place on a+2 and b+2: we strip the codewords */
542 GEN
mulii(GEN a,GEN b)543 mulii(GEN a,GEN b)
544 {
545 long sa,sb;
546 GEN z;
547
548 sa=signe(a); if (!sa) return gen_0;
549 sb=signe(b); if (!sb) return gen_0;
550 if (sb<0) sa = -sa;
551 z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
552 setsigne(z,sa); return z;
553 }
554
555 GEN
remiimul(GEN x,GEN sy)556 remiimul(GEN x, GEN sy)
557 {
558 GEN r, q, y = (GEN)sy[1], invy;
559 long k;
560 pari_sp av = avma;
561
562 k = cmpii(x, y);
563 if (k <= 0) return k? icopy(x): gen_0;
564 invy = (GEN)sy[2];
565 q = mulir(x,invy);
566 q = truncr(q); /* differs from divii(x, y) at most by 1 */
567 r = subii(x, mulii(y,q));
568 if (signe(r) < 0)
569 r = subiispec(y+2,r+2, lgefint(y)-2, lgefint(r)-2); /* y - |r| */
570 else
571 {
572 /* remii(x,y) + y >= r >= remii(x,y) */
573 k = absi_cmp(r, y);
574 if (k >= 0)
575 {
576 if (k == 0) { avma = av; return gen_0; }
577 r = subiispec(r+2, y+2, lgefint(r)-2, lgefint(y)-2);
578 }
579 }
580 #if 0
581 q = subii(r,remii(x,y));
582 if (signe(q))
583 pari_err(talker,"bug in remiimul: x = %Z\ny = %Z\ndif = %Z", x,y,q);
584 #endif
585 return gerepileuptoint(av, r); /* = remii(x, y) */
586 }
587
588 GEN
sqri(GEN a)589 sqri(GEN a) { return sqrispec(a+2, lgefint(a)-2); }
590
591 /* Old cgiv without reference count (which was not used anyway)
592 * Should be a macro.
593 */
594 void
cgiv(GEN x)595 cgiv(GEN x)
596 {
597 if (x == (GEN) avma)
598 avma = (pari_sp) (x+lg(x));
599 }
600
601 /* sqrt()'s result may be off by 1 when a is not representable exactly as a
602 * double [64bit machine] */
603 ulong
usqrtsafe(ulong a)604 usqrtsafe(ulong a)
605 {
606 ulong x = (ulong)sqrt((double)a);
607 #ifdef LONG_IS_64BIT
608 if (x > LOWMASK || x*x > a) x--;
609 #endif
610 return x;
611 }
612
613 /********************************************************************/
614 /** **/
615 /** RANDOM INTEGERS **/
616 /** **/
617 /********************************************************************/
618 static long pari_randseed = 1;
619
620 /* BSD rand gives this: seed = 1103515245*seed + 12345 */
621 /*Return 31 ``random'' bits.*/
622 long
pari_rand31(void)623 pari_rand31(void)
624 {
625 pari_randseed = (1000276549*pari_randseed + 12347) & 0x7fffffff;
626 return pari_randseed;
627 }
628
629 long
setrand(long seed)630 setrand(long seed) { return (pari_randseed = seed); }
631
632 long
getrand(void)633 getrand(void) { return pari_randseed; }
634
635 static ulong
pari_rand(void)636 pari_rand(void)
637 {
638 #define GLUE2(hi, lo) (((hi) << BITS_IN_HALFULONG) | (lo))
639 #if !defined(LONG_IS_64BIT)
640 return GLUE2((pari_rand31()>>12)&LOWMASK,
641 (pari_rand31()>>12)&LOWMASK);
642 #else
643 #define GLUE4(hi1,hi2, lo1,lo2) GLUE2(((hi1)<<16)|(hi2), ((lo1)<<16)|(lo2))
644 # define LOWMASK2 0xffffUL
645 return GLUE4((pari_rand31()>>12)&LOWMASK2,
646 (pari_rand31()>>12)&LOWMASK2,
647 (pari_rand31()>>12)&LOWMASK2,
648 (pari_rand31()>>12)&LOWMASK2);
649 #endif
650 }
651
652 #if 1
653 /* assume N > 0 */
654 GEN
randomi(GEN N)655 randomi(GEN N)
656 {
657 ulong n;
658 long lx, i;
659 GEN x, xMSW, xd, Nd;
660
661 lx = lgefint(N); x = cgeti(lx);
662 x[1] = evalsigne(1) | evallgefint(lx);
663 xMSW = xd = int_MSW(x);
664 for (i=2; i<lx; i++) { *xd = pari_rand(); xd = int_precW(xd); }
665
666 Nd = int_MSW(N); n = (ulong)*Nd;
667 xd = xMSW;
668 if (lx == 3) n--;
669 else
670 for (i=3; i<lx; i++)
671 {
672 xd = int_precW(xd);
673 Nd = int_precW(Nd);
674 if (*xd != *Nd) {
675 if ((ulong)*xd > (ulong)*Nd) n--;
676 break;
677 }
678 }
679 /* MSW needs to be generated between 0 and n */
680 if (n) {
681 LOCAL_HIREMAINDER;
682 (void)mulll((ulong)*xMSW, n + 1);
683 n = hiremainder;
684 }
685 *xMSW = (long)n;
686 if (!n) x = int_normalize(x, 1);
687 return x;
688 }
689 #else
690 /* assume N > 0 */
691 GEN
randomi(GEN N)692 randomi(GEN N)
693 {
694 long i, lx = lgefint(N), shift = bfffo(*int_MSW(N));
695 GEN x = cgeti(lx), xMSW, xd;
696
697 x[1] = evalsigne(1) | evallgefint(lx);
698 xMSW = int_MSW(x);
699 for (xd = xMSW;;) {
700 for (i=2; i<lx; i++) { *xd = pari_rand(); xd = int_precW(xd); }
701 *xMSW = ((ulong)*xMSW) >> shift;
702 x = int_normalize(x, 0);
703 if (absi_cmp(x, N) < 0) return x;
704 }
705 }
706 #endif
707
708 /********************************************************************/
709 /** **/
710 /** EXPONENT / CONVERSION t_REAL --> double **/
711 /** **/
712 /********************************************************************/
713
714 #ifdef LONG_IS_64BIT
715 long
dblexpo(double x)716 dblexpo(double x)
717 {
718 union { double f; ulong i; } fi;
719 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
720 const int exp_mid = 0x3ff;/* exponent bias */
721
722 if (x==0.) return -exp_mid;
723 fi.f = x;
724 return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
725 }
726
727 ulong
dblmantissa(double x)728 dblmantissa(double x)
729 {
730 union { double f; ulong i; } fi;
731 const int expo_len = 11; /* number of bits of exponent */
732
733 if (x==0.) return 0;
734 fi.f = x;
735 return (fi.i << expo_len) | HIGHBIT;
736 }
737
738 GEN
dbltor(double x)739 dbltor(double x)
740 {
741 GEN z;
742 long e;
743 union { double f; ulong i; } fi;
744 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
745 const int exp_mid = 0x3ff;/* exponent bias */
746 const int expo_len = 11; /* number of bits of exponent */
747
748 if (x==0.) return real_0_bit(-exp_mid);
749 fi.f = x; z = cgetr(DEFAULTPREC);
750 {
751 const ulong a = fi.i;
752 ulong A;
753 e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
754 if (e == exp_mid+1) pari_err(talker, "NaN or Infinity in dbltor");
755 A = a << expo_len;
756 if (e == -exp_mid)
757 { /* unnormalized values */
758 int sh = bfffo(A);
759 e -= sh-1;
760 z[2] = A << sh;
761 }
762 else
763 z[2] = HIGHBIT | A;
764 z[1] = evalexpo(e) | evalsigne(x<0? -1: 1);
765 }
766 return z;
767 }
768
769 double
rtodbl(GEN x)770 rtodbl(GEN x)
771 {
772 long ex,s=signe(x);
773 ulong a;
774 union { double f; ulong i; } fi;
775 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
776 const int exp_mid = 0x3ff;/* exponent bias */
777 const int expo_len = 11; /* number of bits of exponent */
778
779 if (typ(x)==t_INT && !s) return 0.0;
780 if (typ(x)!=t_REAL) pari_err(typeer,"rtodbl");
781 if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
782
783 /* start by rounding to closest */
784 a = (x[2] & (HIGHBIT-1)) + 0x400;
785 if (a & HIGHBIT) { ex++; a=0; }
786 if (ex >= exp_mid) pari_err(rtodber);
787 fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
788 if (s<0) fi.i |= HIGHBIT;
789 return fi.f;
790 }
791
792 #else /* LONG_IS_64BIT */
793
794 #if PARI_DOUBLE_FORMAT == 1
795 # define INDEX0 1
796 # define INDEX1 0
797 #elif PARI_DOUBLE_FORMAT == 0
798 # define INDEX0 0
799 # define INDEX1 1
800 #endif
801
802 long
dblexpo(double x)803 dblexpo(double x)
804 {
805 union { double f; ulong i[2]; } fi;
806 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
807 const int exp_mid = 0x3ff;/* exponent bias */
808 const int shift = mant_len-32;
809
810 if (x==0.) return -exp_mid;
811 fi.f = x;
812 {
813 const ulong a = fi.i[INDEX0];
814 return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
815 }
816 }
817
818 ulong
dblmantissa(double x)819 dblmantissa(double x)
820 {
821 union { double f; ulong i[2]; } fi;
822 const int expo_len = 11; /* number of bits of exponent */
823
824 if (x==0.) return 0;
825 fi.f = x;
826 {
827 const ulong a = fi.i[INDEX0];
828 const ulong b = fi.i[INDEX1];
829 return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
830 }
831 }
832
833 GEN
dbltor(double x)834 dbltor(double x)
835 {
836 GEN z;
837 long e;
838 union { double f; ulong i[2]; } fi;
839 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
840 const int exp_mid = 0x3ff;/* exponent bias */
841 const int expo_len = 11; /* number of bits of exponent */
842 const int shift = mant_len-32;
843
844 if (x==0.) return real_0_bit(-exp_mid);
845 fi.f = x; z = cgetr(DEFAULTPREC);
846 {
847 const ulong a = fi.i[INDEX0];
848 const ulong b = fi.i[INDEX1];
849 ulong A, B;
850 e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
851 if (e == exp_mid+1) pari_err(talker, "NaN or Infinity in dbltor");
852 A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
853 B = b << expo_len;
854 if (e == -exp_mid)
855 { /* unnormalized values */
856 int sh;
857 if (A)
858 {
859 sh = bfffo(A);
860 e -= sh-1;
861 z[2] = (A << sh) | (B >> (32-sh));
862 z[3] = B << sh;
863 }
864 else
865 {
866 sh = bfffo(B); /* B != 0 */
867 e -= sh-1 + 32;
868 z[2] = B << sh;
869 z[3] = 0;
870 }
871 }
872 else
873 {
874 z[3] = B;
875 z[2] = HIGHBIT | A;
876 }
877 z[1] = evalexpo(e) | evalsigne(x<0? -1: 1);
878 }
879 return z;
880 }
881
882 double
rtodbl(GEN x)883 rtodbl(GEN x)
884 {
885 long ex,s=signe(x),lx=lg(x);
886 ulong a,b,k;
887 union { double f; ulong i[2]; } fi;
888 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
889 const int exp_mid = 0x3ff;/* exponent bias */
890 const int expo_len = 11; /* number of bits of exponent */
891 const int shift = mant_len-32;
892
893 if (typ(x)==t_INT && !s) return 0.0;
894 if (typ(x)!=t_REAL) pari_err(typeer,"rtodbl");
895 if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
896
897 /* start by rounding to closest */
898 a = x[2] & (HIGHBIT-1);
899 if (lx > 3)
900 {
901 b = x[3] + 0x400UL; if (b < 0x400UL) a++;
902 if (a & HIGHBIT) { ex++; a=0; }
903 }
904 else b = 0;
905 if (ex >= exp_mid) pari_err(rtodber);
906 ex += exp_mid;
907 k = (a >> expo_len) | (ex << shift);
908 if (s<0) k |= HIGHBIT;
909 fi.i[INDEX0] = k;
910 fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
911 return fi.f;
912 }
913 #endif /* LONG_IS_64BIT */
914
915