Line data Source code
1 : #line 2 "../src/kernel/none/mp.c"
2 : /* Copyright (C) 2000-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 : /** MULTIPRECISION KERNEL **/
19 : /** **/
20 : /***********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 : #include "../src/kernel/none/tune-gen.h"
24 :
25 : void
26 779 : pari_kernel_init(void) { }
27 : void
28 777 : pari_kernel_close(void) { }
29 : const char *
30 2 : pari_kernel_version(void) { return ""; }
31 :
32 : /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
33 : * GENs but pairs (long *a, long na) representing a list of digits (in basis
34 : * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
35 : * need to reintroduce codewords ] */
36 :
37 : #define LIMBS(x) ((x)+2)
38 : #define NLIMBS(x) (lgefint(x)-2)
39 :
40 : /* Normalize a nonnegative integer */
41 : GEN
42 850535667 : int_normalize(GEN x, long known_zero_words)
43 : {
44 850535667 : long i, lx = lgefint(x);
45 : GEN x0;
46 850535667 : if (lx == 2) { x[1] = evalsigne(0) | evallgefint(2); return x; }
47 850535667 : if (!known_zero_words && x[2]) return x;
48 3527972817 : for (i = 2+known_zero_words; i < lx; i++)
49 3460516830 : if (x[i]) break;
50 331278855 : x0 = x; i -= 2; x += i;
51 331278855 : if (x0 == (GEN)avma) set_avma((pari_sp)x);
52 198663021 : else stackdummy((pari_sp)(x0+i), (pari_sp)x0);
53 331278855 : lx -= i;
54 331278855 : x[0] = evaltyp(t_INT) | evallg(lx);
55 331278855 : if (lx == 2) x[1] = evalsigne(0) | evallgefint(lx);
56 263822868 : else x[1] = evalsigne(1) | evallgefint(lx);
57 331278855 : return x;
58 : }
59 :
60 : /***********************************************************************/
61 : /** **/
62 : /** ADDITION / SUBTRACTION **/
63 : /** **/
64 : /***********************************************************************/
65 :
66 : GEN
67 2247084 : setloop(GEN a)
68 : {
69 2247084 : pari_sp av = avma;
70 2247084 : (void)cgetg(lgefint(a) + 3, t_VECSMALL);
71 2247084 : return icopy_avma(a, av); /* two cells of extra space before a */
72 : }
73 :
74 : /* we had a = setloop(?), then some incloops. Reset a to b */
75 : GEN
76 130656 : resetloop(GEN a, GEN b) {
77 130656 : long lb = lgefint(b);
78 130656 : a += lgefint(a) - lb;
79 130656 : a[0] = evaltyp(t_INT) | evallg(lb);
80 130656 : affii(b, a); return a;
81 : }
82 :
83 : /* assume a > 0, initialized by setloop. Do a++ */
84 : static GEN
85 31954041 : incpos(GEN a)
86 : {
87 31954041 : long i, l = lgefint(a);
88 31954044 : for (i=l-1; i>1; i--)
89 31954041 : if (++a[i]) return a;
90 3 : l++; a--; /* use extra cell */
91 3 : a[0]=evaltyp(t_INT) | _evallg(l);
92 3 : a[1]=evalsigne(1) | evallgefint(l);
93 3 : a[2]=1; return a;
94 : }
95 :
96 : /* assume a < 0, initialized by setloop. Do a++ */
97 : static GEN
98 50013 : incneg(GEN a)
99 : {
100 50013 : long i, l = lgefint(a)-1;
101 50013 : if (uel(a,l)--)
102 : {
103 50010 : if (l == 2 && !a[2])
104 : {
105 1485 : a++; /* save one cell */
106 1485 : a[0] = evaltyp(t_INT) | _evallg(2);
107 1485 : a[1] = evalsigne(0) | evallgefint(2);
108 : }
109 50010 : return a;
110 : }
111 3 : for (i = l-1;; i--) /* finishes since a[2] != 0 */
112 3 : if (uel(a,i)--) break;
113 3 : if (!a[2])
114 : {
115 3 : a++; /* save one cell */
116 3 : a[0] = evaltyp(t_INT) | _evallg(l);
117 3 : a[1] = evalsigne(-1) | evallgefint(l);
118 : }
119 3 : return a;
120 : }
121 :
122 : /* assume a initialized by setloop. Do a++ */
123 : GEN
124 32256525 : incloop(GEN a)
125 : {
126 32256525 : switch(signe(a))
127 : {
128 252471 : case 0: a--; /* use extra cell */
129 252471 : a[0]=evaltyp(t_INT) | _evallg(3);
130 252471 : a[1]=evalsigne(1) | evallgefint(3);
131 252471 : a[2]=1; return a;
132 50013 : case -1: return incneg(a);
133 31954041 : default: return incpos(a);
134 : }
135 : }
136 :
137 : INLINE GEN
138 2389856328 : adduispec(ulong s, GEN x, long nx)
139 : {
140 2389856328 : GEN xd, zd = (GEN)avma;
141 : long lz;
142 :
143 2389856328 : if (nx == 1) return adduu(s, uel(x,0));
144 877470162 : lz = nx+3; (void)new_chunk(lz);
145 877470162 : xd = x + nx;
146 877470162 : *--zd = (ulong)*--xd + s;
147 877470162 : if ((ulong)*zd < s)
148 : for(;;)
149 : {
150 263390214 : if (xd == x) { *--zd = 1; break; } /* enlarge z */
151 259679952 : *--zd = ((ulong)*--xd) + 1;
152 259679952 : if (*zd) { lz--; break; }
153 : }
154 620763606 : else lz--;
155 2175372852 : while (xd > x) *--zd = *--xd;
156 877470162 : *--zd = evalsigne(1) | evallgefint(lz);
157 877470162 : *--zd = evaltyp(t_INT) | evallg(lz);
158 877470162 : return gc_const((pari_sp)zd, zd);
159 : }
160 :
161 : GEN
162 491446833 : adduispec_offset(ulong s, GEN x, long offset, long nx)
163 : {
164 491446833 : GEN xd = x+lgefint(x)-nx-offset;
165 606252177 : while (nx && *xd==0) {xd++; nx--;}
166 491446833 : if (!nx) return utoi(s);
167 459445197 : return adduispec(s,xd,nx);
168 : }
169 :
170 : static GEN
171 4569102876 : addiispec(GEN x, GEN y, long nx, long ny)
172 : {
173 : GEN xd, yd, zd;
174 4569102876 : long lz, i = -2;
175 : LOCAL_OVERFLOW;
176 :
177 4569102876 : if (nx < ny) swapspec(x,y, nx,ny);
178 4569102876 : if (ny == 1) return adduispec(*y,x,nx);
179 2703838425 : zd = (GEN)avma;
180 2703838425 : lz = nx+3; (void)new_chunk(lz);
181 2703838425 : xd = x + nx;
182 2703838425 : yd = y + ny;
183 2703838425 : zd[-1] = addll(xd[-1], yd[-1]);
184 : #ifdef addllx8
185 2407359283 : for ( ; i-8 > -ny; i-=8)
186 1506079808 : addllx8(xd+i, yd+i, zd+i, overflow);
187 : #endif
188 38073128849 : for ( ; i >= -ny; i--) zd[i] = addllx(xd[i], yd[i]);
189 2703838425 : if (overflow)
190 : for(;;)
191 : {
192 549631152 : if (i < -nx) { zd[i] = 1; i--; break; } /* enlarge z */
193 366695808 : zd[i] = uel(xd,i) + 1;
194 366695808 : if (zd[i]) { i--; lz--; break; }
195 63440619 : i--;
196 : }
197 2217647892 : else lz--;
198 20516724981 : for (; i >= -nx; i--) zd[i] = xd[i];
199 2703838425 : zd += i+1;
200 2703838425 : *--zd = evalsigne(1) | evallgefint(lz);
201 2703838425 : *--zd = evaltyp(t_INT) | evallg(lz);
202 2703838425 : return gc_const((pari_sp)zd, zd);
203 : }
204 :
205 : /* assume x >= s */
206 : INLINE GEN
207 1581491034 : subiuspec(GEN x, ulong s, long nx)
208 : {
209 1581491034 : GEN xd, zd = (GEN)avma;
210 : long lz;
211 : LOCAL_OVERFLOW;
212 :
213 1581491034 : if (nx == 1) return utoi(x[0] - s);
214 :
215 363978273 : lz = nx+2; (void)new_chunk(lz);
216 363978273 : xd = x + nx;
217 363978273 : *--zd = subll(*--xd, s);
218 363978273 : if (overflow)
219 : for(;;)
220 : {
221 161294364 : *--zd = ((ulong)*--xd) - 1;
222 161294364 : if (*xd) break;
223 : }
224 363978273 : if (xd == x)
225 266743287 : while (*zd == 0) { zd++; lz--; } /* shorten z */
226 : else
227 4734387840 : do *--zd = *--xd; while (xd > x);
228 363978273 : *--zd = evalsigne(1) | evallgefint(lz);
229 363978273 : *--zd = evaltyp(t_INT) | evallg(lz);
230 363978273 : return gc_const((pari_sp)zd, zd);
231 : }
232 :
233 : /* assume x > y */
234 : static GEN
235 3374480082 : subiispec(GEN x, GEN y, long nx, long ny)
236 : {
237 : GEN xd,yd,zd;
238 3374480082 : long lz, i = -2;
239 : LOCAL_OVERFLOW;
240 :
241 3374480082 : if (ny==1) return subiuspec(x,*y,nx);
242 1952515821 : zd = (GEN)avma;
243 1952515821 : lz = nx+2; (void)new_chunk(lz);
244 1952515821 : xd = x + nx;
245 1952515821 : yd = y + ny;
246 1952515821 : zd[-1] = subll(xd[-1], yd[-1]);
247 : #ifdef subllx8
248 2180305820 : for ( ; i-8 > -ny; i-=8)
249 1529467213 : subllx8(xd+i, yd+i, zd+i, overflow);
250 : #endif
251 34074567763 : for ( ; i >= -ny; i--) zd[i] = subllx(xd[i], yd[i]);
252 1952515821 : if (overflow)
253 : for(;;)
254 : {
255 989446716 : zd[i] = uel(xd,i) - 1;
256 989446716 : if (xd[i--]) break;
257 : }
258 1952515821 : if (i>=-nx)
259 4535967207 : for (; i >= -nx; i--) zd[i] = xd[i];
260 : else
261 2343175041 : while (zd[i+1] == 0) { i++; lz--; } /* shorten z */
262 1952515821 : zd += i+1;
263 1952515821 : *--zd = evalsigne(1) | evallgefint(lz);
264 1952515821 : *--zd = evaltyp(t_INT) | evallg(lz);
265 1952515821 : return gc_const((pari_sp)zd, zd);
266 : }
267 :
268 : static void
269 443928775 : roundr_up_ip(GEN x, long l)
270 : {
271 443928775 : long i = l;
272 : for(;;)
273 : {
274 444889258 : if (++uel(x,--i)) break;
275 1287795 : if (i == 2) { x[2] = (long)HIGHBIT; shiftr_inplace(x, 1); break; }
276 : }
277 443928775 : }
278 :
279 : void
280 320683815 : affir(GEN x, GEN y)
281 : {
282 320683815 : const long s = signe(x), ly = lg(y);
283 : long lx, sh, i;
284 :
285 320683815 : if (!s)
286 : {
287 31090587 : y[1] = evalexpo(-bit_accuracy(ly));
288 31090587 : return;
289 : }
290 :
291 289593228 : lx = lgefint(x); sh = bfffo(x[2]);
292 289593228 : y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
293 289593228 : if (sh) {
294 284537598 : if (lx <= ly)
295 : {
296 749544981 : for (i=lx; i<ly; i++) y[i]=0;
297 204026241 : shift_left(y,x,2,lx-1, 0,sh);
298 204026241 : return;
299 : }
300 80511357 : shift_left(y,x,2,ly-1, x[ly],sh);
301 : /* lx > ly: round properly */
302 80511357 : if ((uel(x,ly)<<sh) & HIGHBIT) roundr_up_ip(y, ly);
303 : }
304 : else {
305 5055630 : if (lx <= ly)
306 : {
307 4677891 : for (i=2; i<lx; i++) y[i]=x[i];
308 3725493 : for ( ; i<ly; i++) y[i]=0;
309 1271655 : return;
310 : }
311 9521244 : for (i=2; i<ly; i++) y[i]=x[i];
312 : /* lx > ly: round properly */
313 3783975 : if (uel(x,ly) & HIGHBIT) roundr_up_ip(y, ly);
314 : }
315 : }
316 :
317 : INLINE GEN
318 1274740152 : shiftispec(GEN x, long nx, long n)
319 : {
320 : long ny, i, m;
321 : GEN y, yd;
322 1274740152 : if (!n) return icopyspec(x, nx);
323 :
324 1188217158 : if (n > 0)
325 : {
326 731666439 : GEN z = (GEN)avma;
327 731666439 : long d = dvmdsBIL(n, &m);
328 :
329 731666439 : ny = nx+d; y = new_chunk(ny + 2); yd = y + 2;
330 6716376918 : for ( ; d; d--) *--z = 0;
331 1905753735 : if (!m) for (i=0; i<nx; i++) yd[i]=x[i];
332 : else
333 : {
334 710712384 : const ulong sh = BITS_IN_LONG - m;
335 710712384 : shift_left(yd,x, 0,nx-1, 0,m);
336 710712384 : i = uel(x,0) >> sh;
337 : /* Extend y on the left? */
338 710712384 : if (i) { ny++; y = new_chunk(1); y[2] = i; }
339 : }
340 : }
341 : else
342 : {
343 456550719 : ny = nx - dvmdsBIL(-n, &m);
344 456550719 : if (ny<1) return gen_0;
345 455284326 : y = new_chunk(ny + 2); yd = y + 2;
346 455284326 : if (m) {
347 277514970 : shift_right(yd,x, 0,ny, 0,m);
348 277514970 : if (yd[0] == 0)
349 : {
350 33239721 : if (ny==1) return gc_const((pari_sp)(y+3), gen_0);
351 25581135 : ny--; set_avma((pari_sp)(++y));
352 : }
353 : } else {
354 7412777097 : for (i=0; i<ny; i++) yd[i]=x[i];
355 : }
356 : }
357 1179292179 : y[1] = evalsigne(1)|evallgefint(ny + 2);
358 1179292179 : y[0] = evaltyp(t_INT)|evallg(ny + 2); return y;
359 : }
360 :
361 : GEN
362 48739242 : mantissa2nr(GEN x, long n)
363 : { /*This is a kludge since x is not an integer*/
364 48739242 : long s = signe(x);
365 : GEN y;
366 :
367 48739242 : if(s == 0) return gen_0;
368 48738321 : y = shiftispec(x + 2, lg(x) - 2, n);
369 48738321 : if (signe(y)) setsigne(y, s);
370 48738321 : return y;
371 : }
372 :
373 : GEN
374 2663010 : truncr(GEN x)
375 : {
376 : long d,e,i,s,m;
377 : GEN y;
378 :
379 2663010 : if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gen_0;
380 1149627 : d = nbits2lg(e+1); m = remsBIL(e);
381 1149627 : if (d > lg(x)) pari_err_PREC( "truncr (precision loss in truncation)");
382 :
383 1149624 : y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
384 1149624 : if (++m == BITS_IN_LONG)
385 828 : for (i=2; i<d; i++) y[i]=x[i];
386 : else
387 1149273 : shift_right(y,x, 2,d,0, BITS_IN_LONG - m);
388 1149624 : return y;
389 : }
390 :
391 : /* integral part */
392 : GEN
393 5243496 : floorr(GEN x)
394 : {
395 : long d,e,i,lx,m;
396 : GEN y;
397 :
398 5243496 : if (signe(x) >= 0) return truncr(x);
399 3156192 : if ((e=expo(x)) < 0) return gen_m1;
400 2639037 : d = nbits2lg(e+1); m = remsBIL(e);
401 2639037 : lx=lg(x); if (d>lx) pari_err_PREC( "floorr (precision loss in truncation)");
402 2639037 : y = new_chunk(d);
403 2639037 : if (++m == BITS_IN_LONG)
404 : {
405 501 : for (i=2; i<d; i++) y[i]=x[i];
406 210 : i=d; while (i<lx && !x[i]) i++;
407 174 : if (i==lx) goto END;
408 : }
409 : else
410 : {
411 2638863 : shift_right(y,x, 2,d,0, BITS_IN_LONG - m);
412 2638863 : if (uel(x,d-1)<<m == 0)
413 : {
414 313536 : i=d; while (i<lx && !x[i]) i++;
415 84048 : if (i==lx) goto END;
416 : }
417 : }
418 : /* set y:=y+1 */
419 2582442 : for (i=d-1; i>=2; i--) { uel(y,i)++; if (y[i]) goto END; }
420 0 : y=new_chunk(1); y[2]=1; d++;
421 2639037 : END:
422 2639037 : y[1] = evalsigne(-1) | evallgefint(d);
423 2639037 : y[0] = evaltyp(t_INT) | evallg(d); return y;
424 : }
425 :
426 : INLINE int
427 4061637441 : cmpiispec(GEN x, GEN y, long lx, long ly)
428 : {
429 : long i;
430 4061637441 : if (lx < ly) return -1;
431 3770557530 : if (lx > ly) return 1;
432 3750747813 : i = 0; while (i<lx && x[i]==y[i]) i++;
433 3257442813 : if (i==lx) return 0;
434 3036745578 : return (uel(x,i) > uel(y,i))? 1: -1;
435 : }
436 :
437 : INLINE int
438 197853081 : equaliispec(GEN x, GEN y, long lx, long ly)
439 : {
440 : long i;
441 197853081 : if (lx != ly) return 0;
442 361809777 : i = ly-1; while (i>=0 && x[i]==y[i]) i--;
443 197770506 : return i < 0;
444 : }
445 :
446 : /***********************************************************************/
447 : /** **/
448 : /** MULTIPLICATION **/
449 : /** **/
450 : /***********************************************************************/
451 : /* assume ny > 0 */
452 : INLINE GEN
453 4764810786 : muluispec(ulong x, GEN y, long ny)
454 : {
455 4764810786 : GEN yd, z = (GEN)avma;
456 4764810786 : long lz = ny+3;
457 : LOCAL_HIREMAINDER;
458 :
459 4764810786 : (void)new_chunk(lz);
460 4764810786 : yd = y + ny; *--z = mulll(x, *--yd);
461 15358337733 : while (yd > y) *--z = addmul(x,*--yd);
462 4764810786 : if (hiremainder) *--z = hiremainder; else lz--;
463 4764810786 : *--z = evalsigne(1) | evallgefint(lz);
464 4764810786 : *--z = evaltyp(t_INT) | evallg(lz);
465 4764810786 : return gc_const((pari_sp)z, z);
466 : }
467 :
468 : /* a + b*|Y| */
469 : GEN
470 0 : addumului(ulong a, ulong b, GEN Y)
471 : {
472 : GEN yd,y,z;
473 : long ny,lz;
474 : LOCAL_HIREMAINDER;
475 : LOCAL_OVERFLOW;
476 :
477 0 : if (!b || !signe(Y)) return utoi(a);
478 :
479 0 : y = LIMBS(Y); z = (GEN)avma;
480 0 : ny = NLIMBS(Y);
481 0 : lz = ny+3;
482 :
483 0 : (void)new_chunk(lz);
484 0 : yd = y + ny; *--z = addll(a, mulll(b, *--yd));
485 0 : if (overflow) hiremainder++; /* can't overflow */
486 0 : while (yd > y) *--z = addmul(b,*--yd);
487 0 : if (hiremainder) *--z = hiremainder; else lz--;
488 0 : *--z = evalsigne(1) | evallgefint(lz);
489 0 : *--z = evaltyp(t_INT) | evallg(lz);
490 0 : return gc_const((pari_sp)z, z);
491 : }
492 :
493 : /***********************************************************************/
494 : /** **/
495 : /** DIVISION **/
496 : /** **/
497 : /***********************************************************************/
498 :
499 : ulong
500 1408935738 : umodiu(GEN y, ulong x)
501 : {
502 1408935738 : long sy=signe(y),ly,i;
503 : ulong xi;
504 : LOCAL_HIREMAINDER;
505 :
506 1408935738 : if (!x) pari_err_INV("umodiu",gen_0);
507 1408935738 : if (!sy) return 0;
508 1104335865 : ly = lgefint(y);
509 1104335865 : if (x <= uel(y,2))
510 : {
511 333160149 : hiremainder=0;
512 333160149 : if (ly==3)
513 : {
514 302836191 : hiremainder=uel(y,2)%x;
515 302836191 : if (!hiremainder) return 0;
516 254208687 : return (sy > 0)? hiremainder: x - hiremainder;
517 : }
518 : }
519 : else
520 : {
521 771175716 : if (ly==3) return (sy > 0)? uel(y,2): x - uel(y,2);
522 97580139 : hiremainder=uel(y,2); ly--; y++;
523 : }
524 127904097 : xi = get_Fl_red(x);
525 914849400 : for (i=2; i<ly; i++) (void)divll_pre(y[i],x,xi);
526 127904097 : if (!hiremainder) return 0;
527 121929168 : return (sy > 0)? hiremainder: x - hiremainder;
528 : }
529 :
530 : /* return |y| \/ x */
531 : GEN
532 271278447 : absdiviu_rem(GEN y, ulong x, ulong *rem)
533 : {
534 : long ly,i;
535 : GEN z;
536 : ulong xi;
537 : LOCAL_HIREMAINDER;
538 :
539 271278447 : if (!x) pari_err_INV("absdiviu_rem",gen_0);
540 271278447 : if (!signe(y)) { *rem = 0; return gen_0; }
541 :
542 254306082 : ly = lgefint(y);
543 254306082 : if (x <= uel(y,2))
544 : {
545 227138949 : hiremainder=0;
546 227138949 : if (ly==3)
547 : {
548 203531601 : z = cgetipos(3);
549 203531601 : z[2] = divll(uel(y,2),x);
550 203531601 : *rem = hiremainder; return z;
551 : }
552 : }
553 : else
554 : {
555 27167133 : if (ly==3) { *rem = uel(y,2); return gen_0; }
556 6807135 : hiremainder = uel(y,2); ly--; y++;
557 : }
558 30414483 : xi = get_Fl_red(x);
559 30414483 : z = cgetipos(ly);
560 165450273 : for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x,xi);
561 30414483 : *rem = hiremainder; return z;
562 : }
563 :
564 : GEN
565 65137269 : divis_rem(GEN y, long x, long *rem)
566 : {
567 65137269 : long sy=signe(y),ly,s,i;
568 : GEN z;
569 : ulong xi;
570 : LOCAL_HIREMAINDER;
571 :
572 65137269 : if (!x) pari_err_INV("divis_rem",gen_0);
573 65137269 : if (!sy) { *rem=0; return gen_0; }
574 45990102 : if (x<0) { s = -sy; x = -x; } else s = sy;
575 :
576 45990102 : ly = lgefint(y);
577 45990102 : if ((ulong)x <= uel(y,2))
578 : {
579 31730610 : hiremainder=0;
580 31730610 : if (ly==3)
581 : {
582 31423386 : z = cgeti(3); z[1] = evallgefint(3) | evalsigne(s);
583 31423386 : z[2] = divll(uel(y,2),x);
584 31423386 : if (sy<0) hiremainder = - ((long)hiremainder);
585 31423386 : *rem = (long)hiremainder; return z;
586 : }
587 : }
588 : else
589 : {
590 14259492 : if (ly==3) { *rem = itos(y); return gen_0; }
591 257631 : hiremainder = uel(y,2); ly--; y++;
592 : }
593 564855 : xi = get_Fl_red(x);
594 564855 : z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
595 2968047 : for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x,xi);
596 564855 : if (sy<0) hiremainder = - ((long)hiremainder);
597 564855 : *rem = (long)hiremainder; return z;
598 : }
599 :
600 : GEN
601 721872 : divis(GEN y, long x)
602 : {
603 721872 : long sy=signe(y),ly,s,i;
604 : ulong xi;
605 : GEN z;
606 : LOCAL_HIREMAINDER;
607 :
608 721872 : if (!x) pari_err_INV("divis",gen_0);
609 721872 : if (!sy) return gen_0;
610 721836 : if (x<0) { s = -sy; x = -x; } else s = sy;
611 :
612 721836 : ly = lgefint(y);
613 721836 : if ((ulong)x <= uel(y,2))
614 : {
615 712812 : hiremainder=0;
616 712812 : if (ly==3)
617 : {
618 643341 : z = cgeti(3); z[1] = evallgefint(3) | evalsigne(s);
619 643341 : z[2] = divll(y[2],x);
620 643341 : return z;
621 : }
622 : }
623 : else
624 : {
625 9024 : if (ly==3) return gen_0;
626 8790 : hiremainder=y[2]; ly--; y++;
627 : }
628 78261 : xi = get_Fl_red(x);
629 78261 : z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
630 593268 : for (i=2; i<ly; i++) z[i]=divll_pre(y[i],x, xi);
631 78261 : return z;
632 : }
633 :
634 : GEN
635 129584295 : divrr(GEN x, GEN y)
636 : {
637 129584295 : long sx=signe(x), sy=signe(y), lx,ly,lr,e,i,j;
638 : ulong y0,y1;
639 : GEN r, r1;
640 :
641 129584295 : if (!sy) pari_err_INV("divrr",y);
642 129584295 : e = expo(x) - expo(y);
643 129584295 : if (!sx) return real_0_bit(e);
644 129226908 : if (sy<0) sx = -sx;
645 :
646 129226908 : lx=lg(x); ly=lg(y);
647 129226908 : if (ly==3)
648 : {
649 23523993 : ulong k = x[2], l = (lx>3)? x[3]: 0;
650 : LOCAL_HIREMAINDER;
651 23523993 : if (k < uel(y,2)) e--;
652 : else
653 : {
654 6879498 : l >>= 1; if (k&1) l |= HIGHBIT;
655 6879498 : k >>= 1;
656 : }
657 23523993 : hiremainder = k; k = divll(l,y[2]);
658 23523993 : if (hiremainder > (uel(y,2) >> 1) && !++k) { k = HIGHBIT; e++; }
659 23523993 : r = cgetg(3, t_REAL);
660 23523993 : r[1] = evalsigne(sx) | evalexpo(e);
661 23523993 : r[2] = k; return r;
662 : }
663 :
664 105702915 : lr = minss(lx,ly); r = new_chunk(lr);
665 105702915 : r1 = r-1;
666 749262225 : r1[1] = 0; for (i=2; i<lr; i++) r1[i]=x[i];
667 105702915 : r1[lr] = (lx>ly)? x[lr]: 0;
668 105702915 : y0 = y[2]; y1 = y[3];
669 854965140 : for (i=0; i<lr-1; i++)
670 : { /* r1 = r + (i-1), OK up to r1[2] (accesses at most r[lr]) */
671 : ulong k, qp;
672 : LOCAL_HIREMAINDER;
673 : LOCAL_OVERFLOW;
674 :
675 749262225 : if (uel(r1,1) == y0) { qp = ULONG_MAX; k = addll(y0,r1[2]); }
676 : else
677 : {
678 747638946 : if (uel(r1,1) > y0) /* can't happen if i=0 */
679 : {
680 0 : GEN y1 = y+1;
681 0 : j = lr-i; r1[j] = subll(r1[j],y1[j]);
682 0 : for (j--; j>0; j--) r1[j] = subllx(r1[j],y1[j]);
683 0 : j=i; do uel(r,--j)++; while (j && !uel(r,j));
684 : }
685 747638946 : hiremainder = r1[1]; overflow = 0;
686 747638946 : qp = divll(r1[2],y0); k = hiremainder;
687 : }
688 749262225 : j = lr-i+1;
689 749262225 : if (!overflow)
690 : {
691 : long k3, k4;
692 747971937 : k3 = mulll(qp,y1);
693 747971937 : if (j == 3) /* i = lr - 2 maximal, r1[3] undefined -> 0 */
694 105635406 : k4 = subll(hiremainder,k);
695 : else
696 : {
697 642336531 : k3 = subll(k3, r1[3]);
698 642336531 : k4 = subllx(hiremainder,k);
699 : }
700 990968837 : while (!overflow && k4) { qp--; k3 = subll(k3,y1); k4 = subllx(k4,y0); }
701 : }
702 749262225 : if (j<ly) (void)mulll(qp,y[j]); else { hiremainder = 0 ; j = ly; }
703 4961831307 : for (j--; j>1; j--)
704 : {
705 4212569082 : r1[j] = subll(r1[j], addmul(qp,y[j]));
706 4212569082 : hiremainder += overflow;
707 : }
708 749262225 : if (uel(r1,1) != hiremainder)
709 : {
710 609333 : if (uel(r1,1) < hiremainder)
711 : {
712 609333 : qp--;
713 609333 : j = lr-i-(lr-i>=ly); r1[j] = addll(r1[j], y[j]);
714 3398619 : for (j--; j>1; j--) r1[j] = addllx(r1[j], y[j]);
715 : }
716 : else
717 : {
718 0 : r1[1] -= hiremainder;
719 0 : while (r1[1])
720 : {
721 0 : qp++; if (!qp) { j=i; do uel(r,--j)++; while (j && !r[j]); }
722 0 : j = lr-i-(lr-i>=ly); r1[j] = subll(r1[j],y[j]);
723 0 : for (j--; j>1; j--) r1[j] = subllx(r1[j],y[j]);
724 0 : r1[1] -= overflow;
725 : }
726 : }
727 : }
728 749262225 : *++r1 = qp;
729 : }
730 : /* i = lr-1 */
731 : /* round correctly */
732 105702915 : if (uel(r1,1) > (y0>>1))
733 : {
734 51912834 : j=i; do uel(r,--j)++; while (j && !r[j]);
735 : }
736 749262225 : r1 = r-1; for (j=i; j>=2; j--) r[j]=r1[j];
737 105702915 : if (r[0] == 0) e--;
738 45957585 : else if (r[0] == 1) { shift_right(r,r, 2,lr, 1,1); }
739 : else { /* possible only when rounding up to 0x2 0x0 ... */
740 6 : r[2] = (long)HIGHBIT; e++;
741 : }
742 105702915 : r[0] = evaltyp(t_REAL)|evallg(lr);
743 105702915 : r[1] = evalsigne(sx) | evalexpo(e);
744 105702915 : return r;
745 : }
746 :
747 : GEN
748 113307054 : divri(GEN x, GEN y)
749 : {
750 113307054 : long lx, s = signe(y);
751 : pari_sp av;
752 : GEN z;
753 :
754 113307054 : if (!s) pari_err_INV("divri",y);
755 113307054 : if (!signe(x)) return real_0_bit(expo(x) - expi(y));
756 113137116 : if (!is_bigint(y)) {
757 88929720 : GEN z = divru(x, y[2]);
758 88929720 : if (s < 0) togglesign(z);
759 88929720 : return z;
760 : }
761 24207396 : lx = lg(x); z = cgetg(lx, t_REAL); av = avma;
762 24207396 : affrr(divrr(x, itor(y, lg2prec(lx+1))), z);
763 24207396 : return gc_const(av, z);
764 : }
765 :
766 : /* Integer division x / y: such that sign(r) = sign(x)
767 : * if z = ONLY_REM return remainder, otherwise return quotient
768 : * if z != NULL set *z to remainder
769 : * *z is the last object on stack (and thus can be disposed of with cgiv
770 : * instead of gerepile)
771 : * If *z is zero, we put gen_0 here and no copy.
772 : * space needed: lx + ly */
773 : GEN
774 1667385579 : dvmdii(GEN x, GEN y, GEN *z)
775 : {
776 1667385579 : long sx = signe(x), sy = signe(y);
777 1667385579 : long lx, ly = lgefint(y), lz, i, j, sh, lq, lr;
778 : pari_sp av;
779 : ulong y0,y0i,y1, *xd,*rd,*qd;
780 : GEN q, r, r1;
781 :
782 1667385579 : if (!sx)
783 : {
784 52978335 : if (ly < 3) pari_err_INV("dvmdii",gen_0);
785 52978332 : if (!z || z == ONLY_REM) return gen_0;
786 32209836 : *z=gen_0; return gen_0;
787 : }
788 1614407244 : if (ly <= 3)
789 : {
790 : ulong rem;
791 647406522 : if (ly < 3) pari_err_INV("dvmdii",gen_0);
792 647406522 : if (z == ONLY_REM)
793 : {
794 445431168 : rem = umodiu(x,uel(y,2));
795 445431168 : if (!rem) return gen_0;
796 402928149 : return (sx < 0)? utoineg(uel(y,2) - rem): utoipos(rem);
797 : }
798 201975354 : q = absdiviu_rem(x, uel(y,2), &rem);
799 201975354 : if (sx != sy) togglesign(q);
800 201975354 : if (!z) return q;
801 198744756 : if (!rem) *z = gen_0;
802 54402573 : else *z = sx < 0? utoineg(rem): utoipos(rem);
803 198744756 : return q;
804 : }
805 967000722 : lx=lgefint(x);
806 967000722 : lz=lx-ly;
807 967000722 : if (lz <= 0)
808 : {
809 440169801 : if (lz == 0)
810 : {
811 334241115 : for (i=2; i<lx; i++)
812 333614739 : if (x[i] != y[i])
813 : {
814 317271111 : if (uel(x,i) > uel(y,i)) goto DIVIDE;
815 44635410 : goto TRIVIAL;
816 : }
817 626376 : if (z == ONLY_REM) return gen_0;
818 65793 : if (z) *z = gen_0;
819 65793 : if (sx < 0) sy = -sy;
820 65793 : return stoi(sy);
821 : }
822 122272314 : TRIVIAL:
823 166907724 : if (z == ONLY_REM) return icopy(x);
824 2103972 : if (z) *z = icopy(x);
825 2103972 : return gen_0;
826 : }
827 526830921 : DIVIDE: /* quotient is nonzero */
828 799466622 : av=avma; if (sx<0) sy = -sy;
829 799466622 : r1 = new_chunk(lx); sh = bfffo(y[2]);
830 799466622 : if (sh)
831 : { /* normalize so that highbit(y) = 1 (shift left x and y by sh bits)*/
832 790674672 : const ulong m = BITS_IN_LONG - sh;
833 790674672 : r = new_chunk(ly);
834 790674672 : shift_left(r, y,2,ly-1, 0,sh); y = r;
835 790674672 : shift_left(r1,x,2,lx-1, 0,sh);
836 790674672 : r1[1] = uel(x,2) >> m;
837 : }
838 : else
839 : {
840 91013064 : r1[1] = 0; for (j=2; j<lx; j++) r1[j] = x[j];
841 : }
842 799466622 : x = r1;
843 799466622 : y0 = y[2]; y0i = get_Fl_red(y0);
844 799466622 : y1 = y[3];
845 3279770121 : for (i=0; i<=lz; i++)
846 : { /* r1 = x + i */
847 : ulong k, qp;
848 : LOCAL_HIREMAINDER;
849 : LOCAL_OVERFLOW;
850 :
851 2480303499 : if (uel(r1,1) == y0)
852 : {
853 48366 : qp = ULONG_MAX; k = addll(y0,r1[2]);
854 : }
855 : else
856 : {
857 2480255133 : hiremainder = r1[1]; overflow = 0;
858 2480255133 : qp = divll_pre(r1[2],y0,y0i); k = hiremainder;
859 : }
860 2480303499 : if (!overflow)
861 : {
862 2480254226 : long k3 = subll(mulll(qp,y1), r1[3]);
863 2480254226 : long k4 = subllx(hiremainder,k);
864 3027197883 : while (!overflow && k4) { qp--; k3 = subll(k3,y1); k4 = subllx(k4,y0); }
865 : }
866 2480303499 : hiremainder = 0; j = ly;
867 64777576695 : for (j--; j>1; j--)
868 : {
869 62297273196 : r1[j] = subll(r1[j], addmul(qp,y[j]));
870 62297273196 : hiremainder += overflow;
871 : }
872 2480303499 : if (uel(r1,1) < hiremainder)
873 : {
874 5918224 : qp--;
875 5918224 : j = ly-1; r1[j] = addll(r1[j],y[j]);
876 31198063 : for (j--; j>1; j--) r1[j] = addllx(r1[j],y[j]);
877 : }
878 2480303499 : *++r1 = qp;
879 : }
880 :
881 799466622 : lq = lz+2;
882 799466622 : if (!z)
883 : {
884 2811741 : qd = (ulong*)av;
885 2811741 : xd = (ulong*)(x + lq);
886 2811741 : if (x[1]) { lz++; lq++; }
887 34975005 : while (lz--) *--qd = *--xd;
888 2811741 : *--qd = evalsigne(sy) | evallgefint(lq);
889 2811741 : *--qd = evaltyp(t_INT) | evallg(lq);
890 2811741 : return gc_const((pari_sp)qd, (GEN)qd);
891 : }
892 :
893 888466974 : j=lq; while (j<lx && !x[j]) j++;
894 796654881 : lz = lx-j;
895 796654881 : if (z == ONLY_REM)
896 : {
897 517076511 : if (lz==0) return gc_const(av, gen_0);
898 507768447 : rd = (ulong*)av; lr = lz+2;
899 507768447 : xd = (ulong*)(x + lx);
900 539956428 : if (!sh) while (lz--) *--rd = *--xd;
901 : else
902 : { /* shift remainder right by sh bits */
903 499746489 : const ulong shl = BITS_IN_LONG - sh;
904 : ulong l;
905 499746489 : xd--;
906 1503056124 : while (--lz) /* fill r[3..] */
907 : {
908 1003309635 : l = *xd >> sh;
909 1003309635 : *--rd = l | (*--xd << shl);
910 : }
911 499746489 : l = *xd >> sh;
912 499746489 : if (l) *--rd = l; else lr--;
913 : }
914 507768447 : *--rd = evalsigne(sx) | evallgefint(lr);
915 507768447 : *--rd = evaltyp(t_INT) | evallg(lr);
916 507768447 : return gc_const((pari_sp)rd, (GEN)rd);
917 : }
918 :
919 279578370 : lr = lz+2;
920 279578370 : rd = NULL; /* gcc -Wall */
921 279578370 : if (lz)
922 : { /* non zero remainder: initialize rd */
923 274971939 : xd = (ulong*)(x + lx);
924 274971939 : if (!sh)
925 : {
926 569994 : rd = (ulong*)avma; (void)new_chunk(lr);
927 5779566 : while (lz--) *--rd = *--xd;
928 : }
929 : else
930 : { /* shift remainder right by sh bits */
931 274401945 : const ulong shl = BITS_IN_LONG - sh;
932 : ulong l;
933 274401945 : rd = (ulong*)x; /* overwrite shifted y */
934 274401945 : xd--;
935 1224709044 : while (--lz)
936 : {
937 950307099 : l = *xd >> sh;
938 950307099 : *--rd = l | (*--xd << shl);
939 : }
940 274401945 : l = *xd >> sh;
941 274401945 : if (l) *--rd = l; else lr--;
942 : }
943 274971939 : *--rd = evalsigne(sx) | evallgefint(lr);
944 274971939 : *--rd = evaltyp(t_INT) | evallg(lr);
945 274971939 : rd += lr;
946 : }
947 279578370 : qd = (ulong*)av;
948 279578370 : xd = (ulong*)(x + lq);
949 279578370 : if (x[1]) lq++;
950 874857174 : j = lq-2; while (j--) *--qd = *--xd;
951 279578370 : *--qd = evalsigne(sy) | evallgefint(lq);
952 279578370 : *--qd = evaltyp(t_INT) | evallg(lq);
953 279578370 : q = (GEN)qd;
954 279578370 : if (lr==2) *z = gen_0;
955 : else
956 : { /* rd has been properly initialized: we had lz > 0 */
957 1874667801 : while (lr--) *--qd = *--rd;
958 274971939 : *z = (GEN)qd;
959 : }
960 279578370 : return gc_const((pari_sp)qd, q);
961 : }
962 :
963 : /* Montgomery reduction.
964 : * N has k words, assume T >= 0 has less than 2k.
965 : * Return res := T / B^k mod N, where B = 2^BIL
966 : * such that 0 <= res < T/B^k + N and res has less than k words */
967 : GEN
968 39049344 : red_montgomery(GEN T, GEN N, ulong inv)
969 : {
970 : pari_sp av;
971 : GEN Te, Td, Ne, Nd, scratch;
972 39049344 : ulong i, j, m, t, d, k = NLIMBS(N);
973 : int carry;
974 : LOCAL_HIREMAINDER;
975 : LOCAL_OVERFLOW;
976 :
977 39049344 : if (k == 0) return gen_0;
978 39049344 : d = NLIMBS(T); /* <= 2*k */
979 39049344 : if (d == 0) return gen_0;
980 : #ifdef DEBUG
981 : if (d > 2*k) pari_err_BUG("red_montgomery");
982 : #endif
983 39049335 : if (k == 1)
984 : { /* as below, special cased for efficiency */
985 163341 : ulong n = uel(N,2);
986 163341 : if (d == 1) {
987 163194 : hiremainder = uel(T,2);
988 163194 : m = hiremainder * inv;
989 163194 : (void)addmul(m, n); /* t + m*n = 0 */
990 163194 : return utoi(hiremainder);
991 : } else { /* d = 2 */
992 147 : hiremainder = uel(T,3);
993 147 : m = hiremainder * inv;
994 147 : (void)addmul(m, n); /* t + m*n = 0 */
995 147 : t = addll(hiremainder, uel(T,2));
996 147 : if (overflow) t -= n; /* t > n doesn't fit in 1 word */
997 147 : return utoi(t);
998 : }
999 : }
1000 : /* assume k >= 2 */
1001 38885994 : av = avma; scratch = new_chunk(k<<1); /* >= k + 2: result fits */
1002 :
1003 : /* copy T to scratch space (pad with zeroes to 2k words) */
1004 38885994 : Td = (GEN)av;
1005 38885994 : Te = T + (d+2);
1006 855079152 : for (i=0; i < d ; i++) *--Td = *--Te;
1007 64906878 : for ( ; i < (k<<1); i++) *--Td = 0;
1008 :
1009 38885994 : Te = (GEN)av; /* 1 beyond end of current T mantissa (in scratch) */
1010 38885994 : Ne = N + k+2; /* 1 beyond end of N mantissa */
1011 :
1012 38885994 : carry = 0;
1013 459993015 : for (i=0; i<k; i++) /* set T := T/B nod N, k times */
1014 : {
1015 421107021 : Td = Te; /* one beyond end of (new) T mantissa */
1016 421107021 : Nd = Ne;
1017 421107021 : hiremainder = *--Td;
1018 421107021 : m = hiremainder * inv; /* solve T + m N = O(B) */
1019 :
1020 : /* set T := (T + mN) / B */
1021 421107021 : Te = Td;
1022 421107021 : (void)addmul(m, *--Nd); /* = 0 */
1023 6856511649 : for (j=1; j<k; j++)
1024 : {
1025 6435404628 : t = addll(addmul(m, *--Nd), *--Td);
1026 6435404628 : *Td = t;
1027 6435404628 : hiremainder += overflow;
1028 : }
1029 421107021 : t = addll(hiremainder, *--Td); *Td = t + carry;
1030 421107021 : carry = (overflow || (carry && *Td == 0));
1031 : }
1032 38885994 : if (carry)
1033 : { /* Td > N overflows (k+1 words), set Td := Td - N */
1034 373602 : Td = Te;
1035 373602 : Nd = Ne;
1036 373602 : t = subll(*--Td, *--Nd); *Td = t;
1037 6979275 : while (Td > scratch) { t = subllx(*--Td, *--Nd); *Td = t; }
1038 : }
1039 :
1040 : /* copy result */
1041 38885994 : Td = (GEN)av;
1042 42633090 : while (*scratch == 0 && Te > scratch) scratch++; /* strip leading 0s */
1043 456245919 : while (Te > scratch) *--Td = *--Te;
1044 38885994 : k = (GEN)av - Td; if (!k) return gc_const(av, gen_0);
1045 38885994 : k += 2;
1046 38885994 : *--Td = evalsigne(1) | evallgefint(k);
1047 38885994 : *--Td = evaltyp(t_INT) | evallg(k);
1048 : #ifdef DEBUG
1049 : {
1050 : long l = NLIMBS(N), s = BITS_IN_LONG*l;
1051 : GEN R = int2n(s);
1052 : GEN res = remii(mulii(T, Fp_inv(R, N)), N);
1053 : if (k > lgefint(N)
1054 : || !equalii(remii(Td,N),res)
1055 : || cmpii(Td, addii(shifti(T, -s), N)) >= 0) pari_err_BUG("red_montgomery");
1056 : }
1057 : #endif
1058 38885994 : return gc_const((pari_sp)Td, Td);
1059 : }
1060 :
1061 : /* EXACT INTEGER DIVISION */
1062 :
1063 : /* assume xy>0, the division is exact and y is odd. Destroy x */
1064 : static GEN
1065 29497869 : diviuexact_i(GEN x, ulong y)
1066 : {
1067 : long i, lz, lx;
1068 : ulong q, yinv;
1069 : GEN z, z0, x0, x0min;
1070 :
1071 29497869 : if (y == 1) return icopy(x);
1072 23715120 : lx = lgefint(x);
1073 23715120 : if (lx == 3)
1074 : {
1075 853158 : q = uel(x,2) / y;
1076 853158 : if (!q) pari_err_OP("exact division", x, utoi(y));
1077 853158 : return utoipos(q);
1078 : }
1079 22861962 : yinv = invmod2BIL(y);
1080 22861962 : lz = (y <= uel(x,2)) ? lx : lx-1;
1081 22861962 : z = new_chunk(lz);
1082 22861962 : z0 = z + lz;
1083 22861962 : x0 = x + lx; x0min = x + lx-lz+2;
1084 :
1085 82713696 : while (x0 > x0min)
1086 : {
1087 59851734 : *--z0 = q = yinv*uel(--x0,0); /* i-th quotient */
1088 59851734 : if (!q) continue;
1089 : /* x := x - q * y */
1090 : { /* update neither lowest word (could set it to 0) nor highest ones */
1091 59324232 : GEN x1 = x0 - 1;
1092 : LOCAL_HIREMAINDER;
1093 59324232 : (void)mulll(q,y);
1094 59324232 : if (hiremainder)
1095 : {
1096 47602542 : if (uel(x1,0) < hiremainder)
1097 : {
1098 138243 : uel(x1,0) -= hiremainder;
1099 140247 : do uel(--x1,0)--; while (uel(x1,0) == ULONG_MAX);
1100 : }
1101 : else
1102 47464299 : uel(x1,0) -= hiremainder;
1103 : }
1104 : }
1105 : }
1106 22861962 : i=2; while(!z[i]) i++;
1107 22861962 : z += i-2; lz -= i-2;
1108 22861962 : z[0] = evaltyp(t_INT)|evallg(lz);
1109 22861962 : z[1] = evalsigne(1)|evallg(lz);
1110 22861962 : if (lz == 2) pari_err_OP("exact division", x, utoi(y));
1111 22861962 : return gc_const((pari_sp)z, z);
1112 : }
1113 :
1114 : /* assume y != 0 and the division is exact */
1115 : GEN
1116 22234683 : diviuexact(GEN x, ulong y)
1117 : {
1118 : pari_sp av;
1119 22234683 : long lx, vy, s = signe(x);
1120 : GEN z;
1121 :
1122 22234683 : if (!s) return gen_0;
1123 21388224 : if (y == 1) return icopy(x);
1124 18393822 : lx = lgefint(x);
1125 18393822 : if (lx == 3) {
1126 14227632 : ulong q = uel(x,2) / y;
1127 14227632 : if (!q) pari_err_OP("exact division", x, utoi(y));
1128 14227632 : return (s > 0)? utoipos(q): utoineg(q);
1129 : }
1130 4166190 : av = avma; (void)new_chunk(lx); vy = vals(y);
1131 4166190 : if (vy) {
1132 1641192 : y >>= vy;
1133 1641192 : if (y == 1) { set_avma(av); return shifti(x, -vy); }
1134 800880 : x = shifti(x, -vy);
1135 800880 : if (lx == 3) {
1136 0 : ulong q = uel(x,2) / y;
1137 0 : set_avma(av);
1138 0 : if (!q) pari_err_OP("exact division", x, utoi(y));
1139 0 : return (s > 0)? utoipos(q): utoineg(q);
1140 : }
1141 2524998 : } else x = icopy(x);
1142 3325878 : set_avma(av);
1143 3325878 : z = diviuexact_i(x, y);
1144 3325878 : setsigne(z, s); return z;
1145 : }
1146 :
1147 : /* Find z such that x=y*z, knowing that y | x (unchecked)
1148 : * Method: y0 z0 = x0 mod B = 2^BITS_IN_LONG ==> z0 = 1/y0 mod B.
1149 : * Set x := (x - z0 y) / B, updating only relevant words, and repeat */
1150 : GEN
1151 377433558 : diviiexact(GEN x, GEN y)
1152 : {
1153 377433558 : long lx, ly, lz, vy, i, ii, sx = signe(x), sy = signe(y);
1154 : pari_sp av;
1155 : ulong y0inv,q;
1156 : GEN z;
1157 :
1158 377433558 : if (!sy) pari_err_INV("diviiexact",gen_0);
1159 377433558 : if (!sx) return gen_0;
1160 309390417 : lx = lgefint(x);
1161 309390417 : if (lx == 3) {
1162 246217476 : q = uel(x,2) / uel(y,2);
1163 246217476 : if (!q) pari_err_OP("exact division", x, y);
1164 246217476 : return (sx+sy) ? utoipos(q): utoineg(q);
1165 : }
1166 63172941 : vy = vali(y); av = avma;
1167 63172941 : (void)new_chunk(lx); /* enough room for z */
1168 63172941 : if (vy)
1169 : { /* make y odd */
1170 32306574 : y = shifti(y,-vy);
1171 32306574 : x = shifti(x,-vy); lx = lgefint(x);
1172 : }
1173 30866367 : else x = icopy(x); /* necessary because we destroy x */
1174 63172941 : set_avma(av); /* will erase our x,y when exiting */
1175 : /* now y is odd */
1176 63172941 : ly = lgefint(y);
1177 63172941 : if (ly == 3)
1178 : {
1179 26171991 : z = diviuexact_i(x,uel(y,2)); /* x != 0 */
1180 26171991 : setsigne(z, (sx+sy)? 1: -1); return z;
1181 : }
1182 37000950 : y0inv = invmod2BIL(y[ly-1]);
1183 58361577 : i=2; while (i<ly && y[i]==x[i]) i++;
1184 37000950 : lz = (i==ly || uel(y,i) < uel(x,i)) ? lx-ly+3 : lx-ly+2;
1185 37000950 : z = new_chunk(lz);
1186 :
1187 37000950 : y += ly - 1; /* now y[-i] = i-th word of y */
1188 173691783 : for (ii=lx-1,i=lz-1; i>=2; i--,ii--)
1189 : {
1190 : long limj;
1191 : LOCAL_HIREMAINDER;
1192 : LOCAL_OVERFLOW;
1193 :
1194 136690833 : z[i] = q = y0inv*uel(x,ii); /* i-th quotient */
1195 136690833 : if (!q) continue;
1196 :
1197 : /* x := x - q * y */
1198 136559235 : (void)mulll(q,y[0]); limj = maxss(lx - lz, ii+3-ly);
1199 : { /* update neither lowest word (could set it to 0) nor highest ones */
1200 136559235 : GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj;
1201 2318506317 : for (; x0 >= xlim; x0--, y0--)
1202 : {
1203 2181947082 : *x0 = subll(*x0, addmul(q,*y0));
1204 2181947082 : hiremainder += overflow;
1205 : }
1206 136559235 : if (hiremainder && limj != lx - lz)
1207 : {
1208 72305268 : if ((ulong)*x0 < hiremainder)
1209 : {
1210 838125 : *x0 -= hiremainder;
1211 838125 : do (*--x0)--; while ((ulong)*x0 == ULONG_MAX);
1212 : }
1213 : else
1214 71467143 : *x0 -= hiremainder;
1215 : }
1216 : }
1217 : }
1218 37000950 : i=2; while(!z[i]) i++;
1219 37000950 : z += i-2; lz -= (i-2);
1220 37000950 : z[0] = evaltyp(t_INT)|evallg(lz);
1221 37000950 : z[1] = evalsigne((sx+sy)? 1: -1) | evallg(lz);
1222 37000950 : if (lz == 2) pari_err_OP("exact division", x, y);
1223 37000950 : return gc_const((pari_sp)z, z);
1224 : }
1225 :
1226 : /* assume yz != and yz | x */
1227 : GEN
1228 149130 : diviuuexact(GEN x, ulong y, ulong z)
1229 : {
1230 : long tmp[4];
1231 : ulong t;
1232 : LOCAL_HIREMAINDER;
1233 149130 : t = mulll(y, z);
1234 149130 : if (!hiremainder) return diviuexact(x, t);
1235 0 : tmp[0] = evaltyp(t_INT)|_evallg(4);
1236 0 : tmp[1] = evalsigne(1)|evallgefint(4);
1237 0 : tmp[2] = hiremainder;
1238 0 : tmp[3] = t;
1239 0 : return diviiexact(x, tmp);
1240 : }
1241 :
1242 : /********************************************************************/
1243 : /** **/
1244 : /** INTEGER MULTIPLICATION (BASECASE) **/
1245 : /** **/
1246 : /********************************************************************/
1247 : /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
1248 : INLINE GEN
1249 5158520793 : muliispec_basecase(GEN x, GEN y, long nx, long ny)
1250 : {
1251 : GEN z2e,z2d,yd,xd,ye,zd;
1252 : long p1,lz;
1253 : LOCAL_HIREMAINDER;
1254 :
1255 5158520793 : if (ny == 1) return muluispec((ulong)*y, x, nx);
1256 1127779470 : if (ny == 0) return gen_0;
1257 1126559097 : zd = (GEN)avma; lz = nx+ny+2;
1258 1126559097 : (void)new_chunk(lz);
1259 1126559097 : xd = x + nx;
1260 1126559097 : yd = y + ny;
1261 1126559097 : ye = yd; p1 = *--xd;
1262 :
1263 1126559097 : *--zd = mulll(p1, *--yd); z2e = zd;
1264 9444033189 : while (yd > y) *--zd = addmul(p1, *--yd);
1265 1126559097 : *--zd = hiremainder;
1266 :
1267 10812596943 : while (xd > x)
1268 : {
1269 : LOCAL_OVERFLOW;
1270 9686037846 : yd = ye; p1 = *--xd;
1271 :
1272 9686037846 : z2d = --z2e;
1273 9686037846 : *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
1274 >12302*10^7 : while (yd > y)
1275 : {
1276 >11333*10^7 : hiremainder += overflow;
1277 >11333*10^7 : *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
1278 : }
1279 9686037846 : *--zd = hiremainder + overflow;
1280 : }
1281 1126559097 : if (*zd == 0) { zd++; lz--; } /* normalize */
1282 1126559097 : *--zd = evalsigne(1) | evallgefint(lz);
1283 1126559097 : *--zd = evaltyp(t_INT) | evallg(lz);
1284 1126559097 : return gc_const((pari_sp)zd, zd);
1285 : }
1286 :
1287 : INLINE GEN
1288 926242590 : sqrispec_basecase(GEN x, long nx)
1289 : {
1290 : GEN z2e,z2d,yd,xd,zd,x0,z0;
1291 : long p1,lz;
1292 : LOCAL_HIREMAINDER;
1293 : LOCAL_OVERFLOW;
1294 :
1295 926242590 : if (nx == 1) return sqru((ulong)*x);
1296 628644021 : if (nx == 0) return gen_0;
1297 225939060 : zd = (GEN)avma; lz = (nx+1) << 1;
1298 225939060 : z0 = new_chunk(lz);
1299 225939060 : if (nx == 1)
1300 : {
1301 0 : *--zd = mulll(*x, *x);
1302 0 : *--zd = hiremainder; goto END;
1303 : }
1304 225939060 : xd = x + nx;
1305 :
1306 : /* compute double products --> zd */
1307 225939060 : p1 = *--xd; yd = xd; --zd;
1308 225939060 : *--zd = mulll(p1, *--yd); z2e = zd;
1309 1374545859 : while (yd > x) *--zd = addmul(p1, *--yd);
1310 225939060 : *--zd = hiremainder;
1311 :
1312 225939060 : x0 = x+1;
1313 1374545859 : while (xd > x0)
1314 : {
1315 : LOCAL_OVERFLOW;
1316 1148606799 : p1 = *--xd; yd = xd;
1317 :
1318 1148606799 : z2e -= 2; z2d = z2e;
1319 1148606799 : *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
1320 10619931534 : while (yd > x)
1321 : {
1322 9471324735 : hiremainder += overflow;
1323 9471324735 : *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
1324 : }
1325 1148606799 : *--zd = hiremainder + overflow;
1326 : }
1327 : /* multiply zd by 2 (put result in zd - 1) */
1328 225939060 : zd[-1] = ((*zd & HIGHBIT) != 0);
1329 225939060 : shift_left(zd, zd, 0, (nx<<1)-3, 0, 1);
1330 :
1331 : /* add the squares */
1332 225939060 : xd = x + nx; zd = z0 + lz;
1333 225939060 : p1 = *--xd;
1334 225939060 : zd--; *zd = mulll(p1,p1);
1335 225939060 : zd--; *zd = addll(hiremainder, *zd);
1336 1600484919 : while (xd > x)
1337 : {
1338 1374545859 : p1 = *--xd;
1339 1374545859 : zd--; *zd = addll(mulll(p1,p1)+ overflow, *zd);
1340 1374545859 : zd--; *zd = addll(hiremainder + overflow, *zd);
1341 : }
1342 :
1343 225939060 : END:
1344 225939060 : if (*zd == 0) { zd++; lz--; } /* normalize */
1345 225939060 : *--zd = evalsigne(1) | evallgefint(lz);
1346 225939060 : *--zd = evaltyp(t_INT) | evallg(lz);
1347 225939060 : return gc_const((pari_sp)zd, zd);
1348 : }
1349 :
1350 : /********************************************************************/
1351 : /** **/
1352 : /** INTEGER MULTIPLICATION (FFT) **/
1353 : /** **/
1354 : /********************************************************************/
1355 :
1356 : /*
1357 : Compute parameters for FFT:
1358 : len: result length
1359 : k: FFT depth.
1360 : n: number of blocks (2^k)
1361 : bs: block size
1362 : mod: Modulus is M=2^(BIL*mod)+1
1363 : ord: order of 2 in Z/MZ.
1364 : We must have:
1365 : bs*n >= l
1366 : 2^(BIL*mod) > nb*2^(2*BIL*bs)
1367 : 2^k | 2*BIL*mod
1368 : */
1369 : static void
1370 85110 : mulliifft_params(long len, long *k, long *mod, long *bs, long *n, ulong *ord)
1371 : {
1372 : long r;
1373 85110 : *k = expu((3*len)>>2)-3;
1374 : do {
1375 85113 : (*k)--;
1376 85113 : r = *k-(TWOPOTBITS_IN_LONG+2);
1377 85113 : *n = 1L<<*k;
1378 85113 : *bs = (len+*n-1)>>*k;
1379 85113 : *mod= 2**bs+1;
1380 85113 : if (r>0)
1381 5154 : *mod=((*mod+(1L<<r)-1)>>r)<<r;
1382 85113 : } while(*mod>=3**bs);
1383 85110 : *ord= 4**mod*BITS_IN_LONG;
1384 85110 : }
1385 :
1386 : /* Zf_: arithmetic in ring Z/MZ where M= 2^(BITS_IN_LONG*mod)+1
1387 : * for some mod.
1388 : * Do not garbage collect.
1389 : */
1390 :
1391 : static GEN
1392 187122816 : Zf_add(GEN a, GEN b, GEN M)
1393 : {
1394 187122816 : GEN y, z = addii(a,b);
1395 187122816 : long mod = lgefint(M)-3;
1396 187122816 : long l = NLIMBS(z);
1397 187122816 : if (l<=mod) return z;
1398 72497628 : y = subiu(z, 1);
1399 72497628 : if (NLIMBS(y)<=mod) return z;
1400 72497628 : return int_normalize(y,1);
1401 : }
1402 :
1403 : static GEN
1404 190514628 : Zf_sub(GEN a, GEN b, GEN M)
1405 : {
1406 190514628 : GEN z = subii(a,b);
1407 190514628 : return signe(z)>=0? z: addii(M,z);
1408 : }
1409 :
1410 : /* destroy z */
1411 : static GEN
1412 398080593 : Zf_red_destroy(GEN z, GEN M)
1413 : {
1414 398080593 : long mod = lgefint(M)-3;
1415 398080593 : long l = NLIMBS(z);
1416 : GEN y;
1417 398080593 : if (l<=mod) return z;
1418 176688981 : y = shifti(z, -mod*BITS_IN_LONG);
1419 176688981 : z = int_normalize(z, NLIMBS(y));
1420 176688981 : y = Zf_red_destroy(y, M);
1421 176688981 : z = subii(z, y);
1422 176688981 : if (signe(z)<0) z = addii(z, M);
1423 176688981 : return z;
1424 : }
1425 :
1426 : INLINE GEN
1427 205647228 : Zf_shift(GEN a, ulong s, GEN M) { return Zf_red_destroy(shifti(a, s), M); }
1428 :
1429 : /*
1430 : Multiply by sqrt(2)^s
1431 : We use the formula sqrt(2)=z_8*(1-z_4)) && z_8=2^(ord/16) [2^(ord/4)+1]
1432 : */
1433 :
1434 : static GEN
1435 187122816 : Zf_mulsqrt2(GEN a, ulong s, ulong ord, GEN M)
1436 : {
1437 187122816 : ulong hord = ord>>1;
1438 187122816 : if (!signe(a)) return gen_0;
1439 183119220 : if (odd(s)) /* Multiply by 2^(s/2) */
1440 : {
1441 3391812 : GEN az8 = Zf_shift(a, ord>>4, M);
1442 3391812 : GEN az83 = Zf_shift(az8, ord>>3, M);
1443 3391812 : a = Zf_sub(az8, az83, M);
1444 3391812 : s--;
1445 : }
1446 183119220 : if (s < hord)
1447 136033635 : return Zf_shift(a, s>>1, M);
1448 : else
1449 47085585 : return subii(M,Zf_shift(a, (s-hord)>>1, M));
1450 : }
1451 :
1452 : INLINE GEN
1453 448896 : Zf_sqr(GEN a, GEN M) { return Zf_red_destroy(sqri(a), M); }
1454 :
1455 : INLINE GEN
1456 15295488 : Zf_mul(GEN a, GEN b, GEN M) { return Zf_red_destroy(mulii(a,b), M); }
1457 :
1458 : /* In place, bit reversing FFT */
1459 : static void
1460 30871311 : muliifft_dit(ulong o, ulong ord, GEN M, GEN FFT, long d, long step)
1461 : {
1462 30871311 : pari_sp av = avma;
1463 : long i;
1464 30871311 : ulong j, no = (o<<1)%ord;
1465 30871311 : long hstep=step>>1;
1466 154961487 : for (i = d+1, j = 0; i <= d+hstep; ++i, j =(j+o)%ord)
1467 : {
1468 124090176 : GEN a = Zf_add(gel(FFT,i), gel(FFT,i+hstep), M);
1469 124090176 : GEN b = Zf_mulsqrt2(Zf_sub(gel(FFT,i), gel(FFT,i+hstep), M), j, ord, M);
1470 124090176 : affii(a,gel(FFT,i));
1471 124090176 : affii(b,gel(FFT,i+hstep));
1472 124090176 : set_avma(av);
1473 : }
1474 30871311 : if (hstep>1)
1475 : {
1476 15351375 : muliifft_dit(no, ord, M, FFT, d, hstep);
1477 15351375 : muliifft_dit(no, ord, M, FFT, d+hstep, hstep);
1478 : }
1479 30871311 : }
1480 :
1481 : /* In place, bit reversed FFT, inverse of muliifft_dit */
1482 : static void
1483 15659274 : muliifft_dis(ulong o, ulong ord, GEN M, GEN FFT, long d, long step)
1484 : {
1485 15659274 : pari_sp av = avma;
1486 : long i;
1487 15659274 : ulong j, no = (o<<1)%ord;
1488 15659274 : long hstep=step>>1;
1489 15659274 : if (hstep>1)
1490 : {
1491 7787082 : muliifft_dis(no, ord, M, FFT, d, hstep);
1492 7787082 : muliifft_dis(no, ord, M, FFT, d+hstep, hstep);
1493 : }
1494 78691914 : for (i = d+1, j = 0; i <= d+hstep; ++i, j =(j+o)%ord)
1495 : {
1496 63032640 : GEN z = Zf_mulsqrt2(gel(FFT,i+hstep), j, ord, M);
1497 63032640 : GEN a = Zf_add(gel(FFT,i), z, M);
1498 63032640 : GEN b = Zf_sub(gel(FFT,i), z, M);
1499 63032640 : affii(a,gel(FFT,i));
1500 63032640 : affii(b,gel(FFT,i+hstep));
1501 63032640 : set_avma(av);
1502 : }
1503 15659274 : }
1504 :
1505 : static GEN
1506 168561 : muliifft_spliti(GEN a, long na, long bs, long n, long mod)
1507 : {
1508 168561 : GEN ap = a+na-1;
1509 168561 : GEN c = cgetg(n+1, t_VEC);
1510 : long i,j;
1511 31208433 : for(i=1;i<=n;i++)
1512 : {
1513 31039872 : GEN z = cgeti(mod+3);
1514 31039872 : if (na)
1515 : {
1516 15284268 : long m = minss(bs, na), v=0;
1517 15284268 : GEN zp, aa=ap-m+1;
1518 83156736 : while (!*aa && v<m) {aa++; v++;}
1519 15284268 : zp = z+m-v+1;
1520 380476317 : for (j=v; j < m; j++)
1521 365192049 : *zp-- = *ap--;
1522 15284268 : ap -= v; na -= m;
1523 15284268 : z[1] = evalsigne(m!=v) | evallgefint(m-v+2);
1524 : } else
1525 15755604 : z[1] = evalsigne(0) | evallgefint(2);
1526 31039872 : gel(c, i) = z;
1527 : }
1528 168561 : return c;
1529 : }
1530 :
1531 : static GEN
1532 85110 : muliifft_unspliti(GEN V, long bs, long len)
1533 : {
1534 85110 : long s, i, j, l = lg(V);
1535 85110 : GEN a = cgeti(len);
1536 85110 : a[1] = evalsigne(1)|evallgefint(len);
1537 439948050 : for(i=2;i<len;i++)
1538 439862940 : a[i] = 0;
1539 15829494 : for(i=1, s=0; i<l; i++, s+=bs)
1540 : {
1541 15744384 : GEN u = gel(V,i);
1542 15744384 : if (signe(u))
1543 : {
1544 15170145 : GEN ap = int_W(a,s);
1545 15170145 : GEN up = int_LSW(u);
1546 15170145 : long lu = NLIMBS(u);
1547 : LOCAL_OVERFLOW;
1548 15170145 : *ap = addll(*ap, *up--); ap--;
1549 860205717 : for (j=1; j<lu; j++)
1550 845035572 : { *ap = addllx(*ap, *up--); ap--; }
1551 15172701 : while (overflow)
1552 2556 : { *ap = addll(*ap, 1); ap--; }
1553 : }
1554 : }
1555 85110 : return int_normalize(a,0);
1556 : }
1557 :
1558 : static GEN
1559 1659 : sqrispec_fft(GEN a, long na)
1560 : {
1561 1659 : pari_sp av, ltop = avma;
1562 1659 : long len = 2*na;
1563 : long k, mod, bs, n;
1564 : GEN FFT, M;
1565 : long i;
1566 : ulong o, ord;
1567 :
1568 1659 : mulliifft_params(len,&k,&mod,&bs,&n,&ord);
1569 1659 : o = ord>>k;
1570 1659 : M = int2n(mod*BITS_IN_LONG);
1571 1659 : M[2+mod] = 1;
1572 1659 : FFT = muliifft_spliti(a, na, bs, n, mod);
1573 1659 : muliifft_dit(o, ord, M, FFT, 0, n);
1574 1659 : av = avma;
1575 450555 : for(i=1; i<=n; i++)
1576 : {
1577 448896 : affii(Zf_sqr(gel(FFT,i), M), gel(FFT,i));
1578 448896 : set_avma(av);
1579 : }
1580 1659 : muliifft_dis(ord-o, ord, M, FFT, 0, n);
1581 450555 : for(i=1; i<=n; i++)
1582 : {
1583 448896 : affii(Zf_shift(gel(FFT,i), (ord>>1)-k, M), gel(FFT,i));
1584 448896 : set_avma(av);
1585 : }
1586 1659 : return gerepileuptoint(ltop, muliifft_unspliti(FFT,bs,2+len));
1587 : }
1588 :
1589 : static GEN
1590 83451 : muliispec_fft(GEN a, GEN b, long na, long nb)
1591 : {
1592 83451 : pari_sp av, av2, ltop = avma;
1593 83451 : long len = na+nb;
1594 : long k, mod, bs, n;
1595 : GEN FFT, FFTb, M;
1596 : long i;
1597 : ulong o, ord;
1598 :
1599 83451 : mulliifft_params(len,&k,&mod,&bs,&n,&ord);
1600 83451 : o = ord>>k;
1601 83451 : M = int2n(mod*BITS_IN_LONG);
1602 83451 : M[2+mod] = 1;
1603 83451 : FFT = muliifft_spliti(a, na, bs, n, mod);
1604 83451 : av=avma;
1605 83451 : muliifft_dit(o, ord, M, FFT, 0, n);
1606 83451 : FFTb = muliifft_spliti(b, nb, bs, n, mod);
1607 83451 : av2 = avma;
1608 83451 : muliifft_dit(o, ord, M, FFTb, 0, n);
1609 15378939 : for(i=1; i<=n; i++)
1610 : {
1611 15295488 : affii(Zf_mul(gel(FFT,i), gel(FFTb,i), M), gel(FFT,i));
1612 15295488 : set_avma(av2);
1613 : }
1614 83451 : set_avma(av);
1615 83451 : muliifft_dis(ord-o, ord, M, FFT, 0, n);
1616 15378939 : for(i=1; i<=n; i++)
1617 : {
1618 15295488 : affii(Zf_shift(gel(FFT,i),(ord>>1)-k,M), gel(FFT,i));
1619 15295488 : set_avma(av);
1620 : }
1621 83451 : return gerepileuptoint(ltop, muliifft_unspliti(FFT,bs,2+len));
1622 : }
1623 :
1624 : /********************************************************************/
1625 : /** **/
1626 : /** INTEGER MULTIPLICATION (KARATSUBA) **/
1627 : /** **/
1628 : /********************************************************************/
1629 :
1630 : /* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */
1631 : static GEN
1632 703499397 : addshiftw(GEN x, GEN y, long d)
1633 : {
1634 703499397 : GEN z,z0,y0,yd, zd = (GEN)avma;
1635 703499397 : long a,lz,ly = lgefint(y);
1636 :
1637 703499397 : z0 = new_chunk(d);
1638 703499397 : a = ly-2; yd = y+ly;
1639 703499397 : if (a >= d)
1640 : {
1641 12691640283 : y0 = yd-d; while (yd > y0) *--zd = *--yd; /* copy last d words of y */
1642 698904828 : a -= d;
1643 698904828 : if (a)
1644 466572432 : z = addiispec(LIMBS(x), LIMBS(y), NLIMBS(x), a);
1645 : else
1646 232332396 : z = icopy(x);
1647 : }
1648 : else
1649 : {
1650 17039847 : y0 = yd-a; while (yd > y0) *--zd = *--yd; /* copy last a words of y */
1651 69867249 : while (zd > z0) *--zd = 0; /* complete with 0s */
1652 4594569 : z = icopy(x);
1653 : }
1654 703499397 : lz = lgefint(z)+d;
1655 703499397 : z[1] = evalsigne(1) | evallgefint(lz);
1656 703499397 : z[0] = evaltyp(t_INT) | evallg(lz); return z;
1657 : }
1658 :
1659 : /* Fast product (Karatsuba) of integers. a and b are "special" GENs
1660 : * c,c0,c1,c2 are genuine GENs.
1661 : */
1662 : GEN
1663 5385087453 : muliispec(GEN a, GEN b, long na, long nb)
1664 : {
1665 : GEN a0,c,c0;
1666 : long n0, n0a, i;
1667 : pari_sp av;
1668 :
1669 5385087453 : if (na < nb) swapspec(a,b, na,nb);
1670 5385087453 : if (nb < MULII_KARATSUBA_LIMIT) return muliispec_basecase(a,b,na,nb);
1671 226566660 : if (nb >= MULII_FFT_LIMIT) return muliispec_fft(a,b,na,nb);
1672 226483209 : i=(na>>1); n0=na-i; na=i;
1673 226483209 : av=avma; a0=a+na; n0a=n0;
1674 340385085 : while (n0a && !*a0) { a0++; n0a--; }
1675 :
1676 226483209 : if (n0a && nb > n0)
1677 222993813 : { /* nb <= na <= n0 */
1678 : GEN b0,c1,c2;
1679 : long n0b;
1680 :
1681 222993813 : nb -= n0;
1682 222993813 : c = muliispec(a,b,na,nb);
1683 222993813 : b0 = b+nb; n0b = n0;
1684 321657696 : while (n0b && !*b0) { b0++; n0b--; }
1685 222993813 : if (n0b)
1686 : {
1687 222231015 : c0 = muliispec(a0,b0, n0a,n0b);
1688 :
1689 222231015 : c2 = addiispec(a0,a, n0a,na);
1690 222231015 : c1 = addiispec(b0,b, n0b,nb);
1691 222231015 : c1 = muliispec(LIMBS(c1),LIMBS(c2), NLIMBS(c1),NLIMBS(c2));
1692 222231015 : c2 = addiispec(LIMBS(c0),LIMBS(c), NLIMBS(c0),NLIMBS(c));
1693 :
1694 222231015 : c1 = subiispec(LIMBS(c1),LIMBS(c2), NLIMBS(c1),NLIMBS(c2));
1695 : }
1696 : else
1697 : {
1698 762798 : c0 = gen_0;
1699 762798 : c1 = muliispec(a0,b, n0a,nb);
1700 : }
1701 222993813 : c = addshiftw(c,c1, n0);
1702 : }
1703 : else
1704 : {
1705 3489396 : c = muliispec(a,b,na,nb);
1706 3489396 : c0 = muliispec(a0,b,n0a,nb);
1707 : }
1708 226483209 : return gerepileuptoint(av, addshiftw(c,c0, n0));
1709 : }
1710 : GEN
1711 165798 : muluui(ulong x, ulong y, GEN z)
1712 : {
1713 165798 : long t, s = signe(z);
1714 : GEN r;
1715 : LOCAL_HIREMAINDER;
1716 :
1717 165798 : if (!x || !y || !signe(z)) return gen_0;
1718 165501 : t = mulll(x,y);
1719 165501 : if (!hiremainder)
1720 165501 : r = muluispec(t, z+2, lgefint(z)-2);
1721 : else
1722 : {
1723 : long tmp[2];
1724 0 : tmp[0] = hiremainder;
1725 0 : tmp[1] = t;
1726 0 : r = muliispec(z+2,tmp,lgefint(z)-2,2);
1727 : }
1728 165501 : setsigne(r,s); return r;
1729 : }
1730 :
1731 : #define sqrispec_mirror sqrispec
1732 : #define muliispec_mirror muliispec
1733 :
1734 : /* x % (2^n), assuming n >= 0 */
1735 : GEN
1736 18491472 : remi2n(GEN x, long n)
1737 : {
1738 18491472 : long hi,l,k,lx,ly, sx = signe(x);
1739 : GEN z, xd, zd;
1740 :
1741 18491472 : if (!sx || !n) return gen_0;
1742 :
1743 18468837 : k = dvmdsBIL(n, &l);
1744 18468837 : lx = lgefint(x);
1745 18468837 : if (lx < k+3) return icopy(x);
1746 :
1747 18081939 : xd = x + (lx-k-1);
1748 : /* x = |_|...|#|1|...|k| : copy the last l bits of # and the last k words
1749 : * ^--- initial xd */
1750 18081939 : hi = ((ulong)*xd) & ((1UL<<l)-1); /* last l bits of # = top bits of result */
1751 18081939 : if (!hi)
1752 : { /* strip leading zeroes from result */
1753 827559 : xd++; while (k && !*xd) { k--; xd++; }
1754 800646 : if (!k) return gen_0;
1755 602091 : ly = k+2; xd--;
1756 : }
1757 : else
1758 17281293 : ly = k+3;
1759 :
1760 17883384 : zd = z = cgeti(ly);
1761 17883384 : *++zd = evalsigne(sx) | evallgefint(ly);
1762 17883384 : if (hi) *++zd = hi;
1763 102077859 : for ( ;k; k--) *++zd = *++xd;
1764 17883384 : return z;
1765 : }
1766 :
1767 : GEN
1768 936180831 : sqrispec(GEN a, long na)
1769 : {
1770 : GEN a0,c;
1771 : long n0, n0a, i;
1772 : pari_sp av;
1773 :
1774 936180831 : if (na < SQRI_KARATSUBA_LIMIT) return sqrispec_basecase(a,na);
1775 9938241 : if (na >= SQRI_FFT_LIMIT) return sqrispec_fft(a,na);
1776 9936582 : i=(na>>1); n0=na-i; na=i;
1777 9936582 : av=avma; a0=a+na; n0a=n0;
1778 14932635 : while (n0a && !*a0) { a0++; n0a--; }
1779 9936582 : c = sqrispec(a,na);
1780 9936582 : if (n0a)
1781 : {
1782 9927684 : GEN t, c1, c0 = sqrispec(a0,n0a);
1783 : #if 0
1784 : c1 = shifti(muliispec(a0,a, n0a,na),1);
1785 : #else /* faster */
1786 9927684 : t = addiispec(a0,a,n0a,na);
1787 9927684 : t = sqrispec(LIMBS(t),NLIMBS(t));
1788 9927684 : c1= addiispec(LIMBS(c0),LIMBS(c), NLIMBS(c0), NLIMBS(c));
1789 9927684 : c1= subiispec(LIMBS(t),LIMBS(c1), NLIMBS(t), NLIMBS(c1));
1790 : #endif
1791 9927684 : c = addshiftw(c,c1, n0);
1792 9927684 : c = addshiftw(c,c0, n0);
1793 : }
1794 : else
1795 8898 : c = addshiftw(c,gen_0,n0<<1);
1796 9936582 : return gerepileuptoint(av, c);
1797 : }
1798 :
1799 : /********************************************************************/
1800 : /** **/
1801 : /** KARATSUBA SQUARE ROOT **/
1802 : /** adapted from Paul Zimmermann's implementation of **/
1803 : /** his algorithm in GMP (mpn_sqrtrem) **/
1804 : /** **/
1805 : /********************************************************************/
1806 :
1807 : /* Square roots table */
1808 : static const unsigned char approx_tab[192] = {
1809 : 128,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,
1810 : 143,144,144,145,146,147,148,149,150,150,151,152,153,154,155,155,
1811 : 156,157,158,159,160,160,161,162,163,163,164,165,166,167,167,168,
1812 : 169,170,170,171,172,173,173,174,175,176,176,177,178,178,179,180,
1813 : 181,181,182,183,183,184,185,185,186,187,187,188,189,189,190,191,
1814 : 192,192,193,193,194,195,195,196,197,197,198,199,199,200,201,201,
1815 : 202,203,203,204,204,205,206,206,207,208,208,209,209,210,211,211,
1816 : 212,212,213,214,214,215,215,216,217,217,218,218,219,219,220,221,
1817 : 221,222,222,223,224,224,225,225,226,226,227,227,228,229,229,230,
1818 : 230,231,231,232,232,233,234,234,235,235,236,236,237,237,238,238,
1819 : 239,240,240,241,241,242,242,243,243,244,244,245,245,246,246,247,
1820 : 247,248,248,249,249,250,250,251,251,252,252,253,253,254,254,255
1821 : };
1822 :
1823 : /* N[0], assume N[0] >= 2^(BIL-2).
1824 : * Return r,s such that s^2 + r = N, 0 <= r <= 2s */
1825 : static void
1826 94680981 : p_sqrtu1(ulong *N, ulong *ps, ulong *pr)
1827 : {
1828 94680981 : ulong prec, r, s, q, u, n0 = N[0];
1829 :
1830 94680981 : q = n0 >> (BITS_IN_LONG - 8);
1831 : /* 2^6 = 64 <= q < 256 = 2^8 */
1832 94680981 : s = approx_tab[q - 64]; /* 128 <= s < 255 */
1833 94680981 : r = (n0 >> (BITS_IN_LONG - 16)) - s * s; /* r <= 2*s */
1834 94680981 : if (r > (s << 1)) { r -= (s << 1) | 1; s++; }
1835 :
1836 : /* 8-bit approximation from the high 8-bits of N[0] */
1837 94680981 : prec = 8;
1838 94680981 : n0 <<= 2 * prec;
1839 284042943 : while (2 * prec < BITS_IN_LONG)
1840 : { /* invariant: s has prec bits, and r <= 2*s */
1841 189361962 : r = (r << prec) + (n0 >> (BITS_IN_LONG - prec));
1842 189361962 : n0 <<= prec;
1843 189361962 : u = 2 * s;
1844 189361962 : q = r / u; u = r - q * u;
1845 189361962 : s = (s << prec) + q;
1846 189361962 : u = (u << prec) + (n0 >> (BITS_IN_LONG - prec));
1847 189361962 : q = q * q;
1848 189361962 : r = u - q;
1849 189361962 : if (u < q) { s--; r += (s << 1) | 1; }
1850 189361962 : n0 <<= prec;
1851 189361962 : prec = 2 * prec;
1852 : }
1853 94680981 : *ps = s;
1854 94680981 : *pr = r;
1855 94680981 : }
1856 :
1857 : /* N[0..1], assume N[0] >= 2^(BIL-2).
1858 : * Return 1 if remainder overflows, 0 otherwise */
1859 : static int
1860 91978353 : p_sqrtu2(ulong *N, ulong *ps, ulong *pr)
1861 : {
1862 91978353 : ulong cc, qhl, r, s, q, u, n1 = N[1];
1863 : LOCAL_OVERFLOW;
1864 :
1865 91978353 : p_sqrtu1(N, &s, &r); /* r <= 2s */
1866 137594091 : qhl = 0; while (r >= s) { qhl++; r -= s; }
1867 : /* now r < s < 2^(BIL/2) */
1868 91978353 : r = (r << BITS_IN_HALFULONG) | (n1 >> BITS_IN_HALFULONG);
1869 91978353 : u = s << 1;
1870 91978353 : q = r / u; u = r - q * u;
1871 91978353 : q += (qhl & 1) << (BITS_IN_HALFULONG - 1);
1872 91978353 : qhl >>= 1;
1873 : /* (initial r)<<(BIL/2) + n1>>(BIL/2) = (qhl<<(BIL/2) + q) * 2s + u */
1874 91978353 : s = ((s + qhl) << BITS_IN_HALFULONG) + q;
1875 91978353 : cc = u >> BITS_IN_HALFULONG;
1876 91978353 : r = (u << BITS_IN_HALFULONG) | (n1 & LOWMASK);
1877 91978353 : r = subll(r, q * q);
1878 91978353 : cc -= overflow + qhl;
1879 : /* now subtract 2*q*2^(BIL/2) + 2^BIL if qhl is set */
1880 91978353 : if ((long)cc < 0)
1881 : {
1882 23471328 : if (s) {
1883 23417721 : r = addll(r, s);
1884 23417721 : cc += overflow;
1885 23417721 : s--;
1886 : } else {
1887 53607 : cc++;
1888 53607 : s = ~0UL;
1889 : }
1890 23471328 : r = addll(r, s);
1891 23471328 : cc += overflow;
1892 : }
1893 91978353 : *ps = s;
1894 91978353 : *pr = r; return cc;
1895 : }
1896 :
1897 : static void
1898 90916779 : xmpn_zero(GEN x, long n)
1899 : {
1900 695092197 : while (--n >= 0) x[n]=0;
1901 90916779 : }
1902 : static void
1903 1068143847 : xmpn_copy(GEN z, GEN x, long n)
1904 : {
1905 1068143847 : long k = n;
1906 4197216156 : while (--k >= 0) z[k] = x[k];
1907 1068143847 : }
1908 : /* a[0..la-1] * 2^(lb BIL) | b[0..lb-1] */
1909 : static GEN
1910 468316218 : catii(GEN a, long la, GEN b, long lb)
1911 : {
1912 468316218 : long l = la + lb + 2;
1913 468316218 : GEN z = cgetipos(l);
1914 468316218 : xmpn_copy(LIMBS(z), a, la);
1915 468316218 : xmpn_copy(LIMBS(z) + la, b, lb);
1916 468316218 : return int_normalize(z, 0);
1917 : }
1918 :
1919 : /* sqrt n[0..1], assume n normalized */
1920 : static GEN
1921 91703244 : sqrtispec2(GEN n, GEN *pr)
1922 : {
1923 : ulong s, r;
1924 91703244 : int hi = p_sqrtu2((ulong*)n, &s, &r);
1925 91703244 : GEN S = utoi(s);
1926 91703244 : *pr = hi? uutoi(1,r): utoi(r);
1927 91703244 : return S;
1928 : }
1929 :
1930 : /* sqrt n[0], _dont_ assume n normalized */
1931 : static GEN
1932 2702628 : sqrtispec1_sh(GEN n, GEN *pr)
1933 : {
1934 : GEN S;
1935 2702628 : ulong r, s, u0 = uel(n,0);
1936 2702628 : int sh = bfffo(u0) & ~1UL;
1937 2702628 : if (sh) u0 <<= sh;
1938 2702628 : p_sqrtu1(&u0, &s, &r);
1939 : /* s^2 + r = u0, s < 2^(BIL/2). Rescale back:
1940 : * 2^(2k) n = S^2 + R
1941 : * so 2^(2k) n = (S - s0)^2 + (2*S*s0 - s0^2 + R), s0 = S mod 2^k. */
1942 2702628 : if (sh) {
1943 1629945 : int k = sh >> 1;
1944 1629945 : ulong s0 = s & ((1L<<k) - 1);
1945 1629945 : r += s * (s0<<1);
1946 1629945 : s >>= k;
1947 1629945 : r >>= sh;
1948 : }
1949 2702628 : S = utoi(s);
1950 2702628 : if (pr) *pr = utoi(r);
1951 2702628 : return S;
1952 : }
1953 :
1954 : /* sqrt n[0..1], _dont_ assume n normalized */
1955 : static GEN
1956 275109 : sqrtispec2_sh(GEN n, GEN *pr)
1957 : {
1958 : GEN S;
1959 275109 : ulong U[2], r, s, u0 = uel(n,0), u1 = uel(n,1);
1960 275109 : int hi, sh = bfffo(u0) & ~1UL;
1961 275109 : if (sh) {
1962 247146 : u0 = (u0 << sh) | (u1 >> (BITS_IN_LONG-sh));
1963 247146 : u1 <<= sh;
1964 : }
1965 275109 : U[0] = u0;
1966 275109 : U[1] = u1; hi = p_sqrtu2(U, &s, &r);
1967 : /* s^2 + R = u0|u1. Rescale back:
1968 : * 2^(2k) n = S^2 + R
1969 : * so 2^(2k) n = (S - s0)^2 + (2*S*s0 - s0^2 + R), s0 = S mod 2^k. */
1970 275109 : if (sh) {
1971 247146 : int k = sh >> 1;
1972 247146 : ulong s0 = s & ((1L<<k) - 1);
1973 : LOCAL_HIREMAINDER;
1974 : LOCAL_OVERFLOW;
1975 247146 : r = addll(r, mulll(s, (s0<<1)));
1976 247146 : if (overflow) hiremainder++;
1977 247146 : hiremainder += hi; /* + 0 or 1 */
1978 247146 : s >>= k;
1979 247146 : r = (r>>sh) | (hiremainder << (BITS_IN_LONG-sh));
1980 247146 : hi = (hiremainder & (1L<<sh));
1981 : }
1982 275109 : S = utoi(s);
1983 275109 : if (pr) *pr = hi? uutoi(1,r): utoi(r);
1984 275109 : return S;
1985 : }
1986 :
1987 : /* Let N = N[0..2n-1]. Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S
1988 : * Assume N normalized */
1989 : static GEN
1990 325861353 : sqrtispec(GEN N, long n, GEN *r)
1991 : {
1992 : GEN S, R, q, z, u;
1993 : long l, h;
1994 :
1995 325861353 : if (n == 1) return sqrtispec2(N, r);
1996 234158109 : l = n >> 1;
1997 234158109 : h = n - l; /* N = a3(h) | a2(h) | a1(l) | a0(l words) */
1998 234158109 : S = sqrtispec(N, h, &R); /* S^2 + R = a3|a2 */
1999 :
2000 234158109 : z = catii(LIMBS(R), NLIMBS(R), N + 2*h, l); /* = R | a1(l) */
2001 234158109 : q = dvmdii(z, shifti(S,1), &u);
2002 234158109 : z = catii(LIMBS(u), NLIMBS(u), N + n + h, l); /* = u | a0(l) */
2003 :
2004 234158109 : S = addshiftw(S, q, l);
2005 234158109 : R = subii(z, sqri(q));
2006 234158109 : if (signe(R) < 0)
2007 : {
2008 40019820 : GEN S2 = shifti(S,1);
2009 40019820 : R = addis(subiispec(LIMBS(S2),LIMBS(R), NLIMBS(S2),NLIMBS(R)), -1);
2010 40019820 : S = addis(S, -1);
2011 : }
2012 234158109 : *r = R; return S;
2013 : }
2014 :
2015 : /* Return S (and set R) s.t S^2 + R = N, 0 <= R <= 2S.
2016 : * As for dvmdii, R is last on stack and guaranteed to be gen_0 in case the
2017 : * remainder is 0. R = NULL is allowed. */
2018 : GEN
2019 3764757 : sqrtremi(GEN N, GEN *r)
2020 : {
2021 : pari_sp av;
2022 3764757 : GEN S, R, n = N+2;
2023 3764757 : long k, l2, ln = NLIMBS(N);
2024 : int sh;
2025 :
2026 3764757 : if (ln <= 2)
2027 : {
2028 2978292 : if (ln == 2) return sqrtispec2_sh(n, r);
2029 2703183 : if (ln == 1) return sqrtispec1_sh(n, r);
2030 555 : if (r) *r = gen_0;
2031 555 : return gen_0;
2032 : }
2033 786465 : av = avma;
2034 786465 : sh = bfffo(n[0]) >> 1;
2035 786465 : l2 = (ln + 1) >> 1;
2036 786465 : if (sh || (ln & 1)) { /* normalize n, so that n[0] >= 2^BIL / 4 */
2037 785781 : GEN s0, t = new_chunk(ln + 1);
2038 785781 : t[ln] = 0;
2039 785781 : if (sh)
2040 783936 : shift_left(t, n, 0,ln-1, 0, sh << 1);
2041 : else
2042 1845 : xmpn_copy(t, n, ln);
2043 785781 : S = sqrtispec(t, l2, &R); /* t normalized, 2 * l2 words */
2044 : /* Rescale back:
2045 : * 2^(2k) n = S^2 + R, k = sh + (ln & 1)*BIL/2
2046 : * so 2^(2k) n = (S - s0)^2 + (2*S*s0 - s0^2 + R), s0 = S mod 2^k. */
2047 785781 : k = sh + (ln & 1) * (BITS_IN_LONG/2);
2048 785781 : s0 = remi2n(S, k);
2049 785781 : R = addii(shifti(R,-1), mulii(s0, S));
2050 785781 : R = shifti(R, 1 - (k<<1));
2051 785781 : S = shifti(S, -k);
2052 : }
2053 : else
2054 684 : S = sqrtispec(n, l2, &R);
2055 :
2056 786465 : if (!r) { set_avma((pari_sp)S); return gerepileuptoint(av, S); }
2057 722778 : *r = R; return gc_all(av, 2, &S, r);
2058 : }
2059 :
2060 : /* compute sqrt(|a|), assuming a != 0 */
2061 :
2062 : #if 1
2063 : GEN
2064 90916779 : sqrtr_abs(GEN x)
2065 : {
2066 90916779 : long l = lg(x) - 2, e = expo(x), er = e>>1;
2067 90916779 : GEN b, c, res = cgetg(2 + l, t_REAL);
2068 90916779 : res[1] = evalsigne(1) | evalexpo(er);
2069 90916779 : if (e&1) {
2070 40592787 : b = new_chunk(l << 1);
2071 40592787 : xmpn_copy(b, x+2, l);
2072 40592787 : xmpn_zero(b + l,l);
2073 40592787 : b = sqrtispec(b, l, &c);
2074 40592787 : xmpn_copy(res+2, b+2, l);
2075 40592787 : if (cmpii(c, b) > 0) roundr_up_ip(res, l+2);
2076 : } else {
2077 : ulong u;
2078 50323992 : b = new_chunk(2 + (l << 1));
2079 50323992 : shift_left(b+1, x+2, 0,l-1, 0, BITS_IN_LONG-1);
2080 50323992 : b[0] = uel(x,2)>>1;
2081 50323992 : xmpn_zero(b + l+1,l+1);
2082 50323992 : b = sqrtispec(b, l+1, &c);
2083 50323992 : xmpn_copy(res+2, b+2, l);
2084 50323992 : u = uel(b,l+2);
2085 50323992 : if ( u&HIGHBIT || (u == ~HIGHBIT && cmpii(c,b) > 0))
2086 24823770 : roundr_up_ip(res, l+2);
2087 : }
2088 90916779 : return gc_const((pari_sp)res, res);
2089 : }
2090 :
2091 : #else /* use t_REAL: currently much slower (quadratic division) */
2092 :
2093 : #ifdef LONG_IS_64BIT
2094 : /* 64 bits of b = sqrt(a[0] * 2^64 + a[1]) [ up to 1ulp ] */
2095 : static ulong
2096 : sqrtu2(ulong *a)
2097 : {
2098 : ulong c, b = dblmantissa( sqrt((double)a[0]) );
2099 : LOCAL_HIREMAINDER;
2100 : LOCAL_OVERFLOW;
2101 :
2102 : /* > 32 correct bits, 1 Newton iteration to reach 64 */
2103 : if (b <= a[0]) return HIGHBIT | (a[0] >> 1);
2104 : hiremainder = a[0]; c = divll(a[1], b);
2105 : return (addll(c, b) >> 1) | HIGHBIT;
2106 : }
2107 : /* 64 bits of sqrt(a[0] * 2^63) */
2108 : static ulong
2109 : sqrtu2_1(ulong *a)
2110 : {
2111 : ulong t[2];
2112 : t[0] = (a[0] >> 1);
2113 : t[1] = (a[0] << (BITS_IN_LONG-1)) | (a[1] >> 1);
2114 : return sqrtu2(t);
2115 : }
2116 : #else
2117 : /* 32 bits of sqrt(a[0] * 2^32) */
2118 : static ulong
2119 : sqrtu2(ulong *a) { return dblmantissa( sqrt((double)a[0]) ); }
2120 : /* 32 bits of sqrt(a[0] * 2^31) */
2121 : static ulong
2122 : sqrtu2_1(ulong *a) { return dblmantissa( sqrt(2. * a[0]) ); }
2123 : #endif
2124 :
2125 : GEN
2126 : sqrtr_abs(GEN x)
2127 : {
2128 : long l1, i, l = lg(x), ex = expo(x);
2129 : GEN a, t, y = cgetg(l, t_REAL);
2130 : pari_sp av, av0 = avma;
2131 :
2132 : a = rtor(x, lg2prec(l+1));
2133 : t = cgetg(l+1, t_REAL);
2134 : if (ex & 1) { /* odd exponent */
2135 : a[1] = evalsigne(1) | _evalexpo(1);
2136 : t[2] = (long)sqrtu2((ulong*)a + 2);
2137 : } else { /* even exponent */
2138 : a[1] = evalsigne(1) | _evalexpo(0);
2139 : t[2] = (long)sqrtu2_1((ulong*)a + 2);
2140 : }
2141 : t[1] = evalsigne(1) | _evalexpo(0);
2142 : for (i = 3; i <= l; i++) t[i] = 0;
2143 :
2144 : /* |x| = 2^(ex/2) a, t ~ sqrt(a) */
2145 : l--; l1 = 1; av = avma;
2146 : while (l1 < l) { /* let t := (t + a/t)/2 */
2147 : l1 <<= 1; if (l1 > l) l1 = l;
2148 : setlg(a, l1 + 2);
2149 : setlg(t, l1 + 2);
2150 : affrr(addrr(t, divrr(a,t)), t); shiftr_inplace(t, -1);
2151 : set_avma(av);
2152 : }
2153 : affrr(t,y); shiftr_inplace(y, (ex>>1));
2154 : return gc_const(av0, y);
2155 : }
2156 :
2157 : #endif
2158 :
2159 : /*******************************************************************
2160 : * *
2161 : * Base Conversion *
2162 : * *
2163 : *******************************************************************/
2164 :
2165 : static void
2166 733815 : convi_dac(GEN x, ulong l, ulong *res)
2167 : {
2168 733815 : pari_sp ltop=avma;
2169 : ulong m;
2170 : GEN x1,x2;
2171 733815 : if (l==1) { *res=itou(x); return; }
2172 348042 : m=l>>1;
2173 348042 : x1=dvmdii(x,powuu(1000000000UL,m),&x2);
2174 348042 : convi_dac(x1,l-m,res+m);
2175 348042 : convi_dac(x2,m,res);
2176 348042 : set_avma(ltop);
2177 : }
2178 :
2179 : /* convert integer --> base 10^9 [not memory clean] */
2180 : ulong *
2181 316222 : convi(GEN x, long *l)
2182 : {
2183 316222 : long lz, lx = lgefint(x);
2184 : ulong *z;
2185 316222 : if (lx == 3 && uel(x,2) < 1000000000UL) {
2186 278491 : z = (ulong*)new_chunk(1);
2187 278491 : *z = x[2];
2188 278491 : *l = 1; return z+1;
2189 : }
2190 37731 : lz = 1 + (long)bit_accuracy_mul(lx, LOG10_2/9);
2191 37731 : z = (ulong*)new_chunk(lz);
2192 37731 : convi_dac(x,(ulong)lz,z);
2193 67755 : while (z[lz-1]==0) lz--;
2194 37731 : *l=lz; return z+lz;
2195 : }
2196 :
|
