Line data Source code
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
18 60521417 : invmod2BIL(ulong b)
19 : {
20 : static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
21 60521417 : 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 60521417 : x = x*(2-b*x); /* one more pass necessary */
26 : #endif
27 60521417 : x = x*(2-b*x);
28 60521417 : x = x*(2-b*x); return x*(2-b*x);
29 : }
30 :
31 : void
32 1356701086 : affrr(GEN x, GEN y)
33 : {
34 1356701086 : long i, lx, ly = lg(y);
35 1356701086 : if (!signe(x))
36 : {
37 126243 : y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
38 126243 : return;
39 : }
40 1356574843 : y[1] = x[1]; lx = lg(x);
41 1356574843 : if (lx <= ly)
42 : {
43 7299679445 : for (i=2; i<lx; i++) y[i]=x[i];
44 1300240471 : for ( ; i<ly; i++) y[i]=0;
45 1143640100 : return;
46 : }
47 1140168995 : for (i=2; i<ly; i++) y[i]=x[i];
48 : /* lx > ly: round properly */
49 212934743 : if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
50 : }
51 :
52 : GEN
53 62523688 : trunc2nr(GEN x, long n)
54 : {
55 : long ex;
56 62523688 : if (!signe(x)) return gen_0;
57 62085180 : ex = expo(x) + n; if (ex < 0) return gen_0;
58 59533612 : 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
63 52609734 : mantissa_real(GEN x, long *e)
64 : {
65 52609734 : *e = bit_prec(x)-1-expo(x);
66 52609523 : return mantissa2nr(x, 0);
67 : }
68 :
69 : GEN
70 903386888 : mului(ulong x, GEN y)
71 : {
72 903386888 : long s = signe(y);
73 : GEN z;
74 :
75 903386888 : if (!s || !x) return gen_0;
76 748914792 : z = muluispec(x, y+2, lgefint(y)-2);
77 748550530 : setsigne(z,s); return z;
78 : }
79 :
80 : GEN
81 603796325 : mulsi(long x, GEN y)
82 : {
83 603796325 : long s = signe(y);
84 : GEN z;
85 :
86 603796325 : if (!s || !x) return gen_0;
87 360704324 : if (x<0) { s = -s; x = -x; }
88 360704324 : z = muluispec((ulong)x, y+2, lgefint(y)-2);
89 360699881 : setsigne(z,s); return z;
90 : }
91 :
92 : GEN
93 201096682 : mulss(long x, long y)
94 : {
95 : ulong p1;
96 : LOCAL_HIREMAINDER;
97 :
98 201096682 : if (!x || !y) return gen_0;
99 200761242 : if (x<0) {
100 413524 : x = -x;
101 413524 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
102 299122 : p1 = mulll(x,y); return uutoineg(hiremainder, p1);
103 : } else {
104 200347718 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
105 200341384 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
106 : }
107 : }
108 : GEN
109 1498020 : sqrs(long x)
110 : {
111 : ulong p1;
112 : LOCAL_HIREMAINDER;
113 :
114 1498020 : if (!x) return gen_0;
115 1496249 : if (x<0) x = -x;
116 1496249 : p1 = mulll(x,x); return uutoi(hiremainder, p1);
117 : }
118 : GEN
119 4424259578 : muluu(ulong x, ulong y)
120 : {
121 : ulong p1;
122 : LOCAL_HIREMAINDER;
123 :
124 4424259578 : if (!x || !y) return gen_0;
125 4424949077 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
126 : }
127 : GEN
128 606418016 : sqru(ulong x)
129 : {
130 : ulong p1;
131 : LOCAL_HIREMAINDER;
132 :
133 606418016 : if (!x) return gen_0;
134 606084966 : 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
139 701492903 : mulur_2(ulong x, GEN y, long s)
140 : {
141 701492903 : long m, sh, i, lx = lg(y), e = expo(y);
142 701492903 : GEN z = cgetg(lx, t_REAL);
143 : ulong garde;
144 : LOCAL_HIREMAINDER;
145 :
146 700859644 : y--; garde = mulll(x,y[lx]);
147 3674234974 : for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
148 700859644 : z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
149 700859644 : sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
150 700859644 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
151 702142478 : z[1] = evalsigne(s) | evalexpo(m+e);
152 703600139 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
153 703353019 : return z;
154 : }
155 :
156 : INLINE GEN
157 473487 : mul0r(GEN x)
158 : {
159 473487 : long l = realprec(x), e = expo(x);
160 473487 : e = (l > 0)? e - l: (e < 0? 2*e: 0);
161 473487 : return real_0_bit(e);
162 : }
163 : /* lg(x) > 2 */
164 : INLINE GEN
165 186343 : div0r(GEN x) {
166 186343 : long l = realprec(x), e = expo(x);
167 186343 : return real_0_bit(-l - e);
168 : }
169 :
170 : GEN
171 116563547 : mulsr(long x, GEN y)
172 : {
173 : long s;
174 :
175 116563547 : if (!x) return mul0r(y);
176 116563505 : s = signe(y);
177 116563505 : if (!s)
178 : {
179 228716 : if (x < 0) x = -x;
180 228716 : return real_0_bit( expo(y) + expu(x) );
181 : }
182 116334789 : if (x==1) return rcopy(y);
183 107640607 : if (x==-1) return negr(y);
184 104525009 : if (x < 0)
185 34715465 : return mulur_2((ulong)-x, y, -s);
186 : else
187 69809544 : return mulur_2((ulong)x, y, s);
188 : }
189 :
190 : GEN
191 918825042 : mulur(ulong x, GEN y)
192 : {
193 : long s;
194 :
195 918825042 : if (!x) return mul0r(y);
196 918825035 : s = signe(y);
197 918825035 : if (!s) return real_0_bit( expo(y) + expu(x) );
198 911544267 : if (x==1) return rcopy(y);
199 596610226 : return mulur_2(x, y, s);
200 : }
201 :
202 : INLINE void
203 3551155706 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
204 : {
205 : long i;
206 3551155706 : if (hi[2] < 0)
207 : {
208 1644561354 : if (z != hi)
209 1711502039 : for (i=2; i<lz ; i++) z[i] = hi[i];
210 1644561354 : ez++;
211 : }
212 : else
213 : {
214 1906594352 : shift_left(z,hi,2,lz-1, garde, 1);
215 1909846684 : garde <<= 1;
216 : }
217 3554408038 : if (garde & HIGHBIT)
218 : { /* round to nearest */
219 1741063131 : i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
220 1723546082 : if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
221 : }
222 3554408038 : z[1] = evalsigne(sz)|evalexpo(ez);
223 3549735008 : }
224 : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
225 : INLINE void
226 769463112 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
227 : {
228 769463112 : if (z[2] < 0)
229 : { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
230 361108812 : if (garde & HIGHBIT) z[2]++; /* round properly */
231 361108812 : ez++;
232 : }
233 : else
234 : {
235 408354300 : uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
236 408354300 : if (garde & (1UL<<(BITS_IN_LONG-2)))
237 : {
238 167615246 : uel(z,2)++; /* round properly, z2+1 can overflow */
239 167615246 : if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
240 : }
241 : }
242 769463112 : z[1] = evalsigne(sz)|evalexpo(ez);
243 769423123 : }
244 :
245 : /* set z <-- x^2 != 0, floating point multiplication.
246 : * lz = lg(z) = lg(x) */
247 : INLINE void
248 559023094 : sqrz_i(GEN z, GEN x, long lz)
249 : {
250 559023094 : 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 559023094 : if (lz > prec2lg(SQRR_SQRI_LIMIT))
258 : {
259 41200730 : pari_sp av = avma;
260 41200730 : GEN hi = sqrispec_mirror(x+2, lz-2);
261 41512593 : mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
262 41494872 : set_avma(av); return;
263 : }
264 517722175 : if (lz == 3)
265 : {
266 120842059 : garde = mulll(x[2],x[2]);
267 120842059 : z[2] = hiremainder;
268 120842059 : mulrrz_3end(z, 1, ez, garde);
269 120821953 : return;
270 : }
271 :
272 396880116 : lzz = lz-1; p1 = x[lzz];
273 396880116 : if (p1)
274 : {
275 360939803 : (void)mulll(p1,x[3]);
276 360939803 : garde = addmul(p1,x[2]);
277 360939803 : z[lzz] = hiremainder;
278 : }
279 : else
280 : {
281 35940313 : garde = 0;
282 35940313 : z[lzz] = 0;
283 : }
284 1837270072 : for (j=lz-2, x1=x-j; j>=3; j--)
285 : {
286 1440389956 : p1 = x[j]; x1++;
287 1440389956 : if (p1)
288 : {
289 1413379990 : (void)mulll(p1,x1[lz+1]);
290 1413379990 : garde = addll(addmul(p1,x1[lz]), garde);
291 7222074025 : for (i=lzz; i>j; i--)
292 : {
293 5808694035 : hiremainder += overflow;
294 5808694035 : z[i] = addll(addmul(p1,x1[i]), z[i]);
295 : }
296 1413379990 : z[j] = hiremainder+overflow;
297 : }
298 27009966 : else z[j]=0;
299 : }
300 396880116 : p1 = x[2]; x1++;
301 396880116 : garde = addll(mulll(p1,x1[lz]), garde);
302 2239454770 : for (i=lzz; i>2; i--)
303 : {
304 1842574654 : hiremainder += overflow;
305 1842574654 : z[i] = addll(addmul(p1,x1[i]), z[i]);
306 : }
307 396880116 : z[2] = hiremainder+overflow;
308 396880116 : 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
313 54231420 : mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
314 : {
315 54231420 : pari_sp av = avma;
316 54231420 : GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
317 54231420 : mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
318 54231420 : set_avma(av);
319 54231420 : }
320 :
321 : /* lz = 3 */
322 : INLINE void
323 649606530 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
324 : {
325 : ulong garde;
326 : LOCAL_HIREMAINDER;
327 649606530 : if (flag)
328 : {
329 77733385 : (void)mulll(x[2],y[3]);
330 77733385 : garde = addmul(x[2],y[2]);
331 : }
332 : else
333 571873145 : garde = mulll(x[2],y[2]);
334 649606530 : z[2] = hiremainder;
335 649606530 : mulrrz_3end(z, sz, expo(x)+expo(y), garde);
336 649692656 : }
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
342 3697878464 : 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 3697878464 : if (x == y) { sqrz_i(z,x,lz); return; }
351 3697878464 : if (lz > prec2lg(MULRR_MULII_LIMIT)) { mulrrz_int(z,x,y,lz,flag,sz); return; }
352 3646915714 : if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
353 2997308233 : ez = expo(x) + expo(y);
354 2997308233 : if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
355 2997308233 : lzz=lz-1; p1=x[lzz];
356 2997308233 : if (p1)
357 : {
358 2796136897 : (void)mulll(p1,y[3]);
359 2796136897 : garde = addll(addmul(p1,y[2]), garde);
360 2796136897 : z[lzz] = overflow+hiremainder;
361 : }
362 201171336 : else z[lzz]=0;
363 13457646650 : for (j=lz-2, y1=y-j; j>=3; j--)
364 : {
365 10460338417 : p1 = x[j]; y1++;
366 10460338417 : if (p1)
367 : {
368 10161921247 : (void)mulll(p1,y1[lz+1]);
369 10161921247 : garde = addll(addmul(p1,y1[lz]), garde);
370 61430070781 : for (i=lzz; i>j; i--)
371 : {
372 51268149534 : hiremainder += overflow;
373 51268149534 : z[i] = addll(addmul(p1,y1[i]), z[i]);
374 : }
375 10161921247 : z[j] = hiremainder+overflow;
376 : }
377 298417170 : else z[j]=0;
378 : }
379 2997308233 : p1 = x[2]; y1++;
380 2997308233 : garde = addll(mulll(p1,y1[lz]), garde);
381 16516973013 : for (i=lzz; i>2; i--)
382 : {
383 13519664780 : hiremainder += overflow;
384 13519664780 : z[i] = addll(addmul(p1,y1[i]), z[i]);
385 : }
386 2997308233 : z[2] = hiremainder+overflow;
387 2997308233 : mulrrz_end(z, z, lz, sz, ez, garde);
388 : }
389 :
390 : GEN
391 3939736086 : mulrr(GEN x, GEN y)
392 : {
393 : long flag, ly, lz, sx, sy;
394 : GEN z;
395 :
396 3939736086 : if (x == y) return sqrr(x);
397 3938125656 : sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
398 3674423626 : sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
399 3641050875 : if (sy < 0) sx = -sx;
400 3641050875 : lz = lg(x);
401 3641050875 : ly = lg(y);
402 3641050875 : if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
403 3641050875 : z = cgetg(lz, t_REAL);
404 3636358369 : mulrrz_i(z, x,y, lz,flag, sx);
405 3644365549 : return z;
406 : }
407 :
408 : GEN
409 591232412 : sqrr(GEN x)
410 : {
411 591232412 : long lz, sx = signe(x);
412 : GEN z;
413 :
414 591232412 : if (!sx) return real_0_bit(2*expo(x));
415 559643238 : lz = lg(x); z = cgetg(lz, t_REAL);
416 559047833 : sqrz_i(z, x, lz);
417 559890517 : return z;
418 : }
419 :
420 : GEN
421 977518704 : mulir(GEN x, GEN y)
422 : {
423 977518704 : long sx = signe(x), sy;
424 977518704 : if (!sx) return mul0r(y);
425 977045266 : if (lgefint(x) == 3) {
426 866644086 : GEN z = mulur(uel(x,2), y);
427 868049041 : if (sx < 0) togglesign(z);
428 867764392 : return z;
429 : }
430 110401180 : sy = signe(y);
431 110401180 : if (!sy) return real_0_bit(expi(x) + expo(y));
432 109716570 : if (sy < 0) sx = -sx;
433 : {
434 109716570 : long lz = lg(y), lx = lgefint(x);
435 109716570 : GEN hi, z = cgetg(lz, t_REAL);
436 109245127 : pari_sp av = avma;
437 109245127 : if (lx < (lz>>1) || (lx < lz && lz > prec2lg(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 47770457 : x = itor(x, lg2prec(lx));
441 47474714 : hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
442 49056132 : mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
443 : }
444 : else /* dubious: complete x with 0s and call mulrr */
445 61447908 : mulrrz_i(z, itor(x, lg2prec(lz)), y, lz, 0, sx);
446 110901424 : 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
452 90358130 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
453 : {
454 90358130 : long lx = lgefint(x), ly;
455 : pari_sp av;
456 : GEN t;
457 90358130 : if (lx == 2) return mulii(z,y);
458 87582312 : ly = lgefint(y);
459 87582312 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
460 86978534 : av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
461 86979942 : t = mulii(z, y);
462 86972300 : set_avma(av); return addii(t,x);
463 : }
464 : /* x + y*z, lgefint(z) == 3 */
465 : static GEN
466 466334332 : addmulii_lg3(GEN x, GEN y, GEN z)
467 : {
468 466334332 : long s = signe(z), lx, ly;
469 466334332 : ulong w = z[2];
470 : pari_sp av;
471 : GEN t;
472 466334332 : if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
473 384866229 : lx = lgefint(x);
474 384866229 : ly = lgefint(y);
475 384866229 : if (lx == 2)
476 : { /* x = 0 */
477 71320597 : if (ly == 2) return gen_0;
478 39517595 : t = muluispec(w, y+2, ly-2);
479 39518228 : if (signe(y) < 0) s = -s;
480 39518228 : setsigne(t, s); return t;
481 : }
482 313545632 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
483 272914113 : av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
484 273057193 : t = muluispec(w, y+2, ly-2);
485 273008294 : if (signe(y) < 0) s = -s;
486 273008294 : setsigne(t, s);
487 273008294 : set_avma(av); return addii(x,t);
488 : }
489 : /* x + y*z */
490 : GEN
491 167023478 : addmulii(GEN x, GEN y, GEN z)
492 : {
493 167023478 : long lz = lgefint(z);
494 167023478 : switch(lz)
495 : {
496 740447 : case 2: return icopy(x); /* z = 0, wasteful copy */
497 126158826 : case 3: return addmulii_lg3(x, y, z);
498 40124205 : 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
503 1566125672 : addmulii_inplace(GEN x, GEN y, GEN z)
504 : {
505 : long lz;
506 1566125672 : if (lgefint(y) == 2) return x;
507 390220298 : lz = lgefint(z);
508 390220298 : switch(lz)
509 : {
510 17508 : case 2: return x;
511 340207825 : case 3: return addmulii_lg3(x, y, z);
512 49994965 : default:return addmulii_gen(x, y, z, lz);
513 : }
514 : }
515 :
516 : /* written by Bruno Haible following an idea of Robert Harley */
517 : long
518 2202921134 : 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 2202921134 : if (!z) return -1;
526 : #ifdef LONG_IS_64BIT
527 1976955551 : if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
528 : #endif
529 2202921134 : z |= ~z + 1;
530 2202921134 : z += z << 4;
531 2202921134 : z += z << 6;
532 2202921134 : z ^= z << 16; /* or z -= z<<16 */
533 : #ifdef LONG_IS_64BIT
534 1976955551 : return s + tab[(z&0xffffffff)>>26];
535 : #else
536 225965583 : return tab[z>>26];
537 : #endif
538 : }
539 :
540 : GEN
541 161 : divsi(long x, GEN y)
542 : {
543 161 : long p1, s = signe(y);
544 : LOCAL_HIREMAINDER;
545 :
546 161 : if (!s) pari_err_INV("divsi",gen_0);
547 161 : if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
548 161 : hiremainder=0; p1=divll(labs(x),y[2]);
549 161 : if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
550 161 : if (s<0) p1 = -p1;
551 161 : return stoi(p1);
552 : }
553 :
554 : GEN
555 3494459 : divir(GEN x, GEN y)
556 : {
557 : GEN z;
558 3494459 : long ly = lg(y), lx = lgefint(x);
559 : pari_sp av;
560 :
561 3494459 : if (ly == 2) pari_err_INV("divir",y);
562 3494461 : if (lx == 2) return div0r(y);
563 3308118 : if (lx == 3) {
564 1961309 : z = divur(x[2], y);
565 1961305 : if (signe(x) < 0) togglesign(z);
566 1961305 : return z;
567 : }
568 1346809 : z = cgetg(ly, t_REAL); av = avma;
569 1346809 : affrr(divrr(itor(x, lg2prec(ly+1)), y), z);
570 1346809 : set_avma(av); return z;
571 : }
572 :
573 : GEN
574 2493561 : divur(ulong x, GEN y)
575 : {
576 : pari_sp av;
577 2493561 : long p = realprec(y);
578 : GEN z;
579 :
580 2493561 : if (p == 0) pari_err_INV("divur",y);
581 2493561 : if (!x) return div0r(y);
582 2493561 : if (p > INVNEWTON_LIMIT) {
583 183 : av = avma; z = invr(y);
584 183 : if (x == 1) return z;
585 183 : return gerepileuptoleaf(av, mulur(x, z));
586 : }
587 2493378 : z = cgetr(p); av = avma;
588 2493339 : affrr(divrr(utor(x, p + BITS_IN_LONG), y), z);
589 2493372 : set_avma(av); return z;
590 : }
591 :
592 : GEN
593 798 : divsr(long x, GEN y)
594 : {
595 : pari_sp av;
596 798 : long p = realprec(y);
597 : GEN z;
598 :
599 798 : if (p == 0) pari_err_INV("divsr",y);
600 798 : if (!x) return div0r(y);
601 798 : if (p > INVNEWTON_LIMIT) {
602 0 : av = avma; z = invr(y);
603 0 : if (x == 1) return z;
604 0 : if (x ==-1) { togglesign(z); return z; }
605 0 : return gerepileuptoleaf(av, mulsr(x, z));
606 : }
607 798 : z = cgetr(p); av = avma;
608 798 : affrr(divrr(stor(x,p + BITS_IN_LONG), y), z);
609 798 : set_avma(av); return z;
610 : }
611 :
612 : /* returns 1/y, assume y != 0 */
613 : static GEN
614 62383209 : invr_basecase(GEN y)
615 : {
616 62383209 : long p = realprec(y);
617 62383209 : GEN z = cgetr(p);
618 62372803 : pari_sp av = avma;
619 62372803 : affrr(divrr(real_1(p + BITS_IN_LONG), y), z);
620 62382035 : set_avma(av); return z;
621 : }
622 : /* returns 1/b, Newton iteration */
623 : GEN
624 62384167 : invr(GEN b)
625 : {
626 62384167 : const long s = 6;
627 62384167 : long i, p, l = lg(b);
628 : GEN x, a;
629 : ulong mask;
630 :
631 62384167 : if (l <= maxss(prec2lg(INVNEWTON_LIMIT), (1L<<s) + 2)) {
632 62379733 : if (l == 2) pari_err_INV("invr",b);
633 62379733 : return invr_basecase(b);
634 : }
635 3675 : mask = quadratic_prec_mask(l-2);
636 25725 : for(i=0, p=1; i<s; i++) { p <<= 1; if (mask & 1) p--; mask >>= 1; }
637 3675 : x = cgetg(l, t_REAL);
638 3675 : a = rcopy(b); a[1] = _evalexpo(0) | evalsigne(1);
639 3675 : affrr(invr_basecase(rtor(a, lg2prec(p+2))), x);
640 12265 : while (mask > 1)
641 : {
642 8590 : p <<= 1; if (mask & 1) p--;
643 8590 : mask >>= 1;
644 8590 : setlg(a, p + 2);
645 8590 : 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 8590 : affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
649 8590 : set_avma((pari_sp)a);
650 : }
651 3675 : x[1] = (b[1] & SIGNBITS) | evalexpo(expo(x)-expo(b));
652 3675 : set_avma((pari_sp)x); return x;
653 : }
654 :
655 : GEN
656 2570275016 : modii(GEN x, GEN y)
657 : {
658 2570275016 : switch(signe(x))
659 : {
660 485441987 : case 0: return gen_0;
661 1549448334 : case 1: return remii(x,y);
662 535384695 : default:
663 : {
664 535384695 : pari_sp av = avma;
665 535384695 : (void)new_chunk(lgefint(y));
666 539769020 : x = remii(x,y); set_avma(av);
667 539855785 : if (x==gen_0) return x;
668 498479521 : return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
669 : }
670 : }
671 : }
672 :
673 : void
674 781873 : modiiz(GEN x, GEN y, GEN z)
675 : {
676 781873 : const pari_sp av = avma;
677 781873 : affii(modii(x,y),z); set_avma(av);
678 781841 : }
679 :
680 : GEN
681 9806286 : divrs(GEN x, long y)
682 : {
683 : GEN z;
684 9806286 : if (y < 0)
685 : {
686 5754975 : z = divru(x, (ulong)-y);
687 5754948 : togglesign(z);
688 : }
689 : else
690 4051311 : z = divru(x, (ulong)y);
691 9806246 : return z;
692 : }
693 :
694 : GEN
695 1095109773 : divru(GEN x, ulong y)
696 : {
697 1095109773 : long i, lx, sh, e, s = signe(x);
698 : ulong garde;
699 : GEN z;
700 : LOCAL_HIREMAINDER;
701 :
702 1095109773 : if (!y) pari_err_INV("divru",gen_0);
703 1095187859 : if (!s) return real_0_bit(expo(x) - expu(y));
704 1094551212 : if (!(y & (y-1))) /* power of 2 */
705 : {
706 92125516 : if (y == 1) return rcopy(x);
707 89412014 : return shiftr(x, -expu(y));
708 : }
709 1002425696 : e = expo(x);
710 1002425696 : lx = lg(x);
711 1002425696 : z = cgetg(lx, t_REAL);
712 1001565813 : if (lx == 3)
713 : {
714 160135495 : if (y <= uel(x,2))
715 : {
716 160135225 : hiremainder = 0;
717 160135225 : z[2] = divll(x[2],y);
718 : /* we may have hiremainder != 0 ==> garde */
719 160135225 : garde = divll(0,y);
720 : }
721 : else
722 : {
723 270 : hiremainder = x[2];
724 270 : z[2] = divll(0,y);
725 270 : garde = hiremainder;
726 270 : e -= BITS_IN_LONG;
727 : }
728 : }
729 : else
730 : {
731 841430318 : ulong yp = get_Fl_red(y);
732 842200734 : if (y <= uel(x,2))
733 : {
734 842217277 : hiremainder = 0;
735 7262652915 : for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
736 : /* we may have hiremainder != 0 ==> garde */
737 847309112 : garde = divll_pre(0,y,yp);
738 : }
739 : else
740 : {
741 3231 : long l = lx-1;
742 3231 : hiremainder = x[2];
743 56821 : for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
744 3231 : z[i] = divll_pre(0,y,yp);
745 3765 : garde = hiremainder;
746 3765 : e -= BITS_IN_LONG;
747 : }
748 : }
749 1004878206 : sh=bfffo(z[2]); /* z[2] != 0 */
750 1004878206 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
751 1005096075 : z[1] = evalsigne(s) | evalexpo(e-sh);
752 1003940410 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
753 1003834103 : return z;
754 : }
755 :
756 : GEN
757 132091855 : truedvmdii(GEN x, GEN y, GEN *z)
758 : {
759 : pari_sp av;
760 : GEN r, q, *gptr[2];
761 132091855 : if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
762 5278036 : if (z == ONLY_REM) return modii(x,y);
763 :
764 5278036 : av = avma;
765 5278036 : q = dvmdii(x,y,&r); /* assume that r is last on stack */
766 5278036 : switch(signe(r))
767 : {
768 696968 : case 0:
769 696968 : if (z) *z = gen_0;
770 696968 : return q;
771 3551941 : case 1:
772 3551941 : if (z) *z = r; else cgiv(r);
773 3551941 : return q;
774 1029127 : case -1: break;
775 : }
776 1029127 : q = addis(q, -signe(y));
777 1029127 : if (!z) return gerepileuptoint(av, q);
778 :
779 852890 : *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
780 852890 : gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(pari_sp)r,gptr,2);
781 852632 : return q;
782 : }
783 : GEN
784 128968683 : truedvmdis(GEN x, long y, GEN *z)
785 : {
786 128968683 : pari_sp av = avma;
787 : long r;
788 : GEN q;
789 :
790 128968683 : if (z == ONLY_REM) return modis(x, y);
791 128968683 : q = divis_rem(x,y,&r);
792 :
793 128945512 : if (r >= 0)
794 : {
795 117394721 : if (z) *z = utoi(r);
796 117393599 : return q;
797 : }
798 11550791 : q = gerepileuptoint(av, addis(q, (y < 0)? 1: -1));
799 11563770 : if (z) *z = utoi(r + labs(y));
800 11563741 : return q;
801 : }
802 : GEN
803 6201919 : truedvmdsi(long x, GEN y, GEN *z)
804 : {
805 : long q, r;
806 6201919 : if (z == ONLY_REM) return modsi(x, y);
807 6201919 : q = sdivsi_rem(x,y,&r);
808 6201919 : if (r >= 0) {
809 6201919 : if (z) *z = utoi(r);
810 6201919 : return stoi(q);
811 : }
812 0 : q = q - signe(y);
813 0 : if (!z) return stoi(q);
814 :
815 0 : *z = subiuspec(y+2,(ulong)-r, lgefint(y)-2);
816 0 : return stoi(q);
817 : }
818 :
819 : /* 2^n = shifti(gen_1, n) */
820 : GEN
821 102222540 : int2n(long n) {
822 : long i, m, l;
823 : GEN z;
824 102222540 : if (n < 0) return gen_0;
825 102203636 : if (n == 0) return gen_1;
826 :
827 98561927 : l = dvmdsBIL(n, &m) + 3;
828 98560000 : z = cgetipos(l);
829 539665501 : for (i = 2; i < l; i++) z[i] = 0;
830 98534745 : *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
835 25529248 : int2u(ulong n) {
836 : ulong i, m, l;
837 : GEN z;
838 25529248 : if (n == 0) return gen_1;
839 :
840 25529248 : l = dvmduBIL(n, &m) + 3;
841 25529737 : z = cgetipos(l);
842 51437080 : for (i = 2; i < l; i++) z[i] = 0;
843 25528024 : *int_MSW(z) = 1UL << m; return z;
844 : }
845 : /* 2^n - 1 */
846 : GEN
847 7504 : int2um1(ulong n) {
848 : ulong i, m, l;
849 : GEN z;
850 7504 : if (n == 0) return gen_0;
851 :
852 7504 : l = dvmduBIL(n, &m);
853 7504 : l += m? 3: 2;
854 7504 : z = cgetipos(l);
855 15147 : for (i = 2; i < l; i++) z[i] = ~0UL;
856 7504 : if (m) *int_MSW(z) = (1UL << m) - 1;
857 7504 : return z;
858 : }
859 :
860 : GEN
861 2229815619 : shifti(GEN x, long n)
862 : {
863 2229815619 : long s = signe(x);
864 : GEN y;
865 :
866 2229815619 : if(s == 0) return gen_0;
867 1927806685 : y = shiftispec(x + 2, lgefint(x) - 2, n);
868 1927691263 : if (signe(y)) setsigne(y, s);
869 1927691263 : return y;
870 : }
871 :
872 : /* actual operations will take place on a+2 and b+2: we strip the codewords */
873 : GEN
874 20728812132 : mulii(GEN a,GEN b)
875 : {
876 : long sa,sb;
877 : GEN z;
878 :
879 20728812132 : sa=signe(a); if (!sa) return gen_0;
880 14063643642 : sb=signe(b); if (!sb) return gen_0;
881 10344110498 : if (sb<0) sa = -sa;
882 10344110498 : z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
883 10374352815 : setsigne(z,sa); return z;
884 : }
885 :
886 : GEN
887 1832677688 : 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
892 238259239 : usqrt(ulong a)
893 : {
894 238259239 : ulong x = (ulong)sqrt((double)a);
895 : #ifdef LONG_IS_64BIT
896 204135963 : if (x > LOWMASK || x*x > a) x--;
897 : #endif
898 238259239 : return x;
899 : }
900 :
901 : /********************************************************************/
902 : /** **/
903 : /** EXPONENT / CONVERSION t_REAL --> double **/
904 : /** **/
905 : /********************************************************************/
906 :
907 : #ifdef LONG_IS_64BIT
908 : long
909 22287919 : dblexpo(double x)
910 : {
911 : union { double f; ulong i; } fi;
912 22287919 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
913 22287919 : const int exp_mid = 0x3ff;/* exponent bias */
914 :
915 22287919 : if (x==0.) return -exp_mid;
916 22287907 : fi.f = x;
917 22287907 : return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
918 : }
919 :
920 : ulong
921 0 : dblmantissa(double x)
922 : {
923 : union { double f; ulong i; } fi;
924 0 : const int expo_len = 11; /* number of bits of exponent */
925 :
926 0 : if (x==0.) return 0;
927 0 : fi.f = x;
928 0 : return (fi.i << expo_len) | HIGHBIT;
929 : }
930 :
931 : GEN
932 11494598 : dbltor(double x)
933 : {
934 : GEN z;
935 : long e;
936 : union { double f; ulong i; } fi;
937 11494598 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
938 11494598 : const int exp_mid = 0x3ff;/* exponent bias */
939 11494598 : const int expo_len = 11; /* number of bits of exponent */
940 :
941 11494598 : if (x==0.) return real_0_bit(-exp_mid);
942 11489278 : fi.f = x; z = cgetr(DEFAULTPREC);
943 : {
944 11486600 : const ulong a = fi.i;
945 : ulong A;
946 11486600 : e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
947 11486600 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
948 11486614 : A = a << expo_len;
949 11486614 : if (e == -exp_mid)
950 : { /* unnormalized values */
951 0 : int sh = bfffo(A);
952 0 : e -= sh-1;
953 0 : z[2] = A << sh;
954 : }
955 : else
956 11486614 : z[2] = HIGHBIT | A;
957 11486614 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
958 : }
959 11486614 : return z;
960 : }
961 :
962 : double
963 246826147 : rtodbl(GEN x)
964 : {
965 246826147 : long ex,s=signe(x);
966 : ulong a;
967 : union { double f; ulong i; } fi;
968 246826147 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
969 246826147 : const int exp_mid = 0x3ff;/* exponent bias */
970 246826147 : const int expo_len = 11; /* number of bits of exponent */
971 :
972 246826147 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
973 :
974 : /* start by rounding to closest */
975 218338232 : a = (x[2] & (HIGHBIT-1)) + 0x400;
976 218338232 : if (a & HIGHBIT) { ex++; a=0; }
977 218338232 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
978 218340923 : fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
979 218340923 : if (s<0) fi.i |= HIGHBIT;
980 218340923 : 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
994 3745890 : dblexpo(double x)
995 : {
996 : union { double f; ulong i[2]; } fi;
997 3745890 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
998 3745890 : const int exp_mid = 0x3ff;/* exponent bias */
999 3745890 : const int shift = mant_len-32;
1000 :
1001 3745890 : if (x==0.) return -exp_mid;
1002 3745888 : fi.f = x;
1003 : {
1004 3745888 : const ulong a = fi.i[INDEX0];
1005 3745888 : return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
1006 : }
1007 : }
1008 :
1009 : ulong
1010 0 : dblmantissa(double x)
1011 : {
1012 : union { double f; ulong i[2]; } fi;
1013 0 : const int expo_len = 11; /* number of bits of exponent */
1014 :
1015 0 : if (x==0.) return 0;
1016 0 : fi.f = x;
1017 : {
1018 0 : const ulong a = fi.i[INDEX0];
1019 0 : const ulong b = fi.i[INDEX1];
1020 0 : return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1021 : }
1022 : }
1023 :
1024 : GEN
1025 1876811 : dbltor(double x)
1026 : {
1027 : GEN z;
1028 : long e;
1029 : union { double f; ulong i[2]; } fi;
1030 1876811 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1031 1876811 : const int exp_mid = 0x3ff;/* exponent bias */
1032 1876811 : const int expo_len = 11; /* number of bits of exponent */
1033 1876811 : const int shift = mant_len-32;
1034 :
1035 1876811 : if (x==0.) return real_0_bit(-exp_mid);
1036 1875937 : fi.f = x; z = cgetr(DEFAULTPREC);
1037 : {
1038 1875937 : const ulong a = fi.i[INDEX0];
1039 1875937 : const ulong b = fi.i[INDEX1];
1040 : ulong A, B;
1041 1875937 : e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
1042 1875937 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
1043 1875937 : A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1044 1875937 : B = b << expo_len;
1045 1875937 : if (e == -exp_mid)
1046 : { /* unnormalized values */
1047 : int sh;
1048 0 : if (A)
1049 : {
1050 0 : sh = bfffo(A);
1051 0 : e -= sh-1;
1052 0 : z[2] = (A << sh) | (B >> (32-sh));
1053 0 : z[3] = B << sh;
1054 : }
1055 : else
1056 : {
1057 0 : sh = bfffo(B); /* B != 0 */
1058 0 : e -= sh-1 + 32;
1059 0 : z[2] = B << sh;
1060 0 : z[3] = 0;
1061 : }
1062 : }
1063 : else
1064 : {
1065 1875937 : z[3] = B;
1066 1875937 : z[2] = HIGHBIT | A;
1067 : }
1068 1875937 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
1069 : }
1070 1875937 : return z;
1071 : }
1072 :
1073 : double
1074 40483077 : rtodbl(GEN x)
1075 : {
1076 40483077 : long ex,s=signe(x),lx=lg(x);
1077 : ulong a,b,k;
1078 : union { double f; ulong i[2]; } fi;
1079 40483077 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1080 40483077 : const int exp_mid = 0x3ff;/* exponent bias */
1081 40483077 : const int expo_len = 11; /* number of bits of exponent */
1082 40483077 : const int shift = mant_len-32;
1083 :
1084 40483077 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
1085 :
1086 : /* start by rounding to closest */
1087 35665875 : a = x[2] & (HIGHBIT-1);
1088 35665875 : if (lx > 3)
1089 : {
1090 35538836 : b = x[3] + 0x400UL; if (b < 0x400UL) a++;
1091 35538836 : if (a & HIGHBIT) { ex++; a=0; }
1092 : }
1093 127039 : else b = 0;
1094 35665875 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
1095 35665875 : ex += exp_mid;
1096 35665875 : k = (a >> expo_len) | (ex << shift);
1097 35665875 : if (s<0) k |= HIGHBIT;
1098 35665875 : fi.i[INDEX0] = k;
1099 35665875 : fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
1100 35665875 : return fi.f;
1101 : }
1102 : #endif /* LONG_IS_64BIT */
1103 :
|
