OpenXM_contrib/pari-2.2/src/basemath/trans3.c - diff

Return to trans3.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath

Diff for /OpenXM_contrib/pari-2.2/src/basemath/Attic/trans3.c between version 1.1.1.1 and 1.2

version 1.1.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,

/** **/

/********************************************************************/

#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

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;

}

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

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 */

}

GEN

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);

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);

e=cgetr(l1); f=cgetr(l1);

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);

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;

for(;;)

{

Line 131 kbessel(GEN nu, GEN gx, long prec)

Line 453 kbessel(GEN nu, GEN gx, long prec)

}

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 **/

Line 149 GEN

Line 719 GEN

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;

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");

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));

}

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;

{

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

kbessel2(GEN nu, GEN x, long prec)

{

long av = avma,tetpil;

gpmem_t av = avma, tetpil;

GEN p1,p2,x2,a,pitemp;

if (typ(x)==t_REAL) prec = lg(x);

Line 242 GEN

Line 813 GEN

incgam(GEN s, GEN x, long prec)

{

GEN p1,z = cgetr(prec);

long av = avma;

gpmem_t av = avma;

if (gcmp0(x)) return ggamma(s,prec);

if (typ(x)!=t_REAL) { gaffect(x,z); x=z; }

Line 273 GEN

Line 844 GEN

incgam1(GEN a, GEN x, long 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;

if (typ(x) != t_REAL) { gaffect(x,z); x=z; }

Line 296 GEN

Line 868 GEN

incgam2(GEN a, GEN x, long 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;

if (typ(x) != t_REAL) { gaffect(x,z); x=z; }

Line 328 GEN

Line 901 GEN

incgam3(GEN s, GEN x, long 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; }

l=lg(x); n = -bit_accuracy(l)-1;

Line 364 GEN

Line 938 GEN

incgam4(GEN a, GEN x, GEN g, long prec)

{

GEN p1, z = cgetr(prec);

long av = avma;

gpmem_t av = avma;

if (typ(x) != t_REAL) { gaffect(x,z); x=z; }

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

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;

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

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;

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

gerfc(GEN x, long prec)

{

long av=avma;

gpmem_t av=avma;

GEN p1,p2;

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

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;

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

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;

/* treat trivial cases */

Line 689 izeta(long k, long prec)

Line 1267 izeta(long k, long prec)

return gerepileuptoleaf(av, divrs(y,k-1));

}

if (k > bit_accuracy(prec)+1) return realun(prec);

pi2 = mppi(prec); setexpo(pi2,2); /* 2Pi */

pi2 = Pi2n(1, prec);

if ((k&1) == 0)

{

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);

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++)

{

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;

qn = mulrr(qn,q);

if (low_stack(limit,stack_lim(av2,1)))

{

GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;

if (DEBUGMEM>1) err(warnmem,"izeta");

gerepilemany(av2,gptr,2);

gerepileall(av2,2, &z, &qn);

}

setexpo(z,expo(z)+1);

Line 757 izeta(long k, long prec)

Line 1334 izeta(long k, long prec)

qn=mulrr(qn,q);

if (low_stack(limit,stack_lim(av2,1)))

{

GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;

if (DEBUGMEM>1) err(warnmem,"izeta");

gerepilemany(av2,gptr,2);

gerepileall(av2,2, &z, &qn);

}

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 */

static GEN

n_s(long n, GEN *tab)

n_s(ulong n, GEN *tab)

{

byteptr d = diffptr + 2;

GEN x = NULL;

long p, e;

for (p = 3; n > 1; p += *d++)

for (p = 3; n > 1; )

{

e = svaluation(n, p, &n);

if (e)

Line 794 n_s(long n, GEN *tab)

Line 1370 n_s(long n, GEN *tab)

GEN y = tab[pows(p,e)];

if (!x) x = y; else x = gmul(x,y);

}

NEXT_PRIME_VIADIFF_CHECK(p,d);

}

return x;

}

Line 804 GEN czeta(GEN s0, long prec);

Line 1381 GEN czeta(GEN s0, long prec);

static GEN

izeta(long k, long prec)

{

long av = avma;

gpmem_t av = avma;

GEN y,p1,pi2;

/* treat trivial cases */

Line 827 izeta(long k, long prec)

Line 1404 izeta(long k, long 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.

* 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 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;

byteptr d;

if (DEBUGLEVEL) timer2();

if (DEBUGLEVEL>2) (void)timer2();

s = trans_fix_arg(&prec,&s0,&sig,&av,&res);

if (gcmp0(s)) { y = gneg(ghalf); goto END; }

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;

l2 = (lim+ssig/2.-.25);

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);

}

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;

if (typ(s0) == t_INT)

{ /* 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));

NEXT_PRIME_VIADIFF_CHECK(p,d);

}

a = divrr(unr, rpowsi(nn, s0, prec));

}

else

{ /* general case */

ms = gneg(s); p1 = cgetr(prec);

for (p=2; p < nn; p += *d++)

for (p=2; p < nn;)

{

affsr(p, p1);

tab[p] = gexp(gmul(ms, mplog(p1)), prec);

NEXT_PRIME_VIADIFF_CHECK(p,d);

}

affsr(nn,p1);

a = gexp(gmul(ms, mplog(p1)), prec);

}

sqn = (long)sqrt(nn-1.);

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;

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;

for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1;

sim = y = unr;

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++)

sim = gadd(sim, n_s(2*j+1, tab));

sim = gerepileupto(av2, sim);

y = gadd(sim, gmul(tab[2],y));

}

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;

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));

tes = gmulsg(nn, gaddsg(1, u));

}

if (DEBUGLEVEL) msgtimer("Bernoulli sum");

if (DEBUGLEVEL>2) msgtimer("Bernoulli sum");

/* y += tes a / (s-1) */

y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr))));

Line 984 END:

Line 1561 END:

if (funeq)

{

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 */

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

gzetaz(GEN x, GEN y)

{

long av=avma, prec = precision(y);

long prec = precision(y);

gpmem_t av=avma;

if (!prec) err(infprecer,"gzetaz");

gaffect(gzeta(x,prec),y); avma=av;

Line 1038 gzetaz(GEN x, GEN y)

Line 1617 gzetaz(GEN x, GEN y)

static GEN

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;

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

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;

if (m<0) err(talker,"negative index in polylog");

Line 1142 dilog(GEN x, long prec)

Line 1723 dilog(GEN x, long prec)

GEN

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;

m2=m&1; av=avma;

Line 1188 polylogdold(long m, GEN x, long prec)

Line 1770 polylogdold(long m, GEN x, long prec)

GEN

polylogp(long m, GEN x, long prec)

{

long k,l,av,fl,m2;

long k, l, fl, m2;

gpmem_t av;

GEN p1,y;

m2=m&1; av=avma;

Line 1234 polylogp(long m, GEN x, long prec)

Line 1817 polylogp(long m, GEN x, long prec)

GEN

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;

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);

return gerepileupto(av, gneg(p1));

}

if (gcmp0(x)) return gcopy(x);

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);

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

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");

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 (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");

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);

}

Line 1347 inteta(GEN q)

Line 1930 inteta(GEN q)

y=gun; qn=gun; ps=gun;

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(;;)

{

p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));

Line 1358 inteta(GEN q)

Line 1940 inteta(GEN q)

}

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));

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(;;)

{

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 (gexpo(ps)-gexpo(y) < l) return y; }

else

{ if (gval(ps,v)>=precdl) return y; }

{ if (gval(ps,v) >= l) return y; }

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);

}

Line 1391 inteta(GEN q)

Line 1973 inteta(GEN q)

GEN

eta(GEN x, long prec)

{

long av = avma;

gpmem_t av = avma;

GEN q = qq(x,prec);

return gerepileupto(av,inteta(q));

}

Line 1400 eta(GEN x, long prec)

Line 1982 eta(GEN x, long prec)

GEN

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;

if (!is_scalar_t(tx)) err(typeer,"trueeta");

if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0)

err(talker,"argument must belong to upper half-plane");

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) */

unapprox=gsub(gun,gpuigs(stoi(10),-8));

m=gun;

Line 1434 eta0(GEN x, long flag,long prec)

Line 2016 eta0(GEN x, long flag,long prec)

GEN

jell(GEN x, long prec)

{

long av=avma,tetpil,tx = typ(x);

long tx = typ(x);

gpmem_t av=avma, tetpil;

GEN p1;

if (!is_scalar_t(tx) || tx == t_PADIC)

Line 1456 jell(GEN x, long prec)

Line 2039 jell(GEN x, long prec)

GEN

wf2(GEN x, long prec)

{

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

GEN p1,p2;

p1=gsqrt(gdeux,prec);

Line 1468 wf2(GEN x, long prec)

Line 2051 wf2(GEN x, long prec)

GEN

wf1(GEN x, long prec)

{

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

GEN p1,p2;

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

wf(GEN x, long prec)

{

long av=avma,tetpil;

gpmem_t av=avma, tetpil;

GEN p1,p2;

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;

long l,ep;

ulong av;

gpmem_t av;

if (gcmp0(x)) return gcopy(x);

switch(typ(x))

Line 1570 sagm(GEN x, long prec)

Line 2153 sagm(GEN x, long prec)

GEN

agm(GEN x, GEN y, long prec)

{

long av,tetpil, ty=typ(y);

long ty=typ(y);

gpmem_t av, tetpil;

GEN z;

if (is_matvec_t(ty))

Line 1587 agm(GEN x, GEN y, long prec)

Line 2171 agm(GEN x, GEN y, long prec)

GEN

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;

if (typ(q)!=t_REAL) err(typeer,"logagm");

Line 1606 logagm(GEN q)

Line 2191 logagm(GEN q)

GEN

glogagm(GEN x, long prec)

{

long av,tetpil;

gpmem_t av, tetpil;

GEN y,p1,p2;

switch(typ(x))

Line 1644 glogagm(GEN x, long prec)

Line 2229 glogagm(GEN x, long prec)

GEN

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;

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;

p1=realun(prec); z=gmul(p1,z);

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

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;

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));

}

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 */

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>