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 61358106 : invmod2BIL(ulong b)
19 : {
20 : static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
21 61358106 : 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 61358106 : x = x*(2-b*x); /* one more pass necessary */
26 : #endif
27 61358106 : x = x*(2-b*x);
28 61358106 : x = x*(2-b*x); return x*(2-b*x);
29 : }
30 :
31 : void
32 1424430326 : affrr(GEN x, GEN y)
33 : {
34 1424430326 : long i, lx, ly = lg(y);
35 1424430326 : if (!signe(x))
36 : {
37 211479 : y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
38 211479 : return;
39 : }
40 1424218847 : y[1] = x[1]; lx = lg(x);
41 1424218847 : if (lx <= ly)
42 : {
43 7528613622 : for (i=2; i<lx; i++) y[i]=x[i];
44 1363071361 : for ( ; i<ly; i++) y[i]=0;
45 1201096642 : return;
46 : }
47 1169053717 : for (i=2; i<ly; i++) y[i]=x[i];
48 : /* lx > ly: round properly */
49 223122205 : if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
50 : }
51 :
52 : GEN
53 62903004 : trunc2nr(GEN x, long n)
54 : {
55 : long ex;
56 62903004 : if (!signe(x)) return gen_0;
57 62460374 : ex = expo(x) + n; if (ex < 0) return gen_0;
58 59883224 : 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 54676931 : mantissa_real(GEN x, long *e)
64 : {
65 54676931 : *e = bit_prec(x)-1-expo(x);
66 54676931 : return mantissa2nr(x, 0);
67 : }
68 :
69 : GEN
70 1167452871 : mului(ulong x, GEN y)
71 : {
72 1167452871 : long s = signe(y);
73 : GEN z;
74 :
75 1167452871 : if (!s || !x) return gen_0;
76 1008960571 : z = muluispec(x, y+2, lgefint(y)-2);
77 1008960571 : setsigne(z,s); return z;
78 : }
79 :
80 : GEN
81 621291120 : mulsi(long x, GEN y)
82 : {
83 621291120 : long s = signe(y);
84 : GEN z;
85 :
86 621291120 : if (!s || !x) return gen_0;
87 372507168 : if (x<0) { s = -s; x = -x; }
88 372507168 : z = muluispec((ulong)x, y+2, lgefint(y)-2);
89 372507168 : setsigne(z,s); return z;
90 : }
91 :
92 : GEN
93 201188881 : mulss(long x, long y)
94 : {
95 : ulong p1;
96 : LOCAL_HIREMAINDER;
97 :
98 201188881 : if (!x || !y) return gen_0;
99 200852370 : if (x<0) {
100 424620 : x = -x;
101 424620 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoi(hiremainder, p1); }
102 305830 : p1 = mulll(x,y); return uutoineg(hiremainder, p1);
103 : } else {
104 200427750 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
105 200384241 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
106 : }
107 : }
108 : GEN
109 742866 : sqrs(long x)
110 : {
111 : ulong p1;
112 : LOCAL_HIREMAINDER;
113 :
114 742866 : if (!x) return gen_0;
115 741095 : if (x<0) x = -x;
116 741095 : p1 = mulll(x,x); return uutoi(hiremainder, p1);
117 : }
118 : GEN
119 4809824724 : muluu(ulong x, ulong y)
120 : {
121 : ulong p1;
122 : LOCAL_HIREMAINDER;
123 :
124 4809824724 : if (!x || !y) return gen_0;
125 4809786528 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
126 : }
127 : GEN
128 648177240 : sqru(ulong x)
129 : {
130 : ulong p1;
131 : LOCAL_HIREMAINDER;
132 :
133 648177240 : if (!x) return gen_0;
134 647779555 : 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 759151313 : 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 759151313 : if (!(x & (x-1))) { z = shiftr(y, expu(x)); setsigne(z, s); return z; }
147 589274825 : lx = lg(y); z = cgetg(lx, t_REAL); e = expo(y);
148 589274825 : y--; garde = mulll(x,y[lx]);
149 2959495091 : for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
150 589274825 : z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
151 589274825 : sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
152 589274825 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
153 589274825 : z[1] = evalsigne(s) | evalexpo(m+e);
154 589274825 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
155 589274825 : return z;
156 : }
157 :
158 : INLINE GEN
159 481871 : mul0r(GEN x)
160 : {
161 481871 : long l = realprec(x), e = expo(x);
162 481871 : e = (l > 0)? e - l: (e < 0? 2*e: 0);
163 481871 : 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 132688715 : mulsr(long x, GEN y)
174 : {
175 : long s;
176 :
177 132688715 : if (!x) return mul0r(y);
178 132688673 : s = signe(y);
179 132688673 : if (!s)
180 : {
181 247427 : if (x < 0) x = -x;
182 247427 : return real_0_bit( expo(y) + expu(x) );
183 : }
184 132441246 : if (x==1) return rcopy(y);
185 123649502 : if (x==-1) return negr(y);
186 120086821 : if (x < 0)
187 42783324 : return mulur_2((ulong)-x, y, -s);
188 : else
189 77303497 : return mulur_2((ulong)x, y, s);
190 : }
191 :
192 : GEN
193 974583016 : mulur(ulong x, GEN y)
194 : {
195 : long s;
196 :
197 974583016 : if (!x) return mul0r(y);
198 974583009 : s = signe(y);
199 974583009 : if (!s) return real_0_bit( expo(y) + expu(x) );
200 966618285 : if (x==1) return rcopy(y);
201 639064492 : return mulur_2(x, y, s);
202 : }
203 :
204 : INLINE void
205 3775528395 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
206 : {
207 : long i;
208 3775528395 : if (hi[2] < 0)
209 : {
210 1740274190 : if (z != hi)
211 1714926505 : for (i=2; i<lz ; i++) z[i] = hi[i];
212 1740274190 : ez++;
213 : }
214 : else
215 : {
216 2035254205 : shift_left(z,hi,2,lz-1, garde, 1);
217 2035254205 : garde <<= 1;
218 : }
219 3775528395 : if (garde & HIGHBIT)
220 : { /* round to nearest */
221 1850011670 : i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
222 1831722317 : if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
223 : }
224 3775528395 : z[1] = evalsigne(sz)|evalexpo(ez);
225 3775528395 : }
226 : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
227 : INLINE void
228 790185232 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
229 : {
230 790185232 : if (z[2] < 0)
231 : { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
232 368396894 : if (garde & HIGHBIT) z[2]++; /* round properly */
233 368396894 : ez++;
234 : }
235 : else
236 : {
237 421788338 : uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
238 421788338 : if (garde & (1UL<<(BITS_IN_LONG-2)))
239 : {
240 172543743 : uel(z,2)++; /* round properly, z2+1 can overflow */
241 172543743 : if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
242 : }
243 : }
244 790185232 : z[1] = evalsigne(sz)|evalexpo(ez);
245 790185232 : }
246 :
247 : /* set z <-- x^2 != 0, floating point multiplication.
248 : * lz = lg(z) = lg(x) */
249 : INLINE void
250 568767212 : sqrz_i(GEN z, GEN x, long lz)
251 : {
252 568767212 : 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 568767212 : if (lz > prec2lg(SQRR_SQRI_LIMIT))
260 : {
261 40293217 : pari_sp av = avma;
262 40293217 : GEN hi = sqrispec_mirror(x+2, lz-2);
263 40293217 : mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
264 40293217 : set_avma(av); return;
265 : }
266 528473995 : if (lz == 3)
267 : {
268 123028563 : garde = mulll(x[2],x[2]);
269 123028563 : z[2] = hiremainder;
270 123028563 : mulrrz_3end(z, 1, ez, garde);
271 123028563 : return;
272 : }
273 :
274 405445432 : lzz = lz-1; p1 = x[lzz];
275 405445432 : if (p1)
276 : {
277 368930740 : (void)mulll(p1,x[3]);
278 368930740 : garde = addmul(p1,x[2]);
279 368930740 : z[lzz] = hiremainder;
280 : }
281 : else
282 : {
283 36514692 : garde = 0;
284 36514692 : z[lzz] = 0;
285 : }
286 1849334918 : for (j=lz-2, x1=x-j; j>=3; j--)
287 : {
288 1443889486 : p1 = x[j]; x1++;
289 1443889486 : if (p1)
290 : {
291 1415968036 : (void)mulll(p1,x1[lz+1]);
292 1415968036 : garde = addll(addmul(p1,x1[lz]), garde);
293 7187970062 : for (i=lzz; i>j; i--)
294 : {
295 5772002026 : hiremainder += overflow;
296 5772002026 : z[i] = addll(addmul(p1,x1[i]), z[i]);
297 : }
298 1415968036 : z[j] = hiremainder+overflow;
299 : }
300 27921450 : else z[j]=0;
301 : }
302 405445432 : p1 = x[2]; x1++;
303 405445432 : garde = addll(mulll(p1,x1[lz]), garde);
304 2254780350 : for (i=lzz; i>2; i--)
305 : {
306 1849334918 : hiremainder += overflow;
307 1849334918 : z[i] = addll(addmul(p1,x1[i]), z[i]);
308 : }
309 405445432 : z[2] = hiremainder+overflow;
310 405445432 : 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 667156669 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
326 : {
327 : ulong garde;
328 : LOCAL_HIREMAINDER;
329 667156669 : if (flag)
330 : {
331 79529319 : (void)mulll(x[2],y[3]);
332 79529319 : garde = addmul(x[2],y[2]);
333 : }
334 : else
335 587627350 : garde = mulll(x[2],y[2]);
336 667156669 : z[2] = hiremainder;
337 667156669 : mulrrz_3end(z, sz, expo(x)+expo(y), garde);
338 667156669 : }
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 3946707031 : 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 3946707031 : if (x == y) { sqrz_i(z,x,lz); return; }
353 3946707031 : if (lz > prec2lg(MULRR_MULII_LIMIT)) { mulrrz_int(z,x,y,lz,flag,sz); return; }
354 3892604873 : if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
355 3225448204 : ez = expo(x) + expo(y);
356 3225448204 : if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
357 3225448204 : lzz=lz-1; p1=x[lzz];
358 3225448204 : if (p1)
359 : {
360 3010019185 : (void)mulll(p1,y[3]);
361 3010019185 : garde = addll(addmul(p1,y[2]), garde);
362 3010019185 : z[lzz] = overflow+hiremainder;
363 : }
364 215429019 : else z[lzz]=0;
365 14166911236 : for (j=lz-2, y1=y-j; j>=3; j--)
366 : {
367 10941463032 : p1 = x[j]; y1++;
368 10941463032 : if (p1)
369 : {
370 10627532876 : (void)mulll(p1,y1[lz+1]);
371 10627532876 : garde = addll(addmul(p1,y1[lz]), garde);
372 64105261188 : for (i=lzz; i>j; i--)
373 : {
374 53477728312 : hiremainder += overflow;
375 53477728312 : z[i] = addll(addmul(p1,y1[i]), z[i]);
376 : }
377 10627532876 : z[j] = hiremainder+overflow;
378 : }
379 313930156 : else z[j]=0;
380 : }
381 3225448204 : p1 = x[2]; y1++;
382 3225448204 : garde = addll(mulll(p1,y1[lz]), garde);
383 17392359440 : for (i=lzz; i>2; i--)
384 : {
385 14166911236 : hiremainder += overflow;
386 14166911236 : z[i] = addll(addmul(p1,y1[i]), z[i]);
387 : }
388 3225448204 : z[2] = hiremainder+overflow;
389 3225448204 : mulrrz_end(z, z, lz, sz, ez, garde);
390 : }
391 :
392 : GEN
393 4197556848 : mulrr(GEN x, GEN y)
394 : {
395 : long flag, ly, lz, sx, sy;
396 : GEN z;
397 :
398 4197556848 : if (x == y) return sqrr(x);
399 4195937069 : sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
400 3920309101 : sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
401 3882600873 : if (sy < 0) sx = -sx;
402 3882600873 : lz = lg(x);
403 3882600873 : ly = lg(y);
404 3882600873 : if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
405 3882600873 : z = cgetg(lz, t_REAL);
406 3882600873 : mulrrz_i(z, x,y, lz,flag, sx);
407 3882600873 : return z;
408 : }
409 :
410 : GEN
411 602794940 : sqrr(GEN x)
412 : {
413 602794940 : long lz, sx = signe(x);
414 : GEN z;
415 :
416 602794940 : if (!sx) return real_0_bit(2*expo(x));
417 568767212 : lz = lg(x); z = cgetg(lz, t_REAL);
418 568767212 : sqrz_i(z, x, lz);
419 568767212 : return z;
420 : }
421 :
422 : GEN
423 1028139018 : mulir(GEN x, GEN y)
424 : {
425 1028139018 : long sx = signe(x), sy;
426 1028139018 : if (!sx) return mul0r(y);
427 1027657196 : if (lgefint(x) == 3) {
428 912432879 : GEN z = mulur(uel(x,2), y);
429 912432879 : if (sx < 0) togglesign(z);
430 912432879 : return z;
431 : }
432 115224317 : sy = signe(y);
433 115224317 : if (!sy) return real_0_bit(expi(x) + expo(y));
434 114345542 : if (sy < 0) sx = -sx;
435 : {
436 114345542 : long lz = lg(y), lx = lgefint(x);
437 114345542 : GEN hi, z = cgetg(lz, t_REAL);
438 114345542 : pari_sp av = avma;
439 114345542 : 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 50239384 : x = itor(x, lg2prec(lx));
443 50239384 : hi = muliispec_mirror(y+2, x+2, lz-2, lx-2);
444 50239384 : mulrrz_end(z, hi, lz, sx, expo(x)+expo(y), hi[lz]);
445 : }
446 : else /* dubious: complete x with 0s and call mulrr */
447 64106158 : mulrrz_i(z, itor(x, lg2prec(lz)), y, lz, 0, sx);
448 114345542 : 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 92949740 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
455 : {
456 92949740 : long lx = lgefint(x), ly;
457 : pari_sp av;
458 : GEN t;
459 92949740 : if (lx == 2) return mulii(z,y);
460 90128828 : ly = lgefint(y);
461 90128828 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
462 89511175 : av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
463 89511175 : t = mulii(z, y);
464 89511175 : set_avma(av); return addii(t,x);
465 : }
466 : /* x + y*z, lgefint(z) == 3 */
467 : static GEN
468 482952673 : addmulii_lg3(GEN x, GEN y, GEN z)
469 : {
470 482952673 : long s = signe(z), lx, ly;
471 482952673 : ulong w = z[2];
472 : pari_sp av;
473 : GEN t;
474 482952673 : if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
475 396872541 : lx = lgefint(x);
476 396872541 : ly = lgefint(y);
477 396872541 : if (lx == 2)
478 : { /* x = 0 */
479 73369677 : if (ly == 2) return gen_0;
480 40868836 : t = muluispec(w, y+2, ly-2);
481 40868836 : if (signe(y) < 0) s = -s;
482 40868836 : setsigne(t, s); return t;
483 : }
484 323502864 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
485 282685236 : av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
486 282685236 : t = muluispec(w, y+2, ly-2);
487 282685236 : if (signe(y) < 0) s = -s;
488 282685236 : setsigne(t, s);
489 282685236 : set_avma(av); return addii(x,t);
490 : }
491 : /* x + y*z */
492 : GEN
493 169826513 : addmulii(GEN x, GEN y, GEN z)
494 : {
495 169826513 : long lz = lgefint(z);
496 169826513 : switch(lz)
497 : {
498 759297 : case 2: return icopy(x); /* z = 0, wasteful copy */
499 128887740 : case 3: return addmulii_lg3(x, y, z);
500 40179476 : 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 1631435363 : addmulii_inplace(GEN x, GEN y, GEN z)
506 : {
507 : long lz;
508 1631435363 : if (lgefint(y) == 2) return x;
509 406854380 : lz = lgefint(z);
510 406854380 : switch(lz)
511 : {
512 19183 : case 2: return x;
513 354064933 : case 3: return addmulii_lg3(x, y, z);
514 52770264 : 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 2384407099 : 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 2384407099 : if (!z) return -1;
528 : #ifdef LONG_IS_64BIT
529 2143775046 : if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
530 : #endif
531 2384407085 : z |= ~z + 1;
532 2384407085 : z += z << 4;
533 2384407085 : z += z << 6;
534 2384407085 : z ^= z << 16; /* or z -= z<<16 */
535 : #ifdef LONG_IS_64BIT
536 2143775046 : 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 3278537 : divir(GEN x, GEN y)
558 : {
559 : GEN z;
560 3278537 : long ly = lg(y), lx = lgefint(x);
561 : pari_sp av;
562 :
563 3278537 : if (ly == 2) pari_err_INV("divir",y);
564 3278537 : if (lx == 2) return div0r(y);
565 3278537 : if (lx == 3) {
566 1913525 : z = divur(x[2], y);
567 1913525 : if (signe(x) < 0) togglesign(z);
568 1913525 : return z;
569 : }
570 1365012 : z = cgetg(ly, t_REAL); av = avma;
571 1365012 : affrr(divrr(itor(x, lg2prec(ly+1)), y), z);
572 1365012 : set_avma(av); return z;
573 : }
574 :
575 : GEN
576 2456239 : divur(ulong x, GEN y)
577 : {
578 : pari_sp av;
579 2456239 : long p = realprec(y);
580 : GEN z;
581 :
582 2456239 : if (p == 0) pari_err_INV("divur",y);
583 2456239 : if (!x) return div0r(y);
584 2456239 : if (x == 1) return invr(y);
585 2446397 : if (!(x & (x-1))) /* power of 2 */
586 : {
587 913205 : z = invr(y);
588 913205 : shiftr_inplace(z, expu(x)); return z;
589 : }
590 1533192 : 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 1533188 : z = cgetr(p); av = avma;
596 1533188 : affrr(divrr(utor(x, p + BITS_IN_LONG), y), z);
597 1533188 : 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 74743843 : invr_basecase(GEN y)
610 : {
611 74743843 : long p = realprec(y);
612 74743843 : GEN z = cgetr(p);
613 74743843 : pari_sp av = avma;
614 74743843 : affrr(divrr(real_1(p + BITS_IN_LONG), y), z);
615 74743843 : set_avma(av); return z;
616 : }
617 : /* returns 1/b, Newton iteration */
618 : GEN
619 74743843 : invr(GEN b)
620 : {
621 74743843 : const long s = 6;
622 74743843 : long i, p, l = lg(b);
623 : GEN x, a;
624 : ulong mask;
625 :
626 74743843 : if (l <= maxss(prec2lg(INVNEWTON_LIMIT), (1L<<s) + 2)) {
627 74739489 : if (l == 2) pari_err_INV("invr",b);
628 74739489 : 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 3154928971 : modii(GEN x, GEN y)
652 : {
653 3154928971 : switch(signe(x))
654 : {
655 498469264 : case 0: return gen_0;
656 2094802767 : case 1: return remii(x,y);
657 561656940 : default:
658 : {
659 561656940 : pari_sp av = avma;
660 561656940 : (void)new_chunk(lgefint(y));
661 561656940 : x = remii(x,y); set_avma(av);
662 561656940 : if (x==gen_0) return x;
663 519716789 : return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
664 : }
665 : }
666 : }
667 :
668 : GEN
669 20390035 : divrs(GEN x, long y)
670 : {
671 : GEN z;
672 20390035 : if (y < 0)
673 : {
674 5905434 : z = divru(x, (ulong)-y);
675 5905434 : togglesign(z);
676 : }
677 : else
678 14484601 : z = divru(x, (ulong)y);
679 20390035 : return z;
680 : }
681 :
682 : GEN
683 1156631929 : divru(GEN x, ulong y)
684 : {
685 1156631929 : long i, lx, sh, e, s = signe(x);
686 : ulong garde;
687 : GEN z;
688 : LOCAL_HIREMAINDER;
689 :
690 1156631929 : if (!y) pari_err_INV("divru",gen_0);
691 1156631929 : if (!s) return real_0_bit(expo(x) - expu(y));
692 1155937290 : if (!(y & (y-1))) /* power of 2 */
693 : {
694 98861864 : if (y == 1) return rcopy(x);
695 96122364 : return shiftr(x, -expu(y));
696 : }
697 1057075426 : e = expo(x);
698 1057075426 : lx = lg(x);
699 1057075426 : z = cgetg(lx, t_REAL);
700 1057075426 : if (lx == 3)
701 : {
702 175535754 : if (y <= uel(x,2))
703 : {
704 175535484 : hiremainder = 0;
705 175535484 : z[2] = divll(x[2],y);
706 : /* we may have hiremainder != 0 ==> garde */
707 175535484 : 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 881539672 : ulong yp = get_Fl_red(y);
720 881539672 : if (y <= uel(x,2))
721 : {
722 881535906 : hiremainder = 0;
723 7497999035 : for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
724 : /* we may have hiremainder != 0 ==> garde */
725 881535906 : garde = divll_pre(0,y,yp);
726 : }
727 : else
728 : {
729 3766 : long l = lx-1;
730 3766 : hiremainder = x[2];
731 64961 : for (i=2; i<l; i++) z[i] = divll_pre(x[i+1],y,yp);
732 3766 : z[i] = divll_pre(0,y,yp);
733 3766 : garde = hiremainder;
734 3766 : e -= BITS_IN_LONG;
735 : }
736 : }
737 1057075426 : sh=bfffo(z[2]); /* z[2] != 0 */
738 1057075426 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
739 1057075426 : z[1] = evalsigne(s) | evalexpo(e-sh);
740 1057075426 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
741 1057075426 : return z;
742 : }
743 :
744 : GEN
745 142365386 : truedvmdii(GEN x, GEN y, GEN *z)
746 : {
747 : pari_sp av;
748 : GEN r, q;
749 142365386 : if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
750 5488405 : if (z == ONLY_REM) return modii(x,y);
751 :
752 5488405 : av = avma;
753 5488405 : q = dvmdii(x,y,&r); /* assume that r is last on stack */
754 5488405 : switch(signe(r))
755 : {
756 711814 : case 0:
757 711814 : if (z) *z = gen_0;
758 711814 : return q;
759 3735430 : case 1:
760 3735430 : if (z) *z = r; else cgiv(r);
761 3735430 : return q;
762 1041161 : case -1: break;
763 : }
764 1041161 : q = addis(q, -signe(y));
765 1041161 : if (!z) return gc_INT(av, q);
766 :
767 865063 : *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
768 865063 : return gc_all_unsafe(av,(pari_sp)r,2,&q,z);
769 : }
770 : GEN
771 139105733 : truedvmdis(GEN x, long y, GEN *z)
772 : {
773 139105733 : pari_sp av = avma;
774 : long r;
775 : GEN q;
776 :
777 139105733 : if (z == ONLY_REM) return modis(x, y);
778 139105733 : q = divis_rem(x,y,&r);
779 :
780 139105733 : if (r >= 0)
781 : {
782 126211150 : if (z) *z = utoi(r);
783 126211150 : return q;
784 : }
785 12894583 : q = gc_INT(av, addis(q, (y < 0)? 1: -1));
786 12894583 : if (z) *z = utoi(r + labs(y));
787 12894583 : 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 109733767 : int2n(long n) {
809 : long i, m, l;
810 : GEN z;
811 109733767 : if (n < 0) return gen_0;
812 109714755 : if (n == 0) return gen_1;
813 :
814 105960142 : l = dvmdsBIL(n, &m) + 3;
815 105960142 : z = cgetipos(l);
816 568805546 : for (i = 2; i < l; i++) z[i] = 0;
817 105960142 : *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 25996753 : int2u(ulong n) {
823 : ulong i, m, l;
824 : GEN z;
825 25996753 : if (n == 0) return gen_1;
826 :
827 25996753 : l = dvmduBIL(n, &m) + 3;
828 25996753 : z = cgetipos(l);
829 52386955 : for (i = 2; i < l; i++) z[i] = 0;
830 25996753 : *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 2323477028 : shifti(GEN x, long n)
849 : {
850 2323477028 : long s = signe(x);
851 : GEN y;
852 :
853 2323477028 : if(s == 0) return gen_0;
854 2015696809 : y = shiftispec(x + 2, lgefint(x) - 2, n);
855 2015696808 : if (signe(y)) setsigne(y, s);
856 2015696808 : return y;
857 : }
858 :
859 : /* actual operations will take place on a+2 and b+2: we strip the codewords */
860 : GEN
861 21995601330 : mulii(GEN a,GEN b)
862 : {
863 : long sa,sb;
864 : GEN z;
865 :
866 21995601330 : sa=signe(a); if (!sa) return gen_0;
867 14892492511 : sb=signe(b); if (!sb) return gen_0;
868 10914135515 : if (sb<0) sa = -sa;
869 10914135515 : z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
870 10914135514 : setsigne(z,sa); return z;
871 : }
872 :
873 : GEN
874 1963819188 : 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 258620549 : usqrt(ulong a)
880 : {
881 258620549 : ulong x = (ulong)sqrt((double)a);
882 : #ifdef LONG_IS_64BIT
883 224862722 : if (x > LOWMASK || x*x > a) x--;
884 : #endif
885 258620549 : return x;
886 : }
887 :
888 : /********************************************************************/
889 : /** **/
890 : /** EXPONENT / CONVERSION t_REAL --> double **/
891 : /** **/
892 : /********************************************************************/
893 :
894 : #ifdef LONG_IS_64BIT
895 : long
896 23762009 : dblexpo(double x)
897 : {
898 : union { double f; ulong i; } fi;
899 23762009 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
900 23762009 : const int exp_mid = 0x3ff;/* exponent bias */
901 :
902 23762009 : if (x==0.) return -exp_mid;
903 23761997 : fi.f = x;
904 23761997 : 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 11706061 : dbltor(double x)
920 : {
921 : GEN z;
922 : long e;
923 : union { double f; ulong i; } fi;
924 11706061 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
925 11706061 : const int exp_mid = 0x3ff;/* exponent bias */
926 11706061 : const int expo_len = 11; /* number of bits of exponent */
927 :
928 11706061 : if (x==0.) return real_0_bit(-exp_mid);
929 11699739 : fi.f = x; z = cgetr(DEFAULTPREC);
930 : {
931 11699739 : const ulong a = fi.i;
932 : ulong A;
933 11699739 : e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
934 11699739 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
935 11699739 : A = a << expo_len;
936 11699739 : 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 11699739 : z[2] = HIGHBIT | A;
944 11699739 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
945 : }
946 11699739 : return z;
947 : }
948 :
949 : double
950 268113193 : rtodbl(GEN x)
951 : {
952 268113193 : long ex,s=signe(x);
953 : ulong a;
954 : union { double f; ulong i; } fi;
955 268113193 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
956 268113193 : const int exp_mid = 0x3ff;/* exponent bias */
957 268113193 : const int expo_len = 11; /* number of bits of exponent */
958 :
959 268113193 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
960 :
961 : /* start by rounding to closest */
962 236142242 : a = (x[2] & (HIGHBIT-1)) + 0x400;
963 236142242 : if (a & HIGHBIT) { ex++; a=0; }
964 236142242 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
965 236142242 : fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
966 236142242 : if (s<0) fi.i |= HIGHBIT;
967 236142242 : 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 :
|
