Annotation of OpenXM_contrib/pari/src/basemath/arith2.c, Revision 1.1.1.1
1.1 maekawa 1: /*********************************************************************/
2: /** **/
3: /** ARITHMETIC FUNCTIONS **/
4: /** (second part) **/
5: /** **/
6: /*********************************************************************/
7: /* $Id: arith2.c,v 1.1.1.1 1999/09/16 13:47:17 karim Exp $ */
8: #include "pari.h"
9:
10: GEN arith_proto(long f(GEN), GEN x, int do_error);
11: GEN garith_proto(GEN f(GEN), GEN x, int do_error);
12: GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y);
13: static long court_p[]={evaltyp(t_INT)|m_evallg(3),evalsigne(1)|evallgefint(3),0};
14:
15: /***********************************************************************/
16: /** **/
17: /** PRIME NUMBERS **/
18: /** **/
19: /***********************************************************************/
20:
21: GEN
22: prime(long n)
23: {
24: byteptr p = diffptr;
25: long c, prime = 0;
26:
27: while (n--)
28: {
29: c = *p++; if (!c) err(primer1);
30: prime += c;
31: }
32: return stoi(prime);
33: }
34:
35: GEN
36: primes(long n)
37: {
38: byteptr p = diffptr;
39: long c, prime = 0;
40: GEN y,z;
41:
42: z = y = cgetg(n+1,t_VEC);
43: while (n--)
44: {
45: c = *p++; if (!c) err(primer1);
46: prime += c; *++z = lstoi(prime);
47: }
48: return y;
49: }
50:
51: /* Building/Rebuilding the diffptr table. Incorporates Ilya Zakharevich's
52: * block sweep algorithm (see pari-dev-111, 25 April 1998). The actual work
53: * is done by the following two subroutines; the user entry point is the
54: * function initprimes() below. initprimes1() is the old algorithm, called
55: * when maxnum (size) is moderate.
56: */
57: static byteptr
58: initprimes1(long size, long *lenp, long *lastp)
59: {
60: long k;
61: byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2);
62:
63: memset(p, 0, size + 2); fin = p + size;
64: for (r=q=p,k=1; r<=fin; )
65: {
66: do { r+=k; k+=2; r+=k; } while (*++q);
67: for (s=r; s<=fin; s+=k) *s=1;
68: }
69: r=p; *r++=2; *r++=1; /* 2 and 3 */
70: for (s=q=r-1; ; s=q)
71: {
72: do q++; while (*q);
73: if (q>fin) break;
74: *r++ = (q-s) << 1;
75: }
76: *r++=0;
77: *lenp = r - p;
78: *lastp = ((s - p) << 1) + 1;
79: #if 0
80: fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n",
81: q, *lenp, *lastp);
82: #endif
83: return (byteptr) gprealloc(p,r-p,size + 2);
84: }
85:
86: #define PRIME_ARENA (512 * 1024)
87:
88: /* Here's the workhorse. This is recursive, although normally the first
89: recursive call will bottom out and invoke initprimes1() at once.
90: (Not static; might conceivably be useful to someone in library mode) */
91: byteptr
92: initprimes0(ulong maxnum, long *lenp, long *lastp)
93: {
94: long k, size, alloced, asize, psize, rootnum, curlow, last, remains, need;
95: byteptr q,r,s,fin, p, p1, fin1, plast, curdiff;
96:
97: #if 0
98: fprintferr("initprimes0: called for maxnum = %lu\n", maxnum);
99: #endif
100: if (maxnum <= 1ul<<17) /* Arbitrary. */
101: return initprimes1(maxnum>>1, lenp, lastp);
102:
103: maxnum |= 1; /* make it odd. */
104:
105: /* Checked to be enough up to 40e6, attained at 155893 */
106: size = (long) (1.09 * maxnum/log((double)maxnum)) + 145;
107: p1 = (byteptr) gpmalloc(size);
108: rootnum = (long) sqrt((double)maxnum); /* cast it back to a long */
109: rootnum |= 1;
110: #if 0
111: fprintferr("initprimes0: rootnum = %ld\n", rootnum);
112: #endif
113: {
114: byteptr p2 = initprimes0(rootnum, &psize, &last); /* recursive call */
115: memcpy(p1, p2, psize); free(p2);
116: }
117: fin1 = p1 + psize - 1;
118: remains = (maxnum - rootnum) >> 1; /* number of odd numbers to check */
119:
120: need = 100 * rootnum; /* Make % overhead negligeable. */
121: if (need < PRIME_ARENA) need = PRIME_ARENA;
122: if (avma - bot < need>>1) { /* need to do our own allocation */
123: alloced = 1; asize = need;
124: } else { /* scratch area is free part of PARI stack */
125: alloced = 0; asize = avma - bot;
126: }
127: if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */
128: if (alloced)
129: p = (byteptr) gpmalloc(asize);
130: else
131: p = (byteptr) bot;
132: fin = p + asize - 1; /* the 0 sentinel goes at fin. */
133: curlow = rootnum + 2; /* We know all primes up to rootnum (odd). */
134: curdiff = fin1;
135:
136: /* During each iteration p..fin-1 represents a range of odd
137: numbers. plast is a pointer which represents the last prime seen,
138: it may point before p..fin-1. */
139: plast = p - ((rootnum - last) >> 1) - 1;
140: while (remains)
141: {
142: if (asize > remains)
143: {
144: asize = remains + 1;
145: fin = p + asize - 1;
146: }
147: memset(p, 0, asize);
148: /* p corresponds to curlow. q runs over primediffs. */
149: for (q = p1+2, k = 3; q <= fin1; k += *q++)
150: {
151: /* The first odd number which is >= curlow and divisible by p
152: equals (if curlow > p)
153: the last odd number which is <= curlow + 2p - 1 and 0 (mod p)
154: which equals
155: p + the last even number which is <= curlow + p - 1 and 0 (mod p)
156: which equals
157: p + the last even number which is <= curlow + p - 2 and 0 (mod p)
158: which equals
159: p + curlow + p - 2 - (curlow + p - 2)) % 2p.
160: */
161: long k2 = k*k - curlow;
162:
163: if (k2 > 0)
164: {
165: r = p + (k2 >>= 1);
166: if (k2 > remains) break; /* Guard against an address wrap. */
167: } else
168: r = p - (((curlow+k-2) % (2*k)) >> 1) + k - 1;
169: for (s = r; s < fin; s += k) *s = 1;
170: }
171: /* now q runs over addresses corresponding to primes */
172: for (q = p; ; plast = q++)
173: {
174: while (*q) q++; /* use 0-sentinel at end */
175: if (q >= fin) break;
176: *curdiff++ = (q - plast) << 1;
177: }
178: plast -= asize - 1;
179: remains -= asize - 1;
180: curlow += ((asize - 1)<<1);
181: } /* while (remains) */
182: last = curlow - ((p - plast) << 1);
183: *curdiff++ = 0; /* sentinel */
184: *lenp = curdiff - p1;
185: *lastp = last;
186: if (alloced) free(p);
187: return (byteptr) gprealloc(p1, *lenp, size);
188: }
189: #if 0 /* not yet... GN */
190: /* The diffptr table will contain at least 6548 entries (up to and including
191: the 6547th prime, 65557, and the terminal null byte), because the isprime/
192: small-factor-extraction machinery wants to depend on everything up to 65539
193: being in the table, and we might as well go to a multiple of 4 Bytes.--GN */
194:
195: void
196: init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */
197: #endif
198: byteptr
199: initprimes(long maxnum)
200: {
201: long len, last;
202: byteptr p;
203: /* The algorithm must see the next prime beyond maxnum, whence the +257. */
204: /* switch to unsigned: adding the 257 _could_ wrap to a negative number. */
205: ulong maxnum1 = ((maxnum<65302?65302:maxnum)+257ul);
206:
207: if (maxnum1 > 436273000)
208: err(talker,"impossible to have prestored primes > 436273009");
209:
210: p = initprimes0(maxnum1, &len, &last);
211: #if 0 /* not yet... GN */
212: static int build_the_tables = 1;
213: if (build_the_tables) { init_tinyprimes_tridiv(p); build_the_tables=0; }
214: #endif
215: return p;
216: }
217:
218: static void
219: cleanprimetab(void)
220: {
221: long i,j, lp = lg(primetab);
222:
223: for (i=j=1; i < lp; i++)
224: if (primetab[i]) primetab[j++] = primetab[i];
225: setlg(primetab,j);
226: }
227:
228: /* primes is a single element or a row vector with "primes" to add to
229: * primetab. If the "prime" shares a non-trivial factor with another element
230: * in primetab, it is taken into account.
231: */
232: GEN
233: addprimes(GEN prime)
234: {
235: long i,tx, av = avma, lp = lg(primetab);
236: GEN p1,p2;
237:
238: if (!prime) return primetab;
239: tx = typ(prime);
240: if (is_vec_t(tx))
241: {
242: for (i=1; i < lg(prime); i++) addprimes((GEN) prime[i]);
243: return primetab;
244: }
245: if (tx != t_INT) err(typeer,"addprime");
246: if (is_pm1(prime)) return primetab;
247: if (signe(prime) < 0) prime=absi(prime);
248:
249: p1=cgetg(1,t_VEC);
250: for (i=1; i < lp; i++)
251: {
252: p2 = mppgcd((GEN) primetab[i], prime);
253: if (! gcmp1(p2))
254: {
255: if (!egalii(prime,p2)) p1=concatsp(p1,p2);
256: p1 = concatsp(p1,divii((GEN)primetab[i],p2));
257: gunclone((GEN)primetab[i]); primetab[i]=0;
258: }
259: }
260: if (i == NUMPRTBELT+1 && lg(p1) == 1)
261: err(talker,"extra primetable overflows");
262: primetab[i] = lclone(prime); setlg(primetab, lp+1);
263: cleanprimetab();
264: if (lg(p1) > 1) addprimes(p1);
265: avma=av; return primetab;
266: }
267:
268: GEN
269: removeprimes(GEN prime)
270: {
271: long i,tx, av = avma;
272:
273: if (!prime) return primetab;
274: tx = typ(prime);
275: if (is_vec_t(tx))
276: {
277: for (i=1; i < lg(prime); i++) removeprimes((GEN) prime[i]);
278: return primetab;
279: }
280: if (tx != t_INT) err(typeer,"removeprimes");
281: if (is_pm1(prime)) return primetab;
282: if (signe(prime) < 0) prime=absi(prime);
283:
284: for (i=1; i < lg(primetab); i++)
285: if (egalii((GEN) primetab[i], prime))
286: {
287: gunclone((GEN)primetab[i]); primetab[i]=0;
288: cleanprimetab(); avma=av; return primetab;
289: }
290: err(talker,"prime is not in primetable in removeprimes");
291: return NULL; /* not reached */
292: }
293:
294: /***********************************************************************/
295: /** **/
296: /** COMPUTING THE MATRIX OF PRIME DIVISORS AND EXPONENTS **/
297: /** **/
298: /***********************************************************************/
299:
300: /* Configure may/should define this to 1 if MPQS is not installed --GN */
301: #ifndef decomp_default_hint
302: #define decomp_default_hint 0
303: #endif
304:
305: /* Some overkill removed from this (15 spsp for an integer < 2^32 ?!).
306: Should really revert to isprime() once the new primality tester is ready
307: --GN */
308: #define pseudoprime(p) millerrabin(p,3*lgefint(p))
309:
310: /* where to stop trial dividing in factorization */
311:
312: static long
313: tridiv_bound(GEN n, long all)
314: {
315: long size = expi(n) + 1;
316: if (all > 1) /* bounded factoring */
317: return all; /* use the given limit */
318: if (all == 0)
319: return VERYBIGINT; /* smallfact() case */
320: if (size <= 32)
321: return 16384;
322: else if (size <= 512)
323: return (size-16) << 10;
324: return 1L<<19; /* Rho will generally be faster above this */
325: }
326:
327: /********** about to become obsolete --GN **********/
328: #if 0
329: /* If flag != 1, use the new MPQS code: Try to factor N using ECM if flag = 2
330: * and N is not too small, followed by MPQS, or just MPQS otherwise.
331: */
332: static GEN
333: find_fact(GEN N, long *flag)
334: {
335: GEN f;
336:
337: if (*flag == 2)
338: {
339: if ((f = pollardbrent(N))) /* assignment */
340: return f;
341: f = ellfacteur(N, 0);
342: if (!f)
343: {
344: /* *flag = 0; */
345: f = mpqs(N);
346: }
347: }
348: else if (*flag == 1)
349: {
350: if ((f = pollardbrent(N))) /* assignment */
351: return f;
352: f = ellfacteur(N, 1);
353: }
354: else
355: f = mpqs(N); /* may bail out and call ellfacteur(_,1) */
356: /* (this is to be changed) */
357:
358: return f;
359: }
360: #endif
361: /***************************************************/
362:
363: /* function imported from ifactor1.c */
364: long
365: ifac_decomp(GEN n, long hint);
366:
367: static GEN
368: aux_end(GEN n, long nb)
369: {
370: GEN p,p1,p2, z = (GEN)avma;
371: long i;
372:
373: if (n) gunclone(n);
374: p1=cgetg(nb+1,t_COL);
375: p2=cgetg(nb+1,t_COL);
376: for (i=nb; i; i--)
377: {
378: p2[i] = (long)z; z += lg(z);
379: p1[i] = (long)z; z += lg(z);
380: }
381: z[1]=(long)p1;
382: z[2]=(long)p2;
383: if (nb)
384: {
385: long av = avma;
386: GEN p1old = dummycopy(p1), p2old = dummycopy(p2);
387: p=sindexsort(p1);
388: for (i=1;i<=nb; i++)
389: {
390: p1[i]=p1old[p[i]];
391: p2[i]=p2old[p[i]];
392: }
393: avma = av;
394: }
395: return z;
396: }
397:
398: static GEN
399: auxdecomp0(GEN n, long all, long hint)
400: {
401: long pp[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),2,0 };
402: long nb = 0,i,k,lim1,av,lp;
403: byteptr d=diffptr+1;
404:
405: if (typ(n) != t_INT) err(arither1);
406: i=signe(n); if (!i) err(arither2);
407: cgetg(3,t_MAT);
408: if (i<0) { stoi(-1); stoi(1); nb++; }
409: if (is_pm1(n)) return aux_end(NULL,nb);
410:
411: n = gclone(n); setsigne(n,1);
412: i = vali(n);
413: if (i)
414: {
415: stoi(2); stoi(i); nb++;
416: av=avma; affii(shifti(n,-i), n); avma=av;
417: }
418: if (is_pm1(n)) return aux_end(n,nb);
419: lim1 = tridiv_bound(n, all);
420:
421: /* trial division */
422: while (*d && pp[2] <= lim1)
423: {
424: pp[2] += *d++;
425: if (mpdivis(n,pp,n))
426: {
427: nb++; k=1; while (mpdivis(n,pp,n)) k++;
428: icopy(pp); stoi(k);
429: if (is_pm1(n)) return aux_end(n,nb);
430: }
431: }
432:
433: /* pp = square of biggest p tried so far */
434: av=avma; setlg(pp,4); affii(sqri(pp),pp); avma=av;
435: if (cmpii(pp,n) > 0 || pseudoprime(n))
436: {
437: nb++;
438: icopy(n); stoi(1); return aux_end(n,nb);
439: }
440:
441: lp = lg(primetab); k=0;
442: for (i=1; i<lp; i++)
443: if (mpdivis(n,(GEN) primetab[i],n))
444: {
445: nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++;
446: icopy((GEN) primetab[i]); stoi(k);
447: if (is_pm1(n)) return aux_end(n,nb);
448: }
449:
450: /* test primality again, _if_ primetab made a difference */
451: if(k && (cmpii(pp,n) > 0 || pseudoprime(n)))
452: {
453: nb++;
454: icopy(n); stoi(1); return aux_end(n,nb);
455: }
456:
457: /* now we have a large composite. Use hint as is if all==1 */
458: if (all!=1) hint = 15; /* turn off everything except pure powers */
459: nb += ifac_decomp(n, hint);
460:
461: return aux_end(n, nb);
462: }
463:
464: GEN
465: auxdecomp(GEN n, long all)
466: {
467: return auxdecomp0(n,all,decomp_default_hint);
468: }
469:
470: /* see before ifac_crack() in ifactor1.c for current semantics of `hint'
471: (factorint's `flag') */
472: GEN
473: factorint(GEN n, long flag)
474: {
475: return auxdecomp0(n,1,flag);
476: }
477:
478: GEN
479: decomp(GEN n)
480: {
481: return auxdecomp0(n,1,decomp_default_hint);
482: }
483:
484: GEN
485: smallfact(GEN n)
486: {
487: return boundfact(n,0);
488: }
489:
490: GEN
491: gboundfact(GEN n, long lim)
492: {
493: return garith_proto2gs(boundfact,n,lim);
494: }
495:
496: GEN
497: boundfact(GEN n, long lim)
498: {
499: GEN p1,p2,p3,p4,p5,y;
500: long av,tetpil;
501:
502: if (lim<=1) lim=0;
503: switch(typ(n))
504: {
505: case t_INT:
506: return auxdecomp(n,lim);
507: case t_FRACN:
508: av=avma; n=gred(n); /* fall through */
509: case t_FRAC:
510: if (typ(n)==t_FRAC) av=avma;
511: p1=auxdecomp((GEN)n[1],lim);
512: p2=auxdecomp((GEN)n[2],lim);
513: p4=concatsp((GEN)p1[1],(GEN)p2[1]);
514: p5=concatsp((GEN)p1[2],gneg((GEN)p2[2]));
515: p3=indexsort(p4);
516:
517: tetpil=avma; y=cgetg(3,t_MAT);
518: y[1]=(long)extract(p4,p3);
519: y[2]=(long)extract(p5,p3);
520: return gerepile(av,tetpil,y);
521: }
522: err(arither1);
523: return NULL; /* not reached */
524: }
525:
526: /***********************************************************************/
527: /** **/
528: /** BASIC ARITHMETIC FUNCTIONS **/
529: /** **/
530: /***********************************************************************/
531:
532: /* functions imported from ifactor1.c */
533: long
534: ifac_moebius(GEN n, long hint);
535:
536: long
537: ifac_issquarefree(GEN n, long hint);
538:
539: long
540: ifac_omega(GEN n, long hint);
541:
542: long
543: ifac_bigomega(GEN n, long hint);
544:
545: GEN
546: ifac_totient(GEN n, long hint);
547:
548: GEN
549: ifac_numdiv(GEN n, long hint);
550:
551: GEN
552: ifac_sumdiv(GEN n, long hint);
553:
554: GEN
555: ifac_sumdivk(GEN n, long k, long hint);
556:
557: GEN
558: gmu(GEN n)
559: {
560: return arith_proto(mu,n,1);
561: }
562:
563: long
564: mu(GEN n)
565: {
566: byteptr d = diffptr+1; /* point at 3 - 2 */
567: GEN p;
568: long av = avma, lim1, s, v;
569:
570: if (typ(n) != t_INT) err(arither1);
571: if (!signe(n)) err(arither2);
572: if (is_pm1(n)) return 1;
573: v = vali(n);
574: if (v>1) return 0;
575: s = v ? -1 : 1;
576: n=absi(shifti(n,-v)); p=court_p; p[2]=2;
577: if (is_pm1(n)) return s;
578: lim1 = tridiv_bound(n,1);
579:
580: while (*d && p[2] < lim1)
581: {
582: p[2] += *d++;
583: if (mpdivis(n,p,n))
584: {
585: if (divise(n,p)) { avma=av; return 0; }
586: s = -s; if (is_pm1(n)) { avma=av; return s; }
587: }
588: }
589: /* we normally get here with p==nextprime(PRE_TUNE), which will already
590: have been tested against n, so p^2 >= n (instead of p^2 > n) is
591: enough to guarantee n is prime */
592: if (cmpii(sqri(p),n) >= 0 || pseudoprime(n)) { avma=av; return -s; }
593: /* large composite without small factors */
594: v = ifac_moebius(n, decomp_default_hint);
595: avma=av;
596: return (s<0 ? -v : v); /* correct also if v==0 */
597: }
598:
599: GEN
600: gissquarefree(GEN x)
601: {
602: return arith_proto(issquarefree,x,0);
603: }
604:
605: long
606: issquarefree(GEN x)
607: {
608: long av=avma,lim1,tx;
609: GEN p;
610:
611: if (gcmp0(x)) return 0;
612: tx=typ(x);
613: if (tx == t_INT)
614: {
615: long v;
616: byteptr d=diffptr+1;
617: if (is_pm1(x)) return 1;
618: if((v = vali(x)) > 1) return 0;
619: x=absi(shifti(x,-v)); p=court_p; p[2]=2;
620: if (is_pm1(x)) return 1;
621: lim1 = tridiv_bound(x,1);
622:
623: while (*d && p[2] < lim1)
624: {
625: p[2] += *d++;
626: if (mpdivis(x,p,x))
627: {
628: if (divise(x,p)) { avma=av; return 0; }
629: if (is_pm1(x)) { avma=av; return 1; }
630: }
631: }
632: if (cmpii(sqri(p),x) >= 0 || pseudoprime(x)) { avma=av; return 1; }
633: v = ifac_issquarefree(x, decomp_default_hint);
634: avma=av;
635: return v;
636: }
637: if (tx != t_POL) err(typeer,"issquarefree");
638: p = ggcd(x, derivpol(x));
639: avma=av; return (lgef(p)==3);
640: }
641:
642: GEN
643: gomega(GEN n)
644: {
645: return arith_proto(omega,n,1);
646: }
647:
648: long
649: omega(GEN n)
650: {
651: byteptr d=diffptr+1;
652: GEN p;
653: long nb,av=avma,lim1,v;
654:
655: if (typ(n) != t_INT) err(arither1);
656: if (!signe(n)) err(arither2);
657: if (is_pm1(n)) return 0;
658: v=vali(n);
659: nb = v ? 1 : 0;
660: n=absi(shifti(n,-v)); p=court_p; p[2]=2;
661: if (is_pm1(n)) return nb;
662: lim1 = tridiv_bound(n,1);
663:
664: while (*d && p[2] < lim1)
665: {
666: p[2] += *d++;
667: if (mpdivis(n,p,n))
668: {
669: nb++; while (mpdivis(n,p,n)); /* empty */
670: if (is_pm1(n)) { avma = av; return nb; }
671: }
672: }
673: if (cmpii(sqri(p),n) >= 0 || pseudoprime(n)) { avma = av; return nb+1; }
674: /* large composite without small factors */
675: nb += ifac_omega(n, decomp_default_hint);
676: avma=av; return nb;
677: }
678:
679: GEN
680: gbigomega(GEN n)
681: {
682: return arith_proto(bigomega,n,1);
683: }
684:
685: long
686: bigomega(GEN n)
687: {
688: byteptr d=diffptr+1;
689: GEN p;
690: long nb,av=avma,lim1,v;
691:
692: if (typ(n) != t_INT) err(arither1);
693: if (!signe(n)) err(arither2);
694: if (is_pm1(n)) return 0;
695: nb=v=vali(n);
696: n=absi(shifti(n,-v)); p=court_p; p[2]=2;
697: if (is_pm1(n)) { avma = av; return nb; }
698: lim1 = tridiv_bound(n,1);
699:
700: while (*d && p[2] < lim1)
701: {
702: p[2] += *d++;
703: if (mpdivis(n,p,n))
704: {
705: do nb++; while (mpdivis(n,p,n));
706: if (is_pm1(n)) { avma = av; return nb; }
707: }
708: }
709: if (cmpii(sqri(p),n) >= 0 || pseudoprime(n)) { avma = av; return nb+1; }
710: nb += ifac_bigomega(n, decomp_default_hint);
711: avma=av; return nb;
712: }
713:
714: GEN
715: gphi(GEN n)
716: {
717: return garith_proto(phi,n,1);
718: }
719:
720: GEN
721: phi(GEN n)
722: {
723: byteptr d = diffptr+1;
724: GEN p,m;
725: long av = avma, lim1, v;
726:
727: if (typ(n) != t_INT) err(arither1);
728: if (!signe(n)) err(arither2);
729: if (is_pm1(n)) return gun;
730: v=vali(n);
731: n=absi(shifti(n,-v)); p=court_p; p[2]=2;
732: m = (v > 1 ? shifti(gun,v-1) : stoi(1));
733: if (is_pm1(n)) { return gerepileupto(av,m); }
734: lim1 = tridiv_bound(n,1);
735:
736: while (*d && p[2] < lim1)
737: {
738: court_p[2] += *d++;
739: if (mpdivis(n,p,n))
740: {
741: m = mulii(m,addsi(-1,p));
742: while (mpdivis(n,p,n)) m = mulii(m,p);
743: if (is_pm1(n)) { return gerepileupto(av,m); }
744: }
745: }
746: if (cmpii(sqri(p),n) >= 0 || pseudoprime(n))
747: {
748: m = mulii(m,addsi(-1,n));
749: return gerepileupto(av,m);
750: }
751: m = mulii(m, ifac_totient(n, decomp_default_hint));
752: return gerepileupto(av,m);
753: }
754:
755: GEN
756: gnumbdiv(GEN n)
757: {
758: return garith_proto(numbdiv,n,1);
759: }
760:
761: GEN
762: numbdiv(GEN n)
763: {
764: byteptr d=diffptr+1;
765: GEN p,m;
766: long l, av=avma, lim1, v;
767:
768: if (typ(n) != t_INT) err(arither1);
769: if (!signe(n)) err(arither2);
770: if (is_pm1(n)) return gun;
771: v=vali(n);
772: n=absi(shifti(n,-v)); p=court_p; p[2]=2;
773: m = stoi(v+1);
774: if (is_pm1(n)) { return gerepileupto(av,m); }
775: lim1 = tridiv_bound(n,1);
776:
777: while (*d && p[2] < lim1)
778: {
779: p[2] += *d++;
780: l=1; while (mpdivis(n,p,n)) l++;
781: m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); }
782: }
783: if(cmpii(sqri(p),n) >= 0 || pseudoprime(n))
784: {
785: m = shifti(m,1);
786: return gerepileupto(av,m);
787: }
788: m = mulii(m, ifac_numdiv(n, decomp_default_hint));
789: return gerepileupto(av,m);
790: }
791:
792: GEN
793: gsumdiv(GEN n)
794: {
795: return garith_proto(sumdiv,n,1);
796: }
797:
798: GEN
799: sumdiv(GEN n)
800: {
801: byteptr d=diffptr+1;
802: GEN p,m,m1;
803: long av=avma,lim1,v;
804:
805: if (typ(n) != t_INT) err(arither1);
806: if (!signe(n)) err(arither2);
807: if (is_pm1(n)) return gun;
808: v=vali(n);
809: n=absi(shifti(n,-v)); p=court_p; p[2]=2;
810: m = (v ? addsi(-1,shifti(gun,v+1)) : stoi(1));
811: if (is_pm1(n)) { return gerepileupto(av,m); }
812: lim1 = tridiv_bound(n,1);
813:
814: while (*d && p[2] < lim1)
815: {
816: p[2] += *d++;
817: if (mpdivis(n,p,n))
818: {
819: m1=addsi(1,p);
820: while (mpdivis(n,p,n)) m1=addsi(1,mulii(p,m1));
821: m = mulii(m1,m); if (is_pm1(n)) { return gerepileupto(av,m); }
822: }
823: }
824: if(cmpii(sqri(p),n) >= 0 || pseudoprime(n))
825: {
826: m = mulii(m,addsi(1,n));
827: return gerepileupto(av,m);
828: }
829: m = mulii(m, ifac_sumdiv(n, decomp_default_hint));
830: return gerepileupto(av,m);
831: }
832:
833: GEN
834: gsumdivk(GEN n, long k)
835: {
836: return garith_proto2gs(sumdivk,n,k);
837: }
838:
839: GEN
840: sumdivk(GEN n, long k)
841: {
842: byteptr d=diffptr+1;
843: GEN p,p1,m,m1,pk;
844: long av=avma,tetpil,k1,lim1,v;
845:
846: if (!k) return numbdiv(n);
847: if (k==1) return sumdiv(n);
848: if (typ(n) != t_INT) err(arither1);
849: if (!signe(n)) err(arither2);
850: if (is_pm1(n)) return gun;
851: k1=k;
852: if (k==-1) { p1=ginv(n); m=sumdiv(n); goto fin; }
853: if (k<0) { k= -k; p1=gpuigs(n,k1); }
854: v=vali(n);
855: n=absi(shifti(n,-v));
856: p = court_p; p[2]=2; m=stoi(1);
857: while (v--) { m=addsi(1,shifti(m,k)); }
858: if (is_pm1(n)) goto fin;
859: lim1 = tridiv_bound(n,1);
860:
861: while (*d && p[2] < lim1)
862: {
863: p[2] += *d++;
864: if (mpdivis(n,p,n))
865: {
866: pk=gpuigs(p,k); m1=addsi(1,pk);
867: while (mpdivis(n,p,n)) m1=addsi(1,mulii(pk,m1));
868: m = mulii(m1,m); if (is_pm1(n)) goto fin;
869: }
870: }
871: if(cmpii(sqri(p),n) >= 0 || pseudoprime(n))
872: {
873: pk=gpuigs(n,k);
874: m = mulii(m,addsi(1,pk));
875: goto fin;
876: }
877: m = mulii(m, ifac_sumdivk(n, k, decomp_default_hint));
878: fin:
879: if (k1>=0) return gerepileupto(av,m);
880: tetpil=avma; return gerepile(av,tetpil,gmul(p1,m));
881: }
882:
883: /***********************************************************************/
884: /** **/
885: /** MISCELLANEOUS ARITHMETIC FUNCTIONS **/
886: /** (all of these ultimately depend on auxdecomp()) **/
887: /** **/
888: /***********************************************************************/
889:
890: GEN
891: divisors(GEN n)
892: {
893: long tetpil,av=avma,i,j,l;
894: GEN *d,*t,*t1,*t2,*t3, nbdiv,e;
895:
896: if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1);
897:
898: e = (GEN) n[2], n = (GEN) n[1]; l = lg(n);
899: if (l>1 && signe(n[1]) < 0) { e++; n++; l--; } /* skip -1 */
900: nbdiv = gun;
901: for (i=1; i<l; i++)
902: {
903: e[i] = itos((GEN)e[i]);
904: nbdiv = mulis(nbdiv,1+e[i]);
905: }
906: d = t = (GEN*) cgetg(itos(nbdiv)+1,t_VEC);
907: *++d = gun;
908: for (i=1; i<l; i++)
909: for (t1=t,j=e[i]; j; j--,t1=t2)
910: for (t2=d,t3=t1; t3<t2; )
911: *++d = mulii(*++t3, (GEN)n[i]);
912: tetpil=avma; return gerepile(av,tetpil,sort((GEN)t));
913: }
914:
915: GEN
916: core(GEN n)
917: {
918: long av=avma,tetpil,i;
919: GEN fa,p1,p2,res=gun;
920:
921: fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2];
922: for (i=1; i<lg(p1); i++)
923: if (mod2((GEN)p2[i])) res = mulii(res,(GEN)p1[i]);
924: tetpil=avma; return gerepile(av,tetpil,icopy(res));
925: }
926:
927: GEN
928: core2(GEN n)
929: {
930: long av=avma,tetpil,i;
931:
932: GEN fa,p1,p2,p3,res=gun,res2=gun,y;
933:
934: fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2];
935: for (i=1; i<lg(p1); i++)
936: {
937: p3=(GEN)p2[i];
938: if (mod2(p3)) res=mulii(res,(GEN)p1[i]);
939: if (!gcmp1(p3)) res2=mulii(res2,gpui((GEN)p1[i],shifti(p3,-1),0));
940: }
941: tetpil=avma; y=cgetg(3,t_VEC);
942: y[1]=(long)icopy(res); y[2]=(long)icopy(res2);
943: return gerepile(av,tetpil,y);
944: }
945:
946: GEN
947: core0(GEN n,long flag)
948: {
949: return flag? core2(n): core(n);
950: }
951:
952: GEN
953: coredisc(GEN n)
954: {
955: long av=avma,tetpil,r;
956: GEN p1=core(n);
957:
958: r=mod4(p1); if (signe(p1)<0) r = 4-r;
959: if (r==1 || r==4) return p1;
960: tetpil=avma; return gerepile(av,tetpil,shifti(p1,2));
961: }
962:
963: GEN
964: coredisc2(GEN n)
965: {
966: long av=avma,tetpil,r;
967: GEN y,p1,p2=core2(n);
968:
969: p1=(GEN)p2[1]; r=mod4(p1); if (signe(p1)<0) r=4-r;
970: if (r==1 || r==4) return p2;
971: tetpil=avma; y=cgetg(3,t_VEC);
972: y[1]=lshifti(p1,2); y[2]=lmul2n((GEN)p2[2],-1);
973: return gerepile(av,tetpil,y);
974: }
975:
976: GEN
977: coredisc0(GEN n,long flag)
978: {
979: return flag? coredisc2(n): coredisc(n);
980: }
981:
982: /***********************************************************************/
983: /** **/
984: /** BINARY QUADRATIC FORMS **/
985: /** **/
986: /***********************************************************************/
987:
988: GEN
989: qf_disc(GEN x, GEN y, GEN z)
990: {
991: if (!y) { y=(GEN)x[2]; z=(GEN)x[3]; x=(GEN)x[1]; }
992: return subii(sqri(y), shifti(mulii(x,z),2));
993: }
994:
995: static GEN
996: qf_create(GEN x, GEN y, GEN z, long s)
997: {
998: GEN t;
999: if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb");
1000: if (!s)
1001: {
1002: long av=avma; s = signe(qf_disc(x,y,z)); avma=av;
1003: if (!s) err(talker,"zero discriminant in Qfb");
1004: }
1005: t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI);
1006: t[1]=(long)icopy(x); t[2]=(long)icopy(y); t[3]=(long)icopy(z);
1007: return t;
1008: }
1009:
1010: GEN
1011: qfi(GEN x, GEN y, GEN z)
1012: {
1013: return qf_create(x,y,z,-1);
1014: }
1015:
1016: GEN
1017: qfr(GEN x, GEN y, GEN z, GEN d)
1018: {
1019: GEN t = qf_create(x,y,z,1);
1020: if (typ(d) != t_REAL)
1021: err(talker,"Shanks distance should be a t_REAL (in qfr)");
1022: t[4]=lrcopy(d); return t;
1023: }
1024:
1025: GEN
1026: Qfb0(GEN x, GEN y, GEN z, GEN d, long prec)
1027: {
1028: GEN t = qf_create(x,y,z,0);
1029: if (lg(t)==4) return t;
1030: if (!d) d = gzero;
1031: if (typ(d) == t_REAL)
1032: t[4] = lrcopy(d);
1033: else
1034: { t[4]=lgetr(prec); gaffect(d,(GEN)t[4]); }
1035: return t;
1036: }
1037:
1038: /* composition */
1039:
1040: static void
1041: comp_gen(GEN z,GEN x,GEN y)
1042: {
1043: GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,v3,b3,c3,m,p1,r;
1044:
1045: s=shifti(addii((GEN)x[2],(GEN)y[2]),-1);
1046: n=subii((GEN)y[2],s);
1047: d = bezout((GEN)y[1],(GEN)x[1],&y1,&x1);
1048: d1 = bezout(s,d,&x2,&y2);
1049: if (gcmp1(d1))
1050: {
1051: v3 = (GEN)x[1];
1052: v2 = (GEN)y[1];
1053: }
1054: else
1055: {
1056: v1=divii((GEN)x[1],d1);
1057: v2=divii((GEN)y[1],d1);
1058: v3 = mulii(v1, mppgcd(d1,mppgcd((GEN)x[3],mppgcd((GEN)y[3],n))));
1059: }
1060: m = addii(mulii(mulii(y1,y2),n), mulii((GEN)y[3],x2));
1061: setsigne(m,-signe(m));
1062: r=modii(m,v3); p1=mulii(v2,r); b3=shifti(p1,1);
1063: c3=addii(mulii((GEN)y[3],d1),mulii(r,addii((GEN)y[2],p1)));
1064: z[1]=lmulii(v3,v2);
1065: z[2]=laddii((GEN)y[2],b3);
1066: z[3]=ldivii(c3,v3);
1067: }
1068:
1069: static GEN
1070: compimag0(GEN x, GEN y, int raw)
1071: {
1072: long tx = typ(x), av = avma, tetpil;
1073: GEN z;
1074:
1075: if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition");
1076: if (cmpii((GEN)x[1],(GEN)y[1]) > 0) { z=x; x=y; y=z; }
1077: z=cgetg(4,t_QFI); comp_gen(z,x,y); tetpil=avma;
1078: return gerepile(av,tetpil,raw? gcopy(z): redimag(z));
1079: }
1080:
1081: static GEN
1082: compreal0(GEN x, GEN y, int raw)
1083: {
1084: long tx = typ(x), av = avma, tetpil;
1085: GEN z;
1086:
1087: if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition");
1088: z=cgetg(5,t_QFR); comp_gen(z,x,y);
1089: z[4]=laddrr((GEN)x[4],(GEN)y[4]); tetpil=avma;
1090: return gerepile(av,tetpil, raw? gcopy(z): redreal(z));
1091: }
1092:
1093: GEN
1094: compreal(GEN x, GEN y) { return compreal0(x,y,0); }
1095:
1096: GEN
1097: comprealraw(GEN x, GEN y) { return compreal0(x,y,1); }
1098:
1099: GEN
1100: compimag(GEN x, GEN y) { return compimag0(x,y,0); }
1101:
1102: GEN
1103: compimagraw(GEN x, GEN y) { return compimag0(x,y,1); }
1104:
1105: GEN
1106: compraw(GEN x, GEN y)
1107: {
1108: return (typ(x)==t_QFI)? compimagraw(x,y): comprealraw(x,y);
1109: }
1110:
1111: static void
1112: sq_gen(GEN z, GEN x)
1113: {
1114: GEN d1,x2,y2,v1,v3,b3,c3,m,p1,r;
1115:
1116: d1 = bezout((GEN)x[2],(GEN)x[1],&x2,&y2);
1117: if (gcmp1(d1)) v1 = v3 = (GEN)x[1];
1118: else
1119: {
1120: v1 = divii((GEN)x[1],d1);
1121: v3 = mulii(v1,mppgcd(d1,(GEN)x[3]));
1122: }
1123: m=mulii((GEN)x[3],x2); setsigne(m,-signe(m));
1124: r=modii(m,v3); p1=mulii(v1,r); b3=shifti(p1,1);
1125: c3=addii(mulii((GEN)x[3],d1),mulii(r,addii((GEN)x[2],p1)));
1126: z[1]=lmulii(v3,v1);
1127: z[2]=laddii((GEN)x[2],b3);
1128: z[3]=ldivii(c3,v3);
1129: }
1130:
1131: GEN
1132: sqcompimag0(GEN x, long raw)
1133: {
1134: long av = avma, tetpil;
1135: GEN z = cgetg(4,t_QFI);
1136:
1137: if (typ(x)!=t_QFI) err(typeer,"composition");
1138: sq_gen(z,x); tetpil=avma;
1139: return gerepile(av,tetpil,raw? gcopy(z): redimag(z));
1140: }
1141:
1142: static GEN
1143: sqcompreal0(GEN x, long raw)
1144: {
1145: long av = avma, tetpil;
1146: GEN z = cgetg(5,t_QFR);
1147:
1148: if (typ(x)!=t_QFR) err(typeer,"composition");
1149: sq_gen(z,x); z[4]=lshiftr((GEN)x[4],1); tetpil=avma;
1150: return gerepile(av,tetpil,raw? gcopy(z): redreal(z));
1151: }
1152:
1153: GEN
1154: sqcompreal(GEN x) { return sqcompreal0(x,0); }
1155:
1156: GEN
1157: sqcomprealraw(GEN x) { return sqcompreal0(x,1); }
1158:
1159: GEN
1160: sqcompimag(GEN x) { return sqcompimag0(x,0); }
1161:
1162: GEN
1163: sqcompimagraw(GEN x) { return sqcompimag0(x,1); }
1164:
1165: GEN
1166: real_unit_form(GEN x)
1167: {
1168: long av = avma, tetpil,prec;
1169: GEN p1,p2, y = cgetg(5,t_QFR);
1170:
1171: y[1]=un; p1=qf_disc(x,NULL,NULL); p2=racine(p1);
1172: /* we know that p1 and p2 are non-zero */
1173: if (mod2(p1) != mod2(p2)) p2=addsi(-1,p2);
1174: y[2]=(long)p2;
1175: y[3]=lshifti(subii(sqri(p2),p1),-2);
1176: prec = precision((GEN)x[4]);
1177: if (!prec)
1178: err(talker,"not a type real in 4th component of a t_QFR");
1179: y[4]=(long)realzero(prec); tetpil=avma;
1180: return gerepile(av,tetpil,gcopy(y));
1181: }
1182:
1183: GEN
1184: imag_unit_form(GEN x)
1185: {
1186: long av,tetpil;
1187: GEN p1,p2, y = cgetg(4,t_QFI);
1188: y[1]=un;
1189: y[2]=mpodd((GEN)x[2])? un: zero;
1190: av=avma; p1=mulii((GEN)x[1],(GEN)x[3]);
1191: p2=shifti(sqri((GEN)x[2]),-2); tetpil=avma;
1192: y[3]=lpile(av,tetpil,subii(p1,p2));
1193: return y;
1194: }
1195:
1196: GEN
1197: powrealraw(GEN x, long n)
1198: {
1199: long av = avma, m;
1200: GEN y;
1201:
1202: if (typ(x) != t_QFR)
1203: err(talker,"not a real quadratic form in powreal");
1204: if (!n) return real_unit_form(x);
1205: if (n== 1) return gcopy(x);
1206: if (n==-1) return ginv(x);
1207:
1208: y = NULL; m=labs(n);
1209: for (; m>1; m>>=1)
1210: {
1211: if (m&1) y = y? comprealraw(y,x): x;
1212: x=sqcomprealraw(x);
1213: }
1214: y = y? comprealraw(y,x): x;
1215: if (n<0) y = ginv(y);
1216: return gerepileupto(av,y);
1217: }
1218:
1219: GEN
1220: powimagraw(GEN x, long n)
1221: {
1222: long av = avma, m;
1223: GEN y;
1224:
1225: if (typ(x) != t_QFI)
1226: err(talker,"not an imaginary quadratic form in powimag");
1227: if (!n) return imag_unit_form(x);
1228: if (n== 1) return gcopy(x);
1229: if (n==-1) return ginv(x);
1230:
1231: y = NULL; m=labs(n);
1232: for (; m>1; m>>=1)
1233: {
1234: if (m&1) y = y? compimagraw(y,x): x;
1235: x=sqcompimagraw(x);
1236: }
1237: y = y? compimagraw(y,x): x;
1238: if (n<0) y = ginv(y);
1239: return gerepileupto(av,y);
1240: }
1241:
1242: GEN
1243: powraw(GEN x, long n)
1244: {
1245: return (typ(x)==t_QFI)? powimagraw(x,n): powrealraw(x,n);
1246: }
1247:
1248: /* composition: Shanks' NUCOMP & NUDUPL */
1249: /* l = floor((|d|/4)^(1/4)) */
1250: GEN
1251: nucomp(GEN x, GEN y, GEN l)
1252: {
1253: long av=avma,tetpil,cz;
1254: 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;
1255:
1256: if (x==y) return nudupl(x,l);
1257: if (typ(x) != t_QFI || typ(y) != t_QFI)
1258: err(talker,"not an imaginary quadratic form in nucomp");
1259:
1260: if (cmpii((GEN)x[1],(GEN)y[1]) < 0) { s=x; x=y; y=s; }
1261: s=shifti(addii((GEN)x[2],(GEN)y[2]),-1); n=subii((GEN)y[2],s);
1262: a1=(GEN)x[1]; a2=(GEN)y[1]; d=bezout(a2,a1,&u,&v);
1263: if (gcmp1(d)) { a=negi(gmul(u,n)); d1=d; }
1264: else
1265: if (divise(s,d))
1266: {
1267: a=negi(gmul(u,n)); d1=d; a1=divii(a1,d1);
1268: a2=divii(a2,d1); s=divii(s,d1);
1269: }
1270: else
1271: {
1272: d1=bezout(s,d,&u1,&v1);
1273: if (!gcmp1(d1))
1274: {
1275: a1=divii(a1,d1); a2=divii(a2,d1);
1276: s=divii(s,d1); d=divii(d,d1);
1277: }
1278: p1=resii((GEN)x[3],d); p2=resii((GEN)y[3],d);
1279: p3=modii(negi(mulii(u1,addii(mulii(u,p1),mulii(v,p2)))),d);
1280: a=subii(mulii(p3,divii(a1,d)),mulii(u,divii(n,d)));
1281: }
1282: a=modii(a,a1); p1=subii(a1,a); if (cmpii(a,p1)>0) a=negi(p1);
1283: v=gzero; d=a1; v2=gun; v3=a;
1284: for (cz=0; absi_cmp(v3,l) > 0; cz++)
1285: {
1286: p1=dvmdii(d,v3,&t3); t2=subii(v,mulii(p1,v2));
1287: v=v2; d=v3; v2=t2; v3=t3;
1288: }
1289: z=cgetg(4,t_QFI);
1290: if (!cz)
1291: {
1292: q1=mulii(a2,v3); q2=addii(q1,n);
1293: g=divii(addii(mulii(v3,s),(GEN)y[3]),d);
1294: z[1]=lmulii(d,a2);
1295: z[2]=laddii(shifti(q1,1),(GEN)y[2]);
1296: z[3]=laddii(mulii(v3,divii(q2,d)), mulii(g,d1));
1297: }
1298: else
1299: {
1300: if (cz&1) { v3=negi(v3); v2=negi(v2); }
1301: b=divii(addii(mulii(a2,d),mulii(n,v)),a1);
1302: q1=mulii(b,v3); q2=addii(q1,n);
1303: e=divii(addii(mulii(s,d),mulii((GEN)y[3],v)),a1);
1304: q3=mulii(e,v2); q4=subii(q3,s);
1305: g=divii(q4,v); b2=addii(q3,q4);
1306: if (!gcmp1(d1)) { v2=mulii(d1,v2); v=mulii(d1,v); b2=mulii(d1,b2); }
1307: z[1]=laddii(mulii(d,b),mulii(e,v));
1308: z[2]=laddii(b2,addii(q1,q2));
1309: z[3]=laddii(mulii(v3,divii(q2,d)), mulii(g,v2));
1310: }
1311: tetpil=avma; return gerepile(av,tetpil,redimag(z));
1312: }
1313:
1314: GEN
1315: nudupl(GEN x, GEN l)
1316: {
1317: long av=avma,tetpil,cz;
1318: GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g;
1319:
1320: if (typ(x) != t_QFI)
1321: err(talker,"not an imaginary quadratic form in nudupl");
1322: d1=bezout((GEN)x[2],(GEN)x[1],&u,&v);
1323: a=divii((GEN)x[1],d1); b=divii((GEN)x[2],d1);
1324: c=modii(negi(mulii(u,(GEN)x[3])),a); p1=subii(a,c);
1325: if (cmpii(c,p1)>0) c=negi(p1);
1326: v=gzero; d=a; v2=gun; v3=c;
1327: for (cz=0; absi_cmp(v3,l) > 0; cz++)
1328: {
1329: p1=dvmdii(d,v3,&t3); t2=subii(v,mulii(p1,v2));
1330: v=v2; d=v3; v2=t2; v3=t3;
1331: }
1332: z=cgetg(4,t_QFI);
1333: if (!cz)
1334: {
1335: g=divii(addii(mulii(v3,b),(GEN)x[3]),d);
1336: z[1]=(long)sqri(d);
1337: z[2]=laddii((GEN)x[2],shifti(mulii(d,v3),1));
1338: z[3]=laddii(sqri(v3),mulii(g,d1));
1339: }
1340: else
1341: {
1342: if (cz&1) { v=negi(v); d=negi(d); }
1343: e=divii(addii(mulii((GEN)x[3],v),mulii(b,d)),a);
1344: g=divii(subii(mulii(e,v2),b),v);
1345: b2=addii(mulii(e,v2),mulii(v,g));
1346: if (!gcmp1(d1)) { v2=mulii(d1,v2); v=mulii(d1,v); b2=mulii(d1,b2); }
1347: z[1]=laddii(sqri(d),mulii(e,v));
1348: z[2]=laddii(b2,shifti(mulii(d,v3),1));
1349: z[3]=laddii(sqri(v3),mulii(g,v2));
1350: }
1351: tetpil=avma; return gerepile(av,tetpil,redimag(z));
1352: }
1353:
1354: GEN
1355: nupow(GEN x, GEN n)
1356: {
1357: long av,tetpil,i,j;
1358: unsigned long m;
1359: GEN y,l;
1360:
1361: if (typ(n) != t_INT) err(talker,"not an integer exponent in nupow");
1362: if (gcmp1(n)) return gcopy(x);
1363: av=avma; y = imag_unit_form(x);
1364: if (!signe(n)) return y;
1365:
1366: l = racine(shifti(racine((GEN)y[3]),1));
1367: for (i=lgefint(n)-1; i>2; i--)
1368: for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
1369: {
1370: if (m&1) y=nucomp(y,x,l);
1371: x=nudupl(x,l);
1372: }
1373: for (m=n[2]; m>1; m>>=1)
1374: {
1375: if (m&1) y=nucomp(y,x,l);
1376: x=nudupl(x,l);
1377: }
1378: tetpil=avma; y=nucomp(y,x,l);
1379: if (signe(n)<0 && !egalii((GEN)y[1],(GEN)y[2])
1380: && !egalii((GEN)y[1],(GEN)y[3]))
1381: setsigne(y[2],-signe(y[2]));
1382: return gerepile(av,tetpil,y);
1383: }
1384:
1385: /* reduction */
1386:
1387: static GEN
1388: abs_dvmdii(GEN b, GEN p1, GEN *pt, long s)
1389: {
1390: if (s<0) setsigne(b, 1); p1 = dvmdii(b,p1,pt);
1391: if (s<0) setsigne(b,-1); return p1;
1392: }
1393:
1394: static GEN
1395: rhoimag0(GEN x, long *flag)
1396: {
1397: GEN p1,b,d,z;
1398: long fl, s = signe(x[2]);
1399:
1400: fl = cmpii((GEN)x[1], (GEN)x[3]);
1401: if (fl <= 0)
1402: {
1403: long fg = absi_cmp((GEN)x[1], (GEN)x[2]);
1404: if (fg >= 0)
1405: {
1406: *flag = (s<0 && (!fl || !fg))? 2 /* set x[2] = negi(x[2]) in caller */
1407: : 1;
1408: return x;
1409: }
1410: }
1411: p1=shifti((GEN)x[3],1); d = abs_dvmdii((GEN)x[2],p1,&b,s);
1412: if (s>=0)
1413: {
1414: setsigne(d,-signe(d));
1415: if (cmpii(b,(GEN)x[3])<=0) setsigne(b,-signe(b));
1416: else { d=addsi(-1,d); b=subii(p1,b); }
1417: }
1418: else if (cmpii(b,(GEN)x[3])>=0) { d=addsi(1,d); b=subii(b,p1); }
1419:
1420: z=cgetg(4,t_QFI);
1421: z[1] = x[3];
1422: z[3] = laddii((GEN)x[1], mulii(d,shifti(subii((GEN)x[2],b),-1)));
1423: if (signe(b)<0 && !fl) setsigne(b,1);
1424: z[2] = (long)b; *flag = 0; return z;
1425: }
1426:
1427: /* if sqrtD != NULL, compute Shanks' distance as well */
1428: static GEN
1429: rhoreal0(GEN x, GEN D, GEN isqrtD, GEN sqrtD)
1430: {
1431: long s = signe(x[3]);
1432: GEN p1,p2, z = cgetg(5,t_QFR);
1433:
1434: if (!sqrtD)
1435: z[4] = x[4];
1436: else
1437: {
1438: p1 = divrr(addir((GEN)x[2],sqrtD), subir((GEN)x[2],sqrtD));
1439: if (signe(p1)<0) setsigne(p1,1); /* p1 = |(b+sqrtD) / (b-sqrtD)| */
1440: z[4] = laddrr(shiftr(mplog(p1),-1), (GEN) x[4]);
1441: }
1442: z[1] = x[3]; setsigne(z[1],1);
1443: p1=shifti((GEN)z[1],1);
1444: p2 = (cmpii(isqrtD,(GEN)z[1]) >= 0)? isqrtD: (GEN)z[1];
1445: p2 = divii(addii(p2,(GEN)x[2]),p1);
1446: z[2] = lsubii(mulii(p2,p1), (GEN)x[2]);
1447:
1448: setsigne(z[1],s);
1449: p1=shifti(subii(sqri((GEN)z[2]),D),-2);
1450: z[3] = ldivii(p1,(GEN)z[1]); return z;
1451: }
1452:
1453: static GEN
1454: redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
1455: {
1456: long av=avma,l,step, use_d;
1457:
1458: if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal");
1459: switch(flag)
1460: {
1461: case 0: use_d=1; step=0; break;
1462: case 1: use_d=1; step=1; break;
1463: case 2: use_d=0; step=0; break;
1464: case 3: use_d=0; step=1; break;
1465: default: err(flagerr,"qfbred");
1466: }
1467:
1468: if (!D) D = qf_disc(x,NULL,NULL);
1469: else if (typ(D) != t_INT) err(arither1);
1470:
1471: if (!isqrtD || (use_d && !sqrtD))
1472: {
1473: l=max(lg(x[4]),3);
1474: l=max(l,((BITS_IN_LONG-1-expo(x[4]))>>TWOPOTBITS_IN_LONG)+2);
1475: sqrtD=gsqrt(D,l); isqrtD=mptrunc(sqrtD);
1476: }
1477: else
1478: {
1479: if (typ(isqrtD) != t_INT) err(arither1);
1480: if (sqrtD)
1481: {
1482: l=typ(sqrtD); if (!is_intreal_t(l)) err(arither1);
1483: }
1484: }
1485: if (step)
1486: x = rhoreal0(x,D,isqrtD,sqrtD);
1487: else
1488: for(;;)
1489: {
1490: if (signe(x[2]) > 0 && cmpii((GEN)x[2],isqrtD) <= 0 )
1491: {
1492: GEN p1=subii(isqrtD, shifti(absi((GEN)x[1]),1));
1493:
1494: l = absi_cmp((GEN)x[2],p1);
1495: if (l>0 || (l==0 && signe(p1)<0)) break;
1496: }
1497: x = rhoreal0(x,D,isqrtD,sqrtD);
1498: }
1499: l=avma; return gerepile(av,l,gcopy(x));
1500: }
1501:
1502: GEN
1503: redimag(GEN x)
1504: {
1505: long av=avma, tetpil, fl;
1506: do x = rhoimag0(x, &fl); while (fl == 0);
1507: tetpil=avma; x = gerepile(av,tetpil,gcopy(x));
1508: if (fl == 2) setsigne(x[2], -signe(x[2]));
1509: return x;
1510: }
1511:
1512: GEN
1513: redreal(GEN x)
1514: {
1515: return redreal0(x,0,NULL,NULL,NULL);
1516: }
1517:
1518: GEN
1519: rhoreal(GEN x)
1520: {
1521: return redreal0(x,1,NULL,NULL,NULL);
1522: }
1523:
1524: GEN
1525: redrealnod(GEN x, GEN isqrtD)
1526: {
1527: return redreal0(x,2,NULL,isqrtD,NULL);
1528: }
1529:
1530: GEN
1531: rhorealnod(GEN x, GEN isqrtD)
1532: {
1533: return redreal0(x,3,NULL,isqrtD,NULL);
1534: }
1535:
1536: GEN
1537: qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
1538: {
1539: long tx=typ(x),av,tetpil,fl;
1540:
1541: if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred");
1542: if (!D) D = qf_disc(x,NULL,NULL);
1543: switch(signe(D))
1544: {
1545: case 1 :
1546: return redreal0(x,flag,D,isqrtD,sqrtD);
1547:
1548: case -1:
1549: if (!flag) return redimag(x);
1550: if (flag !=1) err(flagerr,"qfbred");
1551: av=avma; x = rhoimag0(x,&fl); tetpil=avma;
1552: x = gerepile(av,tetpil,gcopy(x));
1553: if (fl == 2) setsigne(x[2], -signe(x[2]));
1554: return x;
1555: }
1556: err(redpoler,"qfbred");
1557: return NULL; /* not reached */
1558: }
1559:
1560: GEN
1561: primeform(GEN x, GEN p, long prec)
1562: {
1563: long av,tetpil,s=signe(x);
1564: GEN y,b;
1565:
1566: if (typ(x) != t_INT || !s) err(arither1);
1567: if (typ(p) != t_INT || signe(p) <= 0) err(arither1);
1568: if (s < 0)
1569: {
1570: y = cgetg(4,t_QFI); s = 8-mod8(x);
1571: }
1572: else
1573: {
1574: y = cgetg(5,t_QFR); s = mod8(x);
1575: y[4]=(long)realzero(prec);
1576: }
1577: switch(s&3)
1578: {
1579: case 2: case 3: err(funder2,"primeform");
1580: }
1581: y[1] = (long)icopy(p); av=avma;
1582: if (egalii(p,gdeux))
1583: {
1584: switch(s)
1585: {
1586: case 0: y[2]=zero; break;
1587: case 8: s=0; y[2]=zero; break;
1588: case 1: s=1; y[2]=un; break;
1589: case 4: s=4; y[2]=deux; break;
1590: default: err(sqrter5);
1591: }
1592: b=subsi(s,x); tetpil=avma; b=shifti(b,-3);
1593: }
1594: else
1595: {
1596: b = mpsqrtmod(x,p); if (!b) err(sqrter5);
1597: if (mod2(b) == mod2(x)) y[2] = (long)b;
1598: else
1599: { tetpil = avma; y[2] = lpile(av,tetpil,subii(p,b)); }
1600:
1601: av=avma; b=shifti(subii(sqri((GEN)y[2]),x),-2);
1602: tetpil=avma; b=divii(b,p);
1603: }
1604: y[3]=lpile(av,tetpil,b); return y;
1605: }
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 / src / basemath