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 61358082 : invmod2BIL(ulong b)
19 : {
20 : static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
21 61358082 : 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 61358082 : x = x*(2-b*x); /* one more pass necessary */
26 : #endif
27 61358082 : x = x*(2-b*x);
28 61358082 : x = x*(2-b*x); return x*(2-b*x);
29 : }
30 :
31 : void
32 1417439862 : affrr(GEN x, GEN y)
33 : {
34 1417439862 : long i, lx, ly = lg(y);
35 1417439862 : if (!signe(x))
36 : {
37 211479 : y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
38 211479 : return;
39 : }
40 1417228383 : y[1] = x[1]; lx = lg(x);
41 1417228383 : if (lx <= ly)
42 : {
43 7496604885 : for (i=2; i<lx; i++) y[i]=x[i];
44 1357207207 : for ( ; i<ly; i++) y[i]=0;
45 1195253419 : return;
46 : }
47 1167651175 : for (i=2; i<ly; i++) y[i]=x[i];
48 : /* lx > ly: round properly */
49 221974964 : if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
50 : }
51 :
52 : GEN
53 62864861 : trunc2nr(GEN x, long n)
54 : {
55 : long ex;
56 62864861 : if (!signe(x)) return gen_0;
57 62422104 : ex = expo(x) + n; if (ex < 0) return gen_0;
58 59844996 : 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 54676561 : mantissa_real(GEN x, long *e)
64 : {
65 54676561 : *e = bit_prec(x)-1-expo(x);
66 54676380 : return mantissa2nr(x, 0);
67 : }
68 :
69 : GEN
70 1159533308 : mului(ulong x, GEN y)
71 : {
72 1159533308 : long s = signe(y);
73 : GEN z;
74 :
75 1159533308 : if (!s || !x) return gen_0;
76 1004933485 : z = muluispec(x, y+2, lgefint(y)-2);
77 1004552933 : setsigne(z,s); return z;
78 : }
79 :
80 : GEN
81 621251646 : mulsi(long x, GEN y)
82 : {
83 621251646 : long s = signe(y);
84 : GEN z;
85 :
86 621251646 : if (!s || !x) return gen_0;
87 372499029 : if (x<0) { s = -s; x = -x; }
88 372499029 : z = muluispec((ulong)x, y+2, lgefint(y)-2);
89 372495039 : setsigne(z,s); return z;
90 : }
91 :
92 : GEN
93 201177897 : mulss(long x, long y)
94 : {
95 : ulong p1;
96 : LOCAL_HIREMAINDER;
97 :
98 201177897 : if (!x || !y) return gen_0;
99 200841389 : if (x<0) {
100 413727 : x = -x;
101 413727 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
102 299442 : p1 = mulll(x,y); return uutoineg(hiremainder, p1);
103 : } else {
104 200427662 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
105 200384153 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
106 : }
107 : }
108 : GEN
109 732393 : sqrs(long x)
110 : {
111 : ulong p1;
112 : LOCAL_HIREMAINDER;
113 :
114 732393 : if (!x) return gen_0;
115 730622 : if (x<0) x = -x;
116 730622 : p1 = mulll(x,x); return uutoi(hiremainder, p1);
117 : }
118 : GEN
119 4713164884 : muluu(ulong x, ulong y)
120 : {
121 : ulong p1;
122 : LOCAL_HIREMAINDER;
123 :
124 4713164884 : if (!x || !y) return gen_0;
125 4714142074 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
126 : }
127 : GEN
128 647310612 : sqru(ulong x)
129 : {
130 : ulong p1;
131 : LOCAL_HIREMAINDER;
132 :
133 647310612 : if (!x) return gen_0;
134 646912842 : 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 752629560 : mulur_2(ulong x, GEN y, long s)
140 : {
141 : long m, sh, i, lx, e;
142 : GEN z;
143 : ulong garde;
144 : LOCAL_HIREMAINDER;
145 :
146 752629560 : if (!(x & (x-1))) { z = shiftr(y, expu(x)); setsigne(z, s); return z; }
147 582915357 : lx = lg(y); z = cgetg(lx, t_REAL); e = expo(y);
148 582476816 : y--; garde = mulll(x,y[lx]);
149 2914366137 : for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
150 582476816 : z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
151 582476816 : sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
152 582476816 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
153 583350355 : z[1] = evalsigne(s) | evalexpo(m+e);
154 585073254 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
155 584836619 : return z;
156 : }
157 :
158 : INLINE GEN
159 481861 : mul0r(GEN x)
160 : {
161 481861 : long l = realprec(x), e = expo(x);
162 481861 : e = (l > 0)? e - l: (e < 0? 2*e: 0);
163 481861 : return real_0_bit(e);
164 : }
165 : /* lg(x) > 2 */
166 : INLINE GEN
167 0 : div0r(GEN x) {
168 0 : long l = realprec(x), e = expo(x);
169 0 : return real_0_bit(-l - e);
170 : }
171 :
172 : GEN
173 132258386 : mulsr(long x, GEN y)
174 : {
175 : long s;
176 :
177 132258386 : if (!x) return mul0r(y);
178 132258344 : s = signe(y);
179 132258344 : if (!s)
180 : {
181 247313 : if (x < 0) x = -x;
182 247313 : return real_0_bit( expo(y) + expu(x) );
183 : }
184 132011031 : if (x==1) return rcopy(y);
185 123220174 : if (x==-1) return negr(y);
186 119659141 : if (x < 0)
187 42647034 : return mulur_2((ulong)-x, y, -s);
188 : else
189 77012107 : return mulur_2((ulong)x, y, s);
190 : }
191 :
192 : GEN
193 967800140 : mulur(ulong x, GEN y)
194 : {
195 : long s;
196 :
197 967800140 : if (!x) return mul0r(y);
198 967800133 : s = signe(y);
199 967800133 : if (!s) return real_0_bit( expo(y) + expu(x) );
200 959836024 : if (x==1) return rcopy(y);
201 632478198 : return mulur_2(x, y, s);
202 : }
203 :
204 : INLINE void
205 3768480660 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
206 : {
207 : long i;
208 3768480660 : if (hi[2] < 0)
209 : {
210 1739004226 : if (z != hi)
211 1712055997 : for (i=2; i<lz ; i++) z[i] = hi[i];
212 1739004226 : ez++;
213 : }
214 : else
215 : {
216 2029476434 : shift_left(z,hi,2,lz-1, garde, 1);
217 2032916964 : garde <<= 1;
218 : }
219 3771921190 : if (garde & HIGHBIT)
220 : { /* round to nearest */
221 1848675265 : i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
222 1830238091 : if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
223 : }
224 3771921190 : z[1] = evalsigne(sz)|evalexpo(ez);
225 3767005557 : }
226 : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
227 : INLINE void
228 786285438 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
229 : {
230 786285438 : if (z[2] < 0)
231 : { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
232 367388226 : if (garde & HIGHBIT) z[2]++; /* round properly */
233 367388226 : ez++;
234 : }
235 : else
236 : {
237 418897212 : uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
238 418897212 : if (garde & (1UL<<(BITS_IN_LONG-2)))
239 : {
240 172317953 : uel(z,2)++; /* round properly, z2+1 can overflow */
241 172317953 : if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
242 : }
243 : }
244 786285438 : z[1] = evalsigne(sz)|evalexpo(ez);
245 786287698 : }
246 :
247 : /* set z <-- x^2 != 0, floating point multiplication.
248 : * lz = lg(z) = lg(x) */
249 : INLINE void
250 567252614 : sqrz_i(GEN z, GEN x, long lz)
251 : {
252 567252614 : long ez = 2*expo(x);
253 : long i, j, lzz, p1;
254 : ulong garde;
255 : GEN x1;
256 : LOCAL_HIREMAINDER;
257 : LOCAL_OVERFLOW;
258 :
259 567252614 : if (lz > prec2lg(SQRR_SQRI_LIMIT))
260 : {
261 39973569 : pari_sp av = avma;
262 39973569 : GEN hi = sqrispec_mirror(x+2, lz-2);
263 40278555 : mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
264 40263201 : set_avma(av); return;
265 : }
266 527192778 : if (lz == 3)
267 : {
268 122862996 : garde = mulll(x[2],x[2]);
269 122862996 : z[2] = hiremainder;
270 122862996 : mulrrz_3end(z, 1, ez, garde);
271 122843165 : return;
272 : }
273 :
274 404329782 : lzz = lz-1; p1 = x[lzz];
275 404329782 : if (p1)
276 : {
277 367974935 : (void)mulll(p1,x[3]);
278 367974935 : garde = addmul(p1,x[2]);
279 367974935 : z[lzz] = hiremainder;
280 : }
281 : else
282 : {
283 36354847 : garde = 0;
284 36354847 : z[lzz] = 0;
285 : }
286 1843842595 : for (j=lz-2, x1=x-j; j>=3; j--)
287 : {
288 1439512813 : p1 = x[j]; x1++;
289 1439512813 : if (p1)
290 : {
291 1411592481 : (void)mulll(p1,x1[lz+1]);
292 1411592481 : garde = addll(addmul(p1,x1[lz]), garde);
293 7180180285 : for (i=lzz; i>j; i--)
294 : {
295 5768587804 : hiremainder += overflow;
296 5768587804 : z[i] = addll(addmul(p1,x1[i]), z[i]);
297 : }
298 1411592481 : z[j] = hiremainder+overflow;
299 : }
300 27920332 : else z[j]=0;
301 : }
302 404329782 : p1 = x[2]; x1++;
303 404329782 : garde = addll(mulll(p1,x1[lz]), garde);
304 2253140661 : for (i=lzz; i>2; i--)
305 : {
306 1848810879 : hiremainder += overflow;
307 1848810879 : z[i] = addll(addmul(p1,x1[i]), z[i]);
308 : }
309 404329782 : z[2] = hiremainder+overflow;
310 404329782 : mulrrz_end(z, z, lz, 1, ez, garde);
311 : }
312 :
313 : /* lz "large" = lg(y) = lg(z), lg(x) > lz if flag = 1 and >= if flag = 0 */
314 : INLINE void
315 54102158 : mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
316 : {
317 54102158 : pari_sp av = avma;
318 54102158 : GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
319 54102158 : mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
320 54102158 : set_avma(av);
321 54102158 : }
322 :
323 : /* lz = 3 */
324 : INLINE void
325 664407617 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
326 : {
327 : ulong garde;
328 : LOCAL_HIREMAINDER;
329 664407617 : if (flag)
330 : {
331 79496014 : (void)mulll(x[2],y[3]);
332 79496014 : garde = addmul(x[2],y[2]);
333 : }
334 : else
335 584911603 : garde = mulll(x[2],y[2]);
336 664407617 : z[2] = hiremainder;
337 664407617 : mulrrz_3end(z, sz, expo(x)+expo(y), garde);
338 664514275 : }
339 :
340 : /* set z <-- x*y, floating point multiplication. Trailing 0s for x are
341 : * treated efficiently (important application: mulir).
342 : * lz = lg(z) = lg(x) <= ly <= lg(y), sz = signe(z). flag = lg(x) < lg(y) */
343 : INLINE void
344 3922022647 : mulrrz_i(GEN z, GEN x, GEN y, long lz, long flag, long sz)
345 : {
346 : long ez, i, j, lzz, p1;
347 : ulong garde;
348 : GEN y1;
349 : LOCAL_HIREMAINDER;
350 : LOCAL_OVERFLOW;
351 :
352 3922022647 : if (x == y) { sqrz_i(z,x,lz); return; }
353 3922022647 : if (lz > prec2lg(MULRR_MULII_LIMIT)) { mulrrz_int(z,x,y,lz,flag,sz); return; }
354 3871788269 : if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
355 3207359856 : ez = expo(x) + expo(y);
356 3207359856 : if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
357 3207359856 : lzz=lz-1; p1=x[lzz];
358 3207359856 : if (p1)
359 : {
360 2993990389 : (void)mulll(p1,y[3]);
361 2993990389 : garde = addll(addmul(p1,y[2]), garde);
362 2993990389 : z[lzz] = overflow+hiremainder;
363 : }
364 213369467 : else z[lzz]=0;
365 14050667410 : for (j=lz-2, y1=y-j; j>=3; j--)
366 : {
367 10843307554 : p1 = x[j]; y1++;
368 10843307554 : if (p1)
369 : {
370 10531423194 : (void)mulll(p1,y1[lz+1]);
371 10531423194 : garde = addll(addmul(p1,y1[lz]), garde);
372 63330548957 : for (i=lzz; i>j; i--)
373 : {
374 52799125763 : hiremainder += overflow;
375 52799125763 : z[i] = addll(addmul(p1,y1[i]), z[i]);
376 : }
377 10531423194 : z[j] = hiremainder+overflow;
378 : }
379 311884360 : else z[j]=0;
380 : }
381 3207359856 : p1 = x[2]; y1++;
382 3207359856 : garde = addll(mulll(p1,y1[lz]), garde);
383 17312727111 : for (i=lzz; i>2; i--)
384 : {
385 14105367255 : hiremainder += overflow;
386 14105367255 : z[i] = addll(addmul(p1,y1[i]), z[i]);
387 : }
388 3207359856 : z[2] = hiremainder+overflow;
389 3207359856 : mulrrz_end(z, z, lz, sz, ez, garde);
390 : }
391 :
392 : GEN
393 4178890114 : mulrr(GEN x, GEN y)
394 : {
395 : long flag, ly, lz, sx, sy;
396 : GEN z;
397 :
398 4178890114 : if (x == y) return sqrr(x);
399 4177270404 : sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
400 3902126931 : sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
401 3864443357 : if (sy < 0) sx = -sx;
402 3864443357 : lz = lg(x);
403 3864443357 : ly = lg(y);
404 3864443357 : if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
405 3864443357 : z = cgetg(lz, t_REAL);
406 3859611362 : mulrrz_i(z, x,y, lz,flag, sx);
407 3868208922 : return z;
408 : }
409 :
410 : GEN
411 601847839 : sqrr(GEN x)
412 : {
413 601847839 : long lz, sx = signe(x);
414 : GEN z;
415 :
416 601847839 : if (!sx) return real_0_bit(2*expo(x));
417 567829235 : lz = lg(x); z = cgetg(lz, t_REAL);
418 567267390 : sqrz_i(z, x, lz);
419 568102477 : return z;
420 : }
421 :
422 : GEN
423 1019199562 : mulir(GEN x, GEN y)
424 : {
425 1019199562 : long sx = signe(x), sy;
426 1019199562 : if (!sx) return mul0r(y);
427 1018717750 : if (lgefint(x) == 3) {
428 905901635 : GEN z = mulur(uel(x,2), y);
429 907013438 : if (sx < 0) togglesign(z);
430 906742307 : return z;
431 : }
432 112816115 : sy = signe(y);
433 112816115 : if (!sy) return real_0_bit(expi(x) + expo(y));
434 111937334 : if (sy < 0) sx = -sx;
435 : {
436 111937334 : long lz = lg(y), lx = lgefint(x);
437 111937334 : GEN hi, z = cgetg(lz, t_REAL);
438 111481386 : pari_sp av = avma;
439 111481386 : if (lx < (lz>>1) || (lx < lz && lz > prec2lg(MULRR_MULII_LIMIT)))
440 : { /* size mantissa of x < half size of mantissa z, or lx < lz so large
441 : * that mulrr will call mulii anyway: mulii */
442 48167314 : x = itor(x, lg2prec(lx));
443 47705586 : hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
444 49490183 : mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
445 : }
446 : else /* dubious: complete x with 0s and call mulrr */
447 63279682 : mulrrz_i(z, itor(x, lg2prec(lz)), y, lz, 0, sx);
448 112998054 : set_avma(av); return z;
449 : }
450 : }
451 :
452 : /* x + y*z, generic. If lgefint(z) <= 3, caller should use faster variants */
453 : static GEN
454 92949251 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
455 : {
456 92949251 : long lx = lgefint(x), ly;
457 : pari_sp av;
458 : GEN t;
459 92949251 : if (lx == 2) return mulii(z,y);
460 90127483 : ly = lgefint(y);
461 90127483 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
462 89509830 : av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
463 89511319 : t = mulii(z, y);
464 89502979 : set_avma(av); return addii(t,x);
465 : }
466 : /* x + y*z, lgefint(z) == 3 */
467 : static GEN
468 482686299 : addmulii_lg3(GEN x, GEN y, GEN z)
469 : {
470 482686299 : long s = signe(z), lx, ly;
471 482686299 : ulong w = z[2];
472 : pari_sp av;
473 : GEN t;
474 482686299 : if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
475 396638969 : lx = lgefint(x);
476 396638969 : ly = lgefint(y);
477 396638969 : if (lx == 2)
478 : { /* x = 0 */
479 73358280 : if (ly == 2) return gen_0;
480 40858003 : t = muluispec(w, y+2, ly-2);
481 40858268 : if (signe(y) < 0) s = -s;
482 40858268 : setsigne(t, s); return t;
483 : }
484 323280689 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
485 282465051 : av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
486 282612056 : t = muluispec(w, y+2, ly-2);
487 282563958 : if (signe(y) < 0) s = -s;
488 282563958 : setsigne(t, s);
489 282563958 : set_avma(av); return addii(x,t);
490 : }
491 : /* x + y*z */
492 : GEN
493 169818298 : addmulii(GEN x, GEN y, GEN z)
494 : {
495 169818298 : long lz = lgefint(z);
496 169818298 : switch(lz)
497 : {
498 759297 : case 2: return icopy(x); /* z = 0, wasteful copy */
499 128879453 : case 3: return addmulii_lg3(x, y, z);
500 40179548 : default:return addmulii_gen(x, y, z, lz);
501 : }
502 : }
503 : /* x + y*z, returns x itself and not a copy when y*z = 0 */
504 : GEN
505 1630385136 : addmulii_inplace(GEN x, GEN y, GEN z)
506 : {
507 : long lz;
508 1630385136 : if (lgefint(y) == 2) return x;
509 406342419 : lz = lgefint(z);
510 406342419 : switch(lz)
511 : {
512 19181 : case 2: return x;
513 353837704 : case 3: return addmulii_lg3(x, y, z);
514 52485534 : default:return addmulii_gen(x, y, z, lz);
515 : }
516 : }
517 :
518 : /* written by Bruno Haible following an idea of Robert Harley */
519 : long
520 2374485170 : vals(ulong z)
521 : {
522 : 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};
523 : #ifdef LONG_IS_64BIT
524 : long s;
525 : #endif
526 :
527 2374485170 : if (!z) return -1;
528 : #ifdef LONG_IS_64BIT
529 2133853117 : if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
530 : #endif
531 2374485156 : z |= ~z + 1;
532 2374485156 : z += z << 4;
533 2374485156 : z += z << 6;
534 2374485156 : z ^= z << 16; /* or z -= z<<16 */
535 : #ifdef LONG_IS_64BIT
536 2133853117 : return s + tab[(z&0xffffffff)>>26];
537 : #else
538 240632039 : return tab[z>>26];
539 : #endif
540 : }
541 :
542 : GEN
543 175 : divsi(long x, GEN y)
544 : {
545 175 : long p1, s = signe(y);
546 : LOCAL_HIREMAINDER;
547 :
548 175 : if (!s) pari_err_INV("divsi",gen_0);
549 175 : if (!x || lgefint(y)>3 || ((long)y[2])<0) return gen_0;
550 175 : hiremainder=0; p1=divll(labs(x),y[2]);
551 175 : if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
552 175 : if (s<0) p1 = -p1;
553 175 : return stoi(p1);
554 : }
555 :
556 : GEN
557 3278531 : divir(GEN x, GEN y)
558 : {
559 : GEN z;
560 3278531 : long ly = lg(y), lx = lgefint(x);
561 : pari_sp av;
562 :
563 3278531 : if (ly == 2) pari_err_INV("divir",y);
564 3278531 : if (lx == 2) return div0r(y);
565 3278531 : if (lx == 3) {
566 1913520 : z = divur(x[2], y);
567 1913518 : if (signe(x) < 0) togglesign(z);
568 1913518 : return z;
569 : }
570 1365011 : z = cgetg(ly, t_REAL); av = avma;
571 1365010 : affrr(divrr(itor(x, lg2prec(ly+1)), y), z);
572 1365012 : set_avma(av); return z;
573 : }
574 :
575 : GEN
576 2456204 : divur(ulong x, GEN y)
577 : {
578 : pari_sp av;
579 2456204 : long p = realprec(y);
580 : GEN z;
581 :
582 2456204 : if (p == 0) pari_err_INV("divur",y);
583 2456204 : if (!x) return div0r(y);
584 2456204 : if (x == 1) return invr(y);
585 2446362 : if (!(x & (x-1))) /* power of 2 */
586 : {
587 913185 : z = invr(y);
588 913182 : shiftr_inplace(z, expu(x)); return z;
589 : }
590 1533177 : if (p > INVNEWTON_LIMIT) {
591 4 : av = avma; z = invr(y);
592 4 : if (x == 1) return z;
593 4 : return gc_leaf(av, mulur(x, z));
594 : }
595 1533173 : z = cgetr(p); av = avma;
596 1533182 : affrr(divrr(utor(x, p + BITS_IN_LONG), y), z);
597 1533185 : set_avma(av); return z;
598 : }
599 :
600 : GEN
601 798 : divsr(long x, GEN y)
602 : {
603 798 : if (x >= 0) return divur((ulong)x, y);
604 0 : y = divur((ulong)(-x), y); togglesign(y); return y;
605 : }
606 :
607 : /* returns 1/y, assume y != 0 */
608 : static GEN
609 74731338 : invr_basecase(GEN y)
610 : {
611 74731338 : long p = realprec(y);
612 74731338 : GEN z = cgetr(p);
613 74722330 : pari_sp av = avma;
614 74722330 : affrr(divrr(real_1(p + BITS_IN_LONG), y), z);
615 74730486 : set_avma(av); return z;
616 : }
617 : /* returns 1/b, Newton iteration */
618 : GEN
619 74731522 : invr(GEN b)
620 : {
621 74731522 : const long s = 6;
622 74731522 : long i, p, l = lg(b);
623 : GEN x, a;
624 : ulong mask;
625 :
626 74731522 : if (l <= maxss(prec2lg(INVNEWTON_LIMIT), (1L<<s) + 2)) {
627 74727092 : if (l == 2) pari_err_INV("invr",b);
628 74727092 : return invr_basecase(b);
629 : }
630 4354 : mask = quadratic_prec_mask(l-2);
631 30478 : for(i=0, p=1; i<s; i++) { p <<= 1; if (mask & 1) p--; mask >>= 1; }
632 4354 : x = cgetg(l, t_REAL);
633 4354 : a = rcopy(b); a[1] = _evalexpo(0) | evalsigne(1);
634 4354 : affrr(invr_basecase(rtor(a, lg2prec(p+2))), x);
635 13922 : while (mask > 1)
636 : {
637 9568 : p <<= 1; if (mask & 1) p--;
638 9568 : mask >>= 1;
639 9568 : setlg(a, p + 2);
640 9568 : setlg(x, p + 2);
641 : /* TODO: mulrr(a,x) should be a half product (the higher half is known).
642 : * mulrr(x, ) already is */
643 9568 : affrr(addrr(x, mulrr(x, subsr(1, mulrr(a,x)))), x);
644 9568 : set_avma((pari_sp)a);
645 : }
646 4354 : x[1] = (b[1] & SIGNBITS) | evalexpo(expo(x)-expo(b));
647 4354 : set_avma((pari_sp)x); return x;
648 : }
649 :
650 : GEN
651 3129870857 : modii(GEN x, GEN y)
652 : {
653 3129870857 : switch(signe(x))
654 : {
655 494992597 : case 0: return gen_0;
656 2078862395 : case 1: return remii(x,y);
657 556015865 : default:
658 : {
659 556015865 : pari_sp av = avma;
660 556015865 : (void)new_chunk(lgefint(y));
661 560817867 : x = remii(x,y); set_avma(av);
662 560837515 : if (x==gen_0) return x;
663 518906101 : return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
664 : }
665 : }
666 : }
667 :
668 : GEN
669 20038193 : divrs(GEN x, long y)
670 : {
671 : GEN z;
672 20038193 : if (y < 0)
673 : {
674 5905290 : z = divru(x, (ulong)-y);
675 5905266 : togglesign(z);
676 : }
677 : else
678 14132903 : z = divru(x, (ulong)y);
679 20082825 : return z;
680 : }
681 :
682 : GEN
683 1151967748 : divru(GEN x, ulong y)
684 : {
685 1151967748 : long i, lx, sh, e, s = signe(x);
686 : ulong garde;
687 : GEN z;
688 : LOCAL_HIREMAINDER;
689 :
690 1151967748 : if (!y) pari_err_INV("divru",gen_0);
691 1152063112 : if (!s) return real_0_bit(expo(x) - expu(y));
692 1151368470 : if (!(y & (y-1))) /* power of 2 */
693 : {
694 98770594 : if (y == 1) return rcopy(x);
695 96031167 : return shiftr(x, -expu(y));
696 : }
697 1052597876 : e = expo(x);
698 1052597876 : lx = lg(x);
699 1052597876 : z = cgetg(lx, t_REAL);
700 1051606838 : if (lx == 3)
701 : {
702 175001398 : if (y <= uel(x,2))
703 : {
704 175001128 : hiremainder = 0;
705 175001128 : z[2] = divll(x[2],y);
706 : /* we may have hiremainder != 0 ==> garde */
707 175001128 : garde = divll(0,y);
708 : }
709 : else
710 : {
711 270 : hiremainder = x[2];
712 270 : z[2] = divll(0,y);
713 270 : garde = hiremainder;
714 270 : e -= BITS_IN_LONG;
715 : }
716 : }
717 : else
718 : {
719 876605440 : ulong yp = get_Fl_red(y);
720 877422980 : if (y <= uel(x,2))
721 : {
722 877433822 : hiremainder = 0;
723 7458549880 : for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
724 : /* we may have hiremainder != 0 ==> garde */
725 884928337 : garde = divll_pre(0,y,yp);
726 : }
727 : else
728 : {
729 3232 : long l = lx-1;
730 3232 : hiremainder = x[2];
731 56829 : for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
732 3232 : z[i] = divll_pre(0,y,yp);
733 3766 : garde = hiremainder;
734 3766 : e -= BITS_IN_LONG;
735 : }
736 : }
737 1055466842 : sh=bfffo(z[2]); /* z[2] != 0 */
738 1055466842 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
739 1055353433 : z[1] = evalsigne(s) | evalexpo(e-sh);
740 1053988071 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
741 1053890468 : return z;
742 : }
743 :
744 : GEN
745 142311790 : truedvmdii(GEN x, GEN y, GEN *z)
746 : {
747 : pari_sp av;
748 : GEN r, q;
749 142311790 : if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
750 5488328 : if (z == ONLY_REM) return modii(x,y);
751 :
752 5488328 : av = avma;
753 5488328 : q = dvmdii(x,y,&r); /* assume that r is last on stack */
754 5488326 : switch(signe(r))
755 : {
756 711805 : case 0:
757 711805 : if (z) *z = gen_0;
758 711805 : return q;
759 3735432 : case 1:
760 3735432 : if (z) *z = r; else cgiv(r);
761 3735432 : return q;
762 1041089 : case -1: break;
763 : }
764 1041089 : q = addis(q, -signe(y));
765 1041089 : if (!z) return gc_INT(av, q);
766 :
767 864991 : *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
768 864991 : return gc_all_unsafe(av,(pari_sp)r,2,&q,z);
769 : }
770 : GEN
771 139044627 : truedvmdis(GEN x, long y, GEN *z)
772 : {
773 139044627 : pari_sp av = avma;
774 : long r;
775 : GEN q;
776 :
777 139044627 : if (z == ONLY_REM) return modis(x, y);
778 139044627 : q = divis_rem(x,y,&r);
779 :
780 138973612 : if (r >= 0)
781 : {
782 126121296 : if (z) *z = utoi(r);
783 126120023 : return q;
784 : }
785 12852316 : q = gc_INT(av, addis(q, (y < 0)? 1: -1));
786 12886809 : if (z) *z = utoi(r + labs(y));
787 12886768 : return q;
788 : }
789 : GEN
790 6202318 : truedvmdsi(long x, GEN y, GEN *z)
791 : {
792 : long q, r;
793 6202318 : if (z == ONLY_REM) return modsi(x, y);
794 6202318 : q = sdivsi_rem(x,y,&r);
795 6202318 : if (r >= 0) {
796 6202318 : if (z) *z = utoi(r);
797 6202318 : return stoi(q);
798 : }
799 0 : q = q - signe(y);
800 0 : if (!z) return stoi(q);
801 :
802 0 : *z = subiuspec(y+2,(ulong)-r, lgefint(y)-2);
803 0 : return stoi(q);
804 : }
805 :
806 : /* 2^n = shifti(gen_1, n) */
807 : GEN
808 109704935 : int2n(long n) {
809 : long i, m, l;
810 : GEN z;
811 109704935 : if (n < 0) return gen_0;
812 109685923 : if (n == 0) return gen_1;
813 :
814 105931451 : l = dvmdsBIL(n, &m) + 3;
815 105933097 : z = cgetipos(l);
816 568594122 : for (i = 2; i < l; i++) z[i] = 0;
817 105905492 : *int_MSW(z) = 1UL << m; return z;
818 : }
819 : /* To avoid problems when 2^(BIL-1) < n. Overflow cleanly, where int2n
820 : * returns gen_0 */
821 : GEN
822 25993642 : int2u(ulong n) {
823 : ulong i, m, l;
824 : GEN z;
825 25993642 : if (n == 0) return gen_1;
826 :
827 25993642 : l = dvmduBIL(n, &m) + 3;
828 25994128 : z = cgetipos(l);
829 52378410 : for (i = 2; i < l; i++) z[i] = 0;
830 25992404 : *int_MSW(z) = 1UL << m; return z;
831 : }
832 : /* 2^n - 1 */
833 : GEN
834 8807 : int2um1(ulong n) {
835 : ulong i, m, l;
836 : GEN z;
837 8807 : if (n == 0) return gen_0;
838 :
839 8807 : l = dvmduBIL(n, &m);
840 8807 : l += m? 3: 2;
841 8807 : z = cgetipos(l);
842 17754 : for (i = 2; i < l; i++) z[i] = ~0UL;
843 8807 : if (m) *int_MSW(z) = (1UL << m) - 1;
844 8807 : return z;
845 : }
846 :
847 : GEN
848 2313814404 : shifti(GEN x, long n)
849 : {
850 2313814404 : long s = signe(x);
851 : GEN y;
852 :
853 2313814404 : if(s == 0) return gen_0;
854 2006227739 : y = shiftispec(x + 2, lgefint(x) - 2, n);
855 2005931946 : if (signe(y)) setsigne(y, s);
856 2005931946 : return y;
857 : }
858 :
859 : /* actual operations will take place on a+2 and b+2: we strip the codewords */
860 : GEN
861 21505493369 : mulii(GEN a,GEN b)
862 : {
863 : long sa,sb;
864 : GEN z;
865 :
866 21505493369 : sa=signe(a); if (!sa) return gen_0;
867 14677399201 : sb=signe(b); if (!sb) return gen_0;
868 10723259460 : if (sb<0) sa = -sa;
869 10723259460 : z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
870 10772264103 : setsigne(z,sa); return z;
871 : }
872 :
873 : GEN
874 1955255248 : sqri(GEN a) { return sqrispec(a+2, lgefint(a)-2); }
875 :
876 : /* sqrt()'s result may be off by 1 when a is not representable exactly as a
877 : * double [64bit machine] */
878 : ulong
879 252964488 : usqrt(ulong a)
880 : {
881 252964488 : ulong x = (ulong)sqrt((double)a);
882 : #ifdef LONG_IS_64BIT
883 219206661 : if (x > LOWMASK || x*x > a) x--;
884 : #endif
885 252964488 : return x;
886 : }
887 :
888 : /********************************************************************/
889 : /** **/
890 : /** EXPONENT / CONVERSION t_REAL --> double **/
891 : /** **/
892 : /********************************************************************/
893 :
894 : #ifdef LONG_IS_64BIT
895 : long
896 23735539 : dblexpo(double x)
897 : {
898 : union { double f; ulong i; } fi;
899 23735539 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
900 23735539 : const int exp_mid = 0x3ff;/* exponent bias */
901 :
902 23735539 : if (x==0.) return -exp_mid;
903 23735527 : fi.f = x;
904 23735527 : return ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
905 : }
906 :
907 : ulong
908 0 : dblmantissa(double x)
909 : {
910 : union { double f; ulong i; } fi;
911 0 : const int expo_len = 11; /* number of bits of exponent */
912 :
913 0 : if (x==0.) return 0;
914 0 : fi.f = x;
915 0 : return (fi.i << expo_len) | HIGHBIT;
916 : }
917 :
918 : GEN
919 11703308 : dbltor(double x)
920 : {
921 : GEN z;
922 : long e;
923 : union { double f; ulong i; } fi;
924 11703308 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
925 11703308 : const int exp_mid = 0x3ff;/* exponent bias */
926 11703308 : const int expo_len = 11; /* number of bits of exponent */
927 :
928 11703308 : if (x==0.) return real_0_bit(-exp_mid);
929 11696986 : fi.f = x; z = cgetr(DEFAULTPREC);
930 : {
931 11695181 : const ulong a = fi.i;
932 : ulong A;
933 11695181 : e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
934 11695181 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
935 11695077 : A = a << expo_len;
936 11695077 : if (e == -exp_mid)
937 : { /* unnormalized values */
938 0 : int sh = bfffo(A);
939 0 : e -= sh-1;
940 0 : z[2] = A << sh;
941 : }
942 : else
943 11695077 : z[2] = HIGHBIT | A;
944 11695077 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
945 : }
946 11695077 : return z;
947 : }
948 :
949 : double
950 267907771 : rtodbl(GEN x)
951 : {
952 267907771 : long ex,s=signe(x);
953 : ulong a;
954 : union { double f; ulong i; } fi;
955 267907771 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
956 267907771 : const int exp_mid = 0x3ff;/* exponent bias */
957 267907771 : const int expo_len = 11; /* number of bits of exponent */
958 :
959 267907771 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
960 :
961 : /* start by rounding to closest */
962 236005533 : a = (x[2] & (HIGHBIT-1)) + 0x400;
963 236005533 : if (a & HIGHBIT) { ex++; a=0; }
964 236005533 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
965 236008038 : fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
966 236008038 : if (s<0) fi.i |= HIGHBIT;
967 236008038 : return fi.f;
968 : }
969 :
970 : #else /* LONG_IS_64BIT */
971 :
972 : #if PARI_DOUBLE_FORMAT == 1
973 : # define INDEX0 1
974 : # define INDEX1 0
975 : #elif PARI_DOUBLE_FORMAT == 0
976 : # define INDEX0 0
977 : # define INDEX1 1
978 : #endif
979 :
980 : long
981 3990156 : dblexpo(double x)
982 : {
983 : union { double f; ulong i[2]; } fi;
984 3990156 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
985 3990156 : const int exp_mid = 0x3ff;/* exponent bias */
986 3990156 : const int shift = mant_len-32;
987 :
988 3990156 : if (x==0.) return -exp_mid;
989 3990154 : fi.f = x;
990 : {
991 3990154 : const ulong a = fi.i[INDEX0];
992 3990154 : return ((a & (HIGHBIT-1)) >> shift) - exp_mid;
993 : }
994 : }
995 :
996 : ulong
997 0 : dblmantissa(double x)
998 : {
999 : union { double f; ulong i[2]; } fi;
1000 0 : const int expo_len = 11; /* number of bits of exponent */
1001 :
1002 0 : if (x==0.) return 0;
1003 0 : fi.f = x;
1004 : {
1005 0 : const ulong a = fi.i[INDEX0];
1006 0 : const ulong b = fi.i[INDEX1];
1007 0 : return HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1008 : }
1009 : }
1010 :
1011 : GEN
1012 1897392 : dbltor(double x)
1013 : {
1014 : GEN z;
1015 : long e;
1016 : union { double f; ulong i[2]; } fi;
1017 1897392 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1018 1897392 : const int exp_mid = 0x3ff;/* exponent bias */
1019 1897392 : const int expo_len = 11; /* number of bits of exponent */
1020 1897392 : const int shift = mant_len-32;
1021 :
1022 1897392 : if (x==0.) return real_0_bit(-exp_mid);
1023 1896351 : fi.f = x; z = cgetr(DEFAULTPREC);
1024 : {
1025 1896351 : const ulong a = fi.i[INDEX0];
1026 1896351 : const ulong b = fi.i[INDEX1];
1027 : ulong A, B;
1028 1896351 : e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
1029 1896351 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
1030 1896351 : A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1031 1896351 : B = b << expo_len;
1032 1896351 : if (e == -exp_mid)
1033 : { /* unnormalized values */
1034 : int sh;
1035 0 : if (A)
1036 : {
1037 0 : sh = bfffo(A);
1038 0 : e -= sh-1;
1039 0 : z[2] = (A << sh) | (B >> (32-sh));
1040 0 : z[3] = B << sh;
1041 : }
1042 : else
1043 : {
1044 0 : sh = bfffo(B); /* B != 0 */
1045 0 : e -= sh-1 + 32;
1046 0 : z[2] = B << sh;
1047 0 : z[3] = 0;
1048 : }
1049 : }
1050 : else
1051 : {
1052 1896351 : z[3] = B;
1053 1896351 : z[2] = HIGHBIT | A;
1054 : }
1055 1896351 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
1056 : }
1057 1896351 : return z;
1058 : }
1059 :
1060 : double
1061 44211702 : rtodbl(GEN x)
1062 : {
1063 44211702 : long ex,s=signe(x),lx=lg(x);
1064 : ulong a,b,k;
1065 : union { double f; ulong i[2]; } fi;
1066 44211702 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1067 44211702 : const int exp_mid = 0x3ff;/* exponent bias */
1068 44211702 : const int expo_len = 11; /* number of bits of exponent */
1069 44211702 : const int shift = mant_len-32;
1070 :
1071 44211702 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
1072 :
1073 : /* start by rounding to closest */
1074 38798733 : a = x[2] & (HIGHBIT-1);
1075 38798733 : if (lx > 3)
1076 : {
1077 38671415 : b = x[3] + 0x400UL; if (b < 0x400UL) a++;
1078 38671415 : if (a & HIGHBIT) { ex++; a=0; }
1079 : }
1080 127318 : else b = 0;
1081 38798733 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
1082 38798733 : ex += exp_mid;
1083 38798733 : k = (a >> expo_len) | (ex << shift);
1084 38798733 : if (s<0) k |= HIGHBIT;
1085 38798733 : fi.i[INDEX0] = k;
1086 38798733 : fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
1087 38798733 : return fi.f;
1088 : }
1089 : #endif /* LONG_IS_64BIT */
1090 :
|
