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