/* $Id: bibli2.c,v 1.17.2.2 2001/06/14 19:34:20 karim Exp $
Copyright (C) 2000 The PARI group.
This file is part of the PARI/GP package.
PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.
Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
/********************************************************************/
/** **/
/** BIBLIOTHEQUE MATHEMATIQUE **/
/** (deuxieme partie) **/
/** **/
/********************************************************************/
#include "pari.h"
/********************************************************************/
/** **/
/** DEVELOPPEMENTS LIMITES **/
/** **/
/********************************************************************/
GEN
tayl(GEN x, long v, long precdl)
{
long tetpil,i,vx = gvar9(x), av=avma;
GEN p1,y;
if (v <= vx)
{
long p1[] = { evaltyp(t_SER)|m_evallg(2), 0 };
p1[1] = evalvalp(precdl) | evalvarn(v);
return gadd(p1,x);
}
p1=cgetg(v+2,t_VEC);
for (i=0; i<v; i++) p1[i+1]=lpolx[i];
p1[vx+1]=lpolx[v]; p1[v+1]=lpolx[vx];
y = tayl(changevar(x,p1), vx,precdl); tetpil=avma;
return gerepile(av,tetpil, changevar(y,p1));
}
GEN
grando0(GEN x, long n, long do_clone)
{
long m, v, tx=typ(x);
GEN y;
if (gcmp0(x)) err(talker,"zero argument in O()");
if (tx == t_INT)
{
if (!gcmp1(x)) /* bug 3 + O(1). We suppose x is a truc() */
{
y=cgetg(5,t_PADIC);
y[1] = evalvalp(n) | evalprecp(0);
y[2] = do_clone? lclone(x): licopy(x);
y[3] = un; y[4] = zero; return y;
}
v=0; m=0; /* 1 = x^0 */
}
else
{
if (tx != t_POL && ! is_rfrac_t(tx))
err(talker,"incorrect argument in O()");
v=gvar(x); m=n*gval(x,v);
}
return zeroser(v,m);
}
/*******************************************************************/
/** **/
/** SPECIAL POLYNOMIALS **/
/** **/
/*******************************************************************/
#ifdef LONG_IS_64BIT
# define SQRTVERYBIGINT 3037000500 /* ceil(sqrt(VERYBIGINT)) */
#else
# define SQRTVERYBIGINT 46341
#endif
/* Tchebichev polynomial: T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2)
* T(n) = (n/2) sum_{k=0}^{n/2} a_k x^(n-2k)
* where a_k = (-1)^k 2^(n-2k) (n-k-1)! / k!(n-2k)! is an integer
* and a_0 = 2^(n-1), a_k / a_{k-1} = - (n-2k+2)(n-2k+1) / 4k(n-k) */
GEN
tchebi(long n, long v) /* Assume 4*n < VERYBIGINT */
{
long av,k,l;
GEN q,a,r;
if (v<0) v = 0;
if (n==0) return polun[v];
if (n==1) return polx[v];
q = cgetg(n+3, t_POL); r = q + n+2;
a = shifti(gun, n-1);
*r-- = (long)a;
*r-- = zero;
if (n < SQRTVERYBIGINT)
for (k=1,l=n; l>1; k++,l-=2)
{
av = avma;
a = divis(mulis(a, l*(l-1)), 4*k*(n-k));
a = gerepileuptoint(av, negi(a));
*r-- = (long)a;
*r-- = zero;
}
else
for (k=1,l=n; l>1; k++,l-=2)
{
av = avma;
a = mulis(mulis(a, l), l-1);
a = divis(divis(a, 4*k), n-k);
a = gerepileuptoint(av, negi(a));
*r-- = (long)a;
*r-- = zero;
}
q[1] = evalsigne(1) | evalvarn(v) | evallgef(n+3);
return q;
}
GEN addshiftw(GEN x, GEN y, long d);
/* Legendre polynomial */
/* L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1) */
GEN
legendre(long n, long v)
{
long av,tetpil,m,lim;
GEN p0,p1,p2;
if (v<0) v = 0;
if (n==0) return polun[v];
if (n==1) return polx[v];
p0=polun[v]; av=avma; lim=stack_lim(av,2);
p1=gmul2n(polx[v],1);
for (m=1; m<n; m++)
{
p2 = addshiftw(gmulsg(4*m+2,p1), gmulsg(-4*m,p0), 1);
setvarn(p2,v);
p0 = p1; tetpil=avma; p1 = gdivgs(p2,m+1);
if (low_stack(lim, stack_lim(av,2)))
{
GEN *gptr[2];
if(DEBUGMEM>1) err(warnmem,"legendre");
p0=gcopy(p0); gptr[0]=&p0; gptr[1]=&p1;
gerepilemanysp(av,tetpil,gptr,2);
}
}
tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-n));
}
/* cyclotomic polynomial */
GEN
cyclo(long n, long v)
{
long av=avma,tetpil,d,q,m;
GEN yn,yd;
if (n<=0) err(arither2);
if (v<0) v = 0;
yn = yd = polun[0];
for (d=1; d*d<=n; d++)
{
if (n%d) continue;
q=n/d;
m = mu(stoi(q));
if (m)
{ /* y *= (x^d - 1) */
if (m>0) yn = addshiftw(yn, gneg(yn), d);
else yd = addshiftw(yd, gneg(yd), d);
}
if (q==d) break;
m = mu(stoi(d));
if (m)
{ /* y *= (x^q - 1) */
if (m>0) yn = addshiftw(yn, gneg(yn), q);
else yd = addshiftw(yd, gneg(yd), q);
}
}
tetpil=avma; yn = gerepile(av,tetpil,gdeuc(yn,yd));
setvarn(yn,v); return yn;
}
/* compute prod (L*x ± a[i]) */
GEN
roots_to_pol_intern(GEN L, GEN a, long v, int plus)
{
long i,k,lx = lg(a), code;
GEN p1,p2;
if (lx == 1) return polun[v];
p1 = cgetg(lx, t_VEC);
code = evalsigne(1)|evalvarn(v)|evallgef(5);
for (k=1,i=1; i<lx-1; i+=2)
{
p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
p2[2] = lmul((GEN)a[i],(GEN)a[i+1]);
p2[3] = ladd((GEN)a[i],(GEN)a[i+1]);
if (plus == 0) p2[3] = lneg((GEN)p2[3]);
p2[4] = (long)L; p2[1] = code;
}
if (i < lx)
{
p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
p2[1] = code = evalsigne(1)|evalvarn(v)|evallgef(4);
p2[2] = plus? a[i]: lneg((GEN)a[i]);
p2[3] = (long)L;
}
setlg(p1, k); return divide_conquer_prod(p1, gmul);
}
GEN
roots_to_pol(GEN a, long v)
{
return roots_to_pol_intern(gun,a,v,0);
}
/* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
GEN
roots_to_pol_r1r2(GEN a, long r1, long v)
{
long i,k,lx = lg(a), code;
GEN p1;
if (lx == 1) return polun[v];
p1 = cgetg(lx, t_VEC);
code = evalsigne(1)|evalvarn(v)|evallgef(5);
for (k=1,i=1; i<r1; i+=2)
{
GEN p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
p2[2] = lmul((GEN)a[i],(GEN)a[i+1]);
p2[3] = lneg(gadd((GEN)a[i],(GEN)a[i+1]));
p2[4] = un; p2[1] = code;
}
if (i < r1+1)
p1[k++] = ladd(polx[v], gneg((GEN)a[i]));
for (i=r1+1; i<lx; i++)
{
GEN p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
p2[2] = lnorm((GEN)a[i]);
p2[3] = lneg(gtrace((GEN)a[i]));
p2[4] = un; p2[1] = code;
}
setlg(p1, k); return divide_conquer_prod(p1, gmul);
}
/* finds an equation for the d-th degree subfield of Q(zeta_n).
* (Z/nZ)* must be cyclic.
*/
GEN
subcyclo(GEN nn, GEN dd, int v)
{
long av=avma,tetpil,i,j,k,prec,q,d,p,pp,al,n,ex0,ex,aad,aa;
GEN a,z,pol,fa,powz,alpha;
if (typ(dd)!=t_INT || signe(dd)<=0) err(typeer,"subcyclo");
if (is_bigint(dd)) err(talker,"degree too large in subcyclo");
if (typ(nn)!=t_INT || signe(nn)<=0) err(typeer,"subcyclo");
if (v<0) v = 0;
d=itos(dd); if (d==1) return polx[v];
if (is_bigint(nn)) err(impl,"subcyclo for huge cyclotomic fields");
n = nn[2]; if ((n & 3) == 2) n >>= 1;
if (n == 1) err(talker,"degree does not divide phi(n) in subcyclo");
fa = factor(stoi(n));
p = itos(gmael(fa,1,1));
al= itos(gmael(fa,2,1));
if (lg((GEN)fa[1]) > 2 || (p==2 && al>2))
err(impl,"subcyclo in non-cyclic case");
if (d < n)
{
k = 1 + svaluation(d,p,&i);
if (k<al) { al = k; nn = gpowgs(stoi(p),al); n = nn[2]; }
}
avma=av; q = (n/p)*(p-1); /* = phi(n) */
if (q == d) return cyclo(n,v);
if (q % d) err(talker,"degree does not divide phi(n) in subcyclo");
q /= d;
if (p==2)
{
pol = powgi(polx[v],gdeux); pol[2]=un; /* replace gzero */
return pol; /* = x^2 + 1 */
}
a=gener(stoi(n)); aa = mael(a,2,2);
a=gpowgs(a,d); aad = mael(a,2,2);
#if 1
prec = expi(binome(stoi(d*q-1),d)) + expi(stoi(n));
prec = 2 + (prec>>TWOPOTBITS_IN_LONG);
if (prec<DEFAULTPREC) prec = DEFAULTPREC;
if (DEBUGLEVEL) fprintferr("subcyclo prec = %ld\n",prec);
z = cgetg(3,t_COMPLEX); a=mppi(prec); setexpo(a,2); /* a = 2\pi */
gsincos(divrs(a,n),(GEN*)(z+2),(GEN*)(z+1),prec); /* z = e_n(1) */
powz = cgetg(n,t_VEC); powz[1] = (long)z;
k = (n+3)>>1;
for (i=2; i<k; i++) powz[i] = lmul(z,(GEN)powz[i-1]);
if ((q&1) == 0) /* totally real field, take real part */
{
for (i=1; i<k; i++) powz[i] = mael(powz,i,1);
for ( ; i<n; i++) powz[i] = powz[n-i];
}
else
for ( ; i<n; i++) powz[i] = lconj((GEN)powz[n-i]);
alpha = cgetg(d+1,t_VEC) + 1; pol=gun;
for (ex0=1,k=0; k<d; k++, ex0=(ex0*aa)%n)
{
GEN p1 = gzero;
long av1 = avma; (void)new_chunk(2*prec + 3);
for (ex=ex0,i=0; i<q; i++)
{
for (pp=ex,j=0; j<al; j++)
{
p1 = gadd(p1,(GEN)powz[pp]);
pp = mulssmod(pp,p, n);
}
ex = mulssmod(ex,aad, n);
}
/* p1 = sum z^{p^k*h}, k = 0..al-1, h runs through the subgroup of order
* q = phi(n)/d of (Z/nZ)^* */
avma = av1; alpha[k] = lneg(p1);
}
pol = roots_to_pol_intern(gun,alpha-1,v, 1);
if (q&1) pol=greal(pol); /* already done otherwise */
tetpil=avma; return gerepile(av,tetpil,ground(pol));
#else
{
/* exact computation (much slower) */
GEN p1 = cgetg(n+2,t_POL)+2; for (i=0; i<n; i++) p1[i]=0;
for (ex=1,i=0; i<q; i++, ex=(ex*aad)%n)
for (pp=ex,j=0; j<al; j++, pp=(pp*p)%n) p1[pp]++;
for (i=0; i<n; i++) p1[i] = lstoi(p1[i]);
p1 = normalizepol_i(p1-2,n+2); setvarn(p1,v);
z = cyclo(n,v); a = caract2(z,gres(p1,z),v);
a = gdeuc(a, modulargcd(a,derivpol(a)));
return gerepileupto(av, a);
}
#endif
}
/********************************************************************/
/** **/
/** HILBERT & PASCAL MATRICES **/
/** **/
/********************************************************************/
GEN addshiftpol(GEN x, GEN y, long d);
GEN
mathilbert(long n) /* Hilbert matrix of order n */
{
long i,j;
GEN a,p;
if (n<0) n = 0;
p = cgetg(n+1,t_MAT);
for (j=1; j<=n; j++)
{
p[j]=lgetg(n+1,t_COL);
for (i=1+(j==1); i<=n; i++)
{
a=cgetg(3,t_FRAC); a[1]=un; a[2]=lstoi(i+j-1);
coeff(p,i,j)=(long)a;
}
}
if ( n ) mael(p,1,1)=un;
return p;
}
/* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
GEN
matqpascal(long n, GEN q)
{
long i,j,I, av = avma;
GEN m, *qpow = NULL; /* gcc -Wall */
if (n<0) n = -1;
n++; m = cgetg(n+1,t_MAT);
for (j=1; j<=n; j++) m[j] = lgetg(n+1,t_COL);
if (q)
{
I = (n+1)/2;
if (I > 1) { qpow = (GEN*)new_chunk(I+1); qpow[2]=q; }
for (j=3; j<=I; j++) qpow[j] = gmul(q, qpow[j-1]);
}
for (i=1; i<=n; i++)
{
I = (i+1)/2; coeff(m,i,1)=un;
if (q)
{
for (j=2; j<=I; j++)
coeff(m,i,j) = ladd(gmul(qpow[j],gcoeff(m,i-1,j)), gcoeff(m,i-1,j-1));
}
else
{
for (j=2; j<=I; j++)
coeff(m,i,j) = laddii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
}
for ( ; j<=i; j++) coeff(m,i,j) = coeff(m,i,i+1-j);
for ( ; j<=n; j++) coeff(m,i,j) = zero;
}
return gerepileupto(av, gcopy(m));
}
/********************************************************************/
/** **/
/** LAPLACE TRANSFORM (OF A SERIES) **/
/** **/
/********************************************************************/
GEN
laplace(GEN x)
{
long i,l,ec,av,tetpil;
GEN y,p1;
if (typ(x)!=t_SER) err(talker,"not a series in laplace");
if (gcmp0(x)) return gcopy(x);
av=avma; ec=valp(x);
if (ec<0) err(talker,"negative valuation in laplace");
l=lg(x); y=cgetg(l,t_SER);
p1=mpfact(ec); y[1]=x[1];
for (i=2; i<l; i++)
{
y[i]=lmul(p1,(GEN)x[i]);
ec++; p1=mulsi(ec,p1);
}
tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}
/********************************************************************/
/** **/
/** CONVOLUTION PRODUCT (OF TWO SERIES) **/
/** **/
/********************************************************************/
GEN
convol(GEN x, GEN y)
{
long l,i,j,v, vx=varn(x), lx=lg(x), ly=lg(y), ex=valp(x), ey=valp(y);
GEN z;
if (typ(x) != t_SER || typ(y) != t_SER)
err(talker,"not a series in convol");
if (gcmp0(x) || gcmp0(y))
err(talker,"zero series in convol");
if (varn(y) != vx)
err(talker,"different variables in convol");
v=ex; if (ey>v) v=ey;
l=ex+lx; i=ey+ly; if (i<l) l=i;
l -= v; if (l<3) err(talker,"non significant result in convol");
for (i=v+2; i < v+l; i++)
if (!gcmp0((GEN)x[i-ex]) && !gcmp0((GEN)y[i-ey])) { i++; break; }
if (i == l+v) return zeroser(vx, v+l-2);
z = cgetg(l-i+3+v,t_SER);
z[1] = evalsigne(1) | evalvalp(i-3) | evalvarn(vx);
for (j=i-1; j<l+v; j++) z[j-i+3]=lmul((GEN)x[j-ex],(GEN)y[j-ey]);
return z;
}
/******************************************************************/
/** **/
/** PRECISION CHANGES **/
/** **/
/******************************************************************/
GEN
gprec(GEN x, long l)
{
long tx=typ(x),lx=lg(x),i,pr;
GEN y;
if (l<=0) err(talker,"precision<=0 in gprec");
switch(tx)
{
case t_REAL:
pr = (long) (l*pariK1+3); y=cgetr(pr); affrr(x,y); break;
case t_PADIC:
y=cgetg(lx,tx); copyifstack(x[2], y[2]);
if (!signe(x[4]))
{
y[1]=evalvalp(l+precp(x)) | evalprecp(0);
y[3]=un; y[4]=zero; return y;
}
y[1]=x[1]; setprecp(y,l);
y[3]=lpuigs((GEN)x[2],l);
y[4]=lmodii((GEN)x[4],(GEN)y[3]);
break;
case t_SER:
if (gcmp0(x)) return zeroser(varn(x), l);
y=cgetg(l+2,t_SER); y[1]=x[1]; l++; i=l;
if (l>=lx)
for ( ; i>=lx; i--) y[i]=zero;
for ( ; i>=2; i--) y[i]=lcopy((GEN)x[i]);
break;
case t_POL:
lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=lprec((GEN)x[i],l);
break;
case t_COMPLEX: case t_POLMOD: case t_RFRAC: case t_RFRACN:
case t_VEC: case t_COL: case t_MAT:
y=cgetg(lx,tx);
for (i=1; i<lx; i++) y[i]=lprec((GEN)x[i],l);
break;
default: y=gcopy(x);
}
return y;
}
/* internal: precision given in word length (including codewords) */
GEN
gprec_w(GEN x, long pr)
{
long tx=typ(x),lx,i;
GEN y;
switch(tx)
{
case t_REAL:
y=cgetr(pr); affrr(x,y); break;
case t_POL:
lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
for (i=2; i<lx; i++) y[i]=(long)gprec_w((GEN)x[i],pr);
break;
case t_COMPLEX: case t_POLMOD: case t_RFRAC: case t_RFRACN:
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)gprec_w((GEN)x[i],pr);
break;
default: y=gprec(x,pr);
}
return y;
}
/*******************************************************************/
/** **/
/** RECIPROCAL POLYNOMIAL **/
/** **/
/*******************************************************************/
GEN
polrecip(GEN x)
{
long lx=lgef(x),i,j;
GEN y;
if (typ(x) != t_POL) err(typeer,"polrecip");
y=cgetg(lx,t_POL); y[1]=x[1];
for (i=2,j=lx-1; i<lx; i++,j--) y[i]=lcopy((GEN)x[j]);
return normalizepol_i(y,lx);
}
/* as above. Internal (don't copy or normalize) */
GEN
polrecip_i(GEN x)
{
long lx=lgef(x),i,j;
GEN y;
y=cgetg(lx,t_POL); y[1]=x[1];
for (i=2,j=lx-1; i<lx; i++,j--) y[i]=x[j];
return y;
}
/*******************************************************************/
/** **/
/** BINOMIAL COEFFICIENTS **/
/** **/
/*******************************************************************/
GEN
binome(GEN n, long k)
{
long av,i;
GEN y;
if (k <= 1)
{
if (is_noncalc_t(typ(n))) err(typeer,"binomial");
if (k < 0) return gzero;
if (k == 0) return gun;
return gcopy(n);
}
av = avma; y = n;
if (typ(n) == t_INT)
{
if (signe(n) > 0)
{
GEN z = subis(n,k);
if (cmpis(z,k) < 0) k = itos(z);
avma = av;
if (k <= 1)
{
if (k < 0) return gzero;
if (k == 0) return gun;
return gcopy(n);
}
}
for (i=2; i<=k; i++)
y = gdivgs(gmul(y,addis(n,i-1-k)), i);
}
else
{
for (i=2; i<=k; i++)
y = gdivgs(gmul(y,gaddgs(n,i-1-k)), i);
}
return gerepileupto(av, y);
}
/********************************************************************/
/** **/
/** POLYNOMIAL INTERPOLATION **/
/** **/
/********************************************************************/
/* assume n > 1 */
GEN
polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy)
{
long av = avma,tetpil,i,m, ns=0, tx=typ(x);
GEN den,ho,hp,w,y,c,d,dy;
if (!xa)
{
xa = cgetg(n+1, t_VEC);
for (i=1; i<=n; i++) xa[i] = lstoi(i);
xa++;
}
if (is_scalar_t(tx) && tx != t_INTMOD && tx != t_PADIC && tx != t_POLMOD)
{
GEN dif = NULL, dift;
for (i=0; i<n; i++)
{
dift = gabs(gsub(x,(GEN)xa[i]), MEDDEFAULTPREC);
if (!dif || gcmp(dift,dif)<0) { ns=i; dif=dift; }
}
}
c=new_chunk(n);
d=new_chunk(n); for (i=0; i<n; i++) c[i] = d[i] = ya[i];
y=(GEN)d[ns--];
dy = NULL; tetpil = 0; /* gcc -Wall */
for (m=1; m<n; m++)
{
for (i=0; i<n-m; i++)
{
ho = gsub((GEN)xa[i],x);
hp = gsub((GEN)xa[i+m],x); den = gsub(ho,hp);
if (gcmp0(den)) err(talker,"two abcissas are equal in polint");
w=gsub((GEN)c[i+1],(GEN)d[i]); den = gdiv(w,den);
c[i]=lmul(ho,den);
d[i]=lmul(hp,den);
}
dy = (2*(ns+1) < n-m)? (GEN)c[ns+1]: (GEN)d[ns--];
tetpil=avma; y=gadd(y,dy);
}
if (!ptdy) y = gerepile(av,tetpil,y);
else
{
GEN *gptr[2];
*ptdy=gcopy(dy); gptr[0]=&y; gptr[1]=ptdy;
gerepilemanysp(av,tetpil,gptr,2);
}
return y;
}
GEN
polint(GEN xa, GEN ya, GEN x, GEN *ptdy)
{
long tx=typ(xa), ty, lx=lg(xa);
if (ya) ty = typ(ya); else { ya = xa; ty = tx; xa = NULL; }
if (! is_vec_t(tx) || ! is_vec_t(ty))
err(talker,"not vectors in polinterpolate");
if (lx != lg(ya))
err(talker,"different lengths in polinterpolate");
if (lx <= 2)
{
if (lx == 1) err(talker,"no data in polinterpolate");
ya=gcopy((GEN)ya[1]); if (ptdy) *ptdy = ya;
return ya;
}
if (!x) x = polx[0];
return polint_i(xa? xa+1: xa,ya+1,x,lx-1,ptdy);
}
/***********************************************************************/
/* */
/* SET OPERATIONS */
/* */
/***********************************************************************/
static GEN
gtostr(GEN x)
{
char *s=GENtostr(x);
x = strtoGENstr(s,0); free(s); return x;
}
GEN
gtoset(GEN x)
{
long i,c,av,tetpil,tx,lx;
GEN y;
if (!x) return cgetg(1, t_VEC);
tx = typ(x); lx = lg(x);
if (!is_vec_t(tx))
{
if (tx != t_LIST)
{ y=cgetg(2,t_VEC); y[1]=(long)gtostr(x); return y; }
lx = lgef(x)-1; x++;
}
if (lx==1) return cgetg(1,t_VEC);
av=avma; y=cgetg(lx,t_VEC);
for (i=1; i<lx; i++) y[i]=(long)gtostr((GEN)x[i]);
y = sort(y);
c=1;
for (i=2; i<lx; i++)
if (!gegal((GEN)y[i], (GEN)y[c])) y[++c] = y[i];
tetpil=avma; setlg(y,c+1);
return gerepile(av,tetpil,gcopy(y));
}
long
setisset(GEN x)
{
long lx,i;
if (typ(x)!=t_VEC) return 0;
lx=lg(x)-1; if (!lx) return 1;
for (i=1; i<lx; i++)
if (typ(x[i]) != t_STR || gcmp((GEN)x[i+1],(GEN)x[i])<=0) return 0;
return typ(x[i]) == t_STR;
}
/* looks if y belongs to the set x and returns the index if yes, 0 if no */
long
setsearch(GEN x, GEN y, long flag)
{
long av = avma,lx,j,li,ri,fl, tx = typ(x);
if (tx==t_VEC) lx = lg(x);
else
{
if (tx!=t_LIST) err(talker,"not a set in setsearch");
lx=lgef(x)-1; x++;
}
if (lx==1) return flag? 1: 0;
if (typ(y) != t_STR) y = gtostr(y);
li=1; ri=lx-1;
do
{
j = (ri+li)>>1; fl = gcmp((GEN)x[j],y);
if (!fl) { avma=av; return flag? 0: j; }
if (fl<0) li=j+1; else ri=j-1;
} while (ri>=li);
avma=av; if (!flag) return 0;
return (fl<0)? j+1: j;
}
GEN
setunion(GEN x, GEN y)
{
long av=avma,tetpil;
GEN z;
if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion");
z=concatsp(x,y); tetpil=avma; return gerepile(av,tetpil,gtoset(z));
}
GEN
setintersect(GEN x, GEN y)
{
long av=avma,tetpil,i,lx,c;
GEN z;
if (!setisset(x) || !setisset(y)) err(talker,"not a set in setintersect");
lx=lg(x); z=cgetg(lx,t_VEC); c=1;
for (i=1; i<lx; i++)
if (setsearch(y, (GEN)x[i], 0)) z[c++] = x[i];
tetpil=avma; setlg(z,c);
return gerepile(av,tetpil,gcopy(z));
}
GEN
setminus(GEN x, GEN y)
{
long av=avma,tetpil,i,lx,c;
GEN z;
if (!setisset(x) || !setisset(y)) err(talker,"not a set in setminus");
lx=lg(x); z=cgetg(lx,t_VEC); c=1;
for (i=1; i<lx; i++)
if (setsearch(y, (GEN)x[i], 1)) z[c++] = x[i];
tetpil=avma; setlg(z,c);
return gerepile(av,tetpil,gcopy(z));
}
/***********************************************************************/
/* */
/* OPERATIONS ON DIRICHLET SERIES */
/* */
/***********************************************************************/
/* Addition, subtraction and scalar multiplication of Dirichlet series
are done on the corresponding vectors */
static long
dirval(GEN x)
{
long i=1,lx=lg(x);
while (i<lx && gcmp0((GEN)x[i])) i++;
return i;
}
GEN
dirmul(GEN x, GEN y)
{
long lx,ly,lz,dx,dy,av,tetpil,i,j,k,lim;
GEN z,p1;
if (typ(x)!=t_VEC || typ(y)!=t_VEC) err(talker,"not a dirseries in dirmul");
av=avma; dx=dirval(x); dy=dirval(y); lx=lg(x); ly=lg(y);
if (ly-dy<lx-dx) { z=y; y=x; x=z; lz=ly; ly=lx; lx=lz; lz=dy; dy=dx; dx=lz; }
lz=min(lx*dy,ly*dx);
lim = stack_lim(av,1);
z=cgetg(lz,t_VEC); for (i=1; i<lz; i++) z[i]=zero;
for (j=dx; j<lx; j++)
{
p1=(GEN)x[j];
if (!gcmp0(p1))
{
if (gcmp1(p1))
for (k=dy,i=j*dy; i<lz; i+=j,k++) z[i]=ladd((GEN)z[i],(GEN)y[k]);
else
{
if (gcmp_1(p1))
for (k=dy,i=j*dy; i<lz; i+=j,k++) z[i]=lsub((GEN)z[i],(GEN)y[k]);
else
for (k=dy,i=j*dy; i<lz; i+=j,k++) z[i]=ladd((GEN)z[i],gmul(p1,(GEN)y[k]));
}
}
if (low_stack(lim, stack_lim(av,1)))
{
if (DEBUGLEVEL) fprintferr("doubling stack in dirmul\n");
z = gerepileupto(av,gcopy(z));
}
}
tetpil=avma; return gerepile(av,tetpil,gcopy(z));
}
GEN
dirdiv(GEN x, GEN y)
{
long lx,ly,lz,dx,dy,av,tetpil,i,j;
GEN z,p1;
if (typ(x)!=t_VEC || typ(y)!=t_VEC) err(talker,"not a dirseries in dirmul");
av=avma; dx=dirval(x); dy=dirval(y); lx=lg(x); ly=lg(y);
if (dy!=1) err(talker,"not an invertible dirseries in dirdiv");
lz=min(lx,ly*dx); p1=(GEN)y[1];
if (!gcmp1(p1)) { y=gdiv(y,p1); x=gdiv(x,p1); }
else x=gcopy(x);
z=cgetg(lz,t_VEC); for (i=1; i<dx; i++) z[i]=zero;
for (j=dx; j<lz; j++)
{
p1=(GEN)x[j]; z[j]=(long)p1;
if (!gcmp0(p1))
{
if (gcmp1(p1))
for (i=j+j; i<lz; i+=j) x[i]=lsub((GEN)x[i],(GEN)y[i/j]);
else
{
if (gcmp_1(p1))
for (i=j+j; i<lz; i+=j) x[i]=ladd((GEN)x[i],(GEN)y[i/j]);
else
for (i=j+j; i<lz; i+=j) x[i]=lsub((GEN)x[i],gmul(p1,(GEN)y[i/j]));
}
}
}
tetpil=avma; return gerepile(av,tetpil,gcopy(z));
}
/*************************************************************************/
/** **/
/** RANDOM **/
/** **/
/*************************************************************************/
static long pari_randseed = 1;
/* BSD rand gives this: seed = 1103515245*seed + 12345 */
long
mymyrand(void)
{
#if BITS_IN_RANDOM == 64
pari_randseed = (1000000000000654397*pari_randseed + 12347) & ~HIGHBIT;
#else
pari_randseed = (1000276549*pari_randseed + 12347) & 0x7fffffff;
#endif
return pari_randseed;
}
GEN muluu(ulong x, ulong y);
static ulong
gp_rand(void)
{
#define GLUE2(hi, lo) (((hi) << BITS_IN_HALFULONG) | (lo))
#if !defined(LONG_IS_64BIT) || BITS_IN_RANDOM == 64
return GLUE2((mymyrand()>>12)&LOWMASK,
(mymyrand()>>12)&LOWMASK);
#else
#define GLUE4(hi1,hi2, lo1,lo2) GLUE2(((hi1)<<16)|(hi2), ((lo1)<<16)|(lo2))
# define LOWMASK2 0xffffUL
return GLUE4((mymyrand()>>12)&LOWMASK2,
(mymyrand()>>12)&LOWMASK2,
(mymyrand()>>12)&LOWMASK2,
(mymyrand()>>12)&LOWMASK2);
#endif
}
GEN
genrand(GEN N)
{
long lx,i,nz;
GEN x, p1;
if (!N) return stoi(mymyrand());
if (typ(N)!=t_INT || signe(N)<=0) err(talker,"invalid bound in random");
lx = lgefint(N); x = new_chunk(lx);
nz = lx-1; while (!N[nz]) nz--; /* nz = index of last non-zero word */
for (i=2; i<lx; i++)
{
long n = N[i], r;
if (n == 0) r = 0;
else
{
long av = avma;
if (i < nz) n++; /* allow for equality if we can go down later */
p1 = muluu(n, gp_rand()); /* < n * 2^32, so 0 <= first word < n */
r = (lgefint(p1)<=3)? 0: p1[2]; avma = av;
}
x[i] = r;
if (r < (ulong)N[i]) break;
}
for (i++; i<lx; i++) x[i] = gp_rand();
i=2; while (i<lx && !x[i]) i++;
i -= 2; x += i; lx -= i;
x[1] = evalsigne(lx>2) | evallgefint(lx);
x[0] = evaltyp(t_INT) | evallg(lx);
avma = (long)x; return x;
}
long
setrand(long seed) { return (pari_randseed = seed); }
long
getrand(void) { return pari_randseed; }
long
getstack(void) { return top-avma; }
long
gettime(void) { return timer(); }
/***********************************************************************/
/** **/
/** PERMUTATIONS **/
/** **/
/***********************************************************************/
GEN
permute(long n, GEN x)
{
long av=avma,i,a,r;
GEN v,w;
if (n < 1) err(talker,"n too small (%ld) in numtoperm",n);
if (typ(x) != t_INT) err(arither1);
v = cgetg(n+1, t_VEC);
v[1]=1; av = avma;
if (signe(x) <= 0) x = modii(x, mpfact(n));
for (r=2; r<=n; r++)
{
x = dvmdis(x,r,&w); a = itos(w);
for (i=r; i>=a+2; i--) v[i] = v[i-1];
v[i]=r;
}
avma = av;
for (i=1; i<=n; i++) v[i] = lstoi(v[i]);
return v;
}
GEN
permuteInv(GEN x)
{
long av=avma, lx=lg(x)-1, ini=lx, last, ind, tx = typ(x);
GEN ary,res;
if (!is_vec_t(tx)) err(talker,"not a vector in permuteInv");
ary=cgetg(lx+1,t_VECSMALL);
for (ind=1; ind<=lx; ind++)
{
res = (GEN)*++x;
if (typ(res) != t_INT) err(typeer,"permuteInv");
ary[ind] = itos(res);
}
ary++; res = gzero;
for (last=lx; last>0; last--)
{
lx--; ind=lx;
while (ind>0 && ary[ind] != last) ind--;
res = addis(mulis(res,last), ind);
while (ind++<lx) ary[ind-1] = ary[ind];
}
if (!signe(res)) res=mpfact(ini);
return gerepileuptoint(av, res);
}
/********************************************************************/
/** **/
/** MODREVERSE **/
/** **/
/********************************************************************/
GEN
polymodrecip(GEN x)
{
long v,i,j,n,av,tetpil,lx;
GEN p1,p2,p3,p,phi,y,col;
if (typ(x)!=t_POLMOD) err(talker,"not a polymod in polymodrecip");
p=(GEN)x[1]; phi=(GEN)x[2];
v=varn(p); n=lgef(p)-3; if (n<=0) return gcopy(x);
if (n==1)
{
y=cgetg(3,t_POLMOD);
if (typ(phi)==t_POL) phi = (GEN)phi[2];
p1=cgetg(4,t_POL); p1[1]=p[1]; p1[2]=lneg(phi); p1[3]=un;
y[1]=(long)p1;
if (gcmp0((GEN)p[2])) p1 = zeropol(v);
else
{
p1=cgetg(3,t_POL); av=avma;
p1[1] = evalsigne(1) | evalvarn(n) | evallgef(3);
p2=gdiv((GEN)p[2],(GEN)p[3]); tetpil=avma;
p1[2] = lpile(av,tetpil,gneg(p2));
}
y[2]=(long)p1; return y;
}
if (gcmp0(phi) || typ(phi) != t_POL)
err(talker,"reverse polymod does not exist");
av=avma; y=cgetg(n+1,t_MAT);
y[1]=(long)gscalcol_i(gun,n);
p2=phi;
for (j=2; j<=n; j++)
{
lx=lgef(p2); p1=cgetg(n+1,t_COL); y[j]=(long)p1;
for (i=1; i<=lx-2; i++) p1[i]=p2[i+1];
for ( ; i<=n; i++) p1[i]=zero;
if (j<n) p2 = gmod(gmul(p2,phi), p);
}
col=cgetg(n+1,t_COL); col[1]=zero; col[2]=un;
for (i=3; i<=n; i++) col[i]=zero;
p1=gauss(y,col); p2=gtopolyrev(p1,v); p3=caract(x,v);
tetpil=avma; return gerepile(av,tetpil,gmodulcp(p2,p3));
}
/********************************************************************/
/** **/
/** HEAPSORT **/
/** **/
/********************************************************************/
static GEN vcmp_k;
static int vcmp_lk;
static int (*vcmp_cmp)(GEN,GEN);
int
pari_compare_int(int *a,int *b)
{
return *a - *b;
}
int
pari_compare_long(long *a,long *b)
{
return *a - *b;
}
static int
veccmp(GEN x, GEN y)
{
int i,s;
for (i=1; i<vcmp_lk; i++)
{
s = vcmp_cmp((GEN) x[vcmp_k[i]], (GEN) y[vcmp_k[i]]);
if (s) return s;
}
return 0;
}
static int
longcmp(GEN x, GEN y)
{
return ((long)x > (long)y)? 1: ((x == y)? 0: -1);
}
/* Sort x = vector of elts, using cmp to compare them.
* flag & cmp_IND: indirect sort: return permutation that would sort x
* For private use:
* flag & cmp_C : as cmp_IND, but return permutation as vector of C-longs
*/
GEN
gen_sort(GEN x, int flag, int (*cmp)(GEN,GEN))
{
long i,j,indxt,ir,l,tx=typ(x),lx=lg(x);
GEN q,y,indx;
if (!is_matvec_t(tx) && tx != t_VECSMALL) err(typeer,"gen_sort");
if (flag & cmp_C) tx = t_VECSMALL;
else if (flag & cmp_IND) tx = t_VEC;
y = cgetg(lx, tx);
if (lx==1) return y;
if (lx==2)
{
if (flag & cmp_C)
y[1] = 1;
else if (flag & cmp_IND)
y[1] = un;
else
y[1] = lcopy((GEN)x[1]);
return y;
}
if (!cmp) cmp = &longcmp;
indx = (GEN) gpmalloc(lx*sizeof(long));
for (j=1; j<lx; j++) indx[j]=j;
ir=lx-1; l=(ir>>1)+1;
for(;;)
{
if (l>1)
{ l--; indxt = indx[l]; }
else
{
indxt = indx[ir]; indx[ir]=indx[1]; ir--;
if (ir == 1)
{
indx[1] = indxt;
if (flag & cmp_C)
for (i=1; i<lx; i++) y[i]=indx[i];
else if (flag & cmp_IND)
for (i=1; i<lx; i++) y[i]=lstoi(indx[i]);
else
for (i=1; i<lx; i++) y[i]=lcopy((GEN)x[indx[i]]);
free(indx); return y;
}
}
q = (GEN)x[indxt]; i=l;
if (flag & cmp_REV)
for (j=i<<1; j<=ir; j<<=1)
{
if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) > 0) j++;
if (cmp(q,(GEN)x[indx[j]]) <= 0) break;
indx[i]=indx[j]; i=j;
}
else
for (j=i<<1; j<=ir; j<<=1)
{
if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) < 0) j++;
if (cmp(q,(GEN)x[indx[j]]) >= 0) break;
indx[i]=indx[j]; i=j;
}
indx[i]=indxt;
}
}
#define sort_fun(flag) ((flag & cmp_LEX)? &lexcmp: &gcmp)
GEN
gen_vecsort(GEN x, GEN k, long flag)
{
long i,j,l,t, lx = lg(x), tmp[2];
if (lx<=2) return gen_sort(x,flag,sort_fun(flag));
t = typ(k); vcmp_cmp = sort_fun(flag);
if (t==t_INT)
{
tmp[1] = (long)k; k = tmp;
vcmp_lk = 2;
}
else
{
if (! is_vec_t(t)) err(talker,"incorrect lextype in vecsort");
vcmp_lk = lg(k);
}
l = 0;
vcmp_k = (GEN)gpmalloc(vcmp_lk * sizeof(long));
for (i=1; i<vcmp_lk; i++)
{
j = itos((GEN)k[i]);
if (j<=0) err(talker,"negative index in vecsort");
vcmp_k[i]=j; if (j>l) l=j;
}
t = typ(x);
if (! is_matvec_t(t)) err(typeer,"vecsort");
for (j=1; j<lx; j++)
{
t = typ(x[j]);
if (! is_vec_t(t)) err(typeer,"vecsort");
if (lg((GEN)x[j]) <= l) err(talker,"index too large in vecsort");
}
x = gen_sort(x, flag, veccmp);
free(vcmp_k); return x;
}
GEN
vecsort0(GEN x, GEN k, long flag)
{
if (flag < 0 || flag >= cmp_C) err(flagerr,"vecsort");
return k? gen_vecsort(x,k,flag): gen_sort(x,flag, sort_fun(flag));
}
GEN
vecsort(GEN x, GEN k)
{
return gen_vecsort(x,k, 0);
}
GEN
sindexsort(GEN x)
{
return gen_sort(x, cmp_IND | cmp_C, gcmp);
}
GEN
sindexlexsort(GEN x)
{
return gen_sort(x, cmp_IND | cmp_C, lexcmp);
}
GEN
indexsort(GEN x)
{
return gen_sort(x, cmp_IND, gcmp);
}
GEN
indexlexsort(GEN x)
{
return gen_sort(x, cmp_IND, lexcmp);
}
GEN
sort(GEN x)
{
return gen_sort(x, 0, gcmp);
}
GEN
lexsort(GEN x)
{
return gen_sort(x, 0, lexcmp);
}
/* index of x in table T, 0 otherwise */
long
tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
{
long l=1,u=lg(T)-1,i,s;
while (u>=l)
{
i = (l+u)>>1; s = cmp(x,(GEN)T[i]);
if (!s) return i;
if (s<0) u=i-1; else l=i+1;
}
return 0;
}
/* assume lg(x) = lg(y), x,y in Z^n */
int
cmp_vecint(GEN x, GEN y)
{
long fl,i, lx = lg(x);
for (i=1; i<lx; i++)
if (( fl = cmpii((GEN)x[i], (GEN)y[i]) )) return fl;
return 0;
}
/* assume x and y come from the same primedec call (uniformizer unique) */
int
cmp_prime_over_p(GEN x, GEN y)
{
int k = mael(x,4,2) - mael(y,4,2); /* diff. between residue degree */
return k? ((k > 0)? 1: -1)
: cmp_vecint((GEN)x[2], (GEN)y[2]);
}
int
cmp_prime_ideal(GEN x, GEN y)
{
int k = cmpii((GEN)x[1], (GEN)y[1]);
return k? k: cmp_prime_over_p(x,y);
}