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 1677876463 : nbits2lg(long x) {
22 1677876463 : return (long)(((ulong)x+3*BITS_IN_LONG-1) >> TWOPOTBITS_IN_LONG);
23 : }
24 : INLINE long
25 30275830 : lg2prec(long x){ return (x-2) * BITS_IN_LONG; }
26 :
27 : INLINE long
28 >10009*10^7 : evallg(long x)
29 : {
30 >10009*10^7 : if (x & ~LGBITS) pari_err_OVERFLOW("lg()");
31 >10009*10^7 : return _evallg(x);
32 : }
33 : INLINE long
34 81660145 : evalvalp(long x)
35 : {
36 81660145 : long v = _evalvalp(x);
37 81660145 : if (v & ~VALPBITS) pari_err_OVERFLOW("valp()");
38 81660145 : return v;
39 : }
40 : INLINE long
41 25309350 : evalvalser(long x)
42 : {
43 25309350 : long v = _evalvalser(x);
44 25309350 : if (v & ~VALSERBITS) pari_err_OVERFLOW("valser()");
45 25309350 : return v;
46 : }
47 : INLINE long
48 13971335128 : evalexpo(long x)
49 : {
50 13971335128 : long v = _evalexpo(x);
51 13971335128 : if (v & ~EXPOBITS) pari_err_OVERFLOW("expo()");
52 13971335121 : return v;
53 : }
54 : INLINE long
55 78703542 : evalprecp(long x)
56 : {
57 78703542 : long v = _evalprecp(x);
58 78703542 : if (x & ~((1UL<<(BITS_IN_LONG-VALPnumBITS))-1)) pari_err_OVERFLOW("precp()");
59 78703542 : return v;
60 : }
61 :
62 : INLINE int
63 218513211 : varncmp(long x, long y)
64 : {
65 218513211 : if (varpriority[x] < varpriority[y]) return 1;
66 153653592 : if (varpriority[x] > varpriority[y]) return -1;
67 74173894 : return 0;
68 : }
69 : INLINE long
70 26180 : varnmin(long x, long y)
71 26180 : { return (varpriority[x] <= varpriority[y])? x: y; }
72 : INLINE long
73 203 : varnmax(long x, long y)
74 203 : { return (varpriority[x] >= varpriority[y])? x: y; }
75 :
76 : /* Inhibit some area GC-wise: declare it to be a non recursive
77 : * type, of length l. Thus gc_stack_update won't inspect the zone, just copy it.
78 : * For the following situation:
79 : * z = cgetg(t,a); av = avma; garbage(); ltop = avma;
80 : * for (i=1; i<HUGE; i++) gel(z,i) = blah();
81 : * stackdummy(av,ltop);
82 : * loses (av-ltop) words but saves a GC. */
83 : INLINE void
84 3743990719 : stackdummy(pari_sp av, pari_sp ltop) {
85 3743990719 : long l = ((GEN)av) - ((GEN)ltop);
86 3743990719 : if (l > 0) {
87 1248709976 : GEN z = (GEN)ltop;
88 1248709976 : 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 3743990719 : }
94 : INLINE void
95 106530559 : fixlg(GEN x, long ly) {
96 106530559 : long lx = lg(x), l = lx - ly;
97 106530559 : if (l > 0)
98 : { /* stackdummy(x+lx, x+ly) */
99 49844543 : GEN z = x + ly;
100 49844543 : z[0] = evaltyp(t_VECSMALL) | evallg(l);
101 49844543 : setlg(x, ly);
102 : #ifdef DEBUG
103 : { long i; for (i = 1; i < l; i++) z[i] = 0; }
104 : #endif
105 : }
106 106530559 : }
107 : /* update lg(z) before affrr(y, z) [ to cater for precision loss ]*/
108 : INLINE void
109 56699207 : 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 >13604*10^7 : set_avma(ulong av) { avma = av; }
120 :
121 : INLINE double
122 198969854 : gc_double(pari_sp av, double d) { set_avma(av); return d; }
123 : INLINE long
124 243061012 : gc_long(pari_sp av, long s) { set_avma(av); return s; }
125 : INLINE ulong
126 52329636 : gc_ulong(pari_sp av, ulong s) { set_avma(av); return s; }
127 : INLINE int
128 55261962 : gc_bool(pari_sp av, int s) { set_avma(av); return s; }
129 : INLINE int
130 2818782 : gc_int(pari_sp av, int s) { set_avma(av); return s; }
131 : INLINE GEN
132 7911980 : gc_NULL(pari_sp av) { set_avma(av); return NULL; }
133 : INLINE GEN
134 15930252397 : gc_const(pari_sp av, GEN x) { set_avma(av); return x; }
135 : #define retgc_const(av, x) \
136 : do { set_avma(av); return x; } while(0)
137 :
138 : INLINE GEN
139 155904 : gc_stoi(pari_sp av, long x) { set_avma(av); return stoi(x); }
140 : INLINE GEN
141 467783 : gc_utoi(pari_sp av, ulong x) { set_avma(av); return utoi(x); }
142 : INLINE GEN
143 1155325 : gc_utoipos(pari_sp av, ulong x) { set_avma(av); return utoipos(x); }
144 :
145 : INLINE GEN
146 98658584117 : new_chunk(size_t x) /* x is a number of longs */
147 : {
148 98658584117 : GEN z = ((GEN) avma) - x;
149 : CHECK_CTRLC
150 98658584117 : if (x > (avma-pari_mainstack->bot) / sizeof(long))
151 18 : new_chunk_resize(x);
152 98658584099 : set_avma((pari_sp)z);
153 : #ifdef MEMSTEP
154 : if (DEBUGMEM>1 && pari_mainstack->memused != DISABLE_MEMUSED) {
155 : long d = (long)pari_mainstack->memused - (long)z;
156 : if (labs(d) > 4*MEMSTEP)
157 : {
158 : pari_mainstack->memused = (pari_sp)z;
159 : err_printf("...%4.0lf Mbytes used\n",
160 : (pari_mainstack->top-pari_mainstack->memused)/1048576.);
161 : }
162 : }
163 : #endif
164 98658584096 : return z;
165 : }
166 :
167 : INLINE char *
168 46736502 : stack_malloc(size_t N)
169 : {
170 46736502 : long n = nchar2nlong(N);
171 46736502 : return (char*)new_chunk(n);
172 : }
173 :
174 : INLINE char *
175 55670147 : stack_malloc_align(size_t N, long k)
176 : {
177 55670147 : ulong d = ((ulong)avma) % k, e = ((ulong)N) % k;
178 55670147 : if (d) (void)new_chunk(d/sizeof(long));
179 55670147 : if (e) N += k-e;
180 55670147 : return (char*) new_chunk(nchar2nlong(N));
181 : }
182 :
183 : INLINE char *
184 110628 : stack_calloc(size_t N)
185 : {
186 110628 : char *p = stack_malloc(N);
187 110628 : memset(p, 0, N); return p;
188 : }
189 :
190 : INLINE char *
191 3325 : stack_calloc_align(size_t N, long k)
192 : {
193 3325 : ulong d = ((ulong)avma) % k, e = ((ulong)N) % k;
194 3325 : if (d) (void)new_chunk(d/sizeof(long));
195 3325 : if (e) N += k-e;
196 3325 : return stack_calloc(N);
197 : }
198 :
199 : INLINE GEN
200 857705 : newblock(size_t n)
201 857705 : { return newblock_t(n, 0); }
202 :
203 : /* cgetg(lg(x), typ(x)), set *lx. Implicit unsetisclone() */
204 : INLINE GEN
205 1490907748 : cgetg_copy(GEN x, long *plx) {
206 : GEN y;
207 1490907748 : *plx = lg(x); y = new_chunk((size_t)*plx);
208 1490907748 : y[0] = x[0] & (TYPBITS|LGBITS); return y;
209 : }
210 : INLINE GEN
211 430408 : cgetg_block(long x, long y)
212 : {
213 430408 : GEN z = newblock((size_t)x);
214 430408 : z[0] = CLONEBIT | evaltyp(y) | evallg(x);
215 430408 : return z;
216 : }
217 : INLINE GEN
218 28034801077 : cgetg(long x, long y)
219 : {
220 28034801077 : GEN z = new_chunk((size_t)x);
221 28034801077 : z[0] = evaltyp(y) | evallg(x);
222 28034801077 : return z;
223 : }
224 : INLINE GEN
225 28638383809 : cgeti(long x)
226 : {
227 28638383809 : GEN z = new_chunk((size_t)x);
228 28638383800 : z[0] = evaltyp(t_INT) | evallg(x);
229 28638383795 : return z;
230 : }
231 : INLINE GEN
232 16904123917 : cgetipos(long x)
233 : {
234 16904123917 : GEN z = cgeti(x);
235 16904123917 : z[1] = evalsigne(1) | evallgefint(x);
236 16904123917 : return z;
237 : }
238 : INLINE GEN
239 301591275 : cgetineg(long x)
240 : {
241 301591275 : GEN z = cgeti(x);
242 301591275 : z[1] = evalsigne(-1) | evallgefint(x);
243 301591275 : return z;
244 : }
245 : INLINE GEN
246 44721 : cgetr_block(long x)
247 : {
248 44721 : long l = nbits2lg(x);
249 44721 : GEN z = newblock((size_t)l);
250 44721 : z[0] = CLONEBIT | evaltyp(t_REAL) | evallg(l);
251 44721 : return z;
252 : }
253 : INLINE GEN
254 192135538 : cgetr(long x)
255 : {
256 192135538 : long l = nbits2lg(x);
257 192135538 : GEN z = new_chunk((size_t)l);
258 192135538 : z[0] = evaltyp(t_REAL) | evallg(l);
259 192135538 : return z;
260 : }
261 :
262 : /*******************************************************************/
263 : /* */
264 : /* COPY, NEGATION, ABSOLUTE VALUE */
265 : /* */
266 : /*******************************************************************/
267 : /* cannot do memcpy because sometimes x and y overlap */
268 : INLINE GEN
269 5279206609 : leafcopy(GEN x)
270 : {
271 5279206609 : long lx = lg(x);
272 5279206609 : GEN y = new_chunk(lx); /* can't use cgetg_copy, in case x,y overlap */
273 26997618600 : while (--lx > 0) y[lx] = x[lx];
274 5279206609 : y[0] = x[0] & (TYPBITS|LGBITS); return y;
275 : }
276 : INLINE GEN
277 9395663922 : icopy(GEN x)
278 : {
279 9395663922 : long i = lgefint(x), lx = i;
280 9395663922 : GEN y = new_chunk(lx); /* can't use cgeti, in case x,y overlap */
281 42302986554 : while (--i > 0) y[i] = x[i];
282 9395663915 : y[0] = evaltyp(t_INT) | evallg(lx);
283 9395663915 : return y;
284 : }
285 : INLINE GEN
286 120486428 : icopyspec(GEN x, long nx)
287 : {
288 120486428 : long i = nx+2, lx = i;
289 120486428 : GEN y = new_chunk(lx); /* can't use cgeti, in case x,y overlap */
290 3144588348 : x -= 2; while (--i >= 2) y[i] = x[i];
291 120486428 : y[1] = evalsigne(1) | evallgefint(lx);
292 120486428 : y[0] = evaltyp(t_INT) | evallg(lx);
293 120486428 : return y;
294 : }
295 1101916000 : INLINE GEN rcopy(GEN x) { return leafcopy(x); }
296 665 : INLINE GEN mpcopy(GEN x) { return leafcopy(x); }
297 :
298 : INLINE GEN
299 2118530823 : mpabs(GEN x) { GEN y = leafcopy(x); setabssign(y); return y; }
300 : INLINE GEN
301 13428274 : mpabs_shallow(GEN x) { return signe(x) < 0? mpabs(x): x; }
302 2076460290 : INLINE GEN absi(GEN x) { return mpabs(x); }
303 61597959 : INLINE GEN absi_shallow(GEN x) { return signe(x) < 0? negi(x): x; }
304 20510 : INLINE GEN absr(GEN x) { return mpabs(x); }
305 :
306 : INLINE GEN
307 952047017 : mpneg(GEN x) { GEN y = leafcopy(x); togglesign(y); return y; }
308 661628467 : INLINE GEN negi(GEN x) { return mpneg(x); }
309 3961190 : INLINE GEN negr(GEN x) { return mpneg(x); }
310 :
311 : /* negate in place */
312 : INLINE void
313 1996630297 : togglesign(GEN x) { if (x[1] & SIGNBITS) { x[1] ^= HIGHBIT; } }
314 : INLINE void
315 2187194506 : setabssign(GEN x) { x[1] &= ~HIGHBIT; }
316 : /* negate in place, except universal constants */
317 : INLINE void
318 157546307 : togglesign_safe(GEN *px)
319 : {
320 157546307 : switch(*px - gen_1) /* gen_1, gen_2, gen_m1, gen_m2 */
321 : {
322 3671096 : case 0: *px = gen_m1; break;
323 4 : case 3: *px = gen_m2; break;
324 846651 : case 6: *px = gen_1; break;
325 0 : case 9: *px = gen_2; break;
326 153028556 : default: togglesign(*px);
327 : }
328 157546307 : }
329 : /* setsigne(y, signe(x)) */
330 : INLINE void
331 0 : affectsign(GEN x, GEN y)
332 : {
333 0 : y[1] = (x[1] & SIGNBITS) | (y[1] & ~SIGNBITS);
334 0 : }
335 : /* copies sign in place, except for universal constants */
336 : INLINE void
337 10754891 : affectsign_safe(GEN x, GEN *py)
338 : {
339 10754891 : if (((*py)[1] ^ x[1]) & HIGHBIT) togglesign_safe(py);
340 10754891 : }
341 : /*******************************************************************/
342 : /* */
343 : /* GEN -> LONG, LONG -> GEN */
344 : /* */
345 : /*******************************************************************/
346 : /* assume x != 0, return -x as a t_INT */
347 : INLINE GEN
348 300668514 : utoineg(ulong x) { GEN y = cgetineg(3); y[2] = x; return y; }
349 : /* assume x != 0, return utoi(x) */
350 : INLINE GEN
351 14665164501 : utoipos(ulong x) { GEN y = cgetipos(3); y[2] = x; return y; }
352 : INLINE GEN
353 12431102919 : utoi(ulong x) { return x? utoipos(x): gen_0; }
354 : INLINE GEN
355 776330588 : stoi(long x)
356 : {
357 776330588 : if (!x) return gen_0;
358 553196016 : return x > 0? utoipos((ulong)x): utoineg((ulong)-x);
359 : }
360 :
361 : /* x 2^BIL + y */
362 : INLINE GEN
363 9226859190 : uutoi(ulong x, ulong y)
364 : {
365 : GEN z;
366 9226859190 : if (!x) return utoi(y);
367 884593438 : z = cgetipos(4);
368 884593438 : *int_W_lg(z, 1, 4) = x;
369 884593438 : *int_W_lg(z, 0, 4) = y; return z;
370 : }
371 : /* - (x 2^BIL + y) */
372 : INLINE GEN
373 365874 : uutoineg(ulong x, ulong y)
374 : {
375 : GEN z;
376 365874 : if (!x) return y? utoineg(y): gen_0;
377 16534 : z = cgetineg(4);
378 16534 : *int_W_lg(z, 1, 4) = x;
379 16534 : *int_W_lg(z, 0, 4) = y; return z;
380 : }
381 :
382 : INLINE long
383 477849879 : itos(GEN x)
384 : {
385 477849879 : long s = signe(x);
386 : long u;
387 :
388 477849879 : if (!s) return 0;
389 445948637 : u = x[2];
390 445948637 : if (lgefint(x) > 3 || u < 0)
391 30 : pari_err_OVERFLOW("t_INT-->long assignment");
392 445948607 : return (s>0) ? u : -u;
393 : }
394 : /* as itos, but return 0 if too large. Cf is_bigint */
395 : INLINE long
396 44914236 : itos_or_0(GEN x) {
397 : long n;
398 44914236 : if (lgefint(x) != 3 || (n = x[2]) & HIGHBIT) return 0;
399 38506001 : return signe(x) > 0? n: -n;
400 : }
401 : INLINE ulong
402 180921936 : itou(GEN x)
403 : {
404 180921936 : switch(lgefint(x)) {
405 13844071 : case 2: return 0;
406 167077865 : case 3: return x[2];
407 0 : default:
408 0 : pari_err_OVERFLOW("t_INT-->ulong assignment");
409 : return 0; /* LCOV_EXCL_LINE */
410 : }
411 : }
412 :
413 : /* as itou, but return 0 if too large. Cf is_bigint */
414 : INLINE ulong
415 5015383 : itou_or_0(GEN x) {
416 5015383 : if (lgefint(x) != 3) return 0;
417 4994304 : return (ulong)x[2];
418 : }
419 :
420 : INLINE ulong
421 6349300 : umuluu_or_0(ulong x, ulong y)
422 : {
423 : ulong z;
424 : LOCAL_HIREMAINDER;
425 6349300 : z = mulll(x, y);
426 6349300 : return hiremainder? 0: z;
427 : }
428 : /* return x*y if <= n, else 0. Beware overflow */
429 : INLINE ulong
430 7618608 : umuluu_le(ulong x, ulong y, ulong n)
431 : {
432 : ulong z;
433 : LOCAL_HIREMAINDER;
434 7618608 : z = mulll(x, y);
435 7618608 : return (hiremainder || z > n)? 0: z;
436 : }
437 :
438 : INLINE GEN
439 501081636 : real_0_bit(long bitprec) { GEN x=cgetg(2, t_REAL); x[1]=evalexpo(bitprec); return x; }
440 : INLINE GEN
441 1160900 : real_0(long prec) { return real_0_bit(-prec); }
442 : INLINE GEN
443 4699130 : real_1_bit(long bit) { return real_1(nbits2prec(bit)); }
444 : INLINE GEN
445 133392530 : real_1(long prec) {
446 133392530 : long i, l = nbits2lg(prec);
447 133392530 : GEN x = cgetg(l, t_REAL);
448 133392530 : x[1] = evalsigne(1) | _evalexpo(0);
449 665695529 : x[2] = (long)HIGHBIT; for (i=3; i<l; i++) x[i] = 0;
450 133392530 : return x;
451 : }
452 : INLINE GEN
453 455 : real_m1(long prec) {
454 455 : long i, l = nbits2lg(prec);
455 455 : GEN x = cgetg(l, t_REAL);
456 455 : x[1] = evalsigne(-1) | _evalexpo(0);
457 1761 : x[2] = (long)HIGHBIT; for (i=3; i<l; i++) x[i] = 0;
458 455 : return x;
459 : }
460 :
461 : /* 2.^n */
462 : INLINE GEN
463 1070088 : real2n(long n, long prec) { GEN z = real_1(prec); setexpo(z, n); return z; }
464 : INLINE GEN
465 126 : real_m2n(long n, long prec) { GEN z = real_m1(prec); setexpo(z, n); return z; }
466 :
467 : INLINE GEN
468 509053161 : stor_lg(long s, long l) { GEN z = cgetg(l, t_REAL); affsr(s,z); return z; }
469 : INLINE GEN
470 2140784 : stor(long s, long prec) { return stor_lg(s, nbits2lg(prec)); }
471 :
472 : INLINE GEN
473 13344627 : utor_lg(ulong u, long l){ GEN z = cgetg(l, t_REAL); affur(u,z); return z; }
474 : INLINE GEN
475 13344627 : utor(ulong u, long prec){ return utor_lg(u, nbits2lg(prec)); }
476 :
477 : INLINE GEN
478 739838246 : itor_lg(GEN x, long l) { GEN z = cgetg(l, t_REAL); affir(x,z); return z; }
479 : INLINE GEN
480 451882210 : itor(GEN x, long prec) { return itor_lg(x, nbits2lg(prec)); }
481 :
482 : INLINE GEN
483 308825017 : rtor_lg(GEN x, long l) { GEN z = cgetg(l, t_REAL); affrr(x,z); return z; }
484 : INLINE GEN
485 308820663 : rtor(GEN x, long prec) { return rtor_lg(x, nbits2lg(prec)); }
486 :
487 : INLINE ulong
488 21015730 : int_bit(GEN x, long n)
489 : {
490 21015730 : long r, q = dvmdsBIL(n, &r);
491 21015730 : return q < lgefint(x)-2?((ulong)*int_W(x,q) >> r) & 1UL:0;
492 : }
493 :
494 : /*******************************************************************/
495 : /* */
496 : /* COMPARISON */
497 : /* */
498 : /*******************************************************************/
499 : INLINE int
500 1416201 : cmpss(long a, long b)
501 1416201 : { return a>b? 1: (a<b? -1: 0); }
502 :
503 : INLINE int
504 318396336 : cmpuu(ulong a, ulong b)
505 318396336 : { return a>b? 1: (a<b? -1: 0); }
506 :
507 : INLINE int
508 1730010 : cmpir(GEN x, GEN y)
509 : {
510 : pari_sp av;
511 : GEN z;
512 :
513 1730010 : if (!signe(x)) return -signe(y);
514 421424 : if (!signe(y))
515 : {
516 6182 : if (expo(y) >= expi(x)) return 0;
517 6154 : return signe(x);
518 : }
519 415242 : av=avma; z = itor(x, realprec(y)); set_avma(av);
520 415242 : return cmprr(z,y); /* cmprr does no memory adjustment */
521 : }
522 : INLINE int
523 303856 : cmpri(GEN x, GEN y) { return -cmpir(y,x); }
524 : INLINE int
525 796228 : cmpsr(long x, GEN y)
526 : {
527 : pari_sp av;
528 : GEN z;
529 :
530 796228 : if (!x) return -signe(y);
531 796228 : av=avma; z = stor(x, LOWDEFAULTPREC); set_avma(av);
532 796228 : return cmprr(z,y);
533 : }
534 : INLINE int
535 40996 : cmprs(GEN x, long y) { return -cmpsr(y,x); }
536 : /* compare x and y */
537 : INLINE int
538 13000195 : cmpui(ulong x, GEN y)
539 : {
540 : ulong p;
541 13000195 : if (!x) return -signe(y);
542 13000195 : if (signe(y) <= 0) return 1;
543 12872025 : if (lgefint(y) > 3) return -1;
544 12597467 : p = y[2]; if (p == x) return 0;
545 12493585 : return p < x ? 1 : -1;
546 : }
547 : INLINE int
548 13000195 : cmpiu(GEN x, ulong y) { return -cmpui(y,x); }
549 : /* compare x and |y| */
550 : INLINE int
551 33534405 : abscmpui(ulong x, GEN y)
552 : {
553 33534405 : long l = lgefint(y);
554 : ulong p;
555 :
556 33534405 : if (!x) return (l > 2)? -1: 0;
557 33534384 : if (l == 2) return 1;
558 33355851 : if (l > 3) return -1;
559 33334239 : p = y[2]; if (p == x) return 0;
560 32402665 : return p < x ? 1 : -1;
561 : }
562 : INLINE int
563 33534405 : abscmpiu(GEN x, ulong y) { return -abscmpui(y,x); }
564 : INLINE int
565 6351305 : cmpsi(long x, GEN y)
566 : {
567 : ulong p;
568 :
569 6351305 : if (!x) return -signe(y);
570 :
571 6347092 : if (x > 0)
572 : {
573 6346056 : if (signe(y)<=0) return 1;
574 6335115 : if (lgefint(y)>3) return -1;
575 6329213 : p = y[2]; if (p == (ulong)x) return 0;
576 6259056 : return p < (ulong)x ? 1 : -1;
577 : }
578 :
579 1036 : if (signe(y)>=0) return -1;
580 119 : if (lgefint(y)>3) return 1;
581 119 : p = y[2]; if (p == (ulong)-x) return 0;
582 14 : return p < (ulong)(-x) ? -1 : 1;
583 : }
584 : INLINE int
585 6119306 : cmpis(GEN x, long y) { return -cmpsi(y,x); }
586 : INLINE int
587 2164525 : mpcmp(GEN x, GEN y)
588 : {
589 2164525 : if (typ(x)==t_INT)
590 70406 : return (typ(y)==t_INT) ? cmpii(x,y) : cmpir(x,y);
591 2094119 : return (typ(y)==t_INT) ? -cmpir(y,x) : cmprr(x,y);
592 : }
593 :
594 : /* x == y ? */
595 : INLINE int
596 3042339 : equalui(ulong x, GEN y)
597 : {
598 3042339 : if (!x) return !signe(y);
599 3041646 : if (signe(y) <= 0 || lgefint(y) != 3) return 0;
600 3030083 : return ((ulong)y[2] == (ulong)x);
601 : }
602 : /* x == y ? */
603 : INLINE int
604 1149928 : equalsi(long x, GEN y)
605 : {
606 1149928 : if (!x) return !signe(y);
607 1149928 : if (x > 0)
608 : {
609 1142592 : if (signe(y) <= 0 || lgefint(y) != 3) return 0;
610 1064967 : return ((ulong)y[2] == (ulong)x);
611 : }
612 7336 : if (signe(y) >= 0 || lgefint(y) != 3) return 0;
613 4885 : return ((ulong)y[2] == (ulong)-x);
614 : }
615 : /* x == |y| ? */
616 : INLINE int
617 47655942 : absequalui(ulong x, GEN y)
618 : {
619 47655942 : if (!x) return !signe(y);
620 47655942 : return (lgefint(y) == 3 && (ulong)y[2] == x);
621 : }
622 : INLINE int
623 45906172 : absequaliu(GEN x, ulong y) { return absequalui(y,x); }
624 : INLINE int
625 1149746 : equalis(GEN x, long y) { return equalsi(y,x); }
626 : INLINE int
627 3042339 : equaliu(GEN x, ulong y) { return equalui(y,x); }
628 :
629 : /* assume x != 0, is |x| == 2^n ? */
630 : INLINE int
631 1580223 : absrnz_equal2n(GEN x) {
632 1580223 : if ((ulong)x[2]==HIGHBIT)
633 : {
634 206711 : long i, lx = lg(x);
635 615982 : for (i = 3; i < lx; i++)
636 424920 : if (x[i]) return 0;
637 191062 : return 1;
638 : }
639 1373512 : return 0;
640 : }
641 : /* assume x != 0, is |x| == 1 ? */
642 : INLINE int
643 4559224 : absrnz_equal1(GEN x) { return !expo(x) && absrnz_equal2n(x); }
644 :
645 : INLINE long
646 7491800963 : maxss(long x, long y) { return x>y?x:y; }
647 : INLINE long
648 2044959562 : minss(long x, long y) { return x<y?x:y; }
649 : INLINE long
650 57568285 : minuu(ulong x, ulong y) { return x<y?x:y; }
651 : INLINE long
652 6603110 : maxuu(ulong x, ulong y) { return x>y?x:y; }
653 : INLINE double
654 3218088 : maxdd(double x, double y) { return x>y?x:y; }
655 : INLINE double
656 265750 : mindd(double x, double y) { return x<y?x:y; }
657 :
658 : /*******************************************************************/
659 : /* */
660 : /* ADD / SUB */
661 : /* */
662 : /*******************************************************************/
663 : INLINE GEN
664 25081 : subuu(ulong x, ulong y)
665 : {
666 : ulong z;
667 : LOCAL_OVERFLOW;
668 25081 : z = subll(x, y);
669 25081 : return overflow? utoineg(-z): utoi(z);
670 : }
671 : INLINE GEN
672 3542666893 : adduu(ulong x, ulong y) { ulong t = x+y; return uutoi((t < x), t); }
673 :
674 : INLINE GEN
675 25081 : addss(long x, long y)
676 : {
677 25081 : if (!x) return stoi(y);
678 25081 : if (!y) return stoi(x);
679 25081 : if (x > 0) return y > 0? adduu(x,y): subuu(x, -y);
680 :
681 25081 : if (y > 0) return subuu(y, -x);
682 : else { /* - adduu(-x, -y) */
683 0 : ulong t = (-x)+(-y); return uutoineg((t < (ulong)(-x)), t);
684 : }
685 : }
686 25081 : INLINE GEN subss(long x, long y) { return addss(-y,x); }
687 :
688 : INLINE GEN
689 8008806639 : subii(GEN x, GEN y)
690 : {
691 8008806639 : if (x==y) return gen_0; /* frequent with x = y = gen_0 */
692 6642955318 : return addii_sign(x, signe(x), y, -signe(y));
693 : }
694 : INLINE GEN
695 13115454816 : addii(GEN x, GEN y) { return addii_sign(x, signe(x), y, signe(y)); }
696 : INLINE GEN
697 3032586280 : addrr(GEN x, GEN y) { return addrr_sign(x, signe(x), y, signe(y)); }
698 : INLINE GEN
699 1062988844 : subrr(GEN x, GEN y) { return addrr_sign(x, signe(x), y, -signe(y)); }
700 : INLINE GEN
701 490340776 : addir(GEN x, GEN y) { return addir_sign(x, signe(x), y, signe(y)); }
702 : INLINE GEN
703 3030801 : subir(GEN x, GEN y) { return addir_sign(x, signe(x), y, -signe(y)); }
704 : INLINE GEN
705 1538272 : subri(GEN x, GEN y) { return addir_sign(y, -signe(y), x, signe(x)); }
706 : INLINE GEN
707 320989966 : addsi(long x, GEN y) { return addsi_sign(x, y, signe(y)); }
708 : INLINE GEN
709 107891160 : addui(ulong x, GEN y) { return addui_sign(x, y, signe(y)); }
710 : INLINE GEN
711 6244790 : subsi(long x, GEN y) { return addsi_sign(x, y, -signe(y)); }
712 : INLINE GEN
713 138450476 : subui(ulong x, GEN y) { return addui_sign(x, y, -signe(y)); }
714 :
715 : /*******************************************************************/
716 : /* */
717 : /* MOD, REM, DIV */
718 : /* */
719 : /*******************************************************************/
720 113535923 : INLINE ulong mod2BIL(GEN x) { return *int_LSW(x); }
721 0 : INLINE long mod64(GEN x) { return mod2BIL(x) & 63; }
722 259 : INLINE long mod32(GEN x) { return mod2BIL(x) & 31; }
723 308662 : INLINE long mod16(GEN x) { return mod2BIL(x) & 15; }
724 15673609 : INLINE long mod8(GEN x) { return mod2BIL(x) & 7; }
725 7594956 : INLINE long mod4(GEN x) { return mod2BIL(x) & 3; }
726 65902721 : INLINE long mod2(GEN x) { return mod2BIL(x) & 1; }
727 : INLINE int
728 118265643 : mpodd(GEN x) { return signe(x) && mod2(x); }
729 : /* x mod 2^n, n < BITS_IN_LONG */
730 : INLINE ulong
731 58324899 : umodi2n(GEN x, long n)
732 : {
733 58324899 : long s = signe(x);
734 58324899 : const ulong _2n = 1UL << n;
735 : ulong m;
736 58324899 : if (!s) return 0;
737 56595055 : m = *int_LSW(x) & (_2n - 1);
738 56595055 : if (s < 0 && m) m = _2n - m;
739 56595055 : return m;
740 : }
741 0 : INLINE ulong Mod64(GEN x){ return umodi2n(x,6); }
742 311157 : INLINE ulong Mod32(GEN x){ return umodi2n(x,5); }
743 389809 : INLINE ulong Mod16(GEN x){ return umodi2n(x,4); }
744 3309346 : INLINE ulong Mod8(GEN x) { return umodi2n(x,3); }
745 50718975 : INLINE ulong Mod4(GEN x) { return umodi2n(x,2); }
746 3595514 : INLINE ulong Mod2(GEN x) { return umodi2n(x,1); }
747 :
748 : INLINE GEN
749 54901637 : truedivii(GEN a,GEN b) { return truedvmdii(a,b,NULL); }
750 : INLINE GEN
751 268438 : truedivis(GEN a, long b) { return truedvmdis(a,b,NULL); }
752 : INLINE GEN
753 6202318 : truedivsi(long a, GEN b) { return truedvmdsi(a,b,NULL); }
754 :
755 : INLINE GEN
756 16299032 : divii(GEN a, GEN b) { return dvmdii(a,b,NULL); }
757 : INLINE GEN
758 3104090886 : remii(GEN a, GEN b) { return dvmdii(a,b,ONLY_REM); }
759 :
760 : INLINE GEN
761 0 : divss(long x, long y) { return stoi(x / y); }
762 : INLINE GEN
763 0 : modss(long x, long y) { return utoi(smodss(x, y)); }
764 : INLINE GEN
765 0 : remss(long x, long y) { return stoi(x % y); }
766 : INLINE long
767 13998082 : smodss(long x, long y)
768 : {
769 13998082 : long r = x%y;
770 13998082 : return (r >= 0)? r: labs(y) + r;
771 : }
772 : INLINE ulong
773 734405677 : umodsu(long x, ulong y)
774 : {
775 734405677 : return x>=0 ? x%y: Fl_neg((-x)%y, y);
776 : }
777 :
778 : INLINE long
779 0 : sdivss_rem(long x, long y, long *r)
780 : {
781 : long q;
782 : LOCAL_HIREMAINDER;
783 0 : if (!y) pari_err_INV("sdivss_rem",gen_0);
784 0 : hiremainder = 0; q = divll((ulong)labs(x),(ulong)labs(y));
785 0 : if (x < 0) { hiremainder = -((long)hiremainder); q = -q; }
786 0 : if (y < 0) q = -q;
787 0 : *r = hiremainder; return q;
788 : }
789 : INLINE GEN
790 0 : divss_rem(long x, long y, long *r) { return stoi(sdivss_rem(x,y,r)); }
791 : INLINE ulong
792 2418435 : udivuu_rem(ulong x, ulong y, ulong *r)
793 : {
794 2418435 : if (!y) pari_err_INV("udivuu_rem",gen_0);
795 2418435 : *r = x % y; return x / y;
796 : }
797 : INLINE ulong
798 3916770 : ceildivuu(ulong a, ulong b)
799 : {
800 : ulong c;
801 3916770 : if (!a) return 0;
802 3732768 : c = a / b; return (a % b)? c+1: c;
803 : }
804 :
805 : INLINE ulong
806 128885 : uabsdivui_rem(ulong x, GEN y, ulong *r)
807 : {
808 128885 : long q, s = signe(y);
809 : LOCAL_HIREMAINDER;
810 :
811 128885 : if (!s) pari_err_INV("uabsdivui_rem",gen_0);
812 128885 : if (!x || lgefint(y)>3) { *r = x; return 0; }
813 128178 : hiremainder=0; q = (long)divll(x, (ulong)y[2]);
814 128178 : if (s < 0) q = -q;
815 128178 : *r = hiremainder; return q;
816 : }
817 :
818 : /* assume d != 0 and |n| / d can be represented as an ulong.
819 : * Return |n|/d, set *r = |n| % d */
820 : INLINE ulong
821 26737503 : uabsdiviu_rem(GEN n, ulong d, ulong *r)
822 : {
823 26737503 : switch(lgefint(n))
824 : {
825 0 : case 2: *r = 0; return 0;
826 26737503 : case 3:
827 : {
828 26737503 : ulong nn = n[2];
829 26737503 : *r = nn % d; return nn / d;
830 : }
831 0 : default: /* 4 */
832 : {
833 : ulong n1, n0, q;
834 : LOCAL_HIREMAINDER;
835 0 : n0 = *int_W(n,0);
836 0 : n1 = *int_W(n,1);
837 0 : hiremainder = n1;
838 0 : q = divll(n0, d);
839 0 : *r = hiremainder; return q;
840 : }
841 : }
842 : }
843 :
844 : INLINE long
845 51446757 : sdivsi_rem(long x, GEN y, long *r)
846 : {
847 51446757 : long q, s = signe(y);
848 : LOCAL_HIREMAINDER;
849 :
850 51446757 : if (!s) pari_err_INV("sdivsi_rem",gen_0);
851 51446757 : if (!x || lgefint(y)>3 || ((long)y[2]) < 0) { *r = x; return 0; }
852 49489317 : hiremainder=0; q = (long)divll(labs(x), (ulong)y[2]);
853 49489317 : if (x < 0) { hiremainder = -((long)hiremainder); q = -q; }
854 49489317 : if (s < 0) q = -q;
855 49489317 : *r = hiremainder; return q;
856 : }
857 : INLINE GEN
858 0 : divsi_rem(long s, GEN y, long *r) { return stoi(sdivsi_rem(s,y,r)); }
859 :
860 : INLINE long
861 102165 : sdivsi(long x, GEN y)
862 : {
863 102165 : long q, s = signe(y);
864 :
865 102165 : if (!s) pari_err_INV("sdivsi",gen_0);
866 102165 : if (!x || lgefint(y)>3 || ((long)y[2]) < 0) return 0;
867 102060 : q = labs(x) / y[2];
868 102060 : if (x < 0) q = -q;
869 102060 : if (s < 0) q = -q;
870 102060 : return q;
871 : }
872 :
873 : INLINE GEN
874 0 : dvmdss(long x, long y, GEN *z)
875 : {
876 : long r;
877 0 : GEN q = divss_rem(x,y, &r);
878 0 : *z = stoi(r); return q;
879 : }
880 : INLINE long
881 7387226371 : dvmdsBIL(long n, long *r) { *r = remsBIL(n); return divsBIL(n); }
882 : INLINE ulong
883 170507039 : dvmduBIL(ulong n, ulong *r) { *r = remsBIL(n); return divsBIL(n); }
884 : INLINE GEN
885 0 : dvmdsi(long x, GEN y, GEN *z)
886 : {
887 : long r;
888 0 : GEN q = divsi_rem(x,y, &r);
889 0 : *z = stoi(r); return q;
890 : }
891 : INLINE GEN
892 0 : dvmdis(GEN x, long y, GEN *z)
893 : {
894 : long r;
895 0 : GEN q = divis_rem(x,y, &r);
896 0 : *z = stoi(r); return q;
897 : }
898 :
899 : INLINE long
900 21176116 : smodis(GEN x, long y)
901 : {
902 21176116 : pari_sp av = avma;
903 21176116 : long r; (void)divis_rem(x,y, &r);
904 21176116 : return gc_long(av, (r >= 0)? r: labs(y) + r);
905 : }
906 : INLINE GEN
907 19638526 : modis(GEN x, long y) { return stoi(smodis(x,y)); }
908 : INLINE GEN
909 45242143 : modsi(long x, GEN y) {
910 45242143 : long r; (void)sdivsi_rem(x, y, &r);
911 45242143 : return (r >= 0)? stoi(r): addsi_sign(r, y, 1);
912 : }
913 :
914 : INLINE ulong
915 710986 : umodui(ulong x, GEN y)
916 : {
917 710986 : if (!signe(y)) pari_err_INV("umodui",gen_0);
918 710986 : if (!x || lgefint(y) > 3) return x;
919 414849 : return x % (ulong)y[2];
920 : }
921 :
922 : INLINE ulong
923 10312796 : ugcdiu(GEN x, ulong y) { return ugcd(umodiu(x,y), y); }
924 : INLINE ulong
925 2737 : ugcdui(ulong y, GEN x) { return ugcd(umodiu(x,y), y); }
926 :
927 : INLINE GEN
928 0 : remsi(long x, GEN y)
929 0 : { long r; (void)sdivsi_rem(x,y, &r); return stoi(r); }
930 : INLINE GEN
931 0 : remis(GEN x, long y)
932 : {
933 0 : pari_sp av = avma;
934 : long r;
935 0 : (void)divis_rem(x,y, &r); return gc_stoi(av, r);
936 : }
937 :
938 : INLINE GEN
939 0 : rdivis(GEN x, long y, long prec)
940 : {
941 0 : GEN z = cgetr(prec);
942 0 : pari_sp av = avma;
943 0 : affrr(divrs(itor(x,prec), y),z);
944 0 : set_avma(av); return z;
945 : }
946 : INLINE GEN
947 0 : rdivsi(long x, GEN y, long prec)
948 : {
949 0 : GEN z = cgetr(prec);
950 0 : pari_sp av = avma;
951 0 : affrr(divsr(x, itor(y,prec)), z);
952 0 : set_avma(av); return z;
953 : }
954 : INLINE GEN
955 839647 : rdivss(long x, long y, long prec)
956 : {
957 839647 : GEN z = cgetr(prec);
958 839647 : pari_sp av = avma;
959 839647 : affrr(divrs(stor(x, prec), y), z);
960 839647 : set_avma(av); return z;
961 : }
962 :
963 : INLINE void
964 13173936 : rdiviiz(GEN x, GEN y, GEN z)
965 : {
966 13173936 : long lz = lg(z), lx = lgefint(x), ly = lgefint(y);
967 13173936 : const long EXTRALG64 = EXTRAPREC64 / BITS_IN_LONG;
968 13173936 : if (lx == 2) { affur(0, z); return; }
969 : /* x != 0 */
970 13173936 : if (ly == 3)
971 : {
972 2256176 : affir(x, z); if (signe(y) < 0) togglesign(z);
973 2256176 : affrr(divru(z, y[2]), z);
974 : }
975 10917760 : else if (lx > lz + EXTRALG64 || ly > lz + EXTRALG64)
976 : {
977 5814857 : affir(x,z); affrr(divri(z, y), z);
978 : }
979 : else
980 : { /* bit size of shifti(x,b) < 2*realprec(z) + 64*/
981 5102903 : long b = lg2prec(lz) + expi(y) - expi(x) + 1;
982 5102903 : GEN q = divii(b > 0? shifti(x, b): x, y);
983 5102903 : affir(q, z); if (b > 0) shiftr_inplace(z, -b);
984 : }
985 13173936 : set_avma((ulong)z);
986 : }
987 : INLINE GEN
988 13134309 : rdivii_lg(GEN x, GEN y, long l)
989 13134309 : { GEN z = cgetg(l, t_REAL); rdiviiz(x, y, z); return z; }
990 : INLINE GEN
991 5744744 : rdivii(GEN x, GEN y, long prec) { return rdivii_lg(x, y, nbits2lg(prec)); }
992 : INLINE GEN
993 7389439 : fractor_lg(GEN x, long l) { return rdivii_lg(gel(x,1), gel(x,2), l); }
994 : INLINE GEN
995 7387948 : fractor(GEN x, long prec) { return fractor_lg(x, nbits2lg(prec)); }
996 :
997 : INLINE int
998 17072403 : dvdii(GEN x, GEN y)
999 : {
1000 17072403 : pari_sp av = avma;
1001 : GEN r;
1002 17072403 : if (!signe(x)) return 1;
1003 15635614 : if (!signe(y)) return 0;
1004 15635614 : r = remii(x,y);
1005 15635614 : return gc_bool(av, r == gen_0);
1006 : }
1007 : INLINE int
1008 371 : dvdsi(long x, GEN y)
1009 : {
1010 371 : if (x == 0) return 1;
1011 266 : if (!signe(y)) return 0;
1012 266 : if (lgefint(y) != 3) return 0;
1013 259 : return x % y[2] == 0;
1014 : }
1015 : INLINE int
1016 167545 : dvdui(ulong x, GEN y)
1017 : {
1018 167545 : if (x == 0) return 1;
1019 167545 : if (!signe(y)) return 0;
1020 167545 : if (lgefint(y) != 3) return 0;
1021 156924 : return x % y[2] == 0;
1022 : }
1023 : INLINE int
1024 33765 : dvdis(GEN x, long y)
1025 33765 : { return y? smodis(x, y) == 0: signe(x) == 0; }
1026 : INLINE int
1027 576287 : dvdiu(GEN x, ulong y)
1028 576287 : { return y? umodiu(x, y) == 0: signe(x) == 0; }
1029 :
1030 : INLINE int
1031 0 : dvdisz(GEN x, long y, GEN z)
1032 : {
1033 0 : const pari_sp av = avma;
1034 : long r;
1035 0 : GEN p1 = divis_rem(x,y, &r);
1036 0 : set_avma(av); if (r) return 0;
1037 0 : affii(p1,z); return 1;
1038 : }
1039 : INLINE int
1040 0 : dvdiuz(GEN x, ulong y, GEN z)
1041 : {
1042 0 : const pari_sp av = avma;
1043 : ulong r;
1044 0 : GEN p1 = absdiviu_rem(x,y, &r);
1045 0 : set_avma(av); if (r) return 0;
1046 0 : affii(p1,z); return 1;
1047 : }
1048 : INLINE int
1049 0 : dvdiiz(GEN x, GEN y, GEN z)
1050 : {
1051 0 : const pari_sp av=avma;
1052 0 : GEN p2, p1 = dvmdii(x,y,&p2);
1053 0 : if (signe(p2)) return gc_bool(av,0);
1054 0 : affii(p1,z); return gc_bool(av,1);
1055 : }
1056 :
1057 : INLINE ulong
1058 79167995 : remlll_pre(ulong u2, ulong u1, ulong u0, ulong n, ulong ninv)
1059 : {
1060 79167995 : u1 = remll_pre(u2, u1, n, ninv);
1061 79167995 : return remll_pre(u1, u0, n, ninv);
1062 : }
1063 :
1064 : INLINE ulong
1065 2568103943 : Fl_sqr_pre(ulong a, ulong p, ulong pi)
1066 : {
1067 : ulong x;
1068 : LOCAL_HIREMAINDER;
1069 2568103943 : x = mulll(a,a);
1070 2568103943 : return remll_pre(hiremainder, x, p, pi);
1071 : }
1072 :
1073 : INLINE ulong
1074 4448714363 : Fl_mul_pre(ulong a, ulong b, ulong p, ulong pi)
1075 : {
1076 : ulong x;
1077 : LOCAL_HIREMAINDER;
1078 4448714363 : x = mulll(a,b);
1079 4448714363 : return remll_pre(hiremainder, x, p, pi);
1080 : }
1081 :
1082 : INLINE ulong
1083 8205565349 : Fl_addmul_pre(ulong y0, ulong x0, ulong x1, ulong p, ulong pi)
1084 : {
1085 : ulong l0, h0;
1086 : LOCAL_HIREMAINDER;
1087 8205565349 : hiremainder = y0;
1088 8205565349 : l0 = addmul(x0, x1); h0 = hiremainder;
1089 8205565349 : return remll_pre(h0, l0, p, pi);
1090 : }
1091 :
1092 : INLINE ulong
1093 74211424 : Fl_addmulmul_pre(ulong x0, ulong y0, ulong x1, ulong y1, ulong p, ulong pi)
1094 : {
1095 : ulong l0, l1, h0, h1;
1096 : LOCAL_OVERFLOW;
1097 : LOCAL_HIREMAINDER;
1098 74211424 : l0 = mulll(x0, y0); h0 = hiremainder;
1099 74211424 : l1 = mulll(x1, y1); h1 = hiremainder;
1100 74211424 : l0 = addll(l0, l1); h0 = addllx(h0, h1);
1101 74211424 : return overflow ? remlll_pre(1, h0, l0, p, pi): remll_pre(h0, l0, p, pi);
1102 : }
1103 :
1104 : INLINE ulong
1105 242716 : Fl_ellj_pre(ulong a4, ulong a6, ulong p, ulong pi)
1106 : {
1107 : /* a43 = 4 a4^3 */
1108 242716 : ulong a43 = Fl_double(Fl_double(
1109 : Fl_mul_pre(a4, Fl_sqr_pre(a4, p, pi), p, pi), p), p);
1110 : /* a62 = 27 a6^2 */
1111 242716 : ulong a62 = Fl_mul_pre(Fl_sqr_pre(a6, p, pi), 27 % p, p, pi);
1112 242716 : ulong z1 = Fl_mul_pre(a43, 1728 % p, p, pi);
1113 242716 : ulong z2 = Fl_add(a43, a62, p);
1114 242716 : return Fl_div(z1, z2, p);
1115 : }
1116 :
1117 : /*******************************************************************/
1118 : /* */
1119 : /* MP (INT OR REAL) */
1120 : /* */
1121 : /*******************************************************************/
1122 : INLINE GEN
1123 49 : mptrunc(GEN x) { return typ(x)==t_INT? icopy(x): truncr(x); }
1124 : INLINE GEN
1125 0 : mpfloor(GEN x) { return typ(x)==t_INT? icopy(x): floorr(x); }
1126 : INLINE GEN
1127 0 : mpceil(GEN x) { return typ(x)==t_INT? icopy(x): ceilr(x); }
1128 : INLINE GEN
1129 1216171 : mpround(GEN x) { return typ(x) == t_INT? icopy(x): roundr(x); }
1130 :
1131 : INLINE long
1132 37845523 : mpexpo(GEN x) { return typ(x) == t_INT? expi(x): expo(x); }
1133 :
1134 : INLINE GEN
1135 599834771 : mpadd(GEN x, GEN y)
1136 : {
1137 599834771 : if (typ(x)==t_INT)
1138 16464593 : return (typ(y)==t_INT) ? addii(x,y) : addir(x,y);
1139 583370178 : return (typ(y)==t_INT) ? addir(y,x) : addrr(x,y);
1140 : }
1141 : INLINE GEN
1142 273985453 : mpsub(GEN x, GEN y)
1143 : {
1144 273985453 : if (typ(x)==t_INT)
1145 519589 : return (typ(y)==t_INT) ? subii(x,y) : subir(x,y);
1146 273465864 : return (typ(y)==t_INT) ? subri(x,y) : subrr(x,y);
1147 : }
1148 : INLINE GEN
1149 891408495 : mpmul(GEN x, GEN y)
1150 : {
1151 891408495 : if (typ(x)==t_INT)
1152 38305280 : return (typ(y)==t_INT) ? mulii(x,y) : mulir(x,y);
1153 853103215 : return (typ(y)==t_INT) ? mulir(y,x) : mulrr(x,y);
1154 : }
1155 : INLINE GEN
1156 91858363 : mpsqr(GEN x) { return (typ(x)==t_INT) ? sqri(x) : sqrr(x); }
1157 :
1158 : /* obsolete */
1159 : INLINE GEN
1160 0 : mpdiv(GEN x, GEN y)
1161 : {
1162 0 : if (typ(x) == t_REAL) return divrmp(x,y);
1163 0 : if (typ(y) == t_REAL) return divmpr(x,y);
1164 0 : pari_err_TYPE2("mpdiv",x,y); return NULL;
1165 : }
1166 : INLINE GEN
1167 1203052 : divrmp(GEN x, GEN y) { return (typ(y)==t_INT) ? divri(x,y) : divrr(x,y); }
1168 : INLINE GEN
1169 3759 : divmpr(GEN x, GEN y) { return (typ(x)==t_INT) ? divir(x,y) : divrr(x,y); }
1170 :
1171 : /*******************************************************************/
1172 : /* */
1173 : /* Z/nZ, n ULONG */
1174 : /* */
1175 : /*******************************************************************/
1176 : INLINE ulong
1177 521619649 : Fl_double(ulong a, ulong p)
1178 : {
1179 521619649 : ulong res = a << 1;
1180 521619649 : return (res >= p || res < a) ? res - p : res;
1181 : }
1182 : INLINE ulong
1183 100165434 : Fl_triple(ulong a, ulong p)
1184 : {
1185 100165434 : ulong res = a << 1;
1186 100165434 : if (res >= p || res < a) res -= p;
1187 100165434 : res += a;
1188 100165434 : return (res >= p || res < a)? res - p: res;
1189 : }
1190 : INLINE ulong
1191 20604957 : Fl_halve(ulong a, ulong p)
1192 : {
1193 : ulong ap, ap2;
1194 20604957 : if ((a&1UL)==0) return a>>1;
1195 10393376 : ap = a + p; ap2 = ap>>1;
1196 10393376 : return ap>=a ? ap2: (ap2|HIGHBIT);
1197 : }
1198 :
1199 : INLINE ulong
1200 6561519659 : Fl_add(ulong a, ulong b, ulong p)
1201 : {
1202 6561519659 : ulong res = a + b;
1203 6561519659 : return (res >= p || res < a) ? res - p : res;
1204 : }
1205 : INLINE ulong
1206 741562425 : Fl_neg(ulong x, ulong p) { return x ? p - x: 0; }
1207 :
1208 : INLINE ulong
1209 7709071619 : Fl_sub(ulong a, ulong b, ulong p)
1210 : {
1211 7709071619 : ulong res = a - b;
1212 7709071619 : return (res > a) ? res + p: res;
1213 : }
1214 :
1215 : /* centerlift(u mod p) */
1216 : INLINE long
1217 4490059 : Fl_center(ulong u, ulong p, ulong ps2) { return (long) (u > ps2)? u - p: u; }
1218 :
1219 : INLINE ulong
1220 2551531184 : Fl_mul(ulong a, ulong b, ulong p)
1221 : {
1222 : ulong x;
1223 : LOCAL_HIREMAINDER;
1224 2551531184 : x = mulll(a,b);
1225 2551531184 : if (!hiremainder) return x % p;
1226 431774045 : (void)divll(x,p); return hiremainder;
1227 : }
1228 : INLINE ulong
1229 108115541 : Fl_sqr(ulong a, ulong p)
1230 : {
1231 : ulong x;
1232 : LOCAL_HIREMAINDER;
1233 108115541 : x = mulll(a,a);
1234 108115541 : if (!hiremainder) return x % p;
1235 26959722 : (void)divll(x,p); return hiremainder;
1236 : }
1237 : /* don't assume that p is prime: can't special case a = 0 */
1238 : INLINE ulong
1239 49581353 : Fl_div(ulong a, ulong b, ulong p)
1240 49581353 : { return Fl_mul(a, Fl_inv(b, p), p); }
1241 :
1242 : /*******************************************************************/
1243 : /* */
1244 : /* DEFINED FROM EXISTING ONE EXPLOITING COMMUTATIVITY */
1245 : /* */
1246 : /*******************************************************************/
1247 : INLINE GEN
1248 1113223 : addri(GEN x, GEN y) { return addir(y,x); }
1249 : INLINE GEN
1250 195801248 : addis(GEN x, long s) { return addsi(s,x); }
1251 : INLINE GEN
1252 103816589 : addiu(GEN x, ulong s) { return addui(s,x); }
1253 : INLINE GEN
1254 12387597 : addrs(GEN x, long s) { return addsr(s,x); }
1255 :
1256 : INLINE GEN
1257 133455144 : subiu(GEN x, long y) { GEN z = subui(y, x); togglesign(z); return z; }
1258 : INLINE GEN
1259 171754 : subis(GEN x, long y) { return addsi(-y,x); }
1260 : INLINE GEN
1261 16573235 : subrs(GEN x, long y) { return addsr(-y,x); }
1262 :
1263 : INLINE GEN
1264 470856925 : mulis(GEN x, long s) { return mulsi(s,x); }
1265 : INLINE GEN
1266 378881191 : muliu(GEN x, ulong s) { return mului(s,x); }
1267 : INLINE GEN
1268 2765059 : mulru(GEN x, ulong s) { return mulur(s,x); }
1269 : INLINE GEN
1270 53835790 : mulri(GEN x, GEN s) { return mulir(s,x); }
1271 : INLINE GEN
1272 7181596 : mulrs(GEN x, long s) { return mulsr(s,x); }
1273 :
1274 : /*******************************************************************/
1275 : /* */
1276 : /* VALUATION, EXPONENT, SHIFTS */
1277 : /* */
1278 : /*******************************************************************/
1279 : INLINE long
1280 194748840 : vali(GEN x)
1281 : {
1282 : long i;
1283 : GEN xp;
1284 :
1285 194748840 : if (!signe(x)) return -1;
1286 194663992 : xp=int_LSW(x);
1287 203842709 : for (i=0; !*xp; i++) xp=int_nextW(xp);
1288 194663992 : return vals(*xp) + i * BITS_IN_LONG;
1289 : }
1290 :
1291 : /* assume x > 0 */
1292 : INLINE long
1293 955678938 : expu(ulong x) { return (BITS_IN_LONG-1) - (long)bfffo(x); }
1294 :
1295 : INLINE long
1296 2476536882 : expi(GEN x)
1297 : {
1298 2476536882 : const long lx=lgefint(x);
1299 2476536882 : return lx==2? -(long)HIGHEXPOBIT: bit_accuracy(lx)-(long)bfffo(*int_MSW(x))-1;
1300 : }
1301 :
1302 : INLINE GEN
1303 363948056 : shiftr(GEN x, long n)
1304 : {
1305 363948056 : const long e = evalexpo(expo(x)+n);
1306 363948056 : const GEN y = rcopy(x);
1307 :
1308 363948056 : if (e & ~EXPOBITS) pari_err_OVERFLOW("expo()");
1309 363948056 : y[1] = (y[1]&~EXPOBITS) | e; return y;
1310 : }
1311 : INLINE GEN
1312 153966112 : mpshift(GEN x,long s) { return (typ(x)==t_INT)?shifti(x,s):shiftr(x,s); }
1313 :
1314 : /* FIXME: adapt/use mpn_[lr]shift instead */
1315 : /* z2[imin..imax] := z1[imin..imax].f shifted left sh bits
1316 : * (feeding f from the right). Assume sh > 0 */
1317 : INLINE void
1318 7847729392 : shift_left(GEN z2, GEN z1, long imin, long imax, ulong f, ulong sh)
1319 : {
1320 7847729392 : GEN sb = z1 + imin, se = z1 + imax, te = z2 + imax;
1321 7847729392 : ulong l, m = BITS_IN_LONG - sh, k = f >> m;
1322 50280419983 : while (se > sb) {
1323 42432690591 : l = *se--;
1324 42432690591 : *te-- = (l << sh) | k;
1325 42432690591 : k = l >> m;
1326 : }
1327 7847729392 : *te = (((ulong)*se) << sh) | k;
1328 7847729392 : }
1329 : /* z2[imin..imax] := f.z1[imin..imax-1] shifted right sh bits
1330 : * (feeding f from the left). Assume sh > 0 */
1331 : INLINE void
1332 5947800597 : shift_right(GEN z2, GEN z1, long imin, long imax, ulong f, ulong sh)
1333 : {
1334 5947800597 : GEN sb = z1 + imin, se = z1 + imax, tb = z2 + imin;
1335 5947800597 : ulong k, l = *sb++, m = BITS_IN_LONG - sh;
1336 5947800597 : *tb++ = (l >> sh) | (f << m);
1337 31173357233 : while (sb < se) {
1338 25225556636 : k = l << m;
1339 25225556636 : l = *sb++;
1340 25225556636 : *tb++ = (l >> sh) | k;
1341 : }
1342 5947800597 : }
1343 :
1344 : /* Backward compatibility. Inefficient && unused */
1345 : extern ulong hiremainder;
1346 : INLINE ulong
1347 0 : shiftl(ulong x, ulong y)
1348 0 : { hiremainder = x>>(BITS_IN_LONG-y); return (x<<y); }
1349 :
1350 : INLINE ulong
1351 0 : shiftlr(ulong x, ulong y)
1352 0 : { hiremainder = x<<(BITS_IN_LONG-y); return (x>>y); }
1353 :
1354 : INLINE void
1355 472695776 : shiftr_inplace(GEN z, long d)
1356 : {
1357 472695776 : setexpo(z, expo(z)+d);
1358 472695769 : }
1359 :
1360 : /*******************************************************************/
1361 : /* */
1362 : /* ASSIGNMENT */
1363 : /* */
1364 : /*******************************************************************/
1365 : INLINE void
1366 1002981224 : affii(GEN x, GEN y)
1367 : {
1368 1002981224 : long lx = lgefint(x);
1369 1002981224 : if (lg(y)<lx) pari_err_OVERFLOW("t_INT-->t_INT assignment");
1370 38516368661 : while (--lx) y[lx] = x[lx];
1371 1002981224 : }
1372 : INLINE void
1373 6484770 : affsi(long s, GEN x)
1374 : {
1375 6484770 : if (!s) x[1] = evalsigne(0) | evallgefint(2);
1376 : else
1377 : {
1378 6208291 : if (s > 0) { x[1] = evalsigne( 1) | evallgefint(3); x[2] = s; }
1379 2082248 : else { x[1] = evalsigne(-1) | evallgefint(3); x[2] = -s; }
1380 : }
1381 6484770 : }
1382 : INLINE void
1383 47976646 : affui(ulong u, GEN x)
1384 : {
1385 47976646 : if (!u) x[1] = evalsigne(0) | evallgefint(2);
1386 47937348 : else { x[1] = evalsigne(1) | evallgefint(3); x[2] = u; }
1387 47976646 : }
1388 :
1389 : INLINE void
1390 509053161 : affsr(long x, GEN y)
1391 : {
1392 509053161 : long sh, i, ly = lg(y);
1393 :
1394 509053161 : if (!x)
1395 : {
1396 0 : y[1] = evalexpo(-bit_accuracy(ly));
1397 0 : return;
1398 : }
1399 509053161 : if (x < 0) {
1400 13805 : x = -x; sh = bfffo(x);
1401 13805 : y[1] = evalsigne(-1) | _evalexpo((BITS_IN_LONG-1)-sh);
1402 : }
1403 : else
1404 : {
1405 509039356 : sh = bfffo(x);
1406 509039356 : y[1] = evalsigne(1) | _evalexpo((BITS_IN_LONG-1)-sh);
1407 : }
1408 5220599894 : y[2] = ((ulong)x)<<sh; for (i=3; i<ly; i++) y[i]=0;
1409 : }
1410 :
1411 : INLINE void
1412 13344767 : affur(ulong x, GEN y)
1413 : {
1414 13344767 : long sh, i, ly = lg(y);
1415 :
1416 13344767 : if (!x)
1417 : {
1418 1368994 : y[1] = evalexpo(-bit_accuracy(ly));
1419 1368994 : return;
1420 : }
1421 11975773 : sh = bfffo(x);
1422 11975773 : y[1] = evalsigne(1) | _evalexpo((BITS_IN_LONG-1)-sh);
1423 43194281 : y[2] = x<<sh; for (i=3; i<ly; i++) y[i] = 0;
1424 : }
1425 :
1426 : INLINE long
1427 24078742 : hammingu(ulong x)
1428 : {
1429 : #ifdef LONG_IS_64BIT
1430 20521540 : x = (x & 0x5555555555555555UL) + ((x >> 1) & 0x5555555555555555UL);
1431 20521540 : x = (x & 0x3333333333333333UL) + ((x >> 2) & 0x3333333333333333UL);
1432 20521540 : x = (x & 0x0F0F0F0F0F0F0F0FUL) + ((x >> 4) & 0x0F0F0F0F0F0F0F0FUL);
1433 20521540 : return (x * 0x0101010101010101UL) >> 56;
1434 : #else
1435 3557202 : x = (x & 0x55555555UL) + ((x >> 1) & 0x55555555UL);
1436 3557202 : x = (x & 0x33333333UL) + ((x >> 2) & 0x33333333UL);
1437 3557202 : x = (x & 0x0F0F0F0FUL) + ((x >> 4) & 0x0F0F0F0FUL);
1438 3557202 : return (x * 0x01010101UL) >> 24;
1439 : #endif
1440 : }
1441 :
1442 : INLINE ulong
1443 3311 : thuemorseu(ulong c)
1444 : {
1445 : #ifdef LONG_IS_64BIT
1446 2838 : c ^= c >> 32;
1447 : #endif
1448 3311 : c ^= c >> 16;
1449 3311 : c ^= c >> 8;
1450 3311 : c ^= c >> 4;
1451 3311 : c ^= c >> 2;
1452 3311 : c ^= c >> 1;
1453 3311 : return c & 1;
1454 : }
|
