=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/arith2.c,v retrieving revision 1.1 retrieving revision 1.2 diff -u -p -r1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/basemath/Attic/arith2.c 2001/10/02 11:17:01 1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/arith2.c 2002/09/11 07:26:48 1.2 @@ -1,4 +1,4 @@ -/* $Id: arith2.c,v 1.1 2001/10/02 11:17:01 noro Exp $ +/* $Id: arith2.c,v 1.2 2002/09/11 07:26:48 noro Exp $ Copyright (C) 2000 The PARI group. @@ -25,6 +25,7 @@ extern GEN arith_proto(long f(GEN), GEN x, int do_erro extern GEN arith_proto2(long f(GEN,GEN), GEN x, GEN n); extern GEN garith_proto(GEN f(GEN), GEN x, int do_error); extern GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y); +extern GEN garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z); #define sqru(i) (muluu((i),(i))) @@ -38,13 +39,12 @@ GEN prime(long n) { byteptr p = diffptr; - long c, prime = 0; + long prime = 0; if (n <= 0) err(talker, "n-th prime meaningless if n = %ld",n); while (n--) { - c = *p++; if (!c) err(primer1); - prime += c; + NEXT_PRIME_VIADIFF_CHECK(prime,p); } return stoi(prime); } @@ -56,8 +56,11 @@ pith(long n) ulong prime = 0, res = 0; if (n <= 0) err(talker, "pith meaningless if n = %ld",n); - if (maxprime() <= n) err(primer1); - while (prime<=n) { prime += *p++; res++; } + if (maxprime() <= (ulong)n) err(primer1); + while (prime<=(ulong)n) { + NEXT_PRIME_VIADIFF(prime,p); + res++; + } return stoi(res-1); } @@ -65,15 +68,15 @@ GEN primes(long n) { byteptr p = diffptr; - long c, prime = 0; + long prime = 0; GEN y,z; if (n < 0) n = 0; z = y = cgetg(n+1,t_VEC); while (n--) { - c = *p++; if (!c) err(primer1); - prime += c; *++z = lstoi(prime); + NEXT_PRIME_VIADIFF_CHECK(prime,p); + *++z = lstoi(prime); } return y; } @@ -85,7 +88,7 @@ primes(long n) * when maxnum (size) is moderate. */ static byteptr -initprimes1(long size, long *lenp, long *lastp) +initprimes1(ulong size, long *lenp, long *lastp) { long k; byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2); @@ -110,7 +113,7 @@ initprimes1(long size, long *lenp, long *lastp) fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n", q, *lenp, *lastp); #endif - return (byteptr) gprealloc(p,r-p,size + 2); + return (byteptr) gprealloc(p,r-p); } #define PRIME_ARENA (200 * 1024) /* No slowdown even with 256K level-2 cache */ @@ -119,10 +122,11 @@ initprimes1(long size, long *lenp, long *lastp) recursive call will bottom out and invoke initprimes1() at once. (Not static; might conceivably be useful to someone in library mode) */ byteptr -initprimes0(ulong maxnum, long *lenp, long *lastp) +initprimes0(ulong maxnum, long *lenp, ulong *lastp) { - long k, size, alloced, asize, psize, rootnum, curlow, last, remains; + long k, size, alloced, asize, psize, rootnum; byteptr q,r,s,fin, p, p1, fin1, plast, curdiff; + ulong last, remains, curlow; #if 0 fprintferr("initprimes0: called for maxnum = %lu\n", maxnum); @@ -160,7 +164,7 @@ initprimes0(ulong maxnum, long *lenp, long *lastp) asize = ARENA_IN_ROOTS * rootnum; /* Make % overhead negligeable. */ if (asize < PRIME_ARENA) asize = PRIME_ARENA; if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */ - alloced = (avma - bot < asize>>1); /* enough room on the stack ? */ + alloced = (avma - bot < (size_t)asize>>1); /* enough room on the stack ? */ if (alloced) p = (byteptr) gpmalloc(asize); else @@ -182,6 +186,7 @@ initprimes0(ulong maxnum, long *lenp, long *lastp) } memset(p, 0, asize); /* p corresponds to curlow. q runs over primediffs. */ + /* Don't care about DIFFPTR_SKIP: false positives provide no problem */ for (q = p1+2, k = 3; q <= fin1; k += *q++) { /* The first odd number which is >= curlow and divisible by p @@ -207,9 +212,14 @@ initprimes0(ulong maxnum, long *lenp, long *lastp) /* now q runs over addresses corresponding to primes */ for (q = p; ; plast = q++) { + long d; + while (*q) q++; /* use 0-sentinel at end */ if (q >= fin) break; - *curdiff++ = (q - plast) << 1; + d = (q - plast) << 1; + while (d >= DIFFPTR_SKIP) + *curdiff++ = DIFFPTR_SKIP, d -= DIFFPTR_SKIP; + *curdiff++ = d; } plast -= asize - 1; remains -= asize - 1; @@ -220,7 +230,7 @@ initprimes0(ulong maxnum, long *lenp, long *lastp) *lenp = curdiff - p1; *lastp = last; if (alloced) free(p); - return (byteptr) gprealloc(p1, *lenp, size); + return (byteptr) gprealloc(p1, *lenp); } #if 0 /* not yet... GN */ /* The diffptr table will contain at least 6548 entries (up to and including @@ -235,20 +245,20 @@ init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */ static ulong _maxprime = 0; ulong -maxprime() { return _maxprime; } +maxprime(void) { return _maxprime; } byteptr -initprimes(long maxnum) +initprimes(ulong maxnum) { - long len, last; + long len; + ulong last; byteptr p; - /* The algorithm must see the next prime beyond maxnum, whence the +257. */ - /* switch to unsigned: adding the 257 _could_ wrap to a negative number. */ - ulong maxnum1 = ((maxnum<65302?65302:maxnum)+257ul); + /* The algorithm must see the next prime beyond maxnum, whence the +512. */ + /* switch to unsigned: adding the 512 _could_ wrap to a negative number. */ + ulong maxnum1 = ((maxnum<65302?65302:maxnum)+512ul); - if (maxnum1 > 436273000) - err(talker,"impossible to have prestored primes > 436273009"); - + if ((maxnum>>1) > VERYBIGINT - 1024) + err(talker, "Too large primelimit"); p = initprimes0(maxnum1, &len, &last); #if 0 /* not yet... GN */ static int build_the_tables = 1; @@ -273,7 +283,7 @@ cleanprimetab(void) GEN addprimes(GEN p) { - ulong av; + gpmem_t av; long i,k,tx,lp; GEN L; @@ -281,7 +291,7 @@ addprimes(GEN p) tx = typ(p); if (is_vec_t(tx)) { - for (i=1; i < lg(p); i++) addprimes((GEN) p[i]); + for (i=1; i < lg(p); i++) (void)addprimes((GEN) p[i]); return primetab; } if (tx != t_INT) err(typeer,"addprime"); @@ -302,9 +312,9 @@ addprimes(GEN p) gunclone(n); primetab[i] = 0; } } - primetab = (GEN) gprealloc(primetab, (lp+1)*sizeof(long), lp*sizeof(long)); + primetab = (GEN) gprealloc(primetab, (lp+1)*sizeof(long)); primetab[i] = lclone(p); setlg(primetab, lp+1); - if (k > 1) { cleanprimetab(); setlg(L,k); addprimes(L); } + if (k > 1) { cleanprimetab(); setlg(L,k); (void)addprimes(L); } avma = av; return primetab; } @@ -340,7 +350,7 @@ removeprimes(GEN prime) } else { - for (i=1; i < lg(prime); i++) removeprime((GEN) prime[i]); + for (i=1; i < lg(prime); i++) (void)removeprime((GEN) prime[i]); } return primetab; } @@ -361,7 +371,11 @@ removeprimes(GEN prime) /* Some overkill removed from this (15 spsp for an integer < 2^32 ?!). Should really revert to isprime() once the new primality tester is ready --GN */ +#if 0 #define pseudoprime(p) millerrabin(p,3*lgefint(p)) +#else /* remove further overkill :-) --KB */ +#define pseudoprime(p) BSW_psp(p) +#endif /* where to stop trial dividing in factorization */ @@ -408,44 +422,46 @@ GEN auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state), GEN state, long all, long hint) { - long pp[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),2,0 }; - long nb = 0,i,k,lim1,av,lp; + gpmem_t av; + long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 }; + long nb = 0,i,k,lp,p,lim1; byteptr d=diffptr+1; if (typ(n) != t_INT) err(arither1); i=signe(n); if (!i) err(arither2); - cgetg(3,t_MAT); - if (i<0) { stoi(-1); stoi(1); nb++; } + (void)cgetg(3,t_MAT); + if (i<0) { (void)stoi(-1); (void)stoi(1); nb++; } if (is_pm1(n)) return aux_end(NULL,nb); n = gclone(n); setsigne(n,1); i = vali(n); if (i) { - stoi(2); stoi(i); nb++; + (void)stoi(2); (void)stoi(i); nb++; av=avma; affii(shifti(n,-i), n); avma=av; } if (is_pm1(n)) return aux_end(n,nb); lim1 = tridiv_bound(n, all); /* trial division */ - while (*d && pp[2] <= lim1) + p = 2; + while (*d && p <= lim1) { - pp[2] += *d++; - if (mpdivis(n,pp,n)) + NEXT_PRIME_VIADIFF(p,d); + if (mpdivisis(n,p,n)) { - nb++; k=1; while (mpdivis(n,pp,n)) k++; - icopy(pp); stoi(k); + nb++; k=1; while (mpdivisis(n,p,n)) k++; + (void)utoi(p); (void)stoi(k); if (is_pm1(n)) return aux_end(n,nb); } } /* pp = square of biggest p tried so far */ - av=avma; setlg(pp,4); affii(sqri(pp),pp); avma=av; + av=avma; affii(muluu(p,p),pp); avma=av; if (cmpii(pp,n) > 0) { nb++; - icopy(n); stoi(1); return aux_end(n,nb); + (void)icopy(n); (void)stoi(1); return aux_end(n,nb); } /* trial divide by the "special primes" (usually huge composites...) */ @@ -454,7 +470,7 @@ auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, if (mpdivis(n,(GEN) primetab[i],n)) { nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++; - icopy((GEN) primetab[i]); stoi(k); + (void)icopy((GEN) primetab[i]); (void)stoi(k); if (is_pm1(n)) return aux_end(n,nb); } @@ -462,7 +478,7 @@ auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, if ((k && cmpii(pp,n) > 0) || (all==1 && pseudoprime(n))) { nb++; - icopy(n); stoi(1); return aux_end(n,nb); + (void)icopy(n); (void)stoi(1); return aux_end(n,nb); } /* now we have a large composite. Use hint as is if all==1 */ @@ -484,7 +500,7 @@ auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, static long ifac_break_limit(GEN n, GEN pairs/*unused*/, GEN here, GEN state) { - ulong ltop = avma; + gpmem_t ltop = avma; int res; if (!here) /* initial call */ /* Small prime have been removed since start, n is the new unfactored part. @@ -557,7 +573,7 @@ GEN boundfact(GEN n, long lim) { GEN p1,p2,p3,p4,p5,y; - long av = avma,tetpil; + gpmem_t av = avma,tetpil; if (lim<=1) lim=0; switch(typ(n)) @@ -565,7 +581,10 @@ boundfact(GEN n, long lim) case t_INT: return auxdecomp(n,lim); case t_FRACN: - n=gred(n); /* fall through */ + n=gred(n); + if (typ(n) == t_INT) + return gerepileupto(av, boundfact(n,lim)); + /* fall through */ case t_FRAC: p1=auxdecomp((GEN)n[1],lim); p2=auxdecomp((GEN)n[2],lim); @@ -581,7 +600,51 @@ boundfact(GEN n, long lim) err(arither1); return NULL; /* not reached */ } +/***********************************************************************/ +/** **/ +/** SIMPLE FACTORISATIONS **/ +/** **/ +/***********************************************************************/ +/*Factorize n and output [fp,fe] where + * fp and fe are vecsmall and n=prod{fp[i]^fe[i]} + */ + +GEN +decomp_small(long n) +{ + gpmem_t ltop=avma; + GEN F = factor(stoi(n)); + long i, l=lg(F[1]); + GEN f=cgetg(3,t_VEC); + GEN fp=cgetg(l,t_VECSMALL); + GEN fe=cgetg(l,t_VECSMALL); + f[1] = (long) fp; + f[2] = (long) fe; + for(i = 1; i < l; i++) + { + fp[i]=itos(gcoeff(F,i,1)); + fe[i]=itos(gcoeff(F,i,2)); + } + return gerepileupto(ltop,f); +} + +/*Return the primary factors of a small integer as a vecsmall*/ +GEN +decomp_primary_small(long n) +{ + gpmem_t ltop=avma; + GEN F = factor(stoi(n)); + GEN fc = cgetg(lg(F[1]), t_VECSMALL); + gpmem_t av=avma; + long i; + for (i = 1; i < lg(fc); i++) + fc[i] = itos(powgi(gmael(F,1,i), gmael(F,2,i))); + avma=av; + return gerepileupto(ltop,fc); +} + + /***********************************************************************/ /** **/ /** BASIC ARITHMETIC FUNCTIONS **/ @@ -623,8 +686,9 @@ long mu(GEN n) { byteptr d = diffptr+1; /* point at 3 - 2 */ - ulong p,lim1, av = avma; - long s, v; + gpmem_t av = avma; + ulong p; + long s, v, lim1; if (typ(n) != t_INT) err(arither1); if (!signe(n)) err(arither2); @@ -638,7 +702,7 @@ mu(GEN n) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(n,p,n)) { if (smodis(n,p) == 0) { avma=av; return 0; } @@ -664,8 +728,9 @@ gissquarefree(GEN x) long issquarefree(GEN x) { - ulong av = avma,lim1,p; - long tx; + ulong p; + gpmem_t av = avma; + long tx, lim1; GEN d; if (gcmp0(x)) return 0; @@ -683,7 +748,7 @@ issquarefree(GEN x) p = 2; while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(x,p,x)) { if (smodis(x,p) == 0) { avma = av; return 0; } @@ -709,8 +774,8 @@ long omega(GEN n) { byteptr d=diffptr+1; - ulong p, lim1, av = avma; - long nb,v; + gpmem_t av = avma; + long p,nb,v, lim1; if (typ(n) != t_INT) err(arither1); if (!signe(n)) err(arither2); @@ -723,7 +788,7 @@ omega(GEN n) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(n,p,n)) { nb++; while (mpdivisis(n,p,n)); /* empty */ @@ -746,8 +811,9 @@ long bigomega(GEN n) { byteptr d=diffptr+1; - ulong av = avma, p,lim1; - long nb,v; + ulong p; + gpmem_t av = avma; + long nb,v, lim1; if (typ(n) != t_INT) err(arither1); if (!signe(n)) err(arither2); @@ -759,7 +825,7 @@ bigomega(GEN n) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(n,p,n)) { do nb++; while (mpdivisis(n,p,n)); @@ -782,8 +848,9 @@ phi(GEN n) { byteptr d = diffptr+1; GEN m; - ulong av = avma, lim1, p; - long v; + ulong p; + gpmem_t av = avma; + long v, lim1; if (typ(n) != t_INT) err(arither1); if (!signe(n)) err(arither2); @@ -796,7 +863,7 @@ phi(GEN n) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(n,p,n)) { m = mulis(m, p-1); @@ -824,8 +891,9 @@ numbdiv(GEN n) { byteptr d=diffptr+1; GEN m; - long l, v; - ulong av = avma, p,lim1; + long l, v, lim1; + ulong p; + gpmem_t av = avma; if (typ(n) != t_INT) err(arither1); if (!signe(n)) err(arither2); @@ -838,7 +906,7 @@ numbdiv(GEN n) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); l=1; while (mpdivisis(n,p,n)) l++; m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); } } @@ -862,8 +930,9 @@ sumdiv(GEN n) { byteptr d=diffptr+1; GEN m,m1; - ulong av=avma,lim1,p; - long v; + ulong p; + gpmem_t av=avma; + long v, lim1; if (typ(n) != t_INT) err(arither1); if (!signe(n)) err(arither2); @@ -876,7 +945,7 @@ sumdiv(GEN n) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(n,p,n)) { m1 = utoi(p+1); @@ -904,8 +973,9 @@ sumdivk(GEN n, long k) { byteptr d=diffptr+1; GEN n1,m,m1,pk; - ulong av = avma, p, lim1; - long k1,v; + ulong p; + gpmem_t av = avma; + long k1,v, lim1; if (!k) return numbdiv(n); if (k==1) return sumdiv(n); @@ -924,7 +994,7 @@ sumdivk(GEN n, long k) while (*d && p < lim1) { - p += *d++; + NEXT_PRIME_VIADIFF(p,d); if (mpdivisis(n,p,n)) { pk = gpowgs(utoi(p),k); m1 = addsi(1,pk); @@ -954,7 +1024,8 @@ sumdivk(GEN n, long k) GEN divisors(GEN n) { - long tetpil,av=avma,i,j,l; + gpmem_t tetpil,av=avma; + long i,j,l; GEN *d,*t,*t1,*t2,*t3, nbdiv,e; if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1); @@ -978,72 +1049,78 @@ divisors(GEN n) tetpil=avma; return gerepile(av,tetpil,sort((GEN)t)); } -GEN -core(GEN n) +static GEN +_core(GEN n, int all) { - long av=avma,tetpil,i; - GEN fa,p1,p2,res=gun; + gpmem_t av = avma; + long i; + GEN fa,p1,p2,c = gun; - fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2]; + fa = auxdecomp(n,all); + p1 = (GEN)fa[1]; + p2 = (GEN)fa[2]; for (i=1; i0)? cgetg(5,t_QFR): cgetg(4,t_QFI); @@ -1155,7 +1232,8 @@ comp_gen(GEN z,GEN x,GEN y) static GEN compimag0(GEN x, GEN y, int raw) { - ulong tx = typ(x), av = avma; + gpmem_t av = avma; + long tx = typ(x); GEN z; if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition"); @@ -1168,7 +1246,8 @@ compimag0(GEN x, GEN y, int raw) static GEN compreal0(GEN x, GEN y, int raw) { - ulong tx = typ(x), av = avma; + gpmem_t av = avma; + long tx = typ(x); GEN z; if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition"); @@ -1199,7 +1278,7 @@ compraw(GEN x, GEN y) GEN sqcompimag0(GEN x, long raw) { - long av = avma; + gpmem_t av = avma; GEN z = cgetg(4,t_QFI); if (typ(x)!=t_QFI) err(typeer,"composition"); @@ -1211,7 +1290,7 @@ sqcompimag0(GEN x, long raw) static GEN sqcompreal0(GEN x, long raw) { - long av = avma; + gpmem_t av = avma; GEN z = cgetg(5,t_QFR); if (typ(x)!=t_QFR) err(typeer,"composition"); @@ -1236,7 +1315,7 @@ static GEN real_unit_form_by_disc(GEN D, long prec) { GEN y = cgetg(5,t_QFR), isqrtD; - long av = avma; + gpmem_t av = avma; if (typ(D) != t_INT || signe(D) <= 0) err(typeer,"real_unit_form_by_disc"); switch(mod4(D)) @@ -1249,14 +1328,15 @@ real_unit_form_by_disc(GEN D, long prec) if (mod2(D) != mod2(isqrtD)) isqrtD = gerepileuptoint(av, addsi(-1,isqrtD)); y[2] = (long)isqrtD; av = avma; - y[3] = (long)gerepileuptoint(av, shifti(subii(sqri(isqrtD),D),-2)); + y[3] = lpileuptoint(av, shifti(subii(sqri(isqrtD),D),-2)); y[4] = (long)realzero(prec); return y; } GEN real_unit_form(GEN x) { - long av = avma, prec; + gpmem_t av = avma; + long prec; GEN D; if (typ(x) != t_QFR) err(typeer,"real_unit_form"); prec = precision((GEN)x[4]); @@ -1283,8 +1363,8 @@ imag_unit_form_by_disc(GEN D) y[3] = lshifti(D,-2); setsigne(y[3],1); if (isodd) { - long av = avma; - y[3] = (long)gerepileuptoint(av, addis((GEN)y[3],1)); + gpmem_t av = avma; + y[3] = lpileuptoint(av, addis((GEN)y[3],1)); } return y; } @@ -1293,20 +1373,21 @@ GEN imag_unit_form(GEN x) { GEN p1,p2, y = cgetg(4,t_QFI); - long av; + gpmem_t av; if (typ(x) != t_QFI) err(typeer,"imag_unit_form"); y[1] = un; y[2] = mpodd((GEN)x[2])? un: zero; av = avma; p1 = mulii((GEN)x[1],(GEN)x[3]); p2 = shifti(sqri((GEN)x[2]),-2); - y[3] = (long)gerepileuptoint(av, subii(p1,p2)); + y[3] = lpileuptoint(av, subii(p1,p2)); return y; } GEN powrealraw(GEN x, long n) { - long av = avma, m; + gpmem_t av = avma; + long m; GEN y; if (typ(x) != t_QFR) @@ -1329,7 +1410,8 @@ powrealraw(GEN x, long n) GEN powimagraw(GEN x, long n) { - long av = avma, m; + gpmem_t av = avma; + long m; GEN y; if (typ(x) != t_QFI) @@ -1360,7 +1442,8 @@ powraw(GEN x, long n) GEN nucomp(GEN x, GEN y, GEN l) { - long av=avma,tetpil,cz; + gpmem_t av=avma,tetpil; + long cz; 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; if (x==y) return nudupl(x,l); @@ -1424,7 +1507,8 @@ nucomp(GEN x, GEN y, GEN l) GEN nudupl(GEN x, GEN l) { - long av=avma,tetpil,cz; + gpmem_t av=avma,tetpil; + long cz; GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g; if (typ(x) != t_QFI) @@ -1464,7 +1548,8 @@ nudupl(GEN x, GEN l) GEN nupow(GEN x, GEN n) { - long av,tetpil,i,j; + gpmem_t av,tetpil; + long i,j; unsigned long m; GEN y,l; @@ -1551,14 +1636,13 @@ rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD) GEN p1,p2, y = cgetg(6,t_VEC); GEN b = (GEN)x[2]; GEN c = (GEN)x[3]; - long s = signe(c); y[1] = (long)c; - p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: c; - setsigne(c,1); + p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: absi(c); p1 = shifti(c,1); + if (p1 == gzero) err(talker, "reducible form in rhoreal"); + setsigne(p1,1); /* |2c| */ p2 = divii(addii(p2,b), p1); - setsigne(c,s); y[2] = lsubii(mulii(p2,p1), b); p1 = shifti(subii(sqri((GEN)y[2]),D),-2); @@ -1586,11 +1670,11 @@ GEN codeform5(GEN x, long prec) { GEN y = cgetg(6,t_VEC); - y[1]=x[1]; - y[2]=x[2]; - y[3]=x[3]; - y[4]=zero; - y[5]=(long)realun(prec); return y; + y[1] = x[1]; + y[2] = x[2]; + y[3] = x[3]; + y[4] = zero; + y[5] = (long)realun(prec); return y; } static GEN @@ -1619,7 +1703,7 @@ decodeform(GEN x, GEN d0) p1 = shiftr(p1,-e); p2 = addis(mulsi(EXP220,p2), e); p1 = mplog(p1); - p1 = mpadd(p1, mulir(p2, glog(gdeux, lg(d0)))); + p1 = mpadd(p1, mulir(p2, mplog2(lg(d0)))); } else { /* to avoid loss of precision */ @@ -1665,7 +1749,8 @@ redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD) static GEN redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD) { - long av = avma, prec; + gpmem_t av = avma; + long prec; GEN d0; if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal"); @@ -1703,7 +1788,7 @@ redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrt GEN comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD) { - ulong av = avma; + gpmem_t av = avma; GEN z = cgetg(6,t_VEC); comp_gen(z,x,y); if (x == y) { @@ -1723,7 +1808,8 @@ comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqr GEN powrealform(GEN x, GEN n) { - long av = avma, i,m; + gpmem_t av = avma; + long i,m; GEN y,D,sqrtD,isqrtD,d0; if (typ(x) != t_QFR) @@ -1755,7 +1841,8 @@ powrealform(GEN x, GEN n) GEN redimag(GEN x) { - long av=avma, fl; + gpmem_t av=avma; + long fl; do x = rhoimag0(x, &fl); while (fl == 0); x = gerepilecopy(av,x); if (fl == 2) setsigne(x[2], -signe(x[2])); @@ -1790,7 +1877,7 @@ GEN qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD) { long tx=typ(x),fl; - ulong av; + gpmem_t av; if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred"); if (!D) D = qf_disc(x,NULL,NULL); @@ -1815,51 +1902,55 @@ qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD GEN primeform(GEN x, GEN p, long prec) { - long av,tetpil,s=signe(x); - GEN y,b; + gpmem_t av; + long s, sx = signe(x); + GEN y, b, c; - if (typ(x) != t_INT || !s) err(arither1); + if (typ(x) != t_INT || !sx) err(arither1); if (typ(p) != t_INT || signe(p) <= 0) err(arither1); if (is_pm1(p)) - return s<0? imag_unit_form_by_disc(x) - : real_unit_form_by_disc(x,prec); - if (s < 0) + return sx<0? imag_unit_form_by_disc(x) + : real_unit_form_by_disc(x,prec); + s = mod8(x); + if (sx < 0) { - y = cgetg(4,t_QFI); s = 8-mod8(x); + if (s) s = 8-s; + y = cgetg(4, t_QFI); } else { - y = cgetg(5,t_QFR); s = mod8(x); - y[4]=(long)realzero(prec); + y = cgetg(5, t_QFR); + y[4] = (long)realzero(prec); } switch(s&3) { case 2: case 3: err(funder2,"primeform"); } - y[1] = (long)icopy(p); av = avma; - if (egalii(p,gdeux)) + av = avma; + if (egalii(p, gdeux)) { switch(s) { - case 0: y[2]=zero; break; - case 8: s=0; y[2]=zero; break; - case 1: s=1; y[2]=un; break; - case 4: s=4; y[2]=deux; break; - default: err(sqrter5); + case 0: b = gzero; break; + case 1: b = gun; break; + case 4: b = gdeux; break; + default: err(sqrter5); b = NULL; /* -Wall */ } - b=subsi(s,x); tetpil=avma; b=shifti(b,-3); + c = shifti(subsi(s,x), -3); } else { b = mpsqrtmod(x,p); if (!b) err(sqrter5); - if (mod2(b) == mod2(x)) y[2] = (long)b; - else - { tetpil = avma; y[2] = lpile(av,tetpil,subii(p,b)); } + s &= 1; /* s = x mod 2 */ + /* mod(b) != mod2(x) ? [Warning: we may have b == 0] */ + if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(p,b)); - av=avma; b=shifti(subii(sqri((GEN)y[2]),x),-2); - tetpil=avma; b=divii(b,p); + av = avma; + c = diviiexact(shifti(subii(sqri(b), x), -2), p); } - y[3]=lpile(av,tetpil,b); return y; + y[3] = lpileuptoint(av, c); + y[2] = (long)b; + y[1] = licopy(p); return y; } /*********************************************************************/ @@ -1955,6 +2046,100 @@ bittest(GEN x, long n) return n? 1: 0; } +GEN +bittest_many(GEN x, GEN gn, long c) +{ + long lx = lgefint(x), l1, l2, b1, b2, two_adic = 0, l_res, skip; + ulong *p, *pnew; + long n = itos(gn), extra_words = 0, partial_bits; + GEN res; + + if (c == 1) /* Shortcut */ + return bittest(x,n) ? gun : gzero; + /* Negative n with n+c>0 are too hairy to implement... */ + if (!signe(x) || c == 0) + return gzero; + if (c < 0) { /* Negative x means 2s complemenent */ + c = -c; + if (signe(x) < 0) + two_adic = 1; /* treat x 2-adically... */ + } + if (n < 0) { + gpmem_t oa, na; + + if (n + c <= 0) + return gzero; + oa = avma; + res = bittest_many(x, gzero, two_adic? -(n+c) : n+c); + na = avma; + res = shifti3(res,-n,0); + return gerepile(oa,na,res); + } + partial_bits = (c & (BITS_IN_LONG-1)); + l2 = lx-1 - (n>>TWOPOTBITS_IN_LONG); /* Last=least-significant word */ + l1 = lx-1 - ((n + c - 1)>>TWOPOTBITS_IN_LONG); /* First word */ + b2 = (n & (BITS_IN_LONG-1)); /* Last bit, skip less-signifant */ + b1 = ((n + c - 1) & (BITS_IN_LONG-1)); /* First bit, skip more-significant */ + if (l2 < 2) { /* always: l1 <= l2 */ + if (!two_adic) + return gzero; + /* x is non-zero, so it behaves as -1. */ + return gbitneg(gzero,c); + } + if (l1 < 2) { + if (two_adic) /* If b1 < b2, bits are set by prepend in shift_r */ + extra_words = 2 - l1 - (b1 < b2); + else + partial_bits = b2 ? BITS_IN_LONG - b2 : 0; + l1 = 2; + b1 = (BITS_IN_LONG-1); /* Include all bits in this word */ + } + skip = (b1 < b2); /* Skip the first word of the shift */ + l_res = l2 - l1 + 1 + 2 - skip + extra_words; /* A coarse estimate */ + p = (ulong*) (new_chunk(l_res) + 2 + extra_words); + shift_r(p - skip, (ulong*)x + l1, (ulong*)x + l2 + 1, 0, b2); + if (two_adic) { /* Check the low bits of x */ + int i = lx, nonzero = 0; + + /* Flip the bits */ + pnew = p + l_res - 2 - extra_words; + while (--pnew >= p) + *pnew = MAXULONG - *pnew; + /* Fill the extra words */ + while (extra_words--) + *--p = MAXULONG; + /* The result is one less than 2s-complement if the lower-bits + of x were all 0. */ + while (--i > l2) { + if (x[i]) { + nonzero = 1; + break; + } + } + if (!nonzero && b2) /* Check the partial word too */ + nonzero = x[l2] & ((1<= p) + if (++*pnew) + break; + } + } + if (partial_bits) + *p &= (1<= p - 2 + l_res) + return gzero; + l_res -= (pnew - p); + res = (GEN)(pnew - 2); + res[1] = evalsigne(1)|evallgefint(l_res); + res[0] = evaltyp(t_INT)|evallg(l_res); + return res; +} + static long bittestg(GEN x, GEN n) { @@ -1967,6 +2152,12 @@ gbittest(GEN x, GEN n) return arith_proto2(bittestg,x,n); } +GEN +gbittest3(GEN x, GEN n, long c) +{ + return garith_proto3ggs(bittest_many,x,n,c); +} + /***********************************************************************/ /** **/ /** BITMAP OPS **/ @@ -2226,9 +2417,9 @@ ibitnegimply(GEN x, GEN y) } static GEN -inegate_inplace(GEN z, long ltop) +inegate_inplace(GEN z, gpmem_t ltop) { - long lbot, o; + long o; /* Negate z */ o = incdec(z, 1); /* Can overflow z... */ @@ -2238,16 +2429,14 @@ inegate_inplace(GEN z, long ltop) else if (lgefint(z) == 2) setsigne(z, 0); incdec(z,-1); /* Restore a normalized value */ - lbot = avma; - z = gsub(z,gun); - return gerepile(ltop,lbot,z); + return gerepileupto(ltop, subis(z,1)); } GEN gbitor(GEN x, GEN y) { + gpmem_t ltop; long sx,sy; - long ltop; GEN z; if (typ(x) != t_INT || typ(y) != t_INT) @@ -2278,8 +2467,8 @@ gbitor(GEN x, GEN y) GEN gbitand(GEN x, GEN y) { + gpmem_t ltop; long sx,sy; - long ltop; GEN z; if (typ(x) != t_INT || typ(y) != t_INT) @@ -2312,8 +2501,8 @@ gbitand(GEN x, GEN y) GEN gbitxor(GEN x, GEN y) { + gpmem_t ltop; long sx,sy; - long ltop; GEN z; if (typ(x) != t_INT || typ(y) != t_INT) @@ -2345,8 +2534,8 @@ gbitxor(GEN x, GEN y) GEN gbitnegimply(GEN x, GEN y) /* x & ~y */ { + gpmem_t ltop; long sx,sy; - long ltop; GEN z; if (typ(x) != t_INT || typ(y) != t_INT)