1 #line 2 "../src/kernel/none/level1.h"
2 /* $Id: level1.h 7892 2006-04-19 16:18:26Z kb $
3
4 Copyright (C) 2000 The PARI group.
5
6 This file is part of the PARI/GP package.
7
8 PARI/GP is free software; you can redistribute it and/or modify it under the
9 terms of the GNU General Public License as published by the Free Software
10 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
11 ANY WARRANTY WHATSOEVER.
12
13 Check the License for details. You should have received a copy of it, along
14 with the package; see the file 'COPYING'. If not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
16
17 /* This file defines "level 1" kernel functions.
18 * These functions can be inline; if not they are defined externally in
19 * level1.c, which includes this file and never needs to be changed
20 * The following lines are necessary for level0.c and level1.c */
21
22 #if !defined(INLINE)
23 GEN mkcol(GEN x);
24 GEN mkcol2(GEN x, GEN y);
25 GEN mkcolcopy(GEN x);
26 GEN mkcomplex(GEN x, GEN y);
27 GEN mkfrac(GEN x, GEN y);
28 GEN mkintmod(GEN x, GEN y);
29 GEN mkpolmod(GEN x, GEN y);
30 GEN mkintmodu(ulong x, ulong y);
31 GEN mkmat(GEN x);
32 GEN mkmat2(GEN x, GEN y);
33 GEN mkmatcopy(GEN x);
34 GEN mkrfrac(GEN x, GEN y);
35 GEN mkvec(GEN x);
36 GEN mkvec2(GEN x, GEN y);
37 GEN mkvec2s(long x, long y);
38 GEN mkvec2copy(GEN x, GEN y);
39 GEN mkvec3(GEN x, GEN y, GEN z);
40 GEN mkvec3s(long x, long y, long z);
41 GEN mkvec4(GEN x, GEN y, GEN z, GEN t);
42 GEN mkveccopy(GEN x);
43 GEN mkvecs(long x);
44 GEN mkvecsmall(long x);
45 GEN mkvecsmall2(long x,long y);
46 GEN mkvecsmall3(long x,long y, long z);
47 void affiz(GEN x, GEN y);
48 void affsz(long x, GEN y);
49 GEN addii(GEN x, GEN y);
50 GEN addir(GEN x, GEN y);
51 GEN addrr(GEN x, GEN y);
52 GEN addsi(long x, GEN y);
53 void affii(GEN x, GEN y);
54 void affsi(long s, GEN x);
55 void affsr(long s, GEN x);
56 void affui(long s, GEN x);
57 void affur(ulong s, GEN x);
58 GEN cgetc(long x);
59 GEN cgetg_copy(long lx, GEN x);
60 GEN cgetg(long x, long y);
61 GEN cgeti(long x);
62 GEN cgetr(long x);
63 int cmpir(GEN x, GEN y);
64 int cmpsr(long x, GEN y);
65 GEN constant_term(GEN x);
66 GEN ctofp(GEN x, long prec);
67 void divrrz(GEN x, GEN y, GEN z);
68 GEN divsi_rem(long x, GEN y, long *rem);
69 void divsiz(long x, GEN y, GEN z);
70 GEN divss(long x, long y);
71 GEN divss_rem(long x, long y, long *rem);
72 void divssz(long x, long y, GEN z);
73 ulong Fl_add(ulong a, ulong b, ulong p);
74 long Fl_center(ulong u, ulong p, ulong ps2);
75 ulong Fl_div(ulong a, ulong b, ulong p);
76 ulong Fl_mul(ulong a, ulong b, ulong p);
77 ulong Fl_neg(ulong x, ulong p);
78 ulong Fl_sqr(ulong a, ulong p);
79 ulong Fl_sub(ulong a, ulong b, ulong p);
80 int dvdii(GEN x, GEN y);
81 int dvdiiz(GEN x, GEN y, GEN z);
82 int dvdisz(GEN x, long y, GEN z);
83 int dvdiuz(GEN x, ulong y, GEN z);
84 void dvmdiiz(GEN x, GEN y, GEN z, GEN t);
85 GEN dvmdis(GEN x, long y, GEN *z);
86 void dvmdisz(GEN x, long y, GEN z, GEN t);
87 GEN dvmdsi(long x, GEN y, GEN *z);
88 void dvmdsiz(long x, GEN y, GEN z, GEN t);
89 GEN dvmdss(long x, long y, GEN *z);
90 void dvmdssz(long x, long y, GEN z, GEN t);
91 long evalexpo(long x);
92 long evallg(long x);
93 long evalvalp(long x);
94 long expi(GEN x);
95 GEN fractor(GEN x, long prec);
96 double gtodouble(GEN x);
97 GEN icopy_av(GEN x, GEN y);
98 GEN icopy(GEN x);
99 GEN init_gen_op(GEN x, long tx, long *lx, long *i);
100 GEN itor(GEN x, long prec);
101 long itos(GEN x);
102 long itos_or_0(GEN x);
103 ulong itou(GEN x);
104 ulong itou_or_0(GEN x);
105 GEN leading_term(GEN x);
106 long maxss(long x, long y);
107 long minss(long x, long y);
108 GEN modis(GEN x, long y);
109 GEN modsi(long x, GEN y);
110 GEN modss(long x, long y);
111 GEN mpabs(GEN x);
112 GEN mpadd(GEN x, GEN y);
113 void mpaff(GEN x, GEN y);
114 GEN mpceil(GEN x);
115 int mpcmp(GEN x, GEN y);
116 GEN mpcopy(GEN x);
117 GEN mpdiv(GEN x, GEN y);
118 GEN mpfloor(GEN x);
119 GEN mpmul(GEN x, GEN y);
120 GEN mpneg(GEN x);
121 GEN mpround(GEN x);
122 GEN mpsub(GEN x, GEN y);
123 GEN mptrunc(GEN x);
124 GEN new_chunk(size_t x);
125 long random_bits(long k);
126 GEN rdivii(GEN x, GEN y, long prec);
127 GEN rdiviiz(GEN x, GEN y, GEN z);
128 GEN rdivis(GEN x, long y, long prec);
129 GEN rdivsi(long x, GEN y, long prec);
130 GEN rdivss(long x, long y, long prec);
131 GEN real2n(long n, long prec);
132 GEN real_1(long prec);
133 GEN real_m1(long prec);
134 GEN real_0_bit(long bitprec);
135 GEN real_0(long prec);
136 void remiiz(GEN x, GEN y, GEN z);
137 GEN remis(GEN x, long y);
138 GEN remsi(long x, GEN y);
139 GEN remss(long x, long y);
140 GEN rtor(GEN x, long prec);
141 long sdivsi(long x, GEN y)
142 long sdivsi_rem(long x, GEN y, long *rem);
143 long sdivss_rem(long x, long y, long *rem);
144 void shift_left2(GEN z2, GEN z1, long min, long M, ulong f, ulong sh, ulong m);
145 void shift_right2(GEN z2, GEN z1, long min, long M, ulong f, ulong sh, ulong m);
146 ulong shiftl(ulong x, ulong y);
147 ulong shiftlr(ulong x, ulong y);
148 GEN shiftr(GEN x, long n);
149 long smodis(GEN x, long y);
150 long smodss(long x, long y);
151 void stackdummy(pari_sp av, pari_sp ltop);
152 GEN stoi(long x);
153 GEN stor(long x, long prec);
154 GEN subii(GEN x, GEN y);
155 GEN subir(GEN x, GEN y);
156 GEN subri(GEN x, GEN y);
157 GEN subrr(GEN x, GEN y);
158 GEN subsi(long x, GEN y);
159 ulong umodui(ulong x, GEN y);
160 GEN utoi(ulong x);
161 GEN utoineg(ulong x);
162 GEN utoipos(ulong x);
163 GEN utor(ulong s, long prec);
164 long vali(GEN x);
165 GEN zerocol(long n);
166 GEN zeromat(long m, long n);
167 GEN zeromatcopy(long m, long n);
168 GEN zeropadic(GEN p, long e);
169 GEN zeropol(long v);
170 GEN zeroser(long v, long e);
171 GEN zerovec(long n);
172 GEN col_ei(long n, long i);
173 GEN vec_ei(long n, long i);
174
175 #else
176
177 INLINE long
evallg(long x)178 evallg(long x)
179 {
180 if (x & ~LGBITS) pari_err(errlg);
181 return _evallg(x);
182 }
183
184 INLINE long
evalvalp(long x)185 evalvalp(long x)
186 {
187 const long v = _evalvalp(x);
188 if (v & ~VALPBITS) pari_err(errvalp);
189 return v;
190 }
191
192 INLINE long
evalexpo(long x)193 evalexpo(long x)
194 {
195 const long v = _evalexpo(x);
196 if (v & ~EXPOBITS) pari_err(errexpo);
197 return v;
198 }
199
200 INLINE GEN
constant_term(GEN x)201 constant_term(GEN x) { return signe(x)? gel(x,2): gen_0; }
202 INLINE GEN
leading_term(GEN x)203 leading_term(GEN x) { return lg(x) == 2? gen_0: gel(x,lg(x)-1); }
204
205 /* Inhibit some area gerepile-wise: declare it to be a non recursive
206 * type, of length l. Thus gerepile won't inspect the zone, just copy it.
207 * For the following situation:
208 * z = cgetg(t,a); av = avma; garbage(); ltop = avma;
209 * for (i=1; i<HUGE; i++) gel(z,i) = blah();
210 * stackdummy(av,ltop);
211 * loses (av-ltop) words but save a costly gerepile.
212 */
213 INLINE void
stackdummy(pari_sp av,pari_sp ltop)214 stackdummy(pari_sp av, pari_sp ltop) {
215 long l = ((GEN)av) - ((GEN)ltop);
216 if (l > 0) {
217 GEN z = (GEN)ltop;
218 z[0] = evaltyp(t_VECSMALL) | evallg(l);
219 #ifdef DEBUG
220 { long i; for (i = 1; i < l; i++) z[i] = 0; }
221 #endif
222 }
223 }
224
225 INLINE GEN
new_chunk(size_t x)226 new_chunk(size_t x) /* x is a number of bytes */
227 {
228 const GEN z = ((GEN) avma) - x;
229 if (x > ((avma-bot)>>TWOPOTBYTES_IN_LONG)) pari_err(errpile);
230 #if defined(_WIN32) || defined(__CYGWIN32__)
231 if (win32ctrlc) dowin32ctrlc();
232 #endif
233 avma = (pari_sp)z;
234
235 #ifdef MEMSTEP
236 if (DEBUGMEM && memused != DISABLE_MEMUSED) {
237 long d = (long)memused - (long)z;
238 if (d > 4*MEMSTEP || d < -4*MEMSTEP)
239 {
240 memused = (pari_sp)z;
241 fprintferr("...%4.0lf Mbytes used\n",(top-memused)/1048576.);
242 }
243 }
244 #endif
245 return z;
246 }
247
248 /* cgetg(lg(x), typ(x)), assuming lx = lg(x). Implicit unsetisclone() */
249 INLINE GEN
cgetg_copy(long lx,GEN x)250 cgetg_copy(long lx, GEN x) {
251 GEN y = new_chunk((size_t)lx);
252 y[0] = x[0] & (TYPBITS|LGBITS); return y;
253 }
254 INLINE GEN
init_gen_op(GEN x,long tx,long * lx,long * i)255 init_gen_op(GEN x, long tx, long *lx, long *i) {
256 GEN y;
257 *lx = lg(x); y = cgetg_copy(*lx, x);
258 if (lontyp[tx] == 1) *i = 1; else { y[1] = x[1]; *i = 2; }
259 return y;
260 }
261
262 INLINE GEN
cgetg(long x,long y)263 cgetg(long x, long y)
264 {
265 const GEN z = new_chunk((size_t)x);
266 z[0] = evaltyp(y) | evallg(x);
267 return z;
268 }
269
270 INLINE GEN
cgeti(long x)271 cgeti(long x)
272 {
273 const GEN z = new_chunk((size_t)x);
274 z[0] = evaltyp(t_INT) | evallg(x);
275 return z;
276 }
277
278 INLINE GEN
cgetr(long x)279 cgetr(long x)
280 {
281 const GEN z = new_chunk((size_t)x);
282 z[0] = evaltyp(t_REAL) | evallg(x);
283 return z;
284 }
285 INLINE GEN
cgetc(long l)286 cgetc(long l)
287 {
288 GEN u = cgetg(3,t_COMPLEX); gel(u,1) = cgetr(l); gel(u,2) = cgetr(l);
289 return u;
290 }
291 INLINE GEN
mkintmod(GEN x,GEN y)292 mkintmod(GEN x, GEN y) { GEN v = cgetg(3, t_INTMOD);
293 gel(v,1) = y; gel(v,2) = x; return v; }
294 INLINE GEN
mkpolmod(GEN x,GEN y)295 mkpolmod(GEN x, GEN y) { GEN v = cgetg(3, t_POLMOD);
296 gel(v,1) = y; gel(v,2) = x; return v; }
297 INLINE GEN
mkfrac(GEN x,GEN y)298 mkfrac(GEN x, GEN y) { GEN v = cgetg(3, t_FRAC);
299 gel(v,1) = x; gel(v,2) = y; return v; }
300 INLINE GEN
mkrfrac(GEN x,GEN y)301 mkrfrac(GEN x, GEN y) { GEN v = cgetg(3, t_RFRAC);
302 gel(v,1) = x; gel(v,2) = y; return v; }
303 INLINE GEN
mkcomplex(GEN x,GEN y)304 mkcomplex(GEN x, GEN y) { GEN v = cgetg(3, t_COMPLEX);
305 gel(v,1) = x; gel(v,2) = y; return v; }
306 INLINE GEN
mkvec(GEN x)307 mkvec(GEN x) { GEN v = cgetg(2, t_VEC); gel(v,1) = x; return v; }
308 INLINE GEN
mkvecsmall(long x)309 mkvecsmall(long x) { GEN v = cgetg(2, t_VECSMALL); v[1] = x; return v; }
310 INLINE GEN
mkvecsmall2(long x,long y)311 mkvecsmall2(long x,long y) { GEN v = cgetg(3, t_VECSMALL);
312 v[1]=x; v[2]=y; return v; }
313 INLINE GEN
mkvecsmall3(long x,long y,long z)314 mkvecsmall3(long x,long y, long z) { GEN v = cgetg(4, t_VECSMALL);
315 v[1]=x; v[2]=y; v[3]=z; return v; }
316 INLINE GEN
mkveccopy(GEN x)317 mkveccopy(GEN x) { GEN v = cgetg(2, t_VEC); gel(v,1) = gcopy(x); return v; }
318 INLINE GEN
mkvec2(GEN x,GEN y)319 mkvec2(GEN x, GEN y) {
320 GEN v = cgetg(3,t_VEC); gel(v,1) = x; gel(v,2) = y; return v; }
321 INLINE GEN
mkvec3(GEN x,GEN y,GEN z)322 mkvec3(GEN x, GEN y, GEN z) {
323 GEN v=cgetg(4,t_VEC); gel(v,1) = x; gel(v,2) = y; gel(v,3) = z; return v; }
324 INLINE GEN
mkvec4(GEN x,GEN y,GEN z,GEN t)325 mkvec4(GEN x, GEN y, GEN z, GEN t) {
326 GEN v=cgetg(5,t_VEC); gel(v,1) = x; gel(v,2) = y; gel(v,3) = z; gel(v,4) = t;
327 return v; }
328 INLINE GEN
mkvec2copy(GEN x,GEN y)329 mkvec2copy(GEN x, GEN y) {
330 GEN v = cgetg(3,t_VEC); gel(v,1) = gcopy(x); gel(v,2) = gcopy(y); return v; }
331 INLINE GEN
mkcol(GEN x)332 mkcol(GEN x) { GEN v = cgetg(2, t_COL); gel(v,1) = x; return v; }
333 INLINE GEN
mkcol2(GEN x,GEN y)334 mkcol2(GEN x, GEN y) {
335 GEN v = cgetg(3,t_COL); gel(v,1) = x; gel(v,2) = y; return v; }
336 INLINE GEN
mkcolcopy(GEN x)337 mkcolcopy(GEN x) { GEN v = cgetg(2, t_COL); gel(v,1) = gcopy(x); return v; }
338 INLINE GEN
mkmat(GEN x)339 mkmat(GEN x) { GEN v = cgetg(2, t_MAT); gel(v,1) = x; return v; }
340 INLINE GEN
mkmat2(GEN x,GEN y)341 mkmat2(GEN x, GEN y) { GEN v=cgetg(3,t_MAT); gel(v,1)=x; gel(v,2)=y; return v; }
342 INLINE GEN
mkmatcopy(GEN x)343 mkmatcopy(GEN x) { GEN v = cgetg(2, t_MAT); gel(v,1) = gcopy(x); return v; }
344
345 /*** ZERO ***/
346 /* O(p^e) */
347 INLINE GEN
zeropadic(GEN p,long e)348 zeropadic(GEN p, long e)
349 {
350 GEN y = cgetg(5,t_PADIC);
351 gel(y,4) = gen_0;
352 gel(y,3) = gen_1;
353 copyifstack(p,y[2]);
354 y[1] = evalvalp(e) | evalprecp(0);
355 return y;
356 }
357 /* O(pol_x[v]^e) */
358 INLINE GEN
zeroser(long v,long e)359 zeroser(long v, long e)
360 {
361 GEN x = cgetg(2, t_SER);
362 x[1] = evalvalp(e) | evalvarn(v); return x;
363 }
364 /* 0 * pol_x[v] */
365 INLINE GEN
zeropol(long v)366 zeropol(long v)
367 {
368 GEN x = cgetg(2,t_POL);
369 x[1] = evalvarn(v); return x;
370 }
371 /* vector(n) */
372 INLINE GEN
zerocol(long n)373 zerocol(long n)
374 {
375 GEN y = cgetg(n+1,t_COL);
376 long i; for (i=1; i<=n; i++) gel(y,i) = gen_0;
377 return y;
378 }
379 /* vectorv(n) */
380 INLINE GEN
zerovec(long n)381 zerovec(long n)
382 {
383 GEN y = cgetg(n+1,t_VEC);
384 long i; for (i=1; i<=n; i++) gel(y,i) = gen_0;
385 return y;
386 }
387 /* matrix(m, n) */
388 INLINE GEN
zeromat(long m,long n)389 zeromat(long m, long n)
390 {
391 GEN y = cgetg(n+1,t_MAT);
392 GEN v = zerocol(m);
393 long i; for (i=1; i<=n; i++) gel(y,i) = v;
394 return y;
395 }
396 /* matrix(m, n) */
397 INLINE GEN
zeromatcopy(long m,long n)398 zeromatcopy(long m, long n)
399 {
400 GEN y = cgetg(n+1,t_MAT);
401 long i; for (i=1; i<=n; i++) gel(y,i) = zerocol(m);
402 return y;
403 }
404 /* i-th vector in the standard basis */
405 INLINE GEN
col_ei(long n,long i)406 col_ei(long n, long i) { GEN e = zerocol(n); gel(e,i) = gen_1; return e; }
407 INLINE GEN
vec_ei(long n,long i)408 vec_ei(long n, long i) { GEN e = zerovec(n); gel(e,i) = gen_1; return e; }
409
410 /* cannot do memcpy because sometimes x and y overlap */
411 INLINE GEN
mpcopy(GEN x)412 mpcopy(GEN x)
413 {
414 register long lx = lg(x);
415 const GEN y = cgetg_copy(lx, x);
416 while (--lx > 0) y[lx] = x[lx];
417 return y;
418 }
419
420 INLINE GEN
icopy(GEN x)421 icopy(GEN x)
422 {
423 register long lx = lgefint(x);
424 const GEN y = cgeti(lx);
425 while (--lx > 0) y[lx] = x[lx];
426 return y;
427 }
428
429 /* copy integer x as if we had avma = av */
430 INLINE GEN
icopy_av(GEN x,GEN y)431 icopy_av(GEN x, GEN y)
432 {
433 register long lx = lgefint(x);
434 register long ly = lx;
435 y -= lx;
436 while (--lx > 0) y[lx]=x[lx];
437 y[0] = evaltyp(t_INT)|evallg(ly);
438 return y;
439 }
440
441 INLINE GEN
mpneg(GEN x)442 mpneg(GEN x)
443 {
444 const GEN y=mpcopy(x);
445 setsigne(y,-signe(x)); return y;
446 }
447
448 INLINE GEN
mpabs(GEN x)449 mpabs(GEN x)
450 {
451 const GEN y=mpcopy(x);
452 if (signe(x)<0) setsigne(y,1);
453 return y;
454 }
455
456 INLINE long
smodis(GEN x,long y)457 smodis(GEN x, long y)
458 {
459 long rem;
460 const pari_sp av=avma; (void)divis_rem(x,y, &rem); avma=av;
461 return (rem >= 0) ? rem: labs(y) + rem;
462 }
463
464 INLINE long
smodss(long x,long y)465 smodss(long x, long y)
466 {
467 long rem = x%y;
468 return (rem >= 0)? rem: labs(y) + rem;
469 }
470
471 /* assume x != 0, return -x as a t_INT */
472 INLINE GEN
utoineg(ulong x)473 utoineg(ulong x)
474 {
475 GEN y = cgeti(3);
476 y[1] = evalsigne(-1)| evallgefint(3); y[2] = x; return y;
477 }
478 /* assume x != 0, return utoi(x) */
479 INLINE GEN
utoipos(ulong x)480 utoipos(ulong x)
481 {
482 GEN y = cgeti(3);
483 y[1] = evalsigne(1)| evallgefint(3); y[2] = x; return y;
484 }
485 INLINE GEN
utoi(ulong x)486 utoi(ulong x) { return x? utoipos(x): gen_0; }
487 INLINE GEN
stoi(long x)488 stoi(long x)
489 {
490 if (!x) return gen_0;
491 return x > 0? utoipos((ulong)x): utoineg((ulong)-x);
492 }
493 INLINE GEN
mkvecs(long x)494 mkvecs(long x) { GEN v = cgetg(2, t_VEC); gel(v,1) = stoi(x); return v; }
495 INLINE GEN
mkvec2s(long x,long y)496 mkvec2s(long x, long y) {
497 GEN v = cgetg(3,t_VEC); gel(v,1) = stoi(x); gel(v,2) = stoi(y); return v; }
498 INLINE GEN
mkvec3s(long x,long y,long z)499 mkvec3s(long x, long y, long z) {
500 GEN v=cgetg(4,t_VEC);
501 gel(v,1)=stoi(x); gel(v,2)=stoi(y); gel(v,3)=stoi(z); return v;
502 }
503 INLINE GEN
mkintmodu(ulong x,ulong y)504 mkintmodu(ulong x, ulong y) {
505 GEN v = cgetg(3,t_INTMOD);
506 gel(v,1) = utoipos(y);
507 gel(v,2) = utoi(x); return v;
508 }
509
510 INLINE GEN
stosmall(long x)511 stosmall(long x)
512 {
513 if (labs(x) & SMALL_MASK) return stoi(x);
514 return (GEN) (1 | (x<<1));
515 }
516
517 INLINE long
itos(GEN x)518 itos(GEN x)
519 {
520 const long s = signe(x);
521 long u;
522
523 if (!s) return 0;
524 u = x[2]; if (lgefint(x) > 3 || u < 0) pari_err(affer2);
525 return (s>0) ? u : -u;
526 }
527
528 /* as itos, but return 0 if too large. Cf is_bigint */
529 INLINE long
itos_or_0(GEN x)530 itos_or_0(GEN x) {
531 long n;
532 if (lgefint(x) != 3 || (n = x[2]) & HIGHBIT) return 0;
533 return signe(x) > 0? n: -n;
534 }
535 /* as itou, but return 0 if too large. Cf is_bigint */
536 INLINE ulong
itou_or_0(GEN x)537 itou_or_0(GEN x) {
538 if (lgefint(x) != 3) return 0;
539 return (ulong)x[2];
540 }
541
542 INLINE GEN
modss(long x,long y)543 modss(long x, long y) { return stoi(smodss(x, y)); }
544
545 INLINE GEN
remss(long x,long y)546 remss(long x, long y) { return stoi(x % y); }
547
548 INLINE void
affii(GEN x,GEN y)549 affii(GEN x, GEN y)
550 {
551 long lx;
552
553 if (x==y) return;
554 lx=lgefint(x); if (lg(y)<lx) pari_err(affer3);
555 while (--lx) y[lx]=x[lx];
556 }
557
558 INLINE void
affsi(long s,GEN x)559 affsi(long s, GEN x)
560 {
561 if (!s) x[1] = evalsigne(0) | evallgefint(2);
562 else
563 {
564 if (s > 0) { x[1] = evalsigne( 1) | evallgefint(3); x[2] = s; }
565 else { x[1] = evalsigne(-1) | evallgefint(3); x[2] = -s; }
566 }
567 }
568
569 INLINE void
affsr(long x,GEN y)570 affsr(long x, GEN y)
571 {
572 long sh, i, ly = lg(y);
573
574 if (!x)
575 {
576 y[1] = evalexpo(-bit_accuracy(ly));
577 return;
578 }
579 if (x < 0) {
580 x = -x; sh = bfffo(x);
581 y[1] = evalsigne(-1) | _evalexpo((BITS_IN_LONG-1)-sh);
582 }
583 else
584 {
585 sh = bfffo(x);
586 y[1] = evalsigne(1) | _evalexpo((BITS_IN_LONG-1)-sh);
587 }
588 y[2] = x<<sh; for (i=3; i<ly; i++) y[i]=0;
589 }
590
591 INLINE void
affur(ulong x,GEN y)592 affur(ulong x, GEN y)
593 {
594 long sh, i, ly = lg(y);
595
596 if (!x)
597 {
598 y[1] = evalexpo(-bit_accuracy(ly));
599 return;
600 }
601 sh = bfffo(x);
602 y[1] = evalsigne(1) | _evalexpo((BITS_IN_LONG-1)-sh);
603 y[2] = x<<sh; for (i=3; i<ly; i++) y[i] = 0;
604 }
605
606 INLINE void
affiz(GEN x,GEN y)607 affiz(GEN x, GEN y) { if (typ(y)==t_INT) affii(x,y); else affir(x,y); }
608 INLINE void
affsz(long x,GEN y)609 affsz(long x, GEN y) { if (typ(y)==t_INT) affsi(x,y); else affsr(x,y); }
610 INLINE void
mpaff(GEN x,GEN y)611 mpaff(GEN x, GEN y) { if (typ(x)==t_INT) affiz(x, y); else affrr(x,y); }
612
613 INLINE GEN
real_0_bit(long bitprec)614 real_0_bit(long bitprec) { GEN x=cgetr(2); x[1]=evalexpo(bitprec); return x; }
615 INLINE GEN
real_0(long prec)616 real_0(long prec) { return real_0_bit(-bit_accuracy(prec)); }
617 INLINE GEN
real_1(long prec)618 real_1(long prec) {
619 GEN x = cgetr(prec);
620 long i;
621 x[1] = evalsigne(1) | _evalexpo(0);
622 x[2] = (long)HIGHBIT; for (i=3; i<prec; i++) x[i] = 0;
623 return x;
624 }
625 INLINE GEN
real_m1(long prec)626 real_m1(long prec) {
627 GEN x = cgetr(prec);
628 long i;
629 x[1] = evalsigne(-1) | _evalexpo(0);
630 x[2] = (long)HIGHBIT; for (i=3; i<prec; i++) x[i] = 0;
631 return x;
632 }
633 /* 2.^n */
634 INLINE GEN
real2n(long n,long prec)635 real2n(long n, long prec) { GEN z = real_1(prec); setexpo(z, n); return z; }
636 INLINE GEN
stor(long s,long prec)637 stor(long s, long prec) { GEN z = cgetr(prec); affsr(s,z); return z; }
638 INLINE GEN
utor(ulong s,long prec)639 utor(ulong s, long prec){ GEN z = cgetr(prec); affur(s,z); return z; }
640 INLINE GEN
itor(GEN x,long prec)641 itor(GEN x, long prec) { GEN z = cgetr(prec); affir(x,z); return z; }
642 INLINE GEN
rtor(GEN x,long prec)643 rtor(GEN x, long prec) { GEN z = cgetr(prec); affrr(x,z); return z; }
644 INLINE GEN
ctofp(GEN x,long prec)645 ctofp(GEN x, long prec) { GEN z = cgetg(3,t_COMPLEX);
646 gel(z,1) = gtofp(gel(x,1),prec);
647 gel(z,2) = gtofp(gel(x,2),prec); return z;
648 }
649
650 INLINE GEN
shiftr(GEN x,long n)651 shiftr(GEN x, long n)
652 {
653 const long e = evalexpo(expo(x)+n);
654 const GEN y = rcopy(x);
655
656 if (e & ~EXPOBITS) pari_err(talker,"overflow in real shift");
657 y[1] = (y[1]&~EXPOBITS) | e; return y;
658 }
659
660 INLINE int
cmpir(GEN x,GEN y)661 cmpir(GEN x, GEN y)
662 {
663 pari_sp av;
664 GEN z;
665
666 if (!signe(x)) return -signe(y);
667 if (!signe(y)) return signe(x);
668 av=avma; z = itor(x, lg(y)); avma=av;
669 return cmprr(z,y); /* cmprr does no memory adjustment */
670 }
671
672 INLINE int
cmpsr(long x,GEN y)673 cmpsr(long x, GEN y)
674 {
675 pari_sp av;
676 GEN z;
677
678 if (!x) return -signe(y);
679 av=avma; z = stor(x, 3); avma=av;
680 return cmprr(z,y);
681 }
682
683 INLINE long
maxss(long x,long y)684 maxss(long x, long y) { return x>y?x:y; }
685 INLINE long
minss(long x,long y)686 minss(long x, long y) { return x<y?x:y; }
687
688 INLINE GEN
subii(GEN x,GEN y)689 subii(GEN x, GEN y)
690 {
691 if (x==y) return gen_0; /* frequent with x = y = gen_0 */
692 return addii_sign(x, signe(x), y, -signe(y));
693 }
694 INLINE GEN
addii(GEN x,GEN y)695 addii(GEN x, GEN y) { return addii_sign(x, signe(x), y, signe(y)); }
696 INLINE GEN
addrr(GEN x,GEN y)697 addrr(GEN x, GEN y) { return addrr_sign(x, signe(x), y, signe(y)); }
698 INLINE GEN
subrr(GEN x,GEN y)699 subrr(GEN x, GEN y) { return addrr_sign(x, signe(x), y, -signe(y)); }
700 INLINE GEN
addir(GEN x,GEN y)701 addir(GEN x, GEN y) { return addir_sign(x, signe(x), y, signe(y)); }
702 INLINE GEN
subir(GEN x,GEN y)703 subir(GEN x, GEN y) { return addir_sign(x, signe(x), y, -signe(y)); }
704 INLINE GEN
subri(GEN x,GEN y)705 subri(GEN x, GEN y) { return addir_sign(y, -signe(y), x, signe(x)); }
706 INLINE GEN
addsi(long x,GEN y)707 addsi(long x, GEN y) { return addsi_sign(x, y, signe(y)); }
708 INLINE GEN
subsi(long x,GEN y)709 subsi(long x, GEN y) { return addsi_sign(x, y, -signe(y)); }
710
711 INLINE long
vali(GEN x)712 vali(GEN x)
713 {
714 long i;
715 GEN xp;
716
717 if (!signe(x)) return -1;
718 xp=int_LSW(x);
719 for (i=0; !*xp; i++) xp=int_nextW(xp);
720 return (i<<TWOPOTBITS_IN_LONG) + vals(*xp);
721 }
722
723 INLINE GEN
divss(long x,long y)724 divss(long x, long y) { return stoi(x / y); }
725
726 INLINE long
sdivss_rem(long x,long y,long * rem)727 sdivss_rem(long x, long y, long *rem)
728 {
729 long q;
730 LOCAL_HIREMAINDER;
731 if (!y) pari_err(gdiver);
732 hiremainder = 0; q = divll((ulong)labs(x),(ulong)labs(y));
733 if (x < 0) { hiremainder = -((long)hiremainder); q = -q; }
734 if (y < 0) q = -q;
735 *rem = hiremainder; return q;
736 }
737
738 INLINE GEN
divss_rem(long x,long y,long * rem)739 divss_rem(long x, long y, long *rem) { return stoi(sdivss_rem(x,y,rem)); }
740
741 INLINE GEN
dvmdss(long x,long y,GEN * z)742 dvmdss(long x, long y, GEN *z)
743 {
744 long rem;
745 const GEN q = divss_rem(x,y, &rem);
746 *z = stoi(rem); return q;
747 }
748
749 INLINE long
sdivsi_rem(long x,GEN y,long * rem)750 sdivsi_rem(long x, GEN y, long *rem)
751 {
752 long q, s = signe(y);
753 LOCAL_HIREMAINDER;
754
755 if (!s) pari_err(gdiver);
756 if (!x || lgefint(y)>3 || ((long)y[2]) < 0) { *rem = x; return 0; }
757 hiremainder=0; q = (long)divll(labs(x), (ulong)y[2]);
758 if (x < 0) { hiremainder = -((long)hiremainder); q = -q; }
759 if (s < 0) q = -q;
760 *rem = hiremainder; return q;
761 }
762
763 INLINE long
sdivsi(long x,GEN y)764 sdivsi(long x, GEN y)
765 {
766 long q, s = signe(y);
767
768 if (!s) pari_err(gdiver);
769 if (!x || lgefint(y)>3 || ((long)y[2]) < 0) return 0;
770 q = labs(x) / y[2];
771 if (x < 0) q = -q;
772 if (s < 0) q = -q;
773 return q;
774 }
775
776 INLINE GEN
modsi(long x,GEN y)777 modsi(long x, GEN y) {
778 long r;
779 (void)sdivsi_rem(x, y, &r);
780 return (r >= 0)? stoi(r): addsi_sign(r, y, 1);
781 }
782
783 INLINE GEN
divsi_rem(long s,GEN y,long * rem)784 divsi_rem(long s, GEN y, long *rem) { return stoi(sdivsi_rem(s,y,rem)); }
785
786 INLINE GEN
dvmdsi(long x,GEN y,GEN * z)787 dvmdsi(long x, GEN y, GEN *z)
788 {
789 long rem;
790 const GEN q = divsi_rem(x,y, &rem);
791 *z = stoi(rem); return q;
792 }
793
794 INLINE GEN
dvmdis(GEN x,long y,GEN * z)795 dvmdis(GEN x, long y, GEN *z)
796 {
797 long rem;
798 const GEN q = divis_rem(x,y, &rem);
799 *z = stoi(rem); return q;
800 }
801
802 INLINE void
dvmdssz(long x,long y,GEN z,GEN t)803 dvmdssz(long x, long y, GEN z, GEN t)
804 {
805 long rem;
806 const pari_sp av=avma; affiz(divss_rem(x,y, &rem), z); avma=av;
807 affsi(rem,t);
808 }
809
810 INLINE void
dvmdsiz(long x,GEN y,GEN z,GEN t)811 dvmdsiz(long x, GEN y, GEN z, GEN t)
812 {
813 long rem;
814 const pari_sp av = avma; affiz(divsi_rem(x,y, &rem), z); avma = av;
815 affsi(rem,t);
816 }
817
818 INLINE void
dvmdisz(GEN x,long y,GEN z,GEN t)819 dvmdisz(GEN x, long y, GEN z, GEN t)
820 {
821 long rem;
822 const pari_sp av=avma; affiz(divis_rem(x,y, &rem),z); avma=av;
823 affsz(rem,t);
824 }
825
826 INLINE void
dvmdiiz(GEN x,GEN y,GEN z,GEN t)827 dvmdiiz(GEN x, GEN y, GEN z, GEN t)
828 {
829 const pari_sp av=avma;
830 GEN p;
831
832 affiz(dvmdii(x,y,&p),z); affiz(p,t); avma=av;
833 }
834
835 INLINE GEN
modis(GEN x,long y)836 modis(GEN x, long y)
837 {
838 return stoi(smodis(x,y));
839 }
840
841 INLINE ulong
umodui(ulong x,GEN y)842 umodui(ulong x, GEN y)
843 {
844 LOCAL_HIREMAINDER;
845 if (!signe(y)) pari_err(gdiver);
846 if (!x || lgefint(y) > 3) return x;
847 hiremainder = 0; (void)divll(x, y[2]); return hiremainder;
848 }
849
850 INLINE GEN
remsi(long x,GEN y)851 remsi(long x, GEN y)
852 {
853 long rem;
854 const pari_sp av=avma; (void)divsi_rem(x,y, &rem); avma=av;
855 return stoi(rem);
856 }
857
858 INLINE GEN
remis(GEN x,long y)859 remis(GEN x, long y)
860 {
861 long rem;
862 const pari_sp av=avma; (void)divis_rem(x,y, &rem); avma=av;
863 return stoi(rem);
864 }
865
866 INLINE GEN
rdivis(GEN x,long y,long prec)867 rdivis(GEN x, long y, long prec)
868 {
869 GEN z = cgetr(prec);
870 const pari_sp av = avma;
871 affrr(divrs(itor(x,prec), y),z);
872 avma = av; return z;
873 }
874
875 INLINE void
divsiz(long x,GEN y,GEN z)876 divsiz(long x, GEN y, GEN z)
877 {
878 long junk;
879 affsi(sdivsi_rem(x,y,&junk), z);
880 }
881
882 INLINE void
divssz(long x,long y,GEN z)883 divssz(long x, long y, GEN z) { affsi(x/y, z); }
884
885 INLINE GEN
rdivsi(long x,GEN y,long prec)886 rdivsi(long x, GEN y, long prec)
887 {
888 GEN z = cgetr(prec);
889 const pari_sp av = avma;
890 affrr(divsr(x, itor(y,prec)), z);
891 avma = av; return z;
892 }
893
894 INLINE GEN
rdivss(long x,long y,long prec)895 rdivss(long x, long y, long prec)
896 {
897 GEN z = cgetr(prec);
898 const pari_sp av = avma;
899 affrr(divrs(stor(x, prec), y), z);
900 avma = av; return z;
901 }
902
903 INLINE GEN
rdiviiz(GEN x,GEN y,GEN z)904 rdiviiz(GEN x, GEN y, GEN z)
905 {
906 const pari_sp av = avma;
907 long prec = lg(z);
908 affir(x, z);
909 if (!is_bigint(y)) {
910 affrr(divrs(z, y[2]), z);
911 if (signe(y) < 0) setsigne(z, -signe(z));
912 }
913 else
914 affrr(divrr(z, itor(y,prec)), z);
915 avma = av; return z;
916 }
917
918 INLINE GEN
rdivii(GEN x,GEN y,long prec)919 rdivii(GEN x, GEN y, long prec) { return rdiviiz(x, y, cgetr(prec)); }
920 INLINE GEN
fractor(GEN x,long prec)921 fractor(GEN x, long prec) { return rdivii(gel(x,1), gel(x,2), prec); }
922
923 INLINE void
divrrz(GEN x,GEN y,GEN z)924 divrrz(GEN x, GEN y, GEN z)
925 {
926 const pari_sp av=avma;
927 affrr(divrr(x,y),z); avma=av;
928 }
929
930 INLINE void
remiiz(GEN x,GEN y,GEN z)931 remiiz(GEN x, GEN y, GEN z)
932 {
933 const pari_sp av=avma;
934 affii(remii(x,y),z); avma=av;
935 }
936
937 INLINE int
dvdii(GEN x,GEN y)938 dvdii(GEN x, GEN y)
939 {
940 const pari_sp av=avma;
941 const GEN p1=remii(x,y);
942 avma=av; return p1 == gen_0;
943 }
944
945 INLINE int
mpcmp(GEN x,GEN y)946 mpcmp(GEN x, GEN y)
947 {
948 if (typ(x)==t_INT)
949 return (typ(y)==t_INT) ? cmpii(x,y) : cmpir(x,y);
950 return (typ(y)==t_INT) ? -cmpir(y,x) : cmprr(x,y);
951 }
952
953 INLINE GEN
mptrunc(GEN x)954 mptrunc(GEN x) { return typ(x)==t_INT? icopy(x): truncr(x); }
955 INLINE GEN
mpfloor(GEN x)956 mpfloor(GEN x) { return typ(x)==t_INT? icopy(x): floorr(x); }
957 INLINE GEN
mpceil(GEN x)958 mpceil(GEN x) { return typ(x)==t_INT? icopy(x): ceilr(x); }
959 INLINE GEN
mpround(GEN x)960 mpround(GEN x) { return typ(x) == t_INT? icopy(x): roundr(x); }
961
962 INLINE GEN
mpadd(GEN x,GEN y)963 mpadd(GEN x, GEN y)
964 {
965 if (typ(x)==t_INT)
966 return (typ(y)==t_INT) ? addii(x,y) : addir(x,y);
967 return (typ(y)==t_INT) ? addir(y,x) : addrr(x,y);
968 }
969
970 INLINE GEN
mpsub(GEN x,GEN y)971 mpsub(GEN x, GEN y)
972 {
973 if (typ(x)==t_INT)
974 return (typ(y)==t_INT) ? subii(x,y) : subir(x,y);
975 return (typ(y)==t_INT) ? subri(x,y) : subrr(x,y);
976 }
977
978 INLINE GEN
mpmul(GEN x,GEN y)979 mpmul(GEN x, GEN y)
980 {
981 if (typ(x)==t_INT)
982 return (typ(y)==t_INT) ? mulii(x,y) : mulir(x,y);
983 return (typ(y)==t_INT) ? mulir(y,x) : mulrr(x,y);
984 }
985
986 INLINE GEN
mpdiv(GEN x,GEN y)987 mpdiv(GEN x, GEN y)
988 {
989 if (typ(x)==t_INT)
990 return (typ(y)==t_INT) ? divii(x,y) : divir(x,y);
991 return (typ(y)==t_INT) ? divri(x,y) : divrr(x,y);
992 }
993
994 INLINE int
dvdiiz(GEN x,GEN y,GEN z)995 dvdiiz(GEN x, GEN y, GEN z)
996 {
997 const pari_sp av=avma;
998 GEN p2;
999 const GEN p1=dvmdii(x,y,&p2);
1000
1001 if (signe(p2)) { avma=av; return 0; }
1002 affii(p1,z); avma=av; return 1;
1003 }
1004
1005 /* assume 0 <= k < 32. Return random 0 <= x < (1<<k) */
1006 INLINE long
random_bits(long k)1007 random_bits(long k) { return pari_rand31() >> (31 - k); }
1008
1009 INLINE ulong
itou(GEN x)1010 itou(GEN x)
1011 {
1012 const long s = signe(x);
1013
1014 if (!s) return 0;
1015 if (lgefint(x) > 3) pari_err(affer2);
1016 return x[2];
1017 }
1018
1019 INLINE void
affui(ulong u,GEN x)1020 affui(ulong u, GEN x)
1021 {
1022 if (!u) x[1] = evalsigne(0) | evallgefint(2);
1023 else { x[1] = evalsigne(1) | evallgefint(3); x[2] = u; }
1024 }
1025
1026 INLINE int
dvdisz(GEN x,long y,GEN z)1027 dvdisz(GEN x, long y, GEN z)
1028 {
1029 const pari_sp av = avma;
1030 long rem;
1031 GEN p1 = divis_rem(x,y, &rem);
1032 avma = av; if (rem) return 0;
1033 affii(p1,z); return 1;
1034 }
1035
1036 INLINE int
dvdiuz(GEN x,ulong y,GEN z)1037 dvdiuz(GEN x, ulong y, GEN z)
1038 {
1039 const pari_sp av = avma;
1040 ulong rem;
1041 GEN p1 = diviu_rem(x,y, &rem);
1042 avma = av; if (rem) return 0;
1043 affii(p1,z); return 1;
1044 }
1045
1046 INLINE double
gtodouble(GEN x)1047 gtodouble(GEN x)
1048 {
1049 static long reel4[4]={ evaltyp(t_REAL) | _evallg(4),0,0,0 };
1050
1051 if (typ(x)==t_REAL) return rtodbl(x);
1052 gaffect(x,(GEN)reel4); return rtodbl((GEN)reel4);
1053 }
1054
1055 /* same as Fl_add, assume p <= 2^(BIL - 1), so that overflow can't occur */
1056 INLINE ulong
Fl_add_noofl(ulong a,ulong b,ulong p)1057 Fl_add_noofl(ulong a, ulong b, ulong p)
1058 {
1059 ulong res = a + b;
1060 return (res >= p) ? res - p : res;
1061 }
1062 INLINE ulong
Fl_add(ulong a,ulong b,ulong p)1063 Fl_add(ulong a, ulong b, ulong p)
1064 {
1065 ulong res = a + b;
1066 return (res >= p || res < a) ? res - p : res;
1067 }
1068 INLINE ulong
Fl_neg(ulong x,ulong p)1069 Fl_neg(ulong x, ulong p) { return x ? p - x: 0; }
1070
1071 INLINE ulong
Fl_sub(ulong a,ulong b,ulong p)1072 Fl_sub(ulong a, ulong b, ulong p)
1073 {
1074 ulong res = a - b;
1075 return (res > a) ? res + p: res;
1076 }
1077
1078 /* centerlift(u mod p) */
1079 INLINE long
Fl_center(ulong u,ulong p,ulong ps2)1080 Fl_center(ulong u, ulong p, ulong ps2) { return (long) (u > ps2)? u - p: u; }
1081
1082 INLINE ulong
Fl_mul(ulong a,ulong b,ulong p)1083 Fl_mul(ulong a, ulong b, ulong p)
1084 {
1085 LOCAL_HIREMAINDER;
1086 {
1087 register ulong x = mulll(a,b);
1088 (void)divll(x,p);
1089 }
1090 return hiremainder;
1091 }
1092 INLINE ulong
Fl_sqr(ulong a,ulong p)1093 Fl_sqr(ulong a, ulong p)
1094 {
1095 LOCAL_HIREMAINDER;
1096 {
1097 register ulong x = mulll(a,a);
1098 (void)divll(x,p);
1099 }
1100 return hiremainder;
1101 }
1102 INLINE ulong
Fl_div(ulong a,ulong b,ulong p)1103 Fl_div(ulong a, ulong b, ulong p)
1104 {
1105 return Fl_mul(a, Fl_inv(b, p), p);
1106 }
1107
1108 INLINE long
expi(GEN x)1109 expi(GEN x)
1110 {
1111 const long lx=lgefint(x);
1112 return lx==2? -(long)HIGHEXPOBIT: bit_accuracy(lx)-(long)bfffo(*int_MSW(x))-1;
1113 }
1114
1115 /* z2[imin..imax] := z1[imin..imax].f shifted left sh bits
1116 * (feeding f from the right). Assume sh > 0 */
1117 INLINE void
shift_left2(GEN z2,GEN z1,long imin,long imax,ulong f,ulong sh,ulong m)1118 shift_left2(GEN z2, GEN z1, long imin, long imax, ulong f, ulong sh, ulong m)
1119 {
1120 GEN sb = z1 + imin, se = z1 + imax, te = z2 + imax;
1121 ulong l, k = f >> m;
1122 while (se > sb) {
1123 l = *se--;
1124 *te-- = (l << sh) | k;
1125 k = l >> m;
1126 }
1127 *te = (*se << sh) | k;
1128 }
1129 /* z2[imin..imax] := f.z1[imin..imax-1] shifted right sh bits
1130 * (feeding f from the left). Assume sh > 0 */
1131 INLINE void
shift_right2(GEN z2,GEN z1,long imin,long imax,ulong f,ulong sh,ulong m)1132 shift_right2(GEN z2, GEN z1, long imin, long imax, ulong f, ulong sh, ulong m)
1133 {
1134 GEN sb = z1 + imin, se = z1 + imax, tb = z2 + imin;
1135 ulong k, l = *sb++;
1136 *tb++ = (l >> sh) | (f << m);
1137 while (sb < se) {
1138 k = l << m;
1139 l = *sb++;
1140 *tb++ = (l >> sh) | k;
1141 }
1142 }
1143
1144 /* Backward compatibility. Inefficient && unused */
1145 extern ulong hiremainder;
1146 INLINE ulong
shiftl(ulong x,ulong y)1147 shiftl(ulong x, ulong y)
1148 {
1149 hiremainder=x>>(BITS_IN_LONG-y);
1150 return (x<<y);
1151 }
1152 INLINE ulong
shiftlr(ulong x,ulong y)1153 shiftlr(ulong x, ulong y)
1154 {
1155 hiremainder=x<<(BITS_IN_LONG-y);
1156 return (x>>y);
1157 }
1158 #endif
1159