version 1.1, 2001/10/02 11:17:06 |
version 1.2, 2002/09/11 07:26:53 |
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 rpowsi(ulong a, GEN n, long prec); |
|
extern GEN divrs2_safe(GEN x, long i); |
|
extern void dcxlog(double s, double t, double *a, double *b); |
|
extern double dnorm(double s, double t); |
|
extern GEN trans_fix_arg(long *prec, GEN *s0, GEN *sig, gpmem_t *av, GEN *res); |
|
|
GEN |
GEN |
cgetc(long l) |
cgetc(long l) |
Line 36 fix(GEN x, long l) |
|
Line 41 fix(GEN x, long l) |
|
y = cgetr(l); gaffect(x, y); return y; |
y = cgetr(l); gaffect(x, y); return y; |
} |
} |
|
|
|
long |
|
lgcx(GEN z) |
|
{ |
|
long tz = typ(z); |
|
|
|
switch(tz) |
|
{ |
|
case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: return BIGINT; |
|
case t_REAL: return lg(z); |
|
case t_COMPLEX: return min(lgcx((GEN)z[1]),lgcx((GEN)z[2])); |
|
default: err(typeer,"lgcx"); |
|
} |
|
return 0; |
|
} |
|
|
|
GEN |
|
setlgcx(GEN z, long l) |
|
{ |
|
long tz = typ(z); |
|
GEN p1; |
|
|
|
switch(tz) |
|
{ |
|
case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: return z; |
|
case t_REAL: p1 = cgetr(l); affrr(z,p1); return p1; |
|
case t_COMPLEX: p1 = cgetc(l); gaffect(z,p1); return p1; |
|
default: err(typeer,"setlgcx"); return gzero; /* not reached */ |
|
} |
|
} |
|
|
|
/* force z to be of type real/complex */ |
|
|
|
GEN |
|
setlgcx2(GEN z, long l) |
|
{ |
|
long tz = typ(z); |
|
GEN p1; |
|
|
|
switch(tz) |
|
{ |
|
case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: case t_REAL: |
|
p1 = cgetr(l); gaffect(z,p1); return p1; |
|
case t_COMPLEX: p1 = cgetc(l); gaffect(z,p1); return p1; |
|
default: err(typeer,"setlgcx"); return gzero; /* not reached */ |
|
} |
|
} |
|
|
|
/* a exporter ou ca existe deja ? */ |
|
|
|
long |
|
isint(GEN n, long *ptk) |
|
{ |
|
long tn=typ(n); |
|
GEN p1,p2; |
|
|
|
switch(tn) |
|
{ |
|
case t_INT: |
|
*ptk = itos(n); return 1; |
|
case t_REAL: |
|
p1 = gfloor(n); |
|
if (gcmp(p1,n)==0) {*ptk = itos(p1); return 1;} |
|
else return 0; |
|
case t_FRAC: case t_FRACN: |
|
p1 = dvmdii((GEN)n[1],(GEN)n[2],&p2); |
|
if (!signe(p2)) {*ptk = itos(p1); return 1;} |
|
else return 0; |
|
case t_COMPLEX: |
|
if (gcmp0((GEN)n[2])) return isint((GEN)n[1],ptk); |
|
else return 0; |
|
case t_QUAD: |
|
if (gcmp0((GEN)n[3])) return isint((GEN)n[2],ptk); |
|
else return 0; |
|
default: err(typeer,"isint"); return 0; /* not reached */ |
|
} |
|
} |
|
|
|
double |
|
norml1(GEN n, long prec) |
|
{ |
|
long tn=typ(n); |
|
gpmem_t av; |
|
double res; |
|
|
|
switch(tn) |
|
{ |
|
case t_INT: case t_REAL: case t_FRAC: case t_FRACN: |
|
return gtodouble(gabs(n,prec)); |
|
case t_COMPLEX: |
|
return norml1((GEN)n[1],prec)+norml1((GEN)n[2],prec); |
|
case t_QUAD: |
|
av = avma; res = norml1(gmul(n,realun(prec)),prec); avma = av; |
|
return res; |
|
default: err(typeer,"norml1"); return 0.; /* not reached */ |
|
} |
|
} |
|
|
/***********************************************************************/ |
/***********************************************************************/ |
/** **/ |
/** **/ |
/** K BESSEL FUNCTION **/ |
/** BESSEL FUNCTIONS **/ |
/** **/ |
/** **/ |
/***********************************************************************/ |
/***********************************************************************/ |
|
|
|
/* computes sum_{k=0}^m n!*((-1)^flag*z^2/4)^k/(k!*(k+n)!) */ |
|
|
|
static GEN |
|
_jbessel(GEN n, GEN z, long flag, long m) |
|
{ |
|
long k, limit; |
|
gpmem_t av; |
|
GEN p1,s; |
|
|
|
p1 = gmul2n(gsqr(z),-2); if (flag & 1) p1 = gneg(p1); |
|
if (typ(z) == t_SER) |
|
{ |
|
long v = valp(z); |
|
k = lg(p1)-2 - v; |
|
if (v < 0) err(negexper,"jbessel"); |
|
if (k <= 0) return gadd(gun, zeroser(varn(z), 2*v)); |
|
p1 = gprec(p1, k); |
|
} |
|
s = gun; |
|
av = avma; limit = stack_lim(av,1); |
|
for (k=m; k>=1; k--) |
|
{ |
|
s = gadd(gun,gdiv(gdivgs(gmul(p1,s),k),gaddsg(k,n))); |
|
if (low_stack(limit,stack_lim(av,1))) |
|
{ |
|
if (DEBUGMEM>1) err(warnmem,"jbessel"); |
|
s = gerepilecopy(av, s); |
|
} |
|
} |
|
return s; |
|
} |
|
|
|
static GEN |
|
jbesselintern(GEN n, GEN z, long flag, long prec) |
|
{ |
|
long tz=typ(z), i, lz, lim, k=-1, ki, precnew; |
|
gpmem_t av=avma, tetpil; |
|
double B,N,L,x; |
|
GEN p1,p2,y,znew,nnew; |
|
|
|
switch(tz) |
|
{ |
|
case t_REAL: case t_COMPLEX: |
|
i = precision(z); if (i) prec = i; |
|
if (isint(setlgcx2(n,prec),&ki)) |
|
{ |
|
k = labs(ki); |
|
p2 = gdiv(gpowgs(gmul2n(z,-1),k),mpfact(k)); |
|
if ((flag&1) && (ki<0) && (k&1)) p2 = gneg(p2); |
|
} |
|
else p2 = gdiv(gpow(gmul2n(z,-1),n,prec),ggamma(gaddgs(n,1),prec)); |
|
if (gcmp0(z)) return gerepilecopy(av, p2); |
|
else |
|
{ |
|
x = gtodouble(gabs(z,prec)); |
|
L = x*1.3591409; |
|
B = bit_accuracy(prec)*LOG2/(2*L); |
|
N = 1 + B; |
|
/* 3 Newton iterations are enough except in pathological cases */ |
|
N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); |
|
lim = max((long)(L*N),2); |
|
precnew = prec; |
|
if (x >= 1.0) precnew += 1 + (long)(x/(LOG2*BITS_IN_LONG)); |
|
znew = setlgcx(z,precnew); |
|
if (k >= 0) p1 = setlgcx(_jbessel(stoi(k),znew,flag,lim),prec); |
|
else |
|
{ |
|
i = precision(n); |
|
nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n; |
|
p1 = setlgcx(_jbessel(nnew,znew,flag,lim),prec); |
|
} |
|
tetpil = avma; return gerepile(av,tetpil,gmul(p2,p1)); |
|
} |
|
|
|
case t_SER: |
|
if (isint(setlgcx2(n,prec),&ki)) |
|
{ |
|
k = labs(ki); |
|
p1 = _jbessel(stoi(k),z,flag,lg(z)-2); |
|
} |
|
else p1 = _jbessel(n,z,flag,lg(z)-2); |
|
return gerepilecopy(av,p1); |
|
|
|
case t_VEC: case t_COL: case t_MAT: |
|
lz=lg(z); y=cgetg(lz,typ(z)); |
|
for (i=1; i<lz; i++) |
|
y[i]=(long)jbesselintern(n,(GEN)z[i],flag,prec); |
|
return y; |
|
|
|
case t_INT: case t_FRAC: case t_FRACN: |
|
av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselintern(n,p1,flag,prec)); |
|
|
|
case t_QUAD: |
|
av=avma; p1=gmul(z,realun(prec)); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselintern(n,p1,flag,prec)); |
|
|
|
case t_POL: case t_RFRAC: case t_RFRACN: |
|
av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselintern(n,p1,flag,prec)); |
|
|
|
case t_POLMOD: |
|
av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]); |
|
tetpil=avma; y=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) y[i]=(long)jbesselintern(n,(GEN)p2[i],flag,prec); |
|
return gerepile(av,tetpil,y); |
|
|
|
case t_PADIC: |
|
err(impl,"p-adic jbessel function"); |
|
} |
|
err(typeer,"jbessel"); |
|
return NULL; /* not reached */ |
|
} |
|
|
GEN |
GEN |
kbessel0(GEN nu, GEN gx, long flag, long prec) |
jbessel(GEN n, GEN z, long prec) |
{ |
{ |
switch(flag) |
return jbesselintern(n,z,1,prec); |
|
} |
|
|
|
GEN |
|
ibessel(GEN n, GEN z, long prec) |
|
{ |
|
return jbesselintern(n,z,0,prec); |
|
} |
|
|
|
GEN |
|
_jbesselh(long k, GEN z, long prec) |
|
{ |
|
GEN zinv,s,c,p0,p1,p2; |
|
long i; |
|
|
|
zinv=ginv(z); |
|
gsincos(z,&s,&c,prec); |
|
p1=gmul(zinv,s); |
|
if (k) |
{ |
{ |
case 0: return kbessel(nu,gx,prec); |
p0=p1; p1=gmul(zinv,gsub(p0,c)); |
case 1: return kbessel2(nu,gx,prec); |
for (i=2; i<=k; i++) |
default: err(flagerr,"besselk"); |
{ |
|
p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0); |
|
p0=p1; p1=p2; |
|
} |
} |
} |
|
return p1; |
|
} |
|
|
|
GEN |
|
jbesselh(GEN n, GEN z, long prec) |
|
{ |
|
long gz, k, l, linit, i, lz, tz=typ(z); |
|
gpmem_t av, tetpil; |
|
GEN y,p1,p2; |
|
|
|
if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh"); |
|
k = itos(n); |
|
if (k < 0) |
|
{ |
|
av = avma; n = gadd(ghalf,n); tetpil = avma; |
|
return gerepile(av,tetpil,jbessel(n,z,prec)); |
|
} |
|
|
|
switch(tz) |
|
{ |
|
case t_REAL: case t_COMPLEX: |
|
av = avma; |
|
if (gcmp0(z)) |
|
{ |
|
p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k)); |
|
p1 = gdiv(gmul(mpfact(k),p1),mpfact(2*k+1)); |
|
tetpil = avma; return gerepile(av,tetpil,gmul2n(p1,2*k)); |
|
} |
|
gz = gexpo(z); |
|
linit = lgcx(z); if (linit==BIGINT) linit = prec; |
|
if (gz>=0) l = linit; |
|
else l = lgcx(z) - 1 + ((-2*k*gz)>>TWOPOTBITS_IN_LONG); |
|
if (l>prec) prec = l; |
|
prec += (-gz)>>TWOPOTBITS_IN_LONG; |
|
z = setlgcx(z,prec); |
|
p1 = _jbesselh(k,z,prec); |
|
p1 = gmul(gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec),p1); |
|
tetpil = avma; return gerepile(av,tetpil,setlgcx(p1,linit)); |
|
|
|
case t_SER: |
|
if (gcmp0(z)) return gpowgs(z,k); |
|
av = avma; |
|
if (valp(z) < 0) err(negexper,"jbesselh"); |
|
z = gprec(z, lg(z)-2 + (2*k+1)*valp(z)); |
|
p1 = gdiv(_jbesselh(k,z,prec),gpowgs(z,k)); |
|
for (i=1; i<=k; i++) p1 = gmulgs(p1,2*i+1); |
|
return gerepilecopy(av,p1); |
|
|
|
case t_VEC: case t_COL: case t_MAT: |
|
lz=lg(z); y=cgetg(lz,typ(z)); |
|
for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)z[i],prec); |
|
return y; |
|
|
|
case t_INT: case t_FRAC: case t_FRACN: |
|
av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselh(n,p1,prec)); |
|
|
|
case t_QUAD: |
|
av=avma; p1=gmul(z,realun(prec)); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselh(n,p1,prec)); |
|
|
|
case t_POL: case t_RFRAC: case t_RFRACN: |
|
av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselh(n,p1,prec)); |
|
|
|
case t_POLMOD: |
|
av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]); |
|
tetpil=avma; y=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)p2[i],prec); |
|
return gerepile(av,tetpil,y); |
|
|
|
case t_PADIC: |
|
err(impl,"p-adic jbesselh function"); |
|
} |
|
err(typeer,"jbesselh"); |
return NULL; /* not reached */ |
return NULL; /* not reached */ |
} |
} |
|
|
GEN |
GEN |
kbessel(GEN nu, GEN gx, long prec) |
kbessel(GEN nu, GEN gx, long prec) |
{ |
{ |
GEN x,y,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w; |
GEN x,y,yfin,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w; |
long l,lbin,av,av1,k,k2,l1,n2,n; |
long l, lnew, lbin, k, k2, l1, n2, n, ex, rab=0; |
|
gpmem_t av, av1; |
|
|
if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec); |
if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec); |
l = (typ(gx)==t_REAL)? lg(gx): prec; |
l = (typ(gx)==t_REAL)? lg(gx): prec; |
y=cgetr(l); l1=l+1; |
ex = gexpo(gx); |
av=avma; x = fix(gx, l); |
if (ex < 0) |
|
{ |
|
rab = 1 + ((-ex)>>TWOPOTBITS_IN_LONG); |
|
lnew = l + rab; prec += rab; |
|
} |
|
else lnew = l; |
|
yfin=cgetr(l); l1=lnew+1; |
|
av=avma; x = fix(gx, lnew); y=cgetr(lnew); |
u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1); |
u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1); |
e=cgetr(l1); f=cgetr(l1); |
e=cgetr(l1); f=cgetr(l1); |
nu2=gmulgs(gsqr(nu),-4); |
nu2=gmulgs(gsqr(nu),-4); |
Line 88 kbessel(GEN nu, GEN gx, long prec) |
|
Line 410 kbessel(GEN nu, GEN gx, long prec) |
|
} |
} |
gmulz(s,zf,u); t=gmul2n(t,-1); |
gmulz(s,zf,u); t=gmul2n(t,-1); |
gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1; |
gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1; |
affsr(n2,q=cgetr(l1)); |
q = stor(n2, l1); |
r=gmul2n(x,1); av1=avma; |
r=gmul2n(x,1); av1=avma; |
for(;;) |
for(;;) |
{ |
{ |
Line 131 kbessel(GEN nu, GEN gx, long prec) |
|
Line 453 kbessel(GEN nu, GEN gx, long prec) |
|
} |
} |
gmulz(s,zf,y); |
gmulz(s,zf,y); |
} |
} |
gdivz(y,mpexp(x),y); avma=av; return y; |
gdivz(y,mpexp(x),yfin); |
|
avma=av; return yfin; |
} |
} |
|
|
|
/* computes sum_{k=0}^m ((-1)^flag*z^2/4)^k/(k!*(k+n)!)*(H(k)+H(k+n)) |
|
+sum_{k=0}^{n-1}((-1)^(flag+1)*z^2/4)^(k-n)*(n-k-1)!/k!. Ici n |
|
doit etre un long. Attention, contrairement a _jbessel, pas de n! devant. |
|
Quand flag > 1, calculer exactement les H(k) et factorielles */ |
|
|
|
static GEN |
|
_kbessel(long n, GEN z, long flag, long m, long prec) |
|
{ |
|
long k, limit; |
|
gpmem_t av; |
|
GEN p1,p2,p3,s,*tabh; |
|
|
|
p1 = gmul2n(gsqr(z),-2); if (flag & 1) p1 = gneg(p1); |
|
if (typ(z) == t_SER) |
|
{ |
|
long v = valp(z); |
|
k = lg(p1)-2 - v; |
|
if (v < 0) err(negexper,"kbessel"); |
|
if (k <= 0) return gadd(gun, zeroser(varn(z), 2*v)); |
|
p1 = gprec(p1, k); |
|
} |
|
tabh = (GEN*)cgetg(m+n+2,t_VEC); tabh[1] = gzero; |
|
if (flag <= 1) |
|
{ |
|
s = realun(prec); tabh[2] = s; |
|
for (k=2; k<=m+n; k++) |
|
{ |
|
s = divrs(addsr(1,mulsr(k,s)),k); tabh[k+1] = s; |
|
} |
|
} |
|
else |
|
{ |
|
s = gun; tabh[2] = s; |
|
for (k=2; k<=m+n; k++) |
|
{ |
|
s = gdivgs(gaddsg(1,gmulsg(k,s)),k); tabh[k+1] = s; |
|
} |
|
} |
|
s = gadd(tabh[m+1],tabh[m+n+1]); |
|
av = avma; limit = stack_lim(av,1); |
|
for (k=m; k>=1; k--) |
|
{ |
|
s = gadd(gadd(tabh[k],tabh[k+n]),gdiv(gmul(p1,s),mulss(k,k+n))); |
|
if (low_stack(limit,stack_lim(av,1))) |
|
{ |
|
if (DEBUGMEM>1) err(warnmem,"kbessel"); |
|
s = gerepilecopy(av, s); |
|
} |
|
} |
|
p3 = (flag <= 1) ? mpfactr(n,prec) : mpfact(n); |
|
s = gdiv(s,p3); |
|
if (n) |
|
{ |
|
p1 = gneg(ginv(p1)); |
|
p2 = gmulsg(n,gdiv(p1,p3)); |
|
for (k=n-1; k>=0; k--) |
|
{ |
|
s = gadd(s,p2); |
|
p2 = gmul(p2,gmul(mulss(k,n-k),p1)); |
|
} |
|
} |
|
return s; |
|
} |
|
|
|
static GEN |
|
kbesselintern(GEN n, GEN z, long flag, long prec) |
|
{ |
|
long tz=typ(z), i, k, ki, lz, lim, precnew, fl, fl2, ex; |
|
gpmem_t av=avma, tetpil; |
|
double B,N,L,x,rab; |
|
GEN p1,p2,y,p3,znew,nnew,pplus,pmoins,s,c; |
|
|
|
fl = (flag & 1) == 0; |
|
switch(tz) |
|
{ |
|
case t_REAL: case t_COMPLEX: |
|
if (gcmp0(z)) err(talker,"zero argument in a k/n bessel function"); |
|
i = precision(z); if (i) prec = i; |
|
x = gtodouble(gabs(z,prec)); |
|
/* Experimentally optimal on a PIII 500 Mhz. Your optimum may vary. */ |
|
if ((x > (8*(prec-2)+norml1(n,prec)-3)) && !flag) return kbessel(n,z,prec); |
|
precnew = prec; |
|
if (x >= 1.0) |
|
{ |
|
rab = x/(LOG2*BITS_IN_LONG); if (fl) rab *= 2; |
|
precnew += 1 + (long)rab; |
|
} |
|
znew = setlgcx(z,precnew); |
|
if (isint(setlgcx2(n,precnew),&ki)) |
|
{ |
|
k = labs(ki); |
|
L = x*1.3591409; |
|
B = (bit_accuracy(prec)*LOG2)/(2*L); |
|
if (fl) B += 0.367879; |
|
N = 1 + B; |
|
/* 3 Newton iterations are enough except in pathological cases */ |
|
N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); |
|
lim = max((long)(L*N),2); |
|
p1 = _kbessel(k,znew,flag,lim,precnew); |
|
p1 = gmul(gpowgs(gmul2n(znew,-1),k),p1); |
|
p2 = gadd(mpeuler(precnew),glog(gmul2n(znew,-1),precnew)); |
|
p3 = jbesselintern(stoi(k),znew,flag,precnew); |
|
p2 = gsub(gmul2n(p1,-1),gmul(p2,p3)); |
|
p2 = setlgcx(p2,prec); |
|
if (fl == 0) {p1 = mppi(prec); p2 = gmul2n(gdiv(p2,p1),1);} |
|
fl = (fl && (k&1)) || (!fl && (ki>=0 || (k&1)==0)); |
|
tetpil = avma; return gerepile(av,tetpil,fl ? gneg(p2) : gcopy(p2)); |
|
} |
|
else |
|
{ |
|
i = precision(n); |
|
nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n; |
|
p2 = mppi(precnew); gsincos(gmul(nnew,p2),&s,&c,precnew); |
|
ex = gexpo(s); |
|
if (ex < 0) |
|
{ |
|
rab = (-ex)/(LOG2*BITS_IN_LONG); if (fl) rab *= 2; |
|
precnew += 1 + (long)rab; |
|
} |
|
nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n; |
|
znew = setlgcx(znew,precnew); |
|
p2 = mppi(precnew); gsincos(gmul(nnew,p2),&s,&c,precnew); |
|
pplus = jbesselintern(nnew,znew,flag,precnew); |
|
pmoins = jbesselintern(gneg(nnew),znew,flag,precnew); |
|
if (fl) p1 = gmul(gsub(pmoins,pplus),gdiv(p2,gmul2n(s,1))); |
|
else p1 = gdiv(gsub(gmul(c,pplus),pmoins),s); |
|
tetpil = avma; return gerepile(av,tetpil,setlgcx(p1,prec)); |
|
} |
|
|
|
case t_SER: |
|
if (isint(setlgcx2(n,prec),&ki)) |
|
{ |
|
k = labs(ki); |
|
p1 = _kbessel(k,z,flag+2,lg(z)-2,prec); |
|
return gerepilecopy(av,p1); |
|
} |
|
else |
|
{ |
|
fl2 = isint(setlgcx2(gmul2n(n,1),prec),&ki); |
|
if (!fl2) |
|
err(talker,"cannot give a power series result in k/n bessel function"); |
|
else |
|
{ |
|
k = labs(ki); n = gmul2n(stoi(k),-1); |
|
fl2 = (k&3)==1; |
|
pmoins = jbesselintern(gneg(n),z,flag,prec); |
|
if (fl) |
|
{ |
|
pplus = jbesselintern(n,z,flag,prec); |
|
p2 = gpowgs(z,-k); if (fl2 == 0) p2 = gneg(p2); |
|
p3 = gmul2n(divii(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1)); |
|
p3 = gdivgs(gmul2n(gsqr(p3),1),k); |
|
p2 = gmul(p2,p3); |
|
p1 = gsub(pplus,gmul(p2,pmoins)); |
|
} |
|
else p1 = pmoins; |
|
tetpil = avma; |
|
return gerepile(av,tetpil,fl2 ? gneg(p1) : gcopy(p1)); |
|
} |
|
} |
|
|
|
case t_VEC: case t_COL: case t_MAT: |
|
lz=lg(z); y=cgetg(lz,typ(z)); |
|
for (i=1; i<lz; i++) |
|
y[i]=(long)kbesselintern(n,(GEN)z[i],flag,prec); |
|
return y; |
|
|
|
case t_INT: case t_FRAC: case t_FRACN: |
|
av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma; |
|
return gerepile(av,tetpil,kbesselintern(n,p1,flag,prec)); |
|
|
|
case t_QUAD: |
|
av=avma; p1=gmul(z,realun(prec)); tetpil=avma; |
|
return gerepile(av,tetpil,kbesselintern(n,p1,flag,prec)); |
|
|
|
case t_POL: case t_RFRAC: case t_RFRACN: |
|
av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma; |
|
return gerepile(av,tetpil,kbesselintern(n,p1,flag,prec)); |
|
|
|
case t_POLMOD: |
|
av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]); |
|
tetpil=avma; y=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) y[i]=(long)kbesselintern(n,(GEN)p2[i],flag,prec); |
|
return gerepile(av,tetpil,y); |
|
|
|
case t_PADIC: |
|
err(impl,"p-adic kbessel function"); |
|
} |
|
err(typeer,"kbessel"); |
|
return NULL; /* not reached */ |
|
} |
|
|
|
GEN |
|
kbesselnew(GEN n, GEN z, long prec) |
|
{ |
|
return kbesselintern(n,z,0,prec); |
|
} |
|
|
|
GEN |
|
kbesselnewalways(GEN n, GEN z, long prec) |
|
{ |
|
return kbesselintern(n,z,2,prec); |
|
} |
|
|
|
GEN |
|
nbessel(GEN n, GEN z, long prec) |
|
{ |
|
return kbesselintern(n,z,1,prec); |
|
} |
|
|
|
GEN |
|
hbessel1(GEN n, GEN z, long prec) |
|
{ |
|
gpmem_t av=avma, tetpil; |
|
GEN p1,p2; |
|
|
|
p1 = jbessel(n,z,prec); p2 = gmul(gi,nbessel(n,z,prec)); |
|
tetpil = avma; return gerepile(av,tetpil,gadd(p1,p2)); |
|
} |
|
|
|
GEN |
|
hbessel2(GEN n, GEN z, long prec) |
|
{ |
|
gpmem_t av=avma, tetpil; |
|
GEN p1,p2; |
|
|
|
p1 = jbessel(n,z,prec); p2 = gmul(gi,nbessel(n,z,prec)); |
|
tetpil = avma; return gerepile(av,tetpil,gsub(p1,p2)); |
|
} |
|
|
|
GEN kbessel2(GEN nu, GEN x, long prec); |
|
|
|
GEN |
|
kbessel0(GEN nu, GEN gx, long flag, long prec) |
|
{ |
|
switch(flag) |
|
{ |
|
case 0: return kbesselnew(nu,gx,prec); |
|
case 1: return kbessel(nu,gx,prec); |
|
case 2: return kbessel2(nu,gx,prec); |
|
case 3: return kbesselnewalways(nu,gx,prec); |
|
default: err(flagerr,"besselk"); |
|
} |
|
return NULL; /* not reached */ |
|
} |
|
|
/***********************************************************************/ |
/***********************************************************************/ |
/* **/ |
/* **/ |
/** FONCTION U(a,b,z) GENERALE **/ |
/** FONCTION U(a,b,z) GENERALE **/ |
|
|
hyperu(GEN a, GEN b, GEN gx, long prec) |
hyperu(GEN a, GEN b, GEN gx, long prec) |
{ |
{ |
GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn; |
GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn; |
long l,lbin,av,av1,av2,k,l1,n,ex; |
long l, lbin, k, l1, n, ex; |
|
gpmem_t av, av1, av2; |
|
|
if(gsigne(gx) <= 0) err(talker,"hyperu's third argument must be positive"); |
if(gsigne(gx) <= 0) err(talker,"hyperu's third argument must be positive"); |
ex = (iscomplex(a) || iscomplex(b)); |
ex = (iscomplex(a) || iscomplex(b)); |
Line 181 hyperu(GEN a, GEN b, GEN gx, long prec) |
|
Line 752 hyperu(GEN a, GEN b, GEN gx, long prec) |
|
t=gadd(gaddsg(k,a),gmul(p1,t)); |
t=gadd(gaddsg(k,a),gmul(p1,t)); |
} |
} |
gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1; |
gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1; |
q=cgetr(l1); affsr(n,q); r=x; av1=avma; |
q = stor(n, l1); r=x; av1=avma; |
do |
do |
{ |
{ |
p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf; |
p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf; |
Line 225 hyperu(GEN a, GEN b, GEN gx, long prec) |
|
Line 796 hyperu(GEN a, GEN b, GEN gx, long prec) |
|
GEN |
GEN |
kbessel2(GEN nu, GEN x, long prec) |
kbessel2(GEN nu, GEN x, long prec) |
{ |
{ |
long av = avma,tetpil; |
gpmem_t av = avma, tetpil; |
GEN p1,p2,x2,a,pitemp; |
GEN p1,p2,x2,a,pitemp; |
|
|
if (typ(x)==t_REAL) prec = lg(x); |
if (typ(x)==t_REAL) prec = lg(x); |
|
|
incgam(GEN s, GEN x, long prec) |
incgam(GEN s, GEN x, long prec) |
{ |
{ |
GEN p1,z = cgetr(prec); |
GEN p1,z = cgetr(prec); |
long av = avma; |
gpmem_t av = avma; |
|
|
if (gcmp0(x)) return ggamma(s,prec); |
if (gcmp0(x)) return ggamma(s,prec); |
if (typ(x)!=t_REAL) { gaffect(x,z); x=z; } |
if (typ(x)!=t_REAL) { gaffect(x,z); x=z; } |
|
|
incgam1(GEN a, GEN x, long prec) |
incgam1(GEN a, GEN x, long prec) |
{ |
{ |
GEN p2,p3,y, z = cgetr(prec); |
GEN p2,p3,y, z = cgetr(prec); |
long av=avma, av1, l,n,i; |
long l, n, i; |
|
gpmem_t av=avma, av1; |
double m,mx; |
double m,mx; |
|
|
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
|
|
incgam2(GEN a, GEN x, long prec) |
incgam2(GEN a, GEN x, long prec) |
{ |
{ |
GEN b,p1,p2,p3,y, z = cgetr(prec); |
GEN b,p1,p2,p3,y, z = cgetr(prec); |
long av = avma,av1,l,n,i; |
long l, n, i; |
|
gpmem_t av = avma, av1; |
double m,mx; |
double m,mx; |
|
|
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
|
|
incgam3(GEN s, GEN x, long prec) |
incgam3(GEN s, GEN x, long prec) |
{ |
{ |
GEN b,p1,p2,p3,y, z = cgetr(prec); |
GEN b,p1,p2,p3,y, z = cgetr(prec); |
long av = avma, av1,l,n,i; |
long l, n, i; |
|
gpmem_t av = avma, av1; |
|
|
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
l=lg(x); n = -bit_accuracy(l)-1; |
l=lg(x); n = -bit_accuracy(l)-1; |
|
|
incgam4(GEN a, GEN x, GEN g, long prec) |
incgam4(GEN a, GEN x, GEN g, long prec) |
{ |
{ |
GEN p1, z = cgetr(prec); |
GEN p1, z = cgetr(prec); |
long av = avma; |
gpmem_t av = avma; |
|
|
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
if (typ(x) != t_REAL) { gaffect(x,z); x=z; } |
if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0) |
if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0) |
Line 384 incgam0(GEN a, GEN x, GEN z,long prec) |
|
Line 958 incgam0(GEN a, GEN x, GEN z,long prec) |
|
GEN |
GEN |
eint1(GEN x, long prec) |
eint1(GEN x, long prec) |
{ |
{ |
long av = avma,tetpil,l,n,i; |
long l, n, i; |
|
gpmem_t av = avma, tetpil; |
GEN p1,p2,p3,p4,run,y; |
GEN p1,p2,p3,p4,run,y; |
|
|
if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;} |
if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;} |
Line 454 eint1(GEN x, long prec) |
|
Line 1029 eint1(GEN x, long prec) |
|
GEN |
GEN |
veceint1(GEN C, GEN nmax, long prec) |
veceint1(GEN C, GEN nmax, long prec) |
{ |
{ |
long av,av1,k,n,nstop,i,cd,nmin,G,a, chkpoint; |
long k, n, nstop, i, cd, nmin, G, a, chkpoint; |
|
gpmem_t av, av1; |
GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit; |
GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit; |
|
|
if (!nmax) return eint1(C,prec); |
if (!nmax) return eint1(C,prec); |
Line 527 veceint1(GEN C, GEN nmax, long prec) |
|
Line 1103 veceint1(GEN C, GEN nmax, long prec) |
|
GEN |
GEN |
gerfc(GEN x, long prec) |
gerfc(GEN x, long prec) |
{ |
{ |
long av=avma; |
gpmem_t av=avma; |
GEN p1,p2; |
GEN p1,p2; |
|
|
if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; } |
if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; } |
Line 548 gerfc(GEN x, long prec) |
|
Line 1124 gerfc(GEN x, long prec) |
|
static GEN |
static GEN |
czeta(GEN s, long prec) |
czeta(GEN s, long prec) |
{ |
{ |
long av,n,p,n1,flag1,flag2,i,i2; |
long n, p, n1, flag1, flag2, i, i2; |
|
gpmem_t av; |
double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf; |
double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf; |
GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp; |
GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp; |
|
|
Line 677 next_bin(GEN y, long n, long k) |
|
Line 1254 next_bin(GEN y, long n, long k) |
|
static GEN |
static GEN |
izeta(long k, long prec) |
izeta(long k, long prec) |
{ |
{ |
long av=avma,av2,kk,n,li, limit; |
long kk, n, li; |
|
gpmem_t av=avma, av2, limit; |
GEN y,p1,pi2,qn,z,q,binom; |
GEN y,p1,pi2,qn,z,q,binom; |
|
|
/* treat trivial cases */ |
/* treat trivial cases */ |
Line 689 izeta(long k, long prec) |
|
Line 1267 izeta(long k, long prec) |
|
return gerepileuptoleaf(av, divrs(y,k-1)); |
return gerepileuptoleaf(av, divrs(y,k-1)); |
} |
} |
if (k > bit_accuracy(prec)+1) return realun(prec); |
if (k > bit_accuracy(prec)+1) return realun(prec); |
pi2 = mppi(prec); setexpo(pi2,2); /* 2Pi */ |
pi2 = Pi2n(1, prec); |
if ((k&1) == 0) |
if ((k&1) == 0) |
{ |
{ |
p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec))); |
p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec))); |
Line 715 izeta(long k, long prec) |
|
Line 1293 izeta(long k, long prec) |
|
y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y); |
y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y); |
|
|
av2 = avma; limit = stack_lim(av2,1); |
av2 = avma; limit = stack_lim(av2,1); |
qn = gsqr(q); z = divsr(1,addrs(q,-1)); |
qn = gsqr(q); z = ginv( addrs(q,-1) ); |
for (n=2; ; n++) |
for (n=2; ; n++) |
{ |
{ |
p1 = divsr(1,mulir(gpuigs(stoi(n),k),addrs(qn,-1))); |
p1 = ginv( mulir(gpuigs(stoi(n),k),addrs(qn,-1)) ); |
|
|
z = addrr(z,p1); if (expo(p1)< li) break; |
z = addrr(z,p1); if (expo(p1)< li) break; |
qn = mulrr(qn,q); |
qn = mulrr(qn,q); |
if (low_stack(limit,stack_lim(av2,1))) |
if (low_stack(limit,stack_lim(av2,1))) |
{ |
{ |
GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn; |
|
if (DEBUGMEM>1) err(warnmem,"izeta"); |
if (DEBUGMEM>1) err(warnmem,"izeta"); |
gerepilemany(av2,gptr,2); |
gerepileall(av2,2, &z, &qn); |
} |
} |
} |
} |
setexpo(z,expo(z)+1); |
setexpo(z,expo(z)+1); |
Line 757 izeta(long k, long prec) |
|
Line 1334 izeta(long k, long prec) |
|
qn=mulrr(qn,q); |
qn=mulrr(qn,q); |
if (low_stack(limit,stack_lim(av2,1))) |
if (low_stack(limit,stack_lim(av2,1))) |
{ |
{ |
GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn; |
|
if (DEBUGMEM>1) err(warnmem,"izeta"); |
if (DEBUGMEM>1) err(warnmem,"izeta"); |
gerepilemany(av2,gptr,2); |
gerepileall(av2,2, &z, &qn); |
} |
} |
} |
} |
setexpo(z,expo(z)+1); |
setexpo(z,expo(z)+1); |
Line 780 pows(long x, long n) |
|
Line 1356 pows(long x, long n) |
|
|
|
/* return n^-s, n > 1 odd. tab[q] := q^-s, q prime power */ |
/* return n^-s, n > 1 odd. tab[q] := q^-s, q prime power */ |
static GEN |
static GEN |
n_s(long n, GEN *tab) |
n_s(ulong n, GEN *tab) |
{ |
{ |
byteptr d = diffptr + 2; |
byteptr d = diffptr + 2; |
GEN x = NULL; |
GEN x = NULL; |
long p, e; |
long p, e; |
|
|
for (p = 3; n > 1; p += *d++) |
for (p = 3; n > 1; ) |
{ |
{ |
e = svaluation(n, p, &n); |
e = svaluation(n, p, &n); |
if (e) |
if (e) |
Line 794 n_s(long n, GEN *tab) |
|
Line 1370 n_s(long n, GEN *tab) |
|
GEN y = tab[pows(p,e)]; |
GEN y = tab[pows(p,e)]; |
if (!x) x = y; else x = gmul(x,y); |
if (!x) x = y; else x = gmul(x,y); |
} |
} |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
} |
} |
return x; |
return x; |
} |
} |
Line 804 GEN czeta(GEN s0, long prec); |
|
Line 1381 GEN czeta(GEN s0, long prec); |
|
static GEN |
static GEN |
izeta(long k, long prec) |
izeta(long k, long prec) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN y,p1,pi2; |
GEN y,p1,pi2; |
|
|
/* treat trivial cases */ |
/* treat trivial cases */ |
Line 827 izeta(long k, long prec) |
|
Line 1404 izeta(long k, long prec) |
|
return czeta(stoi(k), prec); |
return czeta(stoi(k), prec); |
} |
} |
|
|
extern GEN rpowsi(ulong a, GEN n, long prec); |
|
extern GEN divrs2_safe(GEN x, long i); |
|
extern void dcxlog(double s, double t, double *a, double *b); |
|
extern double dnorm(double s, double t); |
|
extern GEN trans_fix_arg(long *prec, GEN *s0, GEN *sig, long *av, GEN *res); |
|
|
|
/* s0 a t_INT, t_REAL or t_COMPLEX. |
/* s0 a t_INT, t_REAL or t_COMPLEX. |
* If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) |
* If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) |
* */ |
* */ |
Line 841 czeta(GEN s0, long prec) |
|
Line 1412 czeta(GEN s0, long prec) |
|
{ |
{ |
GEN s, u, a, y, res, tes, sig, invn2, p1, unr; |
GEN s, u, a, y, res, tes, sig, invn2, p1, unr; |
GEN sim, ms, s1, s2, s3, s4, s5, *tab, tabn; |
GEN sim, ms, s1, s2, s3, s4, s5, *tab, tabn; |
long p, i, sqn, nn, lim, lim2, ct, av, av2 = avma, avlim; |
long p, i, sqn, nn, lim, lim2, ct; |
|
gpmem_t av, av2 = avma, avlim; |
int funeq = 0; |
int funeq = 0; |
byteptr d; |
byteptr d; |
|
|
if (DEBUGLEVEL) timer2(); |
if (DEBUGLEVEL>2) (void)timer2(); |
s = trans_fix_arg(&prec,&s0,&sig,&av,&res); |
s = trans_fix_arg(&prec,&s0,&sig,&av,&res); |
if (gcmp0(s)) { y = gneg(ghalf); goto END; } |
if (gcmp0(s)) { y = gneg(ghalf); goto END; } |
if (signe(sig) <= 0 || expo(sig) < -1) |
if (signe(sig) <= 0 || expo(sig) < -1) |
Line 888 czeta(GEN s0, long prec) |
|
Line 1460 czeta(GEN s0, long prec) |
|
lim = (long) ceil(l); if (lim < 2) lim = 2; |
lim = (long) ceil(l); if (lim < 2) lim = 2; |
l2 = (lim+ssig/2.-.25); |
l2 = (lim+ssig/2.-.25); |
nn = 1 + (long)ceil( sqrt(l2*l2 + st*st/4) * la / PI ); |
nn = 1 + (long)ceil( sqrt(l2*l2 + st*st/4) * la / PI ); |
if (DEBUGLEVEL) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn); |
if (DEBUGLEVEL>2) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn); |
if ((ulong)nn >= maxprime()) err(primer1); |
if ((ulong)nn >= maxprime()) err(primer1); |
} |
} |
prec++; unr = realun(prec); /* one extra word of precision */ |
prec++; unr = realun(prec); /* one extra word of precision */ |
Line 897 czeta(GEN s0, long prec) |
|
Line 1469 czeta(GEN s0, long prec) |
|
d = diffptr + 1; |
d = diffptr + 1; |
if (typ(s0) == t_INT) |
if (typ(s0) == t_INT) |
{ /* no explog for 1/p^s */ |
{ /* no explog for 1/p^s */ |
for (p=2; p < nn; p += *d++) |
for (p=2; p < nn;) { |
tab[p] = divrr(unr, rpowsi(p, s0, prec)); |
tab[p] = divrr(unr, rpowsi(p, s0, prec)); |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
|
} |
a = divrr(unr, rpowsi(nn, s0, prec)); |
a = divrr(unr, rpowsi(nn, s0, prec)); |
} |
} |
else |
else |
{ /* general case */ |
{ /* general case */ |
ms = gneg(s); p1 = cgetr(prec); |
ms = gneg(s); p1 = cgetr(prec); |
for (p=2; p < nn; p += *d++) |
for (p=2; p < nn;) |
{ |
{ |
affsr(p, p1); |
affsr(p, p1); |
tab[p] = gexp(gmul(ms, mplog(p1)), prec); |
tab[p] = gexp(gmul(ms, mplog(p1)), prec); |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
} |
} |
affsr(nn,p1); |
affsr(nn,p1); |
a = gexp(gmul(ms, mplog(p1)), prec); |
a = gexp(gmul(ms, mplog(p1)), prec); |
} |
} |
sqn = (long)sqrt(nn-1.); |
sqn = (long)sqrt(nn-1.); |
d = diffptr + 2; /* fill in odd prime powers */ |
d = diffptr + 2; /* fill in odd prime powers */ |
for (p=3; p <= sqn; p += *d++) |
for (p=3; p <= sqn; ) |
{ |
{ |
ulong oldq = p, q = p*p; |
ulong oldq = p, q = p*p; |
while (q<(ulong)nn) { tab[q] = gmul(tab[p], tab[oldq]); oldq = q; q *= p; } |
while (q<(ulong)nn) { tab[q] = gmul(tab[p], tab[oldq]); oldq = q; q *= p; } |
|
NEXT_PRIME_VIADIFF_CHECK(p,d); |
} |
} |
if (DEBUGLEVEL) msgtimer("tab[q^-s] from 1 to N-1"); |
if (DEBUGLEVEL>2) msgtimer("tab[q^-s] from 1 to N-1"); |
|
|
tabn = cgetg(nn, t_VECSMALL); ct = 0; |
tabn = cgetg(nn, t_VECSMALL); ct = 0; |
for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1; |
for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1; |
sim = y = unr; |
sim = y = unr; |
for (i=ct; i > 1; i--) |
for (i=ct; i > 1; i--) |
{ |
{ |
long j, av2 = avma; |
long j; |
|
gpmem_t av2 = avma; |
for (j=tabn[i]+1; j<=tabn[i-1]; j++) |
for (j=tabn[i]+1; j<=tabn[i-1]; j++) |
sim = gadd(sim, n_s(2*j+1, tab)); |
sim = gadd(sim, n_s(2*j+1, tab)); |
sim = gerepileupto(av2, sim); |
sim = gerepileupto(av2, sim); |
y = gadd(sim, gmul(tab[2],y)); |
y = gadd(sim, gmul(tab[2],y)); |
} |
} |
y = gadd(y, gmul2n(a,-1)); |
y = gadd(y, gmul2n(a,-1)); |
if (DEBUGLEVEL) msgtimer("sum from 1 to N-1"); |
if (DEBUGLEVEL>2) msgtimer("sum from 1 to N-1"); |
|
|
invn2 = divrs(unr, nn*nn); lim2 = lim<<1; |
invn2 = divrs(unr, nn*nn); lim2 = lim<<1; |
tes = bernreal(lim2, prec); |
tes = bernreal(lim2, prec); |
Line 976 czeta(GEN s0, long prec) |
|
Line 1553 czeta(GEN s0, long prec) |
|
u = gmul(gmul(tes,invn2), gmul2n(s2, -1)); |
u = gmul(gmul(tes,invn2), gmul2n(s2, -1)); |
tes = gmulsg(nn, gaddsg(1, u)); |
tes = gmulsg(nn, gaddsg(1, u)); |
} |
} |
if (DEBUGLEVEL) msgtimer("Bernoulli sum"); |
if (DEBUGLEVEL>2) msgtimer("Bernoulli sum"); |
/* y += tes a / (s-1) */ |
/* y += tes a / (s-1) */ |
y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr)))); |
y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr)))); |
|
|
|
|
if (funeq) |
if (funeq) |
{ |
{ |
GEN pitemp = mppi(prec); setexpo(pitemp,2); /* 2Pi */ |
GEN pitemp = mppi(prec); setexpo(pitemp,2); /* 2Pi */ |
y = gmul(gmul(y, ggamma(s,prec)), gpui(pitemp,gneg(s),prec)); |
y = gmul(gmul(y, ggamma(gprec_w(s,prec-1),prec)), |
|
gpui(pitemp,gneg(s),prec)); |
setexpo(pitemp,0); /* Pi/2 */ |
setexpo(pitemp,0); /* Pi/2 */ |
y = gmul2n(gmul(y, gcos(gmul(pitemp,s),prec)), 1); |
y = gmul2n(gmul(y, gcos(gmul(pitemp,s),prec)), 1); |
} |
} |
Line 1020 gzeta(GEN x, long prec) |
|
Line 1598 gzeta(GEN x, long prec) |
|
void |
void |
gzetaz(GEN x, GEN y) |
gzetaz(GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gzetaz"); |
if (!prec) err(infprecer,"gzetaz"); |
gaffect(gzeta(x,prec),y); avma=av; |
gaffect(gzeta(x,prec),y); avma=av; |
Line 1038 gzetaz(GEN x, GEN y) |
|
Line 1617 gzetaz(GEN x, GEN y) |
|
static GEN |
static GEN |
cxpolylog(long m, GEN x, long prec) |
cxpolylog(long m, GEN x, long prec) |
{ |
{ |
long av=avma,li,i,n,bern_upto; |
long li, i, n, bern_upto; |
|
gpmem_t av=avma; |
GEN p1,z,h,q,s; |
GEN p1,z,h,q,s; |
|
|
if (gcmp1(x)) return izeta(m,prec); |
if (gcmp1(x)) return izeta(m,prec); |
Line 1070 cxpolylog(long m, GEN x, long prec) |
|
Line 1650 cxpolylog(long m, GEN x, long prec) |
|
GEN |
GEN |
polylog(long m, GEN x, long prec) |
polylog(long m, GEN x, long prec) |
{ |
{ |
long av,av1,limpile,l,e,i,G,sx; |
long l, e, i, G, sx; |
|
gpmem_t av, av1, limpile; |
GEN X,z,p1,p2,n,y,logx; |
GEN X,z,p1,p2,n,y,logx; |
|
|
if (m<0) err(talker,"negative index in polylog"); |
if (m<0) err(talker,"negative index in polylog"); |
Line 1142 dilog(GEN x, long prec) |
|
Line 1723 dilog(GEN x, long prec) |
|
GEN |
GEN |
polylogd0(long m, GEN x, long flag, long prec) |
polylogd0(long m, GEN x, long flag, long prec) |
{ |
{ |
long k,l,av,fl,m2; |
long k, l, fl, m2; |
|
gpmem_t av; |
GEN p1,p2,p3,y; |
GEN p1,p2,p3,y; |
|
|
m2=m&1; av=avma; |
m2=m&1; av=avma; |
Line 1188 polylogdold(long m, GEN x, long prec) |
|
Line 1770 polylogdold(long m, GEN x, long prec) |
|
GEN |
GEN |
polylogp(long m, GEN x, long prec) |
polylogp(long m, GEN x, long prec) |
{ |
{ |
long k,l,av,fl,m2; |
long k, l, fl, m2; |
|
gpmem_t av; |
GEN p1,y; |
GEN p1,y; |
|
|
m2=m&1; av=avma; |
m2=m&1; av=avma; |
Line 1234 polylogp(long m, GEN x, long prec) |
|
Line 1817 polylogp(long m, GEN x, long prec) |
|
GEN |
GEN |
gpolylog(long m, GEN x, long prec) |
gpolylog(long m, GEN x, long prec) |
{ |
{ |
long i,lx,av=avma,v,n; |
long i, lx, v, n; |
|
gpmem_t av=avma; |
GEN y,p1,p2; |
GEN y,p1,p2; |
|
|
if (m<=0) |
if (m<=0) |
Line 1272 gpolylog(long m, GEN x, long prec) |
|
Line 1856 gpolylog(long m, GEN x, long prec) |
|
p1=glog(gsub(gun,x),prec); |
p1=glog(gsub(gun,x),prec); |
return gerepileupto(av, gneg(p1)); |
return gerepileupto(av, gneg(p1)); |
} |
} |
|
if (gcmp0(x)) return gcopy(x); |
if (valp(x)<=0) err(impl,"polylog around a!=0"); |
if (valp(x)<=0) err(impl,"polylog around a!=0"); |
v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2); |
v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2); |
for (i=n; i>=1; i--) |
for (i=n; i>=1; i--) |
Line 1294 gpolylog(long m, GEN x, long prec) |
|
Line 1879 gpolylog(long m, GEN x, long prec) |
|
void |
void |
gpolylogz(long m, GEN x, GEN y) |
gpolylogz(long m, GEN x, GEN y) |
{ |
{ |
long av=avma, prec = precision(y); |
long prec = precision(y); |
|
gpmem_t av=avma; |
|
|
if (!prec) err(infprecer,"gpolylogz"); |
if (!prec) err(infprecer,"gpolylogz"); |
gaffect(gpolylog(m,x,prec),y); avma=av; |
gaffect(gpolylog(m,x,prec),y); avma=av; |
Line 1322 qq(GEN x, long prec) |
|
Line 1908 qq(GEN x, long prec) |
|
if (tx==t_PADIC) return x; |
if (tx==t_PADIC) return x; |
if (is_scalar_t(tx)) |
if (is_scalar_t(tx)) |
{ |
{ |
long l=precision(x); |
long l = precision(x); |
GEN p1,p2,q; |
if (!l) l = prec; |
|
if (tx != t_COMPLEX || gsigne((GEN)x[2]) <= 0) |
if (tx != t_COMPLEX || gcmp((GEN)x[2],gzero)<=0) |
|
err(talker,"argument must belong to upper half-plane"); |
err(talker,"argument must belong to upper half-plane"); |
|
|
if (!l) l=prec; p1=mppi(l); setexpo(p1,2); |
return gexp(gmul(x, PiI2(l)), l); /* e(x) */ |
p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */ |
|
q=gmul(x,p2); return gexp(q,prec); |
|
} |
} |
if (tx != t_POL && tx != t_SER) |
if (tx == t_SER) return x; |
err(talker,"bad argument for modular function"); |
if (tx != t_POL) err(talker,"bad argument for modular function"); |
|
|
return tayl(x,gvar(x),precdl); |
return tayl(x,gvar(x),precdl); |
} |
} |
|
|
y=gun; qn=gun; ps=gun; |
y=gun; qn=gun; ps=gun; |
if (tx==t_PADIC) |
if (tx==t_PADIC) |
{ |
{ |
if (valp(q) <= 0) |
if (valp(q) <= 0) err(talker,"non-positive valuation in eta"); |
err(talker,"non-positive valuation in inteta"); |
|
for(;;) |
for(;;) |
{ |
{ |
p1=gneg_i(gmul(ps,gmul(q,gsqr(qn)))); |
p1=gneg_i(gmul(ps,gmul(q,gsqr(qn)))); |
|
|
} |
} |
else |
else |
{ |
{ |
long l = 0, v = 0; /* gcc -Wall */ |
long l, v = 0; /* gcc -Wall */ |
long av = avma, lim = stack_lim(av,3); |
gpmem_t av = avma, lim = stack_lim(av, 3); |
|
|
if (is_scalar_t(tx)) l = -bit_accuracy(precision(q)); |
if (is_scalar_t(tx)) l = -bit_accuracy(precision(q)); |
else |
else |
{ |
{ |
v=gvar(q); tx = 0; |
v = gvar(q); l = lg(q)-2; tx = 0; |
if (valp(q) <= 0) |
if (valp(q) <= 0) err(talker,"non-positive valuation in eta"); |
err(talker,"non-positive valuation in inteta"); |
|
} |
} |
for(;;) |
for(;;) |
{ |
{ |
p1=gneg_i(gmul(ps,gmul(q,gsqr(qn)))); |
p1 = gneg_i(gmul(ps,gmul(q,gsqr(qn)))); |
y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn); |
/* qn = q^n |
y=gadd(y,ps); |
* ps = (-1)^n q^(n(3n+1)/2) |
|
* p1 = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */ |
|
y = gadd(y,p1); qn = gmul(qn,q); ps = gmul(p1,qn); |
|
y = gadd(y,ps); |
if (tx) |
if (tx) |
{ if (gexpo(ps)-gexpo(y) < l) return y; } |
{ if (gexpo(ps)-gexpo(y) < l) return y; } |
else |
else |
{ if (gval(ps,v)>=precdl) return y; } |
{ if (gval(ps,v) >= l) return y; } |
if (low_stack(lim, stack_lim(av,3))) |
if (low_stack(lim, stack_lim(av,3))) |
{ |
{ |
GEN *gptr[3]; |
if(DEBUGMEM>1) err(warnmem,"eta"); |
if(DEBUGMEM>1) err(warnmem,"inteta"); |
gerepileall(av, 3, &y, &qn, &ps); |
gptr[0]=&y; gptr[1]=&qn; gptr[2]=&ps; |
|
gerepilemany(av,gptr,3); |
|
} |
} |
} |
} |
} |
} |
|
|
GEN |
GEN |
eta(GEN x, long prec) |
eta(GEN x, long prec) |
{ |
{ |
long av = avma; |
gpmem_t av = avma; |
GEN q = qq(x,prec); |
GEN q = qq(x,prec); |
return gerepileupto(av,inteta(q)); |
return gerepileupto(av,inteta(q)); |
} |
} |
Line 1400 eta(GEN x, long prec) |
|
Line 1982 eta(GEN x, long prec) |
|
GEN |
GEN |
trueeta(GEN x, long prec) |
trueeta(GEN x, long prec) |
{ |
{ |
long tx=typ(x), av=avma, tetpil,l; |
long tx=typ(x), l; |
|
gpmem_t av=avma, tetpil; |
GEN p1,p2,q,q24,n,z,m,unapprox; |
GEN p1,p2,q,q24,n,z,m,unapprox; |
|
|
if (!is_scalar_t(tx)) err(typeer,"trueeta"); |
if (!is_scalar_t(tx)) err(typeer,"trueeta"); |
if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0) |
if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0) |
err(talker,"argument must belong to upper half-plane"); |
err(talker,"argument must belong to upper half-plane"); |
l=precision(x); if (l) prec=l; |
l=precision(x); if (l) prec=l; |
p1=mppi(prec); setexpo(p1,2); |
p2 = PiI2(prec); |
p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */ |
|
z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */ |
z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */ |
unapprox=gsub(gun,gpuigs(stoi(10),-8)); |
unapprox=gsub(gun,gpuigs(stoi(10),-8)); |
m=gun; |
m=gun; |
Line 1434 eta0(GEN x, long flag,long prec) |
|
Line 2016 eta0(GEN x, long flag,long prec) |
|
GEN |
GEN |
jell(GEN x, long prec) |
jell(GEN x, long prec) |
{ |
{ |
long av=avma,tetpil,tx = typ(x); |
long tx = typ(x); |
|
gpmem_t av=avma, tetpil; |
GEN p1; |
GEN p1; |
|
|
if (!is_scalar_t(tx) || tx == t_PADIC) |
if (!is_scalar_t(tx) || tx == t_PADIC) |
Line 1456 jell(GEN x, long prec) |
|
Line 2039 jell(GEN x, long prec) |
|
GEN |
GEN |
wf2(GEN x, long prec) |
wf2(GEN x, long prec) |
{ |
{ |
long av=avma,tetpil; |
gpmem_t av=avma, tetpil; |
GEN p1,p2; |
GEN p1,p2; |
|
|
p1=gsqrt(gdeux,prec); |
p1=gsqrt(gdeux,prec); |
Line 1468 wf2(GEN x, long prec) |
|
Line 2051 wf2(GEN x, long prec) |
|
GEN |
GEN |
wf1(GEN x, long prec) |
wf1(GEN x, long prec) |
{ |
{ |
long av=avma,tetpil; |
gpmem_t av=avma, tetpil; |
GEN p1,p2; |
GEN p1,p2; |
|
|
p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec); |
p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec); |
Line 1479 wf1(GEN x, long prec) |
|
Line 2062 wf1(GEN x, long prec) |
|
GEN |
GEN |
wf(GEN x, long prec) |
wf(GEN x, long prec) |
{ |
{ |
long av=avma,tetpil; |
gpmem_t av=avma, tetpil; |
GEN p1,p2; |
GEN p1,p2; |
|
|
p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec)); |
p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec)); |
Line 1506 sagm(GEN x, long prec) |
|
Line 2089 sagm(GEN x, long prec) |
|
{ |
{ |
GEN p1,a,b,a1,b1,y; |
GEN p1,a,b,a1,b1,y; |
long l,ep; |
long l,ep; |
ulong av; |
gpmem_t av; |
|
|
if (gcmp0(x)) return gcopy(x); |
if (gcmp0(x)) return gcopy(x); |
switch(typ(x)) |
switch(typ(x)) |
Line 1570 sagm(GEN x, long prec) |
|
Line 2153 sagm(GEN x, long prec) |
|
GEN |
GEN |
agm(GEN x, GEN y, long prec) |
agm(GEN x, GEN y, long prec) |
{ |
{ |
long av,tetpil, ty=typ(y); |
long ty=typ(y); |
|
gpmem_t av, tetpil; |
GEN z; |
GEN z; |
|
|
if (is_matvec_t(ty)) |
if (is_matvec_t(ty)) |
Line 1587 agm(GEN x, GEN y, long prec) |
|
Line 2171 agm(GEN x, GEN y, long prec) |
|
GEN |
GEN |
logagm(GEN q) |
logagm(GEN q) |
{ |
{ |
long av=avma,prec=lg(q),tetpil,s,n,lim; |
long prec=lg(q), s, n, lim; |
|
gpmem_t av=avma, tetpil; |
GEN y,q4,q1; |
GEN y,q4,q1; |
|
|
if (typ(q)!=t_REAL) err(typeer,"logagm"); |
if (typ(q)!=t_REAL) err(typeer,"logagm"); |
|
|
GEN |
GEN |
glogagm(GEN x, long prec) |
glogagm(GEN x, long prec) |
{ |
{ |
long av,tetpil; |
gpmem_t av, tetpil; |
GEN y,p1,p2; |
GEN y,p1,p2; |
|
|
switch(typ(x)) |
switch(typ(x)) |
Line 1644 glogagm(GEN x, long prec) |
|
Line 2229 glogagm(GEN x, long prec) |
|
GEN |
GEN |
theta(GEN q, GEN z, long prec) |
theta(GEN q, GEN z, long prec) |
{ |
{ |
long av=avma,tetpil,l,n; |
long l, n; |
|
gpmem_t av=avma, tetpil; |
GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold; |
GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold; |
|
|
|
if (!is_scalar_t(typ(q)) || !is_scalar_t(typ(z))) |
|
err(impl,"theta for non-scalar types"); |
|
|
l=precision(q); if (l) prec=l; |
l=precision(q); if (l) prec=l; |
p1=realun(prec); z=gmul(p1,z); |
p1=realun(prec); z=gmul(p1,z); |
if (!l) q=gmul(p1,q); |
if (!l) q=gmul(p1,q); |
Line 1681 theta(GEN q, GEN z, long prec) |
|
Line 2270 theta(GEN q, GEN z, long prec) |
|
GEN |
GEN |
thetanullk(GEN q, long k, long prec) |
thetanullk(GEN q, long k, long prec) |
{ |
{ |
long av=avma,tetpil,l,n; |
long l, n; |
|
gpmem_t av=avma, tetpil; |
GEN p1,ps,qn,y,ps2; |
GEN p1,ps,qn,y,ps2; |
|
|
l=precision(q); |
l=precision(q); |
Line 1703 thetanullk(GEN q, long k, long prec) |
|
Line 2293 thetanullk(GEN q, long k, long prec) |
|
return gerepile(av,tetpil,gmul(p1,y)); |
return gerepile(av,tetpil,gmul(p1,y)); |
} |
} |
|
|
GEN |
|
jbesselh(GEN n, GEN z, long prec) |
|
{ |
|
long av,tetpil,k,l,i,lz; |
|
GEN y,p1,p2,zinv,p0,s,c; |
|
|
|
if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh"); |
|
k=itos(n); if (k<0) err(impl,"ybesselh"); |
|
|
|
switch(typ(z)) |
|
{ |
|
case t_REAL: case t_COMPLEX: |
|
if (gcmp0(z)) return gzero; |
|
av=avma; zinv=ginv(z); |
|
l=precision(z); if (l>prec) prec=l; |
|
gsincos(z,&s,&c,prec); |
|
p1=gmul(zinv,s); |
|
if (k) |
|
{ |
|
p0=p1; p1=gmul(zinv,gsub(p0,c)); |
|
for (i=2; i<=k; i++) |
|
{ |
|
p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0); |
|
p0=p1; p1=p2; |
|
} |
|
} |
|
p2=gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec); |
|
tetpil=avma; return gerepile(av,tetpil,gmul(p2,p1)); |
|
|
|
case t_VEC: case t_COL: case t_MAT: |
|
lz=lg(z); y=cgetg(lz,typ(z)); |
|
for (i=1; i<lz; i++) |
|
y[i]=(long)jbesselh(n,(GEN)z[i],prec); |
|
return y; |
|
|
|
case t_INT: case t_FRAC: case t_FRACN: |
|
av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselh(n,p1,prec)); |
|
|
|
case t_QUAD: |
|
av=avma; p1=gmul(z,realun(prec)); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselh(n,p1,prec)); |
|
|
|
case t_POL: case t_RFRAC: case t_RFRACN: |
|
av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma; |
|
return gerepile(av,tetpil,jbesselh(n,p1,prec)); |
|
|
|
case t_POLMOD: |
|
av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]); |
|
tetpil=avma; y=cgetg(lz,t_COL); |
|
for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)p2[i],prec); |
|
return gerepile(av,tetpil,y); |
|
|
|
case t_PADIC: |
|
err(impl,"p-adic jbessel function"); |
|
case t_SER: |
|
err(impl,"jbessel of power series"); |
|
} |
|
err(typeer,"jbesselh"); |
|
return NULL; /* not reached */ |
|
} |
|