Annotation of OpenXM_contrib/pari-2.2/src/basemath/trans3.c, Revision 1.2
1.2 ! noro 1: /* $Id: trans3.c,v 1.49 2002/07/15 13:30:01 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /********************************************************************/
17: /** **/
18: /** TRANSCENDENTAL FONCTIONS **/
19: /** (part 3) **/
20: /** **/
21: /********************************************************************/
22: #include "pari.h"
1.2 ! noro 23: extern GEN rpowsi(ulong a, GEN n, long prec);
! 24: extern GEN divrs2_safe(GEN x, long i);
! 25: extern void dcxlog(double s, double t, double *a, double *b);
! 26: extern double dnorm(double s, double t);
! 27: extern GEN trans_fix_arg(long *prec, GEN *s0, GEN *sig, gpmem_t *av, GEN *res);
1.1 noro 28:
29: GEN
30: cgetc(long l)
31: {
32: GEN u = cgetg(3,t_COMPLEX); u[1]=lgetr(l); u[2]=lgetr(l);
33: return u;
34: }
35:
36: static GEN
37: fix(GEN x, long l)
38: {
39: GEN y;
40: if (typ(x) == t_REAL) return x;
41: y = cgetr(l); gaffect(x, y); return y;
42: }
43:
1.2 ! noro 44: long
! 45: lgcx(GEN z)
! 46: {
! 47: long tz = typ(z);
! 48:
! 49: switch(tz)
! 50: {
! 51: case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: return BIGINT;
! 52: case t_REAL: return lg(z);
! 53: case t_COMPLEX: return min(lgcx((GEN)z[1]),lgcx((GEN)z[2]));
! 54: default: err(typeer,"lgcx");
! 55: }
! 56: return 0;
! 57: }
! 58:
! 59: GEN
! 60: setlgcx(GEN z, long l)
! 61: {
! 62: long tz = typ(z);
! 63: GEN p1;
! 64:
! 65: switch(tz)
! 66: {
! 67: case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: return z;
! 68: case t_REAL: p1 = cgetr(l); affrr(z,p1); return p1;
! 69: case t_COMPLEX: p1 = cgetc(l); gaffect(z,p1); return p1;
! 70: default: err(typeer,"setlgcx"); return gzero; /* not reached */
! 71: }
! 72: }
! 73:
! 74: /* force z to be of type real/complex */
! 75:
! 76: GEN
! 77: setlgcx2(GEN z, long l)
! 78: {
! 79: long tz = typ(z);
! 80: GEN p1;
! 81:
! 82: switch(tz)
! 83: {
! 84: case t_INT: case t_FRAC: case t_RFRAC: case t_QUAD: case t_REAL:
! 85: p1 = cgetr(l); gaffect(z,p1); return p1;
! 86: case t_COMPLEX: p1 = cgetc(l); gaffect(z,p1); return p1;
! 87: default: err(typeer,"setlgcx"); return gzero; /* not reached */
! 88: }
! 89: }
! 90:
! 91: /* a exporter ou ca existe deja ? */
! 92:
! 93: long
! 94: isint(GEN n, long *ptk)
! 95: {
! 96: long tn=typ(n);
! 97: GEN p1,p2;
! 98:
! 99: switch(tn)
! 100: {
! 101: case t_INT:
! 102: *ptk = itos(n); return 1;
! 103: case t_REAL:
! 104: p1 = gfloor(n);
! 105: if (gcmp(p1,n)==0) {*ptk = itos(p1); return 1;}
! 106: else return 0;
! 107: case t_FRAC: case t_FRACN:
! 108: p1 = dvmdii((GEN)n[1],(GEN)n[2],&p2);
! 109: if (!signe(p2)) {*ptk = itos(p1); return 1;}
! 110: else return 0;
! 111: case t_COMPLEX:
! 112: if (gcmp0((GEN)n[2])) return isint((GEN)n[1],ptk);
! 113: else return 0;
! 114: case t_QUAD:
! 115: if (gcmp0((GEN)n[3])) return isint((GEN)n[2],ptk);
! 116: else return 0;
! 117: default: err(typeer,"isint"); return 0; /* not reached */
! 118: }
! 119: }
! 120:
! 121: double
! 122: norml1(GEN n, long prec)
! 123: {
! 124: long tn=typ(n);
! 125: gpmem_t av;
! 126: double res;
! 127:
! 128: switch(tn)
! 129: {
! 130: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 131: return gtodouble(gabs(n,prec));
! 132: case t_COMPLEX:
! 133: return norml1((GEN)n[1],prec)+norml1((GEN)n[2],prec);
! 134: case t_QUAD:
! 135: av = avma; res = norml1(gmul(n,realun(prec)),prec); avma = av;
! 136: return res;
! 137: default: err(typeer,"norml1"); return 0.; /* not reached */
! 138: }
! 139: }
! 140:
1.1 noro 141: /***********************************************************************/
142: /** **/
1.2 ! noro 143: /** BESSEL FUNCTIONS **/
1.1 noro 144: /** **/
145: /***********************************************************************/
146:
1.2 ! noro 147: /* computes sum_{k=0}^m n!*((-1)^flag*z^2/4)^k/(k!*(k+n)!) */
! 148:
! 149: static GEN
! 150: _jbessel(GEN n, GEN z, long flag, long m)
! 151: {
! 152: long k, limit;
! 153: gpmem_t av;
! 154: GEN p1,s;
! 155:
! 156: p1 = gmul2n(gsqr(z),-2); if (flag & 1) p1 = gneg(p1);
! 157: if (typ(z) == t_SER)
! 158: {
! 159: long v = valp(z);
! 160: k = lg(p1)-2 - v;
! 161: if (v < 0) err(negexper,"jbessel");
! 162: if (k <= 0) return gadd(gun, zeroser(varn(z), 2*v));
! 163: p1 = gprec(p1, k);
! 164: }
! 165: s = gun;
! 166: av = avma; limit = stack_lim(av,1);
! 167: for (k=m; k>=1; k--)
! 168: {
! 169: s = gadd(gun,gdiv(gdivgs(gmul(p1,s),k),gaddsg(k,n)));
! 170: if (low_stack(limit,stack_lim(av,1)))
! 171: {
! 172: if (DEBUGMEM>1) err(warnmem,"jbessel");
! 173: s = gerepilecopy(av, s);
! 174: }
! 175: }
! 176: return s;
! 177: }
! 178:
! 179: static GEN
! 180: jbesselintern(GEN n, GEN z, long flag, long prec)
! 181: {
! 182: long tz=typ(z), i, lz, lim, k=-1, ki, precnew;
! 183: gpmem_t av=avma, tetpil;
! 184: double B,N,L,x;
! 185: GEN p1,p2,y,znew,nnew;
! 186:
! 187: switch(tz)
! 188: {
! 189: case t_REAL: case t_COMPLEX:
! 190: i = precision(z); if (i) prec = i;
! 191: if (isint(setlgcx2(n,prec),&ki))
! 192: {
! 193: k = labs(ki);
! 194: p2 = gdiv(gpowgs(gmul2n(z,-1),k),mpfact(k));
! 195: if ((flag&1) && (ki<0) && (k&1)) p2 = gneg(p2);
! 196: }
! 197: else p2 = gdiv(gpow(gmul2n(z,-1),n,prec),ggamma(gaddgs(n,1),prec));
! 198: if (gcmp0(z)) return gerepilecopy(av, p2);
! 199: else
! 200: {
! 201: x = gtodouble(gabs(z,prec));
! 202: L = x*1.3591409;
! 203: B = bit_accuracy(prec)*LOG2/(2*L);
! 204: N = 1 + B;
! 205: /* 3 Newton iterations are enough except in pathological cases */
! 206: N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1);
! 207: lim = max((long)(L*N),2);
! 208: precnew = prec;
! 209: if (x >= 1.0) precnew += 1 + (long)(x/(LOG2*BITS_IN_LONG));
! 210: znew = setlgcx(z,precnew);
! 211: if (k >= 0) p1 = setlgcx(_jbessel(stoi(k),znew,flag,lim),prec);
! 212: else
! 213: {
! 214: i = precision(n);
! 215: nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n;
! 216: p1 = setlgcx(_jbessel(nnew,znew,flag,lim),prec);
! 217: }
! 218: tetpil = avma; return gerepile(av,tetpil,gmul(p2,p1));
! 219: }
! 220:
! 221: case t_SER:
! 222: if (isint(setlgcx2(n,prec),&ki))
! 223: {
! 224: k = labs(ki);
! 225: p1 = _jbessel(stoi(k),z,flag,lg(z)-2);
! 226: }
! 227: else p1 = _jbessel(n,z,flag,lg(z)-2);
! 228: return gerepilecopy(av,p1);
! 229:
! 230: case t_VEC: case t_COL: case t_MAT:
! 231: lz=lg(z); y=cgetg(lz,typ(z));
! 232: for (i=1; i<lz; i++)
! 233: y[i]=(long)jbesselintern(n,(GEN)z[i],flag,prec);
! 234: return y;
! 235:
! 236: case t_INT: case t_FRAC: case t_FRACN:
! 237: av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma;
! 238: return gerepile(av,tetpil,jbesselintern(n,p1,flag,prec));
! 239:
! 240: case t_QUAD:
! 241: av=avma; p1=gmul(z,realun(prec)); tetpil=avma;
! 242: return gerepile(av,tetpil,jbesselintern(n,p1,flag,prec));
! 243:
! 244: case t_POL: case t_RFRAC: case t_RFRACN:
! 245: av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma;
! 246: return gerepile(av,tetpil,jbesselintern(n,p1,flag,prec));
! 247:
! 248: case t_POLMOD:
! 249: av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL);
! 250: for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]);
! 251: tetpil=avma; y=cgetg(lz,t_COL);
! 252: for (i=1; i<lz; i++) y[i]=(long)jbesselintern(n,(GEN)p2[i],flag,prec);
! 253: return gerepile(av,tetpil,y);
! 254:
! 255: case t_PADIC:
! 256: err(impl,"p-adic jbessel function");
! 257: }
! 258: err(typeer,"jbessel");
! 259: return NULL; /* not reached */
! 260: }
! 261:
! 262: GEN
! 263: jbessel(GEN n, GEN z, long prec)
! 264: {
! 265: return jbesselintern(n,z,1,prec);
! 266: }
! 267:
! 268: GEN
! 269: ibessel(GEN n, GEN z, long prec)
! 270: {
! 271: return jbesselintern(n,z,0,prec);
! 272: }
! 273:
! 274: GEN
! 275: _jbesselh(long k, GEN z, long prec)
! 276: {
! 277: GEN zinv,s,c,p0,p1,p2;
! 278: long i;
! 279:
! 280: zinv=ginv(z);
! 281: gsincos(z,&s,&c,prec);
! 282: p1=gmul(zinv,s);
! 283: if (k)
! 284: {
! 285: p0=p1; p1=gmul(zinv,gsub(p0,c));
! 286: for (i=2; i<=k; i++)
! 287: {
! 288: p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0);
! 289: p0=p1; p1=p2;
! 290: }
! 291: }
! 292: return p1;
! 293: }
! 294:
1.1 noro 295: GEN
1.2 ! noro 296: jbesselh(GEN n, GEN z, long prec)
1.1 noro 297: {
1.2 ! noro 298: long gz, k, l, linit, i, lz, tz=typ(z);
! 299: gpmem_t av, tetpil;
! 300: GEN y,p1,p2;
! 301:
! 302: if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh");
! 303: k = itos(n);
! 304: if (k < 0)
! 305: {
! 306: av = avma; n = gadd(ghalf,n); tetpil = avma;
! 307: return gerepile(av,tetpil,jbessel(n,z,prec));
! 308: }
! 309:
! 310: switch(tz)
1.1 noro 311: {
1.2 ! noro 312: case t_REAL: case t_COMPLEX:
! 313: av = avma;
! 314: if (gcmp0(z))
! 315: {
! 316: p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
! 317: p1 = gdiv(gmul(mpfact(k),p1),mpfact(2*k+1));
! 318: tetpil = avma; return gerepile(av,tetpil,gmul2n(p1,2*k));
! 319: }
! 320: gz = gexpo(z);
! 321: linit = lgcx(z); if (linit==BIGINT) linit = prec;
! 322: if (gz>=0) l = linit;
! 323: else l = lgcx(z) - 1 + ((-2*k*gz)>>TWOPOTBITS_IN_LONG);
! 324: if (l>prec) prec = l;
! 325: prec += (-gz)>>TWOPOTBITS_IN_LONG;
! 326: z = setlgcx(z,prec);
! 327: p1 = _jbesselh(k,z,prec);
! 328: p1 = gmul(gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec),p1);
! 329: tetpil = avma; return gerepile(av,tetpil,setlgcx(p1,linit));
! 330:
! 331: case t_SER:
! 332: if (gcmp0(z)) return gpowgs(z,k);
! 333: av = avma;
! 334: if (valp(z) < 0) err(negexper,"jbesselh");
! 335: z = gprec(z, lg(z)-2 + (2*k+1)*valp(z));
! 336: p1 = gdiv(_jbesselh(k,z,prec),gpowgs(z,k));
! 337: for (i=1; i<=k; i++) p1 = gmulgs(p1,2*i+1);
! 338: return gerepilecopy(av,p1);
! 339:
! 340: case t_VEC: case t_COL: case t_MAT:
! 341: lz=lg(z); y=cgetg(lz,typ(z));
! 342: for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)z[i],prec);
! 343: return y;
! 344:
! 345: case t_INT: case t_FRAC: case t_FRACN:
! 346: av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma;
! 347: return gerepile(av,tetpil,jbesselh(n,p1,prec));
! 348:
! 349: case t_QUAD:
! 350: av=avma; p1=gmul(z,realun(prec)); tetpil=avma;
! 351: return gerepile(av,tetpil,jbesselh(n,p1,prec));
! 352:
! 353: case t_POL: case t_RFRAC: case t_RFRACN:
! 354: av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma;
! 355: return gerepile(av,tetpil,jbesselh(n,p1,prec));
! 356:
! 357: case t_POLMOD:
! 358: av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL);
! 359: for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]);
! 360: tetpil=avma; y=cgetg(lz,t_COL);
! 361: for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)p2[i],prec);
! 362: return gerepile(av,tetpil,y);
! 363:
! 364: case t_PADIC:
! 365: err(impl,"p-adic jbesselh function");
1.1 noro 366: }
1.2 ! noro 367: err(typeer,"jbesselh");
1.1 noro 368: return NULL; /* not reached */
369: }
370:
371: GEN
372: kbessel(GEN nu, GEN gx, long prec)
373: {
1.2 ! noro 374: GEN x,y,yfin,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w;
! 375: long l, lnew, lbin, k, k2, l1, n2, n, ex, rab=0;
! 376: gpmem_t av, av1;
1.1 noro 377:
378: if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec);
379: l = (typ(gx)==t_REAL)? lg(gx): prec;
1.2 ! noro 380: ex = gexpo(gx);
! 381: if (ex < 0)
! 382: {
! 383: rab = 1 + ((-ex)>>TWOPOTBITS_IN_LONG);
! 384: lnew = l + rab; prec += rab;
! 385: }
! 386: else lnew = l;
! 387: yfin=cgetr(l); l1=lnew+1;
! 388: av=avma; x = fix(gx, lnew); y=cgetr(lnew);
1.1 noro 389: u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1);
390: e=cgetr(l1); f=cgetr(l1);
391: nu2=gmulgs(gsqr(nu),-4);
392: n = (long) (bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gnorm(nu)))) / 2;
393: n2=(n<<1); pitemp=mppi(l1);
394: /* this 10 should really be a 5, but then kbessel(15.99) enters oo loop */
395: lbin = 10 - bit_accuracy(l); av1=avma;
396: if (gcmpgs(x,n)<0)
397: {
398: zf=gsqrt(gdivgs(pitemp,n2),prec);
399: zz=cgetr(l1); gaffect(ginv(stoi(n2<<2)), zz);
400: s=gun; t=gzero;
401: for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
402: {
403: if (k2 & ~(MAXHALFULONG>>1))
404: p1 = gadd(mulss(k2,k2),nu2);
405: else
406: p1 = gaddsg(k2*k2,nu2);
407: ak=gdivgs(gmul(p1,zz),-k);
408: s=gadd(gun,gmul(ak,s));
409: t=gaddsg(k2,gmul(ak,t));
410: }
411: gmulz(s,zf,u); t=gmul2n(t,-1);
412: gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1;
1.2 ! noro 413: q = stor(n2, l1);
1.1 noro 414: r=gmul2n(x,1); av1=avma;
415: for(;;)
416: {
417: p1=divsr(5,q); if (expo(p1) >= -1) p1=ghalf;
418: p2=subsr(1,divrr(r,q));
419: if (gcmp(p1,p2)>0) p1=p2;
420: gnegz(p1,c); gaffsg(1,d); affrr(u,e); affrr(v,f);
421: for (k=1; ; k++)
422: {
423: w=gadd(gmul(gsubsg(k,ghalf),u), gmul(subrs(q,k),v));
424: w=gadd(w, gmul(nu,gsub(u,gmul2n(v,1))));
425: gdivgsz(gmul(q,v),k,u);
426: gdivgsz(w,k,v);
427: gmulz(d,c,d);
428: gaddz(e,gmul(d,u),e); p1=gmul(d,v);
429: gaddz(f,p1,f);
430: if (gexpo(p1)-gexpo(f) <= lbin) break;
431: avma=av1;
432: }
433: p1=u; u=e; e=p1;
434: p1=v; v=f; f=p1;
435: gmulz(q,gaddsg(1,c),q);
436: if (expo(subrr(q,r)) <= lbin) break;
437: }
438: gmulz(u,gpui(gdivgs(x,n),nu,prec),y);
439: }
440: else
441: {
442: p2=gmul2n(x,1);
443: zf=gsqrt(gdiv(pitemp,p2),prec);
444: zz=ginv(gmul2n(p2,2)); s=gun;
445: for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
446: {
447: if (k2 & ~(MAXHALFULONG>>1))
448: p1=gadd(mulss(k2,k2),nu2);
449: else
450: p1=gaddsg(k2*k2,nu2);
451: ak=gdivgs(gmul(p1,zz),k);
452: s=gsub(gun,gmul(ak,s));
453: }
454: gmulz(s,zf,y);
455: }
1.2 ! noro 456: gdivz(y,mpexp(x),yfin);
! 457: avma=av; return yfin;
! 458: }
! 459:
! 460: /* computes sum_{k=0}^m ((-1)^flag*z^2/4)^k/(k!*(k+n)!)*(H(k)+H(k+n))
! 461: +sum_{k=0}^{n-1}((-1)^(flag+1)*z^2/4)^(k-n)*(n-k-1)!/k!. Ici n
! 462: doit etre un long. Attention, contrairement a _jbessel, pas de n! devant.
! 463: Quand flag > 1, calculer exactement les H(k) et factorielles */
! 464:
! 465: static GEN
! 466: _kbessel(long n, GEN z, long flag, long m, long prec)
! 467: {
! 468: long k, limit;
! 469: gpmem_t av;
! 470: GEN p1,p2,p3,s,*tabh;
! 471:
! 472: p1 = gmul2n(gsqr(z),-2); if (flag & 1) p1 = gneg(p1);
! 473: if (typ(z) == t_SER)
! 474: {
! 475: long v = valp(z);
! 476: k = lg(p1)-2 - v;
! 477: if (v < 0) err(negexper,"kbessel");
! 478: if (k <= 0) return gadd(gun, zeroser(varn(z), 2*v));
! 479: p1 = gprec(p1, k);
! 480: }
! 481: tabh = (GEN*)cgetg(m+n+2,t_VEC); tabh[1] = gzero;
! 482: if (flag <= 1)
! 483: {
! 484: s = realun(prec); tabh[2] = s;
! 485: for (k=2; k<=m+n; k++)
! 486: {
! 487: s = divrs(addsr(1,mulsr(k,s)),k); tabh[k+1] = s;
! 488: }
! 489: }
! 490: else
! 491: {
! 492: s = gun; tabh[2] = s;
! 493: for (k=2; k<=m+n; k++)
! 494: {
! 495: s = gdivgs(gaddsg(1,gmulsg(k,s)),k); tabh[k+1] = s;
! 496: }
! 497: }
! 498: s = gadd(tabh[m+1],tabh[m+n+1]);
! 499: av = avma; limit = stack_lim(av,1);
! 500: for (k=m; k>=1; k--)
! 501: {
! 502: s = gadd(gadd(tabh[k],tabh[k+n]),gdiv(gmul(p1,s),mulss(k,k+n)));
! 503: if (low_stack(limit,stack_lim(av,1)))
! 504: {
! 505: if (DEBUGMEM>1) err(warnmem,"kbessel");
! 506: s = gerepilecopy(av, s);
! 507: }
! 508: }
! 509: p3 = (flag <= 1) ? mpfactr(n,prec) : mpfact(n);
! 510: s = gdiv(s,p3);
! 511: if (n)
! 512: {
! 513: p1 = gneg(ginv(p1));
! 514: p2 = gmulsg(n,gdiv(p1,p3));
! 515: for (k=n-1; k>=0; k--)
! 516: {
! 517: s = gadd(s,p2);
! 518: p2 = gmul(p2,gmul(mulss(k,n-k),p1));
! 519: }
! 520: }
! 521: return s;
! 522: }
! 523:
! 524: static GEN
! 525: kbesselintern(GEN n, GEN z, long flag, long prec)
! 526: {
! 527: long tz=typ(z), i, k, ki, lz, lim, precnew, fl, fl2, ex;
! 528: gpmem_t av=avma, tetpil;
! 529: double B,N,L,x,rab;
! 530: GEN p1,p2,y,p3,znew,nnew,pplus,pmoins,s,c;
! 531:
! 532: fl = (flag & 1) == 0;
! 533: switch(tz)
! 534: {
! 535: case t_REAL: case t_COMPLEX:
! 536: if (gcmp0(z)) err(talker,"zero argument in a k/n bessel function");
! 537: i = precision(z); if (i) prec = i;
! 538: x = gtodouble(gabs(z,prec));
! 539: /* Experimentally optimal on a PIII 500 Mhz. Your optimum may vary. */
! 540: if ((x > (8*(prec-2)+norml1(n,prec)-3)) && !flag) return kbessel(n,z,prec);
! 541: precnew = prec;
! 542: if (x >= 1.0)
! 543: {
! 544: rab = x/(LOG2*BITS_IN_LONG); if (fl) rab *= 2;
! 545: precnew += 1 + (long)rab;
! 546: }
! 547: znew = setlgcx(z,precnew);
! 548: if (isint(setlgcx2(n,precnew),&ki))
! 549: {
! 550: k = labs(ki);
! 551: L = x*1.3591409;
! 552: B = (bit_accuracy(prec)*LOG2)/(2*L);
! 553: if (fl) B += 0.367879;
! 554: N = 1 + B;
! 555: /* 3 Newton iterations are enough except in pathological cases */
! 556: N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1); N = (N + B)/(log(N)+1);
! 557: lim = max((long)(L*N),2);
! 558: p1 = _kbessel(k,znew,flag,lim,precnew);
! 559: p1 = gmul(gpowgs(gmul2n(znew,-1),k),p1);
! 560: p2 = gadd(mpeuler(precnew),glog(gmul2n(znew,-1),precnew));
! 561: p3 = jbesselintern(stoi(k),znew,flag,precnew);
! 562: p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
! 563: p2 = setlgcx(p2,prec);
! 564: if (fl == 0) {p1 = mppi(prec); p2 = gmul2n(gdiv(p2,p1),1);}
! 565: fl = (fl && (k&1)) || (!fl && (ki>=0 || (k&1)==0));
! 566: tetpil = avma; return gerepile(av,tetpil,fl ? gneg(p2) : gcopy(p2));
! 567: }
! 568: else
! 569: {
! 570: i = precision(n);
! 571: nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n;
! 572: p2 = mppi(precnew); gsincos(gmul(nnew,p2),&s,&c,precnew);
! 573: ex = gexpo(s);
! 574: if (ex < 0)
! 575: {
! 576: rab = (-ex)/(LOG2*BITS_IN_LONG); if (fl) rab *= 2;
! 577: precnew += 1 + (long)rab;
! 578: }
! 579: nnew = (i && (i < precnew)) ? setlgcx(n,precnew) : n;
! 580: znew = setlgcx(znew,precnew);
! 581: p2 = mppi(precnew); gsincos(gmul(nnew,p2),&s,&c,precnew);
! 582: pplus = jbesselintern(nnew,znew,flag,precnew);
! 583: pmoins = jbesselintern(gneg(nnew),znew,flag,precnew);
! 584: if (fl) p1 = gmul(gsub(pmoins,pplus),gdiv(p2,gmul2n(s,1)));
! 585: else p1 = gdiv(gsub(gmul(c,pplus),pmoins),s);
! 586: tetpil = avma; return gerepile(av,tetpil,setlgcx(p1,prec));
! 587: }
! 588:
! 589: case t_SER:
! 590: if (isint(setlgcx2(n,prec),&ki))
! 591: {
! 592: k = labs(ki);
! 593: p1 = _kbessel(k,z,flag+2,lg(z)-2,prec);
! 594: return gerepilecopy(av,p1);
! 595: }
! 596: else
! 597: {
! 598: fl2 = isint(setlgcx2(gmul2n(n,1),prec),&ki);
! 599: if (!fl2)
! 600: err(talker,"cannot give a power series result in k/n bessel function");
! 601: else
! 602: {
! 603: k = labs(ki); n = gmul2n(stoi(k),-1);
! 604: fl2 = (k&3)==1;
! 605: pmoins = jbesselintern(gneg(n),z,flag,prec);
! 606: if (fl)
! 607: {
! 608: pplus = jbesselintern(n,z,flag,prec);
! 609: p2 = gpowgs(z,-k); if (fl2 == 0) p2 = gneg(p2);
! 610: p3 = gmul2n(divii(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
! 611: p3 = gdivgs(gmul2n(gsqr(p3),1),k);
! 612: p2 = gmul(p2,p3);
! 613: p1 = gsub(pplus,gmul(p2,pmoins));
! 614: }
! 615: else p1 = pmoins;
! 616: tetpil = avma;
! 617: return gerepile(av,tetpil,fl2 ? gneg(p1) : gcopy(p1));
! 618: }
! 619: }
! 620:
! 621: case t_VEC: case t_COL: case t_MAT:
! 622: lz=lg(z); y=cgetg(lz,typ(z));
! 623: for (i=1; i<lz; i++)
! 624: y[i]=(long)kbesselintern(n,(GEN)z[i],flag,prec);
! 625: return y;
! 626:
! 627: case t_INT: case t_FRAC: case t_FRACN:
! 628: av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma;
! 629: return gerepile(av,tetpil,kbesselintern(n,p1,flag,prec));
! 630:
! 631: case t_QUAD:
! 632: av=avma; p1=gmul(z,realun(prec)); tetpil=avma;
! 633: return gerepile(av,tetpil,kbesselintern(n,p1,flag,prec));
! 634:
! 635: case t_POL: case t_RFRAC: case t_RFRACN:
! 636: av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma;
! 637: return gerepile(av,tetpil,kbesselintern(n,p1,flag,prec));
! 638:
! 639: case t_POLMOD:
! 640: av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL);
! 641: for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]);
! 642: tetpil=avma; y=cgetg(lz,t_COL);
! 643: for (i=1; i<lz; i++) y[i]=(long)kbesselintern(n,(GEN)p2[i],flag,prec);
! 644: return gerepile(av,tetpil,y);
! 645:
! 646: case t_PADIC:
! 647: err(impl,"p-adic kbessel function");
! 648: }
! 649: err(typeer,"kbessel");
! 650: return NULL; /* not reached */
! 651: }
! 652:
! 653: GEN
! 654: kbesselnew(GEN n, GEN z, long prec)
! 655: {
! 656: return kbesselintern(n,z,0,prec);
! 657: }
! 658:
! 659: GEN
! 660: kbesselnewalways(GEN n, GEN z, long prec)
! 661: {
! 662: return kbesselintern(n,z,2,prec);
! 663: }
! 664:
! 665: GEN
! 666: nbessel(GEN n, GEN z, long prec)
! 667: {
! 668: return kbesselintern(n,z,1,prec);
! 669: }
! 670:
! 671: GEN
! 672: hbessel1(GEN n, GEN z, long prec)
! 673: {
! 674: gpmem_t av=avma, tetpil;
! 675: GEN p1,p2;
! 676:
! 677: p1 = jbessel(n,z,prec); p2 = gmul(gi,nbessel(n,z,prec));
! 678: tetpil = avma; return gerepile(av,tetpil,gadd(p1,p2));
! 679: }
! 680:
! 681: GEN
! 682: hbessel2(GEN n, GEN z, long prec)
! 683: {
! 684: gpmem_t av=avma, tetpil;
! 685: GEN p1,p2;
! 686:
! 687: p1 = jbessel(n,z,prec); p2 = gmul(gi,nbessel(n,z,prec));
! 688: tetpil = avma; return gerepile(av,tetpil,gsub(p1,p2));
! 689: }
! 690:
! 691: GEN kbessel2(GEN nu, GEN x, long prec);
! 692:
! 693: GEN
! 694: kbessel0(GEN nu, GEN gx, long flag, long prec)
! 695: {
! 696: switch(flag)
! 697: {
! 698: case 0: return kbesselnew(nu,gx,prec);
! 699: case 1: return kbessel(nu,gx,prec);
! 700: case 2: return kbessel2(nu,gx,prec);
! 701: case 3: return kbesselnewalways(nu,gx,prec);
! 702: default: err(flagerr,"besselk");
! 703: }
! 704: return NULL; /* not reached */
1.1 noro 705: }
706:
707: /***********************************************************************/
708: /* **/
709: /** FONCTION U(a,b,z) GENERALE **/
710: /** ET CAS PARTICULIERS **/
711: /** **/
712: /***********************************************************************/
713: /* Assume gx > 0 and a,b complex */
714: /* This might one day be extended to handle complex gx */
715: /* see Temme, N. M. "The numerical computation of the confluent */
716: /* hypergeometric function U(a,b,z)" in Numer. Math. 41 (1983), */
717: /* no. 1, 63--82. */
718: GEN
719: hyperu(GEN a, GEN b, GEN gx, long prec)
720: {
721: GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn;
1.2 ! noro 722: long l, lbin, k, l1, n, ex;
! 723: gpmem_t av, av1, av2;
1.1 noro 724:
725: if(gsigne(gx) <= 0) err(talker,"hyperu's third argument must be positive");
726: ex = (iscomplex(a) || iscomplex(b));
727:
728: l = (typ(gx)==t_REAL)? lg(gx): prec;
729: if (ex) y=cgetc(l); else y=cgetr(l);
730: l1=l+1; av=avma; x = fix(gx, l);
731: a1=gaddsg(1,gsub(a,b));
732: if (ex)
733: {
734: u=cgetc(l1); v=cgetc(l1); c=cgetc(l1);
735: d=cgetc(l1); e=cgetc(l1); f=cgetc(l1);
736: }
737: else
738: {
739: u=cgetr(l1); v=cgetr(l1); c=cgetr(l1);
740: d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
741: }
742: n=(long)(bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
743: lbin = 10-bit_accuracy(l); av1=avma;
744: if (gcmpgs(x,n)<0)
745: {
746: gn=stoi(n); zf=gpui(gn,gneg_i(a),l1);
747: zz=gdivsg(-1,gn); s=gun; t=gzero;
748: for (k=n-1; k>=0; k--)
749: {
750: p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
751: s=gaddsg(1,gmul(p1,s));
752: t=gadd(gaddsg(k,a),gmul(p1,t));
753: }
754: gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1;
1.2 ! noro 755: q = stor(n, l1); r=x; av1=avma;
1.1 noro 756: do
757: {
758: p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf;
759: p2=subsr(1,divrr(r,q)); if (gcmp(p1,p2)>0) p1=p2;
760: gnegz(p1,c); avma=av1;
761: k=0; gaffsg(1,d);
762: gaffect(u,e); gaffect(v,f);
763: p3=gsub(q,b); av2=avma;
764: for(;;)
765: {
766: w=gadd(gmul(gaddsg(k,a),u),gmul(gaddsg(-k,p3),v));
767: k++; gdivgsz(gmul(q,v),k,u);
768: gdivgsz(w,k,v);
769: gmulz(d,c,d);
770: gaddz(e,gmul(d,u),e); p1=gmul(d,v);
771: gaddz(f,p1,f);
772: if (gexpo(p1)-gexpo(f) <= lbin) break;
773: avma=av2;
774: }
775: p1=u; u=e; e=p1;
776: p1=v; v=f; f=p1;
777: gmulz(q,gadd(gun,c),q);
778: p1=subrr(q,r); ex=expo(p1); avma=av1;
779: }
780: while (ex>lbin);
781: }
782: else
783: {
784: zf=gpui(x,gneg_i(a),l1);
785: zz=gdivsg(-1,x); s=gun;
786: for (k=n-1; k>=0; k--)
787: {
788: p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
789: s=gadd(gun,gmul(p1,s));
790: }
791: u = gmul(s,zf);
792: }
793: gaffect(u,y); avma=av; return y;
794: }
795:
796: GEN
797: kbessel2(GEN nu, GEN x, long prec)
798: {
1.2 ! noro 799: gpmem_t av = avma, tetpil;
1.1 noro 800: GEN p1,p2,x2,a,pitemp;
801:
802: if (typ(x)==t_REAL) prec = lg(x);
803: x2=gshift(x,1); pitemp=mppi(prec);
804: if (gcmp0(gimag(nu))) a=cgetr(prec); else a=cgetc(prec);
805: gaddz(gun,gshift(nu,1),a);
806: p1=hyperu(gshift(a,-1),a,x2,prec);
807: p1=gmul(gmul(p1,gpui(x2,nu,prec)),mpsqrt(pitemp));
808: p2=gexp(x,prec); tetpil=avma;
809: return gerepile(av,tetpil,gdiv(p1,p2));
810: }
811:
812: GEN
813: incgam(GEN s, GEN x, long prec)
814: {
815: GEN p1,z = cgetr(prec);
1.2 ! noro 816: gpmem_t av = avma;
1.1 noro 817:
818: if (gcmp0(x)) return ggamma(s,prec);
819: if (typ(x)!=t_REAL) { gaffect(x,z); x=z; }
820: if (gcmp(subrs(x,1),s) > 0 || gsigne(greal(s)) <= 0)
821: p1 = incgam2(s,x,prec);
822: else
823: p1 = gsub(ggamma(s,prec), incgam3(s,x,prec));
824: affrr(p1,z); avma = av; return z;
825: }
826:
827: /* = incgam2(0, x, prec). typ(x) = t_REAL. Optimized for eint1 */
828: static GEN
829: incgam2_0(GEN x)
830: {
831: long l,n,i;
832: GEN p1;
833: double m,mx;
834:
835: l = lg(x); mx = rtodbl(x);
836: m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
837: p1 = divsr(-n, addsr(n<<1,x));
838: for (i=n-1; i >= 1; i--)
839: p1 = divsr(-i, addrr(addsr(i<<1,x), mulsr(i,p1)));
840: return mulrr(divrr(mpexp(negr(x)), x), addrr(realun(l),p1));
841: }
842:
843: GEN
844: incgam1(GEN a, GEN x, long prec)
845: {
846: GEN p2,p3,y, z = cgetr(prec);
1.2 ! noro 847: long l, n, i;
! 848: gpmem_t av=avma, av1;
1.1 noro 849: double m,mx;
850:
851: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
852: l=lg(x); mx=rtodbl(x);
853: m=(long) bit_accuracy(l)*LOG2; n=(long)(m/(log(m)-(1+log(mx))));
854: p2 = cgetr(l); affrr(addir(gun,gsub(x,a)), p2);
855: p3 = subrs(p2, n+1); av1 = avma;
856: for (i=n; i>=1; i--)
857: {
858: affrr(addrr(subrs(p2,i), divrr(mulsr(i,x),p3)), p3);
859: avma = av1;
860: }
861: y = mulrr(mpexp(negr(x)),gpui(x,a,prec));
862: affrr(divrr(y,p3), z);
863: avma = av; return z;
864: }
865:
866: /* assume x != 0 */
867: GEN
868: incgam2(GEN a, GEN x, long prec)
869: {
870: GEN b,p1,p2,p3,y, z = cgetr(prec);
1.2 ! noro 871: long l, n, i;
! 872: gpmem_t av = avma, av1;
1.1 noro 873: double m,mx;
874:
875: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
876: if (gcmp0(a)) { affrr(incgam2_0(x), z); avma = av; return z; }
877: l=lg(x); mx=rtodbl(x);
878: m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
879: i = typ(a);
880: if (i == t_REAL) b = addsr(-1,a);
881: else
882: { /* keep b integral : final powering more efficient */
883: gaffect(a,p1=cgetr(prec));
884: b = (i == t_INT)? addsi(-1,a): addsr(-1,p1);
885: a = p1;
886: }
887: p2 = cgetr(l); gaffect(subrr(x,a),p2);
888: p3 = divrr(addsr(-n,a), addrs(p2,n<<1));
889: av1 = avma;
890: for (i=n-1; i>=1; i--)
891: {
892: affrr(divrr(addsr(-i,a), addrr(addrs(p2,i<<1),mulsr(i,p3))), p3);
893: avma = av1;
894: }
895: y = gmul(mpexp(negr(x)), gpui(x,b,prec));
896: affrr(mulrr(y, addsr(1,p3)), z);
897: avma = av; return z;
898: }
899:
900: GEN
901: incgam3(GEN s, GEN x, long prec)
902: {
903: GEN b,p1,p2,p3,y, z = cgetr(prec);
1.2 ! noro 904: long l, n, i;
! 905: gpmem_t av = avma, av1;
1.1 noro 906:
907: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
908: l=lg(x); n = -bit_accuracy(l)-1;
909: p3 = realun(l);
910: p2 = rcopy(p3);
911: i = typ(s);
912: if (i == t_REAL) b = s;
913: else
914: {
915: gaffect(s,p1=cgetr(prec));
916: b = (i == t_INT)? s: p1;
917: s = p1;
918: }
919: if (signe(s) <= 0)
920: {
921: p1 = gcvtoi(s, &i);
922: if (i < 5 - bit_accuracy(prec))
923: err(talker,"negative argument too close to an integer in incgamc");
924: }
925: av1 = avma;
926: for (i=1; expo(p2) >= n; i++)
927: {
928: affrr(divrr(mulrr(x,p2), addsr(i,s)), p2);
929: affrr(addrr(p2,p3), p3); avma = av1;
930: }
931: y = gdiv(mulrr(mpexp(negr(x)), gpui(x,b,prec)), s);
932: affrr(mulrr(y,p3), z);
933: avma = av; return z;
934: }
935:
936: /* Assumes that g=gamma(a,prec). Unchecked */
937: GEN
938: incgam4(GEN a, GEN x, GEN g, long prec)
939: {
940: GEN p1, z = cgetr(prec);
1.2 ! noro 941: gpmem_t av = avma;
1.1 noro 942:
943: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
944: if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0)
945: p1 = incgam2(a,x,prec);
946: else
947: p1 = gsub(g, incgam3(a,x,prec));
948: affrr(p1, z);
949: avma = av; return z;
950: }
951:
952: GEN
953: incgam0(GEN a, GEN x, GEN z,long prec)
954: {
955: return z? incgam4(a,x,z,prec): incgam(a,x,prec);
956: }
957:
958: GEN
959: eint1(GEN x, long prec)
960: {
1.2 ! noro 961: long l, n, i;
! 962: gpmem_t av = avma, tetpil;
1.1 noro 963: GEN p1,p2,p3,p4,run,y;
964:
965: if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;}
966: if (signe(x) >= 0)
967: {
968: if (expo(x) >= 4)
969: return gerepileupto(av, incgam2_0(x));
970:
971: l = lg(x);
972: n = -bit_accuracy(l)-1;
973:
974: run = realun(l);
975: p4 = p3 = p2 = p1 = run;
976: for (i=2; expo(p2)>=n; i++)
977: {
978: p1 = addrr(p1,divrs(run,i)); /* p1 = sum_{i=1} 1/i */
979: p4 = divrs(mulrr(x,p4),i); /* p4 = sum_{i=1} x^(i-1)/i */
980: p2 = mulrr(p4,p1);
981: p3 = addrr(p2,p3);
982: }
983: p3 = mulrr(x,mulrr(mpexp(negr(x)),p3));
984: p1 = addrr(mplog(x), mpeuler(l));
985: return gerepileupto(av, subrr(p3,p1));
986: }
987: else
988: { /* written by Manfred Radimersky */
989: l = lg(x);
990: n = bit_accuracy(l);
991: /* IS: line split to avoid a Workshop cc-5.0 bug (Sun BugID #4254959) */
992: n = 3 * n / 4;
993: y = negr(x);
994: if(gcmpgs(y, n) < 0) {
995: p1 = p2 = p3 = y;
996: p4 = gzero;
997: i = 2;
998: while(gcmp(p3, p4)) {
999: p4 = p3;
1000: p1 = gmul(p1, gdivgs(y, i));
1001: p2 = gdivgs(p1, i);
1002: p3 = gadd(p3, p2);
1003: i++;
1004: }
1005: p1 = gadd(mplog(y), mpeuler(l));
1006: y = gadd(p3, p1);
1007: } else {
1008: p1 = gdivsg(1, y);
1009: p2 = realun(l);
1010: p3 = p2;
1011: p4 = gzero;
1012: i = 1;
1013: while(gcmp(p3, p4)) {
1014: p4 = p3;
1015: p2 = gmulgs(p2, i);
1016: p2 = gmul(p2, p1);
1017: p3 = gadd(p3, p2);
1018: i++;
1019: }
1020: p1 = gdiv(mpexp(y), y);
1021: y = gmul(p3, p1);
1022: }
1023: tetpil = avma;
1024: y = gerepile(av, tetpil, negr(y));
1025: }
1026: return y;
1027: }
1028:
1029: GEN
1030: veceint1(GEN C, GEN nmax, long prec)
1031: {
1.2 ! noro 1032: long k, n, nstop, i, cd, nmin, G, a, chkpoint;
! 1033: gpmem_t av, av1;
1.1 noro 1034: GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit;
1035:
1036: if (!nmax) return eint1(C,prec);
1037:
1038: if (signe(nmax)<=0) return cgetg(1,t_VEC);
1039: if (DEBUGLEVEL>1) fprintferr("Entering veceint1:\n");
1040: if (typ(C) != t_REAL || lg(C) > prec)
1041: { y = cgetr(prec); gaffect(C,y); C = y; }
1042: if (signe(C) <= 0) err(talker,"negative or zero constant in veceint1");
1043:
1044: G = -bit_accuracy(prec);
1045: n=itos(nmax); y=cgetg(n+1,t_VEC);
1046: for(i=1; i<=n; i++) y[i]=lgetr(prec);
1047: av=avma;
1048:
1049: nstop = itos(gceil(divsr(4,C)));
1050: if (nstop<1) nstop=1;
1051: if (nstop>n) nstop=n;
1052: nmin=n-10; if (nmin<nstop) nmin=nstop;
1053: if(DEBUGLEVEL>1) fprintferr("nstop = %ld\n",nstop);
1054:
1055: e1=mpexp(negr(mulsr(n,C)));
1056: e2=mpexp(mulsr(10,C));
1057: unr = realun(prec);
1058: zeror=realzero(prec); deninit=negr(unr);
1059: f=cgetg(3,t_COL); M2=cgetg(3,t_VEC); av1=avma;
1060:
1061: F0=(GEN)y[n]; chkpoint = n;
1062: affrr(eint1(mulsr(n,C),prec), F0);
1063: do
1064: {
1065: if (DEBUGLEVEL>1 && n < chkpoint)
1066: { fprintferr("%ld ",n) ; chkpoint -= (itos(nmax) / 20); }
1067: minvn=divrs(unr,-n); mcn=mulrr(C,minvn);
1068: M2[1] = (long)zeror; M2[2] = lsubrr(minvn,C);
1069: f[1]=(long)zeror; f[2]=ldivrs(e1,-n);
1070: affrr(mulrr(e1,e2), e1);
1071: vdiff=cgetg(2,t_VEC); vdiff[1]=f[2];
1072: for (cd=a=1,n--; n>=nmin; n--,a++)
1073: {
1074: F = F0;
1075: ap = stoi(a); den = deninit;
1076: for (k=1;;)
1077: {
1078: if (k>cd)
1079: {
1080: cd++; p1 = (GEN)f[2];
1081: f[2] = lmul(M2,f);
1082: f[1] = (long)p1;
1083: M2[1] = laddrr((GEN)M2[1],mcn);
1084: M2[2] = laddrr((GEN)M2[2],minvn);
1085: vdiff = concatsp(vdiff,(GEN)f[2]);
1086: }
1087: p1 = mulrr(mulri(den,ap),(GEN)vdiff[k]);
1088: if (expo(p1) < G) { affrr(F,(GEN)y[n]); break; }
1089: F = addrr(F,p1); ap = mulis(ap,a);
1090: k++; den = divrs(den,-k);
1091: }
1092: }
1093: avma=av1; n++; F0=(GEN)y[n]; nmin -= 10;
1094: if (nmin < nstop) nmin=nstop;
1095: }
1096: while(n>nstop);
1097: for(i=1; i<=nstop; i++)
1098: affrr(eint1(mulsr(i,C),prec), (GEN)y[i]);
1099: if (DEBUGLEVEL>1) fprintferr("\n");
1100: avma=av; return y;
1101: }
1102:
1103: GEN
1104: gerfc(GEN x, long prec)
1105: {
1.2 ! noro 1106: gpmem_t av=avma;
1.1 noro 1107: GEN p1,p2;
1108:
1109: if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; }
1110: av = avma; p1 = incgam(ghalf,gsqr(x),prec);
1111: p2 = mpsqrt(mppi(lg(x)));
1112: p1 = divrr(p1,p2);
1113: if (signe(x) < 0) p1 = subsr(2,p1);
1114: return gerepileupto(av,p1);
1115: }
1116:
1117: /***********************************************************************/
1118: /** **/
1119: /** FONCTION ZETA DE RIEMANN **/
1120: /** **/
1121: /***********************************************************************/
1122:
1123: #if 0
1124: static GEN
1125: czeta(GEN s, long prec)
1126: {
1.2 ! noro 1127: long n, p, n1, flag1, flag2, i, i2;
! 1128: gpmem_t av;
1.1 noro 1129: double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf;
1130: GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp;
1131:
1132: i=precision(s); if (i) prec = i;
1133: if (typ(s)==t_COMPLEX)
1134: {
1135: res = cgetc(prec); av=avma;
1136: p1 = cgetc(prec+1);
1137: gaffect(s,p1); s=p1; sig=(GEN)s[1];
1138: }
1139: else
1140: {
1141: res = cgetr(prec); av=avma;
1142: p1 = cgetr(prec+1); affrr(s,p1); sig = s = p1;
1143: }
1144:
1145: if (signe(sig)>0 && expo(sig)>-2) flag1 = 0;
1146: else
1147: {
1148: if (gcmp0(gimag(s)) && gcmp0(gfrac(gmul2n(sig,-1))))
1149: {
1150: if (gcmp0(sig)) gaffect(gneg_i(ghalf),res); else gaffsg(0,res);
1151: avma=av; return res;
1152: }
1153: flag1 = 1; s = gsub(gun,s); sig = greal(s);
1154: }
1155: ssig=rtodbl(sig); st=fabs(rtodbl(gimag(s))); maxbeta = pow(3.0,-2.5);
1156: if (st)
1157: {
1158: ns = ssig*ssig + st*st;
1159: alpha=pariC2*(prec-2)-0.39-0.5*(ssig-1.0)*log(ns)-log(ssig)+ssig*2*pariC1;
1160: beta=(alpha+ssig)/st-atan(ssig/st);
1161: if (beta<=0)
1162: {
1163: if (ssig>=1.0)
1164: {
1165: p=0; sn=sqrt(ns);
1166: n=(long)(ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig)));
1167: }
1168: else
1169: {
1170: p=1; sn=ssig+1; n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
1171: }
1172: }
1173: else
1174: {
1175: if (beta<maxbeta) xinf=beta+pow(3*beta,1.0/3.0);
1176: else
1177: {
1178: double eps=0.0087, x00 = beta+PI/2.0, y00,x11;
1179: for(;;)
1180: {
1181: y00=x00*x00; x11=(beta+atan(x00))*(1+y00)/y00-1/x00;
1182: if (x00-x11 < eps) break;
1183: x00 = x11;
1184: }
1185: xinf=x11;
1186: }
1187: sp=1.0-ssig+st*xinf;
1188: if (sp>0)
1189: {
1190: p=(long)ceil(sp/2.0); sn=ssig+2*p-1;
1191: n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
1192: }
1193: else
1194: {
1195: p=0; sn=sqrt(ns);
1196: n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
1197: }
1198: }
1199: }
1200: else
1201: {
1202: beta=pariC2*(prec-2)+0.61+ssig*2*pariC1-ssig*log(ssig);
1203: if (beta>0)
1204: {
1205: p=(long)ceil(beta/2.0); sn=ssig+2*p-1;
1206: n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
1207: }
1208: else
1209: {
1210: p=0; sn=sqrt(ssig*ssig+st*st);
1211: n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
1212: }
1213: }
1214: if (n < 46340) n1=n*n; else n1=0;
1215: y=gun; ms=gneg_i(s); p1=cgetr(prec+1); p2=gun;
1216: for (i=2; i<=n; i++)
1217: {
1218: affsr(i,p1); p2 = gexp(gmul(ms,mplog(p1)), prec+1);
1219: y = gadd(y,p2);
1220: }
1221: flag2 = (2*p < 46340);
1222: mpbern(p,prec+1); p31=cgetr(prec+1); z=gzero;
1223: for (i=p; i>=1; i--)
1224: {
1225: i2=i<<1;
1226: p1=gmul(gaddsg(i2-1,s),gaddsg(i2,s));
1227: p1=flag2? gdivgs(p1,i2*(i2+1)): gdivgs(gdivgs(p1,i2),i2+1);
1228: p1=n1? gdivgs(p1,n1): gdivgs(gdivgs(p1,n),n);
1229: p3 = bern(i);
1230: if (bernzone[2]>prec+1) { affrr(p3,p31); p3=p31; }
1231: z=gadd(divrs(p3,i),gmul(p1,z));
1232: }
1233: p1=gsub(gdivsg(n,gsubgs(s,1)),ghalf);
1234: z=gmul(gadd(p1,gmul(s,gdivgs(z,n<<1))),p2);
1235: y = gadd(y,z);
1236: if (flag1)
1237: {
1238: pitemp=mppi(prec+1); setexpo(pitemp,2);
1239: y=gmul(gmul(y,ggamma(s,prec+1)),gpui(pitemp,ms,prec+1));
1240: setexpo(pitemp,0);
1241: y=gmul2n(gmul(y,gcos(gmul(pitemp,s),prec+1)),1);
1242: }
1243: gaffect(y,res); avma = av; return res;
1244: }
1245:
1246: /* y = binomial(n,k-2). Return binomial(n,k) */
1247: static GEN
1248: next_bin(GEN y, long n, long k)
1249: {
1250: y = divrs(mulrs(y, n-k+2), k-1);
1251: return divrs(mulrs(y, n-k+1), k);
1252: }
1253:
1254: static GEN
1255: izeta(long k, long prec)
1256: {
1.2 ! noro 1257: long kk, n, li;
! 1258: gpmem_t av=avma, av2, limit;
1.1 noro 1259: GEN y,p1,pi2,qn,z,q,binom;
1260:
1261: /* treat trivial cases */
1262: if (!k) return gneg(ghalf);
1263: if (k < 0)
1264: {
1265: if ((k&1) == 0) return gzero;
1266: y = bernreal(1-k,prec);
1267: return gerepileuptoleaf(av, divrs(y,k-1));
1268: }
1269: if (k > bit_accuracy(prec)+1) return realun(prec);
1.2 ! noro 1270: pi2 = Pi2n(1, prec);
1.1 noro 1271: if ((k&1) == 0)
1272: {
1273: p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec)));
1274: y = divrr(p1, mpfactr(k,prec)); setexpo(y,expo(y)-1);
1275: return gerepileuptoleaf(av, y);
1276: }
1277: /* k > 1 odd */
1278: binom = realun(prec+1);
1279: q = mpexp(pi2); kk = k+1; /* >= 4 */
1280: li = -(1+bit_accuracy(prec));
1281: y = NULL; /* gcc -Wall */
1282: if ((k&3)==3)
1283: {
1284: for (n=0; n <= kk>>1; n+=2)
1285: {
1286: p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
1287: if (n) { binom = next_bin(binom,kk,n); setlg(binom,prec+1); }
1288: p1 = mulrr(binom,p1);
1289: if (n == kk>>1) setexpo(p1, expo(p1)-1);
1290: if ((n>>1)&1) setsigne(p1,-signe(p1));
1291: y = n? addrr(y,p1): p1;
1292: }
1293: y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y);
1294:
1295: av2 = avma; limit = stack_lim(av2,1);
1.2 ! noro 1296: qn = gsqr(q); z = ginv( addrs(q,-1) );
1.1 noro 1297: for (n=2; ; n++)
1298: {
1.2 ! noro 1299: p1 = ginv( mulir(gpuigs(stoi(n),k),addrs(qn,-1)) );
1.1 noro 1300:
1301: z = addrr(z,p1); if (expo(p1)< li) break;
1302: qn = mulrr(qn,q);
1303: if (low_stack(limit,stack_lim(av2,1)))
1304: {
1305: if (DEBUGMEM>1) err(warnmem,"izeta");
1.2 ! noro 1306: gerepileall(av2,2, &z, &qn);
1.1 noro 1307: }
1308: }
1309: setexpo(z,expo(z)+1);
1310: y = addrr(y,z); setsigne(y,-signe(y));
1311: }
1312: else
1313: {
1314: GEN p2 = divrs(pi2, k-1);
1315: for (n=0; n <= k>>1; n+=2)
1316: {
1317: p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
1318: if (n) binom = next_bin(binom,kk,n);
1319: p1 = mulrr(binom,p1);
1320: p1 = mulsr(kk-(n<<1),p1);
1321: if ((n>>1)&1) setsigne(p1,-signe(p1));
1322: y = n? addrr(y,p1): p1;
1323: }
1324: y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y);
1325: y = divrs(y,k-1);
1326: av2 = avma; limit = stack_lim(av2,1);
1327: qn = q; z=gzero;
1328: for (n=1; ; n++)
1329: {
1330: p1=mulir(gpuigs(stoi(n),k),gsqr(addrs(qn,-1)));
1331: p1=divrr(addrs(mulrr(qn,addsr(1,mulsr(n<<1,p2))),-1),p1);
1332:
1333: z=addrr(z,p1); if (expo(p1) < li) break;
1334: qn=mulrr(qn,q);
1335: if (low_stack(limit,stack_lim(av2,1)))
1336: {
1337: if (DEBUGMEM>1) err(warnmem,"izeta");
1.2 ! noro 1338: gerepileall(av2,2, &z, &qn);
1.1 noro 1339: }
1340: }
1341: setexpo(z,expo(z)+1);
1342: y = subrr(y,z);
1343: }
1344: return gerepileuptoleaf(av, y);
1345: }
1346: #else
1347:
1348: /* return x^n, assume n > 0 */
1349: static long
1350: pows(long x, long n)
1351: {
1352: long i, y = x;
1353: for (i=1; i<n; i++) y *= x;
1354: return y;
1355: }
1356:
1357: /* return n^-s, n > 1 odd. tab[q] := q^-s, q prime power */
1358: static GEN
1.2 ! noro 1359: n_s(ulong n, GEN *tab)
1.1 noro 1360: {
1361: byteptr d = diffptr + 2;
1362: GEN x = NULL;
1363: long p, e;
1364:
1.2 ! noro 1365: for (p = 3; n > 1; )
1.1 noro 1366: {
1367: e = svaluation(n, p, &n);
1368: if (e)
1369: {
1370: GEN y = tab[pows(p,e)];
1371: if (!x) x = y; else x = gmul(x,y);
1372: }
1.2 ! noro 1373: NEXT_PRIME_VIADIFF_CHECK(p,d);
1.1 noro 1374: }
1375: return x;
1376: }
1377:
1378: GEN czeta(GEN s0, long prec);
1379:
1380: /* assume k != 1 */
1381: static GEN
1382: izeta(long k, long prec)
1383: {
1.2 ! noro 1384: gpmem_t av = avma;
1.1 noro 1385: GEN y,p1,pi2;
1386:
1387: /* treat trivial cases */
1388: if (!k) { y = realun(prec); setexpo(y,-1); setsigne(y,-1); return y; }
1389: if (k < 0)
1390: {
1391: if ((k&1) == 0) return gzero;
1392: y = bernreal(1-k,prec);
1393: return gerepileuptoleaf(av, divrs(y,k-1));
1394: }
1395: if (k > bit_accuracy(prec)+1) return realun(prec);
1396: if ((k&1) == 0)
1397: {
1398: pi2 = mppi(prec); setexpo(pi2,2); /* 2Pi */
1399: p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec)));
1400: y = divrr(p1, mpfactr(k,prec)); setexpo(y,expo(y)-1);
1401: return gerepileuptoleaf(av, y);
1402: }
1403: /* k > 1 odd */
1404: return czeta(stoi(k), prec);
1405: }
1406:
1407: /* s0 a t_INT, t_REAL or t_COMPLEX.
1408: * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd)
1409: * */
1410: GEN
1411: czeta(GEN s0, long prec)
1412: {
1413: GEN s, u, a, y, res, tes, sig, invn2, p1, unr;
1414: GEN sim, ms, s1, s2, s3, s4, s5, *tab, tabn;
1.2 ! noro 1415: long p, i, sqn, nn, lim, lim2, ct;
! 1416: gpmem_t av, av2 = avma, avlim;
1.1 noro 1417: int funeq = 0;
1418: byteptr d;
1419:
1.2 ! noro 1420: if (DEBUGLEVEL>2) (void)timer2();
1.1 noro 1421: s = trans_fix_arg(&prec,&s0,&sig,&av,&res);
1422: if (gcmp0(s)) { y = gneg(ghalf); goto END; }
1423: if (signe(sig) <= 0 || expo(sig) < -1)
1424: { /* s <--> 1-s */
1425: if (typ(s0) == t_INT)
1426: {
1427: p = itos(s0); avma = av2;
1428: return izeta(p,prec);
1429: }
1430: funeq = 1; s = gsub(gun, s); sig = greal(s);
1431: }
1432: if (gcmp(sig, stoi(bit_accuracy(prec) + 1)) > 0) { y = gun; goto END; }
1433:
1434: { /* find "optimal" parameters [lim, nn] */
1435: double ssig = rtodbl(sig);
1436: double st = rtodbl(gimag(s));
1437: double ns = dnorm(ssig,st), l,l2;
1438: long la = 1;
1439:
1440: if (typ(s0) == t_INT)
1441: {
1442: long ss = itos(s0);
1443: switch(ss)
1444: { /* should depend on prec ? */
1445: case 3: la = 6; break;
1446: default: la = 3; break;
1447: }
1448: }
1449: if (dnorm(ssig-1,st) < 0.1) /* |s - 1| < 0.1 */
1450: l2 = -(ssig - 0.5);
1451: else
1452: { /* l2 = Re( (s - 1/2) log (s-1) ) */
1453: double rlog, ilog; /* log(s-1) */
1454: dcxlog(ssig-1,st, &rlog,&ilog);
1455: l2 = (ssig - 0.5)*rlog - st*ilog;
1456: }
1457: l = (pariC2*(prec-2) - l2 + ssig*2*pariC1) / (2. * (1.+ log((double)la)));
1458: l2 = sqrt(ns)/2;
1459: if (l < l2) l = l2;
1460: lim = (long) ceil(l); if (lim < 2) lim = 2;
1461: l2 = (lim+ssig/2.-.25);
1462: nn = 1 + (long)ceil( sqrt(l2*l2 + st*st/4) * la / PI );
1.2 ! noro 1463: if (DEBUGLEVEL>2) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn);
1.1 noro 1464: if ((ulong)nn >= maxprime()) err(primer1);
1465: }
1466: prec++; unr = realun(prec); /* one extra word of precision */
1467:
1468: tab = (GEN*)cgetg(nn, t_VEC); /* table of q^(-s), q = p^e */
1469: d = diffptr + 1;
1470: if (typ(s0) == t_INT)
1471: { /* no explog for 1/p^s */
1.2 ! noro 1472: for (p=2; p < nn;) {
1.1 noro 1473: tab[p] = divrr(unr, rpowsi(p, s0, prec));
1.2 ! noro 1474: NEXT_PRIME_VIADIFF_CHECK(p,d);
! 1475: }
1.1 noro 1476: a = divrr(unr, rpowsi(nn, s0, prec));
1477: }
1478: else
1479: { /* general case */
1480: ms = gneg(s); p1 = cgetr(prec);
1.2 ! noro 1481: for (p=2; p < nn;)
1.1 noro 1482: {
1483: affsr(p, p1);
1484: tab[p] = gexp(gmul(ms, mplog(p1)), prec);
1.2 ! noro 1485: NEXT_PRIME_VIADIFF_CHECK(p,d);
1.1 noro 1486: }
1487: affsr(nn,p1);
1488: a = gexp(gmul(ms, mplog(p1)), prec);
1489: }
1490: sqn = (long)sqrt(nn-1.);
1491: d = diffptr + 2; /* fill in odd prime powers */
1.2 ! noro 1492: for (p=3; p <= sqn; )
1.1 noro 1493: {
1494: ulong oldq = p, q = p*p;
1495: while (q<(ulong)nn) { tab[q] = gmul(tab[p], tab[oldq]); oldq = q; q *= p; }
1.2 ! noro 1496: NEXT_PRIME_VIADIFF_CHECK(p,d);
1.1 noro 1497: }
1.2 ! noro 1498: if (DEBUGLEVEL>2) msgtimer("tab[q^-s] from 1 to N-1");
1.1 noro 1499:
1500: tabn = cgetg(nn, t_VECSMALL); ct = 0;
1501: for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1;
1502: sim = y = unr;
1503: for (i=ct; i > 1; i--)
1504: {
1.2 ! noro 1505: long j;
! 1506: gpmem_t av2 = avma;
1.1 noro 1507: for (j=tabn[i]+1; j<=tabn[i-1]; j++)
1508: sim = gadd(sim, n_s(2*j+1, tab));
1509: sim = gerepileupto(av2, sim);
1510: y = gadd(sim, gmul(tab[2],y));
1511: }
1512: y = gadd(y, gmul2n(a,-1));
1.2 ! noro 1513: if (DEBUGLEVEL>2) msgtimer("sum from 1 to N-1");
1.1 noro 1514:
1515: invn2 = divrs(unr, nn*nn); lim2 = lim<<1;
1516: tes = bernreal(lim2, prec);
1517: if (typ(s0) == t_INT)
1518: {
1519: av2 = avma; avlim = stack_lim(av2,3);
1520: for (i=lim2-2; i>=2; i-=2)
1521: { /* using single prec (when (s0 + i) < 2^31) not faster (even at \p28) */
1522: u = mulri(mulrr(tes,invn2), mulii(addsi(i,s0), addsi(i-1,s0)));
1523: tes = addrr(bernreal(i,prec), divrs2_safe(u, i+1)); /* u / (i+1)(i+2) */
1524: if (low_stack(avlim,stack_lim(av2,3)))
1525: {
1526: if(DEBUGMEM>1) err(warnmem,"czeta");
1527: tes = gerepileuptoleaf(av2, tes);
1528: }
1529: }
1530: u = gmul(gmul(tes,invn2), gmul2n(mulii(s0, addsi(-1,s0)), -1));
1531: tes = gmulsg(nn, gaddsg(1, u));
1532: }
1533: else /* typ(s0) != t_INT */
1534: {
1535: s1 = gsub(gmul2n(s,1), unr);
1536: s2 = gmul(s, gsub(s,unr));
1537: s3 = gmul2n(invn2,3);
1538: av2 = avma; avlim = stack_lim(av2,3);
1539: s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
1540: s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
1541: for (i = lim2-2; i>=2; i -= 2)
1542: {
1543: s5 = gsub(s5, s4);
1544: s4 = gsub(s4, s3);
1545: tes = gadd(bernreal(i,prec), gdivgs(gmul(s5,tes), (i+1)*(i+2)));
1546: if (low_stack(avlim,stack_lim(av2,3)))
1547: {
1548: GEN *gptr[3]; gptr[0]=&tes; gptr[1]=&s5; gptr[2]=&s4;
1549: if(DEBUGMEM>1) err(warnmem,"czeta");
1550: gerepilemany(av2,gptr,3);
1551: }
1552: }
1553: u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
1554: tes = gmulsg(nn, gaddsg(1, u));
1555: }
1.2 ! noro 1556: if (DEBUGLEVEL>2) msgtimer("Bernoulli sum");
1.1 noro 1557: /* y += tes a / (s-1) */
1558: y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr))));
1559:
1560: END:
1561: if (funeq)
1562: {
1563: GEN pitemp = mppi(prec); setexpo(pitemp,2); /* 2Pi */
1.2 ! noro 1564: y = gmul(gmul(y, ggamma(gprec_w(s,prec-1),prec)),
! 1565: gpui(pitemp,gneg(s),prec));
1.1 noro 1566: setexpo(pitemp,0); /* Pi/2 */
1567: y = gmul2n(gmul(y, gcos(gmul(pitemp,s),prec)), 1);
1568: }
1569: gaffect(y,res); avma = av; return res;
1570: }
1571: #endif
1572:
1573: GEN
1574: gzeta(GEN x, long prec)
1575: {
1576: if (gcmp1(x)) err(talker,"argument equal to one in zeta");
1577: switch(typ(x))
1578: {
1579: case t_INT:
1580: if (is_bigint(x))
1581: {
1582: if (signe(x) > 0) return realun(prec);
1583: if (signe(x) < 0 && mod2(x) == 0) return realzero(prec);
1584: }
1585: return izeta(itos(x),prec);
1586:
1587: case t_REAL: case t_COMPLEX:
1588: return czeta(x,prec);
1589:
1590: case t_INTMOD: case t_PADIC:
1591: err(typeer,"gzeta");
1592: case t_SER:
1593: err(impl,"zeta of power series");
1594: }
1595: return transc(gzeta,x,prec);
1596: }
1597:
1598: void
1599: gzetaz(GEN x, GEN y)
1600: {
1.2 ! noro 1601: long prec = precision(y);
! 1602: gpmem_t av=avma;
1.1 noro 1603:
1604: if (!prec) err(infprecer,"gzetaz");
1605: gaffect(gzeta(x,prec),y); avma=av;
1606: }
1607:
1608: /***********************************************************************/
1609: /** **/
1610: /** FONCTIONS POLYLOGARITHME **/
1611: /** **/
1612: /***********************************************************************/
1613:
1614: /* validity domain contains .005 < |x| < 230
1615: * Li_m(x = e^z) = sum_n=0 zeta(m-n) z^n / n!
1616: * with zeta(1) := H_m - log(-z) */
1617: static GEN
1618: cxpolylog(long m, GEN x, long prec)
1619: {
1.2 ! noro 1620: long li, i, n, bern_upto;
! 1621: gpmem_t av=avma;
1.1 noro 1622: GEN p1,z,h,q,s;
1623:
1624: if (gcmp1(x)) return izeta(m,prec);
1625: z=glog(x,prec); h=gneg_i(glog(gneg_i(z),prec));
1626: for (i=1; i<m; i++) h = gadd(h,ginv(stoi(i)));
1627:
1628: bern_upto=m+50; mpbern(bern_upto,prec);
1629: q=gun; s=izeta(m,prec);
1630: for (n=1; n<=m+1; n++)
1631: {
1632: q=gdivgs(gmul(q,z),n);
1633: s=gadd(s,gmul((n==(m-1))? h: izeta(m-n,prec),q));
1634: }
1635:
1636: n = m+3; z=gsqr(z); li = -(bit_accuracy(prec)+1);
1637:
1638: for(;;)
1639: {
1640: q = gdivgs(gmul(q,z),(n-1)*n);
1641: p1 = gmul(izeta(m-n,prec),q);
1642: s = gadd(s,p1);
1643: if (gexpo(p1) < li) break;
1644: n+=2;
1645: if (n>=bern_upto) { bern_upto += 50; mpbern(bern_upto,prec); }
1646: }
1647: return gerepileupto(av,s);
1648: }
1649:
1650: GEN
1651: polylog(long m, GEN x, long prec)
1652: {
1.2 ! noro 1653: long l, e, i, G, sx;
! 1654: gpmem_t av, av1, limpile;
1.1 noro 1655: GEN X,z,p1,p2,n,y,logx;
1656:
1657: if (m<0) err(talker,"negative index in polylog");
1658: if (!m) return gneg(ghalf);
1659: if (gcmp0(x)) return gcopy(x);
1660: av = avma;
1661: if (m==1)
1662: return gerepileupto(av, gneg(glog(gsub(gun,x), prec)));
1663:
1664: l = precision(x);
1665: if (!l) { l=prec; x=gmul(x, realun(l)); }
1666: e = gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
1667: X = (e > 0)? ginv(x): x;
1668: G = -bit_accuracy(l);
1669: n = icopy(gun);
1670: av1=avma; limpile=stack_lim(av1,1);
1671: y = p1 = X;
1672: for (i=2; ; i++)
1673: {
1674: n[2] = i; p1 = gmul(X,p1); p2 = gdiv(p1,gpuigs(n,m));
1675: y = gadd(y,p2);
1676: if (gexpo(p2) <= G) break;
1677:
1678: if (low_stack(limpile, stack_lim(av1,1)))
1679: { GEN *gptr[2]; gptr[0]=&y; gptr[1]=&p1;
1680: if(DEBUGMEM>1) err(warnmem,"polylog");
1681: gerepilemany(av1,gptr,2);
1682: }
1683: }
1684: if (e < 0) return gerepileupto(av, y);
1685:
1686: sx = gsigne(gimag(x));
1687: if (!sx)
1688: {
1689: if (m&1) sx = gsigne(gsub(gun,greal(x)));
1690: else sx = - gsigne(greal(x));
1691: }
1692: z = cgetg(3,t_COMPLEX);
1693: z[1] = zero;
1694: z[2] = ldivri(mppi(l), mpfact(m-1));
1695: if (sx < 0) z[2] = lnegr((GEN)z[2]);
1696: logx = glog(x,l);
1697:
1698: if (m == 2)
1699: { /* same formula as below, written more efficiently */
1700: y = gneg_i(y);
1701: p1 = gmul2n(gsqr(gsub(logx, z)), -1); /* = (log(-x))^2 / 2 */
1702: p1 = gadd(divrs(gsqr(mppi(l)), 6), p1);
1703: p1 = gneg_i(p1);
1704: }
1705: else
1706: {
1707: GEN logx2 = gsqr(logx); p1 = gneg_i(ghalf);
1708: for (i=m-2; i>=0; i-=2)
1709: p1 = gadd(izeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
1710: if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
1711: p1 = gadd(gmul2n(p1,1), gmul(z,gpuigs(logx,m-1)));
1712: }
1713:
1714: return gerepileupto(av, gadd(y,p1));
1715: }
1716:
1717: GEN
1718: dilog(GEN x, long prec)
1719: {
1720: return gpolylog(2, x, prec);
1721: }
1722:
1723: GEN
1724: polylogd0(long m, GEN x, long flag, long prec)
1725: {
1.2 ! noro 1726: long k, l, fl, m2;
! 1727: gpmem_t av;
1.1 noro 1728: GEN p1,p2,p3,y;
1729:
1730: m2=m&1; av=avma;
1731: if (gcmp0(x)) return gcopy(x);
1732: if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
1733: l=precision(x);
1734: if (!l) { l=prec; x=gmul(x,realun(l)); }
1735: p1=gabs(x,prec); fl=0;
1736: if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
1737:
1738: p1=gneg_i(glog(p1,prec)); p2=gun;
1739: y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
1740: for (k=1; k<m; k++)
1741: {
1742: p2=gdivgs(gmul(p2,p1),k);
1743: p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
1744: y=gadd(y,gmul(p2,p3));
1745: }
1746: if (m2)
1747: {
1748: if (flag)
1749: p2 = gdivgs(gmul(p2,p1),-2*m);
1750: else
1751: p2 = gdivgs(gmul(glog(gnorm(gsub(gun,x)),prec),p2),2*m);
1752: y=gadd(y,p2);
1753: }
1754: if (fl) y = gneg(y);
1755: return gerepileupto(av, y);
1756: }
1757:
1758: GEN
1759: polylogd(long m, GEN x, long prec)
1760: {
1761: return polylogd0(m,x,0,prec);
1762: }
1763:
1764: GEN
1765: polylogdold(long m, GEN x, long prec)
1766: {
1767: return polylogd0(m,x,1,prec);
1768: }
1769:
1770: GEN
1771: polylogp(long m, GEN x, long prec)
1772: {
1.2 ! noro 1773: long k, l, fl, m2;
! 1774: gpmem_t av;
1.1 noro 1775: GEN p1,y;
1776:
1777: m2=m&1; av=avma;
1778: if (gcmp0(x)) return gcopy(x);
1779: if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
1780: l=precision(x);
1781: if (!l) { l=prec; x=gmul(x,realun(l)); }
1782: p1=gabs(x,prec); fl=0;
1783: if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
1784:
1785: p1=gmul2n(glog(p1,prec),1); mpbern(m>>1,prec);
1786: y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
1787:
1788: if (m==1)
1789: {
1790: p1=gmul2n(p1,-2); y=gadd(y,p1);
1791: }
1792: else
1793: {
1794: GEN p2=gun, p3, p4, p5, p51=cgetr(prec);
1795:
1796: for (k=1; k<m; k++)
1797: {
1798: p2=gdivgs(gmul(p2,p1),k);
1799: if (!(k&1) || k==1)
1800: {
1801: if (k!=1)
1802: {
1803: p5=(GEN)bern(k>>1);
1804: if (bernzone[2]>prec) { affrr(p5,p51); p5=p51; }
1805: p4=gmul(p2,p5);
1806: }
1807: else p4=gneg_i(gmul2n(p2,-1));
1808: p3=polylog(m-k,x,prec); p3=m2?greal(p3):gimag(p3);
1809: y=gadd(y,gmul(p4,p3));
1810: }
1811: }
1812: }
1813: if (fl) y = gneg(y);
1814: return gerepileupto(av, y);
1815: }
1816:
1817: GEN
1818: gpolylog(long m, GEN x, long prec)
1819: {
1.2 ! noro 1820: long i, lx, v, n;
! 1821: gpmem_t av=avma;
1.1 noro 1822: GEN y,p1,p2;
1823:
1824: if (m<=0)
1825: {
1826: p1=polx[0]; p2=gsub(gun,p1);
1827: for (i=1; i<=(-m); i++)
1828: p1=gmul(polx[0],gadd(gmul(p2,derivpol(p1)),gmulsg(i,p1)));
1829: p1=gdiv(p1,gpuigs(p2,1-m));
1830: return gerepileupto(av, poleval(p1,x));
1831: }
1832:
1833: switch(typ(x))
1834: {
1835: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
1836: case t_COMPLEX: case t_QUAD:
1837: return polylog(m,x,prec);
1838:
1839: case t_INTMOD: case t_PADIC:
1840: err(impl, "padic polylogarithm");
1841: case t_POLMOD:
1842: p1=roots((GEN)x[1],prec); lx=lg(p1); p2=cgetg(lx,t_COL);
1843: for (i=1; i<lx; i++) p2[i]=lpoleval((GEN)x[2],(GEN)p1[i]);
1844: y=cgetg(lx,t_COL);
1845: for (i=1; i<lx; i++) y[i]=(long)polylog(m,(GEN)p2[i],prec);
1846: return gerepileupto(av, y);
1847:
1848: case t_POL: case t_RFRAC: case t_RFRACN:
1849: p1=tayl(x,gvar(x),precdl);
1850: return gerepileupto(av, gpolylog(m,p1,prec));
1851:
1852: case t_SER:
1853: if (!m) return gneg(ghalf);
1854: if (m==1)
1855: {
1856: p1=glog(gsub(gun,x),prec);
1857: return gerepileupto(av, gneg(p1));
1858: }
1.2 ! noro 1859: if (gcmp0(x)) return gcopy(x);
1.1 noro 1860: if (valp(x)<=0) err(impl,"polylog around a!=0");
1861: v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2);
1862: for (i=n; i>=1; i--)
1863: {
1864: p1=gadd(gpuigs(stoi(i),-m),y);
1865: y=gmul(x,p1);
1866: }
1867: return gerepileupto(av, y);
1868:
1869: case t_VEC: case t_COL: case t_MAT:
1870: lx=lg(x); y=cgetg(lx,typ(x));
1871: for (i=1; i<lx; i++)
1872: y[i]=(long)gpolylog(m,(GEN)x[i],prec);
1873: return y;
1874: }
1875: err(typeer,"gpolylog");
1876: return NULL; /* not reached */
1877: }
1878:
1879: void
1880: gpolylogz(long m, GEN x, GEN y)
1881: {
1.2 ! noro 1882: long prec = precision(y);
! 1883: gpmem_t av=avma;
1.1 noro 1884:
1885: if (!prec) err(infprecer,"gpolylogz");
1886: gaffect(gpolylog(m,x,prec),y); avma=av;
1887: }
1888:
1889: GEN
1890: polylog0(long m, GEN x, long flag, long prec)
1891: {
1892: switch(flag)
1893: {
1894: case 0: return gpolylog(m,x,prec);
1895: case 1: return polylogd(m,x,prec);
1896: case 2: return polylogdold(m,x,prec);
1897: case 3: return polylogp(m,x,prec);
1898: default: err(flagerr,"polylog");
1899: }
1900: return NULL; /* not reached */
1901: }
1902:
1903: static GEN
1904: qq(GEN x, long prec)
1905: {
1906: long tx=typ(x);
1907:
1908: if (tx==t_PADIC) return x;
1909: if (is_scalar_t(tx))
1910: {
1.2 ! noro 1911: long l = precision(x);
! 1912: if (!l) l = prec;
! 1913: if (tx != t_COMPLEX || gsigne((GEN)x[2]) <= 0)
1.1 noro 1914: err(talker,"argument must belong to upper half-plane");
1915:
1.2 ! noro 1916: return gexp(gmul(x, PiI2(l)), l); /* e(x) */
1.1 noro 1917: }
1.2 ! noro 1918: if (tx == t_SER) return x;
! 1919: if (tx != t_POL) err(talker,"bad argument for modular function");
1.1 noro 1920:
1921: return tayl(x,gvar(x),precdl);
1922: }
1923:
1924: static GEN
1925: inteta(GEN q)
1926: {
1927: long tx=typ(q);
1928: GEN p1,ps,qn,y;
1929:
1930: y=gun; qn=gun; ps=gun;
1931: if (tx==t_PADIC)
1932: {
1.2 ! noro 1933: if (valp(q) <= 0) err(talker,"non-positive valuation in eta");
1.1 noro 1934: for(;;)
1935: {
1936: p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
1937: y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
1938: p1=y; y=gadd(y,ps); if (gegal(p1,y)) return y;
1939: }
1940: }
1941: else
1942: {
1.2 ! noro 1943: long l, v = 0; /* gcc -Wall */
! 1944: gpmem_t av = avma, lim = stack_lim(av, 3);
1.1 noro 1945:
1946: if (is_scalar_t(tx)) l = -bit_accuracy(precision(q));
1947: else
1948: {
1.2 ! noro 1949: v = gvar(q); l = lg(q)-2; tx = 0;
! 1950: if (valp(q) <= 0) err(talker,"non-positive valuation in eta");
1.1 noro 1951: }
1952: for(;;)
1953: {
1.2 ! noro 1954: p1 = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
! 1955: /* qn = q^n
! 1956: * ps = (-1)^n q^(n(3n+1)/2)
! 1957: * p1 = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
! 1958: y = gadd(y,p1); qn = gmul(qn,q); ps = gmul(p1,qn);
! 1959: y = gadd(y,ps);
1.1 noro 1960: if (tx)
1961: { if (gexpo(ps)-gexpo(y) < l) return y; }
1962: else
1.2 ! noro 1963: { if (gval(ps,v) >= l) return y; }
1.1 noro 1964: if (low_stack(lim, stack_lim(av,3)))
1965: {
1.2 ! noro 1966: if(DEBUGMEM>1) err(warnmem,"eta");
! 1967: gerepileall(av, 3, &y, &qn, &ps);
1.1 noro 1968: }
1969: }
1970: }
1971: }
1972:
1973: GEN
1974: eta(GEN x, long prec)
1975: {
1.2 ! noro 1976: gpmem_t av = avma;
1.1 noro 1977: GEN q = qq(x,prec);
1978: return gerepileupto(av,inteta(q));
1979: }
1980:
1981: /* returns the true value of eta(x) for Im(x) > 0, using reduction */
1982: GEN
1983: trueeta(GEN x, long prec)
1984: {
1.2 ! noro 1985: long tx=typ(x), l;
! 1986: gpmem_t av=avma, tetpil;
1.1 noro 1987: GEN p1,p2,q,q24,n,z,m,unapprox;
1988:
1989: if (!is_scalar_t(tx)) err(typeer,"trueeta");
1990: if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0)
1991: err(talker,"argument must belong to upper half-plane");
1992: l=precision(x); if (l) prec=l;
1.2 ! noro 1993: p2 = PiI2(prec);
1.1 noro 1994: z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */
1995: unapprox=gsub(gun,gpuigs(stoi(10),-8));
1996: m=gun;
1997: for(;;)
1998: {
1999: n=ground(greal(x));
2000: if (signe(n)) {x=gsub(x,n); m=gmul(m,powgi(z,n));}
2001: if (gcmp(gnorm(x), unapprox)>=0) break;
2002: m=gmul(m,gsqrt(gdiv(gi,x),prec)); x=gdivsg(-1,x);
2003: }
2004: q=gmul(p2,x);
2005: q24=gexp(gdivgs(q,24),prec); q=gpuigs(q24,24);
2006: p1=gmul(m,q24); q = inteta(q); tetpil=avma;
2007: return gerepile(av,tetpil,gmul(p1,q));
2008: }
2009:
2010: GEN
2011: eta0(GEN x, long flag,long prec)
2012: {
2013: return flag? trueeta(x,prec): eta(x,prec);
2014: }
2015:
2016: GEN
2017: jell(GEN x, long prec)
2018: {
1.2 ! noro 2019: long tx = typ(x);
! 2020: gpmem_t av=avma, tetpil;
1.1 noro 2021: GEN p1;
2022:
2023: if (!is_scalar_t(tx) || tx == t_PADIC)
2024: {
2025: GEN p2,q = qq(x,prec);
2026: p1=gdiv(inteta(gsqr(q)), inteta(q));
2027: p1=gmul2n(gsqr(p1),1);
2028: p1=gmul(q,gpuigs(p1,12));
2029: p2=gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
2030: p1=gmulsg(48,p1); tetpil=avma;
2031: return gerepile(av,tetpil,gadd(p2,p1));
2032: }
2033: p1=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
2034: p1=gsqr(gsqr(gsqr(p1)));
2035: p1=gadd(gmul2n(gsqr(p1),8), ginv(p1)); tetpil=avma;
2036: return gerepile(av,tetpil,gpuigs(p1,3));
2037: }
2038:
2039: GEN
2040: wf2(GEN x, long prec)
2041: {
1.2 ! noro 2042: gpmem_t av=avma, tetpil;
1.1 noro 2043: GEN p1,p2;
2044:
2045: p1=gsqrt(gdeux,prec);
2046: p2=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
2047: tetpil=avma;
2048: return gerepile(av,tetpil,gmul(p1,p2));
2049: }
2050:
2051: GEN
2052: wf1(GEN x, long prec)
2053: {
1.2 ! noro 2054: gpmem_t av=avma, tetpil;
1.1 noro 2055: GEN p1,p2;
2056:
2057: p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec);
2058: tetpil=avma;
2059: return gerepile(av,tetpil,gdiv(p1,p2));
2060: }
2061:
2062: GEN
2063: wf(GEN x, long prec)
2064: {
1.2 ! noro 2065: gpmem_t av=avma, tetpil;
1.1 noro 2066: GEN p1,p2;
2067:
2068: p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec));
2069: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=ldivrs(mppi(prec),-24);
2070: p2=gexp(p2,prec); tetpil=avma;
2071: return gerepile(av,tetpil,gmul(p1,p2));
2072: }
2073:
2074: GEN
2075: weber0(GEN x, long flag,long prec)
2076: {
2077: switch(flag)
2078: {
2079: case 0: return wf(x,prec);
2080: case 1: return wf1(x,prec);
2081: case 2: return wf2(x,prec);
2082: default: err(flagerr,"weber");
2083: }
2084: return NULL; /* not reached */
2085: }
2086:
2087: static GEN
2088: sagm(GEN x, long prec)
2089: {
2090: GEN p1,a,b,a1,b1,y;
2091: long l,ep;
1.2 ! noro 2092: gpmem_t av;
1.1 noro 2093:
2094: if (gcmp0(x)) return gcopy(x);
2095: switch(typ(x))
2096: {
2097: case t_REAL:
2098: l = precision(x); y = cgetr(l); av=avma;
2099: a1 = x; b1 = realun(l);
2100: l = 5-bit_accuracy(prec);
2101: do
2102: {
2103: a=a1; b=b1; a1 = addrr(a,b);
2104: setexpo(a1,expo(a1)-1);
2105: b1=mpsqrt(mulrr(a,b));
2106: }
2107: while (expo(subrr(b1,a1))-expo(b1) >= l);
2108: affrr(a1,y); avma=av; return y;
2109:
2110: case t_COMPLEX:
2111: if (gcmp0((GEN)x[2]))
2112: return transc(sagm,(GEN)x[1],prec);
2113: av=avma; l=precision(x); if (l) prec=l;
2114: a1=x; b1=gun; l = 5-bit_accuracy(prec);
2115: do
2116: {
2117: a=a1; b=b1;
2118: a1=gmul2n(gadd(a,b),-1);
2119: b1=gsqrt(gmul(a,b),prec);
2120: }
2121: while (gexpo(gsub(b1,a1))-gexpo(b1) >= l);
2122: return gerepilecopy(av,a1);
2123:
2124: case t_PADIC:
2125: av=avma; a1=x; b1=gun; l=precp(x);
2126: do
2127: {
2128: a=a1; b=b1;
2129: a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),0);
2130: p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
2131: if (ep<=0) { b1=gneg_i(b1); p1=gsub(b1,a1); ep=valp(p1)-valp(b1); }
2132: }
2133: while (ep<l && !gcmp0(p1));
2134: return gerepilecopy(av,a1);
2135:
2136: case t_SER:
2137: av=avma; a1=x; b1=gun; l=lg(x)-2;
2138: do
2139: {
2140: a=a1; b=b1;
2141: a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),prec);
2142: p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
2143: }
2144: while (ep<l && !gcmp0(p1));
2145: return gerepilecopy(av,a1);
2146:
2147: case t_INTMOD:
2148: err(impl,"agm of mod");
2149: }
2150: return transc(sagm,x,prec);
2151: }
2152:
2153: GEN
2154: agm(GEN x, GEN y, long prec)
2155: {
1.2 ! noro 2156: long ty=typ(y);
! 2157: gpmem_t av, tetpil;
1.1 noro 2158: GEN z;
2159:
2160: if (is_matvec_t(ty))
2161: {
2162: ty=typ(x);
2163: if (is_matvec_t(ty)) err(talker,"agm of two vector/matrices");
2164: z=x; x=y; y=z;
2165: }
2166: if (gcmp0(y)) return gcopy(y);
2167: av=avma; z=sagm(gdiv(x,y),prec); tetpil=avma;
2168: return gerepile(av,tetpil,gmul(y,z));
2169: }
2170:
2171: GEN
2172: logagm(GEN q)
2173: {
1.2 ! noro 2174: long prec=lg(q), s, n, lim;
! 2175: gpmem_t av=avma, tetpil;
1.1 noro 2176: GEN y,q4,q1;
2177:
2178: if (typ(q)!=t_REAL) err(typeer,"logagm");
2179: if (signe(q)<=0) err(talker,"non positive argument in logagm");
2180: if (expo(q)<0) s=1; else { q=ginv(q); s=0; }
2181: lim = - (bit_accuracy(prec) >> 1);
2182: q1 = NULL; /* gcc -Wall */
2183: for (n=0; expo(q)>=lim; n++) { q1=q; q=gsqr(q); }
2184: q4=gmul2n(q,2);
2185: if (!n) q1=gsqrt(q,prec);
2186: y=divrr(mppi(prec), agm(addsr(1,q4),gmul2n(q1,2),prec));
2187: tetpil=avma; y=gmul2n(y,-n); if (s) setsigne(y,-1);
2188: return gerepile(av,tetpil,y);
2189: }
2190:
2191: GEN
2192: glogagm(GEN x, long prec)
2193: {
1.2 ! noro 2194: gpmem_t av, tetpil;
1.1 noro 2195: GEN y,p1,p2;
2196:
2197: switch(typ(x))
2198: {
2199: case t_REAL:
2200: if (signe(x)>=0) return logagm(x);
2201:
2202: y=cgetg(3,t_COMPLEX); y[2]=lmppi(lg(x));
2203: setsigne(x,1); y[1]=(long)logagm(x);
2204: setsigne(x,-1); return y;
2205:
2206: case t_COMPLEX:
2207: y=cgetg(3,t_COMPLEX); y[2]=larg(x,prec);
2208: av=avma; p1=glogagm(gnorm(x),prec); tetpil=avma;
2209: y[1]=lpile(av,tetpil,gmul2n(p1,-1));
2210: return y;
2211:
2212: case t_PADIC:
2213: return palog(x);
2214:
2215: case t_SER:
2216: av=avma; if (valp(x)) err(negexper,"glogagm");
2217: p1=gdiv(derivser(x),x);
2218: p1=integ(p1,varn(x));
2219: if (gcmp1((GEN)x[2])) return gerepileupto(av,p1);
2220: p2=glogagm((GEN)x[2],prec); tetpil=avma;
2221: return gerepile(av,tetpil,gadd(p1,p2));
2222:
2223: case t_INTMOD:
2224: err(typeer,"glogagm");
2225: }
2226: return transc(glogagm,x,prec);
2227: }
2228:
2229: GEN
2230: theta(GEN q, GEN z, long prec)
2231: {
1.2 ! noro 2232: long l, n;
! 2233: gpmem_t av=avma, tetpil;
1.1 noro 2234: GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold;
2235:
1.2 ! noro 2236: if (!is_scalar_t(typ(q)) || !is_scalar_t(typ(z)))
! 2237: err(impl,"theta for non-scalar types");
! 2238:
1.1 noro 2239: l=precision(q); if (l) prec=l;
2240: p1=realun(prec); z=gmul(p1,z);
2241: if (!l) q=gmul(p1,q);
2242: if (gexpo(q)>=0) err(thetaer1);
2243: zy = gimag(z);
2244: zold = NULL; /* gcc -Wall */
2245: if (gcmp0(zy)) k=gzero;
2246: else
2247: {
2248: lq=glog(q,prec); k=ground(gdiv(zy,greal(lq)));
2249: if (!gcmp0(k)) { zold=z; z=gadd(z,gdiv(gmul(lq,k),gi)); }
2250: }
2251: y=gsin(z,prec); n=0; qn=gun;
2252: ps2=gsqr(q); ps=gneg_i(ps2);
2253: do
2254: {
2255: n++; p1=gsin(gmulsg(2*n+1,z),prec);
2256: qnold=qn; qn=gmul(qn,ps);
2257: ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
2258: }
2259: while (gexpo(qnold) >= -bit_accuracy(prec));
2260: if (signe(k))
2261: {
2262: y=gmul(y,gmul(gpui(q,gsqr(k),prec),
2263: gexp(gmul2n(gmul(gmul(gi,zold),k),1),prec)));
2264: if (mod2(k)) y=gneg_i(y);
2265: }
2266: p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
2267: return gerepile(av,tetpil,gmul(p1,y));
2268: }
2269:
2270: GEN
2271: thetanullk(GEN q, long k, long prec)
2272: {
1.2 ! noro 2273: long l, n;
! 2274: gpmem_t av=avma, tetpil;
1.1 noro 2275: GEN p1,ps,qn,y,ps2;
2276:
2277: l=precision(q);
2278: if (!l) { l=prec; q=gmul(q,realun(l)); }
2279: if (gexpo(q)>=0) err(thetaer1);
2280:
2281: if (!(k&1)) { avma = av; return gzero; }
2282: ps2=gsqr(q); ps=gneg_i(ps2);
2283: qn = gun; y = gun; n = 0;
2284:
2285: do
2286: {
2287: n++; p1=gpuigs(stoi(n+n+1),k); qn=gmul(qn,ps);
2288: ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
2289: }
2290: while (gexpo(p1) >= -bit_accuracy(l));
2291: p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
2292: if (k&2) { p1=gneg_i(p1); tetpil=avma; }
2293: return gerepile(av,tetpil,gmul(p1,y));
2294: }
2295:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to trans3.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath