1 #line 2 "../src/kernel/none/mp_indep.c"
2 /* Copyright (C) 2000 The PARI group.
3
4 This file is part of the PARI/GP package.
5
6 PARI/GP is free software; you can redistribute it and/or modify it under the
7 terms of the GNU General Public License as published by the Free Software
8 Foundation; either version 2 of the License, or (at your option) any later
9 version. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
15
16 /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */
17 ulong
invmod2BIL(ulong b)18 invmod2BIL(ulong b)
19 {
20 static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
21 ulong x = b + tab[b & 7]; /* b^(-1) mod 2^4 */
22
23 /* Newton applied to 1/x - b = 0 */
24 #ifdef LONG_IS_64BIT
25 x = x*(2-b*x); /* one more pass necessary */
26 #endif
27 x = x*(2-b*x);
28 x = x*(2-b*x); return x*(2-b*x);
29 }
30
31 void
affrr(GEN x,GEN y)32 affrr(GEN x, GEN y)
33 {
34 long i, lx, ly = lg(y);
35 if (!signe(x))
36 {
37 y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
38 return;
39 }
40 y[1] = x[1]; lx = lg(x);
41 if (lx <= ly)
42 {
43 for (i=2; i<lx; i++) y[i]=x[i];
44 for ( ; i<ly; i++) y[i]=0;
45 return;
46 }
47 for (i=2; i<ly; i++) y[i]=x[i];
48 /* lx > ly: round properly */
49 if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
50 }
51
52 GEN
trunc2nr(GEN x,long n)53 trunc2nr(GEN x, long n)
54 {
55 long ex;
56 if (!signe(x)) return gen_0;
57 ex = expo(x) + n; if (ex < 0) return gen_0;
58 return mantissa2nr(x, ex - bit_prec(x) + 1);
59 }
60
61 /* x a t_REAL, x = i/2^e, i a t_INT */
62 GEN
mantissa_real(GEN x,long * e)63 mantissa_real(GEN x, long *e)
64 {
65 *e = bit_prec(x)-1-expo(x);
66 return mantissa2nr(x, 0);
67 }
68
69 GEN
mului(ulong x,GEN y)70 mului(ulong x, GEN y)
71 {
72 long s = signe(y);
73 GEN z;
74
75 if (!s || !x) return gen_0;
76 z = muluispec(x, y+2, lgefint(y)-2);
77 setsigne(z,s); return z;
78 }
79
80 GEN
mulsi(long x,GEN y)81 mulsi(long x, GEN y)
82 {
83 long s = signe(y);
84 GEN z;
85
86 if (!s || !x) return gen_0;
87 if (x<0) { s = -s; x = -x; }
88 z = muluispec((ulong)x, y+2, lgefint(y)-2);
89 setsigne(z,s); return z;
90 }
91
92 GEN
mulss(long x,long y)93 mulss(long x, long y)
94 {
95 long p1;
96 LOCAL_HIREMAINDER;
97
98 if (!x || !y) return gen_0;
99 if (x<0) {
100 x = -x;
101 if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
102 p1 = mulll(x,y); return uutoineg(hiremainder, p1);
103 } else {
104 if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
105 p1 = mulll(x,y); return uutoi(hiremainder, p1);
106 }
107 }
108 GEN
sqrs(long x)109 sqrs(long x)
110 {
111 long p1;
112 LOCAL_HIREMAINDER;
113
114 if (!x) return gen_0;
115 if (x<0) x = -x;
116 p1 = mulll(x,x); return uutoi(hiremainder, p1);
117 }
118 GEN
muluu(ulong x,ulong y)119 muluu(ulong x, ulong y)
120 {
121 long p1;
122 LOCAL_HIREMAINDER;
123
124 if (!x || !y) return gen_0;
125 p1 = mulll(x,y); return uutoi(hiremainder, p1);
126 }
127 GEN
sqru(ulong x)128 sqru(ulong x)
129 {
130 long p1;
131 LOCAL_HIREMAINDER;
132
133 if (!x) return gen_0;
134 p1 = mulll(x,x); return uutoi(hiremainder, p1);
135 }
136
137 /* assume x > 1, y != 0. Return u * y with sign s */
138 static GEN
mulur_2(ulong x,GEN y,long s)139 mulur_2(ulong x, GEN y, long s)
140 {
141 long m, sh, i, lx = lg(y), e = expo(y);
142 GEN z = cgetr(lx);
143 ulong garde;
144 LOCAL_HIREMAINDER;
145
146 y--; garde = mulll(x,y[lx]);
147 for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
148 z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
149 sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
150 if (sh) shift_left(z,z, 2,lx-1, garde,sh);
151 z[1] = evalsigne(s) | evalexpo(m+e);
152 if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
153 return z;
154 }
155
156 INLINE GEN
mul0r(GEN x)157 mul0r(GEN x)
158 {
159 long l = lg(x), e = expo(x);
160 e = (l > 2)? -prec2nbits(l) + e: (e < 0? 2*e: 0);
161 return real_0_bit(e);
162 }
163 /* lg(x) > 2 */
164 INLINE GEN
div0r(GEN x)165 div0r(GEN x) {
166 long l = lg(x), e = expo(x);
167 return real_0_bit(-prec2nbits(l) - e);
168 }
169
170 GEN
mulsr(long x,GEN y)171 mulsr(long x, GEN y)
172 {
173 long s;
174
175 if (!x) return mul0r(y);
176 s = signe(y);
177 if (!s)
178 {
179 if (x < 0) x = -x;
180 return real_0_bit( expo(y) + expu(x) );
181 }
182 if (x==1) return rcopy(y);
183 if (x==-1) return negr(y);
184 if (x < 0)
185 return mulur_2((ulong)-x, y, -s);
186 else
187 return mulur_2((ulong)x, y, s);
188 }
189
190 GEN
mulur(ulong x,GEN y)191 mulur(ulong x, GEN y)
192 {
193 long s;
194
195 if (!x) return mul0r(y);
196 s = signe(y);
197 if (!s) return real_0_bit( expo(y) + expu(x) );
198 if (x==1) return rcopy(y);
199 return mulur_2(x, y, s);
200 }
201
202 INLINE void
mulrrz_end(GEN z,GEN hi,long lz,long sz,long ez,ulong garde)203 mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
204 {
205 long i;
206 if (hi[2] < 0)
207 {
208 if (z != hi)
209 for (i=2; i<lz ; i++) z[i] = hi[i];
210 ez++;
211 }
212 else
213 {
214 shift_left(z,hi,2,lz-1, garde, 1);
215 garde <<= 1;
216 }
217 if (garde & HIGHBIT)
218 { /* round to nearest */
219 i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
220 if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
221 }
222 z[1] = evalsigne(sz)|evalexpo(ez);
223 }
224 /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
225 INLINE void
mulrrz_3end(GEN z,long sz,long ez,ulong garde)226 mulrrz_3end(GEN z, long sz, long ez, ulong garde)
227 {
228 if (z[2] < 0)
229 { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
230 if (garde & HIGHBIT) z[2]++; /* round properly */
231 ez++;
232 }
233 else
234 {
235 uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
236 if (garde & (1UL<<(BITS_IN_LONG-2)))
237 {
238 uel(z,2)++; /* round properly, z2+1 can overflow */
239 if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
240 }
241 }
242 z[1] = evalsigne(sz)|evalexpo(ez);
243 }
244
245 /* set z <-- x^2 != 0, floating point multiplication.
246 * lz = lg(z) = lg(x) */
247 INLINE void
sqrz_i(GEN z,GEN x,long lz)248 sqrz_i(GEN z, GEN x, long lz)
249 {
250 long ez = 2*expo(x);
251 long i, j, lzz, p1;
252 ulong garde;
253 GEN x1;
254 LOCAL_HIREMAINDER;
255 LOCAL_OVERFLOW;
256
257 if (lz > SQRR_SQRI_LIMIT)
258 {
259 pari_sp av = avma;
260 GEN hi = sqrispec_mirror(x+2, lz-2);
261 mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
262 set_avma(av); return;
263 }
264 if (lz == 3)
265 {
266 garde = mulll(x[2],x[2]);
267 z[2] = hiremainder;
268 mulrrz_3end(z, 1, ez, garde);
269 return;
270 }
271
272 lzz = lz-1; p1 = x[lzz];
273 if (p1)
274 {
275 (void)mulll(p1,x[3]);
276 garde = addmul(p1,x[2]);
277 z[lzz] = hiremainder;
278 }
279 else
280 {
281 garde = 0;
282 z[lzz] = 0;
283 }
284 for (j=lz-2, x1=x-j; j>=3; j--)
285 {
286 p1 = x[j]; x1++;
287 if (p1)
288 {
289 (void)mulll(p1,x1[lz+1]);
290 garde = addll(addmul(p1,x1[lz]), garde);
291 for (i=lzz; i>j; i--)
292 {
293 hiremainder += overflow;
294 z[i] = addll(addmul(p1,x1[i]), z[i]);
295 }
296 z[j] = hiremainder+overflow;
297 }
298 else z[j]=0;
299 }
300 p1 = x[2]; x1++;
301 garde = addll(mulll(p1,x1[lz]), garde);
302 for (i=lzz; i>2; i--)
303 {
304 hiremainder += overflow;
305 z[i] = addll(addmul(p1,x1[i]), z[i]);
306 }
307 z[2] = hiremainder+overflow;
308 mulrrz_end(z, z, lz, 1, ez, garde);
309 }
310
311 /* lz "large" = lg(y) = lg(z), lg(x) > lz if flag = 1 and >= if flag = 0 */
312 INLINE void
mulrrz_int(GEN z,GEN x,GEN y,long lz,long flag,long sz)313 mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
314 {
315 pari_sp av = avma;
316 GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
317 mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
318 set_avma(av);
319 }
320
321 /* lz = 3 */
322 INLINE void
mulrrz_3(GEN z,GEN x,GEN y,long flag,long sz)323 mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
324 {
325 ulong garde;
326 LOCAL_HIREMAINDER;
327 if (flag)
328 {
329 (void)mulll(x[2],y[3]);
330 garde = addmul(x[2],y[2]);
331 }
332 else
333 garde = mulll(x[2],y[2]);
334 z[2] = hiremainder;
335 mulrrz_3end(z, sz, expo(x)+expo(y), garde);
336 }
337
338 /* set z <-- x*y, floating point multiplication. Trailing 0s for x are
339 * treated efficiently (important application: mulir).
340 * lz = lg(z) = lg(x) <= ly <= lg(y), sz = signe(z). flag = lg(x) < lg(y) */
341 INLINE void
mulrrz_i(GEN z,GEN x,GEN y,long lz,long flag,long sz)342 mulrrz_i(GEN z, GEN x, GEN y, long lz, long flag, long sz)
343 {
344 long ez, i, j, lzz, p1;
345 ulong garde;
346 GEN y1;
347 LOCAL_HIREMAINDER;
348 LOCAL_OVERFLOW;
349
350 if (x == y) { sqrz_i(z,x,lz); return; }
351 if (lz > MULRR_MULII_LIMIT) { mulrrz_int(z,x,y,lz,flag,sz); return; }
352 if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
353 ez = expo(x) + expo(y);
354 if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
355 lzz=lz-1; p1=x[lzz];
356 if (p1)
357 {
358 (void)mulll(p1,y[3]);
359 garde = addll(addmul(p1,y[2]), garde);
360 z[lzz] = overflow+hiremainder;
361 }
362 else z[lzz]=0;
363 for (j=lz-2, y1=y-j; j>=3; j--)
364 {
365 p1 = x[j]; y1++;
366 if (p1)
367 {
368 (void)mulll(p1,y1[lz+1]);
369 garde = addll(addmul(p1,y1[lz]), garde);
370 for (i=lzz; i>j; i--)
371 {
372 hiremainder += overflow;
373 z[i] = addll(addmul(p1,y1[i]), z[i]);
374 }
375 z[j] = hiremainder+overflow;
376 }
377 else z[j]=0;
378 }
379 p1 = x[2]; y1++;
380 garde = addll(mulll(p1,y1[lz]), garde);
381 for (i=lzz; i>2; i--)
382 {
383 hiremainder += overflow;
384 z[i] = addll(addmul(p1,y1[i]), z[i]);
385 }
386 z[2] = hiremainder+overflow;
387 mulrrz_end(z, z, lz, sz, ez, garde);
388 }
389
390 GEN
mulrr(GEN x,GEN y)391 mulrr(GEN x, GEN y)
392 {
393 long flag, ly, lz, sx, sy;
394 GEN z;
395
396 if (x == y) return sqrr(x);
397 sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
398 sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
399 if (sy < 0) sx = -sx;
400 lz = lg(x);
401 ly = lg(y);
402 if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
403 z = cgetr(lz);
404 mulrrz_i(z, x,y, lz,flag, sx);
405 return z;
406 }
407
408 GEN
sqrr(GEN x)409 sqrr(GEN x)
410 {
411 long lz, sx = signe(x);
412 GEN z;
413
414 if (!sx) return real_0_bit(2*expo(x));
415 lz = lg(x); z = cgetr(lz);
416 sqrz_i(z, x, lz);
417 return z;
418 }
419
420 GEN
mulir(GEN x,GEN y)421 mulir(GEN x, GEN y)
422 {
423 long sx = signe(x), sy;
424 if (!sx) return mul0r(y);
425 if (lgefint(x) == 3) {
426 GEN z = mulur(uel(x,2), y);
427 if (sx < 0) togglesign(z);
428 return z;
429 }
430 sy = signe(y);
431 if (!sy) return real_0_bit(expi(x) + expo(y));
432 if (sy < 0) sx = -sx;
433 {
434 long lz = lg(y), lx = lgefint(x);
435 GEN hi, z = cgetr(lz);
436 pari_sp av = avma;
437 if (lx < (lz>>1) || (lx < lz && lz > MULRR_MULII_LIMIT))
438 { /* size mantissa of x < half size of mantissa z, or lx < lz so large
439 * that mulrr will call mulii anyway: mulii */
440 x = itor(x, lx);
441 hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
442 mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
443 }
444 else /* dubious: complete x with 0s and call mulrr */
445 mulrrz_i(z, itor(x,lz), y, lz, 0, sx);
446 set_avma(av); return z;
447 }
448 }
449
450 /* x + y*z, generic. If lgefint(z) <= 3, caller should use faster variants */
451 static GEN
addmulii_gen(GEN x,GEN y,GEN z,long lz)452 addmulii_gen(GEN x, GEN y, GEN z, long lz)
453 {
454 long lx = lgefint(x), ly;
455 pari_sp av;
456 GEN t;
457 if (lx == 2) return mulii(z,y);
458 ly = lgefint(y);
459 if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
460 av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
461 t = mulii(z, y);
462 set_avma(av); return addii(t,x);
463 }
464 /* x + y*z, lgefint(z) == 3 */
465 static GEN
addmulii_lg3(GEN x,GEN y,GEN z)466 addmulii_lg3(GEN x, GEN y, GEN z)
467 {
468 long s = signe(z), lx, ly;
469 ulong w = z[2];
470 pari_sp av;
471 GEN t;
472 if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
473 lx = lgefint(x);
474 ly = lgefint(y);
475 if (lx == 2)
476 { /* x = 0 */
477 if (ly == 2) return gen_0;
478 t = muluispec(w, y+2, ly-2);
479 if (signe(y) < 0) s = -s;
480 setsigne(t, s); return t;
481 }
482 if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
483 av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
484 t = muluispec(w, y+2, ly-2);
485 if (signe(y) < 0) s = -s;
486 setsigne(t, s);
487 set_avma(av); return addii(x,t);
488 }
489 /* x + y*z */
490 GEN
addmulii(GEN x,GEN y,GEN z)491 addmulii(GEN x, GEN y, GEN z)
492 {
493 long lz = lgefint(z);
494 switch(lz)
495 {
496 case 2: return icopy(x); /* z = 0, wasteful copy */
497 case 3: return addmulii_lg3(x, y, z);
498 default:return addmulii_gen(x, y, z, lz);
499 }
500 }
501 /* x + y*z, returns x itself and not a copy when y*z = 0 */
502 GEN
addmulii_inplace(GEN x,GEN y,GEN z)503 addmulii_inplace(GEN x, GEN y, GEN z)
504 {
505 long lz;
506 if (lgefint(y) == 2) return x;
507 lz = lgefint(z);
508 switch(lz)
509 {
510 case 2: return x;
511 case 3: return addmulii_lg3(x, y, z);
512 default:return addmulii_gen(x, y, z, lz);
513 }
514 }
515
516 /* written by Bruno Haible following an idea of Robert Harley */
517 long
vals(ulong z)518 vals(ulong z)
519 {
520 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};
521 #ifdef LONG_IS_64BIT
522 long s;
523 #endif
524
525 if (!z) return -1;
526 #ifdef LONG_IS_64BIT
527 if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
528 #endif
529 z |= ~z + 1;
530 z += z << 4;
531 z += z << 6;
532 z ^= z << 16; /* or z -= z<<16 */
533 #ifdef LONG_IS_64BIT
534 return s + tab[(z&0xffffffff)>>26];
535 #else
536 return tab[z>>26];
537 #endif
538 }
539
540 GEN
divsi(long x,GEN y)541 divsi(long x, GEN y)
542 {
543 long p1, s = signe(y);
544 LOCAL_HIREMAINDER;
545
546 if (!s) pari_err_INV("divsi",gen_0);
547 if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
548 hiremainder=0; p1=divll(labs(x),y[2]);
549 if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
550 if (s<0) p1 = -p1;
551 return stoi(p1);
552 }
553
554 GEN
divir(GEN x,GEN y)555 divir(GEN x, GEN y)
556 {
557 GEN z;
558 long ly = lg(y), lx = lgefint(x);
559 pari_sp av;
560
561 if (ly == 2) pari_err_INV("divir",y);
562 if (lx == 2) return div0r(y);
563 if (lx == 3) {
564 z = divur(x[2], y);
565 if (signe(x) < 0) togglesign(z);
566 return z;
567 }
568 z = cgetr(ly); av = avma;
569 affrr(divrr(itor(x, ly+1), y), z);
570 set_avma(av); return z;
571 }
572
573 GEN
divur(ulong x,GEN y)574 divur(ulong x, GEN y)
575 {
576 pari_sp av;
577 long ly = lg(y);
578 GEN z;
579
580 if (ly == 2) pari_err_INV("divur",y);
581 if (!x) return div0r(y);
582 if (ly > INVNEWTON_LIMIT) {
583 av = avma; z = invr(y);
584 if (x == 1) return z;
585 return gerepileuptoleaf(av, mulur(x, z));
586 }
587 z = cgetr(ly); av = avma;
588 affrr(divrr(utor(x,ly+1), y), z);
589 set_avma(av); return z;
590 }
591
592 GEN
divsr(long x,GEN y)593 divsr(long x, GEN y)
594 {
595 pari_sp av;
596 long ly = lg(y);
597 GEN z;
598
599 if (ly == 2) pari_err_INV("divsr",y);
600 if (!x) return div0r(y);
601 if (ly > INVNEWTON_LIMIT) {
602 av = avma; z = invr(y);
603 if (x == 1) return z;
604 if (x ==-1) { togglesign(z); return z; }
605 return gerepileuptoleaf(av, mulsr(x, z));
606 }
607 z = cgetr(ly); av = avma;
608 affrr(divrr(stor(x,ly+1), y), z);
609 set_avma(av); return z;
610 }
611
612 /* returns 1/y, assume y != 0 */
613 static GEN
invr_basecase(GEN y)614 invr_basecase(GEN y)
615 {
616 long ly = lg(y);
617 GEN z = cgetr(ly);
618 pari_sp av = avma;
619 affrr(divrr(real_1(ly+1), y), z);
620 set_avma(av); return z;
621 }
622 /* returns 1/b, Newton iteration */
623 GEN
invr(GEN b)624 invr(GEN b)
625 {
626 const long s = 6;
627 long i, p, l = lg(b);
628 GEN x, a;
629 ulong mask;
630
631 if (l <= maxss(INVNEWTON_LIMIT, (1L<<s) + 2)) {
632 if (l == 2) pari_err_INV("invr",b);
633 return invr_basecase(b);
634 }
635 mask = quadratic_prec_mask(l-2);
636 for(i=0, p=1; i<s; i++) { p <<= 1; if (mask & 1) p--; mask >>= 1; }
637 x = cgetr(l);
638 a = rcopy(b); a[1] = _evalexpo(0) | evalsigne(1);
639 affrr(invr_basecase(rtor(a, p+2)), x);
640 while (mask > 1)
641 {
642 p <<= 1; if (mask & 1) p--;
643 mask >>= 1;
644 setlg(a, p + 2);
645 setlg(x, p + 2);
646 /* TODO: mulrr(a,x) should be a half product (the higher half is known).
647 * mulrr(x, ) already is */
648 affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
649 set_avma((pari_sp)a);
650 }
651 x[1] = (b[1] & SIGNBITS) | evalexpo(expo(x)-expo(b));
652 set_avma((pari_sp)x); return x;
653 }
654
655 GEN
modii(GEN x,GEN y)656 modii(GEN x, GEN y)
657 {
658 switch(signe(x))
659 {
660 case 0: return gen_0;
661 case 1: return remii(x,y);
662 default:
663 {
664 pari_sp av = avma;
665 (void)new_chunk(lgefint(y));
666 x = remii(x,y); set_avma(av);
667 if (x==gen_0) return x;
668 return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
669 }
670 }
671 }
672
673 void
modiiz(GEN x,GEN y,GEN z)674 modiiz(GEN x, GEN y, GEN z)
675 {
676 const pari_sp av = avma;
677 affii(modii(x,y),z); set_avma(av);
678 }
679
680 GEN
divrs(GEN x,long y)681 divrs(GEN x, long y)
682 {
683 GEN z;
684 if (y < 0)
685 {
686 z = divru(x, (ulong)-y);
687 togglesign(z);
688 }
689 else
690 z = divru(x, (ulong)y);
691 return z;
692 }
693
694 GEN
divru(GEN x,ulong y)695 divru(GEN x, ulong y)
696 {
697 long i, lx, sh, e, s = signe(x);
698 ulong garde;
699 GEN z;
700 LOCAL_HIREMAINDER;
701
702 if (!y) pari_err_INV("divru",gen_0);
703 if (!s) return real_0_bit(expo(x) - expu(y));
704 if (!(y & (y-1))) /* power of 2 */
705 {
706 if (y == 1) return rcopy(x);
707 return shiftr(x, -expu(y));
708 }
709 e = expo(x);
710 lx = lg(x);
711 z = cgetr(lx);
712 if (lx == 3)
713 {
714 if (y <= uel(x,2))
715 {
716 hiremainder = 0;
717 z[2] = divll(x[2],y);
718 /* we may have hiremainder != 0 ==> garde */
719 garde = divll(0,y);
720 }
721 else
722 {
723 hiremainder = x[2];
724 z[2] = divll(0,y);
725 garde = hiremainder;
726 e -= BITS_IN_LONG;
727 }
728 }
729 else
730 {
731 ulong yp = get_Fl_red(y);
732 if (y <= uel(x,2))
733 {
734 hiremainder = 0;
735 for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
736 /* we may have hiremainder != 0 ==> garde */
737 garde = divll_pre(0,y,yp);
738 }
739 else
740 {
741 long l = lx-1;
742 hiremainder = x[2];
743 for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
744 z[i] = divll_pre(0,y,yp);
745 garde = hiremainder;
746 e -= BITS_IN_LONG;
747 }
748 }
749 sh=bfffo(z[2]); /* z[2] != 0 */
750 if (sh) shift_left(z,z, 2,lx-1, garde,sh);
751 z[1] = evalsigne(s) | evalexpo(e-sh);
752 if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
753 return z;
754 }
755
756 GEN
truedvmdii(GEN x,GEN y,GEN * z)757 truedvmdii(GEN x, GEN y, GEN *z)
758 {
759 pari_sp av;
760 GEN r, q, *gptr[2];
761 if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
762 if (z == ONLY_REM) return modii(x,y);
763
764 av = avma;
765 q = dvmdii(x,y,&r); /* assume that r is last on stack */
766 switch(signe(r))
767 {
768 case 0:
769 if (z) *z = gen_0;
770 return q;
771 case 1:
772 if (z) *z = r; else cgiv(r);
773 return q;
774 case -1: break;
775 }
776 q = addis(q, -signe(y));
777 if (!z) return gerepileuptoint(av, q);
778
779 *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
780 gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(pari_sp)r,gptr,2);
781 return q;
782 }
783 GEN
truedvmdis(GEN x,long y,GEN * z)784 truedvmdis(GEN x, long y, GEN *z)
785 {
786 pari_sp av = avma;
787 long r;
788 GEN q;
789
790 if (z == ONLY_REM) return modis(x, y);
791 q = divis_rem(x,y,&r);
792
793 if (r >= 0)
794 {
795 if (z) *z = utoi(r);
796 return q;
797 }
798 q = gerepileuptoint(av, addis(q, (y < 0)? 1: -1));
799 if (z) *z = utoi(r + labs(y));
800 return q;
801 }
802 GEN
truedvmdsi(long x,GEN y,GEN * z)803 truedvmdsi(long x, GEN y, GEN *z)
804 {
805 long q, r;
806 if (z == ONLY_REM) return modsi(x, y);
807 q = sdivsi_rem(x,y,&r);
808 if (r >= 0) {
809 if (z) *z = utoi(r);
810 return stoi(q);
811 }
812 q = q - signe(y);
813 if (!z) return stoi(q);
814
815 *z = subiuspec(y+2,(ulong)-r, lgefint(y)-2);
816 return stoi(q);
817 }
818
819 /* 2^n = shifti(gen_1, n) */
820 GEN
int2n(long n)821 int2n(long n) {
822 long i, m, l;
823 GEN z;
824 if (n < 0) return gen_0;
825 if (n == 0) return gen_1;
826
827 l = dvmdsBIL(n, &m) + 3;
828 z = cgetipos(l);
829 for (i = 2; i < l; i++) z[i] = 0;
830 *int_MSW(z) = 1UL << m; return z;
831 }
832 /* To avoid problems when 2^(BIL-1) < n. Overflow cleanly, where int2n
833 * returns gen_0 */
834 GEN
int2u(ulong n)835 int2u(ulong n) {
836 ulong i, m, l;
837 GEN z;
838 if (n == 0) return gen_1;
839
840 l = dvmduBIL(n, &m) + 3;
841 z = cgetipos(l);
842 for (i = 2; i < l; i++) z[i] = 0;
843 *int_MSW(z) = 1UL << m; return z;
844 }
845 /* 2^n - 1 */
846 GEN
int2um1(ulong n)847 int2um1(ulong n) {
848 ulong i, m, l;
849 GEN z;
850 if (n == 0) return gen_0;
851
852 l = dvmduBIL(n, &m);
853 l += m? 3: 2;
854 z = cgetipos(l);
855 for (i = 2; i < l; i++) z[i] = ~0UL;
856 if (m) *int_MSW(z) = (1UL << m) - 1;
857 return z;
858 }
859
860 GEN
shifti(GEN x,long n)861 shifti(GEN x, long n)
862 {
863 long s = signe(x);
864 GEN y;
865
866 if(s == 0) return gen_0;
867 y = shiftispec(x + 2, lgefint(x) - 2, n);
868 if (signe(y)) setsigne(y, s);
869 return y;
870 }
871
872 /* actual operations will take place on a+2 and b+2: we strip the codewords */
873 GEN
mulii(GEN a,GEN b)874 mulii(GEN a,GEN b)
875 {
876 long sa,sb;
877 GEN z;
878
879 sa=signe(a); if (!sa) return gen_0;
880 sb=signe(b); if (!sb) return gen_0;
881 if (sb<0) sa = -sa;
882 z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
883 setsigne(z,sa); return z;
884 }
885
886 GEN
sqri(GEN a)887 sqri(GEN a) { return sqrispec(a+2, lgefint(a)-2); }
888
889 /* sqrt()'s result may be off by 1 when a is not representable exactly as a
890 * double [64bit machine] */
891 ulong
usqrt(ulong a)892 usqrt(ulong a)
893 {
894 ulong x = (ulong)sqrt((double)a);
895 #ifdef LONG_IS_64BIT
896 if (x > LOWMASK || x*x > a) x--;
897 #endif
898 return x;
899 }
900
901 /********************************************************************/
902 /** **/
903 /** EXPONENT / CONVERSION t_REAL --> double **/
904 /** **/
905 /********************************************************************/
906
907 #ifdef LONG_IS_64BIT
908 long
dblexpo(double x)909 dblexpo(double x)
910 {
911 union { double f; ulong i; } fi;
912 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
913 const int exp_mid = 0x3ff;/* exponent bias */
914
915 if (x==0.) return -exp_mid;
916 fi.f = x;
917 return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
918 }
919
920 ulong
dblmantissa(double x)921 dblmantissa(double x)
922 {
923 union { double f; ulong i; } fi;
924 const int expo_len = 11; /* number of bits of exponent */
925
926 if (x==0.) return 0;
927 fi.f = x;
928 return (fi.i << expo_len) | HIGHBIT;
929 }
930
931 GEN
dbltor(double x)932 dbltor(double x)
933 {
934 GEN z;
935 long e;
936 union { double f; ulong i; } fi;
937 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
938 const int exp_mid = 0x3ff;/* exponent bias */
939 const int expo_len = 11; /* number of bits of exponent */
940
941 if (x==0.) return real_0_bit(-exp_mid);
942 fi.f = x; z = cgetr(DEFAULTPREC);
943 {
944 const ulong a = fi.i;
945 ulong A;
946 e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
947 if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
948 A = a << expo_len;
949 if (e == -exp_mid)
950 { /* unnormalized values */
951 int sh = bfffo(A);
952 e -= sh-1;
953 z[2] = A << sh;
954 }
955 else
956 z[2] = HIGHBIT | A;
957 z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
958 }
959 return z;
960 }
961
962 double
rtodbl(GEN x)963 rtodbl(GEN x)
964 {
965 long ex,s=signe(x);
966 ulong a;
967 union { double f; ulong i; } fi;
968 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
969 const int exp_mid = 0x3ff;/* exponent bias */
970 const int expo_len = 11; /* number of bits of exponent */
971
972 if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
973
974 /* start by rounding to closest */
975 a = (x[2] & (HIGHBIT-1)) + 0x400;
976 if (a & HIGHBIT) { ex++; a=0; }
977 if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
978 fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
979 if (s<0) fi.i |= HIGHBIT;
980 return fi.f;
981 }
982
983 #else /* LONG_IS_64BIT */
984
985 #if PARI_DOUBLE_FORMAT == 1
986 # define INDEX0 1
987 # define INDEX1 0
988 #elif PARI_DOUBLE_FORMAT == 0
989 # define INDEX0 0
990 # define INDEX1 1
991 #endif
992
993 long
dblexpo(double x)994 dblexpo(double x)
995 {
996 union { double f; ulong i[2]; } fi;
997 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
998 const int exp_mid = 0x3ff;/* exponent bias */
999 const int shift = mant_len-32;
1000
1001 if (x==0.) return -exp_mid;
1002 fi.f = x;
1003 {
1004 const ulong a = fi.i[INDEX0];
1005 return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
1006 }
1007 }
1008
1009 ulong
dblmantissa(double x)1010 dblmantissa(double x)
1011 {
1012 union { double f; ulong i[2]; } fi;
1013 const int expo_len = 11; /* number of bits of exponent */
1014
1015 if (x==0.) return 0;
1016 fi.f = x;
1017 {
1018 const ulong a = fi.i[INDEX0];
1019 const ulong b = fi.i[INDEX1];
1020 return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1021 }
1022 }
1023
1024 GEN
dbltor(double x)1025 dbltor(double x)
1026 {
1027 GEN z;
1028 long e;
1029 union { double f; ulong i[2]; } fi;
1030 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1031 const int exp_mid = 0x3ff;/* exponent bias */
1032 const int expo_len = 11; /* number of bits of exponent */
1033 const int shift = mant_len-32;
1034
1035 if (x==0.) return real_0_bit(-exp_mid);
1036 fi.f = x; z = cgetr(DEFAULTPREC);
1037 {
1038 const ulong a = fi.i[INDEX0];
1039 const ulong b = fi.i[INDEX1];
1040 ulong A, B;
1041 e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
1042 if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
1043 A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1044 B = b << expo_len;
1045 if (e == -exp_mid)
1046 { /* unnormalized values */
1047 int sh;
1048 if (A)
1049 {
1050 sh = bfffo(A);
1051 e -= sh-1;
1052 z[2] = (A << sh) | (B >> (32-sh));
1053 z[3] = B << sh;
1054 }
1055 else
1056 {
1057 sh = bfffo(B); /* B != 0 */
1058 e -= sh-1 + 32;
1059 z[2] = B << sh;
1060 z[3] = 0;
1061 }
1062 }
1063 else
1064 {
1065 z[3] = B;
1066 z[2] = HIGHBIT | A;
1067 }
1068 z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
1069 }
1070 return z;
1071 }
1072
1073 double
rtodbl(GEN x)1074 rtodbl(GEN x)
1075 {
1076 long ex,s=signe(x),lx=lg(x);
1077 ulong a,b,k;
1078 union { double f; ulong i[2]; } fi;
1079 const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1080 const int exp_mid = 0x3ff;/* exponent bias */
1081 const int expo_len = 11; /* number of bits of exponent */
1082 const int shift = mant_len-32;
1083
1084 if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
1085
1086 /* start by rounding to closest */
1087 a = x[2] & (HIGHBIT-1);
1088 if (lx > 3)
1089 {
1090 b = x[3] + 0x400UL; if (b < 0x400UL) a++;
1091 if (a & HIGHBIT) { ex++; a=0; }
1092 }
1093 else b = 0;
1094 if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
1095 ex += exp_mid;
1096 k = (a >> expo_len) | (ex << shift);
1097 if (s<0) k |= HIGHBIT;
1098 fi.i[INDEX0] = k;
1099 fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
1100 return fi.f;
1101 }
1102 #endif /* LONG_IS_64BIT */
1103
1104