version 1.1, 2001/10/02 11:17:02 |
version 1.2, 2002/09/11 07:26:50 |
Line 20 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 20 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
/** **/ |
/** **/ |
/********************************************************************/ |
/********************************************************************/ |
#include "pari.h" |
#include "pari.h" |
|
extern GEN _ei(long n, long i); |
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** **/ |
Line 28 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
Line 29 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
/********************************************************************/ |
/********************************************************************/ |
|
|
GEN |
GEN |
tayl(GEN x, long v, long precdl) |
tayl(GEN x, long v, long precS) |
{ |
{ |
long tetpil,i,vx = gvar9(x), av=avma; |
long i, vx = gvar9(x); |
|
gpmem_t tetpil, av=avma; |
GEN p1,y; |
GEN p1,y; |
|
|
if (v <= vx) |
if (v <= vx) |
{ |
{ |
long p1[] = { evaltyp(t_SER)|m_evallg(2), 0 }; |
long p1[] = { evaltyp(t_SER)|_evallg(2), 0 }; |
p1[1] = evalvalp(precdl) | evalvarn(v); |
p1[1] = evalvalp(precS) | evalvarn(v); |
return gadd(p1,x); |
return gadd(p1,x); |
} |
} |
p1=cgetg(v+2,t_VEC); |
p1=cgetg(v+2,t_VEC); |
for (i=0; i<v; i++) p1[i+1]=lpolx[i]; |
for (i=0; i<v; i++) p1[i+1]=lpolx[i]; |
p1[vx+1]=lpolx[v]; p1[v+1]=lpolx[vx]; |
p1[vx+1]=lpolx[v]; p1[v+1]=lpolx[vx]; |
y = tayl(changevar(x,p1), vx,precdl); tetpil=avma; |
y = tayl(changevar(x,p1), vx,precS); tetpil=avma; |
return gerepile(av,tetpil, changevar(y,p1)); |
return gerepile(av,tetpil, changevar(y,p1)); |
} |
} |
|
|
|
|
grando0(GEN x, long n, long do_clone) |
grando0(GEN x, long n, long do_clone) |
{ |
{ |
long m, v, tx=typ(x); |
long m, v, tx=typ(x); |
GEN y; |
|
|
|
if (gcmp0(x)) err(talker,"zero argument in O()"); |
if (gcmp0(x)) err(talker,"zero argument in O()"); |
if (tx == t_INT) |
if (tx == t_INT) |
{ |
{ |
if (!gcmp1(x)) /* bug 3 + O(1). We suppose x is a truc() */ |
if (!gcmp1(x)) /* bug 3 + O(1). We suppose x is a truc() */ |
{ |
{ |
y=cgetg(5,t_PADIC); |
if (do_clone) x = gclone(x); |
y[1] = evalvalp(n) | evalprecp(0); |
return padiczero(x,n); |
y[2] = do_clone? lclone(x): licopy(x); |
|
y[3] = un; y[4] = zero; return y; |
|
} |
} |
v=0; m=0; /* 1 = x^0 */ |
v=0; m=0; /* 1 = x^0 */ |
} |
} |
Line 68 grando0(GEN x, long n, long do_clone) |
|
Line 67 grando0(GEN x, long n, long do_clone) |
|
{ |
{ |
if (tx != t_POL && ! is_rfrac_t(tx)) |
if (tx != t_POL && ! is_rfrac_t(tx)) |
err(talker,"incorrect argument in O()"); |
err(talker,"incorrect argument in O()"); |
v=gvar(x); m=n*gval(x,v); |
v = gvar(x); if (v > MAXVARN) err(talker,"incorrect object in O()"); |
|
m = n*gval(x,v); |
} |
} |
return zeroser(v,m); |
return zeroser(v,m); |
} |
} |
Line 91 grando0(GEN x, long n, long do_clone) |
|
Line 91 grando0(GEN x, long n, long do_clone) |
|
GEN |
GEN |
tchebi(long n, long v) /* Assume 4*n < VERYBIGINT */ |
tchebi(long n, long v) /* Assume 4*n < VERYBIGINT */ |
{ |
{ |
long av,k,l; |
long k, l; |
|
gpmem_t av; |
GEN q,a,r; |
GEN q,a,r; |
|
|
if (v<0) v = 0; |
if (v<0) v = 0; |
|
if (n < 0) err(talker,"negative degree in tchebi"); |
if (n==0) return polun[v]; |
if (n==0) return polun[v]; |
if (n==1) return polx[v]; |
if (n==1) return polx[v]; |
|
|
Line 131 GEN addshiftw(GEN x, GEN y, long d); |
|
Line 133 GEN addshiftw(GEN x, GEN y, long d); |
|
GEN |
GEN |
legendre(long n, long v) |
legendre(long n, long v) |
{ |
{ |
long av,tetpil,m,lim; |
long m; |
|
gpmem_t av, tetpil, lim; |
GEN p0,p1,p2; |
GEN p0,p1,p2; |
|
|
if (v<0) v = 0; |
if (v<0) v = 0; |
|
if (n < 0) err(talker,"negative degree in legendre"); |
if (n==0) return polun[v]; |
if (n==0) return polun[v]; |
if (n==1) return polx[v]; |
if (n==1) return polx[v]; |
|
|
Line 160 legendre(long n, long v) |
|
Line 164 legendre(long n, long v) |
|
GEN |
GEN |
cyclo(long n, long v) |
cyclo(long n, long v) |
{ |
{ |
long av=avma,tetpil,d,q,m; |
long d, q, m; |
|
gpmem_t av=avma, tetpil; |
GEN yn,yd; |
GEN yn,yd; |
|
|
if (n<=0) err(arither2); |
if (n<=0) err(arither2); |
Line 249 roots_to_pol_r1r2(GEN a, long r1, long v) |
|
Line 254 roots_to_pol_r1r2(GEN a, long r1, long v) |
|
setlg(p1, k); return divide_conquer_prod(p1, gmul); |
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 **/ |
/** HILBERT & PASCAL MATRICES **/ |
Line 363 mathilbert(long n) /* Hilbert matrix of order n */ |
|
Line 276 mathilbert(long n) /* Hilbert matrix of order n */ |
|
coeff(p,i,j)=(long)a; |
coeff(p,i,j)=(long)a; |
} |
} |
} |
} |
if ( n ) mael(p,1,1)=un; |
if ( n ) mael(p,1,1)=un; |
return p; |
return p; |
} |
} |
|
|
Line 371 mathilbert(long n) /* Hilbert matrix of order n */ |
|
Line 284 mathilbert(long n) /* Hilbert matrix of order n */ |
|
GEN |
GEN |
matqpascal(long n, GEN q) |
matqpascal(long n, GEN q) |
{ |
{ |
long i,j,I, av = avma; |
long i, j, I; |
|
gpmem_t av = avma; |
GEN m, *qpow = NULL; /* gcc -Wall */ |
GEN m, *qpow = NULL; /* gcc -Wall */ |
|
|
if (n<0) n = -1; |
if (n<0) n = -1; |
Line 411 matqpascal(long n, GEN q) |
|
Line 325 matqpascal(long n, GEN q) |
|
GEN |
GEN |
laplace(GEN x) |
laplace(GEN x) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long i,l,ec; |
long i,l,ec; |
GEN y,p1; |
GEN y,p1; |
|
|
Line 480 gprec(GEN x, long l) |
|
Line 394 gprec(GEN x, long l) |
|
pr = (long) (l*pariK1+3); y=cgetr(pr); affrr(x,y); break; |
pr = (long) (l*pariK1+3); y=cgetr(pr); affrr(x,y); break; |
|
|
case t_PADIC: |
case t_PADIC: |
y=cgetg(lx,tx); copyifstack(x[2], y[2]); |
|
if (!signe(x[4])) |
if (!signe(x[4])) |
{ |
return padiczero((GEN)x[2], l+precp(x)); |
y[1]=evalvalp(l+precp(x)) | evalprecp(0); |
y=cgetg(lx,tx); copyifstack(x[2], y[2]); |
y[3]=un; y[4]=zero; return y; |
|
} |
|
y[1]=x[1]; setprecp(y,l); |
y[1]=x[1]; setprecp(y,l); |
y[3]=lpuigs((GEN)x[2],l); |
y[3]=lpuigs((GEN)x[2],l); |
y[4]=lmodii((GEN)x[4],(GEN)y[3]); |
y[4]=lmodii((GEN)x[4],(GEN)y[3]); |
Line 580 polrecip_i(GEN x) |
|
Line 491 polrecip_i(GEN x) |
|
GEN |
GEN |
binome(GEN n, long k) |
binome(GEN n, long k) |
{ |
{ |
long av,i; |
long i; |
|
gpmem_t av; |
GEN y; |
GEN y; |
|
|
if (k <= 1) |
if (k <= 1) |
Line 616 binome(GEN n, long k) |
|
Line 528 binome(GEN n, long k) |
|
return gerepileupto(av, y); |
return gerepileupto(av, y); |
} |
} |
|
|
|
/* Assume n >= 1, return bin, bin[k+1] = binomial(n, k) */ |
|
GEN |
|
vecbinome(long n) |
|
{ |
|
long d = (n + 1)/2, k; |
|
GEN bin = cgetg(n+2, t_VEC), *C; |
|
C = (GEN*)(bin + 1); /* C[k] = binomial(n, k) */ |
|
C[0] = gun; |
|
for (k=1; k <= d; k++) |
|
{ |
|
gpmem_t av = avma; |
|
C[k] = gerepileuptoint(av, diviiexact(mulsi(n-k+1, C[k-1]), stoi(k))); |
|
} |
|
for ( ; k <= n; k++) C[k] = C[n - k]; |
|
return bin; |
|
} |
|
|
/********************************************************************/ |
/********************************************************************/ |
/** **/ |
/** **/ |
/** POLYNOMIAL INTERPOLATION **/ |
/** POLYNOMIAL INTERPOLATION **/ |
Line 625 binome(GEN n, long k) |
|
Line 554 binome(GEN n, long k) |
|
GEN |
GEN |
polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy) |
polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy) |
{ |
{ |
long av = avma,tetpil,i,m, ns=0, tx=typ(x); |
long i, m, ns=0, tx=typ(x); |
|
gpmem_t av = avma, tetpil; |
GEN den,ho,hp,w,y,c,d,dy; |
GEN den,ho,hp,w,y,c,d,dy; |
|
|
if (!xa) |
if (!xa) |
|
|
polint(GEN xa, GEN ya, GEN x, GEN *ptdy) |
polint(GEN xa, GEN ya, GEN x, GEN *ptdy) |
{ |
{ |
long tx=typ(xa), ty, lx=lg(xa); |
long tx=typ(xa), ty, lx=lg(xa); |
|
|
if (ya) ty = typ(ya); else { ya = xa; ty = tx; xa = NULL; } |
if (ya) ty = typ(ya); else { ya = xa; ty = tx; xa = NULL; } |
|
|
if (! is_vec_t(tx) || ! is_vec_t(ty)) |
if (! is_vec_t(tx) || ! is_vec_t(ty)) |
|
|
GEN |
GEN |
gtoset(GEN x) |
gtoset(GEN x) |
{ |
{ |
ulong av; |
gpmem_t av; |
long i,c,tx,lx; |
long i,c,tx,lx; |
GEN y; |
GEN y; |
|
|
|
|
|
|
/* looks if y belongs to the set x and returns the index if yes, 0 if no */ |
/* looks if y belongs to the set x and returns the index if yes, 0 if no */ |
long |
long |
setsearch(GEN x, GEN y, long flag) |
gen_search(GEN x, GEN y, int flag, int (*cmp)(GEN,GEN)) |
{ |
{ |
long av = avma,lx,j,li,ri,fl, tx = typ(x); |
long lx,j,li,ri,fl, tx = typ(x); |
|
|
if (tx==t_VEC) lx = lg(x); |
if (tx==t_VEC) lx = lg(x); |
else |
else |
Line 756 setsearch(GEN x, GEN y, long flag) |
|
Line 686 setsearch(GEN x, GEN y, long flag) |
|
} |
} |
if (lx==1) return flag? 1: 0; |
if (lx==1) return flag? 1: 0; |
|
|
if (typ(y) != t_STR) y = gtostr(y); |
|
li=1; ri=lx-1; |
li=1; ri=lx-1; |
do |
do |
{ |
{ |
j = (ri+li)>>1; fl = gcmp((GEN)x[j],y); |
j = (ri+li)>>1; fl = cmp((GEN)x[j],y); |
if (!fl) { avma=av; return flag? 0: j; } |
if (!fl) return flag? 0: j; |
if (fl<0) li=j+1; else ri=j-1; |
if (fl<0) li=j+1; else ri=j-1; |
} while (ri>=li); |
} while (ri>=li); |
avma=av; if (!flag) return 0; |
if (!flag) return 0; |
return (fl<0)? j+1: j; |
return (fl<0)? j+1: j; |
} |
} |
|
|
|
|
|
long |
|
setsearch(GEN x, GEN y, long flag) |
|
{ |
|
gpmem_t av = avma; |
|
long res; |
|
if (typ(y) != t_STR) y = gtostr(y); |
|
res=gen_search(x,y,flag,gcmp); |
|
avma=av; |
|
return res; |
|
} |
|
|
|
#if 0 |
GEN |
GEN |
|
gen_union(GEN x, GEN y, int (*cmp)(GEN,GEN)) |
|
{ |
|
if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion"); |
|
|
|
} |
|
#endif |
|
|
|
GEN |
setunion(GEN x, GEN y) |
setunion(GEN x, GEN y) |
{ |
{ |
long av=avma,tetpil; |
gpmem_t av=avma, tetpil; |
GEN z; |
GEN z; |
|
|
if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion"); |
if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion"); |
Line 781 setunion(GEN x, GEN y) |
|
Line 731 setunion(GEN x, GEN y) |
|
GEN |
GEN |
setintersect(GEN x, GEN y) |
setintersect(GEN x, GEN y) |
{ |
{ |
long av=avma,i,lx,c; |
long i, lx, c; |
|
gpmem_t av=avma; |
GEN z; |
GEN z; |
|
|
if (!setisset(x) || !setisset(y)) err(talker,"not a set in setintersect"); |
if (!setisset(x) || !setisset(y)) err(talker,"not a set in setintersect"); |
Line 792 setintersect(GEN x, GEN y) |
|
Line 743 setintersect(GEN x, GEN y) |
|
} |
} |
|
|
GEN |
GEN |
setminus(GEN x, GEN y) |
gen_setminus(GEN set1, GEN set2, int (*cmp)(GEN,GEN)) |
{ |
{ |
long av=avma,i,lx,c; |
gpmem_t ltop=avma; |
GEN z; |
long find; |
|
long i,j,k; |
|
GEN diff=cgetg(lg(set1),t_VEC); |
|
for(i=1,j=1,k=1; i < lg(set1); i++) |
|
{ |
|
for(find=0; j < lg(set2); j++) |
|
{ |
|
int s=cmp((GEN)set1[i],(GEN)set2[j]); |
|
if (s<0) break ; |
|
if (s>0) continue; |
|
find=1; |
|
} |
|
if (!find) |
|
diff[k++]=set1[i]; |
|
} |
|
setlg(diff,k); |
|
return gerepilecopy(ltop,diff); |
|
} |
|
|
|
GEN |
|
setminus(GEN x, GEN y) |
|
{ |
if (!setisset(x) || !setisset(y)) err(talker,"not a set in setminus"); |
if (!setisset(x) || !setisset(y)) err(talker,"not a set in setminus"); |
lx=lg(x); z=cgetg(lx,t_VEC); c=1; |
return gen_setminus(x,y,gcmp); |
for (i=1; i<lx; i++) |
|
if (setsearch(y, (GEN)x[i], 1)) z[c++] = x[i]; |
|
setlg(z,c); return gerepilecopy(av,z); |
|
} |
} |
|
|
/***********************************************************************/ |
/***********************************************************************/ |
|
|
GEN |
GEN |
dirmul(GEN x, GEN y) |
dirmul(GEN x, GEN y) |
{ |
{ |
ulong av = avma, lim = stack_lim(av,1); |
gpmem_t av = avma, lim = stack_lim(av, 1); |
long lx,ly,lz,dx,dy,i,j,k; |
long lx,ly,lz,dx,dy,i,j,k; |
GEN z,p1; |
GEN z,p1; |
|
|
Line 860 dirmul(GEN x, GEN y) |
|
Line 828 dirmul(GEN x, GEN y) |
|
GEN |
GEN |
dirdiv(GEN x, GEN y) |
dirdiv(GEN x, GEN y) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
long lx,ly,lz,dx,dy,i,j; |
long lx,ly,lz,dx,dy,i,j; |
GEN z,p1; |
GEN z,p1; |
|
|
|
|
ulong n = N[i], r; |
ulong n = N[i], r; |
if (n == 0) r = 0; |
if (n == 0) r = 0; |
else |
else |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
if (i < nz) n++; /* allow for equality if we can go down later */ |
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 */ |
p1 = muluu(n, gp_rand()); /* < n * 2^32, so 0 <= first word < n */ |
r = (lgefint(p1)<=3)? 0: p1[2]; avma = av; |
r = (lgefint(p1)<=3)? 0: p1[2]; avma = av; |
|
|
i -= 2; x += i; lx -= i; |
i -= 2; x += i; lx -= i; |
x[1] = evalsigne(lx>2) | evallgefint(lx); |
x[1] = evalsigne(lx>2) | evallgefint(lx); |
x[0] = evaltyp(t_INT) | evallg(lx); |
x[0] = evaltyp(t_INT) | evallg(lx); |
avma = (long)x; return x; |
avma = (gpmem_t)x; return x; |
} |
} |
|
|
long |
long |
Line 982 gettime(void) { return timer(); } |
|
Line 950 gettime(void) { return timer(); } |
|
GEN |
GEN |
numtoperm(long n, GEN x) |
numtoperm(long n, GEN x) |
{ |
{ |
ulong av; |
gpmem_t av; |
long i,a,r; |
long i,a,r; |
GEN v,w; |
GEN v,w; |
|
|
Line 1005 numtoperm(long n, GEN x) |
|
Line 973 numtoperm(long n, GEN x) |
|
GEN |
GEN |
permtonum(GEN x) |
permtonum(GEN x) |
{ |
{ |
long av=avma, lx=lg(x)-1, n=lx, last, ind, tx = typ(x); |
long lx=lg(x)-1, n=lx, last, ind, tx = typ(x); |
|
gpmem_t av=avma; |
GEN ary,res; |
GEN ary,res; |
|
|
if (!is_vec_t(tx)) err(talker,"not a vector in permtonum"); |
if (!is_vec_t(tx)) err(talker,"not a vector in permtonum"); |
Line 1033 permtonum(GEN x) |
|
Line 1002 permtonum(GEN x) |
|
/** MODREVERSE **/ |
/** MODREVERSE **/ |
/** **/ |
/** **/ |
/********************************************************************/ |
/********************************************************************/ |
|
/* evaluate f(x) mod T */ |
|
GEN |
|
RX_RXQ_compo(GEN f, GEN x, GEN T) |
|
{ |
|
gpmem_t av = avma; |
|
long l; |
|
GEN y; |
|
|
|
if (typ(f) != t_POL) return gcopy(f); |
|
l = lgef(f)-1; y = (GEN)f[l]; |
|
for (l--; l>=2; l--) |
|
y = gres(gadd(gmul(y,x), (GEN)f[l]), T); |
|
return gerepileupto(av, y); |
|
} |
|
|
|
/* return (1,...a^l) mod T. Not memory clean */ |
GEN |
GEN |
|
RXQ_powers(GEN a, GEN T, long l) |
|
{ |
|
long i; |
|
GEN v; |
|
|
|
if (typ(a) != t_POL) err(typeer,"RXQ_powers"); |
|
l += 2; |
|
v = cgetg(l,t_VEC); |
|
v[1] = un; if (l == 2) return v; |
|
|
|
if (degpol(a) >= degpol(T)) a = gres(a, T); |
|
v[2] = (long)a; |
|
for (i=3; i<l; i++) v[i] = lres(gmul((GEN)v[i-1], a), T); |
|
return v; |
|
} |
|
|
|
/* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */ |
|
GEN |
|
modreverse_i(GEN a, GEN T) |
|
{ |
|
gpmem_t av = avma; |
|
long n = degpol(T); |
|
GEN y; |
|
|
|
if (n <= 0) return gcopy(a); |
|
if (n == 1) |
|
return gerepileupto(av, gneg(gdiv((GEN)T[2], (GEN)T[3]))); |
|
if (gcmp0(a) || typ(a) != t_POL) err(talker,"reverse polmod does not exist"); |
|
|
|
y = vecpol_to_mat(RXQ_powers(a,T,n-1), n); |
|
y = gauss(y, _ei(n, 2)); |
|
return gerepilecopy(av, vec_to_pol(y, varn(T))); |
|
} |
|
|
|
GEN |
polymodrecip(GEN x) |
polymodrecip(GEN x) |
{ |
{ |
long v,i,j,n,av,tetpil,lx; |
long v, n; |
GEN p1,p2,p3,p,phi,y,col; |
GEN T, a, y; |
|
|
if (typ(x)!=t_POLMOD) err(talker,"not a polymod in polymodrecip"); |
if (typ(x)!=t_POLMOD) err(talker,"not a polmod in modreverse"); |
p=(GEN)x[1]; phi=(GEN)x[2]; |
T = (GEN)x[1]; |
v=varn(p); n=degpol(p); if (n<=0) return gcopy(x); |
a = (GEN)x[2]; |
if (n==1) |
n = degpol(T); if (n <= 0) return gcopy(x); |
{ |
v = varn(T); |
y=cgetg(3,t_POLMOD); |
y = cgetg(3,t_POLMOD); |
if (typ(phi)==t_POL) phi = (GEN)phi[2]; |
y[1] = (n==1)? lsub(polx[v], a): (long)caract2(T, a, v); |
p1=cgetg(4,t_POL); p1[1]=p[1]; p1[2]=lneg(phi); p1[3]=un; |
y[2] = (long)modreverse_i(a, T); return y; |
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)); |
|
} |
} |
|
|
/********************************************************************/ |
/********************************************************************/ |
Line 1119 longcmp(GEN x, GEN y) |
|
Line 1111 longcmp(GEN x, GEN y) |
|
|
|
/* Sort x = vector of elts, using cmp to compare them. |
/* Sort x = vector of elts, using cmp to compare them. |
* flag & cmp_IND: indirect sort: return permutation that would sort x |
* 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 |
* flag & cmp_C : as cmp_IND, but return permutation as vector of C-longs |
*/ |
*/ |
GEN |
GEN |
Line 1145 gen_sort(GEN x, int flag, int (*cmp)(GEN,GEN)) |
|
Line 1136 gen_sort(GEN x, int flag, int (*cmp)(GEN,GEN)) |
|
} |
} |
if (!cmp) cmp = &longcmp; |
if (!cmp) cmp = &longcmp; |
indx = (GEN) gpmalloc(lx*sizeof(long)); |
indx = (GEN) gpmalloc(lx*sizeof(long)); |
for (j=1; j<lx; j++) indx[j]=j; |
for (j=1; j<lx; j++) indx[j] = j; |
|
|
ir=lx-1; l=(ir>>1)+1; |
ir = lx-1; l = (ir>>1)+1; |
for(;;) |
for(;;) |
{ |
{ |
if (l>1) |
if (l > 1) indxt = indx[--l]; |
{ l--; indxt = indx[l]; } |
|
else |
else |
{ |
{ |
indxt = indx[ir]; indx[ir]=indx[1]; ir--; |
indxt = indx[ir]; indx[ir] = indx[1]; |
if (ir == 1) |
if (--ir == 1) break; /* DONE */ |
{ |
|
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; |
q = (GEN)x[indxt]; i = l; |
if (flag & cmp_REV) |
for (j=i<<1; j<=ir; j<<=1) |
for (j=i<<1; j<=ir; j<<=1) |
{ |
{ |
if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) < 0) j++; |
if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) > 0) j++; |
if (cmp(q,(GEN)x[indx[j]]) >= 0) break; |
if (cmp(q,(GEN)x[indx[j]]) <= 0) break; |
|
|
|
indx[i]=indx[j]; i=j; |
indx[i] = indx[j]; i = j; |
} |
} |
else |
indx[i] = indxt; |
for (j=i<<1; j<=ir; j<<=1) |
} |
{ |
indx[1] = indxt; |
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; |
if (flag & cmp_REV) |
} |
{ /* reverse order */ |
indx[i]=indxt; |
GEN indy = (GEN) gpmalloc(lx*sizeof(long)); |
|
for (j=1; j<lx; j++) indy[j] = indx[lx-j]; |
|
free(indx); indx = indy; |
} |
} |
|
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; |
} |
} |
|
|
#define sort_fun(flag) ((flag & cmp_LEX)? &lexcmp: &gcmp) |
#define sort_fun(flag) ((flag & cmp_LEX)? &lexcmp: &gcmp) |