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