File: [local] / OpenXM_contrib / pari / src / basemath / Attic / arith1.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:30 2000 UTC (24 years, 5 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2
Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/*********************************************************************/
/** **/
/** ARITHMETIC FUNCTIONS **/
/** (first part) **/
/** **/
/*********************************************************************/
/* $Id: arith1.c,v 1.2 1999/09/23 17:50:55 karim Exp $ */
#include "pari.h"
/*********************************************************************/
/** **/
/** ARITHMETIC FUNCTION PROTOTYPES **/
/** **/
/*********************************************************************/
GEN
garith_proto(GEN f(GEN), GEN x, int do_error)
{
long tx=typ(x),lx,i;
GEN y;
if (is_matvec_t(tx))
{
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i] = (long) garith_proto(f,(GEN) x[i], do_error);
return y;
}
if (tx != t_INT && do_error) err(arither1);
return f(x);
}
GEN
arith_proto(long f(GEN), GEN x, int do_error)
{
long tx=typ(x),lx,i;
GEN y;
if (is_matvec_t(tx))
{
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i] = (long) arith_proto(f,(GEN) x[i], do_error);
return y;
}
if (tx != t_INT && do_error) err(arither1);
return stoi(f(x));
}
static GEN
arith_proto2(long f(GEN,GEN), GEN x, GEN n)
{
long l,i,tx = typ(x);
GEN y;
if (is_matvec_t(tx))
{
l=lg(x); y=cgetg(l,tx);
for (i=1; i<l; i++) y[i] = (long) arith_proto2(f,(GEN) x[i],n);
return y;
}
if (tx != t_INT) err(arither1);
tx=typ(n);
if (is_matvec_t(tx))
{
l=lg(n); y=cgetg(l,tx);
for (i=1; i<l; i++) y[i] = (long) arith_proto2(f,x,(GEN) n[i]);
return y;
}
if (tx != t_INT) err(arither1);
return stoi(f(x,n));
}
static GEN
arith_proto2gs(long f(GEN,long), GEN x, long y)
{
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) arith_proto2gs(f,(GEN) x[i],y);
return t;
}
if (tx != t_INT) err(arither1);
return stoi(f(x,y));
}
GEN
garith_proto2gs(GEN f(GEN,long), GEN x, long y)
{
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_proto2gs(f,(GEN) x[i],y);
return t;
}
if (tx != t_INT) err(arither1);
return f(x,y);
}
/*********************************************************************/
/** **/
/** BINARY DECOMPOSITION **/
/** **/
/*********************************************************************/
GEN
binaire(GEN x)
{
ulong m,u;
long i,lx,ex,ly,tx=typ(x);
GEN y,p1,p2;
switch(tx)
{
case t_INT:
lx=lgefint(x);
if (lx==2) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
ly = BITS_IN_LONG+1; m=HIGHBIT; u=x[2];
while (!(m & u)) { m>>=1; ly--; }
y = cgetg(ly+((lx-3)<<TWOPOTBITS_IN_LONG),t_VEC); ly=1;
do { y[ly] = m & u ? un : zero; ly++; } while (m>>=1);
for (i=3; i<lx; i++)
{
m=HIGHBIT; u=x[i];
do { y[ly] = m & u ? un : zero; ly++; } while (m>>=1);
}
break;
case t_REAL:
ex=expo(x);
if (!signe(x))
{
lx=1+max(-ex,0); y=cgetg(lx,t_VEC);
for (i=1; i<lx; i++) y[i]=zero;
return y;
}
lx=lg(x); y=cgetg(3,t_VEC);
if (ex > bit_accuracy(lx)) err(talker,"loss of precision in binary");
p1 = cgetg(max(ex,0)+2,t_VEC);
p2 = cgetg(bit_accuracy(lx)-ex,t_VEC);
y[1] = (long) p1; y[2] = (long) p2;
ly = -ex; ex++;
if (ex<=0)
{
p1[1]=zero; for (i=1; i <= -ex; i++) p2[i]=zero;
i=2; m=HIGHBIT;
}
else
{
ly=1;
for (i=2; i<lx && ly<=ex; i++)
{
m=HIGHBIT; u=x[i];
do
{ p1[ly] = (m & u) ? un : zero; ly++; }
while ((m>>=1) && ly<=ex);
}
ly=1;
if (m) i--; else m=HIGHBIT;
}
for (; i<lx; i++)
{
u=x[i];
do { p2[ly] = m & u ? un : zero; ly++; } while (m>>=1);
m=HIGHBIT;
}
break;
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=(long)binaire((GEN)x[i]);
break;
default: err(typeer,"binaire"); return NULL; /* not reached */
}
return y;
}
/* Renvoie 0 ou 1 selon que le bit numero n de x est a 0 ou 1 */
long
bittest(GEN x, long n)
{
long l;
if (!signe(x) || n<0) return 0;
l = lgefint(x)-1-(n>>TWOPOTBITS_IN_LONG);
if (l<=1) return 0;
n = (1L << (n & (BITS_IN_LONG-1))) & x[l];
return n? 1: 0;
}
static long
bittestg(GEN x, GEN n)
{
return bittest(x,itos(n));
}
GEN
gbittest(GEN x, GEN n)
{
return arith_proto2(bittestg,x,n);
}
/*********************************************************************/
/** **/
/** ORDER of INTEGERMOD x in (Z/nZ)* **/
/** **/
/*********************************************************************/
GEN
order(GEN x)
{
long av=avma,av1,i,e;
GEN o,m,p;
if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2])))
err(talker,"not an element of (Z/nZ)* in order");
o=phi((GEN) x[1]); m=decomp(o);
for (i = lg(m[1])-1; i; i--)
{
p=gcoeff(m,i,1); e=itos(gcoeff(m,i,2));
do
{
GEN o1=divii(o,p), y=powgi(x,o1);
if (!gcmp1((GEN)y[2])) break;
e--; o = o1;
}
while (e);
}
av1=avma; return gerepile(av,av1,icopy(o));
}
/******************************************************************/
/* */
/* GENERATOR of (Z/mZ)* */
/* */
/******************************************************************/
GEN
ggener(GEN m)
{
return garith_proto(gener,m,1);
}
GEN
gener(GEN m)
{
long av=avma,av1,k,i,e;
GEN x,t,q,p;
if (typ(m) != t_INT) err(arither1);
e = signe(m);
if (!e) err(talker,"zero modulus in znprimroot");
if (is_pm1(m)) { avma=av; return gmodulss(0,1); }
if (e<0) m = absi(m);
e = mod4(m);
if (e == 0) /* m = 0 mod 4 */
{
if (cmpis(m,4)) err(generer); /* m != 4, non cyclic */
return gmodulsg(3,m);
}
if (e == 2) /* m = 0 mod 2 */
{
q=shifti(m,-1); x = (GEN) gener(q)[2];
if (!mod2(x)) x = addii(x,q);
av1=avma; return gerepile(av,av1,gmodulcp(x,m));
}
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);
if (e>=2)
{
x = (GEN) gener(p)[2];
if (gcmp1(powmodulo(x,q, sqri(p)))) x = addii(x,p);
av1=avma; return gerepile(av,av1,gmodulcp(x,m));
}
p=(GEN)decomp(q)[1]; k=lg(p)-1; x=stoi(1);
for(;;)
{
x[2]++;
if (gcmp1(mppgcd(m,x)))
{
for (i=k; i; i--)
if (gcmp1(powmodulo(x, divii(q,(GEN)p[i]), m))) break;
if (!i) break;
}
}
av1=avma; return gerepile(av,av1,gmodulcp(x,m));
}
GEN
znstar(GEN n)
{
GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a;
long i,j,nbp,sizeh,av,tetpil;
if (typ(n) != t_INT) err(arither1);
if (!signe(n))
{
z=cgetg(4,t_VEC);
z[1]=deux; p1=cgetg(2,t_VEC);
z[2]=(long)p1; p1[1]=deux; p1=cgetg(2,t_VEC);
z[3]=(long)p1; p1[1]=lneg(gun);
return z;
}
av=avma; n=absi(n);
if (cmpis(n,2)<=0)
{
avma=av; z=cgetg(4,t_VEC);
z[1]=un;
z[2]=lgetg(1,t_VEC);
z[3]=lgetg(1,t_VEC);
return z;
}
list=factor(n); ep=(GEN)list[2]; list=(GEN)list[1]; nbp=lg(list)-1;
h = cgetg(nbp+2,t_VEC);
gen = cgetg(nbp+2,t_VEC);
moduli = cgetg(nbp+2,t_VEC);
switch(mod8(n))
{
case 0:
h[1] = lmul2n(gun,itos((GEN)ep[1])-2); h[2] = deux;
gen[1] = lstoi(5); gen[2] = laddis(gmul2n((GEN)h[1],1),-1);
moduli[1] = moduli[2] = lmul2n(gun,itos((GEN)ep[1]));
sizeh = nbp+1; i=3; j=2; break;
case 4:
h[1] = deux;
gen[1] = lstoi(3);
moduli[1] = lmul2n(gun,itos((GEN)ep[1]));
sizeh = nbp; i=j=2; break;
case 2: case 6:
sizeh = nbp-1; i=1; j=2; break;
case 1: case 3: case 5: case 7:
sizeh = nbp; i=j=1; break;
}
for ( ; j<=nbp; i++,j++)
{
p = (GEN)list[j]; q = gpuigs(p, itos((GEN)ep[j])-1);
h[i] = lmulii(addis(p,-1),q); p1 = mulii(p,q);
gen[i] = gener(p1)[2];
moduli[i] = (long)p1;
}
#if 0
if (nbp == 1 && is_pm1(ep[1]))
gen[1] = lmodulcp((GEN)gen[1],n);
else
#endif
for (i=1; i<=sizeh; i++)
{
q = (GEN)moduli[i]; a = (GEN)gen[i];
u = mpinvmod(q,divii(n,q));
gen[i] = lmodulcp(addii(a,mulii(mulii(subsi(1,a),u),q)), n);
}
for (i=sizeh; i>=2; i--)
for (j=i-1; j>=1; j--)
if (!divise((GEN)h[j],(GEN)h[i]))
{
d=bezout((GEN)h[i],(GEN)h[j],&u,&v);
q=divii((GEN)h[j],d);
h[j]=(long)mulii((GEN)h[i],q); h[i]=(long)d;
gen[j]=ldiv((GEN)gen[j], (GEN)gen[i]);
gen[i]=lmul((GEN)gen[i], powgi((GEN)gen[j], mulii(v,q)));
}
q=gun; for (i=1; i<=sizeh && !gcmp1((GEN)h[i]); i++) q=mulii(q,(GEN)h[i]);
setlg(h,i); setlg(gen,i); z=cgetg(4,t_VEC);
z[1]=(long)q; z[2]=(long)h; z[3]=(long)gen;
tetpil=avma; return gerepile(av,tetpil,gcopy(z));
}
/*********************************************************************/
/** **/
/** INTEGRAL SQUARE ROOT **/
/** **/
/*********************************************************************/
GEN
gracine(GEN a)
{
return garith_proto(racine,a,1); /* hm. --GN */
}
GEN
racine(GEN a)
{
GEN x,y,z;
long av,tetpil,k;
if (typ(a) != t_INT) err(arither1);
switch (signe(a))
{
case 0: return gzero;
case -1:
x=cgetg(3,t_COMPLEX); x[1]=zero;
setsigne(a,1); x[2]= (long) racine(a); setsigne(a,-1);
return x;
case 1:
av=avma; k = (long) sqrt((double)(ulong)a[2]);
x = shifti(stoi(k+1), (lgefint(a)-3)*(BITS_IN_LONG/2));
do
{
z = shifti(addii(x,divii(a,x)), -1);
y = x; x = z;
}
while (cmpii(x,y) < 0);
tetpil=avma; return gerepile(av,tetpil,icopy(y));
}
return NULL; /* not reached */
}
/*********************************************************************/
/** **/
/** PERFECT SQUARE **/
/** **/
/*********************************************************************/
long
carrecomplet(GEN x, GEN *pt)
{
static int carresmod64[]={
1,1,0,0,1,0,0,0,0,1, 0,0,0,0,0,0,1,1,0,0,
0,0,0,0,0,1,0,0,0,0, 0,0,0,1,0,0,1,0,0,0,
0,1,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,1,0,0,
0,0,0,0};
static int carresmod63[]={
1,1,0,0,1,0,0,1,0,1, 0,0,0,0,0,0,1,0,1,0,
0,0,1,0,0,1,0,0,1,0, 0,0,0,0,0,0,1,1,0,0,
0,0,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,1,0,
0,0,0};
static int carresmod65[]={
1,1,0,0,1,0,0,0,0,1, 1,0,0,0,1,0,1,0,0,0,
0,0,0,0,0,1,1,0,0,1, 1,0,0,0,0,1,1,0,0,1,
1,0,0,0,0,0,0,0,0,1, 0,1,0,0,0,1,1,0,0,0,
0,1,0,0,1};
static int carresmod11[]={1,1,0,1,1,1,0,0,0,1,0};
long av;
GEN y;
switch(signe(x))
{
case -1: return 0;
case 0: if (pt) *pt=gzero; return 1;
case 1: if (!carresmod64[mod64(x)]) return 0;
}
if (!carresmod63[smodis(x,63)]) return 0;
if (!carresmod65[smodis(x,65)]) return 0;
if (!carresmod11[smodis(x,11)]) return 0;
av=avma; y=racine(x);
if (!egalii(sqri(y),x)) { avma=av; return 0; }
if (pt) { *pt = y; avma=(long)y; } else avma=av;
return 1;
}
static GEN
polcarrecomplet(GEN x, GEN *pt)
{
long i,l,av,av2;
GEN y,a,b;
if (!signe(x)) return gun;
l=lgef(x); if ((l&1) == 0) return gzero; /* odd degree */
i=2; while (isexactzero((GEN)x[i])) i++;
if (i&1) return gzero;
av2 = avma; a = (GEN)x[i];
switch (typ(a))
{
case t_POL: case t_INT:
y = gcarrecomplet(a,&b); break;
default:
y = gcarreparfait(a); b = NULL; break;
}
if (y == gzero) { avma = av2; return gzero; }
av = avma; x = gdiv(x,a);
y = gtrunc(gsqrt(greffe(x,l,1),0)); av2 = avma;
if (!gegal(gsqr(y), x)) { avma=av; return gzero; }
if (pt)
{
avma = av2;
if (!gcmp1(a))
{
if (!b) b = gsqrt(a,DEFAULTPREC);
y = gmul(b,y);
}
*pt = gerepileupto(av,y);
}
else avma = av;
return gun;
}
GEN
gcarrecomplet(GEN x, GEN *pt)
{
long tx = typ(x),l,i;
GEN z,y,p,t;
if (!pt) return gcarreparfait(x);
if (is_matvec_t(tx))
{
l=lg(x); y=cgetg(l,tx); z=cgetg(l,tx);
for (i=1; i<l; i++)
{
t=gcarrecomplet((GEN)x[i],&p);
y[i] = (long)t;
z[i] = gcmp0(t)? zero: (long)p;
}
*pt=z; return y;
}
if (tx == t_POL) return polcarrecomplet(x,pt);
if (tx != t_INT) err(arither1);
return stoi(carrecomplet(x,pt));
}
GEN
gcarreparfait(GEN x)
{
GEN p1,p2,p3,p4;
long tx=typ(x),l,i,av,v;
switch(tx)
{
case t_INT:
return carreparfait(x)? gun: gzero;
case t_REAL:
return (signe(x)>=0)? gun: gzero;
case t_INTMOD:
if (!signe(x[2])) return gun;
av=avma; p1=factor(absi((GEN)x[1]));
p2=(GEN)p1[1]; l=lg(p2); p3=(GEN)p1[2];
for (i=1; i<l; i++)
{
v = pvaluation((GEN)x[2],(GEN)p2[i],&p4);
if (v < itos((GEN)p3[i]))
{
if (v&1) break;
if (!egalii((GEN)p2[i], gdeux))
{ if (kronecker(p4,(GEN)p2[i]) == -1) break; }
else
{
v = itos((GEN)p3[i]) - v;
if ((v>=3 && mod8(p4) != 1 ) ||
(v==2 && mod4(p4) != 1)) break;
}
}
}
avma=av; return i<l? gzero: gun;
case t_FRAC: case t_FRACN:
av=avma; l=carreparfait(mulii((GEN)x[1],(GEN)x[2]));
avma=av; return l? gun: gzero;
case t_COMPLEX:
return gun;
case t_PADIC:
p4=(GEN)x[4]; if (!signe(p4)) return gun;
if (valp(x)&1) return gzero;
if (cmpis((GEN)x[2],2))
return (kronecker((GEN)x[4],(GEN)x[2])== -1)? gzero: gun;
v=precp(x); /* here p=2 */
if ((v>=3 && mod8(p4) != 1 ) ||
(v==2 && mod4(p4) != 1)) return gzero;
return gun;
case t_POL:
return polcarrecomplet(x,NULL);
case t_SER:
if (!signe(x)) return gun;
if (valp(x)&1) return gzero;
return gcarreparfait((GEN)x[2]);
case t_RFRAC: case t_RFRACN:
av=avma;
l=itos(gcarreparfait(gmul((GEN)x[1],(GEN)x[2])));
avma=av; return stoi(l);
case t_QFR: case t_QFI:
return gcarreparfait((GEN)x[1]);
case t_VEC: case t_COL: case t_MAT:
l=lg(x); p1=cgetg(l,tx);
for (i=1; i<l; i++) p1[i]=(long)gcarreparfait((GEN)x[i]);
return p1;
}
err(impl,"issquare for this type");
return NULL; /* not reached */
}
/*********************************************************************/
/** **/
/** KRONECKER SYMBOL **/
/** **/
/*********************************************************************/
#define ome(t) (labs(mod8(t)-4) == 1)
GEN
gkronecker(GEN x, GEN y)
{
return arith_proto2(kronecker,x,y);
}
long
kronecker(GEN x, GEN y)
{
GEN z;
long av,r,s=1;
av=avma;
switch (signe(y))
{
case -1: y=negi(y); if (signe(x)<0) s= -1; break;
case 0: return is_pm1(x);
}
r=vali(y);
if (r)
{
if (mpodd(x))
{
if (odd(r) && ome(x)) s= -s;
y=shifti(y,-r);
}
else { avma=av; return 0; }
}
x=modii(x,y);
while (signe(x))
{
r=vali(x);
if (r)
{
if (odd(r) && ome(y)) s= -s;
x=shifti(x,-r);
}
/* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
if (y[lgefint(y)-1]&2 && x[lgefint(x)-1]&2) s = -s;
z=resii(y,x); y=x; x=z;
}
avma=av; return is_pm1(y)? s: 0;
}
GEN
gkrogs(GEN x, long y)
{
return arith_proto2gs(krogs,x,y);
}
long
krogs(GEN x, long y)
{
long av,r,s=1,x1,z;
av=avma;
if (y<=0)
{
if (y) { y = -y; if (signe(x)<0) s = -1; }
else return is_pm1(x);
}
r=vals(y);
if (r)
{
if (mpodd(x))
{
if (odd(r) && ome(x)) s= -s;
y>>=r;
}
else return 0;
}
x1=smodis(x,y);
while (x1)
{
r=vals(x1);
if (r)
{
if (odd(r) && (labs((y&7)-4)==1)) s= -s;
x1>>=r;
}
if (y&2 && x1&2) s= -s;
z=y%x1; y=x1; x1=z;
}
avma=av; return (y==1)? s: 0;
}
long
krosg(long s, GEN x)
{
long av = avma, y = kronecker(stoi(s),x);
avma = av; return y;
}
long
kross(long x, long y)
{
long r,s=1,x1,z;
if (y<=0)
{
if (y) { y= -y; if (x<0) s = -1; }
else return (labs(x)==1);
}
r=vals(y);
if (r)
{
if (odd(x))
{
if (odd(r) && labs((x&7)-4) == 1) s = -s;
y>>=r;
}
else return 0;
}
x1=x%y; if (x1<0) x1+=y;
while (x1)
{
r=vals(x1);
if (r)
{
if (odd(r) && labs((y&7)-4) == 1) s= -s;
x1>>=r;
}
if (y&2 && x1&2) s= -s;
z=y%x1; y=x1; x1=z;
}
return (y==1)? s: 0;
}
long
hil0(GEN x, GEN y, GEN p)
{
return p? hil(x,y,p): hil(x,y,gzero);
}
#define eps(t) (((signe(t)*(t[lgefint(t)-1]))&3)==3)
long
hil(GEN x, GEN y, GEN p)
{
long a,b,av,tx,ty,z;
GEN p1,p2,u,v;
if (gcmp0(x) || gcmp0(y)) return 0;
av=avma; tx=typ(x); ty=typ(y);
if (tx>ty) { p1=x; x=y; y=p1; a=tx,tx=ty; ty=a; }
switch(tx)
{
case t_INT:
switch(ty)
{
case t_INT:
if (signe(p)<=0)
return (signe(x)<0 && signe(y)<0)? -1: 1;
a = odd(pvaluation(x,p,&u));
b = odd(pvaluation(y,p,&v));
if (egalii(p,gdeux))
{
z = (eps(u) && eps(v))? -1: 1;
if (a && ome(v)) z= -z;
if (b && ome(u)) z= -z;
}
else
{
z = (a && b && eps(p))? -1: 1;
if (a && kronecker(v,p)<0) z= -z;
if (b && kronecker(u,p)<0) z= -z;
}
avma=av; return z;
case t_REAL:
return (signe(x)<0 && signe(y)<0)? -1: 1;
case t_INTMOD:
if (egalii(gdeux,(GEN)y[1])) err(hiler1);
return hil(x,(GEN)y[2],(GEN)y[1]);
case t_FRAC: case t_FRACN:
p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p);
avma=av; return z;
case t_PADIC:
if (egalii(gdeux,(GEN)y[2]) && precp(y) <= 2) err(hiler1);
p1 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];
z=hil(x,p1,(GEN)y[2]); avma=av; return z;
}
break;
case t_REAL:
if (! is_frac_t(ty)) break;
if (signe(x) > 0) return 1;
return signe(y[1])*signe(y[2]);
case t_INTMOD:
if (egalii(gdeux,(GEN)y[1])) err(hiler1);
switch(ty)
{
case t_INTMOD:
if (!egalii((GEN)x[1],(GEN)y[1])) break;
return hil((GEN)x[2],(GEN)y[2],(GEN)x[1]);
case t_FRAC: case t_FRACN:
return hil((GEN)x[2],y,(GEN)x[1]);
case t_PADIC:
if (!egalii((GEN)x[1],(GEN)y[2])) break;
return hil((GEN)x[2],y,(GEN)x[1]);
}
break;
case t_FRAC: case t_FRACN:
p1=mulii((GEN)x[1],(GEN)x[2]);
switch(ty)
{
case t_FRAC: case t_FRACN:
p2=mulii((GEN)y[1],(GEN)y[2]);
z=hil(p1,p2,p); avma=av; return z;
case t_PADIC:
z=hil(p1,y,NULL); avma=av; return z;
}
break;
case t_PADIC:
if (ty != t_PADIC || !egalii((GEN)x[2],(GEN)y[2])) break;
p1 = odd(valp(x))? mulii((GEN)x[2],(GEN)x[4]): (GEN)x[4];
p2 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];
z=hil(p1,p2,(GEN)x[2]); avma=av; return z;
}
err(talker,"forbidden or incompatible types in hil");
return 0; /* not reached */
}
#undef eps
#undef ome
/*******************************************************************/
/* */
/* SQUARE ROOT MODULO p */
/* */
/*******************************************************************/
#define sqrmod(x,p) (resii(sqri(x),p))
/* Assume p is prime. return NULL if (a,p) = -1 */
GEN
mpsqrtmod(GEN a, GEN p)
{
long av = avma, av1,tetpil,lim,i,k,e;
GEN p1,p2,v,y,w,m1,m;
if (typ(a) != t_INT || typ(p) != t_INT) err(arither1);
p1=addsi(-1,p); e=vali(p1);
if (e == 0) /* p = 2 */
{
avma = av;
if (!signe(a) || !mod2(a)) return gzero;
return gun;
}
p2=shifti(p1,-e);
if (e == 1) y = p1;
else
for (k=2; ; k++)
{
av1 = avma;
m1 = m = powmodulo(stoi(k),p2,p);
for (i=1; i<e; i++)
if (gcmp1(m=sqrmod(m,p))) break;
if (i==e) { y = m1; break; }
avma = av1;
}
/* y contient un generateur de (Z/pZ)* eleve a la puis (p-1)/2^e */
p1=shifti(p2,-1); p1=powmodulo(a,p1,p);
if (!signe(p1)) { avma=av; return gzero; }
p2=mulii(p1,a); v = modii(p2,p);
p1=mulii(sqri(p1),a); w = modii(p1,p);
lim = stack_lim(av,1);
while (!gcmp1(w))
{
/* if p is not prime, next loop will not end */
p1=sqrmod(w,p); for (k=1; !gcmp1(p1); k++) p1 = sqrmod(p1,p);
if (k==e) { avma=av; return NULL; }
p1=y; for (i=1; i<e-k; i++) p1=sqrmod(p1,p);
e = k; v = modii(mulii(p1,v),p);
y = sqrmod(p1,p); w = modii(mulii(y,w),p);
if (low_stack(lim, stack_lim(av,1)))
{
GEN *gptr[3];
if(DEBUGMEM>1) err(warnmem,"mpsqrtmod");
gptr[0]=&y; gptr[1]=&v; gptr[2]=&w;
gerepilemany(av,gptr,3);
}
}
p1=subii(p,v); if (cmpii(v,p1)>0) v = p1;
tetpil=avma; return gerepile(av,tetpil,icopy(v));
}
/*********************************************************************/
/** **/
/** GCD & BEZOUT **/
/** **/
/*********************************************************************/
/* Ultra-fast private ulong gcd for trial divisions. Called with y odd;
x can be arbitrary (but will most of the time be smaller than y).
Will also be used from inside ifactor2.c, so it's `semi-private' really.
--GN */
/* Gotos are Harmful, and Programming is a Science. E.W.Dijkstra. */
ulong
ugcd(ulong x, ulong y) /* assume y&1==1, y > 1 */
{
if (!x) return y;
/* fix up x */
while (!(x&1)) x>>=1;
if (x==1) return 1;
if (x==y) return y;
else if (x>y) goto xislarger;/* will be rare, given how we'll use this */
/* loop invariants: x,y odd and distinct. */
yislarger:
if ((x^y)&2) /* ...01, ...11 or vice versa */
y=(x>>2)+(y>>2)+1; /* ==(x+y)>>2 except it can't overflow */
else /* ...01,...01 or ...11,...11 */
y=(y-x)>>2; /* now y!=0 in either case */
while (!(y&1)) y>>=1; /* kill any windfall-gained powers of 2 */
if (y==1) return 1; /* comparand == return value... */
if (x==y) return y; /* this and the next is just one comparison */
else if (x<y) goto yislarger;/* else fall through to xislarger */
xislarger: /* same as above, seen through a mirror */
if ((x^y)&2)
x=(x>>2)+(y>>2)+1;
else
x=(x-y)>>2; /* x!=0 */
while (!(x&1)) x>>=1;
if (x==1) return 1;
if (x==y) return y;
else if (x>y) goto xislarger;/* again, a decent [g]cc should notice that
it can re-use the comparison */
goto yislarger;
}
/* Gotos are useful, and Programming is an Art. D.E.Knuth. */
/* PS: Of course written with Dijkstra's lessons firmly in mind... --GN */
/* modified right shift binary algorithm with at most one division */
long
cgcd(long a,long b)
{
long v;
a=labs(a); if (!b) return a;
b=labs(b); if (!a) return b;
if (a>b) { a %= b; if (!a) return b; }
else { b %= a; if (!b) return a; }
v=vals(a|b); a>>=v; b>>=v;
if (a==1 || b==1) return 1L<<v;
if (b&1)
return ((long)ugcd((ulong)a, (ulong)b)) << v;
else
return ((long)ugcd((ulong)b, (ulong)a)) << v;
}
/* assume a>b>0, return gcd(a,b) as a GEN (for mppgcd) */
static GEN
mppgcd_cgcd(ulong a, ulong b)
{
GEN r = cgeti(3);
long v;
r[1] = evalsigne(1)|evallgefint(3);
a %= b; if (!a) { r[2]=(long)b; return r; }
v=vals(a|b); a>>=v; b>>=v;
if (a==1 || b==1) { r[2]=(long)(1UL<<v); return r; }
if (b&1)
r[2] = (long)(ugcd((ulong)a, (ulong)b) << v);
else
r[2] = (long)(ugcd((ulong)b, (ulong)a) << v);
return r;
}
void mppgcd_plus_minus(GEN x, GEN y, GEN t);
ulong mppgcd_resiu(GEN y, ulong x);
/* uses modified right-shift binary algorithm now --GN 1998Jul23 */
GEN
mppgcd(GEN a, GEN b)
{
long av,v,w;
GEN t, p1;
if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
switch (absi_cmp(a,b))
{
case 0: return absi(a);
case -1: t=b; b=a; a=t;
}
if (!signe(b)) return absi(a);
/* here |a|>|b|>0. Try single precision first */
if (lgefint(a)==3)
return mppgcd_cgcd((ulong)a[2], (ulong)b[2]);
if (lgefint(b)==3)
{
ulong u = mppgcd_resiu(a,(ulong)b[2]);
if (!u) return absi(b);
return mppgcd_cgcd((ulong)b[2], u);
}
/* larger than gcd: "avma=av" gerepile (erasing t) is valid */
av = avma; new_chunk(lgefint(b));
t = resii(a,b);
if (!signe(t)) { avma=av; return absi(b); }
a = b; b = t;
v = vali(a); a = shifti(a,-v); setsigne(a,1);
w = vali(b); b = shifti(b,-w); setsigne(b,1);
if (w < v) v = w;
switch(absi_cmp(a,b))
{
case 0: avma=av; a=shifti(a,v); return a;
case -1: p1=b; b=a; a=p1;
}
if (is_pm1(b)) { avma=av; return shifti(gun,v); }
/* we have three consecutive memory locations: a,b,t.
* All computations are done in place */
/* a and b are odd now, and a>b>1 */
while (lgefint(a) > 3)
{
/* if a=b mod 4 set t=a-b, otherwise t=a+b, then strip powers of 2 */
/* so that t <= (a+b)/4 < a/2 */
mppgcd_plus_minus(a,b, t);
if (is_pm1(t)) { avma=av; return shifti(gun,v); }
switch(absi_cmp(t,b))
{
case -1: p1 = a; a = b; b = t; t = p1; break;
case 1: p1 = a; a = t; t = p1; break;
case 0: avma = av; b=shifti(b,v); return b;
}
}
{
long r[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
r[2] = (long) ugcd((ulong)b[2], (ulong)a[2]);
avma = av; return shifti(r,v);
}
}
/* Extended bezout. Return d=pgcd(a,b) and &u,&v */
long
cbezout(long a,long b,long *uu,long *vv)
{
long av=avma,v,q,u,t,s, d = labs(a), r = labs(b);
GEN p1;
if (!b)
{
*vv=0;
if (!a) { *uu=1; return 1; }
if (a<0)
{ *uu=-1; return -a; }
else
{ *uu= 1; return a; }
}
u=1; t=0;
for(;;)
{
if (!r)
{
if (a<0) u = -u;
p1 = mulss(a,u); s = signe(p1);
if (!s) v = d / b;
else if (!is_bigint(p1))
{
ulong ul;
long sb = (b<0);
b = labs(b);
if (s < 0)
{
ul = (p1[2]+d);
ul /= b; v = sb? -(long)ul: (long)ul;
}
else
{
ul = p1[2]-d;
ul /= b; v = sb? (long)ul: -(long)ul;
}
}
else v = - itos(divis(addsi(-d,p1), b));
avma=av; *uu = u; *vv = v; return d;
}
v=d; d=r; q = v/r; r = v - q*r;
v=t; t = u - q*t; u=v;
}
}
GEN
bezout(GEN a, GEN b, GEN *ptu, GEN *ptv)
{
GEN u,v,v1,d,d1,q,r, *tmp;
long av;
if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
if (absi_cmp(a,b) < 0)
{
u=b; b=a; a=u;
tmp=ptu; ptu=ptv; ptv=tmp;
}
/* now |a| >= |b| */
if (!signe(b))
{
*ptv = gzero;
switch(signe(a))
{
case 0: *ptu = gun; return gzero;
case 1: *ptu = gun; return icopy(a);
case -1: *ptu = negi(gun); return negi(a);
}
}
if (!is_bigint(a))
{
long uu,vv, dd = cbezout(itos(a),itos(b),&uu,&vv);
*ptu = stoi(uu); *ptv = stoi(vv); return stoi(dd);
}
av = avma;
(void)new_chunk(lgefint(b) + (lgefint(a)<<1)); /* room for d, u and v */
d = a; d1 = b; v = gzero; v1 = gun;
for(;;)
{
q=dvmdii(d,d1,&r);
v=subii(v,mulii(q,v1));
u=v; v=v1; v1=u;
u=r; d=d1; d1=u; if (!signe(d1)) break;
}
u = divii(subii(d,mulii(b,v)),a);
avma = av; d = icopy(d); v = icopy(v); u = icopy(u);
if (signe(d) < 0)
{
setsigne(d,1);
setsigne(u,-signe(u));
setsigne(v,-signe(v));
}
*ptu = u; *ptv = v; return d;
}
/*********************************************************************/
/** **/
/** CHINESE REMAINDERS **/
/** **/
/*********************************************************************/
/* P.M. & M.H.
*
* Chinese Remainder Theorem. x and y must have the same type (integermod,
* polymod, or polynomial/vector/matrix recursively constructed with these
* as coefficients). Creates (with the same type) a z in the same residue
* class as x and the same residue class as y, if it is possible.
*
* We also allow (during recursion) two identical objects even if they are
* not integermod or polymod. For example, if
*
* x = [1. mod(5, 11), mod(X + mod(2, 7), X^2 + 1)]
* y = [1, mod(7, 17), mod(X + mod(0, 3), X^2 + 1)],
*
* then chinois(x, y) returns
*
* [1, mod(16, 187), mod(X + mod(9, 21), X^2 + 1)]
*
* Someone else may want to allow power series, complex numbers, and
* quadratic numbers.
*/
GEN
chinois(GEN x, GEN y)
{
long i,lx,vx,av,tetpil, tx = typ(x);
GEN z,p1,p2,d,u,v;
if (gegal(x,y)) return gcopy(x);
if (tx == typ(y)) switch(tx)
{
case t_POLMOD:
if (gegal((GEN)x[1],(GEN)y[1])) /* same modulus */
{
z=cgetg(3,tx);
z[1]=lcopy((GEN)x[1]);
z[2]=(long)chinois((GEN)x[2],(GEN)y[2]);
return z;
} /* fall through */
case t_INTMOD:
z=cgetg(3,tx); av=avma;
d=gbezout((GEN)x[1],(GEN)y[1],&u,&v);
if (!gegal(gmod((GEN)x[2],d), gmod((GEN)y[2],d))) break;
p1 = gdiv((GEN)x[1],d);
p2 = gadd((GEN)x[2], gmul(gmul(u,p1), gadd((GEN)y[2],gneg((GEN)x[2]))));
tetpil=avma; z[1]=lmul(p1,(GEN)y[1]); z[2]=lmod(p2,(GEN)z[1]);
gerepilemanyvec(av,tetpil,z+1,2); return z;
case t_POL:
lx=lgef(x); vx=varn(x); z=cgetg(lx,tx);
if (lx!=lgef(y) || vx!=varn(y)) break;
z[1]=evalsigne(1)|evallgef(lx)|evalvarn(vx);
for (i=2; i<lx; i++) z[i]=(long)chinois((GEN)x[i],(GEN)y[i]);
return z;
case t_VEC: case t_COL: case t_MAT:
lx=lg(x); z=cgetg(lx,tx); if (lx!=lg(y)) break;
for (i=1; i<lx; i++) z[i]=(long)chinois((GEN)x[i],(GEN)y[i]);
return z;
}
err(talker,"incompatible arguments in chinois");
return NULL; /* not reached */
}
/* return lift(chinois(x2 mod x1, y2 mod y1))
* assume(x1,y1)=1, xi,yi integers and z1 = x1*y1
*/
GEN
chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1)
{
long av = avma;
GEN ax,p1;
(void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1));
ax = mulii(mpinvmod(x1,y1), x1);
p1 = addii(x2, mulii(ax, subii(y2,x2)));
avma = av; return modii(p1,z1);
}
/*********************************************************************/
/** **/
/** INVERSE MODULO b **/
/** **/
/*********************************************************************/
GEN
mpinvmod(GEN a, GEN m)
{
GEN res;
if (invmod(a,m,&res)) return res;
err(talker,"impossible inverse modulo: %Z",gmodulcp(res,m));
return NULL; /* not reached */
}
/*********************************************************************/
/** **/
/** POWERING MODULO (a^n mod m) **/
/** **/
/*********************************************************************/
GEN
powmodulo(GEN a, GEN n, GEN m)
{
GEN y;
long av = avma,av1,lim,*p,man,k,nb;
GEN (*mul)(GEN,GEN) = mulii;
if (typ(a) != t_INT || typ(n) != t_INT || typ(m) != t_INT) err(arither1);
switch (signe(n))
{
case 0:
k = signe(resii(a,m)); avma=av;
return k? gun: gzero;
case -1:
a = mpinvmod(a,m); n = absi(n); /* fall through */
case 1:
y = modii(a,m);
if (!signe(y)) { avma=av; return gzero; }
if (lgefint(y)==3)
switch(y[2])
{
case 1: /* y = 1 */
avma=av; return gun;
case 2: /* y = 2, use shifti not mulii */
mul = (GEN (*)(GEN,GEN))shifti; a = (GEN)1;
}
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--)
{
y=resii(sqri(y), m);
if (man<0) y = modii(mul(y,a), m);
if (low_stack(lim, stack_lim(av1,1)))
{
if (DEBUGMEM>1) err(warnmem,"powmodulo");
y = gerepileupto(av1,y);
}
}
if (--nb == 0) break;
man = *++p, k = BITS_IN_LONG;
}
}
return gerepileupto(av,y);
}
/*********************************************************************/
/** **/
/** NEXT / PRECEDING (PSEUDO) PRIME **/
/** **/
/*********************************************************************/
GEN
gnextprime(GEN n)
{
return garith_proto(nextprime,n,0); /* accept non-integers */
}
GEN
gprecprime(GEN n)
{
return garith_proto(precprime,n,0); /* accept non-integers */
}
GEN
gisprime(GEN x)
{
return arith_proto(isprime,x,1);
}
long
isprime(GEN x)
{
return millerrabin(x,10);
}
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 **/
/** **/
/*********************************************************************/
/* 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
gisfundamental(GEN x)
{
return arith_proto(isfundamental,x,1);
}
long
isfundamental(GEN x)
{
long av,r;
GEN p1;
if (gcmp0(x)) return 0;
r=mod4(x);
if (!r)
{
av=avma; p1=shifti(x,-2);
r=mod4(p1); if (!r) return 0;
if (signe(x)<0) r=4-r;
r = (r==1)? 0: issquarefree(p1);
avma=av; return r;
}
if (signe(x)<0) r=4-r;
return (r==1) ? issquarefree(x) : 0;
}
GEN
quaddisc(GEN x)
{
long av=avma,tetpil=avma,i,r,tx=typ(x);
GEN p1,p2,f,s;
if (tx != t_INT && ! is_frac_t(tx)) err(arither1);
f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2];
s=gun;
for (i=1; i<lg(p1); i++)
if ( odd(mael(p2,i,2)) ) { tetpil=avma; s=gmul(s,(GEN)p1[i]); }
r=mod4(s); if (gsigne(x)<0) r=4-r;
if (r>1) { tetpil=avma; s=shifti(s,2); }
return gerepile(av,tetpil,s);
}
/*********************************************************************/
/** **/
/** FACTORIAL **/
/** **/
/*********************************************************************/
GEN
mpfact(long n)
{
long av,tetpil,lim,k;
GEN f;
if (n<2)
{
if (n<0) err(facter);
return gun;
}
av = avma; lim = stack_lim(av,1); f = gun;
for (k=2; k<n; k++)
if (low_stack(lim, stack_lim(av,1)))
{
if(DEBUGMEM>1) err(warnmem,"mpfact k=%ld",k);
tetpil = avma; f = gerepile(av,tetpil,mulsi(k,f));
}
else f = mulsi(k,f);
tetpil = avma; return gerepile(av,tetpil,mulsi(k,f));
}
GEN
mpfactr(long n, long prec)
{
long av,tetpil,lim,k;
GEN f = realun(prec);
if (n<2)
{
if (n<0) err(facter);
return f;
}
av = avma; lim = stack_lim(av,1);
for (k=2; k<n; k++)
if (low_stack(lim, stack_lim(av,1)))
{
if(DEBUGMEM>1) err(warnmem,"mpfactr k=%ld",k);
tetpil = avma; f = gerepile(av,tetpil,mulsr(k,f));
}
else f = mulsr(k,f);
tetpil = avma; return gerepile(av,tetpil,mulsr(k,f));
}
/*******************************************************************/
/** **/
/** LUCAS & FIBONACCI **/
/** **/
/*******************************************************************/
void
lucas(long n, GEN *ln, GEN *ln1)
{
long taille,av;
GEN z,t;
if (!n) { *ln=stoi(2); *ln1=stoi(1); return; }
taille=(long)(pariC3*(1+labs(n))+3);
*ln=cgeti(taille); *ln1=cgeti(taille);
av=avma; lucas(n/2,&z,&t);
switch(n % 4)
{
case -3:
addsiz(2,sqri(z),*ln1);
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:
subisz(sqri(z),2,*ln1);
subiiz(subis(mulii(z,t),1),*ln1,*ln); 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 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;
}
GEN
fibo(long n)
{
long av = avma;
GEN ln,ln1;
lucas(n-1,&ln,&ln1);
return gerepileupto(av, divis(addii(shifti(ln,1),ln1), 5));
}
/*******************************************************************/
/* */
/* CONTINUED FRACTIONS */
/* */
/*******************************************************************/
static GEN sfcont2(GEN b, GEN x, long k);
GEN
gcf(GEN x)
{
return sfcont(x,x,0);
}
GEN
gcf2(GEN b, GEN x)
{
return sfcont0(x,b,0);
}
GEN
gboundcf(GEN x, long k)
{
return sfcont(x,x,k);
}
GEN
sfcont0(GEN x, GEN b, long flag)
{
long lb,tb,i;
GEN y;
if (!b || gcmp0(b)) return sfcont(x,x,flag);
tb = typ(b);
if (tb == t_INT) return sfcont(x,x,itos(b));
if (! is_matvec_t(tb)) err(typeer,"sfcont0");
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,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
sfcont(GEN x, GEN x1, long k)
{
long av,lx=lg(x),tx=typ(x),e,i,j,l,l1,l2,lx1,tetpil,f;
GEN y,p1,p2,p3,p4,yp;
if (is_scalar_t(tx))
{
if (gcmp0(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
switch(tx)
{
case t_INT:
y=cgetg(2,t_VEC); y[1]=lcopy(x); return y;
case t_REAL:
l=avma; p1=cgetg(3,t_FRACN);
p2 = rcopy(x); settyp(p2,t_INT); setlgefint(p2,lx);
p1[1] = (long) p2; e = bit_accuracy(lx)-1-expo(x);
if (e<0) err(talker,"integral part not significant in scfont");
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:
l=avma; lx1=lgefint(x[2]);
l1 = (long) ((double) BYTES_IN_LONG/4.0*46.093443*(lx1-2)+3);
if (k>0 && l1>k+1) l1=k+1;
if (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]); l2=avma; 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); l2=l2-((l1-i-1)<<TWOPOTBYTES_IN_LONG);
return gerepile(l,l2,y);
}
err(typeer,"sfcont");
}
switch(tx)
{
case t_POL:
y=cgetg(2,t_VEC); y[1]=lcopy(x); break;
case t_SER:
av=avma; p1=gtrunc(x); tetpil=avma;
y=gerepile(av,tetpil,sfcont(p1,p1,k)); break;
case t_RFRAC:
case t_RFRACN:
av=avma; l1=lgef(x[1]); if (lgef(x[2])>l1) l1=lgef(x[2]);
if (k>0 && l1>k+1) l1=k+1;
yp=cgetg(l1,t_VEC); p1=(GEN)x[1]; i=0; p2=(GEN)x[2];
while (!gcmp0(p2) && i<=l1-2)
{
i++; yp[i]=ldivres(p1,p2,&p3);
p1=p2; p2=p3;
}
tetpil=avma; y=cgetg(i+1,t_VEC);
for (j=1; j<=i; j++) y[j]=lcopy((GEN)yp[j]);
y=gerepile(av,tetpil,y); break;
default: err(typeer,"sfcont");
}
return y;
}
static GEN
sfcont2(GEN b, GEN x, long k)
{
long l1 = lg(b), tx = typ(x), i,tetpil, av = avma;
GEN y,p1;
if (k)
{
if (k>=l1) err(typeer,"sfcont2");
l1 = k+1;
}
y=cgetg(l1,t_VEC);
if (l1==1) return y;
if (is_scalar_t(tx))
{
if (!is_intreal_t(tx) && !is_frac_t(tx)) err(typeer,"sfcont2");
}
else if (tx == t_SER) x = gtrunc(x);
if (!gcmp1((GEN)b[1])) x = gmul((GEN)b[1],x);
i = 2; y[1] = lfloor(x); p1 = gsub(x,(GEN)y[1]);
for ( ; i<l1 && !gcmp0(p1); i++)
{
x = gdiv((GEN)b[i],p1);
if (tx == t_REAL)
{
long e = expo(x);
if (e>0 && (e>>TWOPOTBITS_IN_LONG)+3 > lg(x)) break;
}
y[i] = lfloor(x);
p1 = gsub(x,(GEN)y[i]);
}
setlg(y,i); tetpil=avma;
return gerepile(av,tetpil,gcopy(y));
}
GEN
pnqn(GEN x)
{
long av=avma,tetpil,lx,ly,tx=typ(x),i;
GEN y,p0,p1,q0,q1,a,b,p2,q2;
if (! is_matvec_t(tx)) err(typeer,"pnqn");
lx=lg(x); if (lx==1) return idmat(2);
p0=gun; q0=gzero;
if (tx != t_MAT)
{
p1=(GEN)x[1]; q1=gun;
for (i=2; i<lx; i++)
{
a=(GEN)x[i];
p2=gadd(gmul(a,p1),p0); p0=p1; p1=p2;
q2=gadd(gmul(a,q1),q0); q0=q1; q1=q2;
}
}
else
{
ly=lg(x[1]);
if (ly==2)
{
p1=cgetg(lx,t_VEC); for (i=1; i<lx; i++) p1[i]=mael(x,i,1);
tetpil=avma; return gerepile(av,tetpil,pnqn(p1));
}
if (ly!=3) err(talker,"incorrect size in pnqn");
p1=gcoeff(x,2,1); q1=gcoeff(x,1,1);
for (i=2; i<lx; i++)
{
a=gcoeff(x,2,i); b=gcoeff(x,1,i);
p2=gadd(gmul(a,p1),gmul(b,p0)); p0=p1; p1=p2;
q2=gadd(gmul(a,q1),gmul(b,q0)); q0=q1; q1=q2;
}
}
tetpil=avma; y=cgetg(3,t_MAT);
p2=cgetg(3,t_COL); y[1]=(long)p2; p2[1]=lcopy(p1); p2[2]=lcopy(q1);
p2=cgetg(3,t_COL); y[2]=(long)p2; p2[1]=lcopy(p0); p2[2]=lcopy(q0);
return gerepile(av,tetpil,y);
}
GEN
bestappr(GEN x, GEN k)
{
long av=avma,tetpil,tx,tk=typ(k),lx,i,e;
GEN p0,p1,p,q0,q1,q,a,y;
if (tk != t_INT)
{
if (tk != t_REAL && !is_frac_t(tk))
err(talker,"incorrect bound type in bestappr");
k = gcvtoi(k,&e);
}
if (cmpis(k,1) < 0) k=gun;
tx=typ(x); if (tx == t_FRACN) x = gred(x);
switch(tx)
{
case t_INT:
if (av==avma) return icopy(x);
tetpil=avma; return gerepile(av,tetpil,icopy(x));
case t_FRAC:
if (cmpii((GEN)x[2],k) <= 0)
{
if (av==avma) return gcopy(x);
tetpil=avma; return gerepile(av,tetpil,gcopy(x));
}
case t_REAL:
p1=gun; p0=gfloor(x); q1=gzero; q0=gun; a=p0;
while (gcmp(q0,k)<=0)
{
x = gsub(x,a);
if (gcmp0(x)) { p1=p0; q1=q0; break; }
x = ginv(x); a = (gcmp(x,k)<0)? gfloor(x): k;
p = addii(mulii(a,p0),p1); p1=p0; p0=p;
q = addii(mulii(a,q0),q1); q1=q0; q0=q;
}
tetpil=avma; return gerepile(av,tetpil,gdiv(p1,q1));
case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC:
case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
lx = (tx==t_POL)? lgef(x): lg(x); y=cgetg(lx,tx);
for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
for ( ; i<lx; i++) y[i]=(long)bestappr((GEN)x[i],k);
return y;
}
err(typeer,"bestappr");
return NULL; /* not reached */
}
/***********************************************************************/
/** **/
/** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
/** **/
/***********************************************************************/
GEN
gfundunit(GEN x)
{
return garith_proto(fundunit,x,1);
}
static GEN
get_quad(GEN f, GEN pol, long r)
{
GEN y, c=(GEN)f[2], p1=(GEN)c[1], q1=(GEN)c[2];
y=cgetg(4,t_QUAD); y[1]=(long)pol;
y[2]=r? lsubii(p1,q1): (long)p1;
y[3]=(long)q1; return y;
}
/* replace f by f * [a,1; 1,0] */
static void
update_f(GEN f, GEN a)
{
GEN p1;
p1=gcoeff(f,1,1);
coeff(f,1,1)=laddii(mulii(a,p1), gcoeff(f,1,2));
coeff(f,1,2)=(long)p1;
p1=gcoeff(f,2,1);
coeff(f,2,1)=laddii(mulii(a,p1), gcoeff(f,2,2));
coeff(f,2,2)=(long)p1;
}
GEN
fundunit(GEN x)
{
long av = avma, av2,tetpil,lim,r,flp,flq;
GEN pol,y,a,u,v,u1,v1,sqd,f;
if (typ(x) != t_INT) err(arither1);
if (signe(x)<=0) err(arither2);
r=mod4(x); if (r==2 || r==3) err(funder2,"fundunit");
sqd=racine(x); av2=avma; lim=stack_lim(av2,2);
a = shifti(addsi(r,sqd),-1);
f=cgetg(3,t_MAT); f[1]=lgetg(3,t_COL); f[2]=lgetg(3,t_COL);
coeff(f,1,1)=(long)a; coeff(f,1,2)=un;
coeff(f,2,1)=un; coeff(f,2,2)=zero;
v = gdeux; u = stoi(r);
for(;;)
{
u1=subii(mulii(a,v),u); flp=egalii(u,u1); u=u1;
v1=divii(subii(x,sqri(u)),v); flq=egalii(v,v1); v=v1;
if (flq) break; a = divii(addii(sqd,u),v);
if (flp) break; update_f(f,a);
if (low_stack(lim, stack_lim(av2,2)))
{
GEN *gptr[4]; gptr[0]=&a; gptr[1]=&f; gptr[2]=&u; gptr[3]=&v;
if(DEBUGMEM>1) err(warnmem,"fundunit");
gerepilemany(av2,gptr,4);
}
}
pol = quadpoly(x); y = 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);
if (signe(y[3]) < 0) { tetpil=avma; y=gneg(y); }
return gerepile(av,tetpil,y);
}
GEN
gregula(GEN x, long prec)
{
return garith_proto2gs(regula,x,prec);
}
GEN
regula(GEN x, long prec)
{
long av = avma,av2,tetpil,lim,r,fl,rexp;
GEN ln2,reg,reg1,rsqd,y,u,v,a,u1,v1,sqd;
if (typ(x) != t_INT) err(arither1);
if (signe(x)<=0) err(arither2);
r=mod4(x); if (r==2 || r==3) err(funder2,"regula");
sqd=racine(x); rsqd=gsqrt(x,prec);
if (gegal(sqri(sqd),x)) err(talker,"square argument in regula");
rexp=0; reg=cgetr(prec); affsr(2,reg);
ln2 = mplog(reg);
av2=avma; lim = stack_lim(av2,2);
a = shifti(addsi(r,sqd),-1);
v = gdeux; u = stoi(r);
for(;;)
{
u1=subii(mulii(a,v),u);
reg1=divri(addir(u1,rsqd),v);
v1=divii(subii(x,sqri(u1)),v); fl=gegal(v,v1);
if (gegal(u,u1) || fl) break;
v=v1; u=u1; a = divii(addii(sqd,u),v);
reg = mulrr(reg,reg1); rexp += expo(reg); setexpo(reg,0);
if (rexp & ~EXPOBITS) err(muler4);
if (low_stack(lim, stack_lim(av2,2)))
{
GEN *gptr[4]; gptr[0]=&a; gptr[1]=® gptr[2]=&u; gptr[3]=&v;
if(DEBUGMEM>1) err(warnmem,"regula");
gerepilemany(av2,gptr,4);
}
}
reg = gsqr(reg); setexpo(reg,expo(reg)-1);
if (fl) reg = mulrr(reg, divri(addir(u1,rsqd),v));
y = mplog(divri(reg,v)); u1 = mulsr(rexp,ln2);
if (signe(u1)) setexpo(u1, expo(u1)+1);
tetpil=avma; return gerepile(av,tetpil,gadd(y,u1));
}
/*************************************************************************/
/** **/
/** CLASS NUMBER **/
/** **/
/*************************************************************************/
static GEN
gclassno(GEN x)
{
return garith_proto(classno,x,1);
}
static GEN
gclassno2(GEN x)
{
return garith_proto(classno2,x,1);
}
GEN
qfbclassno0(GEN x,long flag)
{
switch(flag)
{
case 0: return gclassno(x);
case 1: return gclassno2(x);
default: err(flagerr,"qfbclassno");
}
return NULL; /* not reached */
}
static GEN
end_classno(GEN h, GEN hin, GEN *formes, long ptforme)
{
long av,i,j,lim,com;
GEN p1,p2,p3,q,hp,fh,fg,ftest, f = formes[0];
if (ptforme>11) ptforme = 11;
p2=factor(h); p1=(GEN)p2[1]; p2=(GEN)p2[2];
for (i=1; i<lg(p1); i++)
{
lim = itos((GEN)p2[i]);
for (j=1; j<=lim; j++)
{
p3=divii(h,(GEN)p1[i]); fh=gpui(f,p3,0);
if (!is_pm1(fh[1])) break;
h = p3;
}
}
q=ground(gdiv(hin,h)); hp=mulii(q,h);
for (i=1; i<ptforme; i++)
{
fg=gpui(formes[i],h,0); fh=gpui(fg,q,0);
if (!is_pm1(fh[1]))
{
av = avma; ftest = fg;
for (com=1; ; com++)
{
if (gegal((GEN) fh[1], (GEN) ftest[1]) &&
absi_equal((GEN) fh[2], (GEN) ftest[2])) break;
ftest = gmul(ftest,fg);
}
if (signe(fh[2]) == signe(ftest[2])) com = -com;
avma = av; q = addsi(com,q);
hp = addii(hp, mulsi(com,h));
}
}
return hp;
}
/* calcul de h(x) pour x<0 par la methode de Shanks */
GEN
classno(GEN x)
{
long av = avma, av2,c,ptforme,k,i,j,j1,lim,com,j2,s;
GEN count,index,tabla,tablb,hash,court,p1,p2,hin,h,f,fh,fg,ftest,formes[100];
byteptr p = diffptr;
if (typ(x) != t_INT) err(arither1);
s=signe(x); if (s>=0) return classno2(x);
k=mod4(x); if (k==1 || k==2) err(funder2,"classno");
if (gcmpgs(x,-12) >= 0) return gun;
p2 = gsqrt(absi(x),DEFAULTPREC);
p1 = divrr(p2,mppi(DEFAULTPREC));
p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1));
s = 1000;
if (signe(p2))
{
s = p2[2];
if (lgefint(p2)>3 || s<0)
err(talker,"discriminant too big in classno");
if (s<1000) s=1000;
}
c=ptforme=0; court = stoi(2);
while (c <= s && *p)
{
c += *p++; k=krogs(x,c);
if (k)
{
av2=avma;
if (k>0)
{
divrsz(mulsr(c,p1),c-1,p1); avma=av2;
if (ptforme<100 && c>2)
{
court[2]=c;
formes[ptforme++]=redimag(primeform(x,court,0));
}
}
else { divrsz(mulsr(c,p1),c+1,p1); avma=av2; }
}
}
s = 2*itos(gceil(gpui(p1,ginv(stoi(4)),DEFAULTPREC)));
if (s>=10000) s=10000;
count = new_chunk(256); for (i=0; i<=255; i++) count[i]=0;
index = new_chunk(257);
tabla = new_chunk(10000);
tablb = new_chunk(10000);
hash = new_chunk(10000);
h = hin = ground(p1); f=formes[0];
p1 = fh = gpui(f,h,0);
for (i=0; i<s; i++)
{
p2=(GEN)p1[1]; tabla[i]=p2[lgefint(p2)-1];
p2=(GEN)p1[2]; tablb[i]=p2[lgefint(p2)-1];
count[tabla[i]&255]++; p1=compimag(p1,f);
}
index[0]=0; for (i=0; i<=254; i++) index[i+1] = index[i]+count[i];
for (i=0; i<s; i++) hash[index[tabla[i]&255]++]=i;
index[0]=0; for (i=0; i<=255; i++) index[i+1] = index[i]+count[i];
fg=gpuigs(f,s); av2=avma; lim=stack_lim(av2,2);
ftest=gpuigs(p1,0);
for (com=0; ; com++)
{
p1=(GEN)ftest[1]; k=p1[lgefint(p1)-1]; j=k&255;
for (j1=index[j]; j1 < index[j+1]; j1++)
{
j2=hash[j1];
if (tabla[j2] == k)
{
p2=(GEN)ftest[2];
if (tablb[j2] == labs(p2[lgefint(p2)-1]))
{
p1 = gmul(gpuigs(f,j2),fh);
if (gegal((GEN) p1[1],(GEN) ftest[1]) &&
absi_equal((GEN) p1[2],(GEN) ftest[2]))
{
/* p1 = ftest or ftest^(-1), we are done */
if (signe(p1[2]) == signe(ftest[2])) com = -com;
h = addii(addis(h,j2), mulss(s,com));
h = end_classno(h,hin,formes,ptforme);
avma = av; return icopy(h); /* safe: stack big enough */
}
}
}
}
ftest = gmul(ftest,fg);
if (is_pm1(ftest[1])) err(impl,"classno with too small order");
if (low_stack(lim, stack_lim(av2,2))) ftest = gerepileupto(av2,ftest);
}
}
/* use Euler products */
GEN
classno2(GEN x)
{
long av=avma,tetpil,n,i,k,s=signe(x),fl2;
GEN p1,p2,p3,p4,p5,p7,p8,pi4,reg,logd,d,fd;
if (typ(x) != t_INT) err(arither1);
if (!s) err(arither2);
if (s<0 && gcmpgs(x,-12) >= 0) return gun;
p1=auxdecomp(absi(x),1); p2=(GEN)p1[2]; p1=(GEN)p1[1];
n=lg(p1); d=gun;
for (i=1; i<n; i++)
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);
p1=gsqrt(gdiv(gmul(d,logd),gmul2n(pi4,1)),DEFAULTPREC);
p4=divri(pi4,d); p7=ginv(mpsqrt(pi4));
if (s>0)
{
reg=regula(d,DEFAULTPREC);
p2=gsubsg(1,gmul2n(gdiv(glog(reg,DEFAULTPREC),logd),1));
p3=gsqrt(gdivsg(2,logd),DEFAULTPREC);
if (gcmp(p2,p3)>=0) p1=gmul(p2,p1);
}
p1 = gtrunc(p1); n=p1[2];
if (lgefint(p1)!=3 || n<0)
err(talker,"discriminant too large in classno");
p1 = gsqrt(d,DEFAULTPREC); p3 = gzero;
if (s>0)
{
for (i=1; i<=n; i++)
{
k=krogs(fd,i);
if (k)
{
p2=mulir(mulss(i,i),p4);
p5=subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
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);
if (!egalii(x,fd)) /* x != p3 */
p8 = gdiv(p8,ground(gdiv(regula(x,DEFAULTPREC),reg)));
}
else
{
p1 = gdiv(p1,pi4);
for (i=1; i<=n; i++)
{
k=krogs(fd,i);
if (k)
{
p2=mulir(mulss(i,i),p4);
p5=subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
p5=addrr(p5,divrr(divrs(p1,i),mpexp(p2)));
p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
}
}
}
p3=ground(p3); tetpil=avma;
return gerepile(av,tetpil,mulii(p8,p3));
}
GEN
hclassno(GEN x)
{
long d,a,b,h,b2,f;
d= -itos(x); if (d>0 || (d & 3) > 1) return gzero;
if (!d) return gdivgs(gun,-12);
if (-d > (VERYBIGINT>>1))
err(talker,"discriminant too big in hclassno. Use quadclassunit");
h=0; b=d&1; b2=(b-d)>>2; f=0;
if (!b)
{
for (a=1; a*a<b2; a++)
if (b2%a == 0) h++;
f = (a*a==b2); b=2; b2=(4-d)>>2;
}
while (b2*3+d<0)
{
if (b2%b == 0) h++;
for (a=b+1; a*a<b2; a++)
if (b2%a == 0) h+=2;
if (a*a==b2) h++;
b+=2; b2=(b*b-d)>>2;
}
if (b2*3+d==0)
{
GEN y = cgetg(3,t_FRAC);
y[1]=lstoi(3*h+1); y[2]=lstoi(3); return y;
}
if (f) return gaddsg(h,ghalf);
return stoi(h);
}
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to arith1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / basemath