version 1.1, 2001/10/02 11:17:01 |
version 1.2, 2002/09/11 07:26:48 |
Line 25 extern GEN arith_proto(long f(GEN), GEN x, int do_erro |
|
Line 25 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 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_proto(GEN f(GEN), GEN x, int do_error); |
extern GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y); |
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))) |
#define sqru(i) (muluu((i),(i))) |
|
|
|
|
prime(long n) |
prime(long n) |
{ |
{ |
byteptr p = diffptr; |
byteptr p = diffptr; |
long c, prime = 0; |
long prime = 0; |
|
|
if (n <= 0) err(talker, "n-th prime meaningless if n = %ld",n); |
if (n <= 0) err(talker, "n-th prime meaningless if n = %ld",n); |
while (n--) |
while (n--) |
{ |
{ |
c = *p++; if (!c) err(primer1); |
NEXT_PRIME_VIADIFF_CHECK(prime,p); |
prime += c; |
|
} |
} |
return stoi(prime); |
return stoi(prime); |
} |
} |
|
|
ulong prime = 0, res = 0; |
ulong prime = 0, res = 0; |
|
|
if (n <= 0) err(talker, "pith meaningless if n = %ld",n); |
if (n <= 0) err(talker, "pith meaningless if n = %ld",n); |
if (maxprime() <= n) err(primer1); |
if (maxprime() <= (ulong)n) err(primer1); |
while (prime<=n) { prime += *p++; res++; } |
while (prime<=(ulong)n) { |
|
NEXT_PRIME_VIADIFF(prime,p); |
|
res++; |
|
} |
return stoi(res-1); |
return stoi(res-1); |
} |
} |
|
|
|
|
primes(long n) |
primes(long n) |
{ |
{ |
byteptr p = diffptr; |
byteptr p = diffptr; |
long c, prime = 0; |
long prime = 0; |
GEN y,z; |
GEN y,z; |
|
|
if (n < 0) n = 0; |
if (n < 0) n = 0; |
z = y = cgetg(n+1,t_VEC); |
z = y = cgetg(n+1,t_VEC); |
while (n--) |
while (n--) |
{ |
{ |
c = *p++; if (!c) err(primer1); |
NEXT_PRIME_VIADIFF_CHECK(prime,p); |
prime += c; *++z = lstoi(prime); |
*++z = lstoi(prime); |
} |
} |
return y; |
return y; |
} |
} |
|
|
* when maxnum (size) is moderate. |
* when maxnum (size) is moderate. |
*/ |
*/ |
static byteptr |
static byteptr |
initprimes1(long size, long *lenp, long *lastp) |
initprimes1(ulong size, long *lenp, long *lastp) |
{ |
{ |
long k; |
long k; |
byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2); |
byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2); |
Line 110 initprimes1(long size, long *lenp, long *lastp) |
|
Line 113 initprimes1(long size, long *lenp, long *lastp) |
|
fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n", |
fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n", |
q, *lenp, *lastp); |
q, *lenp, *lastp); |
#endif |
#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 */ |
#define PRIME_ARENA (200 * 1024) /* No slowdown even with 256K level-2 cache */ |
Line 119 initprimes1(long size, long *lenp, long *lastp) |
|
Line 122 initprimes1(long size, long *lenp, long *lastp) |
|
recursive call will bottom out and invoke initprimes1() at once. |
recursive call will bottom out and invoke initprimes1() at once. |
(Not static; might conceivably be useful to someone in library mode) */ |
(Not static; might conceivably be useful to someone in library mode) */ |
byteptr |
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; |
byteptr q,r,s,fin, p, p1, fin1, plast, curdiff; |
|
ulong last, remains, curlow; |
|
|
#if 0 |
#if 0 |
fprintferr("initprimes0: called for maxnum = %lu\n", maxnum); |
fprintferr("initprimes0: called for maxnum = %lu\n", maxnum); |
Line 160 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
Line 164 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
asize = ARENA_IN_ROOTS * rootnum; /* Make % overhead negligeable. */ |
asize = ARENA_IN_ROOTS * rootnum; /* Make % overhead negligeable. */ |
if (asize < PRIME_ARENA) asize = PRIME_ARENA; |
if (asize < PRIME_ARENA) asize = PRIME_ARENA; |
if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */ |
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) |
if (alloced) |
p = (byteptr) gpmalloc(asize); |
p = (byteptr) gpmalloc(asize); |
else |
else |
Line 182 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
Line 186 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
} |
} |
memset(p, 0, asize); |
memset(p, 0, asize); |
/* p corresponds to curlow. q runs over primediffs. */ |
/* 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++) |
for (q = p1+2, k = 3; q <= fin1; k += *q++) |
{ |
{ |
/* The first odd number which is >= curlow and divisible by p |
/* The first odd number which is >= curlow and divisible by p |
Line 207 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
Line 212 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
/* now q runs over addresses corresponding to primes */ |
/* now q runs over addresses corresponding to primes */ |
for (q = p; ; plast = q++) |
for (q = p; ; plast = q++) |
{ |
{ |
|
long d; |
|
|
while (*q) q++; /* use 0-sentinel at end */ |
while (*q) q++; /* use 0-sentinel at end */ |
if (q >= fin) break; |
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; |
plast -= asize - 1; |
remains -= asize - 1; |
remains -= asize - 1; |
Line 220 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
Line 230 initprimes0(ulong maxnum, long *lenp, long *lastp) |
|
*lenp = curdiff - p1; |
*lenp = curdiff - p1; |
*lastp = last; |
*lastp = last; |
if (alloced) free(p); |
if (alloced) free(p); |
return (byteptr) gprealloc(p1, *lenp, size); |
return (byteptr) gprealloc(p1, *lenp); |
} |
} |
#if 0 /* not yet... GN */ |
#if 0 /* not yet... GN */ |
/* The diffptr table will contain at least 6548 entries (up to and including |
/* The diffptr table will contain at least 6548 entries (up to and including |
Line 235 init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */ |
|
Line 245 init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */ |
|
static ulong _maxprime = 0; |
static ulong _maxprime = 0; |
|
|
ulong |
ulong |
maxprime() { return _maxprime; } |
maxprime(void) { return _maxprime; } |
|
|
byteptr |
byteptr |
initprimes(long maxnum) |
initprimes(ulong maxnum) |
{ |
{ |
long len, last; |
long len; |
|
ulong last; |
byteptr p; |
byteptr p; |
/* The algorithm must see the next prime beyond maxnum, whence the +257. */ |
/* The algorithm must see the next prime beyond maxnum, whence the +512. */ |
/* switch to unsigned: adding the 257 _could_ wrap to a negative number. */ |
/* switch to unsigned: adding the 512 _could_ wrap to a negative number. */ |
ulong maxnum1 = ((maxnum<65302?65302:maxnum)+257ul); |
ulong maxnum1 = ((maxnum<65302?65302:maxnum)+512ul); |
|
|
if (maxnum1 > 436273000) |
if ((maxnum>>1) > VERYBIGINT - 1024) |
err(talker,"impossible to have prestored primes > 436273009"); |
err(talker, "Too large primelimit"); |
|
|
p = initprimes0(maxnum1, &len, &last); |
p = initprimes0(maxnum1, &len, &last); |
#if 0 /* not yet... GN */ |
#if 0 /* not yet... GN */ |
static int build_the_tables = 1; |
static int build_the_tables = 1; |
Line 273 cleanprimetab(void) |
|
Line 283 cleanprimetab(void) |
|
GEN |
GEN |
addprimes(GEN p) |
addprimes(GEN p) |
{ |
{ |
ulong av; |
gpmem_t av; |
long i,k,tx,lp; |
long i,k,tx,lp; |
GEN L; |
GEN L; |
|
|
Line 281 addprimes(GEN p) |
|
Line 291 addprimes(GEN p) |
|
tx = typ(p); |
tx = typ(p); |
if (is_vec_t(tx)) |
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; |
return primetab; |
} |
} |
if (tx != t_INT) err(typeer,"addprime"); |
if (tx != t_INT) err(typeer,"addprime"); |
Line 302 addprimes(GEN p) |
|
Line 312 addprimes(GEN p) |
|
gunclone(n); primetab[i] = 0; |
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); |
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; |
avma = av; return primetab; |
} |
} |
|
|
Line 340 removeprimes(GEN prime) |
|
Line 350 removeprimes(GEN prime) |
|
} |
} |
else |
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; |
return primetab; |
} |
} |
Line 361 removeprimes(GEN prime) |
|
Line 371 removeprimes(GEN prime) |
|
/* Some overkill removed from this (15 spsp for an integer < 2^32 ?!). |
/* Some overkill removed from this (15 spsp for an integer < 2^32 ?!). |
Should really revert to isprime() once the new primality tester is ready |
Should really revert to isprime() once the new primality tester is ready |
--GN */ |
--GN */ |
|
#if 0 |
#define pseudoprime(p) millerrabin(p,3*lgefint(p)) |
#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 */ |
/* where to stop trial dividing in factorization */ |
|
|
|
|
auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state), |
auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, GEN here, GEN state), |
GEN state, long all, long hint) |
GEN state, long all, long hint) |
{ |
{ |
long pp[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),2,0 }; |
gpmem_t av; |
long nb = 0,i,k,lim1,av,lp; |
long pp[] = { evaltyp(t_INT)|_evallg(4), 0,0,0 }; |
|
long nb = 0,i,k,lp,p,lim1; |
byteptr d=diffptr+1; |
byteptr d=diffptr+1; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
i=signe(n); if (!i) err(arither2); |
i=signe(n); if (!i) err(arither2); |
cgetg(3,t_MAT); |
(void)cgetg(3,t_MAT); |
if (i<0) { stoi(-1); stoi(1); nb++; } |
if (i<0) { (void)stoi(-1); (void)stoi(1); nb++; } |
if (is_pm1(n)) return aux_end(NULL,nb); |
if (is_pm1(n)) return aux_end(NULL,nb); |
|
|
n = gclone(n); setsigne(n,1); |
n = gclone(n); setsigne(n,1); |
i = vali(n); |
i = vali(n); |
if (i) |
if (i) |
{ |
{ |
stoi(2); stoi(i); nb++; |
(void)stoi(2); (void)stoi(i); nb++; |
av=avma; affii(shifti(n,-i), n); avma=av; |
av=avma; affii(shifti(n,-i), n); avma=av; |
} |
} |
if (is_pm1(n)) return aux_end(n,nb); |
if (is_pm1(n)) return aux_end(n,nb); |
lim1 = tridiv_bound(n, all); |
lim1 = tridiv_bound(n, all); |
|
|
/* trial division */ |
/* trial division */ |
while (*d && pp[2] <= lim1) |
p = 2; |
|
while (*d && p <= lim1) |
{ |
{ |
pp[2] += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivis(n,pp,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
nb++; k=1; while (mpdivis(n,pp,n)) k++; |
nb++; k=1; while (mpdivisis(n,p,n)) k++; |
icopy(pp); stoi(k); |
(void)utoi(p); (void)stoi(k); |
if (is_pm1(n)) return aux_end(n,nb); |
if (is_pm1(n)) return aux_end(n,nb); |
} |
} |
} |
} |
|
|
/* pp = square of biggest p tried so far */ |
/* 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) |
if (cmpii(pp,n) > 0) |
{ |
{ |
nb++; |
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...) */ |
/* trial divide by the "special primes" (usually huge composites...) */ |
Line 454 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, |
|
Line 470 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, |
|
if (mpdivis(n,(GEN) primetab[i],n)) |
if (mpdivis(n,(GEN) primetab[i],n)) |
{ |
{ |
nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++; |
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); |
if (is_pm1(n)) return aux_end(n,nb); |
} |
} |
|
|
Line 462 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, |
|
Line 478 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, |
|
if ((k && cmpii(pp,n) > 0) || (all==1 && pseudoprime(n))) |
if ((k && cmpii(pp,n) > 0) || (all==1 && pseudoprime(n))) |
{ |
{ |
nb++; |
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 */ |
/* now we have a large composite. Use hint as is if all==1 */ |
Line 484 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, |
|
Line 500 auxdecomp1(GEN n, long (*ifac_break)(GEN n, GEN pairs, |
|
static long |
static long |
ifac_break_limit(GEN n, GEN pairs/*unused*/, GEN here, GEN state) |
ifac_break_limit(GEN n, GEN pairs/*unused*/, GEN here, GEN state) |
{ |
{ |
ulong ltop = avma; |
gpmem_t ltop = avma; |
int res; |
int res; |
if (!here) /* initial call */ |
if (!here) /* initial call */ |
/* Small prime have been removed since start, n is the new unfactored part. |
/* Small prime have been removed since start, n is the new unfactored part. |
|
|
boundfact(GEN n, long lim) |
boundfact(GEN n, long lim) |
{ |
{ |
GEN p1,p2,p3,p4,p5,y; |
GEN p1,p2,p3,p4,p5,y; |
long av = avma,tetpil; |
gpmem_t av = avma,tetpil; |
|
|
if (lim<=1) lim=0; |
if (lim<=1) lim=0; |
switch(typ(n)) |
switch(typ(n)) |
Line 565 boundfact(GEN n, long lim) |
|
Line 581 boundfact(GEN n, long lim) |
|
case t_INT: |
case t_INT: |
return auxdecomp(n,lim); |
return auxdecomp(n,lim); |
case t_FRACN: |
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: |
case t_FRAC: |
p1=auxdecomp((GEN)n[1],lim); |
p1=auxdecomp((GEN)n[1],lim); |
p2=auxdecomp((GEN)n[2],lim); |
p2=auxdecomp((GEN)n[2],lim); |
Line 581 boundfact(GEN n, long lim) |
|
Line 600 boundfact(GEN n, long lim) |
|
err(arither1); |
err(arither1); |
return NULL; /* not reached */ |
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 **/ |
/** BASIC ARITHMETIC FUNCTIONS **/ |
|
|
mu(GEN n) |
mu(GEN n) |
{ |
{ |
byteptr d = diffptr+1; /* point at 3 - 2 */ |
byteptr d = diffptr+1; /* point at 3 - 2 */ |
ulong p,lim1, av = avma; |
gpmem_t av = avma; |
long s, v; |
ulong p; |
|
long s, v, lim1; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) err(arither2); |
if (!signe(n)) err(arither2); |
|
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(n,p,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
if (smodis(n,p) == 0) { avma=av; return 0; } |
if (smodis(n,p) == 0) { avma=av; return 0; } |
Line 664 gissquarefree(GEN x) |
|
Line 728 gissquarefree(GEN x) |
|
long |
long |
issquarefree(GEN x) |
issquarefree(GEN x) |
{ |
{ |
ulong av = avma,lim1,p; |
ulong p; |
long tx; |
gpmem_t av = avma; |
|
long tx, lim1; |
GEN d; |
GEN d; |
|
|
if (gcmp0(x)) return 0; |
if (gcmp0(x)) return 0; |
Line 683 issquarefree(GEN x) |
|
Line 748 issquarefree(GEN x) |
|
p = 2; |
p = 2; |
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(x,p,x)) |
if (mpdivisis(x,p,x)) |
{ |
{ |
if (smodis(x,p) == 0) { avma = av; return 0; } |
if (smodis(x,p) == 0) { avma = av; return 0; } |
|
|
omega(GEN n) |
omega(GEN n) |
{ |
{ |
byteptr d=diffptr+1; |
byteptr d=diffptr+1; |
ulong p, lim1, av = avma; |
gpmem_t av = avma; |
long nb,v; |
long p,nb,v, lim1; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) err(arither2); |
if (!signe(n)) err(arither2); |
|
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(n,p,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
nb++; while (mpdivisis(n,p,n)); /* empty */ |
nb++; while (mpdivisis(n,p,n)); /* empty */ |
|
|
bigomega(GEN n) |
bigomega(GEN n) |
{ |
{ |
byteptr d=diffptr+1; |
byteptr d=diffptr+1; |
ulong av = avma, p,lim1; |
ulong p; |
long nb,v; |
gpmem_t av = avma; |
|
long nb,v, lim1; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) err(arither2); |
if (!signe(n)) err(arither2); |
|
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(n,p,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
do nb++; while (mpdivisis(n,p,n)); |
do nb++; while (mpdivisis(n,p,n)); |
|
|
{ |
{ |
byteptr d = diffptr+1; |
byteptr d = diffptr+1; |
GEN m; |
GEN m; |
ulong av = avma, lim1, p; |
ulong p; |
long v; |
gpmem_t av = avma; |
|
long v, lim1; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) err(arither2); |
if (!signe(n)) err(arither2); |
|
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(n,p,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
m = mulis(m, p-1); |
m = mulis(m, p-1); |
|
|
{ |
{ |
byteptr d=diffptr+1; |
byteptr d=diffptr+1; |
GEN m; |
GEN m; |
long l, v; |
long l, v, lim1; |
ulong av = avma, p,lim1; |
ulong p; |
|
gpmem_t av = avma; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) err(arither2); |
if (!signe(n)) err(arither2); |
|
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
l=1; while (mpdivisis(n,p,n)) l++; |
l=1; while (mpdivisis(n,p,n)) l++; |
m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); } |
m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); } |
} |
} |
|
|
{ |
{ |
byteptr d=diffptr+1; |
byteptr d=diffptr+1; |
GEN m,m1; |
GEN m,m1; |
ulong av=avma,lim1,p; |
ulong p; |
long v; |
gpmem_t av=avma; |
|
long v, lim1; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) err(arither2); |
if (!signe(n)) err(arither2); |
|
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(n,p,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
m1 = utoi(p+1); |
m1 = utoi(p+1); |
Line 904 sumdivk(GEN n, long k) |
|
Line 973 sumdivk(GEN n, long k) |
|
{ |
{ |
byteptr d=diffptr+1; |
byteptr d=diffptr+1; |
GEN n1,m,m1,pk; |
GEN n1,m,m1,pk; |
ulong av = avma, p, lim1; |
ulong p; |
long k1,v; |
gpmem_t av = avma; |
|
long k1,v, lim1; |
|
|
if (!k) return numbdiv(n); |
if (!k) return numbdiv(n); |
if (k==1) return sumdiv(n); |
if (k==1) return sumdiv(n); |
Line 924 sumdivk(GEN n, long k) |
|
Line 994 sumdivk(GEN n, long k) |
|
|
|
while (*d && p < lim1) |
while (*d && p < lim1) |
{ |
{ |
p += *d++; |
NEXT_PRIME_VIADIFF(p,d); |
if (mpdivisis(n,p,n)) |
if (mpdivisis(n,p,n)) |
{ |
{ |
pk = gpowgs(utoi(p),k); m1 = addsi(1,pk); |
pk = gpowgs(utoi(p),k); m1 = addsi(1,pk); |
Line 954 sumdivk(GEN n, long k) |
|
Line 1024 sumdivk(GEN n, long k) |
|
GEN |
GEN |
divisors(GEN n) |
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; |
GEN *d,*t,*t1,*t2,*t3, nbdiv,e; |
|
|
if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1); |
if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1); |
|
Line 1049 divisors(GEN n) |
|
tetpil=avma; return gerepile(av,tetpil,sort((GEN)t)); |
tetpil=avma; return gerepile(av,tetpil,sort((GEN)t)); |
} |
} |
|
|
GEN |
static GEN |
core(GEN n) |
_core(GEN n, int all) |
{ |
{ |
long av=avma,tetpil,i; |
gpmem_t av = avma; |
GEN fa,p1,p2,res=gun; |
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; i<lg(p1); i++) |
for (i=1; i<lg(p1); i++) |
if (mod2((GEN)p2[i])) res = mulii(res,(GEN)p1[i]); |
if (mod2((GEN)p2[i])) c = mulii(c,(GEN)p1[i]); |
tetpil=avma; return gerepile(av,tetpil,icopy(res)); |
return gerepileuptoint(av, c); |
} |
} |
|
|
GEN |
static GEN |
core2(GEN n) |
_core2(GEN n, int all) |
{ |
{ |
long av=avma,tetpil,i; |
gpmem_t av = avma; |
|
long i; |
|
GEN fa,p1,p2,e,c=gun,f=gun,y; |
|
|
GEN fa,p1,p2,p3,res=gun,res2=gun,y; |
fa = auxdecomp(n,all); |
|
p1 = (GEN)fa[1]; |
fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2]; |
p2 = (GEN)fa[2]; |
for (i=1; i<lg(p1); i++) |
for (i=1; i<lg(p1); i++) |
{ |
{ |
p3=(GEN)p2[i]; |
e = (GEN)p2[i]; |
if (mod2(p3)) res=mulii(res,(GEN)p1[i]); |
if (mod2(e)) c = mulii(c, (GEN)p1[i]); |
if (!gcmp1(p3)) res2=mulii(res2,gpui((GEN)p1[i],shifti(p3,-1),0)); |
if (!gcmp1(e)) f = mulii(f, powgi((GEN)p1[i], shifti(e,-1))); |
} |
} |
tetpil=avma; y=cgetg(3,t_VEC); |
y = cgetg(3,t_VEC); |
y[1]=(long)icopy(res); y[2]=(long)icopy(res2); |
y[1] = (long)c; |
return gerepile(av,tetpil,y); |
y[2] = (long)f; return gerepilecopy(av, y); |
} |
} |
|
|
|
GEN core(GEN n) { return _core(n,1); } |
|
GEN corepartial(GEN n) { return _core(n,0); } |
|
GEN core2(GEN n) { return _core2(n,1); } |
|
GEN core2partial(GEN n){ return _core2(n,0); } |
|
|
GEN |
GEN |
core0(GEN n,long flag) |
core0(GEN n,long flag) { return flag? core2(n): core(n); } |
{ |
|
return flag? core2(n): core(n); |
|
} |
|
|
|
|
static long |
|
_mod4(GEN c) { long r = mod4(c); if (signe(c) < 0) r = 4-r; return r; } |
|
|
GEN |
GEN |
coredisc(GEN n) |
coredisc(GEN n) |
{ |
{ |
long av=avma,tetpil,r; |
gpmem_t av = avma; |
GEN p1=core(n); |
GEN c = core(n); |
|
long r = _mod4(c); |
r=mod4(p1); if (signe(p1)<0) r = 4-r; |
if (r==1 || r==4) return c; |
if (r==1 || r==4) return p1; |
return gerepileuptoint(av, shifti(c,2)); |
tetpil=avma; return gerepile(av,tetpil,shifti(p1,2)); |
|
} |
} |
|
|
GEN |
GEN |
coredisc2(GEN n) |
coredisc2(GEN n) |
{ |
{ |
long av=avma,tetpil,r; |
gpmem_t av = avma; |
GEN y,p1,p2=core2(n); |
GEN y = core2(n); |
|
GEN c = (GEN)y[1], f = (GEN)y[2]; |
p1=(GEN)p2[1]; r=mod4(p1); if (signe(p1)<0) r=4-r; |
long r = _mod4(c); |
if (r==1 || r==4) return p2; |
if (r==1 || r==4) return y; |
tetpil=avma; y=cgetg(3,t_VEC); |
y = cgetg(3,t_VEC); |
y[1]=lshifti(p1,2); y[2]=lmul2n((GEN)p2[2],-1); |
y[1] = lshifti(c,2); |
return gerepile(av,tetpil,y); |
y[2] = lmul2n(f,-1); return gerepileupto(av, y); |
} |
} |
|
|
GEN |
GEN |
coredisc0(GEN n,long flag) |
coredisc0(GEN n,long flag) { return flag? coredisc2(n): coredisc(n); } |
{ |
|
return flag? coredisc2(n): coredisc(n); |
|
} |
|
|
|
/***********************************************************************/ |
/***********************************************************************/ |
/** **/ |
/** **/ |
Line 1065 qf_create(GEN x, GEN y, GEN z, long s) |
|
Line 1142 qf_create(GEN x, GEN y, GEN z, long s) |
|
if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb"); |
if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb"); |
if (!s) |
if (!s) |
{ |
{ |
long av=avma; s = signe(qf_disc(x,y,z)); avma=av; |
gpmem_t av=avma; s = signe(qf_disc(x,y,z)); avma=av; |
if (!s) err(talker,"zero discriminant in Qfb"); |
if (!s) err(talker,"zero discriminant in Qfb"); |
} |
} |
t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI); |
t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI); |
Line 1155 comp_gen(GEN z,GEN x,GEN y) |
|
Line 1232 comp_gen(GEN z,GEN x,GEN y) |
|
static GEN |
static GEN |
compimag0(GEN x, GEN y, int raw) |
compimag0(GEN x, GEN y, int raw) |
{ |
{ |
ulong tx = typ(x), av = avma; |
gpmem_t av = avma; |
|
long tx = typ(x); |
GEN z; |
GEN z; |
|
|
if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition"); |
if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition"); |
Line 1168 compimag0(GEN x, GEN y, int raw) |
|
Line 1246 compimag0(GEN x, GEN y, int raw) |
|
static GEN |
static GEN |
compreal0(GEN x, GEN y, int raw) |
compreal0(GEN x, GEN y, int raw) |
{ |
{ |
ulong tx = typ(x), av = avma; |
gpmem_t av = avma; |
|
long tx = typ(x); |
GEN z; |
GEN z; |
|
|
if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition"); |
if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition"); |
Line 1199 compraw(GEN x, GEN y) |
|
Line 1278 compraw(GEN x, GEN y) |
|
GEN |
GEN |
sqcompimag0(GEN x, long raw) |
sqcompimag0(GEN x, long raw) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN z = cgetg(4,t_QFI); |
GEN z = cgetg(4,t_QFI); |
|
|
if (typ(x)!=t_QFI) err(typeer,"composition"); |
if (typ(x)!=t_QFI) err(typeer,"composition"); |
Line 1211 sqcompimag0(GEN x, long raw) |
|
Line 1290 sqcompimag0(GEN x, long raw) |
|
static GEN |
static GEN |
sqcompreal0(GEN x, long raw) |
sqcompreal0(GEN x, long raw) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN z = cgetg(5,t_QFR); |
GEN z = cgetg(5,t_QFR); |
|
|
if (typ(x)!=t_QFR) err(typeer,"composition"); |
if (typ(x)!=t_QFR) err(typeer,"composition"); |
|
|
real_unit_form_by_disc(GEN D, long prec) |
real_unit_form_by_disc(GEN D, long prec) |
{ |
{ |
GEN y = cgetg(5,t_QFR), isqrtD; |
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"); |
if (typ(D) != t_INT || signe(D) <= 0) err(typeer,"real_unit_form_by_disc"); |
switch(mod4(D)) |
switch(mod4(D)) |
Line 1249 real_unit_form_by_disc(GEN D, long prec) |
|
Line 1328 real_unit_form_by_disc(GEN D, long prec) |
|
if (mod2(D) != mod2(isqrtD)) |
if (mod2(D) != mod2(isqrtD)) |
isqrtD = gerepileuptoint(av, addsi(-1,isqrtD)); |
isqrtD = gerepileuptoint(av, addsi(-1,isqrtD)); |
y[2] = (long)isqrtD; av = avma; |
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; |
y[4] = (long)realzero(prec); return y; |
} |
} |
|
|
GEN |
GEN |
real_unit_form(GEN x) |
real_unit_form(GEN x) |
{ |
{ |
long av = avma, prec; |
gpmem_t av = avma; |
|
long prec; |
GEN D; |
GEN D; |
if (typ(x) != t_QFR) err(typeer,"real_unit_form"); |
if (typ(x) != t_QFR) err(typeer,"real_unit_form"); |
prec = precision((GEN)x[4]); |
prec = precision((GEN)x[4]); |
Line 1283 imag_unit_form_by_disc(GEN D) |
|
Line 1363 imag_unit_form_by_disc(GEN D) |
|
y[3] = lshifti(D,-2); setsigne(y[3],1); |
y[3] = lshifti(D,-2); setsigne(y[3],1); |
if (isodd) |
if (isodd) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
y[3] = (long)gerepileuptoint(av, addis((GEN)y[3],1)); |
y[3] = lpileuptoint(av, addis((GEN)y[3],1)); |
} |
} |
return y; |
return y; |
} |
} |
|
|
imag_unit_form(GEN x) |
imag_unit_form(GEN x) |
{ |
{ |
GEN p1,p2, y = cgetg(4,t_QFI); |
GEN p1,p2, y = cgetg(4,t_QFI); |
long av; |
gpmem_t av; |
if (typ(x) != t_QFI) err(typeer,"imag_unit_form"); |
if (typ(x) != t_QFI) err(typeer,"imag_unit_form"); |
y[1] = un; |
y[1] = un; |
y[2] = mpodd((GEN)x[2])? un: zero; |
y[2] = mpodd((GEN)x[2])? un: zero; |
av = avma; p1 = mulii((GEN)x[1],(GEN)x[3]); |
av = avma; p1 = mulii((GEN)x[1],(GEN)x[3]); |
p2 = shifti(sqri((GEN)x[2]),-2); |
p2 = shifti(sqri((GEN)x[2]),-2); |
y[3] = (long)gerepileuptoint(av, subii(p1,p2)); |
y[3] = lpileuptoint(av, subii(p1,p2)); |
return y; |
return y; |
} |
} |
|
|
GEN |
GEN |
powrealraw(GEN x, long n) |
powrealraw(GEN x, long n) |
{ |
{ |
long av = avma, m; |
gpmem_t av = avma; |
|
long m; |
GEN y; |
GEN y; |
|
|
if (typ(x) != t_QFR) |
if (typ(x) != t_QFR) |
Line 1329 powrealraw(GEN x, long n) |
|
Line 1410 powrealraw(GEN x, long n) |
|
GEN |
GEN |
powimagraw(GEN x, long n) |
powimagraw(GEN x, long n) |
{ |
{ |
long av = avma, m; |
gpmem_t av = avma; |
|
long m; |
GEN y; |
GEN y; |
|
|
if (typ(x) != t_QFI) |
if (typ(x) != t_QFI) |
Line 1360 powraw(GEN x, long n) |
|
Line 1442 powraw(GEN x, long n) |
|
GEN |
GEN |
nucomp(GEN x, GEN y, GEN l) |
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; |
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); |
if (x==y) return nudupl(x,l); |
Line 1424 nucomp(GEN x, GEN y, GEN l) |
|
Line 1507 nucomp(GEN x, GEN y, GEN l) |
|
GEN |
GEN |
nudupl(GEN x, GEN l) |
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; |
GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g; |
|
|
if (typ(x) != t_QFI) |
if (typ(x) != t_QFI) |
Line 1464 nudupl(GEN x, GEN l) |
|
Line 1548 nudupl(GEN x, GEN l) |
|
GEN |
GEN |
nupow(GEN x, GEN n) |
nupow(GEN x, GEN n) |
{ |
{ |
long av,tetpil,i,j; |
gpmem_t av,tetpil; |
|
long i,j; |
unsigned long m; |
unsigned long m; |
GEN y,l; |
GEN y,l; |
|
|
Line 1551 rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD) |
|
Line 1636 rhoreal_aux(GEN x, GEN D, GEN sqrtD, GEN isqrtD) |
|
GEN p1,p2, y = cgetg(6,t_VEC); |
GEN p1,p2, y = cgetg(6,t_VEC); |
GEN b = (GEN)x[2]; |
GEN b = (GEN)x[2]; |
GEN c = (GEN)x[3]; |
GEN c = (GEN)x[3]; |
long s = signe(c); |
|
|
|
y[1] = (long)c; |
y[1] = (long)c; |
p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: c; |
p2 = (absi_cmp(isqrtD,c) >= 0)? isqrtD: absi(c); |
setsigne(c,1); |
|
p1 = shifti(c,1); |
p1 = shifti(c,1); |
|
if (p1 == gzero) err(talker, "reducible form in rhoreal"); |
|
setsigne(p1,1); /* |2c| */ |
p2 = divii(addii(p2,b), p1); |
p2 = divii(addii(p2,b), p1); |
setsigne(c,s); |
|
y[2] = lsubii(mulii(p2,p1), b); |
y[2] = lsubii(mulii(p2,p1), b); |
|
|
p1 = shifti(subii(sqri((GEN)y[2]),D),-2); |
p1 = shifti(subii(sqri((GEN)y[2]),D),-2); |
|
|
codeform5(GEN x, long prec) |
codeform5(GEN x, long prec) |
{ |
{ |
GEN y = cgetg(6,t_VEC); |
GEN y = cgetg(6,t_VEC); |
y[1]=x[1]; |
y[1] = x[1]; |
y[2]=x[2]; |
y[2] = x[2]; |
y[3]=x[3]; |
y[3] = x[3]; |
y[4]=zero; |
y[4] = zero; |
y[5]=(long)realun(prec); return y; |
y[5] = (long)realun(prec); return y; |
} |
} |
|
|
static GEN |
static GEN |
Line 1619 decodeform(GEN x, GEN d0) |
|
Line 1703 decodeform(GEN x, GEN d0) |
|
p1 = shiftr(p1,-e); |
p1 = shiftr(p1,-e); |
p2 = addis(mulsi(EXP220,p2), e); |
p2 = addis(mulsi(EXP220,p2), e); |
p1 = mplog(p1); |
p1 = mplog(p1); |
p1 = mpadd(p1, mulir(p2, glog(gdeux, lg(d0)))); |
p1 = mpadd(p1, mulir(p2, mplog2(lg(d0)))); |
} |
} |
else |
else |
{ /* to avoid loss of precision */ |
{ /* to avoid loss of precision */ |
Line 1665 redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD) |
|
Line 1749 redrealform5(GEN x, GEN D, GEN sqrtD, GEN isqrtD) |
|
static GEN |
static GEN |
redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD) |
redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD) |
{ |
{ |
long av = avma, prec; |
gpmem_t av = avma; |
|
long prec; |
GEN d0; |
GEN d0; |
|
|
if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal"); |
if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal"); |
Line 1703 redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrt |
|
Line 1788 redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrt |
|
GEN |
GEN |
comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqrtD) |
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); |
GEN z = cgetg(6,t_VEC); comp_gen(z,x,y); |
if (x == y) |
if (x == y) |
{ |
{ |
Line 1723 comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqr |
|
Line 1808 comprealform5(GEN x, GEN y, GEN D, GEN sqrtD, GEN isqr |
|
GEN |
GEN |
powrealform(GEN x, GEN n) |
powrealform(GEN x, GEN n) |
{ |
{ |
long av = avma, i,m; |
gpmem_t av = avma; |
|
long i,m; |
GEN y,D,sqrtD,isqrtD,d0; |
GEN y,D,sqrtD,isqrtD,d0; |
|
|
if (typ(x) != t_QFR) |
if (typ(x) != t_QFR) |
Line 1755 powrealform(GEN x, GEN n) |
|
Line 1841 powrealform(GEN x, GEN n) |
|
GEN |
GEN |
redimag(GEN x) |
redimag(GEN x) |
{ |
{ |
long av=avma, fl; |
gpmem_t av=avma; |
|
long fl; |
do x = rhoimag0(x, &fl); while (fl == 0); |
do x = rhoimag0(x, &fl); while (fl == 0); |
x = gerepilecopy(av,x); |
x = gerepilecopy(av,x); |
if (fl == 2) setsigne(x[2], -signe(x[2])); |
if (fl == 2) setsigne(x[2], -signe(x[2])); |
|
|
qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD) |
qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD) |
{ |
{ |
long tx=typ(x),fl; |
long tx=typ(x),fl; |
ulong av; |
gpmem_t av; |
|
|
if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred"); |
if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred"); |
if (!D) D = qf_disc(x,NULL,NULL); |
if (!D) D = qf_disc(x,NULL,NULL); |
Line 1815 qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD |
|
Line 1902 qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD |
|
GEN |
GEN |
primeform(GEN x, GEN p, long prec) |
primeform(GEN x, GEN p, long prec) |
{ |
{ |
long av,tetpil,s=signe(x); |
gpmem_t av; |
GEN y,b; |
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 (typ(p) != t_INT || signe(p) <= 0) err(arither1); |
if (is_pm1(p)) |
if (is_pm1(p)) |
return s<0? imag_unit_form_by_disc(x) |
return sx<0? imag_unit_form_by_disc(x) |
: real_unit_form_by_disc(x,prec); |
: real_unit_form_by_disc(x,prec); |
if (s < 0) |
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 |
else |
{ |
{ |
y = cgetg(5,t_QFR); s = mod8(x); |
y = cgetg(5, t_QFR); |
y[4]=(long)realzero(prec); |
y[4] = (long)realzero(prec); |
} |
} |
switch(s&3) |
switch(s&3) |
{ |
{ |
case 2: case 3: err(funder2,"primeform"); |
case 2: case 3: err(funder2,"primeform"); |
} |
} |
y[1] = (long)icopy(p); av = avma; |
av = avma; |
if (egalii(p,gdeux)) |
if (egalii(p, gdeux)) |
{ |
{ |
switch(s) |
switch(s) |
{ |
{ |
case 0: y[2]=zero; break; |
case 0: b = gzero; break; |
case 8: s=0; y[2]=zero; break; |
case 1: b = gun; break; |
case 1: s=1; y[2]=un; break; |
case 4: b = gdeux; break; |
case 4: s=4; y[2]=deux; break; |
default: err(sqrter5); b = NULL; /* -Wall */ |
default: err(sqrter5); |
|
} |
} |
b=subsi(s,x); tetpil=avma; b=shifti(b,-3); |
c = shifti(subsi(s,x), -3); |
} |
} |
else |
else |
{ |
{ |
b = mpsqrtmod(x,p); if (!b) err(sqrter5); |
b = mpsqrtmod(x,p); if (!b) err(sqrter5); |
if (mod2(b) == mod2(x)) y[2] = (long)b; |
s &= 1; /* s = x mod 2 */ |
else |
/* mod(b) != mod2(x) ? [Warning: we may have b == 0] */ |
{ tetpil = avma; y[2] = lpile(av,tetpil,subii(p,b)); } |
if ((!signe(b) && s) || mod2(b) != s) b = gerepileuptoint(av, subii(p,b)); |
|
|
av=avma; b=shifti(subii(sqri((GEN)y[2]),x),-2); |
av = avma; |
tetpil=avma; b=divii(b,p); |
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; |
} |
} |
|
|
/*********************************************************************/ |
/*********************************************************************/ |
Line 1955 bittest(GEN x, long n) |
|
Line 2046 bittest(GEN x, long n) |
|
return n? 1: 0; |
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<<b2) - 1); |
|
if (!nonzero) { /* Increment res. Do not underflow to before p */ |
|
pnew = p + l_res - 2; |
|
while (--pnew >= p) |
|
if (++*pnew) |
|
break; |
|
} |
|
} |
|
if (partial_bits) |
|
*p &= (1<<partial_bits) - 1; |
|
pnew = p; |
|
while ((pnew < p + l_res - 2) && !*pnew) /* Skip 0 words */ |
|
pnew++; |
|
avma = (gpmem_t)(pnew - 2); |
|
if (pnew >= 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 |
static long |
bittestg(GEN x, GEN n) |
bittestg(GEN x, GEN n) |
{ |
{ |
Line 1967 gbittest(GEN x, GEN n) |
|
Line 2152 gbittest(GEN x, GEN n) |
|
return arith_proto2(bittestg,x,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 **/ |
/** BITMAP OPS **/ |
Line 2226 ibitnegimply(GEN x, GEN y) |
|
Line 2417 ibitnegimply(GEN x, GEN y) |
|
} |
} |
|
|
static GEN |
static GEN |
inegate_inplace(GEN z, long ltop) |
inegate_inplace(GEN z, gpmem_t ltop) |
{ |
{ |
long lbot, o; |
long o; |
|
|
/* Negate z */ |
/* Negate z */ |
o = incdec(z, 1); /* Can overflow z... */ |
o = incdec(z, 1); /* Can overflow z... */ |
Line 2238 inegate_inplace(GEN z, long ltop) |
|
Line 2429 inegate_inplace(GEN z, long ltop) |
|
else if (lgefint(z) == 2) |
else if (lgefint(z) == 2) |
setsigne(z, 0); |
setsigne(z, 0); |
incdec(z,-1); /* Restore a normalized value */ |
incdec(z,-1); /* Restore a normalized value */ |
lbot = avma; |
return gerepileupto(ltop, subis(z,1)); |
z = gsub(z,gun); |
|
return gerepile(ltop,lbot,z); |
|
} |
} |
|
|
GEN |
GEN |
gbitor(GEN x, GEN y) |
gbitor(GEN x, GEN y) |
{ |
{ |
|
gpmem_t ltop; |
long sx,sy; |
long sx,sy; |
long ltop; |
|
GEN z; |
GEN z; |
|
|
if (typ(x) != t_INT || typ(y) != t_INT) |
if (typ(x) != t_INT || typ(y) != t_INT) |
Line 2278 gbitor(GEN x, GEN y) |
|
Line 2467 gbitor(GEN x, GEN y) |
|
GEN |
GEN |
gbitand(GEN x, GEN y) |
gbitand(GEN x, GEN y) |
{ |
{ |
|
gpmem_t ltop; |
long sx,sy; |
long sx,sy; |
long ltop; |
|
GEN z; |
GEN z; |
|
|
if (typ(x) != t_INT || typ(y) != t_INT) |
if (typ(x) != t_INT || typ(y) != t_INT) |
Line 2312 gbitand(GEN x, GEN y) |
|
Line 2501 gbitand(GEN x, GEN y) |
|
GEN |
GEN |
gbitxor(GEN x, GEN y) |
gbitxor(GEN x, GEN y) |
{ |
{ |
|
gpmem_t ltop; |
long sx,sy; |
long sx,sy; |
long ltop; |
|
GEN z; |
GEN z; |
|
|
if (typ(x) != t_INT || typ(y) != t_INT) |
if (typ(x) != t_INT || typ(y) != t_INT) |
Line 2345 gbitxor(GEN x, GEN y) |
|
Line 2534 gbitxor(GEN x, GEN y) |
|
GEN |
GEN |
gbitnegimply(GEN x, GEN y) /* x & ~y */ |
gbitnegimply(GEN x, GEN y) /* x & ~y */ |
{ |
{ |
|
gpmem_t ltop; |
long sx,sy; |
long sx,sy; |
long ltop; |
|
GEN z; |
GEN z; |
|
|
if (typ(x) != t_INT || typ(y) != t_INT) |
if (typ(x) != t_INT || typ(y) != t_INT) |