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