Line data Source code
1 : #line 2 "../src/kernel/gmp/mp.c"
2 : /* Copyright (C) 2002-2003 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 : /***********************************************************************/
17 : /** **/
18 : /** GMP KERNEL **/
19 : /** BA2002Sep24 **/
20 : /***********************************************************************/
21 : /* GMP t_INT as just like normal t_INT, just the mantissa is the other way
22 : * round
23 : *
24 : * `How would you like to live in Looking-glass House, Kitty? I
25 : * wonder if they'd give you milk in there? Perhaps Looking-glass
26 : * milk isn't good to drink--But oh, Kitty! now we come to the
27 : * passage. You can just see a little PEEP of the passage in
28 : * Looking-glass House, if you leave the door of our drawing-room
29 : * wide open: and it's very like our passage as far as you can see,
30 : * only you know it may be quite different on beyond. Oh, Kitty!
31 : * how nice it would be if we could only get through into Looking-
32 : * glass House! I'm sure it's got, oh! such beautiful things in it!
33 : *
34 : * Through the Looking Glass, Lewis Carrol
35 : *
36 : * (pityful attempt to beat GN code/comments rate)
37 : * */
38 :
39 : #include <gmp.h>
40 : #include "pari.h"
41 : #include "paripriv.h"
42 : #include "../src/kernel/none/tune-gen.h"
43 :
44 : /*We need PARI invmod renamed to invmod_pari*/
45 : #define INVMOD_PARI
46 :
47 0 : static void *pari_gmp_realloc(void *ptr, size_t old_size, size_t new_size) {
48 0 : (void)old_size; return (void *) pari_realloc(ptr,new_size);
49 : }
50 :
51 1755994 : static void pari_gmp_free(void *ptr, size_t old_size){
52 1755994 : (void)old_size; pari_free(ptr);
53 1755994 : }
54 :
55 : static void *(*old_gmp_malloc)(size_t new_size);
56 : static void *(*old_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
57 : static void (*old_gmp_free)(void *ptr, size_t old_size);
58 :
59 : void
60 1112 : pari_kernel_init(void)
61 : {
62 1112 : mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
63 1112 : mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
64 1112 : }
65 :
66 : const char *
67 4 : pari_kernel_version(void)
68 : {
69 : #ifdef gmp_version
70 4 : return gmp_version;
71 : #else
72 : return "";
73 : #endif
74 : }
75 :
76 : void
77 1104 : pari_kernel_close(void)
78 : {
79 : void *(*new_gmp_malloc)(size_t new_size);
80 : void *(*new_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
81 : void (*new_gmp_free)(void *ptr, size_t old_size);
82 1104 : mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
83 1104 : if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
84 1104 : if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
85 1104 : if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
86 1104 : mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
87 1104 : }
88 :
89 : #define LIMBS(x) ((mp_limb_t *)((x)+2))
90 : #define NLIMBS(x) (lgefint(x)-2)
91 : /*This one is for t_REALs to emphasize they are not t_INTs*/
92 : #define RLIMBS(x) ((mp_limb_t *)((x)+2))
93 : #define RNLIMBS(x) (lg(x)-2)
94 :
95 : INLINE void
96 6890117 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
97 : {
98 57047270 : while (--n >= 0) x[n]=y[n];
99 6890117 : }
100 :
101 : INLINE void
102 587224867 : xmpn_mirror(mp_limb_t *x, long n)
103 : {
104 : long i;
105 3922044559 : for(i=0;i<(n>>1);i++)
106 : {
107 3334819692 : ulong m=x[i];
108 3334819692 : x[i]=x[n-1-i];
109 3334819692 : x[n-1-i]=m;
110 : }
111 587224867 : }
112 :
113 : INLINE void
114 718703914 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
115 : {
116 : long i;
117 9935818011 : for(i=0;i<n;i++)
118 9217114097 : z[i]=x[n-1-i];
119 718703914 : }
120 :
121 : INLINE void
122 237505011 : xmpn_zero(mp_ptr x, mp_size_t n)
123 : {
124 2075607562 : while (--n >= 0) x[n]=0;
125 237505011 : }
126 :
127 : INLINE GEN
128 41266738 : icopy_ef(GEN x, long l)
129 : {
130 41266738 : long lx = lgefint(x);
131 41266738 : const GEN y = cgeti(l);
132 :
133 294269258 : while (--lx > 0) y[lx]=x[lx];
134 41265336 : return y;
135 : }
136 :
137 : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
138 : * GENs but pairs (long *a, long na) representing a list of digits (in basis
139 : * BITS_IN_LONG) : a[0], ..., a[na-1]. In order to facilitate splitting: no
140 : * need to reintroduce codewords. */
141 :
142 : /***********************************************************************/
143 : /** **/
144 : /** ADDITION / SUBTRACTION **/
145 : /** **/
146 : /***********************************************************************/
147 :
148 : GEN
149 2997898 : setloop(GEN a)
150 : {
151 2997898 : pari_sp av = avma - 2 * sizeof(long);
152 2997898 : (void)cgetg(lgefint(a) + 3, t_VECSMALL);
153 2997898 : return icopy_avma(a, av); /* two cells of extra space after a */
154 : }
155 :
156 : /* we had a = setloop(?), then some incloops. Reset a to b */
157 : GEN
158 174328 : resetloop(GEN a, GEN b) {
159 174328 : a[0] = evaltyp(t_INT) | evallg(lgefint(b));
160 174328 : affii(b, a); return a;
161 : }
162 :
163 : /* assume a > 0, initialized by setloop. Do a++ */
164 : static GEN
165 103299962 : incpos(GEN a)
166 : {
167 103299962 : long i, l = lgefint(a);
168 103299967 : for (i=2; i<l; i++)
169 103304919 : if (++uel(a,i)) return a;
170 3 : a[l] = 1; l++;
171 3 : a[0]=evaltyp(t_INT) | _evallg(l);
172 3 : a[1]=evalsigne(1) | evallgefint(l);
173 3 : return a;
174 : }
175 :
176 : /* assume a < 0, initialized by setloop. Do a++ */
177 : static GEN
178 66684 : incneg(GEN a)
179 : {
180 66684 : long i, l = lgefint(a)-1;
181 66684 : if (uel(a,2)--)
182 : {
183 66680 : if (!a[l]) /* implies l = 2 */
184 : {
185 1980 : a[0] = evaltyp(t_INT) | _evallg(2);
186 1980 : a[1] = evalsigne(0) | evallgefint(2);
187 : }
188 66680 : return a;
189 : }
190 5 : for (i=3; i<=l; i++)
191 5 : if (uel(a,i)--) break;
192 4 : if (!a[l])
193 : {
194 4 : a[0] = evaltyp(t_INT) | _evallg(l);
195 4 : a[1] = evalsigne(-1) | evallgefint(l);
196 : }
197 4 : return a;
198 : }
199 :
200 : /* assume a initialized by setloop. Do a++ */
201 : GEN
202 103701952 : incloop(GEN a)
203 : {
204 103701952 : switch(signe(a))
205 : {
206 336652 : case 0:
207 336652 : a[0]=evaltyp(t_INT) | _evallg(3);
208 336652 : a[1]=evalsigne(1) | evallgefint(3);
209 336652 : a[2]=1; return a;
210 66684 : case -1: return incneg(a);
211 103298616 : default: return incpos(a);
212 : }
213 : }
214 :
215 : INLINE GEN
216 2782227608 : adduispec(ulong s, GEN x, long nx)
217 : {
218 : GEN zd;
219 : long lz;
220 :
221 2782227608 : if (nx == 1) return adduu(uel(x,0), s);
222 900714013 : lz = nx+3; zd = cgeti(lz);
223 903728398 : if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
224 5022012 : zd[lz-1]=1;
225 : else
226 898710232 : lz--;
227 903732244 : zd[1] = evalsigne(1) | evallgefint(lz);
228 903732244 : return zd;
229 : }
230 :
231 : GEN
232 764438320 : adduispec_offset(ulong s, GEN x, long offset, long nx)
233 : {
234 764438320 : GEN xd=x+2+offset;
235 996205236 : while (nx && *(xd+nx-1)==0) nx--;
236 764438320 : if (!nx) return utoi(s);
237 719703094 : return adduispec(s,xd,nx);
238 : }
239 :
240 : INLINE GEN
241 3332326573 : addiispec(GEN x, GEN y, long nx, long ny)
242 : {
243 : GEN zd;
244 : long lz;
245 :
246 3332326573 : if (nx < ny) swapspec(x,y, nx,ny);
247 3332326573 : if (ny == 1) return adduispec(*y,x,nx);
248 1356042186 : lz = nx+3; zd = cgeti(lz);
249 :
250 1356088644 : if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
251 30432363 : zd[lz-1]=1;
252 : else
253 1325934775 : lz--;
254 :
255 1356367138 : zd[1] = evalsigne(1) | evallgefint(lz);
256 1356367138 : return zd;
257 : }
258 :
259 : /* assume x >= y */
260 : INLINE GEN
261 1773943631 : subiuspec(GEN x, ulong s, long nx)
262 : {
263 : GEN zd;
264 : long lz;
265 :
266 1773943631 : if (nx == 1) return utoi(x[0] - s);
267 :
268 251788744 : lz = nx + 2; zd = cgeti(lz);
269 252257519 : mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
270 252258233 : if (! zd[lz - 1]) { --lz; }
271 :
272 252258233 : zd[1] = evalsigne(1) | evallgefint(lz);
273 252258233 : return zd;
274 : }
275 :
276 : /* assume x > y */
277 : INLINE GEN
278 3018774916 : subiispec(GEN x, GEN y, long nx, long ny)
279 : {
280 : GEN zd;
281 : long lz;
282 3018774916 : if (ny==1) return subiuspec(x,*y,nx);
283 1275040774 : lz = nx+2; zd = cgeti(lz);
284 :
285 1272766641 : mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
286 1761194383 : while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
287 :
288 1273029560 : zd[1] = evalsigne(1) | evallgefint(lz);
289 1273029560 : return zd;
290 : }
291 :
292 : static void
293 527272288 : roundr_up_ip(GEN x, long l)
294 : {
295 527272288 : long i = l;
296 : for(;;)
297 : {
298 529070909 : if (++((ulong*)x)[--i]) break;
299 2190431 : if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
300 : }
301 527323393 : }
302 :
303 : void
304 403940072 : affir(GEN x, GEN y)
305 : {
306 403940072 : const long s = signe(x), ly = lg(y);
307 : long lx, sh, i;
308 :
309 403940072 : if (!s)
310 : {
311 42185988 : y[1] = evalexpo(-bit_accuracy(ly));
312 42182247 : return;
313 : }
314 361754084 : lx = lgefint(x); sh = bfffo(*int_MSW(x));
315 361754084 : y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
316 361948946 : if (sh) {
317 354296414 : if (lx <= ly)
318 : {
319 961506750 : for (i=lx; i<ly; i++) y[i]=0;
320 249345107 : mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
321 249412273 : xmpn_mirror(LIMBS(y),lx-2);
322 249699600 : return;
323 : }
324 104951307 : mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
325 104952106 : uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
326 104952106 : xmpn_mirror(LIMBS(y),ly-2);
327 : /* lx > ly: round properly */
328 104949662 : if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
329 : }
330 : else {
331 7652532 : GEN xd=int_MSW(x);
332 7652532 : if (lx <= ly)
333 : {
334 9517865 : for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
335 5237161 : for ( ; i<ly; i++) y[i]=0;
336 2363881 : return;
337 : }
338 15117578 : for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
339 : /* lx > ly: round properly */
340 5288651 : if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
341 : }
342 : }
343 :
344 : INLINE GEN
345 715114362 : shiftispec(GEN x, long nx, long n)
346 : {
347 : long ny,m;
348 : GEN yd, y;
349 :
350 715114362 : if (!n) return icopyspec(x, nx);
351 :
352 685821662 : if (n > 0)
353 : {
354 399658059 : long d = dvmdsBIL(n, &m);
355 : long i;
356 :
357 399547790 : ny = nx + d + (m!=0);
358 399547790 : y = cgeti(ny + 2); yd = y + 2;
359 553255224 : for (i=0; i<d; i++) yd[i] = 0;
360 :
361 398604774 : if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
362 : else
363 : {
364 397063102 : ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
365 397126800 : if (carryd) yd[ny - 1] = carryd;
366 373850721 : else ny--;
367 : }
368 : }
369 : else
370 : {
371 286163603 : long d = dvmdsBIL(-n, &m);
372 :
373 288128434 : ny = nx - d;
374 288128434 : if (ny < 1) return gen_0;
375 285120252 : y = cgeti(ny + 2); yd = y + 2;
376 :
377 284805673 : if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
378 : else
379 : {
380 280142362 : mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
381 280155171 : if (yd[ny - 1] == 0)
382 : {
383 61011907 : if (ny == 1) return gc_const((pari_sp)(yd + 1), gen_0);
384 51396808 : ny--;
385 : }
386 : }
387 : }
388 671979110 : y[1] = evalsigne(1)|evallgefint(ny + 2);
389 671979110 : return y;
390 : }
391 :
392 : GEN
393 138195133 : mantissa2nr(GEN x, long n)
394 : {
395 138195133 : long ly, i, m, s = signe(x), lx = lg(x);
396 : GEN y;
397 138195133 : if (!s) return gen_0;
398 138193784 : if (!n)
399 : {
400 30684873 : y = cgeti(lx);
401 30681655 : y[1] = evalsigne(s) | evallgefint(lx);
402 30681655 : xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
403 30681453 : return y;
404 : }
405 107508911 : if (n > 0)
406 : {
407 218413 : GEN z = (GEN)avma;
408 218413 : long d = dvmdsBIL(n, &m);
409 :
410 218413 : ly = lx+d; y = new_chunk(ly);
411 557548 : for ( ; d; d--) *--z = 0;
412 222846 : if (!m) for (i=2; i<lx; i++) y[i]=x[i];
413 : else
414 : {
415 216971 : const ulong sh = BITS_IN_LONG - m;
416 216971 : shift_left(y,x, 2,lx-1, 0,m);
417 216971 : i = uel(x,2) >> sh;
418 : /* Extend y on the left? */
419 216971 : if (i) { ly++; y = new_chunk(1); y[2] = i; }
420 : }
421 : }
422 : else
423 : {
424 107290498 : ly = lx - dvmdsBIL(-n, &m);
425 107286784 : if (ly<3) return gen_0;
426 107286784 : y = new_chunk(ly);
427 107249902 : if (m) {
428 106991425 : shift_right(y,x, 2,ly, 0,m);
429 107028258 : if (y[2] == 0)
430 : {
431 0 : if (ly==3) return gc_const((pari_sp)(y+3), gen_0);
432 0 : ly--; set_avma((pari_sp)(++y));
433 : }
434 : } else {
435 701708 : for (i=2; i<ly; i++) y[i]=x[i];
436 : }
437 : }
438 107505148 : xmpn_mirror(LIMBS(y),ly-2);
439 107553653 : y[1] = evalsigne(s)|evallgefint(ly);
440 107553653 : y[0] = evaltyp(t_INT)|evallg(ly); return y;
441 : }
442 :
443 : GEN
444 3630460 : truncr(GEN x)
445 : {
446 : long s, e, d, m, i;
447 : GEN y;
448 3630460 : if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
449 1514724 : d = nbits2lg(e+1); m = remsBIL(e);
450 1514714 : if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
451 :
452 1514710 : y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
453 1514641 : if (++m == BITS_IN_LONG)
454 95989 : for (i=2; i<d; i++) y[d-i+1]=x[i];
455 : else
456 : {
457 1466905 : GEN z=cgeti(d);
458 3000407 : for (i=2; i<d; i++) z[d-i+1]=x[i];
459 1466805 : mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
460 1466805 : set_avma((pari_sp)y);
461 : }
462 1514508 : return y;
463 : }
464 :
465 : /* integral part */
466 : GEN
467 7132204 : floorr(GEN x)
468 : {
469 : long e, d, m, i, lx;
470 : GEN y;
471 7132204 : if (signe(x) >= 0) return truncr(x);
472 4243620 : if ((e=expo(x)) < 0) return gen_m1;
473 3523938 : d = nbits2lg(e+1); m = remsBIL(e);
474 3522502 : lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
475 3522502 : y = cgeti(d+1);
476 3511845 : if (++m == BITS_IN_LONG)
477 : {
478 3053 : for (i=2; i<d; i++) y[d-i+1]=x[i];
479 1450 : i=d; while (i<lx && !x[i]) i++;
480 906 : if (i==lx) goto END;
481 : }
482 : else
483 : {
484 3510939 : GEN z=cgeti(d);
485 7851596 : for (i=2; i<d; i++) z[d-i+1]=x[i];
486 3497680 : mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
487 3498283 : if (uel(x,d-1)<<m == 0)
488 : {
489 521066 : i=d; while (i<lx && !x[i]) i++;
490 121273 : if (i==lx) goto END;
491 : }
492 : }
493 3421171 : if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
494 0 : y[d++]=1;
495 3421032 : END:
496 3499050 : y[1] = evalsigne(-1) | evallgefint(d);
497 3499050 : return y;
498 : }
499 :
500 : INLINE int
501 3976032183 : cmpiispec(GEN x, GEN y, long lx, long ly)
502 : {
503 3976032183 : if (lx < ly) return -1;
504 3653275161 : if (lx > ly) return 1;
505 3504166981 : return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
506 : }
507 :
508 : INLINE int
509 269588528 : equaliispec(GEN x, GEN y, long lx, long ly)
510 : {
511 269588528 : if (lx != ly) return 0;
512 269447271 : return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
513 : }
514 :
515 : /***********************************************************************/
516 : /** **/
517 : /** MULTIPLICATION **/
518 : /** **/
519 : /***********************************************************************/
520 : /* assume ny > 0 */
521 : INLINE GEN
522 5603777981 : muluispec(ulong x, GEN y, long ny)
523 : {
524 5603777981 : if (ny == 1)
525 4505760259 : return muluu(x, *y);
526 : else
527 : {
528 1098017722 : long lz = ny+3;
529 1098017722 : GEN z = cgeti(lz);
530 1107390853 : ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
531 1108772375 : if (hi) { z[lz - 1] = hi; } else lz--;
532 1108772375 : z[1] = evalsigne(1) | evallgefint(lz);
533 1108772375 : return z;
534 : }
535 : }
536 :
537 : /* a + b*|y| */
538 : GEN
539 0 : addumului(ulong a, ulong b, GEN y)
540 : {
541 : GEN z;
542 : long i, lz;
543 : ulong hi;
544 0 : if (!b || !signe(y)) return utoi(a);
545 0 : lz = lgefint(y)+1;
546 0 : z = cgeti(lz);
547 0 : z[2]=a;
548 0 : for(i=3;i<lz;i++) z[i]=0;
549 0 : hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
550 0 : if (hi) z[lz-1]=hi; else lz--;
551 0 : z[1] = evalsigne(1) | evallgefint(lz);
552 0 : return gc_const((pari_sp)z, z);
553 : }
554 :
555 : /***********************************************************************/
556 : /** **/
557 : /** DIVISION **/
558 : /** **/
559 : /***********************************************************************/
560 :
561 : ulong
562 1324268251 : umodiu(GEN y, ulong x)
563 : {
564 1324268251 : long sy=signe(y);
565 : ulong hi;
566 1324268251 : if (!x) pari_err_INV("umodiu",gen_0);
567 1325239695 : if (!sy) return 0;
568 913865420 : hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
569 913865420 : if (!hi) return 0;
570 897268074 : return (sy > 0)? hi: x - hi;
571 : }
572 :
573 : /* return |y| \/ x */
574 : GEN
575 110615268 : absdiviu_rem(GEN y, ulong x, ulong *rem)
576 : {
577 : long ly;
578 : GEN z;
579 :
580 110615268 : if (!x) pari_err_INV("absdiviu_rem",gen_0);
581 110620933 : if (!signe(y)) { *rem = 0; return gen_0; }
582 :
583 83688431 : ly = lgefint(y);
584 83688431 : if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
585 :
586 68987125 : z = cgeti(ly);
587 68986878 : *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
588 68986955 : if (z [ly - 1] == 0) ly--;
589 68986955 : z[1] = evallgefint(ly) | evalsigne(1);
590 68986955 : return z;
591 : }
592 :
593 : GEN
594 86178643 : divis_rem(GEN y, long x, long *rem)
595 : {
596 86178643 : long sy=signe(y),ly,s;
597 : GEN z;
598 :
599 86178643 : if (!x) pari_err_INV("divis_rem",gen_0);
600 86189167 : if (!sy) { *rem = 0; return gen_0; }
601 60360308 : if (x<0) { s = -sy; x = -x; } else s = sy;
602 :
603 60360308 : ly = lgefint(y);
604 60360308 : if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
605 :
606 41570236 : z = cgeti(ly);
607 41565614 : *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
608 41565751 : if (sy<0) *rem = - *rem;
609 41565751 : if (z[ly - 1] == 0) ly--;
610 41565751 : z[1] = evallgefint(ly) | evalsigne(s);
611 41565751 : return z;
612 : }
613 :
614 : GEN
615 1006045 : divis(GEN y, long x)
616 : {
617 1006045 : long sy=signe(y),ly,s;
618 : GEN z;
619 :
620 1006045 : if (!x) pari_err_INV("divis",gen_0);
621 1006045 : if (!sy) return gen_0;
622 1005997 : if (x<0) { s = -sy; x = -x; } else s=sy;
623 :
624 1005997 : ly = lgefint(y);
625 1005997 : if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
626 :
627 1005685 : z = cgeti(ly);
628 1005682 : (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
629 1005682 : if (z[ly - 1] == 0) ly--;
630 1005682 : z[1] = evallgefint(ly) | evalsigne(s);
631 1005682 : return z;
632 : }
633 :
634 : /* We keep llx bits of x and lly bits of y*/
635 : static GEN
636 76136708 : divrr_with_gmp(GEN x, GEN y)
637 : {
638 76136708 : long lx=RNLIMBS(x),ly=RNLIMBS(y);
639 76136708 : long lw=minss(lx,ly);
640 76137667 : long lly=minss(lw+1,ly);
641 76139144 : GEN w = cgetg(lw+2, t_REAL);
642 76121787 : long lu=lw+lly;
643 76121787 : long llx=minss(lu,lx);
644 76119255 : mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
645 76098692 : mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
646 : mp_limb_t *q,*r;
647 76066218 : long e=expo(x)-expo(y);
648 76066218 : long sx=signe(x),sy=signe(y);
649 76066218 : xmpn_mirrorcopy(z,RLIMBS(y),lly);
650 76079966 : xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
651 76126685 : xmpn_zero(u,lu-llx);
652 76154737 : q = (mp_limb_t *)new_chunk(lw+1);
653 76138220 : r = (mp_limb_t *)new_chunk(lly);
654 :
655 76114453 : mpn_tdiv_qr(q,r,0,u,lu,z,lly);
656 :
657 : /*Round up: This is not exactly correct we should test 2*r>z*/
658 76162860 : if (uel(r,lly-1) > (uel(z,lly-1)>>1))
659 37776133 : mpn_add_1(q,q,lw+1,1);
660 :
661 76162899 : xmpn_mirrorcopy(RLIMBS(w),q,lw);
662 :
663 76161655 : if (q[lw] == 0) e--;
664 42025375 : else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
665 0 : else { w[2] = HIGHBIT; e++; }
666 76160527 : if (sy < 0) sx = -sx;
667 76160527 : w[1] = evalsigne(sx) | evalexpo(e);
668 76153680 : return gc_const((pari_sp)w, w);
669 : }
670 :
671 : /* We keep llx bits of x and lly bits of y*/
672 : static GEN
673 35250553 : divri_with_gmp(GEN x, GEN y)
674 : {
675 35250553 : long llx=RNLIMBS(x),ly=NLIMBS(y);
676 35250553 : long lly=minss(llx+1,ly);
677 35250875 : GEN w = cgetg(llx+2, t_REAL);
678 35246508 : long lu=llx+lly, ld=ly-lly;
679 35246508 : mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
680 35242781 : mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
681 : mp_limb_t *q,*r;
682 35237018 : long sh=bfffo(y[ly+1]);
683 35237018 : long e=expo(x)-expi(y);
684 35237782 : long sx=signe(x),sy=signe(y);
685 35237782 : if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
686 685153 : else xmpn_copy(z,LIMBS(y)+ld,lly);
687 35239542 : xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
688 35246994 : xmpn_zero(u,lu-llx);
689 35253402 : q = (mp_limb_t *)new_chunk(llx+1);
690 35250839 : r = (mp_limb_t *)new_chunk(lly);
691 :
692 35246120 : mpn_tdiv_qr(q,r,0,u,lu,z,lly);
693 :
694 : /*Round up: This is not exactly correct we should test 2*r>z*/
695 35255923 : if (uel(r,lly-1) > (uel(z,lly-1)>>1))
696 16045077 : mpn_add_1(q,q,llx+1,1);
697 :
698 35255924 : xmpn_mirrorcopy(RLIMBS(w),q,llx);
699 :
700 35255579 : if (q[llx] == 0) e--;
701 10628021 : else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
702 0 : else { w[2] = HIGHBIT; e++; }
703 35255458 : if (sy < 0) sx = -sx;
704 35255458 : w[1] = evalsigne(sx) | evalexpo(e);
705 35254211 : return gc_const((pari_sp)w, w);
706 : }
707 :
708 : GEN
709 151397087 : divri(GEN x, GEN y)
710 : {
711 151397087 : long s = signe(y);
712 :
713 151397087 : if (!s) pari_err_INV("divri",gen_0);
714 151397248 : if (!signe(x)) return real_0_bit(expo(x) - expi(y));
715 151160625 : if (!is_bigint(y)) {
716 115915700 : GEN z = divru(x, y[2]);
717 115913471 : if (s < 0) togglesign(z);
718 115913478 : return z;
719 : }
720 35244081 : return divri_with_gmp(x,y);
721 : }
722 :
723 : GEN
724 142776624 : divrr(GEN x, GEN y)
725 : {
726 142776624 : long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
727 : ulong y0,y1;
728 : GEN r, r1;
729 :
730 142776624 : if (!sy) pari_err_INV("divrr",y);
731 142784420 : e = expo(x) - expo(y);
732 142784420 : if (!sx) return real_0_bit(e);
733 142278276 : if (sy<0) sx = -sx;
734 :
735 142278276 : lx=lg(x); ly=lg(y);
736 142278276 : if (ly==3)
737 : {
738 26336966 : ulong k = x[2], l = (lx>3)? x[3]: 0;
739 : LOCAL_HIREMAINDER;
740 26336966 : if (k < uel(y,2)) e--;
741 : else
742 : {
743 8290653 : l >>= 1; if (k&1) l |= HIGHBIT;
744 8290653 : k >>= 1;
745 : }
746 26336966 : hiremainder = k; k = divll(l,y[2]);
747 26336966 : if (hiremainder > (uel(y,2) >> 1) && !++k) { k = HIGHBIT; e++; }
748 26336966 : r = cgetg(3, t_REAL);
749 26335754 : r[1] = evalsigne(sx) | evalexpo(e);
750 26333755 : r[2] = k; return r;
751 : }
752 :
753 115941310 : if (ly >= prec2lg(DIVRR_GMP_LIMIT))
754 76136487 : return divrr_with_gmp(x,y);
755 :
756 39852128 : lr = minss(lx,ly); r = new_chunk(lr);
757 39853632 : r1 = r-1;
758 135488141 : r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
759 39853632 : r1[lr] = (lx>ly)? x[lr]: 0;
760 39853632 : y0 = y[2]; y1 = y[3];
761 175306270 : for (i=0; i<lr-1; i++)
762 : { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
763 : ulong k, qp;
764 : LOCAL_HIREMAINDER;
765 : LOCAL_OVERFLOW;
766 :
767 135452638 : if (uel(r1,1) == y0) { qp = ULONG_MAX; k = addll(y0,r1[2]); }
768 : else
769 : {
770 135199646 : if (uel(r1,1) > y0) /* can't happen if i=0 */
771 : {
772 0 : GEN y1 = y+1;
773 0 : j = lr-i; r1[j] = subll(r1[j],y1[j]);
774 0 : for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
775 0 : j=i; do uel(r,--j)++; while (j && !r[j]);
776 : }
777 135199646 : hiremainder = r1[1]; overflow = 0;
778 135199646 : qp = divll(r1[2],y0); k = hiremainder;
779 : }
780 135452638 : j = lr-i+1;
781 135452638 : if (!overflow)
782 : {
783 : long k3, k4;
784 135294623 : k3 = mulll(qp,y1);
785 135294623 : if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
786 39907830 : k4 = subll(hiremainder,k);
787 : else
788 : {
789 95386793 : k3 = subll(k3, r1[3]);
790 95386793 : k4 = subllx(hiremainder,k);
791 : }
792 173705690 : while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
793 : }
794 135452638 : if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
795 402127263 : for (j--; j>1; j--)
796 : {
797 266674625 : r1[j] = subll(r1[j], addmul(qp,y[j]));
798 266674625 : hiremainder += overflow;
799 : }
800 135452638 : if (uel(r1,1) != hiremainder)
801 : {
802 187473 : if (uel(r1,1) < hiremainder)
803 : {
804 187473 : qp--;
805 187473 : j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
806 534265 : for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
807 : }
808 : else
809 : {
810 0 : uel(r1,1) -= hiremainder;
811 0 : while (r1[1])
812 : {
813 0 : qp++; if (!qp) { j=i; do uel(r,--j)++; while (j && !r[j]); }
814 0 : j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
815 0 : for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
816 0 : uel(r1,1) -= overflow;
817 : }
818 : }
819 : }
820 135452638 : *++r1 = qp;
821 : }
822 : /* i = lr-1 */
823 : /* round correctly */
824 39853632 : if (uel(r1,1) > (y0>>1))
825 : {
826 20132890 : j=i; do uel(r,--j)++; while (j && !r[j]);
827 : }
828 135672195 : r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
829 39853632 : if (r[0] == 0) e--;
830 14281076 : else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
831 : else { /* possible only when rounding up to 0x2 0x0 ... */
832 18 : r[2] = (long)HIGHBIT; e++;
833 : }
834 39853013 : r[0] = evaltyp(t_REAL)|evallg(lr);
835 39912602 : r[1] = evalsigne(sx) | evalexpo(e);
836 39905899 : return r;
837 : }
838 :
839 : /* Integer division x / y: such that sign(r) = sign(x)
840 : * if z = ONLY_REM return remainder, otherwise return quotient
841 : * if z != NULL set *z to remainder
842 : * *z is the last object on stack and can be disposed of with cgiv
843 : * If *z is zero, we put gen_0 here and no copy.
844 : * space needed: lx + ly
845 : */
846 : GEN
847 2196328588 : dvmdii(GEN x, GEN y, GEN *z)
848 : {
849 2196328588 : long sx=signe(x),sy=signe(y);
850 : long lx, ly, lq;
851 : pari_sp av;
852 : GEN r,q;
853 :
854 2196328588 : if (!sy) pari_err_INV("dvmdii",y);
855 2197438748 : if (!sx)
856 : {
857 69917458 : if (!z || z == ONLY_REM) return gen_0;
858 41103097 : *z=gen_0; return gen_0;
859 : }
860 2127521290 : lx=lgefint(x);
861 2127521290 : ly=lgefint(y); lq=lx-ly;
862 2127521290 : if (lq <= 0)
863 : {
864 1260835160 : if (lq == 0)
865 : {
866 1072018017 : long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
867 1072018017 : if (s>0) goto DIVIDE;
868 369649028 : if (s==0)
869 : {
870 33439277 : if (z == ONLY_REM) return gen_0;
871 22327793 : if (z) *z = gen_0;
872 22327793 : if (sx < 0) sy = -sy;
873 22327793 : return stoi(sy);
874 : }
875 : }
876 525026894 : if (z == ONLY_REM) return icopy(x);
877 14332338 : if (z) *z = icopy(x);
878 14332338 : return gen_0;
879 : }
880 866686130 : DIVIDE: /* quotient is nonzero */
881 1569055119 : av=avma; if (sx<0) sy = -sy;
882 1569055119 : if (ly==3)
883 : {
884 574058400 : ulong lq = lx;
885 : ulong si;
886 574058400 : q = cgeti(lq);
887 573618069 : si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
888 574149427 : if (q[lq - 1] == 0) lq--;
889 574149427 : if (z == ONLY_REM)
890 : {
891 333975039 : if (!si) return gc_const(av, gen_0);
892 287653256 : set_avma(av); r = cgeti(3);
893 287106106 : r[1] = evalsigne(sx) | evallgefint(3);
894 287106106 : r[2] = si; return r;
895 : }
896 240174388 : q[1] = evalsigne(sy) | evallgefint(lq);
897 240174388 : if (!z) return q;
898 236066272 : if (!si) { *z=gen_0; return q; }
899 64200936 : r=cgeti(3);
900 64228453 : r[1] = evalsigne(sx) | evallgefint(3);
901 64228453 : r[2] = si; *z=r; return q;
902 : }
903 994996719 : if (z == ONLY_REM)
904 : {
905 971887481 : ulong lr = lgefint(y);
906 971887481 : ulong lq = lgefint(x)-lgefint(y)+3;
907 971887481 : GEN r = cgeti(lr);
908 965305061 : GEN q = cgeti(lq);
909 958823805 : mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
910 972536643 : if (!r[lr - 1])
911 : {
912 1182225156 : while(lr>2 && !r[lr - 1]) lr--;
913 545213151 : if (lr == 2) return gc_const(av, gen_0); /* exact division */
914 : }
915 957599058 : r[1] = evalsigne(sx) | evallgefint(lr);
916 957599058 : return gc_const((pari_sp)r, r);
917 : }
918 : else
919 : {
920 23109238 : ulong lq = lgefint(x)-lgefint(y)+3;
921 23109238 : ulong lr = lgefint(y);
922 23109238 : GEN q = cgeti(lq);
923 28641531 : GEN r = cgeti(lr);
924 28625216 : mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
925 28654173 : if (q[lq - 1] == 0) lq--;
926 28654173 : q[1] = evalsigne(sy) | evallgefint(lq);
927 28654173 : if (!z) return gc_const((pari_sp)q, q);
928 24950234 : if (!r[lr - 1])
929 : {
930 38510854 : while(lr>2 && !r[lr - 1]) lr--;
931 6343283 : if (lr == 2) { *z = gen_0; return gc_const((pari_sp)q, q); } /* exact */
932 : }
933 20169278 : r[1] = evalsigne(sx) | evallgefint(lr);
934 20169278 : *z = gc_const((pari_sp)r, r); return q;
935 : }
936 : }
937 :
938 : /* Montgomery reduction.
939 : * N has k words, assume T >= 0 has less than 2k.
940 : * Return res := T / B^k mod N, where B = 2^BIL
941 : * such that 0 <= res < T/B^k + N and res has less than k words */
942 : GEN
943 35037453 : red_montgomery(GEN T, GEN N, ulong inv)
944 : {
945 : pari_sp av;
946 : GEN Te, Td, Ne, Nd, scratch;
947 35037453 : ulong i, j, m, t, d, k = NLIMBS(N);
948 : int carry;
949 : LOCAL_HIREMAINDER;
950 : LOCAL_OVERFLOW;
951 :
952 35037453 : if (k == 0) return gen_0;
953 35037453 : d = NLIMBS(T); /* <= 2*k */
954 35037453 : if (d == 0) return gen_0;
955 : #ifdef DEBUG
956 : if (d > 2*k) pari_err_BUG("red_montgomery");
957 : #endif
958 35037444 : if (k == 1)
959 : { /* as below, special cased for efficiency */
960 163341 : ulong n = uel(N,2);
961 163341 : if (d == 1) {
962 163194 : hiremainder = uel(T,2);
963 163194 : m = hiremainder * inv;
964 163194 : (void)addmul(m, n); /* t + m*n = 0 */
965 163194 : return utoi(hiremainder);
966 : } else { /* d = 2 */
967 147 : hiremainder = uel(T,2);
968 147 : m = hiremainder * inv;
969 147 : (void)addmul(m, n); /* t + m*n = 0 */
970 147 : t = addll(hiremainder, uel(T,3));
971 147 : if (overflow) t -= n; /* t > n doesn't fit in 1 word */
972 147 : return utoi(t);
973 : }
974 : }
975 : /* assume k >= 2 */
976 34874103 : av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
977 :
978 : /* copy T to scratch space (pad with zeroes to 2k words) */
979 34479071 : Td = scratch;
980 34479071 : Te = T + 2;
981 482615678 : for (i=0; i < d ; i++) *Td++ = *Te++;
982 61663552 : for ( ; i < (k<<1); i++) *Td++ = 0;
983 :
984 34479071 : Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
985 34479071 : Ne = N + 1; /* 1 beyond end of N mantissa */
986 :
987 34479071 : carry = 0;
988 257606754 : for (i=0; i<k; i++) /* set T := T/B nod N, k times */
989 : {
990 223127683 : Td = Te; /* one beyond end of (new) T mantissa */
991 223127683 : Nd = Ne;
992 223127683 : hiremainder = *++Td;
993 223127683 : m = hiremainder * inv; /* solve T + m N = O(B) */
994 :
995 : /* set T := (T + mN) / B */
996 223127683 : Te = Td;
997 223127683 : (void)addmul(m, *++Nd); /* = 0 */
998 1934727282 : for (j=1; j<k; j++)
999 : {
1000 1711599599 : t = addll(addmul(m, *++Nd), *++Td);
1001 1711599599 : *Td = t;
1002 1711599599 : hiremainder += overflow;
1003 : }
1004 223127683 : t = addll(hiremainder, *++Td); *Td = t + carry;
1005 223127683 : carry = (overflow || (carry && *Td == 0));
1006 : }
1007 34479071 : if (carry)
1008 : { /* Td > N overflows (k+1 words), set Td := Td - N */
1009 76306 : GEN NE = N + k+1;
1010 76306 : Td = Te;
1011 76306 : Nd = Ne;
1012 76306 : t = subll(*++Td, *++Nd); *Td = t;
1013 768715 : while (Nd < NE) { t = subllx(*++Td, *++Nd); *Td = t; }
1014 : }
1015 :
1016 : /* copy result */
1017 34479071 : Td = (GEN)av - 1; /* *Td = high word of final result */
1018 38086717 : while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
1019 34479071 : k = Td - Te; if (!k) return gc_const(av, gen_0);
1020 34479071 : Td = (GEN)av - k; /* will write mantissa there */
1021 34479071 : (void)memmove(Td, Te+1, k*sizeof(long));
1022 34479071 : Td -= 2;
1023 34479071 : Td[0] = evaltyp(t_INT) | evallg(k+2);
1024 34788370 : Td[1] = evalsigne(1) | evallgefint(k+2);
1025 : #ifdef DEBUG
1026 : {
1027 : long l = NLIMBS(N), s = BITS_IN_LONG*l;
1028 : GEN R = int2n(s);
1029 : GEN res = remii(mulii(T, Fp_inv(R, N)), N);
1030 : if (k > lgefint(N)
1031 : || !equalii(remii(Td,N),res)
1032 : || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
1033 : }
1034 : #endif
1035 34788370 : return gc_const((pari_sp)Td, Td);
1036 : }
1037 :
1038 : /* EXACT INTEGER DIVISION */
1039 :
1040 : /* use undocumented GMP interface */
1041 : static void
1042 115387697 : GEN2mpz(mpz_t X, GEN x)
1043 : {
1044 115387697 : long l = lgefint(x)-2;
1045 115387697 : X->_mp_alloc = l;
1046 115387697 : X->_mp_size = signe(x) > 0? l: -l;
1047 115387697 : X->_mp_d = LIMBS(x);
1048 115387697 : }
1049 : static void
1050 57695135 : mpz2GEN(GEN z, mpz_t Z)
1051 : {
1052 57695135 : long l = Z->_mp_size;
1053 57695135 : z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
1054 57695135 : }
1055 :
1056 : #ifdef mpn_divexact_1
1057 : static GEN
1058 418957871 : diviuexact_i(GEN x, ulong y)
1059 : {
1060 418957871 : long l = lgefint(x);
1061 418957871 : GEN z = cgeti(l);
1062 418464785 : mpn_divexact_1(LIMBS(z), LIMBS(x), NLIMBS(x), y);
1063 418516219 : if (z[l-1] == 0) l--;
1064 418516219 : z[1] = evallgefint(l) | evalsigne(signe(x));
1065 418516219 : return z;
1066 : }
1067 : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly */
1068 : /* assume y != 0 and the division is exact */
1069 : static GEN
1070 : diviuexact_i(GEN x, ulong y)
1071 : {
1072 : long l = lgefint(x);
1073 : GEN z = cgeti(l);
1074 : mpz_t X, Z;
1075 : GEN2mpz(X, x);
1076 : Z->_mp_alloc = l-2;
1077 : Z->_mp_size = l-2;
1078 : Z->_mp_d = LIMBS(z);
1079 : mpz_divexact_ui(Z, X, y);
1080 : mpz2GEN(z, Z); return z;
1081 : }
1082 : #else
1083 : /* assume y != 0 and the division is exact */
1084 : static GEN
1085 : diviuexact_i(GEN x, ulong y)
1086 : {
1087 : /*TODO: implement true exact division.*/
1088 : return divii(x,utoi(y));
1089 : }
1090 : #endif
1091 :
1092 : GEN
1093 39543611 : diviuexact(GEN x, ulong y)
1094 : {
1095 : GEN z;
1096 39543611 : if (!signe(x)) return gen_0;
1097 38415035 : z = diviuexact_i(x, y);
1098 38411214 : if (lgefint(z) == 2) pari_err_OP("exact division", x, utoi(y));
1099 38411270 : return z;
1100 : }
1101 :
1102 : /* Find z such that x=y*z, knowing that y | x (unchecked) */
1103 : GEN
1104 531472525 : diviiexact(GEN x, GEN y)
1105 : {
1106 : GEN z;
1107 531472525 : if (!signe(y)) pari_err_INV("diviiexact",y);
1108 531648421 : if (!signe(x)) return gen_0;
1109 438443228 : if (lgefint(y) == 3)
1110 : {
1111 380768603 : z = diviuexact_i(x, y[2]);
1112 380315101 : if (signe(y) < 0) togglesign(z);
1113 : }
1114 : else
1115 : {
1116 57674625 : long l = lgefint(x);
1117 : mpz_t X, Y, Z;
1118 57674625 : z = cgeti(l);
1119 57694019 : GEN2mpz(X, x);
1120 57693929 : GEN2mpz(Y, y);
1121 57693801 : Z->_mp_alloc = l-2;
1122 57693801 : Z->_mp_size = l-2;
1123 57693801 : Z->_mp_d = LIMBS(z);
1124 57693801 : mpz_divexact(Z, X, Y);
1125 57695155 : mpz2GEN(z, Z);
1126 : }
1127 438010226 : if (lgefint(z) == 2) pari_err_OP("exact division", x, y);
1128 437986634 : return z;
1129 : }
1130 :
1131 : /* assume yz != and yz | x */
1132 : GEN
1133 200513 : diviuuexact(GEN x, ulong y, ulong z)
1134 : {
1135 : long tmp[4];
1136 : ulong t;
1137 : LOCAL_HIREMAINDER;
1138 200513 : t = mulll(y, z);
1139 200513 : if (!hiremainder) return diviuexact(x, t);
1140 0 : tmp[0] = evaltyp(t_INT)|_evallg(4);
1141 0 : tmp[1] = evalsigne(1)|evallgefint(4);
1142 0 : tmp[2] = t;
1143 0 : tmp[3] = hiremainder;
1144 0 : return diviiexact(x, tmp);
1145 : }
1146 :
1147 : /********************************************************************/
1148 : /** **/
1149 : /** INTEGER MULTIPLICATION **/
1150 : /** **/
1151 : /********************************************************************/
1152 :
1153 : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
1154 : GEN
1155 5938811442 : muliispec(GEN x, GEN y, long nx, long ny)
1156 : {
1157 : GEN zd;
1158 : long lz;
1159 : ulong hi;
1160 :
1161 5938811442 : if (nx < ny) swapspec(x,y, nx,ny);
1162 5938811442 : if (!ny) return gen_0;
1163 5938811442 : if (ny == 1) return muluispec((ulong)*y, x, nx);
1164 :
1165 1057649240 : lz = nx+ny+2;
1166 1057649240 : zd = cgeti(lz);
1167 1060470060 : hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
1168 1067773071 : if (!hi) lz--;
1169 : /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
1170 :
1171 1067773071 : zd[1] = evalsigne(1) | evallgefint(lz);
1172 1067773071 : return zd;
1173 : }
1174 : GEN
1175 222738 : muluui(ulong x, ulong y, GEN z)
1176 : {
1177 222738 : long t, s = signe(z);
1178 : GEN r;
1179 : LOCAL_HIREMAINDER;
1180 :
1181 222738 : if (!x || !y || !signe(z)) return gen_0;
1182 222355 : t = mulll(x,y);
1183 222355 : if (!hiremainder)
1184 222365 : r = muluispec(t, z+2, lgefint(z)-2);
1185 : else
1186 : {
1187 : long tmp[2];
1188 0 : tmp[1] = hiremainder;
1189 0 : tmp[0] = t;
1190 0 : r = muliispec(z+2,tmp, lgefint(z)-2, 2);
1191 : }
1192 222342 : setsigne(r,s); return r;
1193 : }
1194 :
1195 : GEN
1196 1030321959 : sqrispec(GEN x, long nx)
1197 : {
1198 : GEN zd;
1199 : long lz;
1200 :
1201 1030321959 : if (!nx) return gen_0;
1202 493295612 : if (nx==1) return sqru(*x);
1203 :
1204 279456669 : lz = (nx<<1)+2;
1205 279456669 : zd = cgeti(lz);
1206 : #ifdef mpn_sqr
1207 276102895 : mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
1208 : #else
1209 : mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
1210 : #endif
1211 279829452 : if (zd[lz-1]==0) lz--;
1212 :
1213 279829452 : zd[1] = evalsigne(1) | evallgefint(lz);
1214 279829452 : return zd;
1215 : }
1216 :
1217 : INLINE GEN
1218 41408779 : sqrispec_mirror(GEN x, long nx)
1219 : {
1220 41408779 : GEN cx=new_chunk(nx);
1221 : GEN z;
1222 41360699 : xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
1223 41463997 : z=sqrispec(cx, nx);
1224 41534585 : xmpn_mirror(LIMBS(z), NLIMBS(z));
1225 41536077 : return z;
1226 : }
1227 :
1228 : /* leaves garbage on the stack. */
1229 : INLINE GEN
1230 85120199 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
1231 : {
1232 : GEN cx, cy, z;
1233 85120199 : long s = 0;
1234 116259234 : while (nx && x[nx-1]==0) { nx--; s++; }
1235 122875953 : while (ny && y[ny-1]==0) { ny--; s++; }
1236 85120199 : cx=new_chunk(nx); cy=new_chunk(ny);
1237 84651794 : xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
1238 85573392 : xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
1239 86120137 : z = nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
1240 86105665 : if (s)
1241 : {
1242 7856744 : long i, lz = lgefint(z) + s;
1243 7856744 : (void)new_chunk(s);
1244 7856742 : z -= s;
1245 76751528 : for (i=0; i<s; i++) z[2+i]=0;
1246 7856742 : z[1] = evalsigne(1) | evallgefint(lz);
1247 7856742 : z[0] = evaltyp(t_INT) | evallg(lz);
1248 : }
1249 86105662 : xmpn_mirror(LIMBS(z), NLIMBS(z));
1250 86823806 : return z;
1251 : }
1252 :
1253 : /* x % (2^n), assuming n >= 0 */
1254 : GEN
1255 38999084 : remi2n(GEN x, long n)
1256 : {
1257 : ulong hi;
1258 38999084 : long l, k, lx, ly, sx = signe(x);
1259 : GEN z, xd, zd;
1260 :
1261 38999084 : if (!sx || !n) return gen_0;
1262 :
1263 38678155 : k = dvmdsBIL(n, &l);
1264 38683534 : lx = lgefint(x);
1265 38683534 : if (lx < k+3) return icopy(x);
1266 :
1267 37722187 : xd = x + (2 + k);
1268 : /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
1269 : * ^--- initial xd */
1270 37722187 : hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
1271 37722187 : if (!hi)
1272 : { /* strip leading zeroes from result */
1273 4251366 : xd--; while (k && !*xd) { k--; xd--; }
1274 4051925 : if (!k) return gen_0;
1275 2184440 : ly = k+2;
1276 : }
1277 : else
1278 33670262 : ly = k+3;
1279 :
1280 35854702 : zd = z = cgeti(ly);
1281 35814323 : *++zd = evalsigne(sx) | evallgefint(ly);
1282 535795721 : xd = x+1; for ( ;k; k--) *++zd = *++xd;
1283 35814323 : if (hi) *++zd = hi;
1284 35814323 : return z;
1285 : }
1286 :
1287 : /********************************************************************/
1288 : /** **/
1289 : /** INTEGER SQUARE ROOT **/
1290 : /** **/
1291 : /********************************************************************/
1292 :
1293 : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
1294 : * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
1295 : * remainder is 0. R = NULL is allowed. */
1296 : GEN
1297 5182373 : sqrtremi(GEN a, GEN *r)
1298 : {
1299 5182373 : long l, na = NLIMBS(a);
1300 : mp_size_t nr;
1301 : GEN S;
1302 5182373 : if (!na) {
1303 724 : if (r) *r = gen_0;
1304 724 : return gen_0;
1305 : }
1306 5181649 : l = (na + 5) >> 1; /* 2 + ceil(na/2) */
1307 5181649 : S = cgetipos(l);
1308 5181613 : if (r) {
1309 1310197 : GEN R = cgeti(2 + na);
1310 1310197 : nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
1311 1310197 : if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
1312 26517 : else { set_avma((pari_sp)S); R = gen_0; }
1313 1310197 : *r = R;
1314 : }
1315 : else
1316 3871416 : (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
1317 5181617 : return S;
1318 : }
1319 :
1320 : /* compute sqrt(|a|), assuming a != 0 */
1321 : GEN
1322 126388176 : sqrtr_abs(GEN a)
1323 : {
1324 : GEN res;
1325 : mp_limb_t *b, *c;
1326 126388176 : long l = RNLIMBS(a), e = expo(a), er = e>>1;
1327 : long n;
1328 126388176 : res = cgetg(2 + l, t_REAL);
1329 126302798 : res[1] = evalsigne(1) | evalexpo(er);
1330 126379092 : if (e&1)
1331 : {
1332 53329891 : b = (mp_limb_t *) new_chunk(l<<1);
1333 53315200 : xmpn_zero(b,l);
1334 53315829 : xmpn_mirrorcopy(b+l, RLIMBS(a), l);
1335 53329202 : c = (mp_limb_t *) new_chunk(l);
1336 53322600 : n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
1337 53358714 : if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
1338 : }
1339 : else
1340 : {
1341 : ulong u;
1342 73049201 : b = (mp_limb_t *) mantissa2nr(a,-1);
1343 73083991 : b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
1344 73083991 : b = (mp_limb_t *) new_chunk(l);
1345 73059874 : xmpn_zero(b,l+1); /* overwrites the former b[0] */
1346 73062877 : c = (mp_limb_t *) new_chunk(l + 1);
1347 73024727 : n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
1348 73124055 : u = (ulong)*c++;
1349 73124055 : if ( u&HIGHBIT || (u == ~HIGHBIT &&
1350 0 : (n>l || (n==l && mpn_cmp(b,c,l)>0))))
1351 36025221 : mpn_add_1(c,c,l,1);
1352 : }
1353 126491593 : xmpn_mirrorcopy(RLIMBS(res),c,l);
1354 126465681 : return gc_const((pari_sp)res, res);
1355 : }
1356 :
1357 : /* Normalize a nonnegative integer */
1358 : GEN
1359 307430168 : int_normalize(GEN x, long known_zero_words)
1360 : {
1361 307430168 : long i = lgefint(x) - 1 - known_zero_words;
1362 2774633330 : for ( ; i > 1; i--)
1363 2722738304 : if (x[i]) { setlgefint(x, i+1); return x; }
1364 51895026 : x[1] = evalsigne(0) | evallgefint(2); return x;
1365 : }
1366 :
1367 : /********************************************************************
1368 : ** **
1369 : ** Base Conversion **
1370 : ** **
1371 : ********************************************************************/
1372 :
1373 : ulong *
1374 439114 : convi(GEN x, long *l)
1375 : {
1376 439114 : long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
1377 439114 : GEN str = cgetg(n+1, t_VECSMALL);
1378 439114 : unsigned char *res = (unsigned char*) GSTR(str);
1379 439114 : long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
1380 : long lz;
1381 : ulong *z;
1382 : long i, j;
1383 : unsigned char *t;
1384 439114 : while (!*res) {res++; llz--;} /*Strip leading zeros*/
1385 439114 : lz = (8+llz)/9;
1386 439114 : z = (ulong*)new_chunk(1+lz);
1387 439114 : t=res+llz+9;
1388 870532 : for(i=0;i<llz-8;i+=9)
1389 : {
1390 : ulong s;
1391 431418 : t-=18;
1392 431418 : s=*t++;
1393 3882762 : for (j=1; j<9;j++)
1394 3451344 : s=10*s+*t++;
1395 431418 : *z++=s;
1396 : }
1397 439114 : if (i<llz)
1398 : {
1399 435131 : unsigned char *t = res;
1400 435131 : ulong s=*t++;
1401 1229610 : for (j=i+1; j<llz;j++)
1402 794479 : s=10*s+*t++;
1403 435131 : *z++=s;
1404 : }
1405 439114 : *l = lz;
1406 439114 : return z;
1407 : }
|
