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