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