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 1742873 : static void pari_gmp_free(void *ptr, size_t old_size){
52 1742873 : (void)old_size; pari_free(ptr);
53 1742873 : }
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 1100 : pari_kernel_init(void)
61 : {
62 1100 : mp_get_memory_functions (&old_gmp_malloc, &old_gmp_realloc, &old_gmp_free);
63 1100 : mp_set_memory_functions((void *(*)(size_t)) pari_malloc, pari_gmp_realloc, pari_gmp_free);
64 1100 : }
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 1092 : 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 1092 : mp_get_memory_functions (&new_gmp_malloc, &new_gmp_realloc, &new_gmp_free);
83 1092 : if (new_gmp_malloc==pari_malloc) new_gmp_malloc = old_gmp_malloc;
84 1092 : if (new_gmp_realloc==pari_gmp_realloc) new_gmp_realloc = old_gmp_realloc;
85 1092 : if (new_gmp_free==pari_gmp_free) new_gmp_free = old_gmp_free;
86 1092 : mp_set_memory_functions(new_gmp_malloc, new_gmp_realloc, new_gmp_free);
87 1092 : }
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 6874085 : xmpn_copy(mp_limb_t *x, mp_limb_t *y, long n)
97 : {
98 56485090 : while (--n >= 0) x[n]=y[n];
99 6874085 : }
100 :
101 : INLINE void
102 582885403 : xmpn_mirror(mp_limb_t *x, long n)
103 : {
104 : long i;
105 3868893082 : for(i=0;i<(n>>1);i++)
106 : {
107 3286007679 : ulong m=x[i];
108 3286007679 : x[i]=x[n-1-i];
109 3286007679 : x[n-1-i]=m;
110 : }
111 582885403 : }
112 :
113 : INLINE void
114 712408408 : xmpn_mirrorcopy(mp_limb_t *z, mp_limb_t *x, long n)
115 : {
116 : long i;
117 9807461335 : for(i=0;i<n;i++)
118 9095052927 : z[i]=x[n-1-i];
119 712408408 : }
120 :
121 : INLINE void
122 235850290 : xmpn_zero(mp_ptr x, mp_size_t n)
123 : {
124 2055659065 : while (--n >= 0) x[n]=0;
125 235850290 : }
126 :
127 : INLINE GEN
128 41121302 : icopy_ef(GEN x, long l)
129 : {
130 41121302 : long lx = lgefint(x);
131 41121302 : const GEN y = cgeti(l);
132 :
133 293748628 : while (--lx > 0) y[lx]=x[lx];
134 41119685 : 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 ordre to facilitate splitting: no
140 : * need to reintroduce codewords ]
141 : * Use speci(a,na) to visualize the corresponding GEN.
142 : */
143 :
144 : /***********************************************************************/
145 : /** **/
146 : /** ADDITION / SUBTRACTION **/
147 : /** **/
148 : /***********************************************************************/
149 :
150 : GEN
151 2997547 : setloop(GEN a)
152 : {
153 2997547 : pari_sp av = avma - 2 * sizeof(long);
154 2997547 : (void)cgetg(lgefint(a) + 3, t_VECSMALL);
155 2997547 : return icopy_avma(a, av); /* two cells of extra space after a */
156 : }
157 :
158 : /* we had a = setloop(?), then some incloops. Reset a to b */
159 : GEN
160 174328 : resetloop(GEN a, GEN b) {
161 174328 : a[0] = evaltyp(t_INT) | evallg(lgefint(b));
162 174328 : affii(b, a); return a;
163 : }
164 :
165 : /* assume a > 0, initialized by setloop. Do a++ */
166 : static GEN
167 99436923 : incpos(GEN a)
168 : {
169 99436923 : long i, l = lgefint(a);
170 99436928 : for (i=2; i<l; i++)
171 99443104 : if (++uel(a,i)) return a;
172 3 : a[l] = 1; l++;
173 3 : a[0]=evaltyp(t_INT) | _evallg(l);
174 3 : a[1]=evalsigne(1) | evallgefint(l);
175 3 : return a;
176 : }
177 :
178 : /* assume a < 0, initialized by setloop. Do a++ */
179 : static GEN
180 66684 : incneg(GEN a)
181 : {
182 66684 : long i, l = lgefint(a)-1;
183 66684 : if (uel(a,2)--)
184 : {
185 66680 : if (!a[l]) /* implies l = 2 */
186 : {
187 1980 : a[0] = evaltyp(t_INT) | _evallg(2);
188 1980 : a[1] = evalsigne(0) | evallgefint(2);
189 : }
190 66680 : return a;
191 : }
192 5 : for (i=3; i<=l; i++)
193 5 : if (uel(a,i)--) break;
194 4 : if (!a[l])
195 : {
196 4 : a[0] = evaltyp(t_INT) | _evallg(l);
197 4 : a[1] = evalsigne(-1) | evallgefint(l);
198 : }
199 4 : return a;
200 : }
201 :
202 : /* assume a initialized by setloop. Do a++ */
203 : GEN
204 99851458 : incloop(GEN a)
205 : {
206 99851458 : switch(signe(a))
207 : {
208 336632 : case 0:
209 336632 : a[0]=evaltyp(t_INT) | _evallg(3);
210 336632 : a[1]=evalsigne(1) | evallgefint(3);
211 336632 : a[2]=1; return a;
212 66684 : case -1: return incneg(a);
213 99448142 : default: return incpos(a);
214 : }
215 : }
216 :
217 : INLINE GEN
218 2773425053 : adduispec(ulong s, GEN x, long nx)
219 : {
220 : GEN zd;
221 : long lz;
222 :
223 2773425053 : if (nx == 1) return adduu(uel(x,0), s);
224 899896860 : lz = nx+3; zd = cgeti(lz);
225 902970379 : if (mpn_add_1(LIMBS(zd),(mp_limb_t *)x,nx,s))
226 5008503 : zd[lz-1]=1;
227 : else
228 897958913 : lz--;
229 902967416 : zd[1] = evalsigne(1) | evallgefint(lz);
230 902967416 : return zd;
231 : }
232 :
233 : GEN
234 764527742 : adduispec_offset(ulong s, GEN x, long offset, long nx)
235 : {
236 764527742 : GEN xd=x+2+offset;
237 998629753 : while (nx && *(xd+nx-1)==0) nx--;
238 764527742 : if (!nx) return utoi(s);
239 719888127 : return adduispec(s,xd,nx);
240 : }
241 :
242 : INLINE GEN
243 3312559869 : addiispec(GEN x, GEN y, long nx, long ny)
244 : {
245 : GEN zd;
246 : long lz;
247 :
248 3312559869 : if (nx < ny) swapspec(x,y, nx,ny);
249 3312559869 : if (ny == 1) return adduispec(*y,x,nx);
250 1347242565 : lz = nx+3; zd = cgeti(lz);
251 :
252 1347516079 : if (mpn_add(LIMBS(zd),(mp_limb_t *)x,nx,(mp_limb_t *)y,ny))
253 30415717 : zd[lz-1]=1;
254 : else
255 1317408637 : lz--;
256 :
257 1347824354 : zd[1] = evalsigne(1) | evallgefint(lz);
258 1347824354 : return zd;
259 : }
260 :
261 : /* assume x >= y */
262 : INLINE GEN
263 1758992162 : subiuspec(GEN x, ulong s, long nx)
264 : {
265 : GEN zd;
266 : long lz;
267 :
268 1758992162 : if (nx == 1) return utoi(x[0] - s);
269 :
270 250848446 : lz = nx + 2; zd = cgeti(lz);
271 251342742 : mpn_sub_1 (LIMBS(zd), (mp_limb_t *)x, nx, s);
272 251342358 : if (! zd[lz - 1]) { --lz; }
273 :
274 251342358 : zd[1] = evalsigne(1) | evallgefint(lz);
275 251342358 : return zd;
276 : }
277 :
278 : /* assume x > y */
279 : INLINE GEN
280 2994942874 : subiispec(GEN x, GEN y, long nx, long ny)
281 : {
282 : GEN zd;
283 : long lz;
284 2994942874 : if (ny==1) return subiuspec(x,*y,nx);
285 1266044173 : lz = nx+2; zd = cgeti(lz);
286 :
287 1264351302 : mpn_sub (LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
288 1738911992 : while (lz >= 3 && zd[lz - 1] == 0) { lz--; }
289 :
290 1264511427 : zd[1] = evalsigne(1) | evallgefint(lz);
291 1264511427 : return zd;
292 : }
293 :
294 : static void
295 521467876 : roundr_up_ip(GEN x, long l)
296 : {
297 521467876 : long i = l;
298 : for(;;)
299 : {
300 523294797 : if (++((ulong*)x)[--i]) break;
301 2179849 : if (i == 2) { x[2] = HIGHBIT; shiftr_inplace(x, 1); break; }
302 : }
303 521552894 : }
304 :
305 : void
306 400757112 : affir(GEN x, GEN y)
307 : {
308 400757112 : const long s = signe(x), ly = lg(y);
309 : long lx, sh, i;
310 :
311 400757112 : if (!s)
312 : {
313 41735006 : y[1] = evalexpo(-bit_accuracy(ly));
314 41732576 : return;
315 : }
316 359022106 : lx = lgefint(x); sh = bfffo(*int_MSW(x));
317 359022106 : y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
318 359165229 : if (sh) {
319 351649840 : if (lx <= ly)
320 : {
321 952608037 : for (i=lx; i<ly; i++) y[i]=0;
322 247247712 : mpn_lshift(LIMBS(y),LIMBS(x),lx-2,sh);
323 247323098 : xmpn_mirror(LIMBS(y),lx-2);
324 247638798 : return;
325 : }
326 104402128 : mpn_lshift(LIMBS(y),LIMBS(x)+lx-ly,ly-2,sh);
327 104401200 : uel(y,2) |= uel(x,lx-ly+1) >> (BITS_IN_LONG-sh);
328 104401200 : xmpn_mirror(LIMBS(y),ly-2);
329 : /* lx > ly: round properly */
330 104402638 : if ((uel(x,lx-ly+1)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
331 : }
332 : else {
333 7515389 : GEN xd=int_MSW(x);
334 7515389 : if (lx <= ly)
335 : {
336 9403594 : for (i=2; i<lx; i++,xd=int_precW(xd)) y[i]=*xd;
337 5118038 : for ( ; i<ly; i++) y[i]=0;
338 2353664 : return;
339 : }
340 14680137 : for (i=2; i<ly; i++,xd=int_precW(xd)) y[i]=*xd;
341 : /* lx > ly: round properly */
342 5161725 : if (uel(x,lx-ly+1) & HIGHBIT) roundr_up_ip(y, ly);
343 : }
344 : }
345 :
346 : INLINE GEN
347 705697837 : shiftispec(GEN x, long nx, long n)
348 : {
349 : long ny,m;
350 : GEN yd, y;
351 :
352 705697837 : if (!n) return icopyspec(x, nx);
353 :
354 676708801 : if (n > 0)
355 : {
356 396150211 : long d = dvmdsBIL(n, &m);
357 : long i;
358 :
359 396057224 : ny = nx + d + (m!=0);
360 396057224 : y = cgeti(ny + 2); yd = y + 2;
361 549622892 : for (i=0; i<d; i++) yd[i] = 0;
362 :
363 395188269 : if (!m) xmpn_copy((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx);
364 : else
365 : {
366 393651673 : ulong carryd = mpn_lshift((mp_limb_t *) (yd + d), (mp_limb_t *) x, nx, m);
367 393741095 : if (carryd) yd[ny - 1] = carryd;
368 370576172 : else ny--;
369 : }
370 : }
371 : else
372 : {
373 280558590 : long d = dvmdsBIL(-n, &m);
374 :
375 282452017 : ny = nx - d;
376 282452017 : if (ny < 1) return gen_0;
377 279557846 : y = cgeti(ny + 2); yd = y + 2;
378 :
379 279292362 : if (!m) xmpn_copy((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d);
380 : else
381 : {
382 274637913 : mpn_rshift((mp_limb_t *) yd, (mp_limb_t *) (x + d), nx - d, m);
383 274651537 : if (yd[ny - 1] == 0)
384 : {
385 60034727 : if (ny == 1) return gc_const((pari_sp)(yd + 1), gen_0);
386 50536006 : ny--;
387 : }
388 : }
389 : }
390 663259698 : y[1] = evalsigne(1)|evallgefint(ny + 2);
391 663259698 : return y;
392 : }
393 :
394 : GEN
395 137234173 : mantissa2nr(GEN x, long n)
396 : {
397 137234173 : long ly, i, m, s = signe(x), lx = lg(x);
398 : GEN y;
399 137234173 : if (!s) return gen_0;
400 137232822 : if (!n)
401 : {
402 30221201 : y = cgeti(lx);
403 30218335 : y[1] = evalsigne(s) | evallgefint(lx);
404 30218335 : xmpn_mirrorcopy(LIMBS(y),RLIMBS(x),lx-2);
405 30218057 : return y;
406 : }
407 107011621 : if (n > 0)
408 : {
409 218231 : GEN z = (GEN)avma;
410 218231 : long d = dvmdsBIL(n, &m);
411 :
412 218231 : ly = lx+d; y = new_chunk(ly);
413 556814 : for ( ; d; d--) *--z = 0;
414 222664 : if (!m) for (i=2; i<lx; i++) y[i]=x[i];
415 : else
416 : {
417 216787 : const ulong sh = BITS_IN_LONG - m;
418 216787 : shift_left(y,x, 2,lx-1, 0,m);
419 216787 : i = uel(x,2) >> sh;
420 : /* Extend y on the left? */
421 216787 : if (i) { ly++; y = new_chunk(1); y[2] = i; }
422 : }
423 : }
424 : else
425 : {
426 106793390 : ly = lx - dvmdsBIL(-n, &m);
427 106789858 : if (ly<3) return gen_0;
428 106789858 : y = new_chunk(ly);
429 106762255 : if (m) {
430 106503799 : shift_right(y,x, 2,ly, 0,m);
431 106551136 : if (y[2] == 0)
432 : {
433 0 : if (ly==3) return gc_const((pari_sp)(y+3), gen_0);
434 0 : ly--; set_avma((pari_sp)(++y));
435 : }
436 : } else {
437 701385 : for (i=2; i<ly; i++) y[i]=x[i];
438 : }
439 : }
440 107027823 : xmpn_mirror(LIMBS(y),ly-2);
441 107072336 : y[1] = evalsigne(s)|evallgefint(ly);
442 107072336 : y[0] = evaltyp(t_INT)|evallg(ly); return y;
443 : }
444 :
445 : GEN
446 3558885 : truncr(GEN x)
447 : {
448 : long s, e, d, m, i;
449 : GEN y;
450 3558885 : if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
451 1539986 : d = nbits2lg(e+1); m = remsBIL(e);
452 1539968 : if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
453 :
454 1539964 : y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
455 1539905 : if (++m == BITS_IN_LONG)
456 95943 : for (i=2; i<d; i++) y[d-i+1]=x[i];
457 : else
458 : {
459 1492192 : GEN z=cgeti(d);
460 3050973 : for (i=2; i<d; i++) z[d-i+1]=x[i];
461 1492088 : mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
462 1492089 : set_avma((pari_sp)y);
463 : }
464 1539778 : return y;
465 : }
466 :
467 : /* integral part */
468 : GEN
469 6984068 : floorr(GEN x)
470 : {
471 : long e, d, m, i, lx;
472 : GEN y;
473 6984068 : if (signe(x) >= 0) return truncr(x);
474 4200493 : if ((e=expo(x)) < 0) return gen_m1;
475 3520468 : d = nbits2lg(e+1); m = remsBIL(e);
476 3518069 : lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
477 3518069 : y = cgeti(d+1);
478 3509396 : if (++m == BITS_IN_LONG)
479 : {
480 3053 : for (i=2; i<d; i++) y[d-i+1]=x[i];
481 1450 : i=d; while (i<lx && !x[i]) i++;
482 906 : if (i==lx) goto END;
483 : }
484 : else
485 : {
486 3508490 : GEN z=cgeti(d);
487 7846694 : for (i=2; i<d; i++) z[d-i+1]=x[i];
488 3495230 : mpn_rshift(LIMBS(y),LIMBS(z),d-2,BITS_IN_LONG-m);
489 3495743 : if (uel(x,d-1)<<m == 0)
490 : {
491 510271 : i=d; while (i<lx && !x[i]) i++;
492 117390 : if (i==lx) goto END;
493 : }
494 : }
495 3421172 : if (mpn_add_1(LIMBS(y),LIMBS(y),d-2,1))
496 0 : y[d++]=1;
497 3421270 : END:
498 3496747 : y[1] = evalsigne(-1) | evallgefint(d);
499 3496747 : return y;
500 : }
501 :
502 : INLINE int
503 3893946860 : cmpiispec(GEN x, GEN y, long lx, long ly)
504 : {
505 3893946860 : if (lx < ly) return -1;
506 3574228090 : if (lx > ly) return 1;
507 3424903871 : return mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
508 : }
509 :
510 : INLINE int
511 265915486 : equaliispec(GEN x, GEN y, long lx, long ly)
512 : {
513 265915486 : if (lx != ly) return 0;
514 265774147 : return !mpn_cmp((mp_limb_t*)x,(mp_limb_t*)y, lx);
515 : }
516 :
517 : /***********************************************************************/
518 : /** **/
519 : /** MULTIPLICATION **/
520 : /** **/
521 : /***********************************************************************/
522 : /* assume ny > 0 */
523 : INLINE GEN
524 5482020407 : muluispec(ulong x, GEN y, long ny)
525 : {
526 5482020407 : if (ny == 1)
527 4388148518 : return muluu(x, *y);
528 : else
529 : {
530 1093871889 : long lz = ny+3;
531 1093871889 : GEN z = cgeti(lz);
532 1103501935 : ulong hi = mpn_mul_1 (LIMBS(z), (mp_limb_t *)y, ny, x);
533 1104833694 : if (hi) { z[lz - 1] = hi; } else lz--;
534 1104833694 : z[1] = evalsigne(1) | evallgefint(lz);
535 1104833694 : return z;
536 : }
537 : }
538 :
539 : /* a + b*|y| */
540 : GEN
541 0 : addumului(ulong a, ulong b, GEN y)
542 : {
543 : GEN z;
544 : long i, lz;
545 : ulong hi;
546 0 : if (!b || !signe(y)) return utoi(a);
547 0 : lz = lgefint(y)+1;
548 0 : z = cgeti(lz);
549 0 : z[2]=a;
550 0 : for(i=3;i<lz;i++) z[i]=0;
551 0 : hi=mpn_addmul_1(LIMBS(z), LIMBS(y), NLIMBS(y), b);
552 0 : if (hi) z[lz-1]=hi; else lz--;
553 0 : z[1] = evalsigne(1) | evallgefint(lz);
554 0 : return gc_const((pari_sp)z, z);
555 : }
556 :
557 : /***********************************************************************/
558 : /** **/
559 : /** DIVISION **/
560 : /** **/
561 : /***********************************************************************/
562 :
563 : ulong
564 1304508916 : umodiu(GEN y, ulong x)
565 : {
566 1304508916 : long sy=signe(y);
567 : ulong hi;
568 1304508916 : if (!x) pari_err_INV("umodiu",gen_0);
569 1305454700 : if (!sy) return 0;
570 896738473 : hi = mpn_mod_1(LIMBS(y),NLIMBS(y),x);
571 896738473 : if (!hi) return 0;
572 880522834 : return (sy > 0)? hi: x - hi;
573 : }
574 :
575 : /* return |y| \/ x */
576 : GEN
577 105518975 : absdiviu_rem(GEN y, ulong x, ulong *rem)
578 : {
579 : long ly;
580 : GEN z;
581 :
582 105518975 : if (!x) pari_err_INV("absdiviu_rem",gen_0);
583 105524946 : if (!signe(y)) { *rem = 0; return gen_0; }
584 :
585 82837993 : ly = lgefint(y);
586 82837993 : if (ly == 3 && (ulong)x > uel(y,2)) { *rem = uel(y,2); return gen_0; }
587 :
588 68512560 : z = cgeti(ly);
589 68512683 : *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
590 68512866 : if (z [ly - 1] == 0) ly--;
591 68512866 : z[1] = evallgefint(ly) | evalsigne(1);
592 68512866 : return z;
593 : }
594 :
595 : GEN
596 85283162 : divis_rem(GEN y, long x, long *rem)
597 : {
598 85283162 : long sy=signe(y),ly,s;
599 : GEN z;
600 :
601 85283162 : if (!x) pari_err_INV("divis_rem",gen_0);
602 85292311 : if (!sy) { *rem = 0; return gen_0; }
603 59930331 : if (x<0) { s = -sy; x = -x; } else s = sy;
604 :
605 59930331 : ly = lgefint(y);
606 59930331 : if (ly == 3 && (ulong)x > uel(y,2)) { *rem = itos(y); return gen_0; }
607 :
608 41472059 : z = cgeti(ly);
609 41467383 : *rem = mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
610 41467503 : if (sy<0) *rem = - *rem;
611 41467503 : if (z[ly - 1] == 0) ly--;
612 41467503 : z[1] = evallgefint(ly) | evalsigne(s);
613 41467503 : return z;
614 : }
615 :
616 : GEN
617 966689 : divis(GEN y, long x)
618 : {
619 966689 : long sy=signe(y),ly,s;
620 : GEN z;
621 :
622 966689 : if (!x) pari_err_INV("divis",gen_0);
623 966689 : if (!sy) return gen_0;
624 966641 : if (x<0) { s = -sy; x = -x; } else s=sy;
625 :
626 966641 : ly = lgefint(y);
627 966641 : if (ly == 3 && (ulong)x > uel(y,2)) return gen_0;
628 :
629 966329 : z = cgeti(ly);
630 966325 : (void)mpn_divrem_1(LIMBS(z), 0, LIMBS(y), NLIMBS(y), x);
631 966327 : if (z[ly - 1] == 0) ly--;
632 966327 : z[1] = evallgefint(ly) | evalsigne(s);
633 966327 : return z;
634 : }
635 :
636 : /* We keep llx bits of x and lly bits of y*/
637 : static GEN
638 75742126 : divrr_with_gmp(GEN x, GEN y)
639 : {
640 75742126 : long lx=RNLIMBS(x),ly=RNLIMBS(y);
641 75742126 : long lw=minss(lx,ly);
642 75742002 : long lly=minss(lw+1,ly);
643 75743020 : GEN w = cgetg(lw+2, t_REAL);
644 75724027 : long lu=lw+lly;
645 75724027 : long llx=minss(lu,lx);
646 75721798 : mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
647 75704613 : mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
648 : mp_limb_t *q,*r;
649 75675679 : long e=expo(x)-expo(y);
650 75675679 : long sx=signe(x),sy=signe(y);
651 75675679 : xmpn_mirrorcopy(z,RLIMBS(y),lly);
652 75695017 : xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
653 75736557 : xmpn_zero(u,lu-llx);
654 75761684 : q = (mp_limb_t *)new_chunk(lw+1);
655 75747721 : r = (mp_limb_t *)new_chunk(lly);
656 :
657 75726630 : mpn_tdiv_qr(q,r,0,u,lu,z,lly);
658 :
659 : /*Round up: This is not exactly correct we should test 2*r>z*/
660 75769381 : if (uel(r,lly-1) > (uel(z,lly-1)>>1))
661 37534021 : mpn_add_1(q,q,lw+1,1);
662 :
663 75769390 : xmpn_mirrorcopy(RLIMBS(w),q,lw);
664 :
665 75767276 : if (q[lw] == 0) e--;
666 41926416 : else if (q[lw] == 1) { shift_right(w,w, 2,lw+2, 1,1); }
667 0 : else { w[2] = HIGHBIT; e++; }
668 75766127 : if (sy < 0) sx = -sx;
669 75766127 : w[1] = evalsigne(sx) | evalexpo(e);
670 75760028 : return gc_const((pari_sp)w, w);
671 : }
672 :
673 : /* We keep llx bits of x and lly bits of y*/
674 : static GEN
675 35175503 : divri_with_gmp(GEN x, GEN y)
676 : {
677 35175503 : long llx=RNLIMBS(x),ly=NLIMBS(y);
678 35175503 : long lly=minss(llx+1,ly);
679 35175621 : GEN w = cgetg(llx+2, t_REAL);
680 35171250 : long lu=llx+lly, ld=ly-lly;
681 35171250 : mp_limb_t *u=(mp_limb_t *)new_chunk(lu);
682 35168360 : mp_limb_t *z=(mp_limb_t *)new_chunk(lly);
683 : mp_limb_t *q,*r;
684 35163063 : long sh=bfffo(y[ly+1]);
685 35163063 : long e=expo(x)-expi(y);
686 35164169 : long sx=signe(x),sy=signe(y);
687 35164169 : if (sh) mpn_lshift(z,LIMBS(y)+ld,lly,sh);
688 683057 : else xmpn_copy(z,LIMBS(y)+ld,lly);
689 35166285 : xmpn_mirrorcopy(u+lu-llx,RLIMBS(x),llx);
690 35172406 : xmpn_zero(u,lu-llx);
691 35177950 : q = (mp_limb_t *)new_chunk(llx+1);
692 35176002 : r = (mp_limb_t *)new_chunk(lly);
693 :
694 35171743 : mpn_tdiv_qr(q,r,0,u,lu,z,lly);
695 :
696 : /*Round up: This is not exactly correct we should test 2*r>z*/
697 35180518 : if (uel(r,lly-1) > (uel(z,lly-1)>>1))
698 16012908 : mpn_add_1(q,q,llx+1,1);
699 :
700 35180528 : xmpn_mirrorcopy(RLIMBS(w),q,llx);
701 :
702 35180147 : if (q[llx] == 0) e--;
703 10588405 : else if (q[llx] == 1) { shift_right(w,w, 2,llx+2, 1,1); }
704 0 : else { w[2] = HIGHBIT; e++; }
705 35180033 : if (sy < 0) sx = -sx;
706 35180033 : w[1] = evalsigne(sx) | evalexpo(e);
707 35178862 : return gc_const((pari_sp)w, w);
708 : }
709 :
710 : GEN
711 150194486 : divri(GEN x, GEN y)
712 : {
713 150194486 : long s = signe(y);
714 :
715 150194486 : if (!s) pari_err_INV("divri",gen_0);
716 150195141 : if (!signe(x)) return real_0_bit(expo(x) - expi(y));
717 149963319 : if (!is_bigint(y)) {
718 114792526 : GEN z = divru(x, y[2]);
719 114792496 : if (s < 0) togglesign(z);
720 114792503 : return z;
721 : }
722 35169798 : return divri_with_gmp(x,y);
723 : }
724 :
725 : GEN
726 141810397 : divrr(GEN x, GEN y)
727 : {
728 141810397 : long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
729 : ulong y0,y1;
730 : GEN r, r1;
731 :
732 141810397 : if (!sy) pari_err_INV("divrr",y);
733 141830705 : e = expo(x) - expo(y);
734 141830705 : if (!sx) return real_0_bit(e);
735 141326850 : if (sy<0) sx = -sx;
736 :
737 141326850 : lx=lg(x); ly=lg(y);
738 141326850 : if (ly==3)
739 : {
740 26296506 : ulong k = x[2], l = (lx>3)? x[3]: 0;
741 : LOCAL_HIREMAINDER;
742 26296506 : if (k < uel(y,2)) e--;
743 : else
744 : {
745 8275022 : l >>= 1; if (k&1) l |= HIGHBIT;
746 8275022 : k >>= 1;
747 : }
748 26296506 : hiremainder = k; k = divll(l,y[2]);
749 26296506 : if (hiremainder > (uel(y,2) >> 1) && !++k) { k = HIGHBIT; e++; }
750 26296506 : r = cgetg(3, t_REAL);
751 26293873 : r[1] = evalsigne(sx) | evalexpo(e);
752 26291764 : r[2] = k; return r;
753 : }
754 :
755 115030344 : if (ly >= prec2lg(DIVRR_GMP_LIMIT))
756 75742093 : return divrr_with_gmp(x,y);
757 :
758 39339020 : lr = minss(lx,ly); r = new_chunk(lr);
759 39347363 : r1 = r-1;
760 133635944 : r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
761 39347363 : r1[lr] = (lx>ly)? x[lr]: 0;
762 39347363 : y0 = y[2]; y1 = y[3];
763 172960671 : for (i=0; i<lr-1; i++)
764 : { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
765 : ulong k, qp;
766 : LOCAL_HIREMAINDER;
767 : LOCAL_OVERFLOW;
768 :
769 133613308 : if (uel(r1,1) == y0) { qp = ULONG_MAX; k = addll(y0,r1[2]); }
770 : else
771 : {
772 133367055 : if (uel(r1,1) > y0) /* can't happen if i=0 */
773 : {
774 0 : GEN y1 = y+1;
775 0 : j = lr-i; r1[j] = subll(r1[j],y1[j]);
776 0 : for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
777 0 : j=i; do uel(r,--j)++; while (j && !r[j]);
778 : }
779 133367055 : hiremainder = r1[1]; overflow = 0;
780 133367055 : qp = divll(r1[2],y0); k = hiremainder;
781 : }
782 133613308 : j = lr-i+1;
783 133613308 : if (!overflow)
784 : {
785 : long k3, k4;
786 133480749 : k3 = mulll(qp,y1);
787 133480749 : if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
788 39403950 : k4 = subll(hiremainder,k);
789 : else
790 : {
791 94076799 : k3 = subll(k3, r1[3]);
792 94076799 : k4 = subllx(hiremainder,k);
793 : }
794 171234745 : while (!overflow && k4) { qp--; k3=subll(k3,y1); k4=subllx(k4,y0); }
795 : }
796 133613308 : if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
797 396394943 : for (j--; j>1; j--)
798 : {
799 262781635 : r1[j] = subll(r1[j], addmul(qp,y[j]));
800 262781635 : hiremainder += overflow;
801 : }
802 133613308 : if (uel(r1,1) != hiremainder)
803 : {
804 182609 : if (uel(r1,1) < hiremainder)
805 : {
806 182609 : qp--;
807 182609 : j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
808 520274 : for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
809 : }
810 : else
811 : {
812 0 : uel(r1,1) -= hiremainder;
813 0 : while (r1[1])
814 : {
815 0 : qp++; if (!qp) { j=i; do uel(r,--j)++; while (j && !r[j]); }
816 0 : j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
817 0 : for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
818 0 : uel(r1,1) -= overflow;
819 : }
820 : }
821 : }
822 133613308 : *++r1 = qp;
823 : }
824 : /* i = lr-1 */
825 : /* round correctly */
826 39347363 : if (uel(r1,1) > (y0>>1))
827 : {
828 19842311 : j=i; do uel(r,--j)++; while (j && !r[j]);
829 : }
830 133821779 : r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
831 39347363 : if (r[0] == 0) e--;
832 14186618 : else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
833 : else { /* possible only when rounding up to 0x2 0x0 ... */
834 18 : r[2] = (long)HIGHBIT; e++;
835 : }
836 39346255 : r[0] = evaltyp(t_REAL)|evallg(lr);
837 39404974 : r[1] = evalsigne(sx) | evalexpo(e);
838 39394848 : return r;
839 : }
840 :
841 : /* Integer division x / y: such that sign(r) = sign(x)
842 : * if z = ONLY_REM return remainder, otherwise return quotient
843 : * if z != NULL set *z to remainder
844 : * *z is the last object on stack (and thus can be disposed of with cgiv
845 : * instead of gerepile)
846 : * If *z is zero, we put gen_0 here and no copy.
847 : * space needed: lx + ly
848 : */
849 : GEN
850 2169100318 : dvmdii(GEN x, GEN y, GEN *z)
851 : {
852 2169100318 : long sx=signe(x),sy=signe(y);
853 : long lx, ly, lq;
854 : pari_sp av;
855 : GEN r,q;
856 :
857 2169100318 : if (!sy) pari_err_INV("dvmdii",y);
858 2170205351 : if (!sx)
859 : {
860 68116816 : if (!z || z == ONLY_REM) return gen_0;
861 40575281 : *z=gen_0; return gen_0;
862 : }
863 2102088535 : lx=lgefint(x);
864 2102088535 : ly=lgefint(y); lq=lx-ly;
865 2102088535 : if (lq <= 0)
866 : {
867 1243127321 : if (lq == 0)
868 : {
869 1055319071 : long s=mpn_cmp(LIMBS(x),LIMBS(y),NLIMBS(x));
870 1055319071 : if (s>0) goto DIVIDE;
871 360850165 : if (s==0)
872 : {
873 30551123 : if (z == ONLY_REM) return gen_0;
874 20009097 : if (z) *z = gen_0;
875 20009097 : if (sx < 0) sy = -sy;
876 20009097 : return stoi(sy);
877 : }
878 : }
879 518107292 : if (z == ONLY_REM) return icopy(x);
880 13945769 : if (z) *z = icopy(x);
881 13945769 : return gen_0;
882 : }
883 858961214 : DIVIDE: /* quotient is nonzero */
884 1553430120 : av=avma; if (sx<0) sy = -sy;
885 1553430120 : if (ly==3)
886 : {
887 567057451 : ulong lq = lx;
888 : ulong si;
889 567057451 : q = cgeti(lq);
890 566645499 : si = mpn_divrem_1(LIMBS(q), 0, LIMBS(x), NLIMBS(x), y[2]);
891 567169879 : if (q[lq - 1] == 0) lq--;
892 567169879 : if (z == ONLY_REM)
893 : {
894 331823141 : if (!si) return gc_const(av, gen_0);
895 286251981 : set_avma(av); r = cgeti(3);
896 285743666 : r[1] = evalsigne(sx) | evallgefint(3);
897 285743666 : r[2] = si; return r;
898 : }
899 235346738 : q[1] = evalsigne(sy) | evallgefint(lq);
900 235346738 : if (!z) return q;
901 231251698 : if (!si) { *z=gen_0; return q; }
902 59823905 : r=cgeti(3);
903 59848912 : r[1] = evalsigne(sx) | evallgefint(3);
904 59848912 : r[2] = si; *z=r; return q;
905 : }
906 986372669 : if (z == ONLY_REM)
907 : {
908 963202411 : ulong lr = lgefint(y);
909 963202411 : ulong lq = lgefint(x)-lgefint(y)+3;
910 963202411 : GEN r = cgeti(lr);
911 958080444 : GEN q = cgeti(lq);
912 951539088 : mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
913 965344209 : if (!r[lr - 1])
914 : {
915 1180381356 : while(lr>2 && !r[lr - 1]) lr--;
916 544908819 : if (lr == 2) return gc_const(av, gen_0); /* exact division */
917 : }
918 950468105 : r[1] = evalsigne(sx) | evallgefint(lr);
919 950468105 : return gc_const((pari_sp)r, r);
920 : }
921 : else
922 : {
923 23170258 : ulong lq = lgefint(x)-lgefint(y)+3;
924 23170258 : ulong lr = lgefint(y);
925 23170258 : GEN q = cgeti(lq);
926 28366610 : GEN r = cgeti(lr);
927 28352766 : mpn_tdiv_qr(LIMBS(q), LIMBS(r),0, LIMBS(x), NLIMBS(x), LIMBS(y), NLIMBS(y));
928 28379773 : if (q[lq - 1] == 0) lq--;
929 28379773 : q[1] = evalsigne(sy) | evallgefint(lq);
930 28379773 : if (!z) return gc_const((pari_sp)q, q);
931 24682414 : if (!r[lr - 1])
932 : {
933 37789018 : while(lr>2 && !r[lr - 1]) lr--;
934 6229606 : if (lr == 2) { *z = gen_0; return gc_const((pari_sp)q, q); } /* exact */
935 : }
936 20009276 : r[1] = evalsigne(sx) | evallgefint(lr);
937 20009276 : *z = gc_const((pari_sp)r, r); return q;
938 : }
939 : }
940 :
941 : /* Montgomery reduction.
942 : * N has k words, assume T >= 0 has less than 2k.
943 : * Return res := T / B^k mod N, where B = 2^BIL
944 : * such that 0 <= res < T/B^k + N and res has less than k words */
945 : GEN
946 36730312 : red_montgomery(GEN T, GEN N, ulong inv)
947 : {
948 : pari_sp av;
949 : GEN Te, Td, Ne, Nd, scratch;
950 36730312 : ulong i, j, m, t, d, k = NLIMBS(N);
951 : int carry;
952 : LOCAL_HIREMAINDER;
953 : LOCAL_OVERFLOW;
954 :
955 36730312 : if (k == 0) return gen_0;
956 36730312 : d = NLIMBS(T); /* <= 2*k */
957 36730312 : if (d == 0) return gen_0;
958 : #ifdef DEBUG
959 : if (d > 2*k) pari_err_BUG("red_montgomery");
960 : #endif
961 36730303 : if (k == 1)
962 : { /* as below, special cased for efficiency */
963 163341 : ulong n = uel(N,2);
964 163341 : if (d == 1) {
965 163194 : hiremainder = uel(T,2);
966 163194 : m = hiremainder * inv;
967 163194 : (void)addmul(m, n); /* t + m*n = 0 */
968 163194 : return utoi(hiremainder);
969 : } else { /* d = 2 */
970 147 : hiremainder = uel(T,2);
971 147 : m = hiremainder * inv;
972 147 : (void)addmul(m, n); /* t + m*n = 0 */
973 147 : t = addll(hiremainder, uel(T,3));
974 147 : if (overflow) t -= n; /* t > n doesn't fit in 1 word */
975 147 : return utoi(t);
976 : }
977 : }
978 : /* assume k >= 2 */
979 36566962 : av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
980 :
981 : /* copy T to scratch space (pad with zeroes to 2k words) */
982 36273310 : Td = scratch;
983 36273310 : Te = T + 2;
984 527809525 : for (i=0; i < d ; i++) *Td++ = *Te++;
985 64623297 : for ( ; i < (k<<1); i++) *Td++ = 0;
986 :
987 36273310 : Te = scratch - 1; /* 1 beyond end of current T mantissa (in scratch) */
988 36273310 : Ne = N + 1; /* 1 beyond end of N mantissa */
989 :
990 36273310 : carry = 0;
991 276216077 : for (i=0; i<k; i++) /* set T := T/B nod N, k times */
992 : {
993 239942767 : Td = Te; /* one beyond end of (new) T mantissa */
994 239942767 : Nd = Ne;
995 239942767 : hiremainder = *++Td;
996 239942767 : m = hiremainder * inv; /* solve T + m N = O(B) */
997 :
998 : /* set T := (T + mN) / B */
999 239942767 : Te = Td;
1000 239942767 : (void)addmul(m, *++Nd); /* = 0 */
1001 2126595254 : for (j=1; j<k; j++)
1002 : {
1003 1886652487 : t = addll(addmul(m, *++Nd), *++Td);
1004 1886652487 : *Td = t;
1005 1886652487 : hiremainder += overflow;
1006 : }
1007 239942767 : t = addll(hiremainder, *++Td); *Td = t + carry;
1008 239942767 : carry = (overflow || (carry && *Td == 0));
1009 : }
1010 36273310 : if (carry)
1011 : { /* Td > N overflows (k+1 words), set Td := Td - N */
1012 57167 : GEN NE = N + k+1;
1013 57167 : Td = Te;
1014 57167 : Nd = Ne;
1015 57167 : t = subll(*++Td, *++Nd); *Td = t;
1016 514602 : while (Nd < NE) { t = subllx(*++Td, *++Nd); *Td = t; }
1017 : }
1018 :
1019 : /* copy result */
1020 36273310 : Td = (GEN)av - 1; /* *Td = high word of final result */
1021 39789773 : while (*Td == 0 && Te < Td) Td--; /* strip leading 0s */
1022 36273310 : k = Td - Te; if (!k) return gc_const(av, gen_0);
1023 36273310 : Td = (GEN)av - k; /* will write mantissa there */
1024 36273310 : (void)memmove(Td, Te+1, k*sizeof(long));
1025 36273310 : Td -= 2;
1026 36273310 : Td[0] = evaltyp(t_INT) | evallg(k+2);
1027 36467734 : Td[1] = evalsigne(1) | evallgefint(k+2);
1028 : #ifdef DEBUG
1029 : {
1030 : long l = NLIMBS(N), s = BITS_IN_LONG*l;
1031 : GEN R = int2n(s);
1032 : GEN res = remii(mulii(T, Fp_inv(R, N)), N);
1033 : if (k > lgefint(N)
1034 : || !equalii(remii(Td,N),res)
1035 : || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
1036 : }
1037 : #endif
1038 36467734 : return gc_const((pari_sp)Td, Td);
1039 : }
1040 :
1041 : /* EXACT INTEGER DIVISION */
1042 :
1043 : /* use undocumented GMP interface */
1044 : static void
1045 115243836 : GEN2mpz(mpz_t X, GEN x)
1046 : {
1047 115243836 : long l = lgefint(x)-2;
1048 115243836 : X->_mp_alloc = l;
1049 115243836 : X->_mp_size = signe(x) > 0? l: -l;
1050 115243836 : X->_mp_d = LIMBS(x);
1051 115243836 : }
1052 : static void
1053 57623196 : mpz2GEN(GEN z, mpz_t Z)
1054 : {
1055 57623196 : long l = Z->_mp_size;
1056 57623196 : z[1] = evalsigne(l > 0? 1: -1) | evallgefint(labs(l)+2);
1057 57623196 : }
1058 :
1059 : #ifdef mpn_divexact_1
1060 : static GEN
1061 387001640 : diviuexact_i(GEN x, ulong y)
1062 : {
1063 387001640 : long l = lgefint(x);
1064 387001640 : GEN z = cgeti(l);
1065 386549300 : mpn_divexact_1(LIMBS(z), LIMBS(x), NLIMBS(x), y);
1066 386634567 : if (z[l-1] == 0) l--;
1067 386634567 : z[1] = evallgefint(l) | evalsigne(signe(x));
1068 386634567 : return z;
1069 : }
1070 : #elif 1 && !defined(_WIN64) /* mpz_divexact_ui is not LLP64 friendly */
1071 : /* assume y != 0 and the division is exact */
1072 : static GEN
1073 : diviuexact_i(GEN x, ulong y)
1074 : {
1075 : long l = lgefint(x);
1076 : GEN z = cgeti(l);
1077 : mpz_t X, Z;
1078 : GEN2mpz(X, x);
1079 : Z->_mp_alloc = l-2;
1080 : Z->_mp_size = l-2;
1081 : Z->_mp_d = LIMBS(z);
1082 : mpz_divexact_ui(Z, X, y);
1083 : mpz2GEN(z, Z); return z;
1084 : }
1085 : #else
1086 : /* assume y != 0 and the division is exact */
1087 : static GEN
1088 : diviuexact_i(GEN x, ulong y)
1089 : {
1090 : /*TODO: implement true exact division.*/
1091 : return divii(x,utoi(y));
1092 : }
1093 : #endif
1094 :
1095 : GEN
1096 30836211 : diviuexact(GEN x, ulong y)
1097 : {
1098 : GEN z;
1099 30836211 : if (!signe(x)) return gen_0;
1100 29707635 : z = diviuexact_i(x, y);
1101 29704098 : if (lgefint(z) == 2) pari_err_OP("exact division", x, utoi(y));
1102 29703821 : return z;
1103 : }
1104 :
1105 : /* Find z such that x=y*z, knowing that y | x (unchecked) */
1106 : GEN
1107 507243227 : diviiexact(GEN x, GEN y)
1108 : {
1109 : GEN z;
1110 507243227 : if (!signe(y)) pari_err_INV("diviiexact",y);
1111 507414033 : if (!signe(x)) return gen_0;
1112 415128271 : if (lgefint(y) == 3)
1113 : {
1114 357511681 : z = diviuexact_i(x, y[2]);
1115 357137250 : if (signe(y) < 0) togglesign(z);
1116 : }
1117 : else
1118 : {
1119 57616590 : long l = lgefint(x);
1120 : mpz_t X, Y, Z;
1121 57616590 : z = cgeti(l);
1122 57622055 : GEN2mpz(X, x);
1123 57621997 : GEN2mpz(Y, y);
1124 57621881 : Z->_mp_alloc = l-2;
1125 57621881 : Z->_mp_size = l-2;
1126 57621881 : Z->_mp_d = LIMBS(z);
1127 57621881 : mpz_divexact(Z, X, Y);
1128 57623215 : mpz2GEN(z, Z);
1129 : }
1130 414760436 : if (lgefint(z) == 2) pari_err_OP("exact division", x, y);
1131 414739033 : return z;
1132 : }
1133 :
1134 : /* assume yz != and yz | x */
1135 : GEN
1136 199424 : diviuuexact(GEN x, ulong y, ulong z)
1137 : {
1138 : long tmp[4];
1139 : ulong t;
1140 : LOCAL_HIREMAINDER;
1141 199424 : t = mulll(y, z);
1142 199424 : if (!hiremainder) return diviuexact(x, t);
1143 1 : tmp[0] = evaltyp(t_INT)|_evallg(4);
1144 1 : tmp[1] = evalsigne(1)|evallgefint(4);
1145 1 : tmp[2] = t;
1146 1 : tmp[3] = hiremainder;
1147 1 : return diviiexact(x, tmp);
1148 : }
1149 :
1150 : /********************************************************************/
1151 : /** **/
1152 : /** INTEGER MULTIPLICATION **/
1153 : /** **/
1154 : /********************************************************************/
1155 :
1156 : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
1157 : GEN
1158 5823376312 : muliispec(GEN x, GEN y, long nx, long ny)
1159 : {
1160 : GEN zd;
1161 : long lz;
1162 : ulong hi;
1163 :
1164 5823376312 : if (nx < ny) swapspec(x,y, nx,ny);
1165 5823376312 : if (!ny) return gen_0;
1166 5823376312 : if (ny == 1) return muluispec((ulong)*y, x, nx);
1167 :
1168 1045381238 : lz = nx+ny+2;
1169 1045381238 : zd = cgeti(lz);
1170 1048658986 : hi = mpn_mul(LIMBS(zd), (mp_limb_t *)x, nx, (mp_limb_t *)y, ny);
1171 1055761488 : if (!hi) lz--;
1172 : /*else zd[lz-1]=hi; GH tell me it is not necessary.*/
1173 :
1174 1055761488 : zd[1] = evalsigne(1) | evallgefint(lz);
1175 1055761488 : return zd;
1176 : }
1177 : GEN
1178 221641 : muluui(ulong x, ulong y, GEN z)
1179 : {
1180 221641 : long t, s = signe(z);
1181 : GEN r;
1182 : LOCAL_HIREMAINDER;
1183 :
1184 221641 : if (!x || !y || !signe(z)) return gen_0;
1185 221262 : t = mulll(x,y);
1186 221262 : if (!hiremainder)
1187 221289 : r = muluispec(t, z+2, lgefint(z)-2);
1188 : else
1189 : {
1190 : long tmp[2];
1191 0 : tmp[1] = hiremainder;
1192 0 : tmp[0] = t;
1193 0 : r = muliispec(z+2,tmp, lgefint(z)-2, 2);
1194 : }
1195 221261 : setsigne(r,s); return r;
1196 : }
1197 :
1198 : GEN
1199 1015268303 : sqrispec(GEN x, long nx)
1200 : {
1201 : GEN zd;
1202 : long lz;
1203 :
1204 1015268303 : if (!nx) return gen_0;
1205 480612127 : if (nx==1) return sqru(*x);
1206 :
1207 275324011 : lz = (nx<<1)+2;
1208 275324011 : zd = cgeti(lz);
1209 : #ifdef mpn_sqr
1210 272633009 : mpn_sqr(LIMBS(zd), (mp_limb_t *)x, nx);
1211 : #else
1212 : mpn_mul_n(LIMBS(zd), (mp_limb_t *)x, (mp_limb_t *)x, nx);
1213 : #endif
1214 276064617 : if (zd[lz-1]==0) lz--;
1215 :
1216 276064617 : zd[1] = evalsigne(1) | evallgefint(lz);
1217 276064617 : return zd;
1218 : }
1219 :
1220 : INLINE GEN
1221 41389554 : sqrispec_mirror(GEN x, long nx)
1222 : {
1223 41389554 : GEN cx=new_chunk(nx);
1224 : GEN z;
1225 41353846 : xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
1226 41465277 : z=sqrispec(cx, nx);
1227 41517602 : xmpn_mirror(LIMBS(z), NLIMBS(z));
1228 41519297 : return z;
1229 : }
1230 :
1231 : /* leaves garbage on the stack. */
1232 : INLINE GEN
1233 83762575 : muliispec_mirror(GEN x, GEN y, long nx, long ny)
1234 : {
1235 : GEN cx, cy, z;
1236 83762575 : long s = 0;
1237 114828410 : while (nx && x[nx-1]==0) { nx--; s++; }
1238 120978327 : while (ny && y[ny-1]==0) { ny--; s++; }
1239 83762575 : cx=new_chunk(nx); cy=new_chunk(ny);
1240 83337514 : xmpn_mirrorcopy((mp_limb_t *)cx,(mp_limb_t *)x,nx);
1241 84315983 : xmpn_mirrorcopy((mp_limb_t *)cy,(mp_limb_t *)y,ny);
1242 84855364 : z = nx>=ny ? muliispec(cx, cy, nx, ny): muliispec(cy, cx, ny, nx);
1243 84890143 : if (s)
1244 : {
1245 7695728 : long i, lz = lgefint(z) + s;
1246 7695728 : (void)new_chunk(s);
1247 7695723 : z -= s;
1248 75977308 : for (i=0; i<s; i++) z[2+i]=0;
1249 7695723 : z[1] = evalsigne(1) | evallgefint(lz);
1250 7695723 : z[0] = evaltyp(t_INT) | evallg(lz);
1251 : }
1252 84890140 : xmpn_mirror(LIMBS(z), NLIMBS(z));
1253 85558032 : return z;
1254 : }
1255 :
1256 : /* x % (2^n), assuming n >= 0 */
1257 : GEN
1258 36060131 : remi2n(GEN x, long n)
1259 : {
1260 : ulong hi;
1261 36060131 : long l, k, lx, ly, sx = signe(x);
1262 : GEN z, xd, zd;
1263 :
1264 36060131 : if (!sx || !n) return gen_0;
1265 :
1266 35741541 : k = dvmdsBIL(n, &l);
1267 35752312 : lx = lgefint(x);
1268 35752312 : if (lx < k+3) return icopy(x);
1269 :
1270 34874970 : xd = x + (2 + k);
1271 : /* x = |k|...|1|#|... : copy the last l bits of # and the first k words
1272 : * ^--- initial xd */
1273 34874970 : hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
1274 34874970 : if (!hi)
1275 : { /* strip leading zeroes from result */
1276 3150478 : xd--; while (k && !*xd) { k--; xd--; }
1277 2965490 : if (!k) return gen_0;
1278 2040343 : ly = k+2;
1279 : }
1280 : else
1281 31909480 : ly = k+3;
1282 :
1283 33949823 : zd = z = cgeti(ly);
1284 33915148 : *++zd = evalsigne(sx) | evallgefint(ly);
1285 500670461 : xd = x+1; for ( ;k; k--) *++zd = *++xd;
1286 33915148 : if (hi) *++zd = hi;
1287 33915148 : return z;
1288 : }
1289 :
1290 : /********************************************************************/
1291 : /** **/
1292 : /** INTEGER SQUARE ROOT **/
1293 : /** **/
1294 : /********************************************************************/
1295 :
1296 : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
1297 : * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
1298 : * remainder is 0. R = NULL is allowed. */
1299 : GEN
1300 5114919 : sqrtremi(GEN a, GEN *r)
1301 : {
1302 5114919 : long l, na = NLIMBS(a);
1303 : mp_size_t nr;
1304 : GEN S;
1305 5114919 : if (!na) {
1306 724 : if (r) *r = gen_0;
1307 724 : return gen_0;
1308 : }
1309 5114195 : l = (na + 5) >> 1; /* 2 + ceil(na/2) */
1310 5114195 : S = cgetipos(l);
1311 5114134 : if (r) {
1312 1309361 : GEN R = cgeti(2 + na);
1313 1309361 : nr = mpn_sqrtrem(LIMBS(S), LIMBS(R), LIMBS(a), na);
1314 1309361 : if (nr) R[1] = evalsigne(1) | evallgefint(nr+2);
1315 25720 : else { set_avma((pari_sp)S); R = gen_0; }
1316 1309361 : *r = R;
1317 : }
1318 : else
1319 3804773 : (void)mpn_sqrtrem(LIMBS(S), NULL, LIMBS(a), na);
1320 5114135 : return S;
1321 : }
1322 :
1323 : /* compute sqrt(|a|), assuming a != 0 */
1324 : GEN
1325 125159665 : sqrtr_abs(GEN a)
1326 : {
1327 : GEN res;
1328 : mp_limb_t *b, *c;
1329 125159665 : long l = RNLIMBS(a), e = expo(a), er = e>>1;
1330 : long n;
1331 125159665 : res = cgetg(2 + l, t_REAL);
1332 125076490 : res[1] = evalsigne(1) | evalexpo(er);
1333 125156593 : if (e&1)
1334 : {
1335 52580689 : b = (mp_limb_t *) new_chunk(l<<1);
1336 52568016 : xmpn_zero(b,l);
1337 52569197 : xmpn_mirrorcopy(b+l, RLIMBS(a), l);
1338 52587662 : c = (mp_limb_t *) new_chunk(l);
1339 52583062 : n = mpn_sqrtrem(c,b,b,l<<1); /* c <- sqrt; b <- rem */
1340 52612726 : if (n>l || (n==l && mpn_cmp(b,c,l) > 0)) mpn_add_1(c,c,l,1);
1341 : }
1342 : else
1343 : {
1344 : ulong u;
1345 72575904 : b = (mp_limb_t *) mantissa2nr(a,-1);
1346 72617778 : b[1] = uel(a,l+1)<<(BITS_IN_LONG-1);
1347 72617778 : b = (mp_limb_t *) new_chunk(l);
1348 72597759 : xmpn_zero(b,l+1); /* overwrites the former b[0] */
1349 72600714 : c = (mp_limb_t *) new_chunk(l + 1);
1350 72567661 : n = mpn_sqrtrem(c,b,b,(l<<1)+2); /* c <- sqrt; b <- rem */
1351 72653907 : u = (ulong)*c++;
1352 72653907 : if ( u&HIGHBIT || (u == ~HIGHBIT &&
1353 0 : (n>l || (n==l && mpn_cmp(b,c,l)>0))))
1354 35796951 : mpn_add_1(c,c,l,1);
1355 : }
1356 125275046 : xmpn_mirrorcopy(RLIMBS(res),c,l);
1357 125248042 : return gc_const((pari_sp)res, res);
1358 : }
1359 :
1360 : /* Normalize a nonnegative integer */
1361 : GEN
1362 306844965 : int_normalize(GEN x, long known_zero_words)
1363 : {
1364 306844965 : long i = lgefint(x) - 1 - known_zero_words;
1365 2772109882 : for ( ; i > 1; i--)
1366 2720025076 : if (x[i]) { setlgefint(x, i+1); return x; }
1367 52084806 : x[1] = evalsigne(0) | evallgefint(2); return x;
1368 : }
1369 :
1370 : /********************************************************************
1371 : ** **
1372 : ** Base Conversion **
1373 : ** **
1374 : ********************************************************************/
1375 :
1376 : ulong *
1377 437505 : convi(GEN x, long *l)
1378 : {
1379 437505 : long n = nchar2nlong(2 + (long)(NLIMBS(x) * (BITS_IN_LONG * LOG10_2)));
1380 437505 : GEN str = cgetg(n+1, t_VECSMALL);
1381 437505 : unsigned char *res = (unsigned char*) GSTR(str);
1382 437505 : long llz = mpn_get_str(res, 10, LIMBS(icopy(x)), NLIMBS(x));
1383 : long lz;
1384 : ulong *z;
1385 : long i, j;
1386 : unsigned char *t;
1387 437505 : while (!*res) {res++; llz--;} /*Strip leading zeros*/
1388 437505 : lz = (8+llz)/9;
1389 437505 : z = (ulong*)new_chunk(1+lz);
1390 437505 : t=res+llz+9;
1391 867256 : for(i=0;i<llz-8;i+=9)
1392 : {
1393 : ulong s;
1394 429751 : t-=18;
1395 429751 : s=*t++;
1396 3867759 : for (j=1; j<9;j++)
1397 3438008 : s=10*s+*t++;
1398 429751 : *z++=s;
1399 : }
1400 437505 : if (i<llz)
1401 : {
1402 433546 : unsigned char *t = res;
1403 433546 : ulong s=*t++;
1404 1225685 : for (j=i+1; j<llz;j++)
1405 792139 : s=10*s+*t++;
1406 433546 : *z++=s;
1407 : }
1408 437505 : *l = lz;
1409 437505 : return z;
1410 : }
|
