Annotation of OpenXM_contrib/pari-2.2/src/basemath/arith2.c, Revision 1.2
1.2 ! noro 1: /* $Id: arith2.c,v 1.53 2002/09/02 21:14:43 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*********************************************************************/
17: /** **/
18: /** ARITHMETIC FUNCTIONS **/
19: /** (second part) **/
20: /** **/
21: /*********************************************************************/
22: #include "pari.h"
23:
24: extern GEN arith_proto(long f(GEN), GEN x, int do_error);
25: extern GEN arith_proto2(long f(GEN,GEN), GEN x, GEN n);
26: extern GEN garith_proto(GEN f(GEN), GEN x, int do_error);
27: extern GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y);
1.2 ! noro 28: extern GEN garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z);
1.1 noro 29:
30: #define sqru(i) (muluu((i),(i)))
31:
32: /***********************************************************************/
33: /** **/
34: /** PRIME NUMBERS **/
35: /** **/
36: /***********************************************************************/
37:
38: GEN
39: prime(long n)
40: {
41: byteptr p = diffptr;
1.2 ! noro 42: long prime = 0;
1.1 noro 43:
44: if (n <= 0) err(talker, "n-th prime meaningless if n = %ld",n);
45: while (n--)
46: {
1.2 ! noro 47: NEXT_PRIME_VIADIFF_CHECK(prime,p);
1.1 noro 48: }
49: return stoi(prime);
50: }
51:
52: GEN
53: pith(long n)
54: {
55: byteptr p = diffptr;
56: ulong prime = 0, res = 0;
57:
58: if (n <= 0) err(talker, "pith meaningless if n = %ld",n);
1.2 ! noro 59: if (maxprime() <= (ulong)n) err(primer1);
! 60: while (prime<=(ulong)n) {
! 61: NEXT_PRIME_VIADIFF(prime,p);
! 62: res++;
! 63: }
1.1 noro 64: return stoi(res-1);
65: }
66:
67: GEN
68: primes(long n)
69: {
70: byteptr p = diffptr;
1.2 ! noro 71: long prime = 0;
1.1 noro 72: GEN y,z;
73:
74: if (n < 0) n = 0;
75: z = y = cgetg(n+1,t_VEC);
76: while (n--)
77: {
1.2 ! noro 78: NEXT_PRIME_VIADIFF_CHECK(prime,p);
! 79: *++z = lstoi(prime);
1.1 noro 80: }
81: return y;
82: }
83:
84: /* Building/Rebuilding the diffptr table. Incorporates Ilya Zakharevich's
85: * block sweep algorithm (see pari-dev-111, 25 April 1998). The actual work
86: * is done by the following two subroutines; the user entry point is the
87: * function initprimes() below. initprimes1() is the old algorithm, called
88: * when maxnum (size) is moderate.
89: */
90: static byteptr
1.2 ! noro 91: initprimes1(ulong size, long *lenp, long *lastp)
1.1 noro 92: {
93: long k;
94: byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2);
95:
96: memset(p, 0, size + 2); fin = p + size;
97: for (r=q=p,k=1; r<=fin; )
98: {
99: do { r+=k; k+=2; r+=k; } while (*++q);
100: for (s=r; s<=fin; s+=k) *s=1;
101: }
102: r=p; *r++=2; *r++=1; /* 2 and 3 */
103: for (s=q=r-1; ; s=q)
104: {
105: do q++; while (*q);
106: if (q>fin) break;
107: *r++ = (q-s) << 1;
108: }
109: *r++=0;
110: *lenp = r - p;
111: *lastp = ((s - p) << 1) + 1;
112: #if 0
113: fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n",
114: q, *lenp, *lastp);
115: #endif
1.2 ! noro 116: return (byteptr) gprealloc(p,r-p);
1.1 noro 117: }
118:
119: #define PRIME_ARENA (200 * 1024) /* No slowdown even with 256K level-2 cache */
120:
121: /* Here's the workhorse. This is recursive, although normally the first
122: recursive call will bottom out and invoke initprimes1() at once.
123: (Not static; might conceivably be useful to someone in library mode) */
124: byteptr
1.2 ! noro 125: initprimes0(ulong maxnum, long *lenp, ulong *lastp)
1.1 noro 126: {
1.2 ! noro 127: long k, size, alloced, asize, psize, rootnum;
1.1 noro 128: byteptr q,r,s,fin, p, p1, fin1, plast, curdiff;
1.2 ! noro 129: ulong last, remains, curlow;
1.1 noro 130:
131: #if 0
132: fprintferr("initprimes0: called for maxnum = %lu\n", maxnum);
133: #endif
134: if (maxnum <= 1ul<<17) /* Arbitrary. */
135: return initprimes1(maxnum>>1, lenp, lastp);
136:
137: maxnum |= 1; /* make it odd. */
138:
139: /* Checked to be enough up to 40e6, attained at 155893 */
140: size = (long) (1.09 * maxnum/log((double)maxnum)) + 145;
141: p1 = (byteptr) gpmalloc(size);
142: rootnum = (long) sqrt((double)maxnum); /* cast it back to a long */
143: rootnum |= 1;
144: #if 0
145: fprintferr("initprimes0: rootnum = %ld\n", rootnum);
146: #endif
147: {
148: byteptr p2 = initprimes0(rootnum, &psize, &last); /* recursive call */
149: memcpy(p1, p2, psize); free(p2);
150: }
151: fin1 = p1 + psize - 1;
152: remains = (maxnum - rootnum) >> 1; /* number of odd numbers to check */
153:
154: /* ARENA_IN_ROOTS below 12: some slowdown starts to be noticable
155: * when things fit into the cache.
156: * XXX The choice of 10 gives a slowdown of 1-2% on UltraSparcII,
157: * but makes calculations even for (the current) maximum of 436273009
158: * fit into 256K cache (still common for some architectures).
159: *
160: * One may change it when small caches become uncommon, but the gain
161: * is not going to be very noticable... */
162: #define ARENA_IN_ROOTS 10
163:
164: asize = ARENA_IN_ROOTS * rootnum; /* Make % overhead negligeable. */
165: if (asize < PRIME_ARENA) asize = PRIME_ARENA;
166: if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */
1.2 ! noro 167: alloced = (avma - bot < (size_t)asize>>1); /* enough room on the stack ? */
1.1 noro 168: if (alloced)
169: p = (byteptr) gpmalloc(asize);
170: else
171: p = (byteptr) bot;
172: fin = p + asize - 1; /* the 0 sentinel goes at fin. */
173: curlow = rootnum + 2; /* We know all primes up to rootnum (odd). */
174: curdiff = fin1;
175:
176: /* During each iteration p..fin-1 represents a range of odd
177: numbers. plast is a pointer which represents the last prime seen,
178: it may point before p..fin-1. */
179: plast = p - ((rootnum - last) >> 1) - 1;
180: while (remains)
181: {
182: if (asize > remains)
183: {
184: asize = remains + 1;
185: fin = p + asize - 1;
186: }
187: memset(p, 0, asize);
188: /* p corresponds to curlow. q runs over primediffs. */
1.2 ! noro 189: /* Don't care about DIFFPTR_SKIP: false positives provide no problem */
1.1 noro 190: for (q = p1+2, k = 3; q <= fin1; k += *q++)
191: {
192: /* The first odd number which is >= curlow and divisible by p
193: equals (if curlow > p)
194: the last odd number which is <= curlow + 2p - 1 and 0 (mod p)
195: which equals
196: p + the last even number which is <= curlow + p - 1 and 0 (mod p)
197: which equals
198: p + the last even number which is <= curlow + p - 2 and 0 (mod p)
199: which equals
200: p + curlow + p - 2 - (curlow + p - 2)) % 2p.
201: */
202: long k2 = k*k - curlow;
203:
204: if (k2 > 0)
205: {
206: r = p + (k2 >>= 1);
207: if (k2 > remains) break; /* Guard against an address wrap. */
208: } else
209: r = p - (((curlow+k-2) % (2*k)) >> 1) + k - 1;
210: for (s = r; s < fin; s += k) *s = 1;
211: }
212: /* now q runs over addresses corresponding to primes */
213: for (q = p; ; plast = q++)
214: {
1.2 ! noro 215: long d;
! 216:
1.1 noro 217: while (*q) q++; /* use 0-sentinel at end */
218: if (q >= fin) break;
1.2 ! noro 219: d = (q - plast) << 1;
! 220: while (d >= DIFFPTR_SKIP)
! 221: *curdiff++ = DIFFPTR_SKIP, d -= DIFFPTR_SKIP;
! 222: *curdiff++ = d;
1.1 noro 223: }
224: plast -= asize - 1;
225: remains -= asize - 1;
226: curlow += ((asize - 1)<<1);
227: } /* while (remains) */
228: last = curlow - ((p - plast) << 1);
229: *curdiff++ = 0; /* sentinel */
230: *lenp = curdiff - p1;
231: *lastp = last;
232: if (alloced) free(p);
1.2 ! noro 233: return (byteptr) gprealloc(p1, *lenp);
1.1 noro 234: }
235: #if 0 /* not yet... GN */
236: /* The diffptr table will contain at least 6548 entries (up to and including
237: the 6547th prime, 65557, and the terminal null byte), because the isprime/
238: small-factor-extraction machinery wants to depend on everything up to 65539
239: being in the table, and we might as well go to a multiple of 4 Bytes.--GN */
240:
241: void
242: init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */
243: #endif
244:
245: static ulong _maxprime = 0;
246:
247: ulong
1.2 ! noro 248: maxprime(void) { return _maxprime; }
1.1 noro 249:
250: byteptr
1.2 ! noro 251: initprimes(ulong maxnum)
1.1 noro 252: {
1.2 ! noro 253: long len;
! 254: ulong last;
1.1 noro 255: byteptr p;
1.2 ! noro 256: /* The algorithm must see the next prime beyond maxnum, whence the +512. */
! 257: /* switch to unsigned: adding the 512 _could_ wrap to a negative number. */
! 258: ulong maxnum1 = ((maxnum<65302?65302:maxnum)+512ul);
1.1 noro 259:
1.2 ! noro 260: if ((maxnum>>1) > VERYBIGINT - 1024)
! 261: err(talker, "Too large primelimit");
1.1 noro 262: p = initprimes0(maxnum1, &len, &last);
263: #if 0 /* not yet... GN */
264: static int build_the_tables = 1;
265: if (build_the_tables) { init_tinyprimes_tridiv(p); build_the_tables=0; }
266: #endif
267: _maxprime = last; return p;
268: }
269:
270: static void
271: cleanprimetab(void)
272: {
273: long i,j, lp = lg(primetab);
274:
275: for (i=j=1; i < lp; i++)
276: if (primetab[i]) primetab[j++] = primetab[i];
277: setlg(primetab,j);
278: }
279:
280: /* p is a single element or a row vector with "primes" to add to primetab.
281: * If p shares a non-trivial factor with another element in primetab, take it
282: * into account. */
283: GEN
284: addprimes(GEN p)
285: {
1.2 ! noro 286: gpmem_t av;
1.1 noro 287: long i,k,tx,lp;
288: GEN L;
289:
290: if (!p) return primetab;
291: tx = typ(p);
292: if (is_vec_t(tx))
293: {
1.2 ! noro 294: for (i=1; i < lg(p); i++) (void)addprimes((GEN) p[i]);
1.1 noro 295: return primetab;
296: }
297: if (tx != t_INT) err(typeer,"addprime");
298: if (is_pm1(p)) return primetab;
299: av = avma; i = signe(p);
300: if (i == 0) err(talker,"can't accept 0 in addprimes");
301: if (i < 0) p = absi(p);
302:
303: lp = lg(primetab);
304: L = cgetg(2*lp,t_VEC); k = 1;
305: for (i=1; i < lp; i++)
306: {
307: GEN n = (GEN)primetab[i], d = mppgcd(n, p);
308: if (! is_pm1(d))
309: {
310: if (!egalii(p,d)) L[k++] = (long)d;
311: L[k++] = ldivii(n,d);
312: gunclone(n); primetab[i] = 0;
313: }
314: }
1.2 ! noro 315: primetab = (GEN) gprealloc(primetab, (lp+1)*sizeof(long));
1.1 noro 316: primetab[i] = lclone(p); setlg(primetab, lp+1);
1.2 ! noro 317: if (k > 1) { cleanprimetab(); setlg(L,k); (void)addprimes(L); }
1.1 noro 318: avma = av; return primetab;
319: }
320:
321: static GEN
322: removeprime(GEN prime)
323: {
324: long i;
325:
326: if (typ(prime) != t_INT) err(typeer,"removeprime");
327: for (i=lg(primetab) - 1; i; i--)
328: if (absi_equal((GEN) primetab[i], prime))
329: {
330: gunclone((GEN)primetab[i]); primetab[i]=0;
331: cleanprimetab(); break;
332: }
333: if (!i) err(talker,"prime %Z is not in primetable", prime);
334: return primetab;
335: }
336:
337: GEN
338: removeprimes(GEN prime)
339: {
340: long i,tx;
341:
342: if (!prime) return primetab;
343: tx = typ(prime);
344: if (is_vec_t(tx))
345: {
346: if (prime == primetab)
347: {
348: for (i=1; i < lg(prime); i++) gunclone((GEN)prime[i]);
349: setlg(prime, 1);
350: }
351: else
352: {
1.2 ! noro 353: for (i=1; i < lg(prime); i++) (void)removeprime((GEN) prime[i]);
1.1 noro 354: }
355: return primetab;
356: }
357: return removeprime(prime);
358: }
359:
360: /***********************************************************************/
361: /** **/
362: /** COMPUTING THE MATRIX OF PRIME DIVISORS AND EXPONENTS **/
363: /** **/
364: /***********************************************************************/
365:
366: /* Configure may/should define this to 1 if MPQS is not installed --GN */
367: #ifndef decomp_default_hint
368: #define decomp_default_hint 0
369: #endif
370:
371: /* Some overkill removed from this (15 spsp for an integer < 2^32 ?!).
372: Should really revert to isprime() once the new primality tester is ready
373: --GN */
1.2 ! noro 374: #if 0
1.1 noro 375: #define pseudoprime(p) millerrabin(p,3*lgefint(p))
1.2 ! noro 376: #else /* remove further overkill :-) --KB */
! 377: #define pseudoprime(p) BSW_psp(p)
! 378: #endif
1.1 noro 379:
380: /* where to stop trial dividing in factorization */
381:
382: static long
383: tridiv_bound(GEN n, long all)
384: {
385: long size = expi(n) + 1;
386: if (all > 1) /* bounded factoring */
387: return all; /* use the given limit */
388: if (all == 0)
389: return VERYBIGINT; /* smallfact() case */
390: if (size <= 32)
391: return 16384;
392: else if (size <= 512)
393: return (size-16) << 10;
394: return 1L<<19; /* Rho will generally be faster above this */
395: }
396:
397: /* function imported from ifactor1.c */
398: extern long ifac_decomp(GEN n, long hint);
399: extern long ifac_decomp_break(GEN n, long (*ifac_break)(GEN,GEN,GEN,GEN),
400: GEN state, long hint);
401: static GEN
402: aux_end(GEN n, long nb)
403: {
404: GEN p1,p2, z = (GEN)avma;
405: long i;
406:
407: if (n) gunclone(n);
408: p1=cgetg(nb+1,t_COL);
409: p2=cgetg(nb+1,t_COL);
410: for (i=nb; i; i--)
411: {
412: p2[i] = (long)z; z += lg(z);
413: p1[i] = (long)z; z += lg(z);
414: }
415: z[1] = (long)p1;
416: z[2] = (long)p2;
417: if (nb) (void)sort_factor_gen(z,cmpii);
418: return z;
419: }
420:
421: GEN
422: auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state),
423: GEN state, long all, long hint)
424: {
1.2 ! noro 425: gpmem_t av;
! 426: long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 };
! 427: long nb = 0,i,k,lp,p,lim1;
1.1 noro 428: byteptr d=diffptr+1;
429:
430: if (typ(n) != t_INT) err(arither1);
431: i=signe(n); if (!i) err(arither2);
1.2 ! noro 432: (void)cgetg(3,t_MAT);
! 433: if (i<0) { (void)stoi(-1); (void)stoi(1); nb++; }
1.1 noro 434: if (is_pm1(n)) return aux_end(NULL,nb);
435:
436: n = gclone(n); setsigne(n,1);
437: i = vali(n);
438: if (i)
439: {
1.2 ! noro 440: (void)stoi(2); (void)stoi(i); nb++;
1.1 noro 441: av=avma; affii(shifti(n,-i), n); avma=av;
442: }
443: if (is_pm1(n)) return aux_end(n,nb);
444: lim1 = tridiv_bound(n, all);
445:
446: /* trial division */
1.2 ! noro 447: p = 2;
! 448: while (*d && p <= lim1)
1.1 noro 449: {
1.2 ! noro 450: NEXT_PRIME_VIADIFF(p,d);
! 451: if (mpdivisis(n,p,n))
1.1 noro 452: {
1.2 ! noro 453: nb++; k=1; while (mpdivisis(n,p,n)) k++;
! 454: (void)utoi(p); (void)stoi(k);
1.1 noro 455: if (is_pm1(n)) return aux_end(n,nb);
456: }
457: }
458:
459: /* pp = square of biggest p tried so far */
1.2 ! noro 460: av=avma; affii(muluu(p,p),pp); avma=av;
1.1 noro 461: if (cmpii(pp,n) > 0)
462: {
463: nb++;
1.2 ! noro 464: (void)icopy(n); (void)stoi(1); return aux_end(n,nb);
1.1 noro 465: }
466:
467: /* trial divide by the "special primes" (usually huge composites...) */
468: lp = lg(primetab); k=0;
469: for (i=1; i<lp; i++)
470: if (mpdivis(n,(GEN) primetab[i],n))
471: {
472: nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++;
1.2 ! noro 473: (void)icopy((GEN) primetab[i]); (void)stoi(k);
1.1 noro 474: if (is_pm1(n)) return aux_end(n,nb);
475: }
476:
477: /* test primality (unless all != 1 (i.e smallfact))*/
478: if ((k && cmpii(pp,n) > 0) || (all==1 && pseudoprime(n)))
479: {
480: nb++;
1.2 ! noro 481: (void)icopy(n); (void)stoi(1); return aux_end(n,nb);
1.1 noro 482: }
483:
484: /* now we have a large composite. Use hint as is if all==1 */
485: if (all!=1) hint = 15; /* turn off everything except pure powers */
486: if (ifac_break && (*ifac_break)(n,NULL,NULL,state))
487: /*Initialize ifac_break*/
488: {
489: if (DEBUGLEVEL >= 3)
490: fprintferr("IFAC: (Partial fact.) Initial stop requested.\n");
491: }
492: else
493: nb += ifac_decomp_break(n, ifac_break, state, hint);
494:
495: return aux_end(n, nb);
496: }
497:
498: /* state[1]: current unfactored part.
499: * state[2]: limit. */
500: static long
501: ifac_break_limit(GEN n, GEN pairs/*unused*/, GEN here, GEN state)
502: {
1.2 ! noro 503: gpmem_t ltop = avma;
1.1 noro 504: int res;
505: if (!here) /* initial call */
506: /* Small prime have been removed since start, n is the new unfactored part.
507: * Result is affect()ed to state[1] to preserve stack. */
508: affii(n, (GEN)state[1]);
509: else
510: {
511: GEN q = powgi((GEN)here[0],(GEN)here[1]); /* primary factor found.*/
512: if (DEBUGLEVEL>2) fprintferr("IFAC: Stop: Primary factor: %Z\n",q);
513: /* divide unfactored part by q and assign the result to state[1] */
514: diviiz((GEN)state[1],q, (GEN)state[1]);
515: }
516: if (DEBUGLEVEL>=3) fprintferr("IFAC: Stop: remaining %Z\n",state[1]);
517: /* check the stopping criterion, then restore stack */
518: res = cmpii((GEN)state[1],(GEN)state[2]) <= 0;
519: avma = ltop; return res;
520: }
521:
522: static GEN
523: auxdecomp0(GEN n, long all, long hint)
524: {
525: return auxdecomp1(n, NULL, gzero, all, hint);
526: }
527:
528: GEN
529: auxdecomp(GEN n, long all)
530: {
531: return auxdecomp0(n,all,decomp_default_hint);
532: }
533:
534: /* see before ifac_crack() in ifactor1.c for current semantics of `hint'
535: (factorint's `flag') */
536: GEN
537: factorint(GEN n, long flag)
538: {
539: return auxdecomp0(n,1,flag);
540: }
541:
542: GEN
543: decomp(GEN n)
544: {
545: return auxdecomp0(n,1,decomp_default_hint);
546: }
547:
548: /* Factor until the unfactored part is smaller than limit. */
549: GEN
550: decomp_limit(GEN n, GEN limit)
551: {
552: GEN state = cgetg(3,t_VEC);
553: /* icopy is mainly done to allocate memory for affect().
554: * Currently state[1] value is discarded in initial call to ifac_break_limit */
555: state[1] = licopy(n);
556: state[2] = lcopy(limit);
557: return auxdecomp1(n, &ifac_break_limit, state, 1, decomp_default_hint);
558: }
559:
560: GEN
561: smallfact(GEN n)
562: {
563: return boundfact(n,0);
564: }
565:
566: GEN
567: gboundfact(GEN n, long lim)
568: {
569: return garith_proto2gs(boundfact,n,lim);
570: }
571:
572: GEN
573: boundfact(GEN n, long lim)
574: {
575: GEN p1,p2,p3,p4,p5,y;
1.2 ! noro 576: gpmem_t av = avma,tetpil;
1.1 noro 577:
578: if (lim<=1) lim=0;
579: switch(typ(n))
580: {
581: case t_INT:
582: return auxdecomp(n,lim);
583: case t_FRACN:
1.2 ! noro 584: n=gred(n);
! 585: if (typ(n) == t_INT)
! 586: return gerepileupto(av, boundfact(n,lim));
! 587: /* fall through */
1.1 noro 588: case t_FRAC:
589: p1=auxdecomp((GEN)n[1],lim);
590: p2=auxdecomp((GEN)n[2],lim);
591: p4=concatsp((GEN)p1[1],(GEN)p2[1]);
592: p5=concatsp((GEN)p1[2],gneg((GEN)p2[2]));
593: p3=indexsort(p4);
594:
595: tetpil=avma; y=cgetg(3,t_MAT);
596: y[1]=(long)extract(p4,p3);
597: y[2]=(long)extract(p5,p3);
598: return gerepile(av,tetpil,y);
599: }
600: err(arither1);
601: return NULL; /* not reached */
602: }
1.2 ! noro 603: /***********************************************************************/
! 604: /** **/
! 605: /** SIMPLE FACTORISATIONS **/
! 606: /** **/
! 607: /***********************************************************************/
! 608:
! 609: /*Factorize n and output [fp,fe] where
! 610: * fp and fe are vecsmall and n=prod{fp[i]^fe[i]}
! 611: */
! 612:
! 613: GEN
! 614: decomp_small(long n)
! 615: {
! 616: gpmem_t ltop=avma;
! 617: GEN F = factor(stoi(n));
! 618: long i, l=lg(F[1]);
! 619: GEN f=cgetg(3,t_VEC);
! 620: GEN fp=cgetg(l,t_VECSMALL);
! 621: GEN fe=cgetg(l,t_VECSMALL);
! 622: f[1] = (long) fp;
! 623: f[2] = (long) fe;
! 624: for(i = 1; i < l; i++)
! 625: {
! 626: fp[i]=itos(gcoeff(F,i,1));
! 627: fe[i]=itos(gcoeff(F,i,2));
! 628: }
! 629: return gerepileupto(ltop,f);
! 630: }
! 631:
! 632: /*Return the primary factors of a small integer as a vecsmall*/
! 633: GEN
! 634: decomp_primary_small(long n)
! 635: {
! 636: gpmem_t ltop=avma;
! 637: GEN F = factor(stoi(n));
! 638: GEN fc = cgetg(lg(F[1]), t_VECSMALL);
! 639: gpmem_t av=avma;
! 640: long i;
! 641: for (i = 1; i < lg(fc); i++)
! 642: fc[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i)));
! 643: avma=av;
! 644: return gerepileupto(ltop,fc);
! 645: }
! 646:
1.1 noro 647:
648: /***********************************************************************/
649: /** **/
650: /** BASIC ARITHMETIC FUNCTIONS **/
651: /** **/
652: /***********************************************************************/
653:
654: /* functions imported from ifactor1.c */
655: long
656: ifac_moebius(GEN n, long hint);
657:
658: long
659: ifac_issquarefree(GEN n, long hint);
660:
661: long
662: ifac_omega(GEN n, long hint);
663:
664: long
665: ifac_bigomega(GEN n, long hint);
666:
667: GEN
668: ifac_totient(GEN n, long hint);
669:
670: GEN
671: ifac_numdiv(GEN n, long hint);
672:
673: GEN
674: ifac_sumdiv(GEN n, long hint);
675:
676: GEN
677: ifac_sumdivk(GEN n, long k, long hint);
678:
679: GEN
680: gmu(GEN n)
681: {
682: return arith_proto(mu,n,1);
683: }
684:
685: long
686: mu(GEN n)
687: {
688: byteptr d = diffptr+1; /* point at 3 - 2 */
1.2 ! noro 689: gpmem_t av = avma;
! 690: ulong p;
! 691: long s, v, lim1;
1.1 noro 692:
693: if (typ(n) != t_INT) err(arither1);
694: if (!signe(n)) err(arither2);
695: if (is_pm1(n)) return 1;
696: v = vali(n);
697: if (v>1) return 0;
698: s = v ? -1 : 1;
699: n=absi(shifti(n,-v)); p = 2;
700: if (is_pm1(n)) return s;
701: lim1 = tridiv_bound(n,1);
702:
703: while (*d && p < lim1)
704: {
1.2 ! noro 705: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 706: if (mpdivisis(n,p,n))
707: {
708: if (smodis(n,p) == 0) { avma=av; return 0; }
709: s = -s; if (is_pm1(n)) { avma=av; return s; }
710: }
711: }
712: /* we normally get here with p==nextprime(PRE_TUNE), which will already
713: have been tested against n, so p^2 >= n (instead of p^2 > n) is
714: enough to guarantee n is prime */
715: if (cmpii(sqru(p),n) >= 0 || pseudoprime(n)) { avma=av; return -s; }
716: /* large composite without small factors */
717: v = ifac_moebius(n, decomp_default_hint);
718: avma=av;
719: return (s<0 ? -v : v); /* correct also if v==0 */
720: }
721:
722: GEN
723: gissquarefree(GEN x)
724: {
725: return arith_proto(issquarefree,x,0);
726: }
727:
728: long
729: issquarefree(GEN x)
730: {
1.2 ! noro 731: ulong p;
! 732: gpmem_t av = avma;
! 733: long tx, lim1;
1.1 noro 734: GEN d;
735:
736: if (gcmp0(x)) return 0;
737: tx = typ(x);
738: if (tx == t_INT)
739: {
740: long v;
741: byteptr d=diffptr+1;
742: if (is_pm1(x)) return 1;
743: if((v = vali(x)) > 1) return 0;
744: x = absi(shifti(x,-v));
745: if (is_pm1(x)) return 1;
746: lim1 = tridiv_bound(x,1);
747:
748: p = 2;
749: while (*d && p < lim1)
750: {
1.2 ! noro 751: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 752: if (mpdivisis(x,p,x))
753: {
754: if (smodis(x,p) == 0) { avma = av; return 0; }
755: if (is_pm1(x)) { avma = av; return 1; }
756: }
757: }
758: if (cmpii(sqru(p),x) >= 0 || pseudoprime(x)) { avma = av; return 1; }
759: v = ifac_issquarefree(x, decomp_default_hint);
760: avma = av; return v;
761: }
762: if (tx != t_POL) err(typeer,"issquarefree");
763: d = ggcd(x, derivpol(x));
764: avma = av; return (lgef(d) == 3);
765: }
766:
767: GEN
768: gomega(GEN n)
769: {
770: return arith_proto(omega,n,1);
771: }
772:
773: long
774: omega(GEN n)
775: {
776: byteptr d=diffptr+1;
1.2 ! noro 777: gpmem_t av = avma;
! 778: long p,nb,v, lim1;
1.1 noro 779:
780: if (typ(n) != t_INT) err(arither1);
781: if (!signe(n)) err(arither2);
782: if (is_pm1(n)) return 0;
783: v=vali(n);
784: nb = v ? 1 : 0;
785: n = absi(shifti(n,-v)); p = 2;
786: if (is_pm1(n)) return nb;
787: lim1 = tridiv_bound(n,1);
788:
789: while (*d && p < lim1)
790: {
1.2 ! noro 791: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 792: if (mpdivisis(n,p,n))
793: {
794: nb++; while (mpdivisis(n,p,n)); /* empty */
795: if (is_pm1(n)) { avma = av; return nb; }
796: }
797: }
798: if (cmpii(sqru(p),n) >= 0 || pseudoprime(n)) { avma = av; return nb+1; }
799: /* large composite without small factors */
800: nb += ifac_omega(n, decomp_default_hint);
801: avma=av; return nb;
802: }
803:
804: GEN
805: gbigomega(GEN n)
806: {
807: return arith_proto(bigomega,n,1);
808: }
809:
810: long
811: bigomega(GEN n)
812: {
813: byteptr d=diffptr+1;
1.2 ! noro 814: ulong p;
! 815: gpmem_t av = avma;
! 816: long nb,v, lim1;
1.1 noro 817:
818: if (typ(n) != t_INT) err(arither1);
819: if (!signe(n)) err(arither2);
820: if (is_pm1(n)) return 0;
821: nb=v=vali(n);
822: n=absi(shifti(n,-v)); p = 2;
823: if (is_pm1(n)) { avma = av; return nb; }
824: lim1 = tridiv_bound(n,1);
825:
826: while (*d && p < lim1)
827: {
1.2 ! noro 828: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 829: if (mpdivisis(n,p,n))
830: {
831: do nb++; while (mpdivisis(n,p,n));
832: if (is_pm1(n)) { avma = av; return nb; }
833: }
834: }
835: if (cmpii(sqru(p),n) >= 0 || pseudoprime(n)) { avma = av; return nb+1; }
836: nb += ifac_bigomega(n, decomp_default_hint);
837: avma=av; return nb;
838: }
839:
840: GEN
841: gphi(GEN n)
842: {
843: return garith_proto(phi,n,1);
844: }
845:
846: GEN
847: phi(GEN n)
848: {
849: byteptr d = diffptr+1;
850: GEN m;
1.2 ! noro 851: ulong p;
! 852: gpmem_t av = avma;
! 853: long v, lim1;
1.1 noro 854:
855: if (typ(n) != t_INT) err(arither1);
856: if (!signe(n)) err(arither2);
857: if (is_pm1(n)) return gun;
858: v = vali(n);
859: n = absi(shifti(n,-v)); p = 2;
860: m = (v > 1 ? shifti(gun,v-1) : stoi(1));
861: if (is_pm1(n)) { return gerepileupto(av,m); }
862: lim1 = tridiv_bound(n,1);
863:
864: while (*d && p < lim1)
865: {
1.2 ! noro 866: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 867: if (mpdivisis(n,p,n))
868: {
869: m = mulis(m, p-1);
870: while (mpdivisis(n,p,n)) m = mulis(m,p);
871: if (is_pm1(n)) { return gerepileupto(av,m); }
872: }
873: }
874: if (cmpii(sqru(p),n) >= 0 || pseudoprime(n))
875: {
876: m = mulii(m,addsi(-1,n));
877: return gerepileupto(av,m);
878: }
879: m = mulii(m, ifac_totient(n, decomp_default_hint));
880: return gerepileupto(av,m);
881: }
882:
883: GEN
884: gnumbdiv(GEN n)
885: {
886: return garith_proto(numbdiv,n,1);
887: }
888:
889: GEN
890: numbdiv(GEN n)
891: {
892: byteptr d=diffptr+1;
893: GEN m;
1.2 ! noro 894: long l, v, lim1;
! 895: ulong p;
! 896: gpmem_t av = avma;
1.1 noro 897:
898: if (typ(n) != t_INT) err(arither1);
899: if (!signe(n)) err(arither2);
900: if (is_pm1(n)) return gun;
901: v = vali(n);
902: n = absi(shifti(n,-v)); p = 2;
903: m = stoi(v+1);
904: if (is_pm1(n)) return gerepileupto(av,m);
905: lim1 = tridiv_bound(n,1);
906:
907: while (*d && p < lim1)
908: {
1.2 ! noro 909: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 910: l=1; while (mpdivisis(n,p,n)) l++;
911: m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); }
912: }
913: if(cmpii(sqru(p),n) >= 0 || pseudoprime(n))
914: {
915: m = shifti(m,1);
916: return gerepileupto(av,m);
917: }
918: m = mulii(m, ifac_numdiv(n, decomp_default_hint));
919: return gerepileupto(av,m);
920: }
921:
922: GEN
923: gsumdiv(GEN n)
924: {
925: return garith_proto(sumdiv,n,1);
926: }
927:
928: GEN
929: sumdiv(GEN n)
930: {
931: byteptr d=diffptr+1;
932: GEN m,m1;
1.2 ! noro 933: ulong p;
! 934: gpmem_t av=avma;
! 935: long v, lim1;
1.1 noro 936:
937: if (typ(n) != t_INT) err(arither1);
938: if (!signe(n)) err(arither2);
939: if (is_pm1(n)) return gun;
940: v = vali(n);
941: n = absi(shifti(n,-v)); p = 2;
942: m = (v ? addsi(-1,shifti(gun,v+1)) : stoi(1));
943: if (is_pm1(n)) { return gerepileupto(av,m); }
944: lim1 = tridiv_bound(n,1);
945:
946: while (*d && p < lim1)
947: {
1.2 ! noro 948: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 949: if (mpdivisis(n,p,n))
950: {
951: m1 = utoi(p+1);
952: while (mpdivisis(n,p,n)) m1 = addsi(1, mului(p,m1));
953: m = mulii(m1,m); if (is_pm1(n)) { return gerepileupto(av,m); }
954: }
955: }
956: if(cmpii(sqru(p),n) >= 0 || pseudoprime(n))
957: {
958: m = mulii(m,addsi(1,n));
959: return gerepileupto(av,m);
960: }
961: m = mulii(m, ifac_sumdiv(n, decomp_default_hint));
962: return gerepileupto(av,m);
963: }
964:
965: GEN
966: gsumdivk(GEN n, long k)
967: {
968: return garith_proto2gs(sumdivk,n,k);
969: }
970:
971: GEN
972: sumdivk(GEN n, long k)
973: {
974: byteptr d=diffptr+1;
975: GEN n1,m,m1,pk;
1.2 ! noro 976: ulong p;
! 977: gpmem_t av = avma;
! 978: long k1,v, lim1;
1.1 noro 979:
980: if (!k) return numbdiv(n);
981: if (k==1) return sumdiv(n);
982: if (typ(n) != t_INT) err(arither1);
983: if (!signe(n)) err(arither2);
984: if (is_pm1(n)) return gun;
985: k1 = k; n1 = n;
986: if (k==-1) { m=sumdiv(n); goto fin; }
987: if (k<0) k = -k;
988: v=vali(n);
989: n=absi(shifti(n,-v));
990: p = 2; m = stoi(1);
991: while (v--) m = addsi(1,shifti(m,k));
992: if (is_pm1(n)) goto fin;
993: lim1 = tridiv_bound(n,1);
994:
995: while (*d && p < lim1)
996: {
1.2 ! noro 997: NEXT_PRIME_VIADIFF(p,d);
1.1 noro 998: if (mpdivisis(n,p,n))
999: {
1000: pk = gpowgs(utoi(p),k); m1 = addsi(1,pk);
1001: while (mpdivisis(n,p,n)) m1 = addsi(1, mulii(pk,m1));
1002: m = mulii(m1,m); if (is_pm1(n)) goto fin;
1003: }
1004: }
1005: if(cmpii(sqru(p),n) >= 0 || pseudoprime(n))
1006: {
1007: pk = gpuigs(n,k);
1008: m = mulii(m,addsi(1,pk));
1009: goto fin;
1010: }
1011: m = mulii(m, ifac_sumdivk(n, k, decomp_default_hint));
1012: fin:
1013: if (k1 < 0) m = gdiv(m, gpowgs(n1,k));
1014: return gerepileupto(av,m);
1015: }
1016:
1017: /***********************************************************************/
1018: /** **/
1019: /** MISCELLANEOUS ARITHMETIC FUNCTIONS **/
1020: /** (all of these ultimately depend on auxdecomp()) **/
1021: /** **/
1022: /***********************************************************************/
1023:
1024: GEN
1025: divisors(GEN n)
1026: {
1.2 ! noro 1027: gpmem_t tetpil,av=avma;
! 1028: long i,j,l;
1.1 noro 1029: GEN *d,*t,*t1,*t2,*t3, nbdiv,e;
1030:
1031: if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1);
1032:
1033: e = (GEN) n[2], n = (GEN) n[1]; l = lg(n);
1034: if (l>1 && signe(n[1]) < 0) { e++; n++; l--; } /* skip -1 */
1035: nbdiv = gun;
1036: for (i=1; i<l; i++)
1037: {
1038: e[i] = itos((GEN)e[i]);
1039: nbdiv = mulis(nbdiv,1+e[i]);
1040: }
1041: if (is_bigint(nbdiv) || (itos(nbdiv) & ~LGBITS))
1042: err(talker, "too many divisors (more than %ld)", LGBITS - 1);
1043: d = t = (GEN*) cgetg(itos(nbdiv)+1,t_VEC);
1044: *++d = gun;
1045: for (i=1; i<l; i++)
1046: for (t1=t,j=e[i]; j; j--,t1=t2)
1047: for (t2=d,t3=t1; t3<t2; )
1048: *++d = mulii(*++t3, (GEN)n[i]);
1049: tetpil=avma; return gerepile(av,tetpil,sort((GEN)t));
1050: }
1051:
1.2 ! noro 1052: static GEN
! 1053: _core(GEN n, int all)
1.1 noro 1054: {
1.2 ! noro 1055: gpmem_t av = avma;
! 1056: long i;
! 1057: GEN fa,p1,p2,c = gun;
1.1 noro 1058:
1.2 ! noro 1059: fa = auxdecomp(n,all);
! 1060: p1 = (GEN)fa[1];
! 1061: p2 = (GEN)fa[2];
1.1 noro 1062: for (i=1; i<lg(p1); i++)
1.2 ! noro 1063: if (mod2((GEN)p2[i])) c = mulii(c,(GEN)p1[i]);
! 1064: return gerepileuptoint(av, c);
1.1 noro 1065: }
1066:
1.2 ! noro 1067: static GEN
! 1068: _core2(GEN n, int all)
1.1 noro 1069: {
1.2 ! noro 1070: gpmem_t av = avma;
! 1071: long i;
! 1072: GEN fa,p1,p2,e,c=gun,f=gun,y;
1.1 noro 1073:
1.2 ! noro 1074: fa = auxdecomp(n,all);
! 1075: p1 = (GEN)fa[1];
! 1076: p2 = (GEN)fa[2];
1.1 noro 1077: for (i=1; i<lg(p1); i++)
1078: {
1.2 ! noro 1079: e = (GEN)p2[i];
! 1080: if (mod2(e)) c = mulii(c, (GEN)p1[i]);
! 1081: if (!gcmp1(e)) f = mulii(f, powgi((GEN)p1[i], shifti(e,-1)));
1.1 noro 1082: }
1.2 ! noro 1083: y = cgetg(3,t_VEC);
! 1084: y[1] = (long)c;
! 1085: y[2] = (long)f; return gerepilecopy(av, y);
1.1 noro 1086: }
1087:
1.2 ! noro 1088: GEN core(GEN n) { return _core(n,1); }
! 1089: GEN corepartial(GEN n) { return _core(n,0); }
! 1090: GEN core2(GEN n) { return _core2(n,1); }
! 1091: GEN core2partial(GEN n){ return _core2(n,0); }
! 1092:
1.1 noro 1093: GEN
1.2 ! noro 1094: core0(GEN n,long flag) { return flag? core2(n): core(n); }
! 1095:
! 1096: static long
! 1097: _mod4(GEN c) { long r = mod4(c); if (signe(c) < 0) r = 4-r; return r; }
1.1 noro 1098:
1099: GEN
1100: coredisc(GEN n)
1101: {
1.2 ! noro 1102: gpmem_t av = avma;
! 1103: GEN c = core(n);
! 1104: long r = _mod4(c);
! 1105: if (r==1 || r==4) return c;
! 1106: return gerepileuptoint(av, shifti(c,2));
1.1 noro 1107: }
1108:
1109: GEN
1110: coredisc2(GEN n)
1111: {
1.2 ! noro 1112: gpmem_t av = avma;
! 1113: GEN y = core2(n);
! 1114: GEN c = (GEN)y[1], f = (GEN)y[2];
! 1115: long r = _mod4(c);
! 1116: if (r==1 || r==4) return y;
! 1117: y = cgetg(3,t_VEC);
! 1118: y[1] = lshifti(c,2);
! 1119: y[2] = lmul2n(f,-1); return gerepileupto(av, y);
1.1 noro 1120: }
1121:
1122: GEN
1.2 ! noro 1123: coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); }
1.1 noro 1124:
1125: /***********************************************************************/
1126: /** **/
1127: /** BINARY QUADRATIC FORMS **/
1128: /** **/
1129: /***********************************************************************/
1130:
1131: GEN
1132: qf_disc(GEN x, GEN y, GEN z)
1133: {
1134: if (!y) { y=(GEN)x[2]; z=(GEN)x[3]; x=(GEN)x[1]; }
1135: return subii(sqri(y), shifti(mulii(x,z),2));
1136: }
1137:
1138: static GEN
1139: qf_create(GEN x, GEN y, GEN z, long s)
1140: {
1141: GEN t;
1142: if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb");
1143: if (!s)
1144: {
1.2 ! noro 1145: gpmem_t av=avma; s = signe(qf_disc(x,y,z)); avma=av;
1.1 noro 1146: if (!s) err(talker,"zero discriminant in Qfb");
1147: }
1148: t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI);
1149: t[1]=(long)icopy(x); t[2]=(long)icopy(y); t[3]=(long)icopy(z);
1150: return t;
1151: }
1152:
1153: GEN
1154: qfi(GEN x, GEN y, GEN z)
1155: {
1156: return qf_create(x,y,z,-1);
1157: }
1158:
1159: GEN
1160: qfr(GEN x, GEN y, GEN z, GEN d)
1161: {
1162: GEN t = qf_create(x,y,z,1);
1163: if (typ(d) != t_REAL)
1164: err(talker,"Shanks distance should be a t_REAL (in qfr)");
1165: t[4]=lrcopy(d); return t;
1166: }
1167:
1168: GEN
1169: Qfb0(GEN x, GEN y, GEN z, GEN d, long prec)
1170: {
1171: GEN t = qf_create(x,y,z,0);
1172: if (lg(t)==4) return t;
1173: if (!d) d = gzero;
1174: if (typ(d) == t_REAL)
1175: t[4] = lrcopy(d);
1176: else
1177: { t[4]=lgetr(prec); gaffect(d,(GEN)t[4]); }
1178: return t;
1179: }
1180:
1181: /* composition */
1182: static void
1183: sq_gen(GEN z, GEN x)
1184: {
1185: GEN d1,x2,y2,v1,v2,c3,m,p1,r;
1186:
1187: d1 = bezout((GEN)x[2],(GEN)x[1],&x2,&y2);
1188: if (gcmp1(d1)) v1 = v2 = (GEN)x[1];
1189: else
1190: {
1191: v1 = divii((GEN)x[1],d1);
1192: v2 = mulii(v1,mppgcd(d1,(GEN)x[3]));
1193: }
1194: m = mulii((GEN)x[3],x2); setsigne(m,-signe(m));
1195: r = modii(m,v2); p1 = mulii(v1,r);
1196: c3 = addii(mulii((GEN)x[3],d1), mulii(r,addii((GEN)x[2],p1)));
1197: z[1] = lmulii(v1,v2);
1198: z[2] = laddii((GEN)x[2], shifti(p1,1));
1199: z[3] = ldivii(c3,v2);
1200: }
1201:
1202: void
1203: comp_gen(GEN z,GEN x,GEN y)
1204: {
1205: GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,c3,m,p1,r;
1206:
1207: if (x == y) { sq_gen(z,x); return; }
1208: s = shifti(addii((GEN)x[2],(GEN)y[2]),-1);
1209: n = subii((GEN)y[2],s);
1210: d = bezout((GEN)y[1],(GEN)x[1],&y1,&x1);
1211: d1 = bezout(s,d,&x2,&y2);
1212: if (gcmp1(d1))
1213: {
1214: v1 = (GEN)x[1];
1215: v2 = (GEN)y[1];
1216: }
1217: else
1218: {
1219: v1 = divii((GEN)x[1],d1);
1220: v2 = divii((GEN)y[1],d1);
1221: v1 = mulii(v1, mppgcd(d1,mppgcd((GEN)x[3],mppgcd((GEN)y[3],n))));
1222: }
1223: m = addii(mulii(mulii(y1,y2),n), mulii((GEN)y[3],x2));
1224: setsigne(m,-signe(m));
1225: r = modii(m,v1); p1 = mulii(v2,r);
1226: c3 = addii(mulii((GEN)y[3],d1), mulii(r,addii((GEN)y[2],p1)));
1227: z[1] = lmulii(v1,v2);
1228: z[2] = laddii((GEN)y[2], shifti(p1,1));
1229: z[3] = ldivii(c3,v1);
1230: }
1231:
1232: static GEN
1233: compimag0(GEN x, GEN y, int raw)
1234: {
1.2 ! noro 1235: gpmem_t av = avma;
! 1236: long tx = typ(x);
1.1 noro 1237: GEN z;
1238:
1239: if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition");
1240: if (cmpii((GEN)x[1],(GEN)y[1]) > 0) { z=x; x=y; y=z; }
1241: z=cgetg(4,t_QFI); comp_gen(z,x,y);
1242: if (raw) return gerepilecopy(av,z);
1243: return gerepileupto(av, redimag(z));
1244: }
1245:
1246: static GEN
1247: compreal0(GEN x, GEN y, int raw)
1248: {
1.2 ! noro 1249: gpmem_t av = avma;
! 1250: long tx = typ(x);
1.1 noro 1251: GEN z;
1252:
1253: if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition");
1254: z=cgetg(5,t_QFR); comp_gen(z,x,y);
1255: z[4]=laddrr((GEN)x[4],(GEN)y[4]);
1256: if (raw) return gerepilecopy(av,z);
1257: return gerepileupto(av,redreal(z));
1258: }
1259:
1260: GEN
1261: compreal(GEN x, GEN y) { return compreal0(x,y,0); }
1262:
1263: GEN
1264: comprealraw(GEN x, GEN y) { return compreal0(x,y,1); }
1265:
1266: GEN
1267: compimag(GEN x, GEN y) { return compimag0(x,y,0); }
1268:
1269: GEN
1270: compimagraw(GEN x, GEN y) { return compimag0(x,y,1); }
1271:
1272: GEN
1273: compraw(GEN x, GEN y)
1274: {
1275: return (typ(x)==t_QFI)? compimagraw(x,y): comprealraw(x,y);
1276: }
1277:
1278: GEN
1279: sqcompimag0(GEN x, long raw)
1280: {
1.2 ! noro 1281: gpmem_t av = avma;
1.1 noro 1282: GEN z = cgetg(4,t_QFI);
1283:
1284: if (typ(x)!=t_QFI) err(typeer,"composition");
1285: sq_gen(z,x);
1286: if (raw) return gerepilecopy(av,z);
1287: return gerepileupto(av,redimag(z));
1288: }
1289:
1290: static GEN
1291: sqcompreal0(GEN x, long raw)
1292: {
1.2 ! noro 1293: gpmem_t av = avma;
1.1 noro 1294: GEN z = cgetg(5,t_QFR);
1295:
1296: if (typ(x)!=t_QFR) err(typeer,"composition");
1297: sq_gen(z,x); z[4]=lshiftr((GEN)x[4],1);
1298: if (raw) return gerepilecopy(av,z);
1299: return gerepileupto(av,redreal(z));
1300: }
1301:
1302: GEN
1303: sqcompreal(GEN x) { return sqcompreal0(x,0); }
1304:
1305: GEN
1306: sqcomprealraw(GEN x) { return sqcompreal0(x,1); }
1307:
1308: GEN
1309: sqcompimag(GEN x) { return sqcompimag0(x,0); }
1310:
1311: GEN
1312: sqcompimagraw(GEN x) { return sqcompimag0(x,1); }
1313:
1314: static GEN
1315: real_unit_form_by_disc(GEN D, long prec)
1316: {
1317: GEN y = cgetg(5,t_QFR), isqrtD;
1.2 ! noro 1318: gpmem_t av = avma;
1.1 noro 1319:
1320: if (typ(D) != t_INT || signe(D) <= 0) err(typeer,"real_unit_form_by_disc");
1321: switch(mod4(D))
1322: {
1323: case 2:
1324: case 3: err(funder2,"real_unit_form_by_disc");
1325: }
1326: y[1]=un; isqrtD = racine(D);
1327: /* we know that D and isqrtD are non-zero */
1328: if (mod2(D) != mod2(isqrtD))
1329: isqrtD = gerepileuptoint(av, addsi(-1,isqrtD));
1330: y[2] = (long)isqrtD; av = avma;
1.2 ! noro 1331: y[3] = lpileuptoint(av, shifti(subii(sqri(isqrtD),D),-2));
1.1 noro 1332: y[4] = (long)realzero(prec); return y;
1333: }
1334:
1335: GEN
1336: real_unit_form(GEN x)
1337: {
1.2 ! noro 1338: gpmem_t av = avma;
! 1339: long prec;
1.1 noro 1340: GEN D;
1341: if (typ(x) != t_QFR) err(typeer,"real_unit_form");
1342: prec = precision((GEN)x[4]);
1343: if (!prec) err(talker,"not a t_REAL in 4th component of a t_QFR");
1344: D = qf_disc(x,NULL,NULL);
1345: return gerepileupto(av, real_unit_form_by_disc(D,prec));
1346: }
1347:
1348: static GEN
1349: imag_unit_form_by_disc(GEN D)
1350: {
1351: GEN y = cgetg(4,t_QFI);
1352: long isodd;
1353:
1354: if (typ(D) != t_INT || signe(D) >= 0) err(typeer,"real_unit_form_by_disc");
1355: switch(4 - mod4(D))
1356: {
1357: case 2:
1358: case 3: err(funder2,"imag_unit_form_by_disc");
1359: }
1360: y[1] = un; isodd = mpodd(D);
1361: y[2] = isodd? un: zero;
1362: /* y[3] = (1-D) / 4 or -D / 4, whichever is an integer */
1363: y[3] = lshifti(D,-2); setsigne(y[3],1);
1364: if (isodd)
1365: {
1.2 ! noro 1366: gpmem_t av = avma;
! 1367: y[3] = lpileuptoint(av, addis((GEN)y[3],1));
1.1 noro 1368: }
1369: return y;
1370: }
1371:
1372: GEN
1373: imag_unit_form(GEN x)
1374: {
1375: GEN p1,p2, y = cgetg(4,t_QFI);
1.2 ! noro 1376: gpmem_t av;
1.1 noro 1377: if (typ(x) != t_QFI) err(typeer,"imag_unit_form");
1378: y[1] = un;
1379: y[2] = mpodd((GEN)x[2])? un: zero;
1380: av = avma; p1 = mulii((GEN)x[1],(GEN)x[3]);
1381: p2 = shifti(sqri((GEN)x[2]),-2);
1.2 ! noro 1382: y[3] = lpileuptoint(av, subii(p1,p2));
1.1 noro 1383: return y;
1384: }
1385:
1386: GEN
1387: powrealraw(GEN x, long n)
1388: {
1.2 ! noro 1389: gpmem_t av = avma;
! 1390: long m;
1.1 noro 1391: GEN y;
1392:
1393: if (typ(x) != t_QFR)
1394: err(talker,"not a real quadratic form in powrealraw");
1395: if (!n) return real_unit_form(x);
1396: if (n== 1) return gcopy(x);
1397: if (n==-1) return ginv(x);
1398:
1399: y = NULL; m=labs(n);
1400: for (; m>1; m>>=1)
1401: {
1402: if (m&1) y = y? comprealraw(y,x): x;
1403: x=sqcomprealraw(x);
1404: }
1405: y = y? comprealraw(y,x): x;
1406: if (n<0) y = ginv(y);
1407: return gerepileupto(av,y);
1408: }
1409:
1410: GEN
1411: powimagraw(GEN x, long n)
1412: {
1.2 ! noro 1413: gpmem_t av = avma;
! 1414: long m;
1.1 noro 1415: GEN y;
1416:
1417: if (typ(x) != t_QFI)
1418: err(talker,"not an imaginary quadratic form in powimag");
1419: if (!n) return imag_unit_form(x);
1420: if (n== 1) return gcopy(x);
1421: if (n==-1) return ginv(x);
1422:
1423: y = NULL; m=labs(n);
1424: for (; m>1; m>>=1)
1425: {
1426: if (m&1) y = y? compimagraw(y,x): x;
1427: x=sqcompimagraw(x);
1428: }
1429: y = y? compimagraw(y,x): x;
1430: if (n<0) y = ginv(y);
1431: return gerepileupto(av,y);
1432: }
1433:
1434: GEN
1435: powraw(GEN x, long n)
1436: {
1437: return (typ(x)==t_QFI)? powimagraw(x,n): powrealraw(x,n);
1438: }
1439:
1440: /* composition: Shanks' NUCOMP & NUDUPL */
1441: /* l = floor((|d|/4)^(1/4)) */
1442: GEN
1443: nucomp(GEN x, GEN y, GEN l)
1444: {
1.2 ! noro 1445: gpmem_t av=avma,tetpil;
! 1446: long cz;
1.1 noro 1447: GEN a,a1,a2,b,b2,d,d1,e,g,n,p1,p2,p3,q1,q2,q3,q4,s,t2,t3,u,u1,v,v1,v2,v3,z;
1448:
1449: if (x==y) return nudupl(x,l);
1450: if (typ(x) != t_QFI || typ(y) != t_QFI)
1451: err(talker,"not an imaginary quadratic form in nucomp");
1452:
1453: if (cmpii((GEN)x[1],(GEN)y[1]) < 0) { s=x; x=y; y=s; }
1454: s=shifti(addii((GEN)x[2],(GEN)y[2]),-1); n=subii((GEN)y[2],s);
1455: a1=(GEN)x[1]; a2=(GEN)y[1]; d=bezout(a2,a1,&u,&v);
1456: if (gcmp1(d)) { a=negi(gmul(u,n)); d1=d; }
1457: else
1458: if (divise(s,d))
1459: {
1460: a=negi(gmul(u,n)); d1=d; a1=divii(a1,d1);
1461: a2=divii(a2,d1); s=divii(s,d1);
1462: }
1463: else
1464: {
1465: d1=bezout(s,d,&u1,&v1);
1466: if (!gcmp1(d1))
1467: {
1468: a1=divii(a1,d1); a2=divii(a2,d1);
1469: s=divii(s,d1); d=divii(d,d1);
1470: }
1471: p1=resii((GEN)x[3],d); p2=resii((GEN)y[3],d);
1472: p3=modii(negi(mulii(u1,addii(mulii(u,p1),mulii(v,p2)))),d);
1473: a=subii(mulii(p3,divii(a1,d)),mulii(u,divii(n,d)));
1474: }
1475: a=modii(a,a1); p1=subii(a1,a); if (cmpii(a,p1)>0) a=negi(p1);
1476: v=gzero; d=a1; v2=gun; v3=a;
1477: for (cz=0; absi_cmp(v3,l) > 0; cz++)
1478: {
1479: p1=dvmdii(d,v3,&t3); t2=subii(v,mulii(p1,v2));
1480: v=v2; d=v3; v2=t2; v3=t3;
1481: }
1482: z=cgetg(4,t_QFI);
1483: if (!cz)
1484: {
1485: q1=mulii(a2,v3); q2=addii(q1,n);
1486: g=divii(addii(mulii(v3,s),(GEN)y[3]),d);
1487: z[1]=lmulii(d,a2);
1488: z[2]=laddii(shifti(q1,1),(GEN)y[2]);
1489: z[3]=laddii(mulii(v3,divii(q2,d)), mulii(g,d1));
1490: }
1491: else
1492: {
1493: if (cz&1) { v3=negi(v3); v2=negi(v2); }
1494: b=divii(addii(mulii(a2,d),mulii(n,v)),a1);
1495: q1=mulii(b,v3); q2=addii(q1,n);
1496: e=divii(addii(mulii(s,d),mulii((GEN)y[3],v)),a1);
1497: q3=mulii(e,v2); q4=subii(q3,s);
1498: g=divii(q4,v); b2=addii(q3,q4);
1499: if (!gcmp1(d1)) { v2=mulii(d1,v2); v=mulii(d1,v); b2=mulii(d1,b2); }
1500: z[1]=laddii(mulii(d,b),mulii(e,v));
1501: z[2]=laddii(b2,addii(q1,q2));
1502: z[3]=laddii(mulii(v3,divii(q2,d)), mulii(g,v2));
1503: }
1504: tetpil=avma; return gerepile(av,tetpil,redimag(z));
1505: }
1506:
1507: GEN
1508: nudupl(GEN x, GEN l)
1509: {
1.2 ! noro 1510: gpmem_t av=avma,tetpil;
! 1511: long cz;
1.1 noro 1512: GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g;
1513:
1514: if (typ(x) != t_QFI)
1515: err(talker,"not an imaginary quadratic form in nudupl");
1516: d1=bezout((GEN)x[2],(GEN)x[1],&u,&v);
1517: a=divii((GEN)x[1],d1); b=divii((GEN)x[2],d1);
1518: c=modii(negi(mulii(u,(GEN)x[3])),a); p1=subii(a,c);
1519: if (cmpii(c,p1)>0) c=negi(p1);
1520: v=gzero; d=a; v2=gun; v3=c;
1521: for (cz=0; absi_cmp(v3,l) > 0; cz++)
1522: {
1523: p1=dvmdii(d,v3,&t3); t2=subii(v,mulii(p1,v2));
1524: v=v2; d=v3; v2=t2; v3=t3;
1525: }
1526: z=cgetg(4,t_QFI);
1527: if (!cz)
1528: {
1529: g=divii(addii(mulii(v3,b),(GEN)x[3]),d);
1530: z[1]=(long)sqri(d);
1531: z[2]=laddii((GEN)x[2],shifti(mulii(d,v3),1));
1532: z[3]=laddii(sqri(v3),mulii(g,d1));
1533: }
1534: else
1535: {
1536: if (cz&1) { v=negi(v); d=negi(d); }
1537: e=divii(addii(mulii((GEN)x[3],v),mulii(b,d)),a);
1538: g=divii(subii(mulii(e,v2),b),v);
1539: b2=addii(mulii(e,v2),mulii(v,g));
1540: if (!gcmp1(d1)) { v2=mulii(d1,v2); v=mulii(d1,v); b2=mulii(d1,b2); }
1541: z[1]=laddii(sqri(d),mulii(e,v));
1542: z[2]=laddii(b2,shifti(mulii(d,v3),1));
1543: z[3]=laddii(sqri(v3),mulii(g,v2));
1544: }
1545: tetpil=avma; return gerepile(av,tetpil,redimag(z));
1546: }
1547:
1548: GEN
1549: nupow(GEN x, GEN n)
1550: {
1.2 ! noro 1551: gpmem_t av,tetpil;
! 1552: long i,j;
1.1 noro 1553: unsigned long m;
1554: GEN y,l;
1555:
1556: if (typ(n) != t_INT) err(talker,"not an integer exponent in nupow");
1557: if (gcmp1(n)) return gcopy(x);
1558: av=avma; y = imag_unit_form(x);
1559: if (!signe(n)) return y;
1560:
1561: l = racine(shifti(racine((GEN)y[3]),1));
1562: for (i=lgefint(n)-1; i>2; i--)
1563: for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
1564: {
1565: if (m&1) y=nucomp(y,x,l);
1566: x=nudupl(x,l);
1567: }
1568: for (m=n[2]; m>1; m>>=1)
1569: {
1570: if (m&1) y=nucomp(y,x,l);
1571: x=nudupl(x,l);
1572: }
1573: tetpil=avma; y=nucomp(y,x,l);
1574: if (signe(n)<0 && !egalii((GEN)y[1],(GEN)y[2])
1575: && !egalii((GEN)y[1],(GEN)y[3]))
1576: setsigne(y[2],-signe(y[2]));
1577: return gerepile(av,tetpil,y);
1578: }
1579:
1580: /* reduction */
1581:
1582: static GEN
1583: abs_dvmdii(GEN b, GEN p1, GEN *pt, long s)
1584: {
1585: if (s<0) setsigne(b, 1); p1 = dvmdii(b,p1,pt);
1586: if (s<0) setsigne(b,-1); return p1;
1587: }
1588:
1589: static GEN
1590: rhoimag0(GEN x, long *flag)
1591: {
1592: GEN p1,b,d,z;
1593: long fl, s = signe(x[2]);
1594:
1595: fl = cmpii((GEN)x[1], (GEN)x[3]);
1596: if (fl <= 0)
1597: {
1598: long fg = absi_cmp((GEN)x[1], (GEN)x[2]);
1599: if (fg >= 0)
1600: {
1601: *flag = (s<0 && (!fl || !fg))? 2 /* set x[2] = negi(x[2]) in caller */
1602: : 1;
1603: return x;
1604: }
1605: }
1606: p1=shifti((GEN)x[3],1); d = abs_dvmdii((GEN)x[2],p1,&b,s);
1607: if (s>=0)
1608: {
1609: setsigne(d,-signe(d));
1610: if (cmpii(b,(GEN)x[3])<=0) setsigne(b,-signe(b));
1611: else { d=addsi(-1,d); b=subii(p1,b); }
1612: }
1613: else if (cmpii(b,(GEN)x[3])>=0) { d=addsi(1,d); b=subii(b,p1); }
1614:
1615: z=cgetg(4,t_QFI);
1616: z[1] = x[3];
1617: z[3] = laddii((GEN)x[1], mulii(d,shifti(subii((GEN)x[2],b),-1)));
1618: if (signe(b)<0 && !fl) setsigne(b,1);
1619: z[2] = (long)b; *flag = 0; return z;
1620: }
1621:
1622: static void
1623: fix_expo(GEN x)
1624: {
1625: long e = expo(x[5]);
1626: if (e >= EXP220)
1627: {
1628: x[4] = laddsi(1,(GEN)x[4]);
1629: setexpo(x[5], e - EXP220);
1630: }
1631: }
1632:
1633: GEN
1634: rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD)
1635: {
1636: GEN p1,p2, y = cgetg(6,t_VEC);
1637: GEN b = (GEN)x[2];
1638: GEN c = (GEN)x[3];
1639:
1640: y[1] = (long)c;
1.2 ! noro 1641: p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: absi(c);
1.1 noro 1642: p1 = shifti(c,1);
1.2 ! noro 1643: if (p1 == gzero) err(talker, "reducible form in rhoreal");
! 1644: setsigne(p1,1); /* |2c| */
1.1 noro 1645: p2 = divii(addii(p2,b), p1);
1646: y[2] = lsubii(mulii(p2,p1), b);
1647:
1648: p1 = shifti(subii(sqri((GEN)y[2]),D),-2);
1649: y[3] = ldivii(p1,(GEN)y[1]);
1650:
1651: if (lg(x) <= 5) setlg(y,4);
1652: else
1653: {
1654: y[4] = x[4];
1655: y[5] = x[5];
1656: if (signe(b))
1657: {
1658: p1 = divrr(addir(b,sqrtD), subir(b,sqrtD));
1659: y[5] = lmulrr(p1, (GEN)y[5]);
1660: fix_expo(y);
1661: }
1662: }
1663: return y;
1664: }
1665:
1666: #define qf_NOD 2
1667: #define qf_STEP 1
1668:
1669: GEN
1670: codeform5(GEN x, long prec)
1671: {
1672: GEN y = cgetg(6,t_VEC);
1.2 ! noro 1673: y[1] = x[1];
! 1674: y[2] = x[2];
! 1675: y[3] = x[3];
! 1676: y[4] = zero;
! 1677: y[5] = (long)realun(prec); return y;
1.1 noro 1678: }
1679:
1680: static GEN
1681: add_distance(GEN x, GEN d0)
1682: {
1683: GEN y = cgetg(5, t_QFR);
1684: y[1] = licopy((GEN)x[1]);
1685: y[2] = licopy((GEN)x[2]);
1686: y[3] = licopy((GEN)x[3]);
1687: y[4] = lcopy(d0); return y;
1688: }
1689:
1690: /* d0 = initial distance, assume |n| > 1 */
1691: static GEN
1692: decodeform(GEN x, GEN d0)
1693: {
1694: GEN p1,p2;
1695:
1696: if (lg(x) < 6) return add_distance(x,d0);
1697: /* x = (a,b,c, expo(d), d), d = exp(2*distance) */
1698: p1 = absr((GEN)x[5]);
1699: p2 = (GEN)x[4];
1700: if (signe(p2))
1701: {
1702: long e = expo(p1);
1703: p1 = shiftr(p1,-e);
1704: p2 = addis(mulsi(EXP220,p2), e);
1705: p1 = mplog(p1);
1.2 ! noro 1706: p1 = mpadd(p1, mulir(p2, mplog2(lg(d0))));
1.1 noro 1707: }
1708: else
1709: { /* to avoid loss of precision */
1710: p1 = gcmp1(p1)? NULL: mplog(p1);
1711: }
1712: if (p1)
1713: d0 = addrr(d0, shiftr(p1,-1));
1714: return add_distance(x, d0);
1715: }
1716:
1717: static long
1718: get_prec(GEN d)
1719: {
1720: long k = lg(d);
1721: long l = ((BITS_IN_LONG-1-expo(d))>>TWOPOTBITS_IN_LONG)+2;
1722: if (l < k) l = k;
1723: if (l < 3) l = 3;
1724: return l;
1725: }
1726:
1727: static int
1728: real_isreduced(GEN x, GEN isqrtD)
1729: {
1730: GEN a = (GEN)x[1];
1731: GEN b = (GEN)x[2];
1732: if (signe(b) > 0 && cmpii(b,isqrtD) <= 0 )
1733: {
1734: GEN p1 = subii(isqrtD, shifti(absi(a),1));
1735: long l = absi_cmp(b, p1);
1736: if (l > 0 || (l == 0 && signe(p1) < 0)) return 1;
1737: }
1738: return 0;
1739: }
1740:
1741: GEN
1742: redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD)
1743: {
1744: while (!real_isreduced(x,isqrtD))
1745: x = rhoreal_aux(x,D,sqrtD,isqrtD);
1746: return x;
1747: }
1748:
1749: static GEN
1750: redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
1751: {
1.2 ! noro 1752: gpmem_t av = avma;
! 1753: long prec;
1.1 noro 1754: GEN d0;
1755:
1756: if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal");
1757:
1758: if (!D)
1759: D = qf_disc(x,NULL,NULL);
1760: else
1761: if (typ(D) != t_INT) err(arither1);
1762:
1763: d0 = (GEN)x[4]; prec = get_prec(d0);
1764: x = codeform5(x,prec);
1765: if ((flag & qf_NOD)) setlg(x,4);
1766: else
1767: {
1768: if (!sqrtD)
1769: sqrtD = gsqrt(D,prec);
1770: else
1771: {
1772: long l = typ(sqrtD);
1773: if (!is_intreal_t(l)) err(arither1);
1774: }
1775: }
1776: if (!isqrtD)
1777: isqrtD = sqrtD? mptrunc(sqrtD): racine(D);
1778: else
1779: if (typ(isqrtD) != t_INT) err(arither1);
1780:
1781: if (flag & qf_STEP)
1782: x = rhoreal_aux(x,D,sqrtD,isqrtD);
1783: else
1784: x = redrealform5(x,D,sqrtD,isqrtD);
1785: return gerepileupto(av, decodeform(x,d0));
1786: }
1787:
1788: GEN
1789: comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD)
1790: {
1.2 ! noro 1791: gpmem_t av = avma;
1.1 noro 1792: GEN z = cgetg(6,t_VEC); comp_gen(z,x,y);
1793: if (x == y)
1794: {
1795: z[4] = lshifti((GEN)x[4],1);
1796: z[5] = lsqr((GEN)x[5]);
1797: }
1798: else
1799: {
1800: z[4] = laddii((GEN)x[4],(GEN)y[4]);
1801: z[5] = lmulrr((GEN)x[5],(GEN)y[5]);
1802: }
1803: fix_expo(z); z = redrealform5(z,D,sqrtD,isqrtD);
1804: return gerepilecopy(av,z);
1805: }
1806:
1807: /* assume n!=0 */
1808: GEN
1809: powrealform(GEN x, GEN n)
1810: {
1.2 ! noro 1811: gpmem_t av = avma;
! 1812: long i,m;
1.1 noro 1813: GEN y,D,sqrtD,isqrtD,d0;
1814:
1815: if (typ(x) != t_QFR)
1816: err(talker,"not a real quadratic form in powreal");
1817: if (gcmp1(n)) return gcopy(x);
1818: if (gcmp_1(n)) return ginv(x);
1819:
1820: d0 = (GEN)x[4];
1821: D = qf_disc(x,NULL,NULL);
1822: sqrtD = gsqrt(D, get_prec(d0));
1823: isqrtD = mptrunc(sqrtD);
1824: if (signe(n) < 0) { x = ginv(x); d0 = (GEN)x[4]; }
1825: n = absi(n);
1826: x = codeform5(x, lg(d0)); y = NULL;
1827: for (i=lgefint(n)-1; i>1; i--)
1828: {
1829: m = n[i];
1830: for (; m; m>>=1)
1831: {
1832: if (m&1) y = y? comprealform5(y,x,D,sqrtD,isqrtD): x;
1833: if (m == 1 && i == 2) break;
1834: x = comprealform5(x,x,D,sqrtD,isqrtD);
1835: }
1836: }
1837: d0 = mulri(d0,n);
1838: return gerepileupto(av, decodeform(y,d0));
1839: }
1840:
1841: GEN
1842: redimag(GEN x)
1843: {
1.2 ! noro 1844: gpmem_t av=avma;
! 1845: long fl;
1.1 noro 1846: do x = rhoimag0(x, &fl); while (fl == 0);
1847: x = gerepilecopy(av,x);
1848: if (fl == 2) setsigne(x[2], -signe(x[2]));
1849: return x;
1850: }
1851:
1852: GEN
1853: redreal(GEN x)
1854: {
1855: return redreal0(x,0,NULL,NULL,NULL);
1856: }
1857:
1858: GEN
1859: rhoreal(GEN x)
1860: {
1861: return redreal0(x,qf_STEP,NULL,NULL,NULL);
1862: }
1863:
1864: GEN
1865: redrealnod(GEN x, GEN isqrtD)
1866: {
1867: return redreal0(x,qf_NOD,NULL,isqrtD,NULL);
1868: }
1869:
1870: GEN
1871: rhorealnod(GEN x, GEN isqrtD)
1872: {
1873: return redreal0(x,qf_STEP|qf_NOD,NULL,isqrtD,NULL);
1874: }
1875:
1876: GEN
1877: qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
1878: {
1879: long tx=typ(x),fl;
1.2 ! noro 1880: gpmem_t av;
1.1 noro 1881:
1882: if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred");
1883: if (!D) D = qf_disc(x,NULL,NULL);
1884: switch(signe(D))
1885: {
1886: case 1 :
1887: return redreal0(x,flag,D,isqrtD,sqrtD);
1888:
1889: case -1:
1890: if (!flag) return redimag(x);
1891: if (flag !=1) err(flagerr,"qfbred");
1892: av=avma; x = rhoimag0(x,&fl);
1893: x = gerepilecopy(av,x);
1894: if (fl == 2) setsigne(x[2], -signe(x[2]));
1895: return x;
1896: }
1897: err(redpoler,"qfbred");
1898: return NULL; /* not reached */
1899: }
1900:
1901: /* special case: p = 1 return unit form */
1902: GEN
1903: primeform(GEN x, GEN p, long prec)
1904: {
1.2 ! noro 1905: gpmem_t av;
! 1906: long s, sx = signe(x);
! 1907: GEN y, b, c;
1.1 noro 1908:
1.2 ! noro 1909: if (typ(x) != t_INT || !sx) err(arither1);
1.1 noro 1910: if (typ(p) != t_INT || signe(p) <= 0) err(arither1);
1911: if (is_pm1(p))
1.2 ! noro 1912: return sx<0? imag_unit_form_by_disc(x)
! 1913: : real_unit_form_by_disc(x,prec);
! 1914: s = mod8(x);
! 1915: if (sx < 0)
1.1 noro 1916: {
1.2 ! noro 1917: if (s) s = 8-s;
! 1918: y = cgetg(4, t_QFI);
1.1 noro 1919: }
1920: else
1921: {
1.2 ! noro 1922: y = cgetg(5, t_QFR);
! 1923: y[4] = (long)realzero(prec);
1.1 noro 1924: }
1925: switch(s&3)
1926: {
1927: case 2: case 3: err(funder2,"primeform");
1928: }
1.2 ! noro 1929: av = avma;
! 1930: if (egalii(p, gdeux))
1.1 noro 1931: {
1932: switch(s)
1933: {
1.2 ! noro 1934: case 0: b = gzero; break;
! 1935: case 1: b = gun; break;
! 1936: case 4: b = gdeux; break;
! 1937: default: err(sqrter5); b = NULL; /* -Wall */
1.1 noro 1938: }
1.2 ! noro 1939: c = shifti(subsi(s,x), -3);
1.1 noro 1940: }
1941: else
1942: {
1943: b = mpsqrtmod(x,p); if (!b) err(sqrter5);
1.2 ! noro 1944: s &= 1; /* s = x mod 2 */
! 1945: /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */
! 1946: if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(p,b));
! 1947:
! 1948: av = avma;
! 1949: c = diviiexact(shifti(subii(sqri(b), x), -2), p);
! 1950: }
! 1951: y[3] = lpileuptoint(av, c);
! 1952: y[2] = (long)b;
! 1953: y[1] = licopy(p); return y;
1.1 noro 1954: }
1955:
1956: /*********************************************************************/
1957: /** **/
1958: /** BINARY DECOMPOSITION **/
1959: /** **/
1960: /*********************************************************************/
1961:
1962: GEN
1963: binaire(GEN x)
1964: {
1965: ulong m,u;
1966: long i,lx,ex,ly,tx=typ(x);
1967: GEN y,p1,p2;
1968:
1969: switch(tx)
1970: {
1971: case t_INT:
1972: lx=lgefint(x);
1973: if (lx==2) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
1974: ly = BITS_IN_LONG+1; m=HIGHBIT; u=x[2];
1975: while (!(m & u)) { m>>=1; ly--; }
1976: y = cgetg(ly+((lx-3)<<TWOPOTBITS_IN_LONG),t_VEC); ly=1;
1977: do { y[ly] = m & u ? un : zero; ly++; } while (m>>=1);
1978: for (i=3; i<lx; i++)
1979: {
1980: m=HIGHBIT; u=x[i];
1981: do { y[ly] = m & u ? un : zero; ly++; } while (m>>=1);
1982: }
1983: break;
1984:
1985: case t_REAL:
1986: ex=expo(x);
1987: if (!signe(x))
1988: {
1989: lx=1+max(-ex,0); y=cgetg(lx,t_VEC);
1990: for (i=1; i<lx; i++) y[i]=zero;
1991: return y;
1992: }
1993:
1994: lx=lg(x); y=cgetg(3,t_VEC);
1995: if (ex > bit_accuracy(lx)) err(precer,"binary");
1996: p1 = cgetg(max(ex,0)+2,t_VEC);
1997: p2 = cgetg(bit_accuracy(lx)-ex,t_VEC);
1998: y[1] = (long) p1; y[2] = (long) p2;
1999: ly = -ex; ex++; m = HIGHBIT;
2000: if (ex<=0)
2001: {
2002: p1[1]=zero; for (i=1; i <= -ex; i++) p2[i]=zero;
2003: i=2;
2004: }
2005: else
2006: {
2007: ly=1;
2008: for (i=2; i<lx && ly<=ex; i++)
2009: {
2010: m=HIGHBIT; u=x[i];
2011: do
2012: { p1[ly] = (m & u) ? un : zero; ly++; }
2013: while ((m>>=1) && ly<=ex);
2014: }
2015: ly=1;
2016: if (m) i--; else m=HIGHBIT;
2017: }
2018: for (; i<lx; i++)
2019: {
2020: u=x[i];
2021: do { p2[ly] = m & u ? un : zero; ly++; } while (m>>=1);
2022: m=HIGHBIT;
2023: }
2024: break;
2025:
2026: case t_VEC: case t_COL: case t_MAT:
2027: lx=lg(x); y=cgetg(lx,tx);
2028: for (i=1; i<lx; i++) y[i]=(long)binaire((GEN)x[i]);
2029: break;
2030: default: err(typeer,"binaire");
2031: return NULL; /* not reached */
2032: }
2033: return y;
2034: }
2035:
2036: /* return 1 if bit n of x is set, 0 otherwise */
2037: long
2038: bittest(GEN x, long n)
2039: {
2040: long l;
2041:
2042: if (!signe(x) || n<0) return 0;
2043: l = lgefint(x)-1 - (n>>TWOPOTBITS_IN_LONG);
2044: if (l < 2) return 0;
2045: n = (1L << (n & (BITS_IN_LONG-1))) & x[l];
2046: return n? 1: 0;
2047: }
2048:
1.2 ! noro 2049: GEN
! 2050: bittest_many(GEN x, GEN gn, long c)
! 2051: {
! 2052: long lx = lgefint(x), l1, l2, b1, b2, two_adic = 0, l_res, skip;
! 2053: ulong *p, *pnew;
! 2054: long n = itos(gn), extra_words = 0, partial_bits;
! 2055: GEN res;
! 2056:
! 2057: if (c == 1) /* Shortcut */
! 2058: return bittest(x,n) ? gun : gzero;
! 2059: /* Negative n with n+c>0 are too hairy to implement... */
! 2060: if (!signe(x) || c == 0)
! 2061: return gzero;
! 2062: if (c < 0) { /* Negative x means 2s complemenent */
! 2063: c = -c;
! 2064: if (signe(x) < 0)
! 2065: two_adic = 1; /* treat x 2-adically... */
! 2066: }
! 2067: if (n < 0) {
! 2068: gpmem_t oa, na;
! 2069:
! 2070: if (n + c <= 0)
! 2071: return gzero;
! 2072: oa = avma;
! 2073: res = bittest_many(x, gzero, two_adic? -(n+c) : n+c);
! 2074: na = avma;
! 2075: res = shifti3(res,-n,0);
! 2076: return gerepile(oa,na,res);
! 2077: }
! 2078: partial_bits = (c & (BITS_IN_LONG-1));
! 2079: l2 = lx-1 - (n>>TWOPOTBITS_IN_LONG); /* Last=least-significant word */
! 2080: l1 = lx-1 - ((n + c - 1)>>TWOPOTBITS_IN_LONG); /* First word */
! 2081: b2 = (n & (BITS_IN_LONG-1)); /* Last bit, skip less-signifant */
! 2082: b1 = ((n + c - 1) & (BITS_IN_LONG-1)); /* First bit, skip more-significant */
! 2083: if (l2 < 2) { /* always: l1 <= l2 */
! 2084: if (!two_adic)
! 2085: return gzero;
! 2086: /* x is non-zero, so it behaves as -1. */
! 2087: return gbitneg(gzero,c);
! 2088: }
! 2089: if (l1 < 2) {
! 2090: if (two_adic) /* If b1 < b2, bits are set by prepend in shift_r */
! 2091: extra_words = 2 - l1 - (b1 < b2);
! 2092: else
! 2093: partial_bits = b2 ? BITS_IN_LONG - b2 : 0;
! 2094: l1 = 2;
! 2095: b1 = (BITS_IN_LONG-1); /* Include all bits in this word */
! 2096: }
! 2097: skip = (b1 < b2); /* Skip the first word of the shift */
! 2098: l_res = l2 - l1 + 1 + 2 - skip + extra_words; /* A coarse estimate */
! 2099: p = (ulong*) (new_chunk(l_res) + 2 + extra_words);
! 2100: shift_r(p - skip, (ulong*)x + l1, (ulong*)x + l2 + 1, 0, b2);
! 2101: if (two_adic) { /* Check the low bits of x */
! 2102: int i = lx, nonzero = 0;
! 2103:
! 2104: /* Flip the bits */
! 2105: pnew = p + l_res - 2 - extra_words;
! 2106: while (--pnew >= p)
! 2107: *pnew = MAXULONG - *pnew;
! 2108: /* Fill the extra words */
! 2109: while (extra_words--)
! 2110: *--p = MAXULONG;
! 2111: /* The result is one less than 2s-complement if the lower-bits
! 2112: of x were all 0. */
! 2113: while (--i > l2) {
! 2114: if (x[i]) {
! 2115: nonzero = 1;
! 2116: break;
! 2117: }
! 2118: }
! 2119: if (!nonzero && b2) /* Check the partial word too */
! 2120: nonzero = x[l2] & ((1<<b2) - 1);
! 2121: if (!nonzero) { /* Increment res. Do not underflow to before p */
! 2122: pnew = p + l_res - 2;
! 2123: while (--pnew >= p)
! 2124: if (++*pnew)
! 2125: break;
! 2126: }
! 2127: }
! 2128: if (partial_bits)
! 2129: *p &= (1<<partial_bits) - 1;
! 2130: pnew = p;
! 2131: while ((pnew < p + l_res - 2) && !*pnew) /* Skip 0 words */
! 2132: pnew++;
! 2133: avma = (gpmem_t)(pnew - 2);
! 2134: if (pnew >= p - 2 + l_res)
! 2135: return gzero;
! 2136: l_res -= (pnew - p);
! 2137: res = (GEN)(pnew - 2);
! 2138: res[1] = evalsigne(1)|evallgefint(l_res);
! 2139: res[0] = evaltyp(t_INT)|evallg(l_res);
! 2140: return res;
! 2141: }
! 2142:
1.1 noro 2143: static long
2144: bittestg(GEN x, GEN n)
2145: {
2146: return bittest(x,itos(n));
2147: }
2148:
2149: GEN
2150: gbittest(GEN x, GEN n)
2151: {
2152: return arith_proto2(bittestg,x,n);
2153: }
2154:
1.2 ! noro 2155: GEN
! 2156: gbittest3(GEN x, GEN n, long c)
! 2157: {
! 2158: return garith_proto3ggs(bittest_many,x,n,c);
! 2159: }
! 2160:
1.1 noro 2161: /***********************************************************************/
2162: /** **/
2163: /** BITMAP OPS **/
2164: /** x & y (and), x | y (or), x ^ y (xor), ~x (neg), x & ~y (negimply) **/
2165: /** **/
2166: /***********************************************************************/
2167:
2168: /* Normalize a non-negative integer. */
2169: static void
2170: inormalize(GEN x, long known_zero_words)
2171: {
2172: int xl = lgefint(x);
2173: int i, j;
2174:
2175: /* Normalize */
2176: i = 2 + known_zero_words;
2177: while (i < xl) {
2178: if (x[i])
2179: break;
2180: i++;
2181: }
2182: j = 2;
2183: while (i < xl)
2184: x[j++] = x[i++];
2185: xl -= i - j;
2186: setlgefint(x, xl);
2187: if (xl == 2)
2188: setsigne(x,0);
2189: }
2190:
2191: /* Truncate a non-negative integer to a number of bits. */
2192: static void
2193: ibittrunc(GEN x, long bits, long normalized)
2194: {
2195: int xl = lgefint(x);
2196: int len_out = ((bits + BITS_IN_LONG - 1) >> TWOPOTBITS_IN_LONG) + 2;
2197: int known_zero_words, i = 2 + xl - len_out;
2198:
2199: if (xl < len_out && normalized)
2200: return;
2201: /* Check whether mask is trivial */
2202: if (!(bits & (BITS_IN_LONG - 1))) {
2203: if (xl == len_out && normalized)
2204: return;
2205: } else if (len_out <= xl) {
2206: /* Non-trival mask is given by a formula, if x is not
2207: normalized, this works even in the exceptional case */
2208: x[i] = x[i] & ((1 << (bits & (BITS_IN_LONG - 1))) - 1);
2209: if (x[i] && xl == len_out)
2210: return;
2211: }
2212: /* Normalize */
2213: if (xl <= len_out) /* Not normalized */
2214: known_zero_words = 0;
2215: else
2216: known_zero_words = xl - len_out;
2217: inormalize(x, known_zero_words);
2218: }
2219:
2220: /* Increment/decrement absolute value of non-zero integer in place.
2221: It is assumed that the increment will not increase the length.
2222: Decrement may result in a non-normalized value.
2223: Returns 1 if the increment overflows (thus the result is 0). */
2224: static int
2225: incdec(GEN x, long incdec)
2226: {
2227: long *xf = x + 2, *xl;
2228: long len = lgefint(x);
2229: const ulong uzero = 0;
2230:
2231: xl = x + len;
2232: if (incdec == 1) {
2233: while (--xl >= xf) {
2234: if ((ulong)*xl != ~uzero) {
2235: (*xl)++;
2236: return 0;
2237: }
2238: *xl = 0;
2239: }
2240: return 1;
2241: } else {
2242: while (--xl >= xf) {
2243: if (*xl != 0) {
2244: (*xl)--;
2245: return 0;
2246: }
2247: *xl = (long)~uzero;
2248: }
2249: return 0;
2250: }
2251: }
2252:
2253: GEN
2254: gbitneg(GEN x, long bits)
2255: {
2256: long xl, len_out, i;
2257: const ulong uzero = 0;
2258:
2259: if (typ(x) != t_INT)
2260: err(typeer, "bitwise negation");
2261: if (bits < -1)
2262: err(talker, "negative exponent in bitwise negation");
2263: if (bits == -1)
2264: return gsub(gneg(gun),x);
2265: if (bits == 0)
2266: return gzero;
2267: if (signe(x) == -1) { /* Consider as if mod big power of 2 */
2268: x = gcopy(x);
2269: setsigne(x, 1);
2270: incdec(x, -1);
2271: /* Now truncate this! */
2272: ibittrunc(x, bits, x[2]);
2273: return x;
2274: }
2275: xl = lgefint(x);
2276: len_out = ((bits + BITS_IN_LONG - 1) >> TWOPOTBITS_IN_LONG) + 2;
2277: if (len_out > xl) { /* Need to grow */
2278: GEN out = cgeti(len_out);
2279: int j = 2;
2280:
2281: if (!(bits & (BITS_IN_LONG - 1)))
2282: out[2] = ~uzero;
2283: else
2284: out[2] = (1 << (bits & (BITS_IN_LONG - 1))) - 1;
2285: for (i = 3; i < len_out - xl + 2; i++)
2286: out[i] = ~uzero;
2287: while (i < len_out)
2288: out[i++] = ~x[j++];
2289: setlgefint(out, len_out);
2290: setsigne(out,1);
2291: return out;
2292: }
2293: x = gcopy(x);
2294: for (i = 2; i < xl; i++)
2295: x[i] = ~x[i];
2296: ibittrunc(x, bits, x[2]);
2297: return x;
2298: }
2299:
2300: /* bitwise 'and' of two positive integers (any integers, but we ignore sign).
2301: * Inputs are not necessary normalized. */
2302: static GEN
2303: ibitand(GEN x, GEN y)
2304: {
2305: long lx, ly, lout;
2306: long *xp, *yp, *outp, *xlim;
2307: GEN out;
2308:
2309: lx=lgefint(x); ly=lgefint(y);
2310: if (lx > ly)
2311: lout = ly;
2312: else
2313: lout = lx;
2314: xlim = x + lx;
2315: xp = xlim + 2 - lout;
2316: yp = y + 2 + ly - lout;
2317: out = cgeti(lout);
2318: outp = out + 2;
2319: while (xp < xlim)
2320: *outp++ = (*xp++) & (*yp++);
2321: setsigne(out,1);
2322: setlgefint(out,lout);
2323: if (lout == 2)
2324: setsigne(out,0);
2325: else if ( !out[2] )
2326: inormalize(out, 1);
2327: return out;
2328: }
2329:
2330: #define swaplen(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
2331:
2332: /* bitwise 'or' of two positive integers (any integers, but we ignore sign).
2333: * Inputs are not necessary normalized. */
2334: static GEN
2335: ibitor(GEN x, GEN y, long exclusive)
2336: {
2337: long lx, ly, lout;
2338: long *xp, *yp, *outp, *xlim, *xprep;
2339: GEN out;
2340:
2341: lx=lgefint(x); ly=lgefint(y);
2342: if (lx < ly)
2343: swaplen(x,y,lx,ly);
2344: lout = lx;
2345: xlim = x + lx;
2346: xp = xlim + 2 - ly;
2347: yp = y + 2;
2348: out = cgeti(lout);
2349: outp = out + 2;
2350: if (lx > ly) {
2351: xprep = x + 2;
2352: while (xprep < xp)
2353: *outp++ = *xprep++;
2354: }
2355: if (exclusive) {
2356: while (xp < xlim)
2357: *outp++ = (*xp++) ^ (*yp++);
2358: } else {
2359: while (xp < xlim)
2360: *outp++ = (*xp++) | (*yp++);
2361: }
2362: setsigne(out,1);
2363: setlgefint(out,lout);
2364: if (lout == 2)
2365: setsigne(out,0);
2366: else if ( !out[2] )
2367: inormalize(out, 1);
2368: return out;
2369: }
2370:
2371: /* bitwise negated 'implies' of two positive integers (any integers, but we
2372: * ignore sign). "Neg-Implies" is x & ~y unless "negated".
2373: * Inputs are not necessary normalized. */
2374: static GEN
2375: ibitnegimply(GEN x, GEN y)
2376: {
2377: long lx, ly, lout, inverted = 0;
2378: long *xp, *yp, *outp, *xlim, *xprep;
2379: GEN out;
2380:
2381: lx=lgefint(x); ly=lgefint(y);
2382: if (lx < ly) {
2383: inverted = 1;
2384: swaplen(x,y,lx,ly);
2385: }
2386: /* x is longer than y */
2387: lout = lx;
2388: xlim = x + lx;
2389: xp = xlim + 2 - ly;
2390: yp = y + 2;
2391: out = cgeti(lout);
2392: outp = out + 2;
2393: if (lx > ly) {
2394: xprep = x + 2;
2395: if (!inverted) { /* x & ~y */
2396: while (xprep < xp)
2397: *outp++ = *xprep++;
2398: } else { /* ~x & y */
2399: while (xprep++ < xp)
2400: *outp++ = 0;
2401: }
2402: }
2403: if (inverted) { /* ~x & y */
2404: while (xp < xlim)
2405: *outp++ = ~(*xp++) & (*yp++);
2406: } else {
2407: while (xp < xlim)
2408: *outp++ = (*xp++) & ~(*yp++);
2409: }
2410: setsigne(out,1);
2411: setlgefint(out,lout);
2412: if (lout == 2)
2413: setsigne(out,0);
2414: else if ( !out[2] )
2415: inormalize(out, 1);
2416: return out;
2417: }
2418:
2419: static GEN
1.2 ! noro 2420: inegate_inplace(GEN z, gpmem_t ltop)
1.1 noro 2421: {
1.2 ! noro 2422: long o;
1.1 noro 2423:
2424: /* Negate z */
2425: o = incdec(z, 1); /* Can overflow z... */
2426: setsigne(z, -1);
2427: if (!o)
2428: return z;
2429: else if (lgefint(z) == 2)
2430: setsigne(z, 0);
2431: incdec(z,-1); /* Restore a normalized value */
1.2 ! noro 2432: return gerepileupto(ltop, subis(z,1));
1.1 noro 2433: }
2434:
2435: GEN
2436: gbitor(GEN x, GEN y)
2437: {
1.2 ! noro 2438: gpmem_t ltop;
1.1 noro 2439: long sx,sy;
2440: GEN z;
2441:
2442: if (typ(x) != t_INT || typ(y) != t_INT)
2443: err(typeer, "bitwise or");
2444: sx=signe(x); if (!sx) return icopy(y);
2445: sy=signe(y); if (!sy) return icopy(x);
2446: if (sx == 1) {
2447: if (sy == 1)
2448: return ibitor(x,y,0);
2449: goto posneg;
2450: } else if (sy == -1) {
2451: ltop = avma;
2452: incdec(x, -1); incdec(y, -1);
2453: z = ibitand(x,y);
2454: incdec(x, 1); incdec(y, 1); /* Restore x and y... */
2455: } else {
2456: z = x; x = y; y = z;
2457: /* x is positive, y is negative. The result will be negative. */
2458: posneg:
2459: ltop = avma;
2460: incdec(y, -1);
2461: z = ibitnegimply(y,x); /* ~x & y */
2462: incdec(y, 1);
2463: }
2464: return inegate_inplace(z, ltop);
2465: }
2466:
2467: GEN
2468: gbitand(GEN x, GEN y)
2469: {
1.2 ! noro 2470: gpmem_t ltop;
1.1 noro 2471: long sx,sy;
2472: GEN z;
2473:
2474: if (typ(x) != t_INT || typ(y) != t_INT)
2475: err(typeer, "bitwise and");
2476: sx=signe(x); if (!sx) return gzero;
2477: sy=signe(y); if (!sy) return gzero;
2478: if (sx == 1) {
2479: if (sy == 1)
2480: return ibitand(x,y);
2481: goto posneg;
2482: } else if (sy == -1) {
2483: ltop = avma;
2484: incdec(x, -1); incdec(y, -1);
2485: z = ibitor(x,y,0);
2486: incdec(x, 1); incdec(y, 1); /* Restore x and y... */
2487: return inegate_inplace(z, ltop);
2488: } else {
2489: z = x; x = y; y = z;
2490: /* x is positive, y is negative. The result will be positive. */
2491: posneg:
2492: ltop = avma;
2493: incdec(y, -1);
2494: /* x & ~y */
2495: z = ibitnegimply(x,y); /* x & ~y */
2496: incdec(y, 1);
2497: return z;
2498: }
2499: }
2500:
2501: GEN
2502: gbitxor(GEN x, GEN y)
2503: {
1.2 ! noro 2504: gpmem_t ltop;
1.1 noro 2505: long sx,sy;
2506: GEN z;
2507:
2508: if (typ(x) != t_INT || typ(y) != t_INT)
2509: err(typeer, "bitwise xor");
2510: sx=signe(x); if (!sx) return icopy(y);
2511: sy=signe(y); if (!sy) return icopy(x);
2512: if (sx == 1) {
2513: if (sy == 1)
2514: return ibitor(x,y,1);
2515: goto posneg;
2516: } else if (sy == -1) {
2517: incdec(x, -1); incdec(y, -1);
2518: z = ibitor(x,y,1);
2519: incdec(x, 1); incdec(y, 1); /* Restore x and y... */
2520: return z;
2521: } else {
2522: z = x; x = y; y = z;
2523: /* x is positive, y is negative. The result will be negative. */
2524: posneg:
2525: ltop = avma;
2526: incdec(y, -1);
2527: /* ~(x ^ ~y) == x ^ y */
2528: z = ibitor(x,y,1);
2529: incdec(y, 1);
2530: }
2531: return inegate_inplace(z, ltop);
2532: }
2533:
2534: GEN
2535: gbitnegimply(GEN x, GEN y) /* x & ~y */
2536: {
1.2 ! noro 2537: gpmem_t ltop;
1.1 noro 2538: long sx,sy;
2539: GEN z;
2540:
2541: if (typ(x) != t_INT || typ(y) != t_INT)
2542: err(typeer, "bitwise negated imply");
2543: sx=signe(x); if (!sx) return gzero;
2544: sy=signe(y); if (!sy) return icopy(x);
2545: if (sx == 1) {
2546: if (sy == 1)
2547: return ibitnegimply(x,y);
2548: /* x is positive, y is negative. The result will be positive. */
2549: incdec(y, -1);
2550: z = ibitand(x,y);
2551: incdec(y, 1);
2552: return z;
2553: } else if (sy == -1) {
2554: /* both negative. The result will be positive. */
2555: incdec(x, -1); incdec(y, -1);
2556: /* ((~x) & ~(~y)) == ~x & y */
2557: z = ibitnegimply(y,x);
2558: incdec(x, 1); incdec(y, 1); /* Restore x and y... */
2559: return z;
2560: } else {
2561: /* x is negative, y is positive. The result will be negative. */
2562: ltop = avma;
2563: incdec(x, -1);
2564: /* ~((~x) & ~y) == x | y */
2565: z = ibitor(x,y,0);
2566: incdec(x, 1);
2567: }
2568: return inegate_inplace(z, ltop);
2569: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to arith2.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath