version 1.1, 2001/10/02 11:17:01 |
version 1.2, 2002/09/11 07:26:48 |
Line 115 garith_proto2gs(GEN f(GEN,long), GEN x, long y) |
|
Line 115 garith_proto2gs(GEN f(GEN,long), GEN x, long y) |
|
} |
} |
|
|
GEN |
GEN |
|
garith_proto3ggs(GEN f(GEN,GEN,long), GEN x, GEN y, long z) |
|
{ |
|
long l,i,tx = typ(x); |
|
GEN t; |
|
|
|
if (is_matvec_t(tx)) |
|
{ |
|
l=lg(x); t=cgetg(l,tx); |
|
for (i=1; i<l; i++) t[i]= (long) garith_proto3ggs(f,(GEN) x[i],y,z); |
|
return t; |
|
} |
|
if (tx != t_INT) err(arither1); |
|
tx = typ(y); |
|
if (is_matvec_t(tx)) |
|
{ |
|
l=lg(y); t=cgetg(l,tx); |
|
for (i=1; i<l; i++) t[i]= (long) garith_proto3ggs(f,x,(GEN) y[i],z); |
|
return t; |
|
} |
|
if (tx != t_INT) err(arither1); |
|
return f(x,y,z); |
|
} |
|
|
|
GEN |
gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y) |
gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y) |
{ |
{ |
int tx=typ(x); |
int tx=typ(x); |
if (!y) |
if (!y) |
{ |
{ |
ulong av=avma; |
gpmem_t av = avma; |
if (tx!=t_VEC && tx!=t_COL) |
if (tx!=t_VEC && tx!=t_COL) |
err(typeer,"association"); |
err(typeer,"association"); |
return gerepileupto(av,divide_conquer_prod(x,f)); |
return gerepileupto(av,divide_conquer_prod(x,f)); |
Line 137 gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y) |
|
Line 161 gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y) |
|
GEN |
GEN |
order(GEN x) |
order(GEN x) |
{ |
{ |
long av=avma,av1,i,e; |
gpmem_t av = avma,av1; |
|
long i,e; |
GEN o,m,p; |
GEN o,m,p; |
|
|
if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2]))) |
if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2]))) |
|
|
GEN |
GEN |
gener(GEN m) |
gener(GEN m) |
{ |
{ |
long av=avma,av1,k,i,e; |
gpmem_t av=avma,av1; |
|
long k,i,e; |
GEN x,t,q,p; |
GEN x,t,q,p; |
|
|
if (typ(m) != t_INT) err(arither1); |
if (typ(m) != t_INT) err(arither1); |
e = signe(m); |
e = signe(m); |
if (!e) err(talker,"zero modulus in znprimroot"); |
if (!e) err(talker,"zero modulus in znprimroot"); |
if (is_pm1(m)) { avma=av; return gmodulss(0,1); } |
if (is_pm1(m)) return gmodulss(0,1); |
if (e<0) m = absi(m); |
if (e < 0) m = absi(m); |
|
|
e = mod4(m); |
e = mod4(m); |
if (e == 0) /* m = 0 mod 4 */ |
if (e == 0) /* m = 0 mod 4 */ |
|
|
|
|
t=decomp(m); if (lg(t[1]) != 2) err(generer); |
t=decomp(m); if (lg(t[1]) != 2) err(generer); |
p=gcoeff(t,1,1); e=itos(gcoeff(t,1,2)); q=subis(p,1); |
p=gcoeff(t,1,1); e=itos(gcoeff(t,1,2)); q=subis(p,1); |
if (e>=2) |
if (e >= 2) |
{ |
{ |
x = (GEN) gener(p)[2]; |
x = (GEN)gener(p)[2]; |
if (gcmp1(powmodulo(x,q, sqri(p)))) x = addii(x,p); |
if (gcmp1(powmodulo(x,q, sqri(p)))) x = addii(x,p); |
av1=avma; return gerepile(av,av1,gmodulcp(x,m)); |
av1=avma; return gerepile(av,av1,gmodulcp(x,m)); |
} |
} |
|
|
p=(GEN)decomp(q)[1]; k=lg(p)-1; x=stoi(1); |
p=(GEN)decomp(q)[1]; k=lg(p)-1; x=stoi(1); |
|
for (i=1; i<=k; i++) p[i] = (long)diviiexact(q, (GEN)p[i]); |
for(;;) |
for(;;) |
{ |
{ |
x[2]++; |
x[2]++; |
if (gcmp1(mppgcd(m,x))) |
if (gcmp1(mppgcd(m,x))) |
{ |
{ |
for (i=k; i; i--) |
for (i=k; i; i--) |
if (gcmp1(powmodulo(x, divii(q,(GEN)p[i]), m))) break; |
if (gcmp1(powmodulo(x, (GEN)p[i], m))) break; |
if (!i) break; |
if (!i) break; |
} |
} |
} |
} |
av1=avma; return gerepile(av,av1,gmodulcp(x,m)); |
av1=avma; return gerepile(av,av1,gmodulcp(x,m)); |
} |
} |
|
|
|
/* assume p odd prime */ |
|
ulong |
|
u_gener(ulong p) |
|
{ |
|
const gpmem_t av = avma; |
|
const long q = p - 1; |
|
const GEN L = (GEN)decomp(utoi(q))[1]; |
|
const int k = lg(L) - 1; |
|
int i,x; |
|
|
|
for (x=2;;x++) |
|
if (x % p) |
|
{ |
|
for (i=k; i; i--) |
|
if (powuumod(x, q/itos((GEN)L[i]), p) == 1) break; |
|
if (!i) break; |
|
} |
|
avma = av; return x; |
|
} |
|
|
GEN |
GEN |
znstar(GEN n) |
znstar(GEN n) |
{ |
{ |
GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a; |
GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a; |
long i,j,nbp,sizeh; |
long i,j,nbp,sizeh; |
ulong av; |
gpmem_t av; |
|
|
if (typ(n) != t_INT) err(arither1); |
if (typ(n) != t_INT) err(arither1); |
if (!signe(n)) |
if (!signe(n)) |
|
|
static GEN |
static GEN |
racine_r(GEN a, long l) |
racine_r(GEN a, long l) |
{ |
{ |
long av,s; |
gpmem_t av; |
|
long s; |
GEN x,y,z; |
GEN x,y,z; |
|
|
if (l <= 4) return utoi(mpsqrtl(a)); |
if (l <= 4) return utoi(mpsqrtl(a)); |
Line 331 racine_r(GEN a, long l) |
|
Line 379 racine_r(GEN a, long l) |
|
y = x; x = z; |
y = x; x = z; |
} |
} |
while (cmpii(x,y) < 0); |
while (cmpii(x,y) < 0); |
avma = (long)y; |
avma = (gpmem_t)y; |
return gerepileuptoint(av,y); |
return gerepileuptoint(av,y); |
} |
} |
|
|
GEN |
GEN |
racine_i(GEN a, int roundup) |
racine_i(GEN a, int roundup) |
{ |
{ |
long k,m,l = lgefint(a), av = avma; |
gpmem_t av = avma; |
|
long k,m,l = lgefint(a); |
GEN x = racine_r(a, l); |
GEN x = racine_r(a, l); |
if (roundup && signe(x)) |
if (roundup && signe(x)) |
{ |
{ |
m = modBIL(x); |
m = modBIL(x); |
k = (m * m != a[l-1] || !egalii(sqri(x),a)); |
k = (m * m != a[l-1] || !egalii(sqri(x),a)); |
avma = (long)x; |
avma = (gpmem_t)x; |
if (k) x = gerepileuptoint(av, addis(x,1)); |
if (k) x = gerepileuptoint(av, addis(x,1)); |
} |
} |
return x; |
return x; |
Line 409 ucarrecomplet(ulong A) |
|
Line 458 ucarrecomplet(ulong A) |
|
long |
long |
carrecomplet(GEN x, GEN *pt) |
carrecomplet(GEN x, GEN *pt) |
{ |
{ |
long av; |
gpmem_t av; |
GEN y; |
GEN y; |
|
|
switch(signe(x)) |
switch(signe(x)) |
Line 427 carrecomplet(GEN x, GEN *pt) |
|
Line 476 carrecomplet(GEN x, GEN *pt) |
|
if (!carremod((ulong)smodis(x, 64*63*65*11))) return 0; |
if (!carremod((ulong)smodis(x, 64*63*65*11))) return 0; |
av=avma; y = racine(x); |
av=avma; y = racine(x); |
if (!egalii(sqri(y),x)) { avma=av; return 0; } |
if (!egalii(sqri(y),x)) { avma=av; return 0; } |
if (pt) { *pt = y; avma=(long)y; } else avma=av; |
if (pt) { *pt = y; avma=(gpmem_t)y; } else avma=av; |
return 1; |
return 1; |
} |
} |
|
|
static GEN |
static GEN |
polcarrecomplet(GEN x, GEN *pt) |
polcarrecomplet(GEN x, GEN *pt) |
{ |
{ |
long i,l,av,av2; |
gpmem_t av,av2; |
|
long i,l; |
GEN y,a,b; |
GEN y,a,b; |
|
|
if (!signe(x)) return gun; |
if (!signe(x)) return gun; |
Line 493 gcarrecomplet(GEN x, GEN *pt) |
|
Line 543 gcarrecomplet(GEN x, GEN *pt) |
|
GEN |
GEN |
gcarreparfait(GEN x) |
gcarreparfait(GEN x) |
{ |
{ |
|
gpmem_t av; |
GEN p1,a,p; |
GEN p1,a,p; |
long tx=typ(x),l,i,av,v; |
long tx=typ(x),l,i,v; |
|
|
switch(tx) |
switch(tx) |
{ |
{ |
Line 506 gcarreparfait(GEN x) |
|
Line 557 gcarreparfait(GEN x) |
|
|
|
case t_INTMOD: |
case t_INTMOD: |
{ |
{ |
GEN b, q, e; |
GEN b, q; |
long w; |
long w; |
a = (GEN)x[2]; if (!signe(a)) return gun; |
a = (GEN)x[2]; if (!signe(a)) return gun; |
av = avma; |
av = avma; |
Line 533 gcarreparfait(GEN x) |
|
Line 584 gcarreparfait(GEN x) |
|
if (i==0) |
if (i==0) |
{ |
{ |
GEN d = mppgcd(a,q); |
GEN d = mppgcd(a,q); |
p1 = factor(d); |
p = (GEN)factor(d)[1]; l = lg(p); |
p = (GEN)p1[1]; |
|
e = (GEN)p1[2]; l = lg(e); |
|
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
{ |
{ |
v = pvaluation(a,(GEN)p[i],&p1); |
v = pvaluation(a,(GEN)p[i],&p1); |
w = itos((GEN)e[i]); |
w = pvaluation(q,(GEN)p[i], &q); |
if (v < w && (v&1 || kronecker(p1,(GEN)p[i]) == -1)) |
if (v < w && (v&1 || kronecker(p1,(GEN)p[i]) == -1)) |
{ avma = av; return gzero; } |
{ avma = av; return gzero; } |
q = diviiexact(q, gpowgs((GEN)p[i], w)); |
|
} |
} |
if (kronecker(a,q) == -1) { avma = av; return gzero; } |
if (kronecker(a,q) == -1) { avma = av; return gzero; } |
} |
} |
/* kro(a,q) = 1, q odd: need to factor q */ |
/* kro(a,q) = 1, q odd: need to factor q */ |
p1 = factor(q); |
p = (GEN)factor(q)[1]; l = lg(p) - 1; |
p = (GEN)p1[1]; |
|
e = (GEN)p1[2]; l = lg(e) - 1; |
|
/* kro(a,q) = 1, check all p|q but the last (product formula) */ |
/* kro(a,q) = 1, check all p|q but the last (product formula) */ |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
if (kronecker(a,(GEN)p[i]) == -1) { avma = av; return gzero; } |
if (kronecker(a,(GEN)p[i]) == -1) { avma = av; return gzero; } |
Line 618 gkronecker(GEN x, GEN y) |
|
Line 664 gkronecker(GEN x, GEN y) |
|
long |
long |
kronecker(GEN x, GEN y) |
kronecker(GEN x, GEN y) |
{ |
{ |
|
const gpmem_t av = avma; |
GEN z; |
GEN z; |
long av,r,s=1; |
long r,s=1; |
|
|
av=avma; |
|
switch (signe(y)) |
switch (signe(y)) |
{ |
{ |
case -1: y=negi(y); if (signe(x)<0) s= -1; break; |
case -1: y=negi(y); if (signe(x)<0) s= -1; break; |
Line 662 gkrogs(GEN x, long y) |
|
Line 708 gkrogs(GEN x, long y) |
|
long |
long |
krogs(GEN x, long y) |
krogs(GEN x, long y) |
{ |
{ |
long av,r,s=1,x1,z; |
const gpmem_t av = avma; |
|
long r,s=1,x1,z; |
|
|
av=avma; |
|
if (y<=0) |
if (y<=0) |
{ |
{ |
if (y) { y = -y; if (signe(x)<0) s = -1; } |
if (y) { y = -y; if (signe(x)<0) s = -1; } |
Line 698 krogs(GEN x, long y) |
|
Line 744 krogs(GEN x, long y) |
|
long |
long |
krosg(long s, GEN x) |
krosg(long s, GEN x) |
{ |
{ |
long av = avma, y = kronecker(stoi(s),x); |
gpmem_t av = avma; |
|
long y = kronecker(stoi(s),x); |
avma = av; return y; |
avma = av; return y; |
} |
} |
|
|
Line 747 hil0(GEN x, GEN y, GEN p) |
|
Line 794 hil0(GEN x, GEN y, GEN p) |
|
long |
long |
hil(GEN x, GEN y, GEN p) |
hil(GEN x, GEN y, GEN p) |
{ |
{ |
long a,b,av,tx,ty,z; |
gpmem_t av; |
|
long a,b,tx,ty,z; |
GEN p1,p2,u,v; |
GEN p1,p2,u,v; |
|
|
if (gcmp0(x) || gcmp0(y)) return 0; |
if (gcmp0(x) || gcmp0(y)) return 0; |
Line 779 hil(GEN x, GEN y, GEN p) |
|
Line 827 hil(GEN x, GEN y, GEN p) |
|
case t_REAL: |
case t_REAL: |
return (signe(x)<0 && signe(y)<0)? -1: 1; |
return (signe(x)<0 && signe(y)<0)? -1: 1; |
case t_INTMOD: |
case t_INTMOD: |
if (egalii(gdeux,(GEN)y[1])) err(hiler1); |
p = (GEN)y[1]; |
return hil(x,(GEN)y[2],(GEN)y[1]); |
if (egalii(gdeux,p)) err(hiler1); |
|
return hil(x,(GEN)y[2],p); |
case t_FRAC: case t_FRACN: |
case t_FRAC: case t_FRACN: |
p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p); |
p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p); |
avma=av; return z; |
avma=av; return z; |
case t_PADIC: |
case t_PADIC: |
if (egalii(gdeux,(GEN)y[2]) && precp(y) <= 2) err(hiler1); |
p = (GEN)y[2]; |
p1 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4]; |
if (egalii(gdeux,p) && precp(y) <= 1) err(hiler1); |
z=hil(x,p1,(GEN)y[2]); avma=av; return z; |
p1 = odd(valp(y))? mulii(p,(GEN)y[4]): (GEN)y[4]; |
|
z=hil(x,p1,p); avma=av; return z; |
} |
} |
break; |
break; |
|
|
Line 797 hil(GEN x, GEN y, GEN p) |
|
Line 847 hil(GEN x, GEN y, GEN p) |
|
return signe(y[1])*signe(y[2]); |
return signe(y[1])*signe(y[2]); |
|
|
case t_INTMOD: |
case t_INTMOD: |
if (egalii(gdeux,(GEN)y[1])) err(hiler1); |
p = (GEN)x[1]; |
|
if (egalii(gdeux,p)) err(hiler1); |
switch(ty) |
switch(ty) |
{ |
{ |
case t_INTMOD: |
case t_INTMOD: |
if (!egalii((GEN)x[1],(GEN)y[1])) break; |
if (!egalii(p, (GEN)y[1])) break; |
return hil((GEN)x[2],(GEN)y[2],(GEN)x[1]); |
return hil((GEN)x[2],(GEN)y[2],p); |
case t_FRAC: case t_FRACN: |
case t_FRAC: case t_FRACN: |
return hil((GEN)x[2],y,(GEN)x[1]); |
return hil((GEN)x[2],y,p); |
case t_PADIC: |
case t_PADIC: |
if (!egalii((GEN)x[1],(GEN)y[2])) break; |
if (!egalii(p, (GEN)y[2])) break; |
return hil((GEN)x[2],y,(GEN)x[1]); |
return hil((GEN)x[2],y,p); |
} |
} |
break; |
break; |
|
|
Line 824 hil(GEN x, GEN y, GEN p) |
|
Line 875 hil(GEN x, GEN y, GEN p) |
|
break; |
break; |
|
|
case t_PADIC: |
case t_PADIC: |
if (ty != t_PADIC || !egalii((GEN)x[2],(GEN)y[2])) break; |
p = (GEN)x[2]; |
p1 = odd(valp(x))? mulii((GEN)x[2],(GEN)x[4]): (GEN)x[4]; |
if (ty != t_PADIC || !egalii(p,(GEN)y[2])) break; |
p2 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4]; |
if (egalii(p, gdeux) && (precp(x) <= 1 || precp(y) <= 1)) err(hiler1); |
z=hil(p1,p2,(GEN)x[2]); avma=av; return z; |
p1 = odd(valp(x))? mulii(p,(GEN)x[4]): (GEN)x[4]; |
|
p2 = odd(valp(y))? mulii(p,(GEN)y[4]): (GEN)y[4]; |
|
z=hil(p1,p2,p); avma=av; return z; |
} |
} |
err(talker,"forbidden or incompatible types in hil"); |
err(talker,"forbidden or incompatible types in hil"); |
return 0; /* not reached */ |
return 0; /* not reached */ |
Line 844 hil(GEN x, GEN y, GEN p) |
|
Line 897 hil(GEN x, GEN y, GEN p) |
|
|
|
#define sqrmod(x,p) (resii(sqri(x),p)) |
#define sqrmod(x,p) (resii(sqri(x),p)) |
|
|
/* Assume p is prime. return NULL if (a,p) = -1 */ |
static GEN ffsqrtmod(GEN a, GEN p); |
|
|
|
/* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */ |
GEN |
GEN |
mpsqrtmod(GEN a, GEN p) |
mpsqrtmod(GEN a, GEN p) |
{ |
{ |
long av = avma, av1,lim,i,k,e; |
gpmem_t av = avma, av1,lim; |
|
long i,k,e; |
GEN p1,q,v,y,w,m; |
GEN p1,q,v,y,w,m; |
|
|
if (typ(a) != t_INT || typ(p) != t_INT) err(arither1); |
if (typ(a) != t_INT || typ(p) != t_INT) err(arither1); |
if (signe(p) <= 0 || is_pm1(p)) err(talker,"not a prime in mpsqrtmod"); |
if (signe(p) <= 0 || is_pm1(p)) err(talker,"not a prime in mpsqrtmod"); |
p1 = addsi(-1,p); e = vali(p1); |
p1 = addsi(-1,p); e = vali(p1); |
|
|
|
/* If `e' is "too big", use Cipolla algorithm ! [GTL] */ |
|
if (e*(e-1) > 20 + 8 * bit_accuracy(lgefint(p))) |
|
{ |
|
v = ffsqrtmod(a,p); |
|
if (!v) { avma = av; return NULL; } |
|
return gerepileuptoint(av,v); |
|
} |
|
|
if (e == 0) /* p = 2 */ |
if (e == 0) /* p = 2 */ |
{ |
{ |
avma = av; |
avma = av; |
Line 911 mpsqrtmod(GEN a, GEN p) |
|
Line 976 mpsqrtmod(GEN a, GEN p) |
|
return gerepileuptoint(av, v); |
return gerepileuptoint(av, v); |
} |
} |
|
|
|
/* Cipolla's algorithm is better when e = v_2(p-1) is "too big". |
|
* Otherwise, is a constant times worse than the above one. |
|
* For p = 3 (mod 4), is about 3 times worse, and in average |
|
* is about 2 or 2.5 times worse. |
|
* |
|
* But try both algorithms for e.g. S(n)=(2^n+3)^2-8 with |
|
* n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc. |
|
* |
|
* If X^2 = t^2 - a is not a square in F_p, then |
|
* |
|
* (t+X)^(p+1) = (t+X)(t-X) = a |
|
* |
|
* so we get sqrt(a) in F_p^2 by |
|
* |
|
* sqrt(a) = (t+X)^((p+1)/2) |
|
* |
|
* If (a|p)=1, then sqrt(a) is in F_p. |
|
* |
|
* cf: LNCS 2286, pp 430-434 (2002) [Gonzalo Tornaria] */ |
|
static GEN |
|
ffsqrtmod(GEN a, GEN p) |
|
{ |
|
gpmem_t av = avma, av1, lim; |
|
long e, t, man, k, nb; |
|
GEN q, n, u, v, d, d2, b, u2, v2, qptr; |
|
|
|
if (kronecker(a, p) < 0) return NULL; |
|
|
|
av1 = avma; |
|
for(t=1; ; t++) |
|
{ |
|
n = subsi(t*t, a); |
|
if (kronecker(n, p) < 0) break; |
|
avma = av1; |
|
} |
|
|
|
u = stoi(t); v = gun; /* u+vX = t+X */ |
|
e = vali(addsi(-1,p)); q = shifti(p, -e); |
|
/* p = 2^e q + 1 and (p+1)/2 = 2^(e-1)q + 1 */ |
|
|
|
/* Compute u+vX := (t+X)^q */ |
|
av1 = avma; lim = stack_lim(av1, 1); |
|
qptr = q+2; man = *qptr; |
|
k = 1+bfffo(man); man<<=k; k=BITS_IN_LONG-k; |
|
for (nb=lgefint(q)-2;nb; nb--) |
|
{ |
|
for (; k; man<<=1, k--) |
|
{ |
|
if (man < 0) |
|
{ /* u+vX := (u+vX)^2 (t+X) */ |
|
d = addii(u, mulsi(t,v)); |
|
d2 = sqri(d); |
|
b = resii(mulii(a,v), p); |
|
u = modii(subii(mulsi(t,d2), mulii(b,addii(u,d))), p); |
|
v = modii(subii(d2, mulii(b,v)), p); |
|
} |
|
else |
|
{ /* u+vX := (u+vX)^2 */ |
|
u2 = sqri(u); |
|
v2 = sqri(v); |
|
v = modii(subii(sqri(addii(v,u)), addii(u2,v2)), p); |
|
u = modii(addii(u2, mulii(v2,n)), p); |
|
} |
|
} |
|
|
|
if (low_stack(lim, stack_lim(av1, 1))) |
|
{ |
|
GEN *gptr[2]; gptr[0]=&u; gptr[1]=&v; |
|
if (DEBUGMEM>1) err(warnmem, "ffsqrtmod"); |
|
gerepilemany(av1,gptr,2); |
|
} |
|
|
|
man = *++qptr; k = BITS_IN_LONG; |
|
} |
|
|
|
while (--e) |
|
{ /* u+vX := (u+vX)^2 */ |
|
u2 = sqri(u); |
|
v2 = sqri(v); |
|
v = modii(subii(sqri(addii(v,u)), addii(u2,v2)), p); |
|
u = modii(addii(u2, mulii(v2,n)), p); |
|
|
|
if (low_stack(lim, stack_lim(av1, 1))) |
|
{ |
|
GEN *gptr[2]; gptr[0]=&u; gptr[1]=&v; |
|
if (DEBUGMEM>1) err(warnmem, "ffsqrtmod"); |
|
gerepilemany(av1,gptr,2); |
|
} |
|
} |
|
|
|
/* Now u+vX = (t+X)^(2^(e-1)q); thus |
|
* |
|
* (u+vX)(t+X) = sqrt(a) + 0 X |
|
* |
|
* Whence, |
|
* |
|
* sqrt(a) = (u+vt)t - v*a |
|
* 0 = (u+vt) |
|
* |
|
* Thus a square root is v*a */ |
|
|
|
v = modii(mulii(v,a), p); |
|
|
|
u = subii(p,v); if (cmpii(v,u) > 0) v = u; |
|
return gerepileuptoint(av,v); |
|
} |
|
|
/*******************************************************************/ |
/*******************************************************************/ |
/* */ |
/* */ |
/* n-th ROOT MODULO p */ |
/* n-th ROOT MODULO p */ |
Line 930 mpsqrtmod(GEN a, GEN p) |
|
Line 1102 mpsqrtmod(GEN a, GEN p) |
|
static GEN |
static GEN |
mplgenmod(GEN l, long e, GEN r,GEN p,GEN *zeta) |
mplgenmod(GEN l, long e, GEN r,GEN p,GEN *zeta) |
{ |
{ |
ulong av1; |
const gpmem_t av1 = avma; |
GEN m,m1; |
GEN m,m1; |
long k,i; |
long k,i; |
av1 = avma; |
|
for (k=1; ; k++) |
for (k=1; ; k++) |
{ |
{ |
m1 = m = powmodulo(stoi(k+1),r,p); |
m1 = m = powmodulo(stoi(k+1),r,p); |
|
if (gcmp1(m)) {avma=av1; continue;} |
for (i=1; i<e; i++) |
for (i=1; i<e; i++) |
if (gcmp1(m=powmodulo(m,l,p))) break; |
if (gcmp1(m=powmodulo(m,l,p))) break; |
if (i==e) break; |
if (i==e) break; |
Line 960 GEN Fp_shanks(GEN x,GEN g0,GEN p, GEN q); |
|
Line 1132 GEN Fp_shanks(GEN x,GEN g0,GEN p, GEN q); |
|
static GEN |
static GEN |
mpsqrtlmod(GEN a, GEN l, GEN p, GEN q,long e, GEN r, GEN y, GEN m) |
mpsqrtlmod(GEN a, GEN l, GEN p, GEN q,long e, GEN r, GEN y, GEN m) |
{ |
{ |
ulong av = avma, tetpil,lim; |
gpmem_t av = avma, tetpil,lim; |
long k; |
long k; |
GEN p1,u1,u2,v,w,z,dl; |
GEN p1,u1,u2,v,w,z,dl; |
/* y contient un generateur de (Z/pZ)^* eleve a la puis (p-1)/(l^e) */ |
/* y contient un generateur de (Z/pZ)^* eleve a la puis (p-1)/(l^e) */ |
bezout(r,l,&u1,&u2); |
(void)bezout(r,l,&u1,&u2); |
v=powmodulo(a,u2,p); |
v=powmodulo(a,u2,p); |
w=powmodulo(a,modii(mulii(negi(u1),r),q),p); |
w=powmodulo(a,modii(mulii(negi(u1),r),q),p); |
lim = stack_lim(av,1); |
lim = stack_lim(av,1); |
Line 1015 If a=0 ,return 0 and if zetan is not NULL zetan is set |
|
Line 1187 If a=0 ,return 0 and if zetan is not NULL zetan is set |
|
GEN |
GEN |
mpsqrtnmod(GEN a, GEN n, GEN p, GEN *zetan) |
mpsqrtnmod(GEN a, GEN n, GEN p, GEN *zetan) |
{ |
{ |
ulong ltop=avma,lbot=0,av1,lim; |
gpmem_t ltop=avma,lbot=0,av1,lim; |
GEN m,u1,u2; |
GEN m,u1,u2; |
GEN q,z; |
GEN q,z; |
GEN *gptr[2]; |
GEN *gptr[2]; |
Line 1183 ulong mppgcd_resiu(GEN y, ulong x); |
|
Line 1355 ulong mppgcd_resiu(GEN y, ulong x); |
|
GEN |
GEN |
mppgcd(GEN a, GEN b) |
mppgcd(GEN a, GEN b) |
{ |
{ |
long av,v,w; |
long v, w; |
|
gpmem_t av; |
GEN t, p1; |
GEN t, p1; |
|
|
if (typ(a) != t_INT || typ(b) != t_INT) err(arither1); |
if (typ(a) != t_INT || typ(b) != t_INT) err(arither1); |
Line 1204 mppgcd(GEN a, GEN b) |
|
Line 1377 mppgcd(GEN a, GEN b) |
|
} |
} |
|
|
/* larger than gcd: "avma=av" gerepile (erasing t) is valid */ |
/* larger than gcd: "avma=av" gerepile (erasing t) is valid */ |
av = avma; new_chunk(lgefint(b)); /* HACK */ |
av = avma; (void)new_chunk(lgefint(b)); /* HACK */ |
t = resii(a,b); |
t = resii(a,b); |
if (!signe(t)) { avma=av; return absi(b); } |
if (!signe(t)) { avma=av; return absi(b); } |
|
|
Line 1237 mppgcd(GEN a, GEN b) |
|
Line 1410 mppgcd(GEN a, GEN b) |
|
} |
} |
} |
} |
{ |
{ |
long r[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0}; |
long r[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0}; |
r[2] = (long) ugcd((ulong)b[2], (ulong)a[2]); |
r[2] = (long) ugcd((ulong)b[2], (ulong)a[2]); |
avma = av; return shifti(r,v); |
avma = av; return shifti(r,v); |
} |
} |
Line 1246 mppgcd(GEN a, GEN b) |
|
Line 1419 mppgcd(GEN a, GEN b) |
|
GEN |
GEN |
mpppcm(GEN x, GEN y) |
mpppcm(GEN x, GEN y) |
{ |
{ |
ulong av; |
gpmem_t av; |
GEN p1,p2; |
GEN p1,p2; |
if (typ(x) != t_INT || typ(y) != t_INT) err(arither1); |
if (typ(x) != t_INT || typ(y) != t_INT) err(arither1); |
if (!signe(x)) return gzero; |
if (!signe(x)) return gzero; |
Line 1292 GEN chinese(GEN x, GEN y) |
|
Line 1465 GEN chinese(GEN x, GEN y) |
|
GEN |
GEN |
chinois(GEN x, GEN y) |
chinois(GEN x, GEN y) |
{ |
{ |
ulong av,tetpil, tx = typ(x); |
gpmem_t av,tetpil; |
long i,lx,vx; |
long i,lx,vx, tx = typ(x); |
GEN z,p1,p2,d,u,v; |
GEN z,p1,p2,d,u,v; |
|
|
if (gegal(x,y)) return gcopy(x); |
if (gegal(x,y)) return gcopy(x); |
Line 1339 chinois(GEN x, GEN y) |
|
Line 1512 chinois(GEN x, GEN y) |
|
GEN |
GEN |
chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1) |
chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN ax,p1; |
GEN ax,p1; |
(void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1)); /* HACK */ |
(void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1)); /* HACK */ |
ax = mulii(mpinvmod(x1,y1), x1); |
ax = mulii(mpinvmod(x1,y1), x1); |
Line 1364 mpinvmod(GEN a, GEN m) |
|
Line 1537 mpinvmod(GEN a, GEN m) |
|
|
|
/*********************************************************************/ |
/*********************************************************************/ |
/** **/ |
/** **/ |
/** POWERING MODULO (A^N mod M) **/ |
/** MODULAR EXPONENTIATION **/ |
/** **/ |
/** **/ |
/*********************************************************************/ |
/*********************************************************************/ |
GEN resiimul(GEN x, GEN y); |
extern ulong invrev(ulong b); |
GEN resmod2n(GEN x, long y); |
extern GEN resiimul(GEN x, GEN y); |
GEN _resii(GEN x, GEN y) { return resii(x,y); } |
|
|
|
|
static GEN _resii(GEN x, GEN y) { return resii(x,y); } |
|
|
|
/* Montgomery reduction */ |
|
|
|
typedef struct { |
|
GEN N; |
|
ulong inv; /* inv = -N^(-1) mod B, */ |
|
} montdata; |
|
|
|
void |
|
init_montdata(GEN N, montdata *s) |
|
{ |
|
s->N = N; |
|
s->inv = (ulong) -invrev(modBIL(N)); |
|
} |
|
|
GEN |
GEN |
init_remainder(GEN M) |
init_remainder(GEN M) |
{ |
{ |
GEN sM = cgetg(3, t_VEC); |
GEN sM = cgetg(3, t_VEC); |
GEN Mr = cgetr(lgefint(M) + 1); |
|
affir(M, Mr); |
|
sM[1] = (long)M; |
sM[1] = (long)M; |
sM[2] = (long)linv(Mr); |
sM[2] = (long)linv( itor(M, lgefint(M) + 1) ); /* 1. / M */ |
return sM; |
return sM; |
} |
} |
|
|
/* optimal on UltraSPARC */ |
/* optimal on UltraSPARC */ |
static long RESIIMUL_LIMIT = 150; |
static long RESIIMUL_LIMIT = 150; |
|
static long MONTGOMERY_LIMIT = 32; |
|
|
void |
void |
setresiilimit(long n) { RESIIMUL_LIMIT = n; } |
setresiilimit(long n) { RESIIMUL_LIMIT = n; } |
|
void |
|
setmontgomerylimit(long n) { MONTGOMERY_LIMIT = n; } |
|
|
|
typedef struct { |
|
GEN N; |
|
GEN (*res)(GEN,GEN); |
|
GEN (*mulred)(GEN,GEN,GEN); |
|
} muldata; |
|
|
|
/* reduction for multiplication by 2 */ |
|
static GEN |
|
_redsub(GEN x, GEN N) |
|
{ |
|
return (cmpii(x,N) >= 0)? subii(x,N): x; |
|
} |
|
/* Montgomery reduction */ |
|
extern GEN red_montgomery(GEN T, GEN N, ulong inv); |
|
static GEN |
|
montred(GEN x, GEN N) |
|
{ |
|
return red_montgomery(x, ((montdata*)N)->N, ((montdata*)N)->inv); |
|
} |
|
/* 2x mod N */ |
|
static GEN |
|
_muli2red(GEN x, GEN y/* ignored */, GEN N) { |
|
(void)y; return _redsub(shifti(x,1), N); |
|
} |
|
static GEN |
|
_muli2montred(GEN x, GEN y/* ignored */, GEN N) { |
|
GEN n = ((montdata*)N)->N; |
|
GEN z = _muli2red(x,y, n); |
|
long l = lgefint(n); |
|
while (lgefint(z) > l) z = subii(z,n); |
|
return z; |
|
} |
|
static GEN |
|
_muli2invred(GEN x, GEN y/* ignored */, GEN N) { |
|
return _muli2red(x,y, (GEN)N[1]); |
|
} |
|
/* xy mod N */ |
|
static GEN |
|
_muliired(GEN x, GEN y, GEN N) { return resii(mulii(x,y), N); } |
|
static GEN |
|
_muliimontred(GEN x, GEN y, GEN N) { return montred(mulii(x,y), N); } |
|
static GEN |
|
_muliiinvred(GEN x, GEN y, GEN N) { return resiimul(mulii(x,y), N); } |
|
|
|
static GEN |
|
_mul(void *data, GEN x, GEN y) |
|
{ |
|
muldata *D = (muldata *)data; |
|
return D->mulred(x,y,D->N); |
|
} |
|
static GEN |
|
_sqr(void *data, GEN x) |
|
{ |
|
muldata *D = (muldata *)data; |
|
return D->res(sqri(x), D->N); |
|
} |
|
|
|
/* A^k mod N */ |
GEN |
GEN |
powmodulo(GEN A, GEN N, GEN M) |
powmodulo(GEN A, GEN k, GEN N) |
{ |
{ |
|
gpmem_t av = avma; |
|
long t,s, lN; |
|
int base_is_2, use_montgomery; |
GEN y; |
GEN y; |
long av = avma,av1,lim,man,k,nb,s, *p; |
muldata D; |
GEN (*mul)(GEN,GEN) = mulii; |
montdata S; |
GEN (*res)(GEN,GEN) = _resii; |
|
|
|
if (typ(A) != t_INT || typ(N) != t_INT || typ(M) != t_INT) err(arither1); |
if (typ(A) != t_INT || typ(k) != t_INT || typ(N) != t_INT) err(arither1); |
s = signe(N); |
s = signe(k); |
if (!s) |
if (!s) |
{ |
{ |
k = signe(resii(A,M)); avma=av; |
t = signe(resii(A,N)); avma = av; |
return k? gun: gzero; |
return t? gun: gzero; |
} |
} |
if (s < 0) { A = mpinvmod(A,M); N = absi(N); } |
if (s < 0) y = mpinvmod(A,N); |
else |
else |
{ |
{ |
A = modii(A,M); |
y = modii(A,N); |
if (!signe(A)) { avma=av; return gzero; } |
if (!signe(y)) { avma = av; return gzero; } |
} |
} |
y = A; |
|
if (lgefint(y)==3) switch(y[2]) |
base_is_2 = 0; |
|
if (lgefint(y) == 3) switch(y[2]) |
{ |
{ |
case 1: /* y = 1 */ |
case 1: avma = av; return gun; |
avma=av; return gun; |
case 2: base_is_2 = 1; break; |
case 2: /* y = 2, use shifti not mulii */ |
|
mul = (GEN (*)(GEN,GEN))shifti; A = (GEN)1L; |
|
} |
} |
|
|
/* TODO: Move this out of here and use for general modular computations */ |
/* TODO: Move this out of here and use for general modular computations */ |
if ((k = vali(M)) && k == expi(M)) |
lN = lgefint(N); |
{ /* M is a power of 2 */ |
use_montgomery = mod2(N) && lN < MONTGOMERY_LIMIT; |
M = (GEN)k; |
if (use_montgomery) |
res = (GEN(*)(GEN,GEN))resmod2n; |
{ |
|
init_montdata(N, &S); |
|
y = resii(shifti(y, bit_accuracy(lN)), N); |
|
D.mulred = base_is_2? &_muli2montred: &_muliimontred; |
|
D.res = &montred; |
|
D.N = (GEN)&S; |
} |
} |
else if (lgefint(M) > RESIIMUL_LIMIT && (lgefint(N) > 3 || N[2] > 4)) |
else if (lN > RESIIMUL_LIMIT |
{ /* compute x % M using multiplication by 1./M */ |
&& (lgefint(k) > 3 || (((double)k[2])*expi(A)) > 2 + expi(N))) |
M = init_remainder(M); |
{ |
res = (GEN(*)(GEN,GEN))resiimul; |
D.mulred = base_is_2? &_muli2invred: &_muliiinvred; |
|
D.res = &resiimul; |
|
D.N = init_remainder(N); |
} |
} |
|
else |
p = N+2; man = *p; /* see puissii */ |
|
k = 1+bfffo(man); man<<=k; k = BITS_IN_LONG-k; |
|
av1=avma; lim=stack_lim(av1,1); |
|
for (nb=lgefint(N)-2;;) |
|
{ |
{ |
for (; k; man<<=1,k--) |
D.mulred = base_is_2? &_muli2red: &_muliired; |
{ |
D.res = &_resii; |
y = res(sqri(y), M); |
D.N = N; |
if (man < 0) y = res(mul(y,A), M); |
|
if (low_stack(lim, stack_lim(av1,1))) |
|
{ |
|
if (DEBUGMEM>1) err(warnmem,"powmodulo"); |
|
y = gerepileuptoint(av1,y); |
|
} |
|
} |
|
if (--nb == 0) break; |
|
man = *++p, k = BITS_IN_LONG; |
|
} |
} |
return gerepileupto(av,y); |
y = leftright_pow(y, k, (void*)&D, &_sqr, &_mul); |
|
if (use_montgomery) |
|
{ |
|
y = montred(y, (GEN)&S); |
|
if (cmpii(y,N) >= 0) y = subii(y,N); |
|
} |
|
return gerepileuptoint(av,y); |
} |
} |
|
|
/*********************************************************************/ |
/*********************************************************************/ |
Line 1457 powmodulo(GEN A, GEN N, GEN M) |
|
Line 1705 powmodulo(GEN A, GEN N, GEN M) |
|
/** **/ |
/** **/ |
/*********************************************************************/ |
/*********************************************************************/ |
GEN |
GEN |
gnextprime(GEN n) |
gnextprime(GEN n) { return garith_proto(nextprime,n,0); } |
{ |
|
return garith_proto(nextprime,n,0); |
|
} |
|
|
|
GEN |
GEN |
gprecprime(GEN n) |
gprecprime(GEN n) { return garith_proto(precprime,n,0); } |
{ |
|
return garith_proto(precprime,n,0); |
|
} |
|
|
|
GEN |
GEN |
gisprime(GEN x, long flag) |
gisprime(GEN x, long flag) |
{ |
{ |
switch (flag) |
switch (flag) |
{ |
{ |
case 0: |
case 0: return arith_proto(isprime,x,1); |
return arith_proto(isprime,x,1); |
case 1: return garith_proto2gs(plisprime,x,1); |
case 1: |
case 2: return arith_proto(isprimeAPRCL,x,1); |
return garith_proto2gs(plisprime,x,0); |
|
case 2: |
|
return garith_proto2gs(plisprime,x,1); |
|
} |
} |
err(flagerr,"gisprime"); |
err(flagerr,"gisprime"); |
return 0; |
return 0; |
} |
} |
|
|
long |
long |
isprime(GEN x) |
isprimeSelfridge(GEN x) { return (plisprime(x,0)==gun); } |
|
|
|
/* assume x BSW pseudoprime. Check whether it's small enough to be certified |
|
* prime */ |
|
int |
|
BSW_isprime_small(GEN x) |
{ |
{ |
return millerrabin(x,10); |
long l = lgefint(x); |
|
if (l < 4) return 1; |
|
if (l == 4) |
|
{ |
|
gpmem_t av = avma; |
|
long t = cmpii(x, u2toi(0x918UL, 0x4e72a000UL)); /* 10^13 */ |
|
avma = av; |
|
if (t < 0) return 1; |
|
} |
|
return 0; |
} |
} |
|
|
GEN |
/* assume x a BSW pseudoprime */ |
gispsp(GEN x) |
int |
|
BSW_isprime(GEN x) |
{ |
{ |
return arith_proto(ispsp,x,1); |
gpmem_t av = avma; |
|
long l, res; |
|
GEN F, p; |
|
|
|
if (BSW_isprime_small(x)) return 1; |
|
F = (GEN)auxdecomp(subis(x,1), 0)[1]; |
|
l = lg(F); p = (GEN)F[l-1]; |
|
if (BSW_psp(p)) |
|
{ /* smooth */ |
|
GEN z = cgetg(3, t_VEC); |
|
z[1] = (long)x; |
|
z[2] = (long)F; res = isprimeSelfridge(z); |
|
} |
|
else |
|
res = isprimeAPRCL(x); |
|
avma = av; return res; |
} |
} |
|
|
long |
long |
ispsp(GEN x) |
isprime(GEN x) |
{ |
{ |
return millerrabin(x,1); |
return BSW_psp(x) && BSW_isprime(x); |
} |
} |
|
|
GEN |
GEN |
gmillerrabin(GEN n, long k) |
gispseudoprime(GEN x, long flag) |
{ |
{ |
return arith_proto2gs(millerrabin,n,k); |
if (flag == 0) return arith_proto(BSW_psp,x,1); |
|
return gmillerrabin(x, flag); |
} |
} |
|
|
|
long |
|
ispseudoprime(GEN x, long flag) |
|
{ |
|
if (flag == 0) return BSW_psp(x); |
|
return millerrabin(x, flag); |
|
} |
|
|
|
GEN |
|
gispsp(GEN x) { return arith_proto(ispsp,x,1); } |
|
|
|
long |
|
ispsp(GEN x) { return millerrabin(x,1); } |
|
|
|
GEN |
|
gmillerrabin(GEN n, long k) { return arith_proto2gs(millerrabin,n,k); } |
|
|
/*********************************************************************/ |
/*********************************************************************/ |
/** **/ |
/** **/ |
/** FUNDAMENTAL DISCRIMINANTS **/ |
/** FUNDAMENTAL DISCRIMINANTS **/ |
/** **/ |
/** **/ |
/*********************************************************************/ |
/*********************************************************************/ |
|
|
/* gissquarefree, issquarefree moved into arith2.c with the other |
|
basic arithmetic functions (issquarefree is the square of the |
|
Moebius mu function, after all...) --GN */ |
|
|
|
GEN |
GEN |
gisfundamental(GEN x) |
gisfundamental(GEN x) { return arith_proto(isfundamental,x,1); } |
{ |
|
return arith_proto(isfundamental,x,1); |
|
} |
|
|
|
long |
long |
isfundamental(GEN x) |
isfundamental(GEN x) |
{ |
{ |
long av,r; |
long r; |
GEN p1; |
GEN p1; |
|
|
if (gcmp0(x)) return 0; |
if (gcmp0(x)) return 0; |
r=mod4(x); |
r=mod4(x); |
if (!r) |
if (!r) |
{ |
{ |
av=avma; p1=shifti(x,-2); |
const gpmem_t av = avma; |
|
p1=shifti(x,-2); |
r=mod4(p1); if (!r) return 0; |
r=mod4(p1); if (!r) return 0; |
if (signe(x)<0) r=4-r; |
if (signe(x)<0) r=4-r; |
r = (r==1)? 0: issquarefree(p1); |
r = (r==1)? 0: issquarefree(p1); |
Line 1547 isfundamental(GEN x) |
|
Line 1826 isfundamental(GEN x) |
|
GEN |
GEN |
quaddisc(GEN x) |
quaddisc(GEN x) |
{ |
{ |
long av=avma,tetpil=avma,i,r,tx=typ(x); |
const gpmem_t av = avma; |
|
long i,r,tx=typ(x); |
GEN p1,p2,f,s; |
GEN p1,p2,f,s; |
|
|
if (tx != t_INT && ! is_frac_t(tx)) err(arither1); |
if (tx != t_INT && ! is_frac_t(tx)) err(arither1); |
f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2]; |
f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2]; |
s=gun; |
s = gun; |
for (i=1; i<lg(p1); i++) |
for (i=1; i<lg(p1); i++) |
if ( odd(mael(p2,i,2)) ) { tetpil=avma; s=gmul(s,(GEN)p1[i]); } |
if (odd(mael(p2,i,2))) s = gmul(s,(GEN)p1[i]); |
r=mod4(s); if (gsigne(x)<0) r=4-r; |
r=mod4(s); if (gsigne(x)<0) r=4-r; |
if (r>1) { tetpil=avma; s=shifti(s,2); } |
if (r>1) s = shifti(s,2); |
return gerepile(av,tetpil,s); |
return gerepileuptoint(av, s); |
} |
} |
|
|
/*********************************************************************/ |
/*********************************************************************/ |
Line 1565 quaddisc(GEN x) |
|
Line 1845 quaddisc(GEN x) |
|
/** FACTORIAL **/ |
/** FACTORIAL **/ |
/** **/ |
/** **/ |
/*********************************************************************/ |
/*********************************************************************/ |
GEN muluu(ulong x, ulong y); |
|
|
|
GEN |
GEN |
mpfact(long n) |
mpfact(long n) |
{ |
{ |
long lx,k,l, av = avma; |
const gpmem_t av = avma; |
|
long lx,k,l; |
GEN x; |
GEN x; |
|
|
if (n<2) |
if (n<2) |
|
|
GEN |
GEN |
mpfactr(long n, long prec) |
mpfactr(long n, long prec) |
{ |
{ |
long av = avma, lim,k; |
gpmem_t av = avma, lim; |
GEN f = cgetr(prec); |
long k; |
|
GEN f = realun(prec); |
|
|
affsr(1,f); |
|
if (n<2) |
if (n<2) |
{ |
{ |
if (n<0) err(facter); |
if (n<0) err(facter); |
Line 1629 mpfactr(long n, long prec) |
|
Line 1908 mpfactr(long n, long prec) |
|
void |
void |
lucas(long n, GEN *ln, GEN *ln1) |
lucas(long n, GEN *ln, GEN *ln1) |
{ |
{ |
long taille,av; |
gpmem_t av; |
|
long taille; |
GEN z,t; |
GEN z,t; |
|
|
if (!n) { *ln=stoi(2); *ln1=stoi(1); return; } |
if (!n) { *ln = stoi(2); *ln1 = stoi(1); return; } |
|
|
taille=(long)(pariC3*(1+labs(n))+3); |
taille = 3 + (long)(pariC3 * (1+labs(n))); |
*ln=cgeti(taille); *ln1=cgeti(taille); |
*ln = cgeti(taille); |
av=avma; lucas(n/2,&z,&t); |
*ln1= cgeti(taille); |
|
av = avma; lucas(n/2, &z, &t); |
switch(n % 4) |
switch(n % 4) |
{ |
{ |
case -3: |
case -3: |
addsiz(2,sqri(z),*ln1); |
addsiz(2,sqri(z), *ln1); |
subiiz(addsi(1,mulii(z,t)),*ln1,*ln); break; |
subiiz(addsi(1,mulii(z,t)),*ln1, *ln); break; |
case -2: |
|
addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break; |
|
case -1: |
case -1: |
subisz(sqri(z),2,*ln1); |
subisz(sqri(z),2, *ln1); |
subiiz(subis(mulii(z,t),1),*ln1,*ln); break; |
subiiz(subis(mulii(z,t),1),*ln1, *ln); break; |
case 0: subisz(sqri(z),2,*ln); subisz(mulii(z,t),1,*ln1); break; |
case 0: subisz(sqri(z),2, *ln); subisz(mulii(z,t),1, *ln1); break; |
case 1: subisz(mulii(z,t),1,*ln); addsiz(2,sqri(t),*ln1); break; |
case 1: subisz(mulii(z,t),1, *ln); addsiz(2,sqri(t), *ln1); break; |
case 2: addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break; |
case -2: |
case 3: addsiz(1,mulii(z,t),*ln); subisz(sqri(t),2,*ln1); |
case 2: addsiz(2,sqri(z), *ln); addsiz(1,mulii(z,t), *ln1); break; |
|
case 3: addsiz(1,mulii(z,t), *ln); subisz(sqri(t),2, *ln1); |
} |
} |
avma=av; |
avma = av; |
} |
} |
|
|
GEN |
GEN |
fibo(long n) |
fibo(long n) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN ln,ln1; |
GEN ln,ln1; |
|
|
lucas(n-1,&ln,&ln1); |
lucas(n-1,&ln,&ln1); |
|
|
/* CONTINUED FRACTIONS */ |
/* CONTINUED FRACTIONS */ |
/* */ |
/* */ |
/*******************************************************************/ |
/*******************************************************************/ |
static GEN sfcont2(GEN b, GEN x, long k); |
static GEN |
|
icopy_lg(GEN x, long l) |
GEN |
|
gcf(GEN x) |
|
{ |
{ |
return sfcont(x,x,0); |
long lx = lgefint(x); |
|
GEN y; |
|
|
|
if (lx >= l) return icopy(x); |
|
y = cgeti(l); affii(x, y); return y; |
} |
} |
|
|
GEN |
/* if y != NULL, stop as soon as partial quotients differ from y */ |
gcf2(GEN b, GEN x) |
static GEN |
|
Qsfcont(GEN x, GEN y, long k) |
{ |
{ |
return contfrac0(x,b,0); |
GEN z, p1, p2, p3; |
} |
long i, l, ly; |
|
|
GEN |
ly = lgefint(x[2]); |
gboundcf(GEN x, long k) |
l = 3 + (long) ((ly-2) / pariC3); |
{ |
if (k > 0 && ++k > 0 && l > k) l = k; /* beware overflow */ |
return sfcont(x,x,k); |
if ((ulong)l > LGBITS) l = LGBITS; |
} |
|
|
|
GEN |
p1 = (GEN)x[1]; |
contfrac0(GEN x, GEN b, long flag) |
p2 = (GEN)x[2]; |
{ |
z = cgetg(l,t_VEC); |
long lb,tb,i; |
l--; |
GEN y; |
if (y) { |
|
gpmem_t av = avma; |
if (!b || gcmp0(b)) return sfcont(x,x,flag); |
if (l >= lg(y)) l = lg(y)-1; |
tb = typ(b); |
for (i = 1; i <= l; i++) |
if (tb == t_INT) return sfcont(x,x,itos(b)); |
{ |
if (! is_matvec_t(tb)) err(typeer,"contfrac0"); |
z[i] = y[i]; |
|
p3 = p2; if (!gcmp1((GEN)z[i])) p3 = mulii((GEN)z[i], p2); |
lb=lg(b); if (lb==1) return cgetg(1,t_VEC); |
p3 = subii(p1, p3); |
if (tb != t_MAT) return sfcont2(b,x,flag); |
if (signe(p3) < 0) |
if (lg(b[1])==1) return sfcont(x,x,flag); |
{ /* partial quotient too large */ |
y = (GEN) gpmalloc(lb*sizeof(long)); |
p3 = addii(p3, p2); |
for (i=1; i<lb; i++) y[i]=coeff(b,1,i); |
if (signe(p3) >= 0) i++; /* by 1 */ |
x = sfcont2(y,x,flag); free(y); return x; |
break; |
|
} |
|
if (cmpii(p3, p2) >= 0) |
|
{ /* partial quotient too small */ |
|
p3 = subii(p3, p2); |
|
if (cmpii(p3, p2) < 0) i++; /* by 1 */ |
|
break; |
|
} |
|
if ((i & 0xff) == 0) gerepileall(av, 2, &p2, &p3); |
|
p1 = p2; p2 = p3; |
|
} |
|
} else { |
|
p1 = icopy_lg(p1, ly); |
|
p2 = icopy(p2); |
|
for (i = 1; i <= l; i++) |
|
{ |
|
z[i] = (long)truedvmdii(p1,p2,&p3); |
|
if (p3 == gzero) { i++; break; } |
|
affii(p3, p1); cgiv(p3); p3 = p1; |
|
p1 = p2; p2 = p3; |
|
} |
|
} |
|
i--; |
|
if (i > 2 && gcmp1((GEN)z[i])) |
|
{ |
|
cgiv((GEN)z[i]); --i; |
|
addsiz(1,(GEN)z[i], (GEN)z[i]); |
|
} |
|
setlg(z,i+1); return z; |
} |
} |
|
|
GEN |
static GEN |
sfcont(GEN x, GEN x1, long k) |
sfcont(GEN x, long k) |
{ |
{ |
long av,lx=lg(x),tx=typ(x),e,i,j,l,l1,lx1,tetpil,f; |
gpmem_t av; |
GEN y,p1,p2,p3,p4,yp; |
long lx,tx=typ(x),e,i,l; |
|
GEN y,p1,p2,p3; |
|
|
if (is_scalar_t(tx)) |
if (is_scalar_t(tx)) |
{ |
{ |
if (gcmp0(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; } |
if (gcmp0(x)) return _vec(gzero); |
switch(tx) |
switch(tx) |
{ |
{ |
case t_INT: |
case t_INT: y = cgetg(2,t_VEC); y[1] = (long)icopy(x); return y; |
y=cgetg(2,t_VEC); y[1]=lcopy(x); return y; |
|
case t_REAL: |
case t_REAL: |
l=avma; p1=cgetg(3,t_FRACN); |
av = avma; lx = lg(x); |
p2 = rcopy(x); settyp(p2,t_INT); setlgefint(p2,lx); |
p2 = rcopy(x); settyp(p2,t_INT); setlgefint(p2,lx); |
p1[1] = (long) p2; e = bit_accuracy(lx)-1-expo(x); |
e = bit_accuracy(lx)-1-expo(x); |
if (e<0) err(talker,"integral part not significant in scfont"); |
if (e < 0) err(talker,"integral part not significant in scfont"); |
|
p1 = cgetg(3, t_FRACN); |
|
p1[1] = (long)p2; |
|
p1[2] = lshifti(gun, e); |
|
|
|
p3 = cgetg(3, t_FRACN); |
|
p3[1] = laddsi(signe(x), p2); |
|
p3[2] = p1[2]; |
|
p1 = Qsfcont(p1,NULL,k); |
|
return gerepilecopy(av, Qsfcont(p3,p1,k)); |
|
|
l1 = (e>>TWOPOTBITS_IN_LONG)+3; p2=cgeti(l1); |
|
p2[1]= evalsigne(1)|evallgefint(l1); |
|
p2[2]= (1L<<(e&(BITS_IN_LONG-1))); |
|
for (i=3; i<l1; i++) p2[i]=0; |
|
p1[2]=(long) p2; p3=cgetg(3,t_FRACN); |
|
p3[2]=lcopy(p2); |
|
p3[1]=laddsi(signe(x),(GEN)p1[1]); |
|
p1=sfcont(p1,p1,k); tetpil=avma; |
|
return gerepile(l,tetpil,sfcont(p3,p1,k)); |
|
|
|
case t_FRAC: case t_FRACN: |
case t_FRAC: case t_FRACN: |
l=avma; lx1=lgefint(x[2]); |
av = avma; |
l1 = (long) ((double) BYTES_IN_LONG/4.0*46.093443*(lx1-2)+3); |
return gerepileupto(av, Qsfcont(x, NULL, k)); |
if (k>0 && ++k > 0 && l1 > k) l1 = k; /* beware overflow */ |
|
if ((ulong)l1>LGBITS) l1=LGBITS; |
|
if (lgefint(x[1]) >= lx1) p1=icopy((GEN)x[1]); |
|
else affii((GEN)x[1],p1=cgeti(lx1)); |
|
p2=icopy((GEN)x[2]); lx1=lg(x1); |
|
y=cgetg(l1,t_VEC); f=(x!=x1); i=0; |
|
while (!gcmp0(p2) && i<=l1-2) |
|
{ |
|
i++; y[i]=ldvmdii(p1,p2,&p3); |
|
if (signe(p3)>=0) affii(p3,p1); |
|
else |
|
{ |
|
p4=addii(p3,p2); affii(p4,p1); cgiv(p4); |
|
y[1]=laddsi(-1,(GEN)y[1]); |
|
} |
|
cgiv(p3); p4=p1; p1=p2; p2=p4; |
|
if (f) |
|
{ |
|
if (i>=lx1) { i--; break; } |
|
if (!egalii((GEN)y[i],(GEN)x1[i])) |
|
{ |
|
av=avma; |
|
if (gcmp1(absi(subii((GEN)x1[i],(GEN)y[i])))) |
|
{ |
|
if (i>=lx1-1 || !gcmp1((GEN)x1[i+1])) |
|
affii((GEN)x1[i],(GEN)y[i]); |
|
} |
|
else i--; |
|
avma=av; break; |
|
} |
|
} |
|
} |
|
if (i>=2 && gcmp1((GEN)y[i])) |
|
{ |
|
cgiv((GEN)y[i]); --i; cgiv((GEN)y[i]); |
|
y[i]=laddsi(1,(GEN)y[i]); |
|
} |
|
setlg(y,i+1); return gerepileupto(l, y); |
|
} |
} |
err(typeer,"sfcont"); |
err(typeer,"sfcont"); |
} |
} |
|
|
switch(tx) |
switch(tx) |
{ |
{ |
case t_POL: |
case t_POL: return _veccopy(x); |
y=cgetg(2,t_VEC); y[1]=lcopy(x); break; |
|
case t_SER: |
case t_SER: |
av=avma; p1=gtrunc(x); tetpil=avma; |
av = avma; p1 = gtrunc(x); |
y=gerepile(av,tetpil,sfcont(p1,p1,k)); break; |
return gerepileupto(av,sfcont(p1,k)); |
case t_RFRAC: |
case t_RFRAC: |
case t_RFRACN: |
case t_RFRACN: |
av=avma; l1=lgef(x[1]); if (lgef(x[2])>l1) l1=lgef(x[2]); |
av = avma; |
if (k>0 && l1>k+1) l1=k+1; |
l = typ(x[1]) == t_POL? lgef(x[1]): 3; |
yp=cgetg(l1,t_VEC); p1=(GEN)x[1]; i=0; p2=(GEN)x[2]; |
if (lgef(x[2]) > l) l = lgef(x[2]); |
while (!gcmp0(p2) && i<=l1-2) |
if (k > 0 && l > k+1) l = k+1; |
|
y = cgetg(l,t_VEC); |
|
p1 = (GEN)x[1]; |
|
p2 = (GEN)x[2]; |
|
for (i=1; i<l; i++) |
{ |
{ |
i++; yp[i]=ldivres(p1,p2,&p3); |
y[i] = ldivres(p1,p2,&p3); |
p1=p2; p2=p3; |
if (gcmp0(p3)) { i++; break; } |
|
p1 = p2; p2 = p3; |
} |
} |
tetpil=avma; y=cgetg(i+1,t_VEC); |
setlg(y, i); return gerepilecopy(av, y); |
for (j=1; j<=i; j++) y[j]=lcopy((GEN)yp[j]); |
|
y=gerepile(av,tetpil,y); break; |
|
default: err(typeer,"sfcont"); |
|
return NULL; /* not reached */ |
|
} |
} |
return y; |
err(typeer,"sfcont"); |
|
return NULL; /* not reached */ |
} |
} |
|
|
static GEN |
static GEN |
sfcont2(GEN b, GEN x, long k) |
sfcont2(GEN b, GEN x, long k) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long l1 = lg(b), tx = typ(x), i; |
long l1 = lg(b), tx = typ(x), i; |
GEN y,p1; |
GEN y,p1; |
|
|
Line 1846 sfcont2(GEN b, GEN x, long k) |
|
Line 2118 sfcont2(GEN b, GEN x, long k) |
|
return gerepilecopy(av,y); |
return gerepilecopy(av,y); |
} |
} |
|
|
|
|
GEN |
GEN |
|
gcf(GEN x) |
|
{ |
|
return sfcont(x,0); |
|
} |
|
|
|
GEN |
|
gcf2(GEN b, GEN x) |
|
{ |
|
return contfrac0(x,b,0); |
|
} |
|
|
|
GEN |
|
gboundcf(GEN x, long k) |
|
{ |
|
return sfcont(x,k); |
|
} |
|
|
|
GEN |
|
contfrac0(GEN x, GEN b, long flag) |
|
{ |
|
long lb,tb,i; |
|
GEN y; |
|
|
|
if (!b || gcmp0(b)) return sfcont(x,flag); |
|
tb = typ(b); |
|
if (tb == t_INT) return sfcont(x,itos(b)); |
|
if (! is_matvec_t(tb)) err(typeer,"contfrac0"); |
|
|
|
lb=lg(b); if (lb==1) return cgetg(1,t_VEC); |
|
if (tb != t_MAT) return sfcont2(b,x,flag); |
|
if (lg(b[1])==1) return sfcont(x,flag); |
|
y = (GEN) gpmalloc(lb*sizeof(long)); |
|
for (i=1; i<lb; i++) y[i]=coeff(b,1,i); |
|
x = sfcont2(y,x,flag); free(y); return x; |
|
} |
|
|
|
GEN |
pnqn(GEN x) |
pnqn(GEN x) |
{ |
{ |
long av=avma,tetpil,lx,ly,tx=typ(x),i; |
gpmem_t av=avma,tetpil; |
|
long lx,ly,tx=typ(x),i; |
GEN y,p0,p1,q0,q1,a,b,p2,q2; |
GEN y,p0,p1,q0,q1,a,b,p2,q2; |
|
|
if (! is_matvec_t(tx)) err(typeer,"pnqn"); |
if (! is_matvec_t(tx)) err(typeer,"pnqn"); |
Line 1899 bestappr_mod(GEN x, GEN A, GEN B) |
|
Line 2210 bestappr_mod(GEN x, GEN A, GEN B) |
|
{ |
{ |
case t_INTMOD: |
case t_INTMOD: |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN a,b,d, t = cgetg(3, t_FRAC); |
GEN a,b,d, t = cgetg(3, t_FRAC); |
if (! ratlift((GEN)x[2], (GEN)x[1], &a,&b,A,B)) return NULL; |
if (! ratlift((GEN)x[2], (GEN)x[1], &a,&b,A,B)) return NULL; |
if (is_pm1(b)) return icopy_av(a, (GEN)av); |
if (is_pm1(b)) return icopy_av(a, (GEN)av); |
Line 1928 bestappr_mod(GEN x, GEN A, GEN B) |
|
Line 2239 bestappr_mod(GEN x, GEN A, GEN B) |
|
GEN |
GEN |
bestappr(GEN x, GEN k) |
bestappr(GEN x, GEN k) |
{ |
{ |
ulong av = avma, tetpil; |
gpmem_t av = avma, tetpil; |
long tx,tk=typ(k),lx,i,e; |
long tx,tk=typ(k),lx,i,e; |
GEN p0,p1,p,q0,q1,q,a,y; |
GEN p0,p1,p,q0,q1,q,a,y; |
|
|
Line 1939 bestappr(GEN x, GEN k) |
|
Line 2250 bestappr(GEN x, GEN k) |
|
k = gcvtoi(k,&e); |
k = gcvtoi(k,&e); |
} |
} |
if (signe(k) <= 0) k=gun; |
if (signe(k) <= 0) k=gun; |
tx=typ(x); if (tx == t_FRACN) x = gred(x); |
tx=typ(x); if (tx == t_FRACN) { x = gred(x); tx = typ(x); } |
switch(tx) |
switch(tx) |
{ |
{ |
case t_INT: |
case t_INT: |
Line 1954 bestappr(GEN x, GEN k) |
|
Line 2265 bestappr(GEN x, GEN k) |
|
} |
} |
|
|
case t_REAL: |
case t_REAL: |
p1=gun; p0=gfloor(x); q1=gzero; q0=gun; a=p0; |
p1=gun; a=p0=gfloor(x); q1=gzero; q0=gun; |
while (gcmp(q0,k)<=0) |
while (cmpii(q0,k)<=0) |
{ |
{ |
x = gsub(x,a); |
x = gsub(x,a); |
if (gcmp0(x)) { p1=p0; q1=q0; break; } |
if (gcmp0(x)) { p1=p0; q1=q0; break; } |
Line 1980 bestappr(GEN x, GEN k) |
|
Line 2291 bestappr(GEN x, GEN k) |
|
GEN |
GEN |
bestappr0(GEN x, GEN a, GEN b) |
bestappr0(GEN x, GEN a, GEN b) |
{ |
{ |
ulong av; |
gpmem_t av; |
GEN t; |
GEN t; |
if (!b) return bestappr(x,a); |
if (!b) return bestappr(x,a); |
av = avma; |
av = avma; |
Line 2027 update_f(GEN f, GEN a) |
|
Line 2338 update_f(GEN f, GEN a) |
|
GEN |
GEN |
fundunit(GEN x) |
fundunit(GEN x) |
{ |
{ |
long av = avma, av2,tetpil,lim,r,flp,flq; |
gpmem_t av = avma, av2, lim; |
|
long r,flp,flq; |
GEN pol,y,a,u,v,u1,v1,sqd,f; |
GEN pol,y,a,u,v,u1,v1,sqd,f; |
|
|
if (typ(x) != t_INT) err(arither1); |
if (typ(x) != t_INT) err(arither1); |
Line 2056 fundunit(GEN x) |
|
Line 2368 fundunit(GEN x) |
|
pol = quadpoly(x); y = get_quad(f,pol,r); |
pol = quadpoly(x); y = get_quad(f,pol,r); |
if (!flq) v1 = y; else { update_f(f,a); v1 = get_quad(f,pol,r); } |
if (!flq) v1 = y; else { update_f(f,a); v1 = get_quad(f,pol,r); } |
|
|
u1=gconj(y); tetpil=avma; y=gdiv(v1,u1); |
u1=gconj(y); y=gdiv(v1,u1); |
if (signe(y[3]) < 0) { tetpil=avma; y=gneg(y); } |
if (signe(y[3]) < 0) y = gneg(y); |
return gerepile(av,tetpil,y); |
return gerepileupto(av, y); |
} |
} |
|
|
GEN |
GEN |
Line 2070 gregula(GEN x, long prec) |
|
Line 2382 gregula(GEN x, long prec) |
|
GEN |
GEN |
regula(GEN x, long prec) |
regula(GEN x, long prec) |
{ |
{ |
long av = avma,av2,lim,r,fl,rexp; |
gpmem_t av = avma,av2,lim; |
|
long r,fl,rexp; |
GEN reg,rsqd,y,u,v,u1,v1, sqd = racine(x); |
GEN reg,rsqd,y,u,v,u1,v1, sqd = racine(x); |
|
|
if (signe(x)<=0) err(arither2); |
if (signe(x)<=0) err(arither2); |
Line 2079 regula(GEN x, long prec) |
|
Line 2392 regula(GEN x, long prec) |
|
rsqd = gsqrt(x,prec); |
rsqd = gsqrt(x,prec); |
if (egalii(sqri(sqd),x)) err(talker,"square argument in regula"); |
if (egalii(sqri(sqd),x)) err(talker,"square argument in regula"); |
|
|
rexp=0; reg=cgetr(prec); affsr(2,reg); |
rexp=0; reg = stor(2, prec); |
av2=avma; lim = stack_lim(av2,2); |
av2=avma; lim = stack_lim(av2,2); |
u = stoi(r); v = gdeux; |
u = stoi(r); v = gdeux; |
for(;;) |
for(;;) |
Line 2103 regula(GEN x, long prec) |
|
Line 2416 regula(GEN x, long prec) |
|
y = mplog(divri(reg,v)); |
y = mplog(divri(reg,v)); |
if (rexp) |
if (rexp) |
{ |
{ |
u1 = mulsr(rexp, glog(gdeux, prec)); |
u1 = mulsr(rexp, mplog2(prec)); |
setexpo(u1, expo(u1)+1); |
setexpo(u1, expo(u1)+1); |
y = addrr(y,u1); |
y = addrr(y,u1); |
} |
} |
Line 2174 end_classno(GEN h, GEN hin, GEN forms, long lform) |
|
Line 2487 end_classno(GEN h, GEN hin, GEN forms, long lform) |
|
q = ground(gdiv(hin,h)); /* approximate order of G/H */ |
q = ground(gdiv(hin,h)); /* approximate order of G/H */ |
for (i=1; i < lform; i++) |
for (i=1; i < lform; i++) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
fg = powgi((GEN)forms[i], h); |
fg = powgi((GEN)forms[i], h); |
fh = powgi(fg, q); |
fh = powgi(fg, q); |
a = (GEN)fh[1]; |
a = (GEN)fh[1]; |
Line 2198 to_form(GEN x, long f) |
|
Line 2511 to_form(GEN x, long f) |
|
return redimag(primeform(x, stoi(f), 0)); |
return redimag(primeform(x, stoi(f), 0)); |
} |
} |
|
|
|
static GEN |
|
conductor_part(GEN x, GEN *ptD, GEN *ptreg, GEN *ptfa) |
|
{ |
|
long n,i,k,s=signe(x),fl2; |
|
GEN e,p,H,d,D,fa,reg; |
|
|
|
fa = auxdecomp(absi(x),1); |
|
e = gtovecsmall((GEN)fa[2]); |
|
fa = (GEN)fa[1]; |
|
n = lg(fa); d = gun; |
|
for (i=1; i<n; i++) |
|
if (e[i] & 1) d = mulii(d,(GEN)fa[i]); |
|
if (mod4(d) == 2-s) fl2 = 0; |
|
else |
|
{ |
|
fl2 = (mod4(x)==0); |
|
if (!fl2) err(funder2,"classno"); |
|
d = shifti(d,2); |
|
} |
|
H = gun; D = (s<0)? negi(d): d; /* d = abs(D) */ |
|
/* f \prod_{p|f} [ 1 - (D/p) p^-1 ] */ |
|
for (i=1; i<n; i++) |
|
{ |
|
p = (GEN)fa[i]; |
|
k = e[i]; |
|
if (fl2 && i==1) k -= 2; /* p = 2 */ |
|
if (k >= 2) |
|
{ |
|
H = mulii(H, subis(p, kronecker(D,p))); |
|
if (k>=4) H = mulii(H,gpowgs(p,(k>>1)-1)); |
|
} |
|
} |
|
|
|
/* divide by [ O^* : O_K^* ] */ |
|
if (s < 0) |
|
{ |
|
reg = NULL; |
|
if (!is_bigint(d)) |
|
switch(itos(d)) |
|
{ |
|
case 4: H = divis(H,2); break; |
|
case 3: H = divis(H,3); break; |
|
} |
|
} else { |
|
reg = regula(D,DEFAULTPREC); |
|
if (!egalii(x,D)) |
|
H = divii(H, ground(gdiv(regula(x,DEFAULTPREC), reg))); |
|
} |
|
*ptfa = fa; *ptD = D; if (ptreg) *ptreg = reg; |
|
return H; |
|
} |
|
|
|
|
static long |
static long |
two_rank(GEN x) |
two_rank(GEN x) |
{ |
{ |
Line 2215 two_rank(GEN x) |
|
Line 2581 two_rank(GEN x) |
|
} |
} |
|
|
#define MAXFORM 11 |
#define MAXFORM 11 |
#define _low(x) (__x=(GEN)x, modBIL(__x)) |
#if 0 /* gcc-ism */ |
|
# define _low(x) ({GEN __x=(GEN)x; modBIL(__x);}) |
|
#else |
|
# define _low(x) (__x = (GEN)(x), modBIL(__x)) |
|
#endif |
|
|
/* h(x) for x<0 using Baby Step/Giant Step. |
/* h(x) for x<0 using Baby Step/Giant Step. |
* Assumes G is not too far from being cyclic. |
* Assumes G is not too far from being cyclic. |
Line 2224 two_rank(GEN x) |
|
Line 2594 two_rank(GEN x) |
|
GEN |
GEN |
classno(GEN x) |
classno(GEN x) |
{ |
{ |
long av = avma, av2,c,lforms,k,l,i,j,j1,j2,lim,com,s, forms[MAXFORM]; |
gpmem_t av = avma, av2, lim; |
long r2; |
long r2,c,lforms,k,l,i,j,j1,j2,com,s, forms[MAXFORM]; |
GEN a,b,count,index,tabla,tablb,hash,p1,p2,hin,h,f,fh,fg,ftest; |
GEN a,b,count,index,tabla,tablb,hash,p1,p2,hin,h,f,fh,fg,ftest; |
GEN __x; |
GEN Hf,D,fa; |
byteptr p = diffptr; |
byteptr p = diffptr; |
|
GEN __x; |
|
|
if (typ(x) != t_INT) err(arither1); |
if (typ(x) != t_INT) err(arither1); |
s=signe(x); if (s>=0) return classno2(x); |
s=signe(x); if (s>=0) return classno2(x); |
|
|
k=mod4(x); if (k==1 || k==2) err(funder2,"classno"); |
k=mod4(x); if (k==1 || k==2) err(funder2,"classno"); |
if (gcmpgs(x,-12) >= 0) return gun; |
if (cmpis(x,-12) >= 0) return gun; |
|
|
p2 = gsqrt(absi(x),DEFAULTPREC); |
Hf = conductor_part(x,&D,NULL,&fa); |
|
if (cmpis(D,-12) >= 0) return gerepilecopy(av, Hf); |
|
|
|
p2 = gsqrt(absi(D),DEFAULTPREC); |
p1 = divrr(p2,mppi(DEFAULTPREC)); |
p1 = divrr(p2,mppi(DEFAULTPREC)); |
p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1)); |
p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1)); |
s = 1000; |
s = 1000; |
|
|
if (is_bigint(p2)) err(talker,"discriminant too big in classno"); |
if (is_bigint(p2)) err(talker,"discriminant too big in classno"); |
s = itos(p2); if (s < 1000) s = 1000; |
s = itos(p2); if (s < 1000) s = 1000; |
} |
} |
r2 = two_rank(x); |
r2 = two_rank(D); |
c=lforms=0; |
c=lforms=0; |
while (c <= s && *p) |
while (c <= s && *p) |
{ |
{ |
c += *p++; k = krogs(x,c); |
NEXT_PRIME_VIADIFF(c,p); |
|
k = krogs(D,c); |
if (!k) continue; |
if (!k) continue; |
|
|
av2 = avma; |
av2 = avma; |
|
|
tabla = new_chunk(10000); |
tabla = new_chunk(10000); |
tablb = new_chunk(10000); |
tablb = new_chunk(10000); |
hash = new_chunk(10000); |
hash = new_chunk(10000); |
f = gsqr(to_form(x, forms[0])); |
f = gsqr(to_form(D, forms[0])); |
p1 = fh = powgi(f, h); |
p1 = fh = powgi(f, h); |
for (i=0; i<s; i++) |
for (i=0; i<s; i++) |
{ |
{ |
|
|
h = addii(addis(h,j2), mulss(s,com)); |
h = addii(addis(h,j2), mulss(s,com)); |
forms[0] = (long)f; |
forms[0] = (long)f; |
for (i=1; i<lforms; i++) |
for (i=1; i<lforms; i++) |
forms[i] = (long)gsqr(to_form(x, forms[i])); |
forms[i] = (long)gsqr(to_form(D, forms[i])); |
h = end_classno(h,hin,forms,lforms); |
h = end_classno(h,hin,forms,lforms); |
|
h = mulii(h,Hf); |
return gerepileuptoint(av, shifti(h, r2)); |
return gerepileuptoint(av, shifti(h, r2)); |
} |
} |
} |
} |
|
|
GEN |
GEN |
classno2(GEN x) |
classno2(GEN x) |
{ |
{ |
long av=avma,tetpil,n,i,k,s=signe(x),fl2; |
gpmem_t av = avma; |
GEN p1,p2,p3,p4,p5,p7,p8,pi4,reg,logd,d,fd; |
long n,i,k,s = signe(x); |
|
GEN p1,p2,p3,p4,p5,p7,Hf,Pi,reg,logd,d,D; |
|
|
if (typ(x) != t_INT) err(arither1); |
if (typ(x) != t_INT) err(arither1); |
if (!s) err(arither2); |
if (!s) err(arither2); |
if (s<0 && gcmpgs(x,-12) >= 0) return gun; |
if (s < 0 && cmpis(x,-12) >= 0) return gun; |
|
|
p1=auxdecomp(absi(x),1); p2=(GEN)p1[2]; p1=(GEN)p1[1]; |
Hf = conductor_part(x, &D, ®, &p1); |
n=lg(p1); d=gun; |
if (s < 0 && cmpis(D,-12) >= 0) |
for (i=1; i<n; i++) |
return gerepilecopy(av, Hf); /* |D| < 12*/ |
if (mod2((GEN)p2[i])) d=mulii(d,(GEN)p1[i]); |
|
fl2 = (mod4(x)==0); /* 4 | x */ |
|
if (mod4(d) != 2-s) |
|
{ |
|
if (!fl2) err(funder2,"classno2"); |
|
d = shifti(d,2); |
|
} |
|
else fl2=0; |
|
p8=gun; fd = (s<0)? negi(d): d; /* d = abs(fd) */ |
|
for (i=1; i<n; i++) |
|
{ |
|
k=itos((GEN)p2[i]); p4=(GEN)p1[i]; |
|
if (fl2 && i==1) k -= 2; /* p4 = 2 */ |
|
if (k>=2) |
|
{ |
|
p8=mulii(p8,subis(p4,kronecker(fd,p4))); |
|
if (k>=4) p8=mulii(p8,gpuigs(p4,(k>>1)-1)); |
|
} |
|
} |
|
if (s<0 && lgefint(d)==3) |
|
{ |
|
switch(d[2]) |
|
{ |
|
case 4: p8=gdivgs(p8,2); break; |
|
case 3: p8=gdivgs(p8,3); break; |
|
} |
|
if (d[2] < 12) /* |fd| < 12*/ |
|
{ tetpil=avma; return gerepile(av,tetpil,icopy(p8)); } |
|
} |
|
|
|
pi4 = mppi(DEFAULTPREC); logd = glog(d,DEFAULTPREC); |
Pi = mppi(DEFAULTPREC); |
p1 = gsqrt(gdiv(gmul(d,logd),gmul2n(pi4,1)),DEFAULTPREC); |
d = absi(D); logd = glog(d,DEFAULTPREC); |
p4 = divri(pi4,d); p7 = ginv(mpsqrt(pi4)); |
p1 = mpsqrt(gdiv(gmul(d,logd), gmul2n(Pi,1))); |
if (s > 0) |
if (s > 0) |
{ |
{ |
reg = regula(d,DEFAULTPREC); |
p2 = subsr(1, gmul2n(divrr(mplog(reg),logd),1)); |
p2 = gsubsg(1, gmul2n(gdiv(glog(reg,DEFAULTPREC),logd),1)); |
if (gcmp(gsqr(p2), divsr(2,logd)) >= 0) p1 = gmul(p2,p1); |
p3 = gsqrt(gdivsg(2,logd),DEFAULTPREC); |
|
if (gcmp(p2,p3)>=0) p1 = gmul(p2,p1); |
|
} |
} |
else reg = NULL; /* for gcc -Wall */ |
p1 = gtrunc(p1); |
p1 = gtrunc(p1); n=p1[2]; |
if (is_bigint(p1)) err(talker,"discriminant too large in classno"); |
if (lgefint(p1)!=3 || n<0) |
n = itos(p1); |
err(talker,"discriminant too large in classno"); |
|
|
|
|
p4 = divri(Pi,d); p7 = ginv(mpsqrt(Pi)); |
p1 = gsqrt(d,DEFAULTPREC); p3 = gzero; |
p1 = gsqrt(d,DEFAULTPREC); p3 = gzero; |
if (s > 0) |
if (s > 0) |
{ |
{ |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
{ |
{ |
k=krogs(fd,i); |
k = krogs(D,i); if (!k) continue; |
if (k) |
p2 = mulir(mulss(i,i),p4); |
{ |
p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC))); |
p2 = mulir(mulss(i,i),p4); |
p5 = addrr(divrs(mulrr(p1,p5),i), eint1(p2,DEFAULTPREC)); |
p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC))); |
p3 = (k>0)? addrr(p3,p5): subrr(p3,p5); |
p5 = addrr(divrs(mulrr(p1,p5),i),eint1(p2,DEFAULTPREC)); |
|
p3 = (k>0)? addrr(p3,p5): subrr(p3,p5); |
|
} |
|
} |
} |
p3 = shiftr(divrr(p3,reg),-1); |
p3 = shiftr(divrr(p3,reg),-1); |
if (!egalii(x,fd)) /* x != p3 */ |
|
p8 = gdiv(p8,ground(gdiv(regula(x,DEFAULTPREC),reg))); |
|
} |
} |
else |
else |
{ |
{ |
p1 = gdiv(p1,pi4); |
p1 = gdiv(p1,Pi); |
for (i=1; i<=n; i++) |
for (i=1; i<=n; i++) |
{ |
{ |
k=krogs(fd,i); |
k = krogs(D,i); if (!k) continue; |
if (k) |
p2 = mulir(mulss(i,i),p4); |
{ |
p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC))); |
p2=mulir(mulss(i,i),p4); |
p5 = addrr(p5, divrr(divrs(p1,i),mpexp(p2))); |
p5=subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC))); |
p3 = (k>0)? addrr(p3,p5): subrr(p3,p5); |
p5=addrr(p5,divrr(divrs(p1,i),mpexp(p2))); |
|
p3 = (k>0)? addrr(p3,p5): subrr(p3,p5); |
|
} |
|
} |
} |
} |
} |
p3=ground(p3); tetpil=avma; |
return gerepileuptoint(av, mulii(Hf,ground(p3))); |
return gerepile(av,tetpil,mulii(p8,p3)); |
|
} |
} |
|
|
GEN |
GEN |