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. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /***********************************************************************/
16 : /** **/
17 : /** GMP KERNEL **/
18 : /** BA2002Sep24 **/
19 : /***********************************************************************/
20 : /* GMP t_INT as just like normal t_INT, just the mantissa is the other way
21 : * round
22 : *
23 : * `How would you like to live in Looking-glass House, Kitty? I
24 : * wonder if they'd give you milk in there? Perhaps Looking-glass
25 : * milk isn't good to drink--But oh, Kitty! now we come to the
26 : * passage. You can just see a little PEEP of the passage in
27 : * Looking-glass House, if you leave the door of our drawing-room
28 : * wide open: and it's very like our passage as far as you can see,
29 : * only you know it may be quite different on beyond. Oh, Kitty!
30 : * how nice it would be if we could only get through into Looking-
31 : * glass House! I'm sure it's got, oh! such beautiful things in it!
32 : *
33 : * Through the Looking Glass, Lewis Carrol
34 : *
35 : * (pityful attempt to beat GN code/comments rate)
36 : * */
37 :
38 : #include <gmp.h>
39 : #include "pari.h"
40 : #include "paripriv.h"
41 : #include "../src/kernel/none/tune-gen.h"
42 :
43 : /*We need PARI invmod renamed to invmod_pari*/
44 : #define INVMOD_PARI
45 :
46 0 : static void *pari_gmp_realloc(void *ptr, size_t old_size, size_t new_size) {
47 0 : (void)old_size; return (void *) pari_realloc(ptr,new_size);
48 : }
49 :
50 268727 : static void pari_gmp_free(void *ptr, size_t old_size){
51 268727 : (void)old_size; pari_free(ptr);
52 268727 : }
53 :
54 : static void *(*old_gmp_malloc)(size_t new_size);
55 : static void *(*old_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
56 : static void (*old_gmp_free)(void *ptr, size_t old_size);
57 :
58 : void
59 916 : pari_kernel_init(void)
60 : {
61 916 : mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
62 916 : mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
63 916 : }
64 :
65 : void
66 916 : pari_kernel_close(void)
67 : {
68 : void *(*new_gmp_malloc)(size_t new_size);
69 : void *(*new_gmp_realloc)(void *ptr, size_t old_size, size_t new_size);
70 : void (*new_gmp_free)(void *ptr, size_t old_size);
71 916 : mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
72 916 : if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
73 916 : if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
74 916 : if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
75 916 : mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
76 916 : }
77 :
78 : #define LIMBS(x) ((mp_limb_t *)((x)+2))
79 : #define NLIMBS(x) (lgefint(x)-2)
80 : /*This one is for t_REALs to emphasize they are not t_INTs*/
81 : #define RLIMBS(x) ((mp_limb_t *)((x)+2))
82 : #define RNLIMBS(x) (lg(x)-2)
83 :
84 : INLINE void
85 1874705 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
86 : {
87 1874705 : while (--n >= 0) x[n]=y[n];
88 1874705 : }
89 :
90 : INLINE void
91 204376146 : xmpn_mirror(mp_limb_t *x, long n)
92 : {
93 : long i;
94 1188177846 : for(i=0;i<(n>>1);i++)
95 : {
96 983801700 : ulong m=x[i];
97 983801700 : x[i]=x[n-1-i];
98 983801700 : x[n-1-i]=m;
99 : }
100 204376146 : }
101 :
102 : INLINE void
103 302485109 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
104 : {
105 : long i;
106 2879244422 : for(i=0;i<n;i++)
107 2576759313 : z[i]=x[n-1-i];
108 302485109 : }
109 :
110 : INLINE void
111 63794049 : xmpn_zero(mp_ptr x, mp_size_t n)
112 : {
113 63794049 : while (--n >= 0) x[n]=0;
114 63794049 : }
115 :
116 : INLINE GEN
117 14130960 : icopy_ef(GEN x, long l)
118 : {
119 14130960 : register long lx = lgefint(x);
120 14130960 : const GEN y = cgeti(l);
121 :
122 14130937 : while (--lx > 0) y[lx]=x[lx];
123 14130937 : return y;
124 : }
125 :
126 : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
127 : * GENs but pairs (long *a, long na) representing a list of digits (in basis
128 : * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
129 : * need to reintroduce codewords ]
130 : * Use speci(a,na) to visualize the corresponding GEN.
131 : */
132 :
133 : /***********************************************************************/
134 : /** **/
135 : /** ADDITION / SUBTRACTION **/
136 : /** **/
137 : /***********************************************************************/
138 :
139 : GEN
140 2897370 : setloop(GEN a)
141 : {
142 2897370 : pari_sp av = avma - 2 * sizeof(long);
143 2897370 : (void)cgetg(lgefint(a) + 3, t_VECSMALL);
144 2897373 : return icopy_avma(a, av); /* two cells of extra space after a */
145 : }
146 :
147 : /* we had a = setloop(?), then some incloops. Reset a to b */
148 : GEN
149 170732 : resetloop(GEN a, GEN b) {
150 170732 : a[0] = evaltyp(t_INT) | evallg(lgefint(b));
151 170732 : affii(b, a); return a;
152 : }
153 :
154 : /* assume a > 0, initialized by setloop. Do a++ */
155 : static GEN
156 69122201 : incpos(GEN a)
157 : {
158 69122201 : long i, l = lgefint(a);
159 69122206 : for (i=2; i<l; i++)
160 69122202 : if (++uel(a,i)) return a;
161 4 : a[l] = 1; l++;
162 4 : a[0]=evaltyp(t_INT) | _evallg(l);
163 4 : a[1]=evalsigne(1) | evallgefint(l);
164 4 : return a;
165 : }
166 :
167 : /* assume a < 0, initialized by setloop. Do a++ */
168 : static GEN
169 6268 : incneg(GEN a)
170 : {
171 6268 : long i, l = lgefint(a)-1;
172 6268 : if (uel(a,2)--)
173 : {
174 6264 : if (!a[l]) /* implies l = 2 */
175 : {
176 552 : a[0] = evaltyp(t_INT) | _evallg(2);
177 552 : a[1] = evalsigne(0) | evallgefint(2);
178 : }
179 6264 : return a;
180 : }
181 5 : for (i=3; i<=l; i++)
182 5 : if (uel(a,i)--) break;
183 4 : if (!a[l])
184 : {
185 4 : a[0] = evaltyp(t_INT) | _evallg(l);
186 4 : a[1] = evalsigne(-1) | evallgefint(l);
187 : }
188 4 : return a;
189 : }
190 :
191 : /* assume a initialized by setloop. Do a++ */
192 : GEN
193 69435851 : incloop(GEN a)
194 : {
195 69435851 : switch(signe(a))
196 : {
197 : case 0:
198 331620 : a[0]=evaltyp(t_INT) | _evallg(3);
199 331620 : a[1]=evalsigne(1) | evallgefint(3);
200 331620 : a[2]=1; return a;
201 6268 : case -1: return incneg(a);
202 69097963 : default: return incpos(a);
203 : }
204 : }
205 :
206 : INLINE GEN
207 1162269746 : adduispec(ulong s, GEN x, long nx)
208 : {
209 : GEN zd;
210 : long lz;
211 :
212 1162269746 : if (nx == 1) return adduu(uel(x,0), s);
213 234519245 : lz = nx+3; zd = cgeti(lz);
214 233886327 : if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
215 4916099 : zd[lz-1]=1;
216 : else
217 229065235 : lz--;
218 233981334 : zd[1] = evalsigne(1) | evallgefint(lz);
219 233981334 : return zd;
220 : }
221 :
222 : GEN
223 263051777 : adduispec_offset(ulong s, GEN x, long offset, long nx)
224 : {
225 263051777 : GEN xd=x+2+offset;
226 263051777 : while (nx && *(xd+nx-1)==0) nx--;
227 263051777 : if (!nx) return utoi(s);
228 246550718 : return adduispec(s,xd,nx);
229 : }
230 :
231 : INLINE GEN
232 1637008522 : addiispec(GEN x, GEN y, long nx, long ny)
233 : {
234 : GEN zd;
235 : long lz;
236 :
237 1637008522 : if (nx < ny) swapspec(x,y, nx,ny);
238 1637008522 : if (ny == 1) return adduispec(*y,x,nx);
239 777013296 : lz = nx+3; zd = cgeti(lz);
240 :
241 775072173 : if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
242 14202047 : zd[lz-1]=1;
243 : else
244 762036591 : lz--;
245 :
246 776238638 : zd[1] = evalsigne(1) | evallgefint(lz);
247 776238638 : return zd;
248 : }
249 :
250 : /* assume x >= y */
251 : INLINE GEN
252 829423602 : subiuspec(GEN x, ulong s, long nx)
253 : {
254 : GEN zd;
255 : long lz;
256 :
257 829423602 : if (nx == 1) return utoi(x[0] - s);
258 :
259 96301248 : lz = nx + 2; zd = cgeti(lz);
260 96283617 : mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
261 96284423 : if (! zd[lz - 1]) { --lz; }
262 :
263 96284423 : zd[1] = evalsigne(1) | evallgefint(lz);
264 96284423 : return zd;
265 : }
266 :
267 : /* assume x > y */
268 : INLINE GEN
269 1656994678 : subiispec(GEN x, GEN y, long nx, long ny)
270 : {
271 : GEN zd;
272 : long lz;
273 1656994678 : if (ny==1) return subiuspec(x,*y,nx);
274 850369056 : lz = nx+2; zd = cgeti(lz);
275 :
276 844428491 : mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
277 846314416 : while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
278 :
279 846314416 : zd[1] = evalsigne(1) | evallgefint(lz);
280 846314416 : return zd;
281 : }
282 :
283 : static void
284 307070248 : roundr_up_ip(GEN x, long l)
285 : {
286 307070248 : long i = l;
287 : for(;;)
288 : {
289 308089848 : if (++((ulong*)x)[--i]) break;
290 1127441 : if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
291 1019600 : }
292 307070248 : }
293 :
294 : void
295 120386636 : affir(GEN x, GEN y)
296 : {
297 120386636 : const long s = signe(x), ly = lg(y);
298 : long lx, sh, i;
299 :
300 120386636 : if (!s)
301 : {
302 4661644 : y[1] = evalexpo(-bit_accuracy(ly));
303 4661644 : return;
304 : }
305 115724992 : lx = lgefint(x); sh = bfffo(*int_MSW(x));
306 115724992 : y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
307 115724992 : if (sh) {
308 114421785 : if (lx <= ly)
309 : {
310 72606489 : for (i=lx; i<ly; i++) y[i]=0;
311 72606489 : mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
312 72606489 : xmpn_mirror(LIMBS(y),lx-2);
313 72606489 : return;
314 : }
315 41815296 : mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
316 41815296 : uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
317 41815296 : xmpn_mirror(LIMBS(y),ly-2);
318 : /* lx > ly: round properly */
319 41815296 : if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
320 : }
321 : else {
322 1303207 : GEN xd=int_MSW(x);
323 1303207 : if (lx <= ly)
324 : {
325 709222 : for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
326 709222 : for ( ; i<ly; i++) y[i]=0;
327 709222 : return;
328 : }
329 593985 : for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
330 : /* lx > ly: round properly */
331 593985 : if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
332 : }
333 : }
334 :
335 : INLINE GEN
336 400629705 : shiftispec(GEN x, long nx, long n)
337 : {
338 : long ny,m;
339 : GEN yd, y;
340 :
341 400629705 : if (!n) return icopyspec(x, nx);
342 :
343 392142204 : if (n > 0)
344 : {
345 253774504 : long d = dvmdsBIL(n, &m);
346 : long i;
347 :
348 253742484 : ny = nx + d + (m!=0);
349 253742484 : y = cgeti(ny + 2); yd = y + 2;
350 253505432 : for (i=0; i<d; i++) yd[i] = 0;
351 :
352 253505432 : if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
353 : else
354 : {
355 252711816 : ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
356 252910147 : if (carryd) yd[ny - 1] = carryd;
357 240279972 : else ny--;
358 : }
359 : }
360 : else
361 : {
362 138367700 : long d = dvmdsBIL(-n, &m);
363 :
364 138367624 : ny = nx - d;
365 138367624 : if (ny < 1) return gen_0;
366 137710568 : y = cgeti(ny + 2); yd = y + 2;
367 :
368 137709896 : if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
369 : else
370 : {
371 136709350 : mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
372 136709414 : if (yd[ny - 1] == 0)
373 : {
374 15462024 : if (ny == 1) { avma = (pari_sp)(yd + 1); return gen_0; }
375 14170661 : ny--;
376 : }
377 : }
378 : }
379 390086685 : y[1] = evalsigne(1)|evallgefint(ny + 2);
380 390086685 : return y;
381 : }
382 :
383 : GEN
384 25930787 : mantissa2nr(GEN x, long n)
385 : {
386 25930787 : long ly, i, m, s = signe(x), lx = lg(x);
387 : GEN y;
388 25930787 : if (!s) return gen_0;
389 25929479 : if (!n)
390 : {
391 9391253 : y = cgeti(lx);
392 9391253 : y[1] = evalsigne(s) | evallgefint(lx);
393 9391253 : xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
394 9391253 : return y;
395 : }
396 16538226 : if (n > 0)
397 : {
398 314414 : GEN z = (GEN)avma;
399 314414 : long d = dvmdsBIL(n, &m);
400 :
401 314414 : ly = lx+d; y = new_chunk(ly);
402 314414 : for ( ; d; d--) *--z = 0;
403 314414 : if (!m) for (i=2; i<lx; i++) y[i]=x[i];
404 : else
405 : {
406 314019 : register const ulong sh = BITS_IN_LONG - m;
407 314019 : shift_left(y,x, 2,lx-1, 0,m);
408 314019 : i = uel(x,2) >> sh;
409 : /* Extend y on the left? */
410 314019 : if (i) { ly++; y = new_chunk(1); y[2] = i; }
411 : }
412 : }
413 : else
414 : {
415 16223812 : ly = lx - dvmdsBIL(-n, &m);
416 16223812 : if (ly<3) return gen_0;
417 16223812 : y = new_chunk(ly);
418 16223812 : if (m) {
419 16181487 : shift_right(y,x, 2,ly, 0,m);
420 16181487 : if (y[2] == 0)
421 : {
422 0 : if (ly==3) { avma = (pari_sp)(y+3); return gen_0; }
423 0 : ly--; avma = (pari_sp)(++y);
424 : }
425 : } else {
426 42325 : for (i=2; i<ly; i++) y[i]=x[i];
427 : }
428 : }
429 16538226 : xmpn_mirror(LIMBS(y),ly-2);
430 16538226 : y[1] = evalsigne(s)|evallgefint(ly);
431 16538226 : y[0] = evaltyp(t_INT)|evallg(ly); return y;
432 : }
433 :
434 : GEN
435 1745514 : truncr(GEN x)
436 : {
437 : long s, e, d, m, i;
438 : GEN y;
439 1745514 : if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
440 1531204 : d = nbits2lg(e+1); m = remsBIL(e);
441 1531204 : if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
442 :
443 1531204 : y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
444 1531204 : if (++m == BITS_IN_LONG)
445 152 : for (i=2; i<d; i++) y[d-i+1]=x[i];
446 : else
447 : {
448 1531052 : GEN z=cgeti(d);
449 1531052 : for (i=2; i<d; i++) z[d-i+1]=x[i];
450 1531052 : mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
451 1531052 : avma=(pari_sp)y;
452 : }
453 1531204 : return y;
454 : }
455 :
456 : /* integral part */
457 : GEN
458 512194 : floorr(GEN x)
459 : {
460 : long e, d, m, i, lx;
461 : GEN y;
462 512194 : if (signe(x) >= 0) return truncr(x);
463 209376 : if ((e=expo(x)) < 0) return gen_m1;
464 78969 : d = nbits2lg(e+1); m = remsBIL(e);
465 78969 : lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
466 78969 : y = cgeti(d+1);
467 78969 : if (++m == BITS_IN_LONG)
468 : {
469 40 : for (i=2; i<d; i++) y[d-i+1]=x[i];
470 40 : i=d; while (i<lx && !x[i]) i++;
471 40 : if (i==lx) goto END;
472 : }
473 : else
474 : {
475 78929 : GEN z=cgeti(d);
476 78929 : for (i=2; i<d; i++) z[d-i+1]=x[i];
477 78929 : mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
478 78929 : if (uel(x,d-1)<<m == 0)
479 : {
480 17010 : i=d; while (i<lx && !x[i]) i++;
481 17010 : if (i==lx) goto END;
482 : }
483 : }
484 71699 : if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
485 0 : y[d++]=1;
486 : END:
487 78969 : y[1] = evalsigne(-1) | evallgefint(d);
488 78969 : return y;
489 : }
490 :
491 : INLINE int
492 2065666237 : cmpiispec(GEN x, GEN y, long lx, long ly)
493 : {
494 2065666237 : if (lx < ly) return -1;
495 1977509454 : if (lx > ly) return 1;
496 1882670929 : return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
497 : }
498 :
499 : INLINE int
500 112309703 : equaliispec(GEN x, GEN y, long lx, long ly)
501 : {
502 112309703 : if (lx != ly) return 0;
503 112307194 : return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
504 : }
505 :
506 : /***********************************************************************/
507 : /** **/
508 : /** MULTIPLICATION **/
509 : /** **/
510 : /***********************************************************************/
511 : /* assume ny > 0 */
512 : INLINE GEN
513 2293514785 : muluispec(ulong x, GEN y, long ny)
514 : {
515 2293514785 : if (ny == 1)
516 1785712716 : return muluu(x, *y);
517 : else
518 : {
519 507802069 : long lz = ny+3;
520 507802069 : GEN z = cgeti(lz);
521 506414287 : ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
522 506776820 : if (hi) { z[lz - 1] = hi; } else lz--;
523 506776820 : z[1] = evalsigne(1) | evallgefint(lz);
524 506776820 : return z;
525 : }
526 : }
527 :
528 : /* a + b*|y| */
529 : GEN
530 1012036 : addumului(ulong a, ulong b, GEN y)
531 : {
532 : GEN z;
533 : long i, lz;
534 : ulong hi;
535 1012036 : if (!b || !signe(y)) return utoi(a);
536 1011954 : lz = lgefint(y)+1;
537 1011954 : z = cgeti(lz);
538 1011954 : z[2]=a;
539 1011954 : for(i=3;i<lz;i++) z[i]=0;
540 1011954 : hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
541 1011954 : if (hi) z[lz-1]=hi; else lz--;
542 1011954 : z[1] = evalsigne(1) | evallgefint(lz);
543 1011954 : avma=(pari_sp)z; return z;
544 : }
545 :
546 : /***********************************************************************/
547 : /** **/
548 : /** DIVISION **/
549 : /** **/
550 : /***********************************************************************/
551 :
552 : ulong
553 928848223 : umodiu(GEN y, ulong x)
554 : {
555 928848223 : long sy=signe(y);
556 : ulong hi;
557 928848223 : if (!x) pari_err_INV("umodiu",gen_0);
558 928882078 : if (!sy) return 0;
559 796113463 : hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
560 796142970 : if (!hi) return 0;
561 786932471 : return (sy > 0)? hi: x - hi;
562 : }
563 :
564 : /* return |y| \/ x */
565 : GEN
566 127843470 : diviu_rem(GEN y, ulong x, ulong *rem)
567 : {
568 : long ly;
569 : GEN z;
570 :
571 127843470 : if (!x) pari_err_INV("diviu_rem",gen_0);
572 127843497 : if (!signe(y)) { *rem = 0; return gen_0; }
573 :
574 127295908 : ly = lgefint(y);
575 127295908 : if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
576 :
577 127228590 : z = cgeti(ly);
578 127200338 : *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
579 127230314 : if (z [ly - 1] == 0) ly--;
580 127230314 : z[1] = evallgefint(ly) | evalsigne(1);
581 127230314 : return z;
582 : }
583 :
584 : GEN
585 42745813 : divis_rem(GEN y, long x, long *rem)
586 : {
587 42745813 : long sy=signe(y),ly,s;
588 : GEN z;
589 :
590 42745813 : if (!x) pari_err_INV("divis_rem",gen_0);
591 42745813 : if (!sy) { *rem = 0; return gen_0; }
592 31693163 : if (x<0) { s = -sy; x = -x; } else s = sy;
593 :
594 31693163 : ly = lgefint(y);
595 31693163 : if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
596 :
597 20707306 : z = cgeti(ly);
598 20707306 : *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
599 20707306 : if (sy<0) *rem = - *rem;
600 20707306 : if (z[ly - 1] == 0) ly--;
601 20707306 : z[1] = evallgefint(ly) | evalsigne(s);
602 20707306 : return z;
603 : }
604 :
605 : GEN
606 416827 : divis(GEN y, long x)
607 : {
608 416827 : long sy=signe(y),ly,s;
609 : GEN z;
610 :
611 416827 : if (!x) pari_err_INV("divis",gen_0);
612 416827 : if (!sy) return gen_0;
613 416827 : if (x<0) { s = -sy; x = -x; } else s=sy;
614 :
615 416827 : ly = lgefint(y);
616 416827 : if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
617 :
618 416827 : z = cgeti(ly);
619 416827 : (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
620 416827 : if (z[ly - 1] == 0) ly--;
621 416827 : z[1] = evallgefint(ly) | evalsigne(s);
622 416827 : return z;
623 : }
624 :
625 : /* We keep llx bits of x and lly bits of y*/
626 : static GEN
627 36572152 : divrr_with_gmp(GEN x, GEN y)
628 : {
629 36572152 : long lx=RNLIMBS(x),ly=RNLIMBS(y);
630 36572152 : long lw=minss(lx,ly);
631 36572152 : long lly=minss(lw+1,ly);
632 36572152 : GEN w=cgetr(lw+2);
633 36572152 : long lu=lw+lly;
634 36572152 : long llx=minss(lu,lx);
635 36572152 : mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
636 36572152 : mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
637 : mp_limb_t *q,*r;
638 36572152 : long e=expo(x)-expo(y);
639 36572152 : long sx=signe(x),sy=signe(y);
640 36572152 : xmpn_mirrorcopy(z,RLIMBS(y),lly);
641 36572152 : xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
642 36572152 : xmpn_zero(u,lu-llx);
643 36572152 : q = (mp_limb_t *)new_chunk(lw+1);
644 36572152 : r = (mp_limb_t *)new_chunk(lly);
645 :
646 36572152 : mpn_tdiv_qr(q,r,0,u,lu,z,lly);
647 :
648 : /*Round up: This is not exactly correct we should test 2*r>z*/
649 36572152 : if (uel(r,lly-1) > (uel(z,lly-1)>>1))
650 18081973 : mpn_add_1(q,q,lw+1,1);
651 :
652 36572152 : xmpn_mirrorcopy(RLIMBS(w),q,lw);
653 :
654 36572152 : if (q[lw] == 0) e--;
655 14288860 : else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
656 18 : else { w[2] = HIGHBIT; e++; }
657 36572152 : if (sy < 0) sx = -sx;
658 36572152 : w[1] = evalsigne(sx) | evalexpo(e);
659 36572152 : avma=(pari_sp) w;
660 36572152 : return w;
661 : }
662 :
663 : /* We keep llx bits of x and lly bits of y*/
664 : static GEN
665 3388807 : divri_with_gmp(GEN x, GEN y)
666 : {
667 3388807 : long llx=RNLIMBS(x),ly=NLIMBS(y);
668 3388807 : long lly=minss(llx+1,ly);
669 3388807 : GEN w=cgetr(llx+2);
670 3388807 : long lu=llx+lly, ld=ly-lly;
671 3388807 : mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
672 3388807 : mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
673 : mp_limb_t *q,*r;
674 3388807 : long sh=bfffo(y[ly+1]);
675 3388807 : long e=expo(x)-expi(y);
676 3388807 : long sx=signe(x),sy=signe(y);
677 3388807 : if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
678 80531 : else xmpn_copy(z,LIMBS(y)+ld,lly);
679 3388807 : xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
680 3388807 : xmpn_zero(u,lu-llx);
681 3388807 : q = (mp_limb_t *)new_chunk(llx+1);
682 3388807 : r = (mp_limb_t *)new_chunk(lly);
683 :
684 3388807 : mpn_tdiv_qr(q,r,0,u,lu,z,lly);
685 :
686 : /*Round up: This is not exactly correct we should test 2*r>z*/
687 3388807 : if (uel(r,lly-1) > (uel(z,lly-1)>>1))
688 1434111 : mpn_add_1(q,q,llx+1,1);
689 :
690 3388807 : xmpn_mirrorcopy(RLIMBS(w),q,llx);
691 :
692 3388807 : if (q[llx] == 0) e--;
693 2849756 : else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
694 0 : else { w[2] = HIGHBIT; e++; }
695 3388807 : if (sy < 0) sx = -sx;
696 3388807 : w[1] = evalsigne(sx) | evalexpo(e);
697 3388807 : avma=(pari_sp) w;
698 3388807 : return w;
699 : }
700 :
701 : GEN
702 19022350 : divri(GEN x, GEN y)
703 : {
704 19022350 : long s = signe(y);
705 :
706 19022350 : if (!s) pari_err_INV("divri",gen_0);
707 19022350 : if (!signe(x)) return real_0_bit(expo(x) - expi(y));
708 18890704 : if (!is_bigint(y)) {
709 15501897 : GEN z = divru(x, y[2]);
710 15501897 : if (s < 0) togglesign(z);
711 15501897 : return z;
712 : }
713 3388807 : return divri_with_gmp(x,y);
714 : }
715 :
716 : GEN
717 90721310 : divrr(GEN x, GEN y)
718 : {
719 90721310 : long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
720 : ulong y0,y1;
721 : GEN r, r1;
722 :
723 90721310 : if (!sy) pari_err_INV("divrr",y);
724 90721310 : e = expo(x) - expo(y);
725 90721310 : if (!sx) return real_0_bit(e);
726 89150659 : if (sy<0) sx = -sx;
727 :
728 89150659 : lx=lg(x); ly=lg(y);
729 89150659 : if (ly==3)
730 : {
731 52578507 : ulong k = x[2], l = (lx>3)? x[3]: 0;
732 : LOCAL_HIREMAINDER;
733 52578507 : if (k < uel(y,2)) e--;
734 : else
735 : {
736 24237438 : l >>= 1; if (k&1) l |= HIGHBIT;
737 24237438 : k >>= 1;
738 : }
739 52578507 : hiremainder = k; k = divll(l,y[2]);
740 52578507 : if (hiremainder & HIGHBIT)
741 : {
742 14549346 : k++;
743 14549346 : if (!k) { k = HIGHBIT; e++; }
744 : }
745 52578507 : r = cgetr(3);
746 52578507 : r[1] = evalsigne(sx) | evalexpo(e);
747 52578507 : r[2] = k; return r;
748 : }
749 :
750 36572152 : if (ly>=DIVRR_GMP_LIMIT)
751 36572152 : return divrr_with_gmp(x,y);
752 :
753 0 : lr = minss(lx,ly); r = new_chunk(lr);
754 0 : r1 = r-1;
755 0 : r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
756 0 : r1[lr] = (lx>ly)? x[lr]: 0;
757 0 : y0 = y[2]; y1 = y[3];
758 0 : for (i=0; i<lr-1; i++)
759 : { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
760 : ulong k, qp;
761 : LOCAL_HIREMAINDER;
762 : LOCAL_OVERFLOW;
763 :
764 0 : if (uel(r1,1) == y0)
765 : {
766 0 : qp = ULONG_MAX; k = addll(y0,r1[2]);
767 : }
768 : else
769 : {
770 0 : 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 0 : hiremainder = r1[1]; overflow = 0;
778 0 : qp = divll(r1[2],y0); k = hiremainder;
779 : }
780 0 : j = lr-i+1;
781 0 : if (!overflow)
782 : {
783 : long k3, k4;
784 0 : k3 = mulll(qp,y1);
785 0 : if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
786 0 : k4 = subll(hiremainder,k);
787 : else
788 : {
789 0 : k3 = subll(k3, r1[3]);
790 0 : k4 = subllx(hiremainder,k);
791 : }
792 0 : while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
793 : }
794 0 : if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
795 0 : for (j--; j>1; j--)
796 : {
797 0 : r1[j] = subll(r1[j], addmul(qp,y[j]));
798 0 : hiremainder += overflow;
799 : }
800 0 : if (uel(r1,1) != hiremainder)
801 : {
802 0 : if (uel(r1,1) < hiremainder)
803 : {
804 0 : qp--;
805 0 : j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
806 0 : 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 0 : *++r1 = qp;
821 : }
822 : /* i = lr-1 */
823 : /* round correctly */
824 0 : if (uel(r1,1) > (y0>>1))
825 : {
826 0 : j=i; do uel(r,--j)++; while (j && !r[j]);
827 : }
828 0 : r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
829 0 : if (r[0] == 0) e--;
830 0 : else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
831 : else { /* possible only when rounding up to 0x2 0x0 ... */
832 0 : r[2] = (long)HIGHBIT; e++;
833 : }
834 0 : r[0] = evaltyp(t_REAL)|evallg(lr);
835 0 : r[1] = evalsigne(sx) | evalexpo(e);
836 0 : 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 thus can be disposed of with cgiv
843 : * instead of gerepile)
844 : * If *z is zero, we put gen_0 here and no copy.
845 : * space needed: lx + ly
846 : */
847 : GEN
848 974053169 : dvmdii(GEN x, GEN y, GEN *z)
849 : {
850 974053169 : long sx=signe(x),sy=signe(y);
851 : long lx, ly, lq;
852 : pari_sp av;
853 : GEN r,q;
854 :
855 974053169 : if (!sy) pari_err_INV("dvmdii",y);
856 975053692 : if (!sx)
857 : {
858 20695888 : if (!z || z == ONLY_REM) return gen_0;
859 8017010 : *z=gen_0; return gen_0;
860 : }
861 954357804 : lx=lgefint(x);
862 954357804 : ly=lgefint(y); lq=lx-ly;
863 954357804 : if (lq <= 0)
864 : {
865 514037394 : if (lq == 0)
866 : {
867 491013776 : long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
868 491482604 : if (s>0) goto DIVIDE;
869 179348269 : if (s==0)
870 : {
871 7887705 : if (z == ONLY_REM) return gen_0;
872 3857929 : if (z) *z = gen_0;
873 3857929 : if (sx < 0) sy = -sy;
874 3857929 : return stoi(sy);
875 : }
876 : }
877 194484182 : if (z == ONLY_REM) return icopy(x);
878 5744380 : if (z) *z = icopy(x);
879 5744380 : return gen_0;
880 : }
881 : DIVIDE: /* quotient is non-zero */
882 752454745 : av=avma; if (sx<0) sy = -sy;
883 752454745 : if (ly==3)
884 : {
885 257972494 : ulong lq = lx;
886 : ulong si;
887 257972494 : q = cgeti(lq);
888 257848723 : si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
889 258004188 : if (q[lq - 1] == 0) lq--;
890 258004188 : if (z == ONLY_REM)
891 : {
892 217031386 : avma=av; if (!si) return gen_0;
893 191867708 : r=cgeti(3);
894 191755680 : r[1] = evalsigne(sx) | evallgefint(3);
895 191755680 : r[2] = si; return r;
896 : }
897 40972802 : q[1] = evalsigne(sy) | evallgefint(lq);
898 40972802 : if (!z) return q;
899 37403665 : if (!si) { *z=gen_0; return q; }
900 15290151 : r=cgeti(3);
901 15290139 : r[1] = evalsigne(sx) | evallgefint(3);
902 15290139 : r[2] = si; *z=r; return q;
903 : }
904 494482251 : if (z == ONLY_REM)
905 : {
906 488720844 : ulong lr = lgefint(y);
907 488720844 : ulong lq = lgefint(x)-lgefint(y)+3;
908 488720844 : GEN r = cgeti(lr);
909 481714674 : GEN q = cgeti(lq);
910 479291176 : mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
911 489988207 : if (!r[lr - 1])
912 : {
913 155331135 : while(lr>2 && !r[lr - 1]) lr--;
914 155331135 : if (lr == 2) {avma=av; return gen_0;} /* exact division */
915 : }
916 485397258 : r[1] = evalsigne(sx) | evallgefint(lr);
917 485397258 : avma = (pari_sp) r; return r;
918 : }
919 : else
920 : {
921 5761407 : ulong lq = lgefint(x)-lgefint(y)+3;
922 5761407 : ulong lr = lgefint(y);
923 5761407 : GEN q = cgeti(lq);
924 5761407 : GEN r = cgeti(lr);
925 5761407 : mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
926 5761407 : if (q[lq - 1] == 0) lq--;
927 5761407 : q[1] = evalsigne(sy) | evallgefint(lq);
928 5761407 : if (!z) { avma = (pari_sp)q; return q; }
929 5137554 : if (!r[lr - 1])
930 : {
931 1828941 : while(lr>2 && !r[lr - 1]) lr--;
932 1828941 : if (lr == 2) {avma=(pari_sp) q; *z=gen_0; return q;} /* exact division */
933 : }
934 3760927 : r[1] = evalsigne(sx) | evallgefint(lr);
935 3760927 : avma = (pari_sp) r; *z = r; return q;
936 : }
937 : }
938 :
939 : /* Montgomery reduction.
940 : * N has k words, assume T >= 0 has less than 2k.
941 : * Return res := T / B^k mod N, where B = 2^BIL
942 : * such that 0 <= res < T/B^k + N and res has less than k words */
943 : GEN
944 13205886 : red_montgomery(GEN T, GEN N, ulong inv)
945 : {
946 : pari_sp av;
947 : GEN Te, Td, Ne, Nd, scratch;
948 13205886 : ulong i, j, m, t, d, k = NLIMBS(N);
949 : int carry;
950 : LOCAL_HIREMAINDER;
951 : LOCAL_OVERFLOW;
952 :
953 13205886 : if (k == 0) return gen_0;
954 13205886 : d = NLIMBS(T); /* <= 2*k */
955 13205886 : if (d == 0) return gen_0;
956 : #ifdef DEBUG
957 : if (d > 2*k) pari_err_BUG("red_montgomery");
958 : #endif
959 13205877 : if (k == 1)
960 : { /* as below, special cased for efficiency */
961 492 : ulong n = uel(N,2);
962 492 : if (d == 1) {
963 492 : hiremainder = uel(T,2);
964 492 : m = hiremainder * inv;
965 492 : (void)addmul(m, n); /* t + m*n = 0 */
966 492 : return utoi(hiremainder);
967 : } else { /* d = 2 */
968 0 : hiremainder = uel(T,2);
969 0 : m = hiremainder * inv;
970 0 : (void)addmul(m, n); /* t + m*n = 0 */
971 0 : t = addll(hiremainder, uel(T,3));
972 0 : if (overflow) t -= n; /* t > n doesn't fit in 1 word */
973 0 : return utoi(t);
974 : }
975 : }
976 : /* assume k >= 2 */
977 13205385 : av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
978 :
979 : /* copy T to scratch space (pad with zeroes to 2k words) */
980 13219926 : Td = scratch;
981 13219926 : Te = T + 2;
982 13219926 : for (i=0; i < d ; i++) *Td++ = *Te++;
983 13219926 : for ( ; i < (k<<1); i++) *Td++ = 0;
984 :
985 13219926 : Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
986 13219926 : Ne = N + 1; /* 1 beyond end of N mantissa */
987 :
988 13219926 : carry = 0;
989 100012035 : for (i=0; i<k; i++) /* set T := T/B nod N, k times */
990 : {
991 86792109 : Td = Te; /* one beyond end of (new) T mantissa */
992 86792109 : Nd = Ne;
993 86792109 : hiremainder = *++Td;
994 86792109 : m = hiremainder * inv; /* solve T + m N = O(B) */
995 :
996 : /* set T := (T + mN) / B */
997 86792109 : Te = Td;
998 86792109 : (void)addmul(m, *++Nd); /* = 0 */
999 809214084 : for (j=1; j<k; j++)
1000 : {
1001 722421975 : t = addll(addmul(m, *++Nd), *++Td);
1002 722421975 : *Td = t;
1003 722421975 : hiremainder += overflow;
1004 : }
1005 86792109 : t = addll(hiremainder, *++Td); *Td = t + carry;
1006 86792109 : carry = (overflow || (carry && *Td == 0));
1007 : }
1008 13219926 : if (carry)
1009 : { /* Td > N overflows (k+1 words), set Td := Td - N */
1010 34743 : Td = Te;
1011 34743 : Nd = Ne;
1012 34743 : t = subll(*++Td, *++Nd); *Td = t;
1013 34743 : while (Td < (GEN)av) { t = subllx(*++Td, *++Nd); *Td = t; }
1014 : }
1015 :
1016 : /* copy result */
1017 13219926 : Td = (GEN)av - 1; /* *Td = high word of final result */
1018 13219926 : while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
1019 13219926 : k = Td - Te; if (!k) { avma = av; return gen_0; }
1020 13219926 : Td = (GEN)av - k; /* will write mantissa there */
1021 13219926 : (void)memmove(Td, Te+1, k*sizeof(long));
1022 13219926 : Td -= 2;
1023 13219926 : Td[0] = evaltyp(t_INT) | evallg(k+2);
1024 13201365 : 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 13201365 : avma = (pari_sp)Td; return Td;
1036 : }
1037 :
1038 : /* EXACT INTEGER DIVISION */
1039 :
1040 : #if 1 && !defined(_WIN64) /* use undocumented GMP interface */
1041 : /* mpz_divexact_ui is not LLP64 friendly */
1042 : static void
1043 276507922 : GEN2mpz(mpz_t X, GEN x)
1044 : {
1045 276507922 : long l = lgefint(x)-2;
1046 276507922 : X->_mp_alloc = l;
1047 276507922 : X->_mp_size = signe(x) > 0? l: -l;
1048 276507922 : X->_mp_d = LIMBS(x);
1049 276507922 : }
1050 : static void
1051 264999449 : mpz2GEN(GEN z, mpz_t Z)
1052 : {
1053 264999449 : long l = Z->_mp_size;
1054 264999449 : z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
1055 264999449 : }
1056 :
1057 : /* Find z such that x=y*z, knowing that y | x (unchecked) */
1058 : GEN
1059 273317446 : diviiexact(GEN x, GEN y)
1060 : {
1061 273317446 : if (!signe(y)) pari_err_INV("diviiexact",y);
1062 273317652 : if (lgefint(y) == 3)
1063 : {
1064 258671388 : GEN z = diviuexact(x, y[2]);
1065 258671066 : if (signe(y) < 0) togglesign(z);
1066 258671513 : return z;
1067 : }
1068 14646264 : if (!signe(x)) return gen_0;
1069 : {
1070 11518353 : long l = lgefint(x);
1071 : mpz_t X, Y, Z;
1072 11518353 : GEN z = cgeti(l);
1073 11518353 : GEN2mpz(X, x);
1074 11518353 : GEN2mpz(Y, y);
1075 11518353 : Z->_mp_alloc = l-2;
1076 11518353 : Z->_mp_size = l-2;
1077 11518353 : Z->_mp_d = LIMBS(z);
1078 11518353 : mpz_divexact(Z, X, Y);
1079 11518353 : mpz2GEN(z, Z); return z;
1080 : }
1081 : }
1082 : #else
1083 : /* Find z such that x=y*z, knowing that y | x (unchecked)
1084 : * Method: y0 z0 = x0 mod B = 2^BITS_IN_LONG ==> z0 = 1/y0 mod B.
1085 : * Set x := (x - z0 y) / B, updating only relevant words, and repeat */
1086 : GEN
1087 : diviiexact(GEN x, GEN y)
1088 : {
1089 : /*TODO: use mpn_bdivmod instead*/
1090 : if (!signe(y)) pari_err_INV("diviiexact",y);
1091 : return divii(x,y);
1092 : }
1093 : #endif
1094 :
1095 : #ifdef mpn_divexact_1
1096 : GEN
1097 : diviuexact(GEN y, ulong x)
1098 : {
1099 : long ly;
1100 : GEN z;
1101 : if (!signe(y)) return gen_0;
1102 : ly = lgefint(y);
1103 : z = cgeti(ly);
1104 : mpn_divexact_1(LIMBS(z), LIMBS(y), NLIMBS(y), x);
1105 : if (z [ly - 1] == 0) ly--;
1106 : z[1] = evallgefint(ly) | evalsigne(signe(y));
1107 : return z;
1108 : }
1109 : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly */
1110 : /* assume y != 0 and the division is exact */
1111 : GEN
1112 291488286 : diviuexact(GEN x, ulong y)
1113 : {
1114 291488286 : if (!signe(x)) return gen_0;
1115 : {
1116 253479968 : long l = lgefint(x);
1117 : mpz_t X, Z;
1118 253479968 : GEN z = cgeti(l);
1119 253466744 : GEN2mpz(X, x);
1120 253467407 : Z->_mp_alloc = l-2;
1121 253467407 : Z->_mp_size = l-2;
1122 253467407 : Z->_mp_d = LIMBS(z);
1123 253467407 : mpz_divexact_ui(Z, X, y);
1124 253479333 : mpz2GEN(z, Z); return z;
1125 : }
1126 : }
1127 : #else
1128 : /* assume y != 0 and the division is exact */
1129 : GEN
1130 : diviuexact(GEN x, ulong y)
1131 : {
1132 : /*TODO: implement true exact division.*/
1133 : return divii(x,utoi(y));
1134 : }
1135 : #endif
1136 :
1137 : /* assume yz != and yz | x */
1138 : GEN
1139 544603 : diviuuexact(GEN x, ulong y, ulong z)
1140 : {
1141 : long tmp[4];
1142 : ulong t;
1143 : LOCAL_HIREMAINDER;
1144 544603 : t = mulll(y, z);
1145 544603 : if (!hiremainder) return diviuexact(x, t);
1146 0 : tmp[0] = evaltyp(t_INT)|_evallg(4);
1147 0 : tmp[1] = evalsigne(1)|evallgefint(4);
1148 0 : tmp[2] = t;
1149 0 : tmp[3] = hiremainder;
1150 0 : return diviiexact(x, tmp);
1151 : }
1152 :
1153 :
1154 : /********************************************************************/
1155 : /** **/
1156 : /** INTEGER MULTIPLICATION **/
1157 : /** **/
1158 : /********************************************************************/
1159 :
1160 : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
1161 : GEN
1162 2322853032 : muliispec(GEN x, GEN y, long nx, long ny)
1163 : {
1164 : GEN zd;
1165 : long lz;
1166 : ulong hi;
1167 :
1168 2322853032 : if (nx < ny) swapspec(x,y, nx,ny);
1169 2322853032 : if (!ny) return gen_0;
1170 2322853032 : if (ny == 1) return muluispec((ulong)*y, x, nx);
1171 :
1172 678536775 : lz = nx+ny+2;
1173 678536775 : zd = cgeti(lz);
1174 676523435 : hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
1175 679847840 : if (!hi) lz--;
1176 : /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
1177 :
1178 679847840 : zd[1] = evalsigne(1) | evallgefint(lz);
1179 679847840 : return zd;
1180 : }
1181 : GEN
1182 566811 : muluui(ulong x, ulong y, GEN z)
1183 : {
1184 566811 : long t, s = signe(z);
1185 : GEN r;
1186 : LOCAL_HIREMAINDER;
1187 :
1188 566811 : if (!x || !y || !signe(z)) return gen_0;
1189 566415 : t = mulll(x,y);
1190 566415 : if (!hiremainder)
1191 566415 : r = muluispec(t, z+2, lgefint(z)-2);
1192 : else
1193 : {
1194 : long tmp[2];
1195 0 : tmp[1] = hiremainder;
1196 0 : tmp[0] = t;
1197 0 : r = muliispec(z+2,tmp, lgefint(z)-2, 2);
1198 : }
1199 566415 : setsigne(r,s); return r;
1200 : }
1201 :
1202 : GEN
1203 788037787 : sqrispec(GEN x, long nx)
1204 : {
1205 : GEN zd;
1206 : long lz;
1207 :
1208 788037787 : if (!nx) return gen_0;
1209 263504804 : if (nx==1) return sqru(*x);
1210 :
1211 194243780 : lz = (nx<<1)+2;
1212 194243780 : zd = cgeti(lz);
1213 : #ifdef mpn_sqr
1214 191771461 : mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
1215 : #else
1216 : mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
1217 : #endif
1218 195434937 : if (zd[lz-1]==0) lz--;
1219 :
1220 195434937 : zd[1] = evalsigne(1) | evallgefint(lz);
1221 195434937 : return zd;
1222 : }
1223 :
1224 : INLINE GEN
1225 5864762 : sqrispec_mirror(GEN x, long nx)
1226 : {
1227 5864762 : GEN cx=new_chunk(nx);
1228 : GEN z;
1229 5864762 : xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
1230 5864762 : z=sqrispec(cx, nx);
1231 5864762 : xmpn_mirror(LIMBS(z), NLIMBS(z));
1232 5864762 : return z;
1233 : }
1234 :
1235 : /* leaves garbage on the stack. */
1236 : INLINE GEN
1237 67551373 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
1238 : {
1239 : GEN cx, cy, z;
1240 67551373 : long s = 0;
1241 67551373 : while (nx && x[nx-1]==0) { nx--; s++; }
1242 67551373 : while (ny && y[ny-1]==0) { ny--; s++; }
1243 67551373 : cx=new_chunk(nx); cy=new_chunk(ny);
1244 67551373 : xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
1245 67551373 : xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
1246 67551373 : z = nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
1247 67551373 : if (s)
1248 : {
1249 5452099 : long i, lz = lgefint(z) + s;
1250 5452099 : (void)new_chunk(s);
1251 5452099 : z -= s;
1252 5452099 : for (i=0; i<s; i++) z[2+i]=0;
1253 5452099 : z[1] = evalsigne(1) | evallgefint(lz);
1254 5452099 : z[0] = evaltyp(t_INT) | evallg(lz);
1255 : }
1256 67551373 : xmpn_mirror(LIMBS(z), NLIMBS(z));
1257 67551373 : return z;
1258 : }
1259 :
1260 : /* x % (2^n), assuming n >= 0 */
1261 : GEN
1262 17518273 : remi2n(GEN x, long n)
1263 : {
1264 : ulong hi;
1265 17518273 : long l, k, lx, ly, sx = signe(x);
1266 : GEN z, xd, zd;
1267 :
1268 17518273 : if (!sx || !n) return gen_0;
1269 :
1270 17155868 : k = dvmdsBIL(n, &l);
1271 17174716 : lx = lgefint(x);
1272 17174716 : if (lx < k+3) return icopy(x);
1273 :
1274 16832580 : xd = x + (2 + k);
1275 : /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
1276 : * ^--- initial xd */
1277 16832580 : hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
1278 16832580 : if (!hi)
1279 : { /* strip leading zeroes from result */
1280 2188157 : xd--; while (k && !*xd) { k--; xd--; }
1281 2188157 : if (!k) return gen_0;
1282 1146760 : ly = k+2;
1283 : }
1284 : else
1285 14644423 : ly = k+3;
1286 :
1287 15791183 : zd = z = cgeti(ly);
1288 15649165 : *++zd = evalsigne(sx) | evallgefint(ly);
1289 15649165 : xd = x+1; for ( ;k; k--) *++zd = *++xd;
1290 15649165 : if (hi) *++zd = hi;
1291 15649165 : return z;
1292 : }
1293 :
1294 : /********************************************************************/
1295 : /** **/
1296 : /** INTEGER SQUARE ROOT **/
1297 : /** **/
1298 : /********************************************************************/
1299 :
1300 : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
1301 : * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
1302 : * remainder is 0. R = NULL is allowed. */
1303 : GEN
1304 3555622 : sqrtremi(GEN a, GEN *r)
1305 : {
1306 3555622 : long l, na = NLIMBS(a);
1307 : mp_size_t nr;
1308 : GEN S;
1309 3555622 : if (!na) {
1310 4 : if (r) *r = gen_0;
1311 4 : return gen_0;
1312 : }
1313 3555618 : l = (na + 5) >> 1; /* 2 + ceil(na/2) */
1314 3555618 : S = cgetipos(l);
1315 3555618 : if (r) {
1316 348142 : GEN R = cgeti(2 + na);
1317 348142 : nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
1318 348142 : if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
1319 10400 : else { avma = (pari_sp)S; R = gen_0; }
1320 348142 : *r = R;
1321 : }
1322 : else
1323 3207476 : (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
1324 3555618 : return S;
1325 : }
1326 :
1327 : /* compute sqrt(|a|), assuming a != 0 */
1328 : GEN
1329 23833090 : sqrtr_abs(GEN a)
1330 : {
1331 : GEN res;
1332 : mp_limb_t *b, *c;
1333 23833090 : long l = RNLIMBS(a), e = expo(a), er = e>>1;
1334 : long n;
1335 23833090 : res = cgetr(2 + l);
1336 23833090 : res[1] = evalsigne(1) | evalexpo(er);
1337 23833090 : if (e&1)
1338 : {
1339 11799188 : b = (mp_limb_t *) new_chunk(l<<1);
1340 11799188 : xmpn_zero(b,l);
1341 11799188 : xmpn_mirrorcopy(b+l, RLIMBS(a), l);
1342 11799188 : c = (mp_limb_t *) new_chunk(l);
1343 11799188 : n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
1344 11799188 : if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
1345 : }
1346 : else
1347 : {
1348 : ulong u;
1349 12033902 : b = (mp_limb_t *) mantissa2nr(a,-1);
1350 12033902 : b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
1351 12033902 : b = (mp_limb_t *) new_chunk(l);
1352 12033902 : xmpn_zero(b,l+1); /* overwrites the former b[0] */
1353 12033902 : c = (mp_limb_t *) new_chunk(l + 1);
1354 12033902 : n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
1355 12033902 : u = (ulong)*c++;
1356 12033902 : if ( u&HIGHBIT || (u == ~HIGHBIT &&
1357 0 : (n>l || (n==l && mpn_cmp(b,c,l)>0))))
1358 5863205 : mpn_add_1(c,c,l,1);
1359 : }
1360 23833090 : xmpn_mirrorcopy(RLIMBS(res),c,l);
1361 23833090 : avma = (pari_sp)res; return res;
1362 : }
1363 :
1364 : /* Normalize a non-negative integer */
1365 : GEN
1366 92841895 : int_normalize(GEN x, long known_zero_words)
1367 : {
1368 92841895 : long i = lgefint(x) - 1 - known_zero_words;
1369 1102013432 : for ( ; i > 1; i--)
1370 1069305804 : if (x[i]) { setlgefint(x, i+1); return x; }
1371 32707628 : x[1] = evalsigne(0) | evallgefint(2); return x;
1372 : }
1373 :
1374 : /********************************************************************
1375 : ** **
1376 : ** Base Conversion **
1377 : ** **
1378 : ********************************************************************/
1379 :
1380 : ulong *
1381 859292 : convi(GEN x, long *l)
1382 : {
1383 859292 : long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
1384 859292 : GEN str = cgetg(n+1, t_VECSMALL);
1385 859292 : unsigned char *res = (unsigned char*) GSTR(str);
1386 859292 : long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
1387 : long lz;
1388 : ulong *z;
1389 : long i, j;
1390 : unsigned char *t;
1391 859292 : while (!*res) {res++; llz--;} /*Strip leading zeros*/
1392 859292 : lz = (8+llz)/9;
1393 859292 : z = (ulong*)new_chunk(1+lz);
1394 859292 : t=res+llz+9;
1395 1937304 : for(i=0;i<llz-8;i+=9)
1396 : {
1397 : ulong s;
1398 1078012 : t-=18;
1399 1078012 : s=*t++;
1400 9702108 : for (j=1; j<9;j++)
1401 8624096 : s=10*s+*t++;
1402 1078012 : *z++=s;
1403 : }
1404 859292 : if (i<llz)
1405 : {
1406 848914 : unsigned char *t = res;
1407 848914 : ulong s=*t++;
1408 4357800 : for (j=i+1; j<llz;j++)
1409 3508886 : s=10*s+*t++;
1410 848914 : *z++=s;
1411 : }
1412 859292 : *l = lz;
1413 859292 : return z;
1414 : }
|
