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