Annotation of OpenXM_contrib/pari/src/basemath/arith2.c, Revision 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