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 61368762 : invmod2BIL(ulong b)
19 : {
20 : static int tab[] = { 0, 0, 0, 8, 0, 8, 0, 0 };
21 61368762 : 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 61368762 : x = x*(2-b*x); /* one more pass necessary */
26 : #endif
27 61368762 : x = x*(2-b*x);
28 61368762 : x = x*(2-b*x); return x*(2-b*x);
29 : }
30 :
31 : void
32 1424436239 : affrr(GEN x, GEN y)
33 : {
34 1424436239 : long i, lx, ly = lg(y);
35 1424436239 : if (!signe(x))
36 : {
37 211493 : y[1] = evalexpo(minss(expo(x), -bit_accuracy(ly)));
38 211493 : return;
39 : }
40 1424224746 : y[1] = x[1]; lx = lg(x);
41 1424224746 : if (lx <= ly)
42 : {
43 7528635134 : for (i=2; i<lx; i++) y[i]=x[i];
44 1363073287 : for ( ; i<ly; i++) y[i]=0;
45 1201098698 : return;
46 : }
47 1169067088 : for (i=2; i<ly; i++) y[i]=x[i];
48 : /* lx > ly: round properly */
49 223126048 : if (x[ly] & HIGHBIT) roundr_up_ip(y, ly);
50 : }
51 :
52 : GEN
53 62903940 : trunc2nr(GEN x, long n)
54 : {
55 : long ex;
56 62903940 : if (!signe(x)) return gen_0;
57 62461178 : ex = expo(x) + n; if (ex < 0) return gen_0;
58 59883911 : 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 54677500 : mantissa_real(GEN x, long *e)
64 : {
65 54677500 : *e = bit_prec(x)-1-expo(x);
66 54677500 : return mantissa2nr(x, 0);
67 : }
68 :
69 : GEN
70 1168026446 : mului(ulong x, GEN y)
71 : {
72 1168026446 : long s = signe(y);
73 : GEN z;
74 :
75 1168026446 : if (!s || !x) return gen_0;
76 1009534263 : z = muluispec(x, y+2, lgefint(y)-2);
77 1009534262 : setsigne(z,s); return z;
78 : }
79 :
80 : GEN
81 621294299 : mulsi(long x, GEN y)
82 : {
83 621294299 : long s = signe(y);
84 : GEN z;
85 :
86 621294299 : if (!s || !x) return gen_0;
87 372525456 : if (x<0) { s = -s; x = -x; }
88 372525456 : z = muluispec((ulong)x, y+2, lgefint(y)-2);
89 372525456 : setsigne(z,s); return z;
90 : }
91 :
92 : GEN
93 201188895 : mulss(long x, long y)
94 : {
95 : ulong p1;
96 : LOCAL_HIREMAINDER;
97 :
98 201188895 : if (!x || !y) return gen_0;
99 200852384 : 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 200427764 : if (y<0) { y = -y; p1 = mulll(x,y); return uutoineg(hiremainder, p1); }
105 200384255 : 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 4810117312 : muluu(ulong x, ulong y)
120 : {
121 : ulong p1;
122 : LOCAL_HIREMAINDER;
123 :
124 4810117312 : if (!x || !y) return gen_0;
125 4810079116 : p1 = mulll(x,y); return uutoi(hiremainder, p1);
126 : }
127 : GEN
128 648196261 : sqru(ulong x)
129 : {
130 : ulong p1;
131 : LOCAL_HIREMAINDER;
132 :
133 648196261 : if (!x) return gen_0;
134 647798799 : 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 759159789 : 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 759159789 : if (!(x & (x-1))) { z = shiftr(y, expu(x)); setsigne(z, s); return z; }
147 589285301 : lx = lg(y); z = cgetg(lx, t_REAL); e = expo(y);
148 589285301 : y--; garde = mulll(x,y[lx]);
149 2959578384 : for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
150 589285301 : z[2]=hiremainder; /* != 0 since y normalized and |x| > 1 */
151 589285301 : sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
152 589285301 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
153 589285301 : z[1] = evalsigne(s) | evalexpo(m+e);
154 589285301 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
155 589285301 : return z;
156 : }
157 :
158 : INLINE GEN
159 481881 : mul0r(GEN x)
160 : {
161 481881 : long l = realprec(x), e = expo(x);
162 481881 : e = (l > 0)? e - l: (e < 0? 2*e: 0);
163 481881 : 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 132688777 : mulsr(long x, GEN y)
174 : {
175 : long s;
176 :
177 132688777 : if (!x) return mul0r(y);
178 132688735 : s = signe(y);
179 132688735 : if (!s)
180 : {
181 247247 : if (x < 0) x = -x;
182 247247 : return real_0_bit( expo(y) + expu(x) );
183 : }
184 132441488 : if (x==1) return rcopy(y);
185 123649947 : if (x==-1) return negr(y);
186 120087267 : if (x < 0)
187 42783339 : return mulur_2((ulong)-x, y, -s);
188 : else
189 77303928 : return mulur_2((ulong)x, y, s);
190 : }
191 :
192 : GEN
193 974587249 : mulur(ulong x, GEN y)
194 : {
195 : long s;
196 :
197 974587249 : if (!x) return mul0r(y);
198 974587242 : s = signe(y);
199 974587242 : if (!s) return real_0_bit( expo(y) + expu(x) );
200 966623144 : if (x==1) return rcopy(y);
201 639072522 : return mulur_2(x, y, s);
202 : }
203 :
204 : INLINE void
205 3775417690 : mulrrz_end(GEN z, GEN hi, long lz, long sz, long ez, ulong garde)
206 : {
207 : long i;
208 3775417690 : if (hi[2] < 0)
209 : {
210 1740324721 : if (z != hi)
211 1714925936 : for (i=2; i<lz ; i++) z[i] = hi[i];
212 1740324721 : ez++;
213 : }
214 : else
215 : {
216 2035092969 : shift_left(z,hi,2,lz-1, garde, 1);
217 2035092969 : garde <<= 1;
218 : }
219 3775417690 : if (garde & HIGHBIT)
220 : { /* round to nearest */
221 1850059951 : i = lz; do ((ulong*)z)[--i]++; while (i>1 && z[i]==0);
222 1831771054 : if (i == 1) { z[2] = (long)HIGHBIT; ez++; }
223 : }
224 3775417690 : z[1] = evalsigne(sz)|evalexpo(ez);
225 3775417690 : }
226 : /* mulrrz_end for lz = 3, minor simplifications. z[2]=hiremainder from mulll */
227 : INLINE void
228 790144359 : mulrrz_3end(GEN z, long sz, long ez, ulong garde)
229 : {
230 790144359 : if (z[2] < 0)
231 : { /* z2 < (2^BIL-1)^2 / 2^BIL, hence z2+1 != 0 */
232 368400084 : if (garde & HIGHBIT) z[2]++; /* round properly */
233 368400084 : ez++;
234 : }
235 : else
236 : {
237 421744275 : uel(z,2) = (uel(z,2)<<1) | (garde>>(BITS_IN_LONG-1));
238 421744275 : if (garde & (1UL<<(BITS_IN_LONG-2)))
239 : {
240 172538979 : uel(z,2)++; /* round properly, z2+1 can overflow */
241 172538979 : if (!uel(z,2)) { uel(z,2) = HIGHBIT; ez++; }
242 : }
243 : }
244 790144359 : z[1] = evalsigne(sz)|evalexpo(ez);
245 790144359 : }
246 :
247 : /* set z <-- x^2 != 0, floating point multiplication.
248 : * lz = lg(z) = lg(x) */
249 : INLINE void
250 568775532 : sqrz_i(GEN z, GEN x, long lz)
251 : {
252 568775532 : 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 568775532 : if (lz > prec2lg(SQRR_SQRI_LIMIT))
260 : {
261 40293313 : pari_sp av = avma;
262 40293313 : GEN hi = sqrispec_mirror(x+2, lz-2);
263 40293313 : mulrrz_end(z, hi, lz, 1, ez, hi[lz]);
264 40293313 : set_avma(av); return;
265 : }
266 528482219 : if (lz == 3)
267 : {
268 123035282 : garde = mulll(x[2],x[2]);
269 123035282 : z[2] = hiremainder;
270 123035282 : mulrrz_3end(z, 1, ez, garde);
271 123035282 : return;
272 : }
273 :
274 405446937 : lzz = lz-1; p1 = x[lzz];
275 405446937 : if (p1)
276 : {
277 368932166 : (void)mulll(p1,x[3]);
278 368932166 : garde = addmul(p1,x[2]);
279 368932166 : z[lzz] = hiremainder;
280 : }
281 : else
282 : {
283 36514771 : garde = 0;
284 36514771 : z[lzz] = 0;
285 : }
286 1849335337 : for (j=lz-2, x1=x-j; j>=3; j--)
287 : {
288 1443888400 : p1 = x[j]; x1++;
289 1443888400 : if (p1)
290 : {
291 1415967371 : (void)mulll(p1,x1[lz+1]);
292 1415967371 : garde = addll(addmul(p1,x1[lz]), garde);
293 7187963202 : for (i=lzz; i>j; i--)
294 : {
295 5771995831 : hiremainder += overflow;
296 5771995831 : z[i] = addll(addmul(p1,x1[i]), z[i]);
297 : }
298 1415967371 : z[j] = hiremainder+overflow;
299 : }
300 27921029 : else z[j]=0;
301 : }
302 405446937 : p1 = x[2]; x1++;
303 405446937 : garde = addll(mulll(p1,x1[lz]), garde);
304 2254782274 : for (i=lzz; i>2; i--)
305 : {
306 1849335337 : hiremainder += overflow;
307 1849335337 : z[i] = addll(addmul(p1,x1[i]), z[i]);
308 : }
309 405446937 : z[2] = hiremainder+overflow;
310 405446937 : 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 54102156 : mulrrz_int(GEN z, GEN x, GEN y, long lz, long flag, long sz)
316 : {
317 54102156 : pari_sp av = avma;
318 54102156 : GEN hi = muliispec_mirror(y+2, x+2, lz+flag-2, lz-2);
319 54102156 : mulrrz_end(z, hi, lz, sz, expo(x)+expo(y), hi[lz]);
320 54102156 : set_avma(av);
321 54102156 : }
322 :
323 : /* lz = 3 */
324 : INLINE void
325 667109077 : mulrrz_3(GEN z, GEN x, GEN y, long flag, long sz)
326 : {
327 : ulong garde;
328 : LOCAL_HIREMAINDER;
329 667109077 : if (flag)
330 : {
331 79497588 : (void)mulll(x[2],y[3]);
332 79497588 : garde = addmul(x[2],y[2]);
333 : }
334 : else
335 587611489 : garde = mulll(x[2],y[2]);
336 667109077 : z[2] = hiremainder;
337 667109077 : mulrrz_3end(z, sz, expo(x)+expo(y), garde);
338 667109077 : }
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 3946547133 : 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 3946547133 : if (x == y) { sqrz_i(z,x,lz); return; }
353 3946547133 : if (lz > prec2lg(MULRR_MULII_LIMIT)) { mulrrz_int(z,x,y,lz,flag,sz); return; }
354 3892444977 : if (lz == 3) { mulrrz_3(z,x,y,flag,sz); return; }
355 3225335900 : ez = expo(x) + expo(y);
356 3225335900 : if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde = 0;
357 3225335900 : lzz=lz-1; p1=x[lzz];
358 3225335900 : if (p1)
359 : {
360 3009902320 : (void)mulll(p1,y[3]);
361 3009902320 : garde = addll(addmul(p1,y[2]), garde);
362 3009902320 : z[lzz] = overflow+hiremainder;
363 : }
364 215433580 : else z[lzz]=0;
365 14166590152 : for (j=lz-2, y1=y-j; j>=3; j--)
366 : {
367 10941254252 : p1 = x[j]; y1++;
368 10941254252 : if (p1)
369 : {
370 10627329786 : (void)mulll(p1,y1[lz+1]);
371 10627329786 : garde = addll(addmul(p1,y1[lz]), garde);
372 64104783108 : for (i=lzz; i>j; i--)
373 : {
374 53477453322 : hiremainder += overflow;
375 53477453322 : z[i] = addll(addmul(p1,y1[i]), z[i]);
376 : }
377 10627329786 : z[j] = hiremainder+overflow;
378 : }
379 313924466 : else z[j]=0;
380 : }
381 3225335900 : p1 = x[2]; y1++;
382 3225335900 : garde = addll(mulll(p1,y1[lz]), garde);
383 17391926052 : for (i=lzz; i>2; i--)
384 : {
385 14166590152 : hiremainder += overflow;
386 14166590152 : z[i] = addll(addmul(p1,y1[i]), z[i]);
387 : }
388 3225335900 : z[2] = hiremainder+overflow;
389 3225335900 : mulrrz_end(z, z, lz, sz, ez, garde);
390 : }
391 :
392 : GEN
393 4197059676 : mulrr(GEN x, GEN y)
394 : {
395 : long flag, ly, lz, sx, sy;
396 : GEN z;
397 :
398 4197059676 : if (x == y) return sqrr(x);
399 4195439905 : sx = signe(x); if (!sx) return real_0_bit(expo(x) + expo(y));
400 3919887071 : sy = signe(y); if (!sy) return real_0_bit(expo(x) + expo(y));
401 3882437384 : if (sy < 0) sx = -sx;
402 3882437384 : lz = lg(x);
403 3882437384 : ly = lg(y);
404 3882437384 : if (lz > ly) { lz = ly; swap(x, y); flag = 1; } else flag = (lz != ly);
405 3882437384 : z = cgetg(lz, t_REAL);
406 3882437384 : mulrrz_i(z, x,y, lz,flag, sx);
407 3882437384 : return z;
408 : }
409 :
410 : GEN
411 602801893 : sqrr(GEN x)
412 : {
413 602801893 : long lz, sx = signe(x);
414 : GEN z;
415 :
416 602801893 : if (!sx) return real_0_bit(2*expo(x));
417 568775532 : lz = lg(x); z = cgetg(lz, t_REAL);
418 568775532 : sqrz_i(z, x, lz);
419 568775532 : return z;
420 : }
421 :
422 : GEN
423 1028146591 : mulir(GEN x, GEN y)
424 : {
425 1028146591 : long sx = signe(x), sy;
426 1028146591 : if (!sx) return mul0r(y);
427 1027664759 : if (lgefint(x) == 3) {
428 912436845 : GEN z = mulur(uel(x,2), y);
429 912436845 : if (sx < 0) togglesign(z);
430 912436845 : return z;
431 : }
432 115227914 : sy = signe(y);
433 115227914 : if (!sy) return real_0_bit(expi(x) + expo(y));
434 114349133 : if (sy < 0) sx = -sx;
435 : {
436 114349133 : long lz = lg(y), lx = lgefint(x);
437 114349133 : GEN hi, z = cgetg(lz, t_REAL);
438 114349133 : pari_sp av = avma;
439 114349133 : 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 64109749 : mulrrz_i(z, itor(x, lg2prec(lz)), y, lz, 0, sx);
448 114349133 : 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 92960915 : addmulii_gen(GEN x, GEN y, GEN z, long lz)
455 : {
456 92960915 : long lx = lgefint(x), ly;
457 : pari_sp av;
458 : GEN t;
459 92960915 : if (lx == 2) return mulii(z,y);
460 90139175 : ly = lgefint(y);
461 90139175 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
462 89521124 : av = avma; (void)new_chunk(lx+ly+lz); /*HACK*/
463 89521124 : t = mulii(z, y);
464 89521124 : set_avma(av); return addii(t,x);
465 : }
466 : /* x + y*z, lgefint(z) == 3 */
467 : static GEN
468 482972290 : addmulii_lg3(GEN x, GEN y, GEN z)
469 : {
470 482972290 : long s = signe(z), lx, ly;
471 482972290 : ulong w = z[2];
472 : pari_sp av;
473 : GEN t;
474 482972290 : if (w == 1) return (s > 0)? addii(x,y): subii(x,y); /* z = +- 1 */
475 396894271 : lx = lgefint(x);
476 396894271 : ly = lgefint(y);
477 396894271 : if (lx == 2)
478 : { /* x = 0 */
479 73372860 : if (ly == 2) return gen_0;
480 40869433 : t = muluispec(w, y+2, ly-2);
481 40869433 : if (signe(y) < 0) s = -s;
482 40869433 : setsigne(t, s); return t;
483 : }
484 323521411 : if (ly == 2) return icopy(x); /* y = 0, wasteful copy */
485 282696800 : av = avma; (void)new_chunk(1+lx+ly);/*HACK*/
486 282696800 : t = muluispec(w, y+2, ly-2);
487 282696800 : if (signe(y) < 0) s = -s;
488 282696800 : setsigne(t, s);
489 282696800 : set_avma(av); return addii(x,t);
490 : }
491 : /* x + y*z */
492 : GEN
493 169849648 : addmulii(GEN x, GEN y, GEN z)
494 : {
495 169849648 : long lz = lgefint(z);
496 169849648 : switch(lz)
497 : {
498 759297 : case 2: return icopy(x); /* z = 0, wasteful copy */
499 128908862 : case 3: return addmulii_lg3(x, y, z);
500 40181489 : 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 1631473218 : addmulii_inplace(GEN x, GEN y, GEN z)
506 : {
507 : long lz;
508 1631473218 : if (lgefint(y) == 2) return x;
509 406862035 : lz = lgefint(z);
510 406862035 : switch(lz)
511 : {
512 19181 : case 2: return x;
513 354063428 : case 3: return addmulii_lg3(x, y, z);
514 52779426 : 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 2384242349 : 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 2384242349 : if (!z) return -1;
528 : #ifdef LONG_IS_64BIT
529 2143610294 : if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
530 : #endif
531 2384242335 : z |= ~z + 1;
532 2384242335 : z += z << 4;
533 2384242335 : z += z << 6;
534 2384242335 : z ^= z << 16; /* or z -= z<<16 */
535 : #ifdef LONG_IS_64BIT
536 2143610294 : return s + tab[(z&0xffffffff)>>26];
537 : #else
538 240632041 : 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 2456238 : divur(ulong x, GEN y)
577 : {
578 : pari_sp av;
579 2456238 : long p = realprec(y);
580 : GEN z;
581 :
582 2456238 : if (p == 0) pari_err_INV("divur",y);
583 2456238 : if (!x) return div0r(y);
584 2456238 : if (x == 1) return invr(y);
585 2446396 : if (!(x & (x-1))) /* power of 2 */
586 : {
587 913204 : z = invr(y);
588 913204 : 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 74746357 : invr_basecase(GEN y)
610 : {
611 74746357 : long p = realprec(y);
612 74746357 : GEN z = cgetr(p);
613 74746357 : pari_sp av = avma;
614 74746357 : affrr(divrr(real_1(p + BITS_IN_LONG), y), z);
615 74746357 : set_avma(av); return z;
616 : }
617 : /* returns 1/b, Newton iteration */
618 : GEN
619 74746357 : invr(GEN b)
620 : {
621 74746357 : const long s = 6;
622 74746357 : long i, p, l = lg(b);
623 : GEN x, a;
624 : ulong mask;
625 :
626 74746357 : if (l <= maxss(prec2lg(INVNEWTON_LIMIT), (1L<<s) + 2)) {
627 74742003 : if (l == 2) pari_err_INV("invr",b);
628 74742003 : 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 3154711214 : modii(GEN x, GEN y)
652 : {
653 3154711214 : switch(signe(x))
654 : {
655 498554740 : case 0: return gen_0;
656 2094456063 : case 1: return remii(x,y);
657 561700411 : default:
658 : {
659 561700411 : pari_sp av = avma;
660 561700411 : (void)new_chunk(lgefint(y));
661 561700411 : x = remii(x,y); set_avma(av);
662 561700411 : if (x==gen_0) return x;
663 519760449 : return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
664 : }
665 : }
666 : }
667 :
668 : GEN
669 20390042 : divrs(GEN x, long y)
670 : {
671 : GEN z;
672 20390042 : if (y < 0)
673 : {
674 5905441 : z = divru(x, (ulong)-y);
675 5905441 : togglesign(z);
676 : }
677 : else
678 14484601 : z = divru(x, (ulong)y);
679 20390042 : return z;
680 : }
681 :
682 : GEN
683 1156629128 : divru(GEN x, ulong y)
684 : {
685 1156629128 : long i, lx, sh, e, s = signe(x);
686 : ulong garde;
687 : GEN z;
688 : LOCAL_HIREMAINDER;
689 :
690 1156629128 : if (!y) pari_err_INV("divru",gen_0);
691 1156629128 : if (!s) return real_0_bit(expo(x) - expu(y));
692 1155934454 : if (!(y & (y-1))) /* power of 2 */
693 : {
694 98861342 : if (y == 1) return rcopy(x);
695 96121896 : return shiftr(x, -expu(y));
696 : }
697 1057073112 : e = expo(x);
698 1057073112 : lx = lg(x);
699 1057073112 : z = cgetg(lx, t_REAL);
700 1057073112 : if (lx == 3)
701 : {
702 175536379 : if (y <= uel(x,2))
703 : {
704 175536109 : hiremainder = 0;
705 175536109 : z[2] = divll(x[2],y);
706 : /* we may have hiremainder != 0 ==> garde */
707 175536109 : 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 881536733 : ulong yp = get_Fl_red(y);
720 881536733 : if (y <= uel(x,2))
721 : {
722 881532967 : hiremainder = 0;
723 7497985593 : for (i=2; i<lx; i++) z[i] = divll_pre(x[i],y,yp);
724 : /* we may have hiremainder != 0 ==> garde */
725 881532967 : 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 1057073112 : sh=bfffo(z[2]); /* z[2] != 0 */
738 1057073112 : if (sh) shift_left(z,z, 2,lx-1, garde,sh);
739 1057073112 : z[1] = evalsigne(s) | evalexpo(e-sh);
740 1057073112 : if ((garde << sh) & HIGHBIT) roundr_up_ip(z, lx);
741 1057073112 : return z;
742 : }
743 :
744 : GEN
745 142350569 : truedvmdii(GEN x, GEN y, GEN *z)
746 : {
747 : pari_sp av;
748 : GEN r, q;
749 142350569 : if (!is_bigint(y)) return truedvmdis(x, itos(y), z);
750 5488260 : if (z == ONLY_REM) return modii(x,y);
751 :
752 5488260 : av = avma;
753 5488260 : q = dvmdii(x,y,&r); /* assume that r is last on stack */
754 5488260 : switch(signe(r))
755 : {
756 711836 : case 0:
757 711836 : if (z) *z = gen_0;
758 711836 : return q;
759 3735402 : case 1:
760 3735402 : if (z) *z = r; else cgiv(r);
761 3735402 : return q;
762 1041022 : case -1: break;
763 : }
764 1041022 : q = addis(q, -signe(y));
765 1041022 : if (!z) return gc_INT(av, q);
766 :
767 864924 : *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
768 864924 : return gc_all_unsafe(av,(pari_sp)r,2,&q,z);
769 : }
770 : GEN
771 139090958 : truedvmdis(GEN x, long y, GEN *z)
772 : {
773 139090958 : pari_sp av = avma;
774 : long r;
775 : GEN q;
776 :
777 139090958 : if (z == ONLY_REM) return modis(x, y);
778 139090958 : q = divis_rem(x,y,&r);
779 :
780 139090958 : if (r >= 0)
781 : {
782 126200771 : if (z) *z = utoi(r);
783 126200771 : return q;
784 : }
785 12890187 : q = gc_INT(av, addis(q, (y < 0)? 1: -1));
786 12890187 : if (z) *z = utoi(r + labs(y));
787 12890187 : 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 109744236 : int2n(long n) {
809 : long i, m, l;
810 : GEN z;
811 109744236 : if (n < 0) return gen_0;
812 109725224 : if (n == 0) return gen_1;
813 :
814 105970646 : l = dvmdsBIL(n, &m) + 3;
815 105970646 : z = cgetipos(l);
816 568853121 : for (i = 2; i < l; i++) z[i] = 0;
817 105970646 : *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 25996714 : int2u(ulong n) {
823 : ulong i, m, l;
824 : GEN z;
825 25996714 : if (n == 0) return gen_1;
826 :
827 25996714 : l = dvmduBIL(n, &m) + 3;
828 25996714 : z = cgetipos(l);
829 52386777 : for (i = 2; i < l; i++) z[i] = 0;
830 25996714 : *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 2323599674 : shifti(GEN x, long n)
849 : {
850 2323599674 : long s = signe(x);
851 : GEN y;
852 :
853 2323599674 : if(s == 0) return gen_0;
854 2015816677 : y = shiftispec(x + 2, lgefint(x) - 2, n);
855 2015816675 : if (signe(y)) setsigne(y, s);
856 2015816675 : return y;
857 : }
858 :
859 : /* actual operations will take place on a+2 and b+2: we strip the codewords */
860 : GEN
861 21996170510 : mulii(GEN a,GEN b)
862 : {
863 : long sa,sb;
864 : GEN z;
865 :
866 21996170510 : sa=signe(a); if (!sa) return gen_0;
867 14892789591 : sb=signe(b); if (!sb) return gen_0;
868 10914187628 : if (sb<0) sa = -sa;
869 10914187628 : z = muliispec(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
870 10914187628 : setsigne(z,sa); return z;
871 : }
872 :
873 : GEN
874 1963621213 : 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 258734435 : usqrt(ulong a)
880 : {
881 258734435 : ulong x = (ulong)sqrt((double)a);
882 : #ifdef LONG_IS_64BIT
883 224976589 : if (x > LOWMASK || x*x > a) x--;
884 : #endif
885 258734435 : return x;
886 : }
887 :
888 : /********************************************************************/
889 : /** **/
890 : /** EXPONENT / CONVERSION t_REAL --> double **/
891 : /** **/
892 : /********************************************************************/
893 :
894 : #ifdef LONG_IS_64BIT
895 : long
896 23761967 : dblexpo(double x)
897 : {
898 : union { double f; ulong i; } fi;
899 23761967 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
900 23761967 : const int exp_mid = 0x3ff;/* exponent bias */
901 :
902 23761967 : if (x==0.) return -exp_mid;
903 23761955 : fi.f = x;
904 23761955 : 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 11706747 : dbltor(double x)
920 : {
921 : GEN z;
922 : long e;
923 : union { double f; ulong i; } fi;
924 11706747 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
925 11706747 : const int exp_mid = 0x3ff;/* exponent bias */
926 11706747 : const int expo_len = 11; /* number of bits of exponent */
927 :
928 11706747 : if (x==0.) return real_0_bit(-exp_mid);
929 11700413 : fi.f = x; z = cgetr(DEFAULTPREC);
930 : {
931 11700413 : const ulong a = fi.i;
932 : ulong A;
933 11700413 : e = ((a & (HIGHBIT-1)) >> mant_len) - exp_mid;
934 11700413 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
935 11700413 : A = a << expo_len;
936 11700413 : 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 11700413 : z[2] = HIGHBIT | A;
944 11700413 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
945 : }
946 11700413 : return z;
947 : }
948 :
949 : double
950 268317732 : rtodbl(GEN x)
951 : {
952 268317732 : long ex,s=signe(x);
953 : ulong a;
954 : union { double f; ulong i; } fi;
955 268317732 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
956 268317732 : const int exp_mid = 0x3ff;/* exponent bias */
957 268317732 : const int expo_len = 11; /* number of bits of exponent */
958 :
959 268317732 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
960 :
961 : /* start by rounding to closest */
962 236293255 : a = (x[2] & (HIGHBIT-1)) + 0x400;
963 236293255 : if (a & HIGHBIT) { ex++; a=0; }
964 236293255 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
965 236293255 : fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
966 236293255 : if (s<0) fi.i |= HIGHBIT;
967 236293255 : 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 1897516 : dbltor(double x)
1013 : {
1014 : GEN z;
1015 : long e;
1016 : union { double f; ulong i[2]; } fi;
1017 1897516 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1018 1897516 : const int exp_mid = 0x3ff;/* exponent bias */
1019 1897516 : const int expo_len = 11; /* number of bits of exponent */
1020 1897516 : const int shift = mant_len-32;
1021 :
1022 1897516 : if (x==0.) return real_0_bit(-exp_mid);
1023 1896473 : fi.f = x; z = cgetr(DEFAULTPREC);
1024 : {
1025 1896473 : const ulong a = fi.i[INDEX0];
1026 1896473 : const ulong b = fi.i[INDEX1];
1027 : ulong A, B;
1028 1896473 : e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
1029 1896473 : if (e == exp_mid+1) pari_err_OVERFLOW("dbltor [NaN or Infinity]");
1030 1896473 : A = b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
1031 1896473 : B = b << expo_len;
1032 1896473 : 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 1896473 : z[3] = B;
1053 1896473 : z[2] = HIGHBIT | A;
1054 : }
1055 1896473 : z[1] = _evalexpo(e) | evalsigne(x<0? -1: 1);
1056 : }
1057 1896473 : return z;
1058 : }
1059 :
1060 : double
1061 44247131 : rtodbl(GEN x)
1062 : {
1063 44247131 : long ex,s=signe(x),lx=lg(x);
1064 : ulong a,b,k;
1065 : union { double f; ulong i[2]; } fi;
1066 44247131 : const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
1067 44247131 : const int exp_mid = 0x3ff;/* exponent bias */
1068 44247131 : const int expo_len = 11; /* number of bits of exponent */
1069 44247131 : const int shift = mant_len-32;
1070 :
1071 44247131 : if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
1072 :
1073 : /* start by rounding to closest */
1074 38825885 : a = x[2] & (HIGHBIT-1);
1075 38825885 : if (lx > 3)
1076 : {
1077 38698567 : b = x[3] + 0x400UL; if (b < 0x400UL) a++;
1078 38698567 : if (a & HIGHBIT) { ex++; a=0; }
1079 : }
1080 127318 : else b = 0;
1081 38825885 : if (ex >= exp_mid) pari_err_OVERFLOW("t_REAL->double conversion");
1082 38825885 : ex += exp_mid;
1083 38825885 : k = (a >> expo_len) | (ex << shift);
1084 38825885 : if (s<0) k |= HIGHBIT;
1085 38825885 : fi.i[INDEX0] = k;
1086 38825885 : fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
1087 38825885 : return fi.f;
1088 : }
1089 : #endif /* LONG_IS_64BIT */
1090 :
|
