1 #line 2 "../src/kernel/none/level1.h"
2 /* Copyright (C) 2000  The PARI group.
3 
4 This file is part of the PARI/GP package.
5 
6 PARI/GP is free software; you can redistribute it and/or modify it under the
7 terms of the GNU General Public License as published by the Free Software
8 Foundation; either version 2 of the License, or (at your option) any later
9 version. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11 
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
15 
16 /* This file defines "level 1" kernel functions.
17  * These functions can be inline; they are also defined externally in
18  * mpinl.c, which includes this file and never needs to be changed */
19 
20 INLINE long
evallg(long x)21 evallg(long x)
22 {
23   if (x & ~LGBITS) pari_err_OVERFLOW("lg()");
24   return _evallg(x);
25 }
26 INLINE long
evalvalp(long x)27 evalvalp(long x)
28 {
29   long v = _evalvalp(x);
30   if (v & ~VALPBITS) pari_err_OVERFLOW("valp()");
31   return v;
32 }
33 INLINE long
evalexpo(long x)34 evalexpo(long x)
35 {
36   long v = _evalexpo(x);
37   if (v & ~EXPOBITS) pari_err_OVERFLOW("expo()");
38   return v;
39 }
40 INLINE long
evalprecp(long x)41 evalprecp(long x)
42 {
43   long v = _evalprecp(x);
44   if (x & ~((1UL<<(BITS_IN_LONG-VALPnumBITS))-1)) pari_err_OVERFLOW("precp()");
45   return v;
46 }
47 
48 INLINE int
varncmp(long x,long y)49 varncmp(long x, long y)
50 {
51   if (varpriority[x] < varpriority[y]) return  1;
52   if (varpriority[x] > varpriority[y]) return -1;
53   return 0;
54 }
55 INLINE long
varnmin(long x,long y)56 varnmin(long x, long y)
57 { return (varpriority[x] <= varpriority[y])? x: y; }
58 INLINE long
varnmax(long x,long y)59 varnmax(long x, long y)
60 { return (varpriority[x] >= varpriority[y])? x: y; }
61 
62 /* Inhibit some area gerepile-wise: declare it to be a non recursive
63  * type, of length l. Thus gerepile won't inspect the zone, just copy it.
64  * For the following situation:
65  *   z = cgetg(t,a); av = avma; garbage(); ltop = avma;
66  *   for (i=1; i<HUGE; i++) gel(z,i) = blah();
67  *   stackdummy(av,ltop);
68  * loses (av-ltop) words but save a costly gerepile. */
69 INLINE void
stackdummy(pari_sp av,pari_sp ltop)70 stackdummy(pari_sp av, pari_sp ltop) {
71   long l = ((GEN)av) - ((GEN)ltop);
72   if (l > 0) {
73     GEN z = (GEN)ltop;
74     z[0] = evaltyp(t_VECSMALL) | evallg(l);
75 #ifdef DEBUG
76     { long i; for (i = 1; i < l; i++) z[i] = 0; }
77 #endif
78   }
79 }
80 INLINE void
fixlg(GEN x,long ly)81 fixlg(GEN x, long ly) {
82   long lx = lg(x), l = lx - ly;
83   if (l > 0)
84   { /* stackdummy(x+lx, x+ly) */
85     GEN z = x + ly;
86     z[0] = evaltyp(t_VECSMALL) | evallg(l);
87     setlg(x, ly);
88 #ifdef DEBUG
89     { long i; for (i = 1; i < l; i++) z[i] = 0; }
90 #endif
91   }
92 }
93 /* update lg(z) before affrr(y, z)  [ to cater for precision loss ]*/
94 INLINE void
affrr_fixlg(GEN y,GEN z)95 affrr_fixlg(GEN y, GEN z) { fixlg(z, lg(y)); affrr(y, z); }
96 
97 /*******************************************************************/
98 /*                                                                 */
99 /*                       ALLOCATE ON STACK                         */
100 /*                                                                 */
101 /*******************************************************************/
102 INLINE void
set_avma(ulong av)103 set_avma(ulong av) { avma = av; }
104 
105 INLINE double
gc_double(pari_sp av,double d)106 gc_double(pari_sp av, double d) { set_avma(av); return d; }
107 INLINE long
gc_long(pari_sp av,long s)108 gc_long(pari_sp av, long s) { set_avma(av); return s; }
109 INLINE ulong
gc_ulong(pari_sp av,ulong s)110 gc_ulong(pari_sp av, ulong s) { set_avma(av); return s; }
111 INLINE int
gc_bool(pari_sp av,int s)112 gc_bool(pari_sp av, int s) { set_avma(av); return s; }
113 INLINE int
gc_int(pari_sp av,int s)114 gc_int(pari_sp av, int s) { set_avma(av); return s; }
115 INLINE GEN
gc_NULL(pari_sp av)116 gc_NULL(pari_sp av) { set_avma(av); return NULL; }
117 INLINE GEN
gc_const(pari_sp av,GEN x)118 gc_const(pari_sp av, GEN x) { set_avma(av); return x; }
119 
120 INLINE GEN
new_chunk(size_t x)121 new_chunk(size_t x) /* x is a number of longs */
122 {
123   GEN z = ((GEN) avma) - x;
124   CHECK_CTRLC
125   if (x > (avma-pari_mainstack->bot) / sizeof(long))
126     new_chunk_resize(x);
127   set_avma((pari_sp)z);
128 #ifdef MEMSTEP
129   if (DEBUGMEM>1 && pari_mainstack->memused != DISABLE_MEMUSED) {
130     long d = (long)pari_mainstack->memused - (long)z;
131     if (labs(d) > 4*MEMSTEP)
132     {
133       pari_mainstack->memused = (pari_sp)z;
134       err_printf("...%4.0lf Mbytes used\n",
135                 (pari_mainstack->top-pari_mainstack->memused)/1048576.);
136     }
137   }
138 #endif
139   return z;
140 }
141 
142 INLINE char *
stack_malloc(size_t N)143 stack_malloc(size_t N)
144 {
145   long n = nchar2nlong(N);
146   return (char*)new_chunk(n);
147 }
148 
149 INLINE char *
stack_malloc_align(size_t N,long k)150 stack_malloc_align(size_t N, long k)
151 {
152   ulong d = ((ulong)avma) % k, e = ((ulong)N) % k;
153   if (d) (void)new_chunk(d/sizeof(long));
154   if (e) N += k-e;
155   return (char*) new_chunk(nchar2nlong(N));
156 }
157 
158 INLINE char *
stack_calloc(size_t N)159 stack_calloc(size_t N)
160 {
161   char *p = stack_malloc(N);
162   memset(p, 0, N); return p;
163 }
164 
165 INLINE char *
stack_calloc_align(size_t N,long k)166 stack_calloc_align(size_t N, long k)
167 {
168   ulong d = ((ulong)avma) % k, e = ((ulong)N) % k;
169   if (d) (void)new_chunk(d/sizeof(long));
170   if (e) N += k-e;
171   return stack_calloc(N);
172 }
173 
174 /* cgetg(lg(x), typ(x)), set *lx. Implicit unsetisclone() */
175 INLINE GEN
cgetg_copy(GEN x,long * plx)176 cgetg_copy(GEN x, long *plx) {
177   GEN y;
178   *plx = lg(x); y = new_chunk((size_t)*plx);
179   y[0] = x[0] & (TYPBITS|LGBITS); return y;
180 }
181 INLINE GEN
cgetg_block(long x,long y)182 cgetg_block(long x, long y)
183 {
184   GEN z = newblock((size_t)x);
185   z[0] = CLONEBIT | evaltyp(y) | evallg(x);
186   return z;
187 }
188 INLINE GEN
cgetg(long x,long y)189 cgetg(long x, long y)
190 {
191   GEN z = new_chunk((size_t)x);
192   z[0] = evaltyp(y) | evallg(x);
193   return z;
194 }
195 INLINE GEN
cgeti(long x)196 cgeti(long x)
197 {
198   GEN z = new_chunk((size_t)x);
199   z[0] = evaltyp(t_INT) | evallg(x);
200   return z;
201 }
202 INLINE GEN
cgetipos(long x)203 cgetipos(long x)
204 {
205   GEN z = cgeti(x);
206   z[1] = evalsigne(1) | evallgefint(x);
207   return z;
208 }
209 INLINE GEN
cgetineg(long x)210 cgetineg(long x)
211 {
212   GEN z = cgeti(x);
213   z[1] = evalsigne(-1) | evallgefint(x);
214   return z;
215 }
216 INLINE GEN
cgetr_block(long x)217 cgetr_block(long x)
218 {
219   GEN z = newblock((size_t)x);
220   z[0] = CLONEBIT | evaltyp(t_REAL) | evallg(x);
221   return z;
222 }
223 INLINE GEN
cgetr(long x)224 cgetr(long x)
225 {
226   GEN z = new_chunk((size_t)x);
227   z[0] = evaltyp(t_REAL) | evallg(x);
228   return z;
229 }
230 
231 /*******************************************************************/
232 /*                                                                 */
233 /*                     COPY, NEGATION, ABSOLUTE VALUE              */
234 /*                                                                 */
235 /*******************************************************************/
236 /* cannot do memcpy because sometimes x and y overlap */
237 INLINE GEN
leafcopy(GEN x)238 leafcopy(GEN x)
239 {
240   register long lx = lg(x);
241   GEN y = new_chunk(lx); /* can't use cgetg_copy, in case x,y overlap */
242   while (--lx > 0) y[lx] = x[lx];
243   y[0] = x[0] & (TYPBITS|LGBITS); return y;
244 }
245 INLINE GEN
icopy(GEN x)246 icopy(GEN x)
247 {
248   long i = lgefint(x), lx = i;
249   GEN y = new_chunk(lx); /* can't use cgeti, in case x,y overlap */
250   while (--i > 0) y[i] = x[i];
251   y[0] = evaltyp(t_INT) | evallg(lx);
252   return y;
253 }
254 INLINE GEN
icopyspec(GEN x,long nx)255 icopyspec(GEN x, long nx)
256 {
257   long i = nx+2, lx = i;
258   GEN y = new_chunk(lx); /* can't use cgeti, in case x,y overlap */
259   x -= 2; while (--i >= 2) y[i] = x[i];
260   y[1] = evalsigne(1) | evallgefint(lx);
261   y[0] = evaltyp(t_INT) | evallg(lx);
262   return y;
263 }
rcopy(GEN x)264 INLINE GEN rcopy(GEN x) { return leafcopy(x); }
mpcopy(GEN x)265 INLINE GEN mpcopy(GEN x) { return leafcopy(x); }
266 
267 INLINE GEN
mpabs(GEN x)268 mpabs(GEN x) { GEN y = leafcopy(x); setabssign(y); return y; }
269 INLINE GEN
mpabs_shallow(GEN x)270 mpabs_shallow(GEN x) { return signe(x) < 0? mpabs(x): x; }
absi(GEN x)271 INLINE GEN absi(GEN x) { return mpabs(x); }
absi_shallow(GEN x)272 INLINE GEN absi_shallow(GEN x) { return signe(x) < 0? negi(x): x; }
absr(GEN x)273 INLINE GEN absr(GEN x) { return mpabs(x); }
274 
275 INLINE GEN
mpneg(GEN x)276 mpneg(GEN x) { GEN y = leafcopy(x); togglesign(y); return y; }
negi(GEN x)277 INLINE GEN negi(GEN x) { return mpneg(x); }
negr(GEN x)278 INLINE GEN negr(GEN x) { return mpneg(x); }
279 
280 /* negate in place */
281 INLINE void
togglesign(GEN x)282 togglesign(GEN x) { if (x[1] & SIGNBITS) { x[1] ^= HIGHBIT; } }
283 INLINE void
setabssign(GEN x)284 setabssign(GEN x) { x[1] &= ~HIGHBIT; }
285 /* negate in place, except universal constants */
286 INLINE void
togglesign_safe(GEN * px)287 togglesign_safe(GEN *px)
288 {
289   switch(*px - gen_1) /* gen_1, gen_2, gen_m1, gen_m2 */
290   {
291     case 0: *px = gen_m1; break;
292     case 3: *px = gen_m2;  break;
293     case 6: *px = gen_1; break;
294     case 9: *px = gen_2;  break;
295     default: togglesign(*px);
296   }
297 }
298 /* setsigne(y, signe(x)) */
299 INLINE void
affectsign(GEN x,GEN y)300 affectsign(GEN x, GEN y)
301 {
302   y[1] = (x[1] & SIGNBITS) | (y[1] & ~SIGNBITS);
303 }
304 /* copies sign in place, except for universal constants */
305 INLINE void
affectsign_safe(GEN x,GEN * py)306 affectsign_safe(GEN x, GEN *py)
307 {
308   if (((*py)[1] ^ x[1]) & HIGHBIT) togglesign_safe(py);
309 }
310 /*******************************************************************/
311 /*                                                                 */
312 /*                     GEN -> LONG, LONG -> GEN                    */
313 /*                                                                 */
314 /*******************************************************************/
315 /* assume x != 0, return -x as a t_INT */
316 INLINE GEN
utoineg(ulong x)317 utoineg(ulong x) { GEN y = cgetineg(3); y[2] = x; return y; }
318 /* assume x != 0, return utoi(x) */
319 INLINE GEN
utoipos(ulong x)320 utoipos(ulong x) { GEN y = cgetipos(3); y[2] = x; return y; }
321 INLINE GEN
utoi(ulong x)322 utoi(ulong x) { return x? utoipos(x): gen_0; }
323 INLINE GEN
stoi(long x)324 stoi(long x)
325 {
326   if (!x) return gen_0;
327   return x > 0? utoipos((ulong)x): utoineg((ulong)-x);
328 }
329 
330 /* x 2^BIL + y */
331 INLINE GEN
uutoi(ulong x,ulong y)332 uutoi(ulong x, ulong y)
333 {
334   GEN z;
335   if (!x) return utoi(y);
336   z = cgetipos(4);
337   *int_W_lg(z, 1, 4) = x;
338   *int_W_lg(z, 0, 4) = y; return z;
339 }
340 /* - (x 2^BIL + y) */
341 INLINE GEN
uutoineg(ulong x,ulong y)342 uutoineg(ulong x, ulong y)
343 {
344   GEN z;
345   if (!x) return y? utoineg(y): gen_0;
346   z = cgetineg(4);
347   *int_W_lg(z, 1, 4) = x;
348   *int_W_lg(z, 0, 4) = y; return z;
349 }
350 
351 INLINE long
itos(GEN x)352 itos(GEN x)
353 {
354   long s = signe(x);
355   long u;
356 
357   if (!s) return 0;
358   u = x[2];
359   if (lgefint(x) > 3 || u < 0)
360     pari_err_OVERFLOW("t_INT-->long assignment");
361   return (s>0) ? u : -u;
362 }
363 /* as itos, but return 0 if too large. Cf is_bigint */
364 INLINE long
itos_or_0(GEN x)365 itos_or_0(GEN x) {
366   long n;
367   if (lgefint(x) != 3 || (n = x[2]) & HIGHBIT) return 0;
368   return signe(x) > 0? n: -n;
369 }
370 INLINE ulong
itou(GEN x)371 itou(GEN x)
372 {
373   switch(lgefint(x)) {
374     case 2: return 0;
375     case 3: return x[2];
376     default:
377       pari_err_OVERFLOW("t_INT-->ulong assignment");
378       return 0; /* LCOV_EXCL_LINE */
379   }
380 }
381 
382 /* as itou, but return 0 if too large. Cf is_bigint */
383 INLINE ulong
itou_or_0(GEN x)384 itou_or_0(GEN x) {
385   if (lgefint(x) != 3) return 0;
386   return (ulong)x[2];
387 }
388 
389 INLINE ulong
umuluu_or_0(ulong x,ulong y)390 umuluu_or_0(ulong x, ulong y)
391 {
392   ulong z;
393   LOCAL_HIREMAINDER;
394   z = mulll(x, y);
395   return hiremainder? 0: z;
396 }
397 /* return x*y if <= n, else 0. Beware overflow */
398 INLINE ulong
umuluu_le(ulong x,ulong y,ulong n)399 umuluu_le(ulong x, ulong y, ulong n)
400 {
401   ulong z;
402   LOCAL_HIREMAINDER;
403   z = mulll(x, y);
404   return (hiremainder || z > n)? 0: z;
405 }
406 
407 INLINE GEN
real_0_bit(long bitprec)408 real_0_bit(long bitprec) { GEN x=cgetr(2); x[1]=evalexpo(bitprec); return x; }
409 INLINE GEN
real_0(long prec)410 real_0(long prec) { return real_0_bit(-prec2nbits(prec)); }
411 INLINE GEN
real_1_bit(long bit)412 real_1_bit(long bit) { return real_1(nbits2prec(bit)); }
413 INLINE GEN
real_1(long prec)414 real_1(long prec) {
415   GEN x = cgetr(prec);
416   long i;
417   x[1] = evalsigne(1) | _evalexpo(0);
418   x[2] = (long)HIGHBIT; for (i=3; i<prec; i++) x[i] = 0;
419   return x;
420 }
421 INLINE GEN
real_m1(long prec)422 real_m1(long prec) {
423   GEN x = cgetr(prec);
424   long i;
425   x[1] = evalsigne(-1) | _evalexpo(0);
426   x[2] = (long)HIGHBIT; for (i=3; i<prec; i++) x[i] = 0;
427   return x;
428 }
429 
430 /* 2.^n */
431 INLINE GEN
real2n(long n,long prec)432 real2n(long n, long prec) { GEN z = real_1(prec); setexpo(z, n); return z; }
433 INLINE GEN
real_m2n(long n,long prec)434 real_m2n(long n, long prec) { GEN z = real_m1(prec); setexpo(z, n); return z; }
435 INLINE GEN
stor(long s,long prec)436 stor(long s, long prec) { GEN z = cgetr(prec); affsr(s,z); return z; }
437 INLINE GEN
utor(ulong s,long prec)438 utor(ulong s, long prec){ GEN z = cgetr(prec); affur(s,z); return z; }
439 INLINE GEN
itor(GEN x,long prec)440 itor(GEN x, long prec) { GEN z = cgetr(prec); affir(x,z); return z; }
441 INLINE GEN
rtor(GEN x,long prec)442 rtor(GEN x, long prec) { GEN z = cgetr(prec); affrr(x,z); return z; }
443 
int_bit(GEN x,long n)444 INLINE ulong int_bit(GEN x, long n)
445 {
446   long r, q = dvmdsBIL(n, &r);
447   return q < lgefint(x)-2?((ulong)*int_W(x,q) >> r) & 1UL:0;
448 }
449 
450 /*******************************************************************/
451 /*                                                                 */
452 /*                           COMPARISON                            */
453 /*                                                                 */
454 /*******************************************************************/
455 INLINE int
cmpss(long a,long b)456 cmpss(long a, long b)
457 { return a>b? 1: (a<b? -1: 0); }
458 
459 INLINE int
cmpuu(ulong a,ulong b)460 cmpuu(ulong a, ulong b)
461 { return a>b? 1: (a<b? -1: 0); }
462 
463 INLINE int
cmpir(GEN x,GEN y)464 cmpir(GEN x, GEN y)
465 {
466   pari_sp av;
467   GEN z;
468 
469   if (!signe(x)) return -signe(y);
470   if (!signe(y))
471   {
472     if (expo(y) >= expi(x)) return 0;
473     return signe(x);
474   }
475   av=avma; z = itor(x, realprec(y)); set_avma(av);
476   return cmprr(z,y); /* cmprr does no memory adjustment */
477 }
478 INLINE int
cmpri(GEN x,GEN y)479 cmpri(GEN x, GEN y) { return -cmpir(y,x); }
480 INLINE int
cmpsr(long x,GEN y)481 cmpsr(long x, GEN y)
482 {
483   pari_sp av;
484   GEN z;
485 
486   if (!x) return -signe(y);
487   av=avma; z = stor(x, LOWDEFAULTPREC); set_avma(av);
488   return cmprr(z,y);
489 }
490 INLINE int
cmprs(GEN x,long y)491 cmprs(GEN x, long y) { return -cmpsr(y,x); }
492 /* compare x and y */
493 INLINE int
cmpui(ulong x,GEN y)494 cmpui(ulong x, GEN y)
495 {
496   ulong p;
497   if (!x) return -signe(y);
498   if (signe(y) <= 0) return 1;
499   if (lgefint(y) > 3) return -1;
500   p = y[2]; if (p == x) return 0;
501   return p < x ? 1 : -1;
502 }
503 INLINE int
cmpiu(GEN x,ulong y)504 cmpiu(GEN x, ulong y) { return -cmpui(y,x); }
505 /* compare x and |y| */
506 INLINE int
abscmpui(ulong x,GEN y)507 abscmpui(ulong x, GEN y)
508 {
509   long l = lgefint(y);
510   ulong p;
511 
512   if (!x) return (l > 2)? -1: 0;
513   if (l == 2) return 1;
514   if (l > 3) return -1;
515   p = y[2]; if (p == x) return 0;
516   return p < x ? 1 : -1;
517 }
518 INLINE int
abscmpiu(GEN x,ulong y)519 abscmpiu(GEN x, ulong y) { return -abscmpui(y,x); }
520 INLINE int
cmpsi(long x,GEN y)521 cmpsi(long x, GEN y)
522 {
523   ulong p;
524 
525   if (!x) return -signe(y);
526 
527   if (x > 0)
528   {
529     if (signe(y)<=0) return 1;
530     if (lgefint(y)>3) return -1;
531     p = y[2]; if (p == (ulong)x) return 0;
532     return p < (ulong)x ? 1 : -1;
533   }
534 
535   if (signe(y)>=0) return -1;
536   if (lgefint(y)>3) return 1;
537   p = y[2]; if (p == (ulong)-x) return 0;
538   return p < (ulong)(-x) ? -1 : 1;
539 }
540 INLINE int
cmpis(GEN x,long y)541 cmpis(GEN x, long y) { return -cmpsi(y,x); }
542 INLINE int
mpcmp(GEN x,GEN y)543 mpcmp(GEN x, GEN y)
544 {
545   if (typ(x)==t_INT)
546     return (typ(y)==t_INT) ? cmpii(x,y) : cmpir(x,y);
547   return (typ(y)==t_INT) ? -cmpir(y,x) : cmprr(x,y);
548 }
549 
550 /* x == y ? */
551 INLINE int
equalui(ulong x,GEN y)552 equalui(ulong x, GEN y)
553 {
554   if (!x) return !signe(y);
555   if (signe(y) <= 0 || lgefint(y) != 3) return 0;
556   return ((ulong)y[2] == (ulong)x);
557 }
558 /* x == y ? */
559 INLINE int
equalsi(long x,GEN y)560 equalsi(long x, GEN y)
561 {
562   if (!x) return !signe(y);
563   if (x > 0)
564   {
565     if (signe(y) <= 0 || lgefint(y) != 3) return 0;
566     return ((ulong)y[2] == (ulong)x);
567   }
568   if (signe(y) >= 0 || lgefint(y) != 3) return 0;
569   return ((ulong)y[2] == (ulong)-x);
570 }
571 /* x == |y| ? */
572 INLINE int
absequalui(ulong x,GEN y)573 absequalui(ulong x, GEN y)
574 {
575   if (!x) return !signe(y);
576   return (lgefint(y) == 3 && (ulong)y[2] == x);
577 }
578 INLINE int
absequaliu(GEN x,ulong y)579 absequaliu(GEN x, ulong y) { return absequalui(y,x); }
580 INLINE int
equalis(GEN x,long y)581 equalis(GEN x, long y) { return equalsi(y,x); }
582 INLINE int
equaliu(GEN x,ulong y)583 equaliu(GEN x, ulong y) { return equalui(y,x); }
584 
585 /* assume x != 0, is |x| == 2^n ? */
586 INLINE int
absrnz_equal2n(GEN x)587 absrnz_equal2n(GEN x) {
588   if ((ulong)x[2]==HIGHBIT)
589   {
590     long i, lx = lg(x);
591     for (i = 3; i < lx; i++)
592       if (x[i]) return 0;
593     return 1;
594   }
595   return 0;
596 }
597 /* assume x != 0, is |x| == 1 ? */
598 INLINE int
absrnz_equal1(GEN x)599 absrnz_equal1(GEN x) { return !expo(x) && absrnz_equal2n(x); }
600 
601 INLINE long
maxss(long x,long y)602 maxss(long x, long y) { return x>y?x:y; }
603 INLINE long
minss(long x,long y)604 minss(long x, long y) { return x<y?x:y; }
605 INLINE long
minuu(ulong x,ulong y)606 minuu(ulong x, ulong y) { return x<y?x:y; }
607 INLINE long
maxuu(ulong x,ulong y)608 maxuu(ulong x, ulong y) { return x>y?x:y; }
609 INLINE double
maxdd(double x,double y)610 maxdd(double x, double y) { return x>y?x:y; }
611 INLINE double
mindd(double x,double y)612 mindd(double x, double y) { return x<y?x:y; }
613 
614 /*******************************************************************/
615 /*                                                                 */
616 /*                             ADD / SUB                           */
617 /*                                                                 */
618 /*******************************************************************/
619 INLINE GEN
subuu(ulong x,ulong y)620 subuu(ulong x, ulong y)
621 {
622   ulong z;
623   LOCAL_OVERFLOW;
624   z = subll(x, y);
625   return overflow? utoineg(-z): utoi(z);
626 }
627 INLINE GEN
adduu(ulong x,ulong y)628 adduu(ulong x, ulong y) { ulong t = x+y; return uutoi((t < x), t); }
629 
630 INLINE GEN
addss(long x,long y)631 addss(long x, long y)
632 {
633   if (!x) return stoi(y);
634   if (!y) return stoi(x);
635   if (x > 0) return y > 0? adduu(x,y): subuu(x, -y);
636 
637   if (y > 0) return subuu(y, -x);
638   else { /* - adduu(-x, -y) */
639     ulong t = (-x)+(-y); return uutoineg((t < (ulong)(-x)), t);
640   }
641 }
subss(long x,long y)642 INLINE GEN subss(long x, long y) { return addss(-y,x); }
643 
644 INLINE GEN
subii(GEN x,GEN y)645 subii(GEN x, GEN y)
646 {
647   if (x==y) return gen_0; /* frequent with x = y = gen_0 */
648   return addii_sign(x, signe(x), y, -signe(y));
649 }
650 INLINE GEN
addii(GEN x,GEN y)651 addii(GEN x, GEN y) { return addii_sign(x, signe(x), y, signe(y)); }
652 INLINE GEN
addrr(GEN x,GEN y)653 addrr(GEN x, GEN y) { return addrr_sign(x, signe(x), y, signe(y)); }
654 INLINE GEN
subrr(GEN x,GEN y)655 subrr(GEN x, GEN y) { return addrr_sign(x, signe(x), y, -signe(y)); }
656 INLINE GEN
addir(GEN x,GEN y)657 addir(GEN x, GEN y) { return addir_sign(x, signe(x), y, signe(y)); }
658 INLINE GEN
subir(GEN x,GEN y)659 subir(GEN x, GEN y) { return addir_sign(x, signe(x), y, -signe(y)); }
660 INLINE GEN
subri(GEN x,GEN y)661 subri(GEN x, GEN y) { return addir_sign(y, -signe(y), x, signe(x)); }
662 INLINE GEN
addsi(long x,GEN y)663 addsi(long x, GEN y) { return addsi_sign(x, y, signe(y)); }
664 INLINE GEN
addui(ulong x,GEN y)665 addui(ulong x, GEN y) { return addui_sign(x, y, signe(y)); }
666 INLINE GEN
subsi(long x,GEN y)667 subsi(long x, GEN y) { return addsi_sign(x, y, -signe(y)); }
668 INLINE GEN
subui(ulong x,GEN y)669 subui(ulong x, GEN y) { return addui_sign(x, y, -signe(y)); }
670 
671 /*******************************************************************/
672 /*                                                                 */
673 /*                           MOD, REM, DIV                         */
674 /*                                                                 */
675 /*******************************************************************/
mod2BIL(GEN x)676 INLINE ulong mod2BIL(GEN x) { return *int_LSW(x); }
mod64(GEN x)677 INLINE long mod64(GEN x) { return mod2BIL(x) & 63; }
mod32(GEN x)678 INLINE long mod32(GEN x) { return mod2BIL(x) & 31; }
mod16(GEN x)679 INLINE long mod16(GEN x) { return mod2BIL(x) & 15; }
mod8(GEN x)680 INLINE long mod8(GEN x)  { return mod2BIL(x) & 7; }
mod4(GEN x)681 INLINE long mod4(GEN x)  { return mod2BIL(x) & 3; }
mod2(GEN x)682 INLINE long mod2(GEN x)  { return mod2BIL(x) & 1; }
683 INLINE int
mpodd(GEN x)684 mpodd(GEN x) { return signe(x) && mod2(x); }
685 /* x mod 2^n, n < BITS_IN_LONG */
686 INLINE ulong
umodi2n(GEN x,long n)687 umodi2n(GEN x, long n)
688 {
689   long s = signe(x);
690   const ulong _2n = 1UL << n;
691   ulong m;
692   if (!s) return 0;
693   m = *int_LSW(x) & (_2n - 1);
694   if (s < 0 && m) m = _2n - m;
695   return m;
696 }
Mod64(GEN x)697 INLINE ulong Mod64(GEN x){ return umodi2n(x,6); }
Mod32(GEN x)698 INLINE ulong Mod32(GEN x){ return umodi2n(x,5); }
Mod16(GEN x)699 INLINE ulong Mod16(GEN x){ return umodi2n(x,4); }
Mod8(GEN x)700 INLINE ulong Mod8(GEN x) { return umodi2n(x,3); }
Mod4(GEN x)701 INLINE ulong Mod4(GEN x) { return umodi2n(x,2); }
Mod2(GEN x)702 INLINE ulong Mod2(GEN x) { return umodi2n(x,1); }
703 
704 INLINE GEN
truedivii(GEN a,GEN b)705 truedivii(GEN a,GEN b) { return truedvmdii(a,b,NULL); }
706 INLINE GEN
truedivis(GEN a,long b)707 truedivis(GEN a, long b) { return truedvmdis(a,b,NULL); }
708 INLINE GEN
truedivsi(long a,GEN b)709 truedivsi(long a, GEN b) { return truedvmdsi(a,b,NULL); }
710 
711 INLINE GEN
divii(GEN a,GEN b)712 divii(GEN a, GEN b) { return dvmdii(a,b,NULL); }
713 INLINE GEN
remii(GEN a,GEN b)714 remii(GEN a, GEN b) { return dvmdii(a,b,ONLY_REM); }
715 
716 INLINE GEN
divss(long x,long y)717 divss(long x, long y) { return stoi(x / y); }
718 INLINE GEN
modss(long x,long y)719 modss(long x, long y) { return utoi(smodss(x, y)); }
720 INLINE GEN
remss(long x,long y)721 remss(long x, long y) { return stoi(x % y); }
722 INLINE long
smodss(long x,long y)723 smodss(long x, long y)
724 {
725   long r = x%y;
726   return (r >= 0)? r: labs(y) + r;
727 }
728 INLINE ulong
umodsu(long x,ulong y)729 umodsu(long x, ulong y)
730 {
731   return x>=0 ? x%y: Fl_neg((-x)%y, y);
732 }
733 
734 INLINE long
sdivss_rem(long x,long y,long * r)735 sdivss_rem(long x, long y, long *r)
736 {
737   long q;
738   LOCAL_HIREMAINDER;
739   if (!y) pari_err_INV("sdivss_rem",gen_0);
740   hiremainder = 0; q = divll((ulong)labs(x),(ulong)labs(y));
741   if (x < 0) { hiremainder = -((long)hiremainder); q = -q; }
742   if (y < 0) q = -q;
743   *r = hiremainder; return q;
744 }
745 INLINE GEN
divss_rem(long x,long y,long * r)746 divss_rem(long x, long y, long *r) { return stoi(sdivss_rem(x,y,r)); }
747 INLINE ulong
udivuu_rem(ulong x,ulong y,ulong * r)748 udivuu_rem(ulong x, ulong y, ulong *r)
749 {
750   if (!y) pari_err_INV("udivuu_rem",gen_0);
751   *r = x % y; return x / y;
752 }
753 INLINE ulong
ceildivuu(ulong a,ulong b)754 ceildivuu(ulong a, ulong b)
755 {
756   ulong c = a/b;
757   return (a%b)? c+1: c;
758 }
759 
760 INLINE ulong
uabsdivui_rem(ulong x,GEN y,ulong * r)761 uabsdivui_rem(ulong x, GEN y, ulong *r)
762 {
763   long q, s = signe(y);
764   LOCAL_HIREMAINDER;
765 
766   if (!s) pari_err_INV("uabsdivui_rem",gen_0);
767   if (!x || lgefint(y)>3) { *r = x; return 0; }
768   hiremainder=0; q = (long)divll(x, (ulong)y[2]);
769   if (s < 0) q = -q;
770   *r = hiremainder; return q;
771 }
772 
773 /* assume d != 0 and |n| / d can be represented as an ulong.
774  * Return |n|/d, set *r = |n| % d */
775 INLINE ulong
uabsdiviu_rem(GEN n,ulong d,ulong * r)776 uabsdiviu_rem(GEN n, ulong d, ulong *r)
777 {
778   switch(lgefint(n))
779   {
780     case 2: *r = 0; return 0;
781     case 3:
782     {
783       ulong nn = n[2];
784       *r = nn % d; return nn / d;
785     }
786     default: /* 4 */
787     {
788       ulong n1, n0, q;
789       LOCAL_HIREMAINDER;
790       n0 = *int_W(n,0);
791       n1 = *int_W(n,1);
792       hiremainder = n1;
793       q = divll(n0, d);
794       *r = hiremainder; return q;
795     }
796   }
797 }
798 
799 INLINE long
sdivsi_rem(long x,GEN y,long * r)800 sdivsi_rem(long x, GEN y, long *r)
801 {
802   long q, s = signe(y);
803   LOCAL_HIREMAINDER;
804 
805   if (!s) pari_err_INV("sdivsi_rem",gen_0);
806   if (!x || lgefint(y)>3 || ((long)y[2]) < 0) { *r = x; return 0; }
807   hiremainder=0; q = (long)divll(labs(x), (ulong)y[2]);
808   if (x < 0) { hiremainder = -((long)hiremainder); q = -q; }
809   if (s < 0) q = -q;
810   *r = hiremainder; return q;
811 }
812 INLINE GEN
divsi_rem(long s,GEN y,long * r)813 divsi_rem(long s, GEN y, long *r) { return stoi(sdivsi_rem(s,y,r)); }
814 
815 INLINE long
sdivsi(long x,GEN y)816 sdivsi(long x, GEN y)
817 {
818   long q, s = signe(y);
819 
820   if (!s) pari_err_INV("sdivsi",gen_0);
821   if (!x || lgefint(y)>3 || ((long)y[2]) < 0) return 0;
822   q = labs(x) / y[2];
823   if (x < 0) q = -q;
824   if (s < 0) q = -q;
825   return q;
826 }
827 
828 INLINE GEN
dvmdss(long x,long y,GEN * z)829 dvmdss(long x, long y, GEN *z)
830 {
831   long r;
832   GEN q = divss_rem(x,y, &r);
833   *z = stoi(r); return q;
834 }
835 INLINE long
dvmdsBIL(long n,long * r)836 dvmdsBIL(long n, long *r) { *r = remsBIL(n); return divsBIL(n); }
837 INLINE ulong
dvmduBIL(ulong n,ulong * r)838 dvmduBIL(ulong n, ulong *r) { *r = remsBIL(n); return divsBIL(n); }
839 INLINE GEN
dvmdsi(long x,GEN y,GEN * z)840 dvmdsi(long x, GEN y, GEN *z)
841 {
842   long r;
843   GEN q = divsi_rem(x,y, &r);
844   *z = stoi(r); return q;
845 }
846 INLINE GEN
dvmdis(GEN x,long y,GEN * z)847 dvmdis(GEN x, long y, GEN *z)
848 {
849   long r;
850   GEN q = divis_rem(x,y, &r);
851   *z = stoi(r); return q;
852 }
853 
854 INLINE long
smodis(GEN x,long y)855 smodis(GEN x, long y)
856 {
857   pari_sp av = avma;
858   long r; (void)divis_rem(x,y, &r);
859   return gc_long(av, (r >= 0)? r: labs(y) + r);
860 }
861 INLINE GEN
modis(GEN x,long y)862 modis(GEN x, long y) { return stoi(smodis(x,y)); }
863 INLINE GEN
modsi(long x,GEN y)864 modsi(long x, GEN y) {
865   long r; (void)sdivsi_rem(x, y, &r);
866   return (r >= 0)? stoi(r): addsi_sign(r, y, 1);
867 }
868 
869 INLINE ulong
umodui(ulong x,GEN y)870 umodui(ulong x, GEN y)
871 {
872   if (!signe(y)) pari_err_INV("umodui",gen_0);
873   if (!x || lgefint(y) > 3) return x;
874   return x % (ulong)y[2];
875 }
876 
877 INLINE ulong
ugcdiu(GEN x,ulong y)878 ugcdiu(GEN x, ulong y) { return ugcd(umodiu(x,y), y); }
879 INLINE ulong
ugcdui(ulong y,GEN x)880 ugcdui(ulong y, GEN x) { return ugcd(umodiu(x,y), y); }
881 
882 INLINE GEN
remsi(long x,GEN y)883 remsi(long x, GEN y)
884 { long r; (void)sdivsi_rem(x,y, &r); return stoi(r); }
885 INLINE GEN
remis(GEN x,long y)886 remis(GEN x, long y)
887 {
888   pari_sp av = avma;
889   long r;
890   (void)divis_rem(x,y, &r); set_avma(av); return stoi(r);
891 }
892 
893 INLINE GEN
rdivis(GEN x,long y,long prec)894 rdivis(GEN x, long y, long prec)
895 {
896   GEN z = cgetr(prec);
897   pari_sp av = avma;
898   affrr(divrs(itor(x,prec), y),z);
899   set_avma(av); return z;
900 }
901 INLINE GEN
rdivsi(long x,GEN y,long prec)902 rdivsi(long x, GEN y, long prec)
903 {
904   GEN z = cgetr(prec);
905   pari_sp av = avma;
906   affrr(divsr(x, itor(y,prec)), z);
907   set_avma(av); return z;
908 }
909 INLINE GEN
rdivss(long x,long y,long prec)910 rdivss(long x, long y, long prec)
911 {
912   GEN z = cgetr(prec);
913   pari_sp av = avma;
914   affrr(divrs(stor(x, prec), y), z);
915   set_avma(av); return z;
916 }
917 
918 INLINE void
rdiviiz(GEN x,GEN y,GEN z)919 rdiviiz(GEN x, GEN y, GEN z)
920 {
921   long prec = realprec(z), lx = lgefint(x), ly = lgefint(y);
922   if (lx == 2) { affur(0, z); return; }
923   if (ly == 3)
924   {
925     affir(x, z); if (signe(y) < 0) togglesign(z);
926     affrr(divru(z, y[2]), z);
927   }
928   else if (lx > prec + 1 || ly > prec + 1)
929   {
930     affir(x,z); affrr(divri(z, y), z);
931   }
932   else
933   {
934     long b = bit_accuracy(prec) + expi(y) - expi(x) + 1;
935     GEN q = divii(b > 0? shifti(x, b): x, y);
936     affir(q, z); if (b > 0) shiftr_inplace(z, -b);
937   }
938   set_avma((ulong)z);
939 }
940 INLINE GEN
rdivii(GEN x,GEN y,long prec)941 rdivii(GEN x, GEN y, long prec)
942 { GEN z = cgetr(prec); rdiviiz(x, y, z); return z; }
943 INLINE GEN
fractor(GEN x,long prec)944 fractor(GEN x, long prec)
945 { return rdivii(gel(x,1), gel(x,2), prec); }
946 
947 INLINE int
dvdii(GEN x,GEN y)948 dvdii(GEN x, GEN y)
949 {
950   pari_sp av = avma;
951   GEN r;
952   if (!signe(x)) return 1;
953   if (!signe(y)) return 0;
954   r = remii(x,y);
955   return gc_bool(av, r == gen_0);
956 }
957 INLINE int
dvdsi(long x,GEN y)958 dvdsi(long x, GEN y)
959 {
960   if (x == 0) return 1;
961   if (!signe(y)) return 0;
962   if (lgefint(y) != 3) return 0;
963   return x % y[2] == 0;
964 }
965 INLINE int
dvdui(ulong x,GEN y)966 dvdui(ulong x, GEN y)
967 {
968   if (x == 0) return 1;
969   if (!signe(y)) return 0;
970   if (lgefint(y) != 3) return 0;
971   return x % y[2] == 0;
972 }
973 INLINE int
dvdis(GEN x,long y)974 dvdis(GEN x, long y)
975 { return y? smodis(x, y) == 0: signe(x) == 0; }
976 INLINE int
dvdiu(GEN x,ulong y)977 dvdiu(GEN x, ulong y)
978 { return y? umodiu(x, y) == 0: signe(x) == 0; }
979 
980 INLINE int
dvdisz(GEN x,long y,GEN z)981 dvdisz(GEN x, long y, GEN z)
982 {
983   const pari_sp av = avma;
984   long r;
985   GEN p1 = divis_rem(x,y, &r);
986   set_avma(av); if (r) return 0;
987   affii(p1,z); return 1;
988 }
989 INLINE int
dvdiuz(GEN x,ulong y,GEN z)990 dvdiuz(GEN x, ulong y, GEN z)
991 {
992   const pari_sp av = avma;
993   ulong r;
994   GEN p1 = absdiviu_rem(x,y, &r);
995   set_avma(av); if (r) return 0;
996   affii(p1,z); return 1;
997 }
998 INLINE int
dvdiiz(GEN x,GEN y,GEN z)999 dvdiiz(GEN x, GEN y, GEN z)
1000 {
1001   const pari_sp av=avma;
1002   GEN p2, p1 = dvmdii(x,y,&p2);
1003   if (signe(p2)) return gc_bool(av,0);
1004   affii(p1,z); return gc_bool(av,1);
1005 }
1006 
1007 INLINE ulong
remlll_pre(ulong u2,ulong u1,ulong u0,ulong n,ulong ninv)1008 remlll_pre(ulong u2, ulong u1, ulong u0, ulong n, ulong ninv)
1009 {
1010   u1 = remll_pre(u2, u1, n, ninv);
1011   return remll_pre(u1, u0, n, ninv);
1012 }
1013 
1014 INLINE ulong
Fl_sqr_pre(ulong a,ulong p,ulong pi)1015 Fl_sqr_pre(ulong a, ulong p, ulong pi)
1016 {
1017   register ulong x;
1018   LOCAL_HIREMAINDER;
1019   x = mulll(a,a);
1020   return remll_pre(hiremainder, x, p, pi);
1021 }
1022 
1023 INLINE ulong
Fl_mul_pre(ulong a,ulong b,ulong p,ulong pi)1024 Fl_mul_pre(ulong a, ulong b, ulong p, ulong pi)
1025 {
1026   register ulong x;
1027   LOCAL_HIREMAINDER;
1028   x = mulll(a,b);
1029   return remll_pre(hiremainder, x, p, pi);
1030 }
1031 
1032 INLINE ulong
Fl_addmul_pre(ulong y0,ulong x0,ulong x1,ulong p,ulong pi)1033 Fl_addmul_pre(ulong y0, ulong x0, ulong x1, ulong p, ulong pi)
1034 {
1035   ulong l0, h0;
1036   LOCAL_HIREMAINDER;
1037   hiremainder = y0;
1038   l0 = addmul(x0, x1); h0 = hiremainder;
1039   return remll_pre(h0, l0, p, pi);
1040 }
1041 
1042 INLINE ulong
Fl_addmulmul_pre(ulong x0,ulong y0,ulong x1,ulong y1,ulong p,ulong pi)1043 Fl_addmulmul_pre(ulong x0, ulong y0, ulong x1, ulong y1, ulong p, ulong pi)
1044 {
1045   ulong l0, l1, h0, h1;
1046   LOCAL_OVERFLOW;
1047   LOCAL_HIREMAINDER;
1048   l0 = mulll(x0, y0); h0 = hiremainder;
1049   l1 = mulll(x1, y1); h1 = hiremainder;
1050   l0 = addll(l0, l1); h0 = addllx(h0, h1);
1051   return overflow ? remlll_pre(1, h0, l0, p, pi): remll_pre(h0, l0, p, pi);
1052 }
1053 
1054 INLINE ulong
Fl_ellj_pre(ulong a4,ulong a6,ulong p,ulong pi)1055 Fl_ellj_pre(ulong a4, ulong a6, ulong p, ulong pi)
1056 {
1057   /* a43 = 4 a4^3 */
1058   ulong a43 = Fl_double(Fl_double(
1059               Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi), p), p);
1060   /* a62 = 27 a6^2 */
1061   ulong a62 = Fl_mul_pre(Fl_sqr_pre(a6, p, pi), 27 % p, p, pi);
1062   ulong z1 = Fl_mul_pre(a43, 1728 % p, p, pi);
1063   ulong z2 = Fl_add(a43, a62, p);
1064   return Fl_div(z1, z2, p);
1065 }
1066 
1067 /*******************************************************************/
1068 /*                                                                 */
1069 /*                        MP (INT OR REAL)                         */
1070 /*                                                                 */
1071 /*******************************************************************/
1072 INLINE GEN
mptrunc(GEN x)1073 mptrunc(GEN x) { return typ(x)==t_INT? icopy(x): truncr(x); }
1074 INLINE GEN
mpfloor(GEN x)1075 mpfloor(GEN x) { return typ(x)==t_INT? icopy(x): floorr(x); }
1076 INLINE GEN
mpceil(GEN x)1077 mpceil(GEN x) { return typ(x)==t_INT? icopy(x): ceilr(x); }
1078 INLINE GEN
mpround(GEN x)1079 mpround(GEN x) { return typ(x) == t_INT? icopy(x): roundr(x); }
1080 
1081 INLINE long
mpexpo(GEN x)1082 mpexpo(GEN x) { return typ(x) == t_INT? expi(x): expo(x); }
1083 
1084 INLINE GEN
mpadd(GEN x,GEN y)1085 mpadd(GEN x, GEN y)
1086 {
1087   if (typ(x)==t_INT)
1088     return (typ(y)==t_INT) ? addii(x,y) : addir(x,y);
1089   return (typ(y)==t_INT) ? addir(y,x) : addrr(x,y);
1090 }
1091 INLINE GEN
mpsub(GEN x,GEN y)1092 mpsub(GEN x, GEN y)
1093 {
1094   if (typ(x)==t_INT)
1095     return (typ(y)==t_INT) ? subii(x,y) : subir(x,y);
1096   return (typ(y)==t_INT) ? subri(x,y) : subrr(x,y);
1097 }
1098 INLINE GEN
mpmul(GEN x,GEN y)1099 mpmul(GEN x, GEN y)
1100 {
1101   if (typ(x)==t_INT)
1102     return (typ(y)==t_INT) ? mulii(x,y) : mulir(x,y);
1103   return (typ(y)==t_INT) ? mulir(y,x) : mulrr(x,y);
1104 }
1105 INLINE GEN
mpsqr(GEN x)1106 mpsqr(GEN x) { return (typ(x)==t_INT) ? sqri(x) : sqrr(x); }
1107 INLINE GEN
mpdiv(GEN x,GEN y)1108 mpdiv(GEN x, GEN y)
1109 {
1110   if (typ(x)==t_INT)
1111     return (typ(y)==t_INT) ? divii(x,y) : divir(x,y);
1112   return (typ(y)==t_INT) ? divri(x,y) : divrr(x,y);
1113 }
1114 
1115 /*******************************************************************/
1116 /*                                                                 */
1117 /*                          Z/nZ, n ULONG                          */
1118 /*                                                                 */
1119 /*******************************************************************/
1120 INLINE ulong
Fl_double(ulong a,ulong p)1121 Fl_double(ulong a, ulong p)
1122 {
1123   ulong res = a << 1;
1124   return (res >= p || res < a) ? res - p : res;
1125 }
1126 INLINE ulong
Fl_triple(ulong a,ulong p)1127 Fl_triple(ulong a, ulong p)
1128 {
1129   ulong res = a << 1;
1130   if (res >= p || res < a) res -= p;
1131   res += a;
1132   return (res >= p || res < a)? res - p: res;
1133 }
1134 INLINE ulong
Fl_halve(ulong a,ulong p)1135 Fl_halve(ulong a, ulong p)
1136 {
1137   ulong ap, ap2;
1138   if ((a&1UL)==0) return a>>1;
1139   ap = a + p; ap2 = ap>>1;
1140   return ap>=a ? ap2: (ap2|HIGHBIT);
1141 }
1142 
1143 INLINE ulong
Fl_add(ulong a,ulong b,ulong p)1144 Fl_add(ulong a, ulong b, ulong p)
1145 {
1146   ulong res = a + b;
1147   return (res >= p || res < a) ? res - p : res;
1148 }
1149 INLINE ulong
Fl_neg(ulong x,ulong p)1150 Fl_neg(ulong x, ulong p) { return x ? p - x: 0; }
1151 
1152 INLINE ulong
Fl_sub(ulong a,ulong b,ulong p)1153 Fl_sub(ulong a, ulong b, ulong p)
1154 {
1155   ulong res = a - b;
1156   return (res > a) ? res + p: res;
1157 }
1158 
1159 /* centerlift(u mod p) */
1160 INLINE long
Fl_center(ulong u,ulong p,ulong ps2)1161 Fl_center(ulong u, ulong p, ulong ps2) { return (long) (u > ps2)? u - p: u; }
1162 
1163 INLINE ulong
Fl_mul(ulong a,ulong b,ulong p)1164 Fl_mul(ulong a, ulong b, ulong p)
1165 {
1166   register ulong x;
1167   LOCAL_HIREMAINDER;
1168   x = mulll(a,b);
1169   if (!hiremainder) return x % p;
1170   (void)divll(x,p); return hiremainder;
1171 }
1172 INLINE ulong
Fl_sqr(ulong a,ulong p)1173 Fl_sqr(ulong a, ulong p)
1174 {
1175   register ulong x;
1176   LOCAL_HIREMAINDER;
1177   x = mulll(a,a);
1178   if (!hiremainder) return x % p;
1179   (void)divll(x,p); return hiremainder;
1180 }
1181 /* don't assume that p is prime: can't special case a = 0 */
1182 INLINE ulong
Fl_div(ulong a,ulong b,ulong p)1183 Fl_div(ulong a, ulong b, ulong p)
1184 { return Fl_mul(a, Fl_inv(b, p), p); }
1185 
1186 /*******************************************************************/
1187 /*                                                                 */
1188 /*        DEFINED FROM EXISTING ONE EXPLOITING COMMUTATIVITY       */
1189 /*                                                                 */
1190 /*******************************************************************/
1191 INLINE GEN
addri(GEN x,GEN y)1192 addri(GEN x, GEN y) { return addir(y,x); }
1193 INLINE GEN
addis(GEN x,long s)1194 addis(GEN x, long s) { return addsi(s,x); }
1195 INLINE GEN
addiu(GEN x,ulong s)1196 addiu(GEN x, ulong s) { return addui(s,x); }
1197 INLINE GEN
addrs(GEN x,long s)1198 addrs(GEN x, long s) { return addsr(s,x); }
1199 
1200 INLINE GEN
subiu(GEN x,long y)1201 subiu(GEN x, long y) { GEN z = subui(y, x); togglesign(z); return z; }
1202 INLINE GEN
subis(GEN x,long y)1203 subis(GEN x, long y) { return addsi(-y,x); }
1204 INLINE GEN
subrs(GEN x,long y)1205 subrs(GEN x, long y) { return addsr(-y,x); }
1206 
1207 INLINE GEN
mulis(GEN x,long s)1208 mulis(GEN x, long s) { return mulsi(s,x); }
1209 INLINE GEN
muliu(GEN x,ulong s)1210 muliu(GEN x, ulong s) { return mului(s,x); }
1211 INLINE GEN
mulru(GEN x,ulong s)1212 mulru(GEN x, ulong s) { return mulur(s,x); }
1213 INLINE GEN
mulri(GEN x,GEN s)1214 mulri(GEN x, GEN s) { return mulir(s,x); }
1215 INLINE GEN
mulrs(GEN x,long s)1216 mulrs(GEN x, long s) { return mulsr(s,x); }
1217 
1218 /*******************************************************************/
1219 /*                                                                 */
1220 /*                  VALUATION, EXPONENT, SHIFTS                    */
1221 /*                                                                 */
1222 /*******************************************************************/
1223 INLINE long
vali(GEN x)1224 vali(GEN x)
1225 {
1226   long i;
1227   GEN xp;
1228 
1229   if (!signe(x)) return -1;
1230   xp=int_LSW(x);
1231   for (i=0; !*xp; i++) xp=int_nextW(xp);
1232   return vals(*xp) + i * BITS_IN_LONG;
1233 }
1234 
1235 /* assume x > 0 */
1236 INLINE long
expu(ulong x)1237 expu(ulong x) { return (BITS_IN_LONG-1) - (long)bfffo(x); }
1238 
1239 INLINE long
expi(GEN x)1240 expi(GEN x)
1241 {
1242   const long lx=lgefint(x);
1243   return lx==2? -(long)HIGHEXPOBIT: bit_accuracy(lx)-(long)bfffo(*int_MSW(x))-1;
1244 }
1245 
1246 INLINE GEN
shiftr(GEN x,long n)1247 shiftr(GEN x, long n)
1248 {
1249   const long e = evalexpo(expo(x)+n);
1250   const GEN y = rcopy(x);
1251 
1252   if (e & ~EXPOBITS) pari_err_OVERFLOW("expo()");
1253   y[1] = (y[1]&~EXPOBITS) | e; return y;
1254 }
1255 INLINE GEN
mpshift(GEN x,long s)1256 mpshift(GEN x,long s) { return (typ(x)==t_INT)?shifti(x,s):shiftr(x,s); }
1257 
1258 /* FIXME: adapt/use mpn_[lr]shift instead */
1259 /* z2[imin..imax] := z1[imin..imax].f shifted left sh bits
1260  * (feeding f from the right). Assume sh > 0 */
1261 INLINE void
shift_left(GEN z2,GEN z1,long imin,long imax,ulong f,ulong sh)1262 shift_left(GEN z2, GEN z1, long imin, long imax, ulong f,  ulong sh)
1263 {
1264   GEN sb = z1 + imin, se = z1 + imax, te = z2 + imax;
1265   ulong l, m = BITS_IN_LONG - sh, k = f >> m;
1266   while (se > sb) {
1267     l     = *se--;
1268     *te-- = (l << sh) | k;
1269     k     = l >> m;
1270   }
1271   *te = (((ulong)*se) << sh) | k;
1272 }
1273 /* z2[imin..imax] := f.z1[imin..imax-1] shifted right sh bits
1274  * (feeding f from the left). Assume sh > 0 */
1275 INLINE void
shift_right(GEN z2,GEN z1,long imin,long imax,ulong f,ulong sh)1276 shift_right(GEN z2, GEN z1, long imin, long imax, ulong f, ulong sh)
1277 {
1278   GEN sb = z1 + imin, se = z1 + imax, tb = z2 + imin;
1279   ulong k, l = *sb++, m = BITS_IN_LONG - sh;
1280   *tb++ = (l >> sh) | (f << m);
1281   while (sb < se) {
1282     k     = l << m;
1283     l     = *sb++;
1284     *tb++ = (l >> sh) | k;
1285   }
1286 }
1287 
1288 /* Backward compatibility. Inefficient && unused */
1289 extern ulong hiremainder;
1290 INLINE ulong
shiftl(ulong x,ulong y)1291 shiftl(ulong x, ulong y)
1292 { hiremainder = x>>(BITS_IN_LONG-y); return (x<<y); }
1293 
1294 INLINE ulong
shiftlr(ulong x,ulong y)1295 shiftlr(ulong x, ulong y)
1296 { hiremainder = x<<(BITS_IN_LONG-y); return (x>>y); }
1297 
1298 INLINE void
shiftr_inplace(GEN z,long d)1299 shiftr_inplace(GEN z, long d)
1300 {
1301   setexpo(z, expo(z)+d);
1302 }
1303 
1304 /*******************************************************************/
1305 /*                                                                 */
1306 /*                           ASSIGNMENT                            */
1307 /*                                                                 */
1308 /*******************************************************************/
1309 INLINE void
affii(GEN x,GEN y)1310 affii(GEN x, GEN y)
1311 {
1312   long lx = lgefint(x);
1313   if (lg(y)<lx) pari_err_OVERFLOW("t_INT-->t_INT assignment");
1314   while (--lx) y[lx] = x[lx];
1315 }
1316 INLINE void
affsi(long s,GEN x)1317 affsi(long s, GEN x)
1318 {
1319   if (!s) x[1] = evalsigne(0) | evallgefint(2);
1320   else
1321   {
1322     if (s > 0) { x[1] = evalsigne( 1) | evallgefint(3); x[2] =  s; }
1323     else       { x[1] = evalsigne(-1) | evallgefint(3); x[2] = -s; }
1324   }
1325 }
1326 INLINE void
affui(ulong u,GEN x)1327 affui(ulong u, GEN x)
1328 {
1329   if (!u) x[1] = evalsigne(0) | evallgefint(2);
1330   else  { x[1] = evalsigne(1) | evallgefint(3); x[2] = u; }
1331 }
1332 
1333 INLINE void
affsr(long x,GEN y)1334 affsr(long x, GEN y)
1335 {
1336   long sh, i, ly = lg(y);
1337 
1338   if (!x)
1339   {
1340     y[1] = evalexpo(-prec2nbits(ly));
1341     return;
1342   }
1343   if (x < 0) {
1344     x = -x; sh = bfffo(x);
1345     y[1] = evalsigne(-1) | _evalexpo((BITS_IN_LONG-1)-sh);
1346   }
1347   else
1348   {
1349     sh = bfffo(x);
1350     y[1] = evalsigne(1) | _evalexpo((BITS_IN_LONG-1)-sh);
1351   }
1352   y[2] = ((ulong)x)<<sh; for (i=3; i<ly; i++) y[i]=0;
1353 }
1354 
1355 INLINE void
affur(ulong x,GEN y)1356 affur(ulong x, GEN y)
1357 {
1358   long sh, i, ly = lg(y);
1359 
1360   if (!x)
1361   {
1362     y[1] = evalexpo(-prec2nbits(ly));
1363     return;
1364   }
1365   sh = bfffo(x);
1366   y[1] = evalsigne(1) | _evalexpo((BITS_IN_LONG-1)-sh);
1367   y[2] = x<<sh; for (i=3; i<ly; i++) y[i] = 0;
1368 }
1369 
1370 INLINE void
affiz(GEN x,GEN y)1371 affiz(GEN x, GEN y) { if (typ(y)==t_INT) affii(x,y); else affir(x,y); }
1372 INLINE void
affsz(long x,GEN y)1373 affsz(long x, GEN y) { if (typ(y)==t_INT) affsi(x,y); else affsr(x,y); }
1374 INLINE void
mpaff(GEN x,GEN y)1375 mpaff(GEN x, GEN y) { if (typ(x)==t_INT) affiz(x, y); else affrr(x,y); }
1376 
1377 /*******************************************************************/
1378 /*                                                                 */
1379 /*                    OPERATION + ASSIGNMENT                       */
1380 /*                                                                 */
1381 /*******************************************************************/
1382 
addiiz(GEN x,GEN y,GEN z)1383 INLINE void addiiz(GEN x, GEN y, GEN z)
1384 { pari_sp av = avma; affii(addii(x,y),z); set_avma(av); }
addirz(GEN x,GEN y,GEN z)1385 INLINE void addirz(GEN x, GEN y, GEN z)
1386 { pari_sp av = avma; affrr(addir(x,y),z); set_avma(av); }
addriz(GEN x,GEN y,GEN z)1387 INLINE void addriz(GEN x, GEN y, GEN z)
1388 { pari_sp av = avma; affrr(addri(x,y),z); set_avma(av); }
addrrz(GEN x,GEN y,GEN z)1389 INLINE void addrrz(GEN x, GEN y, GEN z)
1390 { pari_sp av = avma; affrr(addrr(x,y),z); set_avma(av); }
addsiz(long s,GEN y,GEN z)1391 INLINE void addsiz(long s, GEN y, GEN z)
1392 { pari_sp av = avma; affii(addsi(s,y),z); set_avma(av); }
addsrz(long s,GEN y,GEN z)1393 INLINE void addsrz(long s, GEN y, GEN z)
1394 { pari_sp av = avma; affrr(addsr(s,y),z); set_avma(av); }
addssz(long s,long y,GEN z)1395 INLINE void addssz(long s, long y, GEN z)
1396 { pari_sp av = avma; affii(addss(s,y),z); set_avma(av); }
1397 
diviiz(GEN x,GEN y,GEN z)1398 INLINE void diviiz(GEN x, GEN y, GEN z)
1399 { pari_sp av = avma; affii(divii(x,y),z); set_avma(av); }
divirz(GEN x,GEN y,GEN z)1400 INLINE void divirz(GEN x, GEN y, GEN z)
1401 { pari_sp av = avma; mpaff(divir(x,y),z); set_avma(av); }
divisz(GEN x,long y,GEN z)1402 INLINE void divisz(GEN x, long y, GEN z)
1403 { pari_sp av = avma; affii(divis(x,y),z); set_avma(av); }
divriz(GEN x,GEN y,GEN z)1404 INLINE void divriz(GEN x, GEN y, GEN z)
1405 { pari_sp av = avma; affrr(divri(x,y),z); set_avma(av); }
divrrz(GEN x,GEN y,GEN z)1406 INLINE void divrrz(GEN x, GEN y, GEN z)
1407 { pari_sp av = avma; affrr(divrr(x,y),z); set_avma(av); }
divrsz(GEN y,long s,GEN z)1408 INLINE void divrsz(GEN y, long s, GEN z)
1409 { pari_sp av = avma; affrr(divrs(y,s),z); set_avma(av); }
divsiz(long x,GEN y,GEN z)1410 INLINE void divsiz(long x, GEN y, GEN z)
1411 { long junk; affsi(sdivsi_rem(x,y,&junk), z); }
divsrz(long s,GEN y,GEN z)1412 INLINE void divsrz(long s, GEN y, GEN z)
1413 { pari_sp av = avma; mpaff(divsr(s,y),z); set_avma(av); }
divssz(long x,long y,GEN z)1414 INLINE void divssz(long x, long y, GEN z)
1415 { affsi(x/y, z); }
1416 
modisz(GEN y,long s,GEN z)1417 INLINE void modisz(GEN y, long s, GEN z)
1418 { affsi(smodis(y,s),z); }
modsiz(long s,GEN y,GEN z)1419 INLINE void modsiz(long s, GEN y, GEN z)
1420 { pari_sp av = avma; affii(modsi(s,y),z); set_avma(av); }
modssz(long s,long y,GEN z)1421 INLINE void modssz(long s, long y, GEN z)
1422 { affsi(smodss(s,y),z); }
1423 
mpaddz(GEN x,GEN y,GEN z)1424 INLINE void mpaddz(GEN x, GEN y, GEN z)
1425 { pari_sp av = avma; mpaff(mpadd(x,y),z); set_avma(av); }
mpsubz(GEN x,GEN y,GEN z)1426 INLINE void mpsubz(GEN x, GEN y, GEN z)
1427 { pari_sp av = avma; mpaff(mpsub(x,y),z); set_avma(av); }
mpmulz(GEN x,GEN y,GEN z)1428 INLINE void mpmulz(GEN x, GEN y, GEN z)
1429 { pari_sp av = avma; mpaff(mpmul(x,y),z); set_avma(av); }
1430 
muliiz(GEN x,GEN y,GEN z)1431 INLINE void muliiz(GEN x, GEN y, GEN z)
1432 { pari_sp av = avma; affii(mulii(x,y),z); set_avma(av); }
mulirz(GEN x,GEN y,GEN z)1433 INLINE void mulirz(GEN x, GEN y, GEN z)
1434 { pari_sp av = avma; mpaff(mulir(x,y),z); set_avma(av); }
mulriz(GEN x,GEN y,GEN z)1435 INLINE void mulriz(GEN x, GEN y, GEN z)
1436 { pari_sp av = avma; mpaff(mulri(x,y),z); set_avma(av); }
mulrrz(GEN x,GEN y,GEN z)1437 INLINE void mulrrz(GEN x, GEN y, GEN z)
1438 { pari_sp av = avma; affrr(mulrr(x,y),z); set_avma(av); }
mulsiz(long s,GEN y,GEN z)1439 INLINE void mulsiz(long s, GEN y, GEN z)
1440 { pari_sp av = avma; affii(mulsi(s,y),z); set_avma(av); }
mulsrz(long s,GEN y,GEN z)1441 INLINE void mulsrz(long s, GEN y, GEN z)
1442 { pari_sp av = avma; mpaff(mulsr(s,y),z); set_avma(av); }
mulssz(long s,long y,GEN z)1443 INLINE void mulssz(long s, long y, GEN z)
1444 { pari_sp av = avma; affii(mulss(s,y),z); set_avma(av); }
1445 
remiiz(GEN x,GEN y,GEN z)1446 INLINE void remiiz(GEN x, GEN y, GEN z)
1447 { pari_sp av = avma; affii(remii(x,y),z); set_avma(av); }
remisz(GEN y,long s,GEN z)1448 INLINE void remisz(GEN y, long s, GEN z)
1449 { pari_sp av = avma; affii(remis(y,s),z); set_avma(av); }
remsiz(long s,GEN y,GEN z)1450 INLINE void remsiz(long s, GEN y, GEN z)
1451 { pari_sp av = avma; affii(remsi(s,y),z); set_avma(av); }
remssz(long s,long y,GEN z)1452 INLINE void remssz(long s, long y, GEN z)
1453 { pari_sp av = avma; affii(remss(s,y),z); set_avma(av); }
1454 
subiiz(GEN x,GEN y,GEN z)1455 INLINE void subiiz(GEN x, GEN y, GEN z)
1456 { pari_sp av = avma; affii(subii(x,y),z); set_avma(av); }
subirz(GEN x,GEN y,GEN z)1457 INLINE void subirz(GEN x, GEN y, GEN z)
1458 { pari_sp av = avma; affrr(subir(x,y),z); set_avma(av); }
subisz(GEN y,long s,GEN z)1459 INLINE void subisz(GEN y, long s, GEN z)
1460 { pari_sp av = avma; affii(addsi(-s,y),z); set_avma(av); }
subriz(GEN x,GEN y,GEN z)1461 INLINE void subriz(GEN x, GEN y, GEN z)
1462 { pari_sp av = avma; affrr(subri(x,y),z); set_avma(av); }
subrrz(GEN x,GEN y,GEN z)1463 INLINE void subrrz(GEN x, GEN y, GEN z)
1464 { pari_sp av = avma; affrr(subrr(x,y),z); set_avma(av); }
subrsz(GEN y,long s,GEN z)1465 INLINE void subrsz(GEN y, long s, GEN z)
1466 { pari_sp av = avma; affrr(addsr(-s,y),z); set_avma(av); }
subsiz(long s,GEN y,GEN z)1467 INLINE void subsiz(long s, GEN y, GEN z)
1468 { pari_sp av = avma; affii(subsi(s,y),z); set_avma(av); }
subsrz(long s,GEN y,GEN z)1469 INLINE void subsrz(long s, GEN y, GEN z)
1470 { pari_sp av = avma; affrr(subsr(s,y),z); set_avma(av); }
subssz(long x,long y,GEN z)1471 INLINE void subssz(long x, long y, GEN z) { addssz(x,-y,z); }
1472 
1473 INLINE void
dvmdssz(long x,long y,GEN z,GEN t)1474 dvmdssz(long x, long y, GEN z, GEN t) {
1475   pari_sp av = avma;
1476   long r;
1477   affii(divss_rem(x,y, &r), z); set_avma(av); affsi(r,t);
1478 }
1479 INLINE void
dvmdsiz(long x,GEN y,GEN z,GEN t)1480 dvmdsiz(long x, GEN y, GEN z, GEN t) {
1481   pari_sp av = avma;
1482   long r;
1483   affii(divsi_rem(x,y, &r), z); set_avma(av); affsi(r,t);
1484 }
1485 INLINE void
dvmdisz(GEN x,long y,GEN z,GEN t)1486 dvmdisz(GEN x, long y, GEN z, GEN t) {
1487   pari_sp av = avma;
1488   long r;
1489   affii(divis_rem(x,y, &r),z); set_avma(av); affsi(r,t);
1490 }
1491 INLINE void
dvmdiiz(GEN x,GEN y,GEN z,GEN t)1492 dvmdiiz(GEN x, GEN y, GEN z, GEN t) {
1493   pari_sp av = avma;
1494   GEN r;
1495   affii(dvmdii(x,y,&r),z); affii(r,t); set_avma(av);
1496 }
1497