Annotation of OpenXM_contrib/pari/src/basemath/trans3.c, Revision 1.1
1.1 ! maekawa 1: /********************************************************************/
! 2: /** **/
! 3: /** TRANSCENDENTAL FONCTIONS **/
! 4: /** (part 3) **/
! 5: /** **/
! 6: /********************************************************************/
! 7: /* $Id: trans3.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: /***********************************************************************/
! 11: /** **/
! 12: /** K BESSEL FUNCTION **/
! 13: /** **/
! 14: /***********************************************************************/
! 15:
! 16: GEN
! 17: kbessel0(GEN nu, GEN gx, long flag, long prec)
! 18: {
! 19: switch(flag)
! 20: {
! 21: case 0: return kbessel(nu,gx,prec);
! 22: case 1: return kbessel2(nu,gx,prec);
! 23: default: err(flagerr,"besselk");
! 24: }
! 25: return NULL; /* not reached */
! 26: }
! 27:
! 28: GEN
! 29: kbessel(GEN nu, GEN gx, long prec)
! 30: {
! 31: GEN x,y,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w;
! 32: long l,lbin,av,av1,k,k2,l1,n2,n;
! 33:
! 34: if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec);
! 35: if (typ(gx)!=t_REAL)
! 36: { l=prec; k=1; }
! 37: else
! 38: { l=lg(gx); k=0; x=gx; }
! 39: y=cgetr(l); l1=l+1;
! 40: av=avma; if (k) gaffect(gx,x=cgetr(l));
! 41: u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1);
! 42: e=cgetr(l1); f=cgetr(l1);
! 43: nu2=gmulgs(gsqr(nu),-4);
! 44: n = (long) (bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gnorm(nu)))) / 2;
! 45: n2=(n<<1); pitemp=mppi(l1);
! 46: /* this 10 should really be a 5, but then kbessel(15.99) enters oo loop */
! 47: lbin = 10 - bit_accuracy(l); av1=avma;
! 48: if (gcmpgs(x,n)<0)
! 49: {
! 50: zf=gsqrt(gdivgs(pitemp,n2),prec);
! 51: zz=cgetr(l1); gaffect(ginv(stoi(n2<<2)), zz);
! 52: s=gun; t=gzero;
! 53: for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
! 54: {
! 55: if (k2 & ~(MAXHALFULONG>>1))
! 56: p1 = gadd(mulss(k2,k2),nu2);
! 57: else
! 58: p1 = gaddsg(k2*k2,nu2);
! 59: ak=gdivgs(gmul(p1,zz),-k);
! 60: s=gadd(gun,gmul(ak,s));
! 61: t=gaddsg(k2,gmul(ak,t));
! 62: }
! 63: gmulz(s,zf,u); t=gmul2n(t,-1);
! 64: gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1;
! 65: affsr(n2,q=cgetr(l1));
! 66: r=gmul2n(x,1); av1=avma;
! 67: for(;;)
! 68: {
! 69: p1=divsr(5,q); if (expo(p1) >= -1) p1=ghalf;
! 70: p2=subsr(1,divrr(r,q));
! 71: if (gcmp(p1,p2)>0) p1=p2;
! 72: gnegz(p1,c); gaffsg(1,d); affrr(u,e); affrr(v,f);
! 73: for (k=1; ; k++)
! 74: {
! 75: w=gadd(gmul(gsubsg(k,ghalf),u), gmul(subrs(q,k),v));
! 76: w=gadd(w, gmul(nu,gsub(u,gmul2n(v,1))));
! 77: gdivgsz(gmul(q,v),k,u);
! 78: gdivgsz(w,k,v);
! 79: gmulz(d,c,d);
! 80: gaddz(e,gmul(d,u),e); p1=gmul(d,v);
! 81: gaddz(f,p1,f);
! 82: if (gexpo(p1)-gexpo(f) <= lbin) break;
! 83: avma=av1;
! 84: }
! 85: p1=u; u=e; e=p1;
! 86: p1=v; v=f; f=p1;
! 87: gmulz(q,gaddsg(1,c),q);
! 88: if (expo(subrr(q,r)) <= lbin) break;
! 89: }
! 90: gmulz(u,gpui(gdivgs(x,n),nu,prec),y);
! 91: }
! 92: else
! 93: {
! 94: p2=gmul2n(x,1);
! 95: zf=gsqrt(gdiv(pitemp,p2),prec);
! 96: zz=ginv(gmul2n(p2,2)); s=gun;
! 97: for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
! 98: {
! 99: if (k2 & ~(MAXHALFULONG>>1))
! 100: p1=gadd(mulss(k2,k2),nu2);
! 101: else
! 102: p1=gaddsg(k2*k2,nu2);
! 103: ak=gdivgs(gmul(p1,zz),k);
! 104: s=gsub(gun,gmul(ak,s));
! 105: }
! 106: gmulz(s,zf,y);
! 107: }
! 108: gdivz(y,mpexp(x),y); avma=av; return y;
! 109: }
! 110:
! 111: /***********************************************************************/
! 112: /* **/
! 113: /** FONCTION U(a,b,z) GENERALE **/
! 114: /** ET CAS PARTICULIERS **/
! 115: /** **/
! 116: /***********************************************************************/
! 117:
! 118: /* Assume gx > 0 and a,b real */
! 119: GEN
! 120: hyperu(GEN a, GEN b, GEN gx, long prec)
! 121: {
! 122: GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn;
! 123: long l,lbin,av,av1,av2,k,l1,n,ex;
! 124:
! 125: if (typ(gx)!=t_REAL)
! 126: { l=prec; k=1; }
! 127: else
! 128: { l=lg(gx); k=0; x=gx; }
! 129: ex = (iscomplex(a) || iscomplex(b));
! 130: if (ex) { y=cgetg(3,t_COMPLEX); y[1]=lgetr(l); y[2]=lgetr(l); }
! 131: else y=cgetr(l);
! 132: l1=l+1; av=avma; if (k) gaffect(gx,x=cgetr(l));
! 133: a1=gaddsg(1,gsub(a,b));
! 134: if (ex)
! 135: {
! 136: u=cgetg(3,t_COMPLEX); u[1]=lgetr(l1); u[2]=lgetr(l1);
! 137: v=cgetg(3,t_COMPLEX); v[1]=lgetr(l1); v[2]=lgetr(l1);
! 138: c=cgetg(3,t_COMPLEX); c[1]=lgetr(l1); c[2]=lgetr(l1);
! 139: d=cgetg(3,t_COMPLEX); d[1]=lgetr(l1); d[2]=lgetr(l1);
! 140: e=cgetg(3,t_COMPLEX); e[1]=lgetr(l1); e[2]=lgetr(l1);
! 141: f=cgetg(3,t_COMPLEX); f[1]=lgetr(l1); f[2]=lgetr(l1);
! 142: }
! 143: else
! 144: {
! 145: u=cgetr(l1); v=cgetr(l1); c=cgetr(l1);
! 146: d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
! 147: }
! 148: n=(long)(bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
! 149: lbin = 10-bit_accuracy(l); av1=avma;
! 150: if (gcmpgs(x,n)<0)
! 151: {
! 152: gn=stoi(n); zf=gpui(gn,gneg_i(a),l1);
! 153: zz=gdivsg(-1,gn); s=gun; t=gzero;
! 154: for (k=n-1; k>=0; k--)
! 155: {
! 156: p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
! 157: s=gaddsg(1,gmul(p1,s));
! 158: t=gadd(gaddsg(k,a),gmul(p1,t));
! 159: }
! 160: gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1;
! 161: q=cgetr(l1); affsr(n,q); r=x; av1=avma;
! 162: do
! 163: {
! 164: p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf;
! 165: p2=subsr(1,divrr(r,q)); if (gcmp(p1,p2)>0) p1=p2;
! 166: gnegz(p1,c); avma=av1;
! 167: k=0; gaffsg(1,d);
! 168: gaffect(u,e); gaffect(v,f);
! 169: p3=gsub(q,b); av2=avma;
! 170: for(;;)
! 171: {
! 172: w=gadd(gmul(gaddsg(k,a),u),gmul(gaddsg(-k,p3),v));
! 173: k++; gdivgsz(gmul(q,v),k,u);
! 174: gdivgsz(w,k,v);
! 175: gmulz(d,c,d);
! 176: gaddz(e,gmul(d,u),e); p1=gmul(d,v);
! 177: gaddz(f,p1,f);
! 178: if (gexpo(p1)-gexpo(f) <= lbin) break;
! 179: avma=av2;
! 180: }
! 181: p1=u; u=e; e=p1;
! 182: p1=v; v=f; f=p1;
! 183: gmulz(q,gadd(gun,c),q);
! 184: p1=subrr(q,r); ex=expo(p1); avma=av1;
! 185: }
! 186: while (ex>lbin);
! 187: }
! 188: else
! 189: {
! 190: zf=gpui(x,gneg_i(a),l1);
! 191: zz=gdivsg(-1,x); s=gun;
! 192: for (k=n-1; k>=0; k--)
! 193: {
! 194: p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
! 195: s=gadd(gun,gmul(p1,s));
! 196: }
! 197: u = gmul(s,zf);
! 198: }
! 199: gaffect(u,y); avma=av; return y;
! 200: }
! 201:
! 202: GEN
! 203: kbessel2(GEN nu, GEN x, long prec)
! 204: {
! 205: long av = avma,tetpil;
! 206: GEN p1,p2,x2,a,pitemp;
! 207:
! 208: if (typ(x)==t_REAL) prec = lg(x);
! 209: x2=gshift(x,1); pitemp=mppi(prec);
! 210: if (gcmp0(gimag(nu))) a=cgetr(prec);
! 211: else
! 212: { a=cgetg(3,t_COMPLEX); a[1]=lgetr(prec); a[2]=lgetr(prec); }
! 213: gaddz(gun,gshift(nu,1),a);
! 214: p1=hyperu(gshift(a,-1),a,x2,prec);
! 215: p1=gmul(gmul(p1,gpui(x2,nu,prec)),mpsqrt(pitemp));
! 216: p2=gexp(x,prec); tetpil=avma;
! 217: return gerepile(av,tetpil,gdiv(p1,p2));
! 218: }
! 219:
! 220: GEN
! 221: incgam(GEN a, GEN x, long prec)
! 222: {
! 223: GEN p1,z = cgetr(prec);
! 224: long av = avma;
! 225:
! 226: if (typ(x)!=t_REAL) { gaffect(x,z); x=z; }
! 227: if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0)
! 228: p1 = incgam2(a,x,prec);
! 229: else
! 230: p1 = gsub(ggamma(a,prec), incgam3(a,x,prec));
! 231: affrr(p1,z); avma = av; return z;
! 232: }
! 233:
! 234: /* = incgam2(0, x, prec). typ(x) = t_REAL. Optimized for eint1 */
! 235: static GEN
! 236: incgam2_0(GEN x)
! 237: {
! 238: long l,n,i;
! 239: GEN p1;
! 240: double m,mx;
! 241:
! 242: l = lg(x); mx = rtodbl(x);
! 243: m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
! 244: p1 = divsr(-n, addsr(n<<1,x));
! 245: for (i=n-1; i >= 1; i--)
! 246: p1 = divsr(-i, addrr(addsr(i<<1,x), mulsr(i,p1)));
! 247: return mulrr(divrr(mpexp(negr(x)), x), addrr(realun(l),p1));
! 248: }
! 249:
! 250: GEN
! 251: incgam1(GEN a, GEN x, long prec)
! 252: {
! 253: GEN p2,p3,y, z = cgetr(prec);
! 254: long av=avma, av1, l,n,i;
! 255: double m,mx;
! 256:
! 257: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
! 258: l=lg(x); mx=rtodbl(x);
! 259: m=(long) bit_accuracy(l)*LOG2; n=(long)(m/(log(m)-(1+log(mx))));
! 260: p2 = cgetr(l); affrr(addir(gun,gsub(x,a)), p2);
! 261: p3 = subrs(p2, n+1); av1 = avma;
! 262: for (i=n; i>=1; i--)
! 263: {
! 264: affrr(addrr(subrs(p2,i), divrr(mulsr(i,x),p3)), p3);
! 265: avma = av1;
! 266: }
! 267: y = mulrr(mpexp(negr(x)),gpui(x,a,prec));
! 268: affrr(divrr(y,p3), z);
! 269: avma = av; return z;
! 270: }
! 271:
! 272: GEN
! 273: incgam2(GEN a, GEN x, long prec)
! 274: {
! 275: GEN b,p1,p2,p3,y, z = cgetr(prec);
! 276: long av = avma,av1,l,n,i;
! 277: double m,mx;
! 278:
! 279: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
! 280: if (gcmp0(a)) { affrr(incgam2_0(x), z); avma = av; return z; }
! 281: l=lg(x); mx=rtodbl(x);
! 282: m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
! 283: i = typ(a);
! 284: if (i == t_REAL) b = addsr(-1,a);
! 285: else
! 286: { /* keep b integral : final powering more efficient */
! 287: gaffect(a,p1=cgetr(prec));
! 288: b = (i == t_INT)? addsi(-1,a): addsr(-1,p1);
! 289: a = p1;
! 290: }
! 291: p2 = cgetr(l); gaffect(subrr(x,a),p2);
! 292: p3 = divrr(addsr(-n,a), addrs(p2,n<<1));
! 293: av1 = avma;
! 294: for (i=n-1; i>=1; i--)
! 295: {
! 296: affrr(divrr(addsr(-i,a), addrr(addrs(p2,i<<1),mulsr(i,p3))), p3);
! 297: avma = av1;
! 298: }
! 299: y = gmul(mpexp(negr(x)), gpui(x,b,prec));
! 300: affrr(mulrr(y, addsr(1,p3)), z);
! 301: avma = av; return z;
! 302: }
! 303:
! 304: GEN
! 305: incgam3(GEN a, GEN x, long prec)
! 306: {
! 307: GEN b,p1,p2,p3,y, z = cgetr(prec);
! 308: long av = avma, av1,l,n,i;
! 309:
! 310: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
! 311: l=lg(x); n = -bit_accuracy(l)-1;
! 312: p3 = realun(l);
! 313: p2 = rcopy(p3);
! 314: i = typ(a);
! 315: if (i == t_REAL) b = a;
! 316: else
! 317: {
! 318: gaffect(a,p1=cgetr(prec));
! 319: b = (i == t_INT)? a: p1;
! 320: a = p1;
! 321: }
! 322: av1 = avma;
! 323: for (i=1; expo(p2) >= n; i++)
! 324: {
! 325: affrr(divrr(mulrr(x,p2), addsr(i,a)), p2);
! 326: affrr(addrr(p2,p3), p3); avma = av1;
! 327: }
! 328: y = gdiv(mulrr(mpexp(negr(x)), gpui(x,b,prec)), a);
! 329: affrr(mulrr(y,p3), z);
! 330: avma = av; return z;
! 331: }
! 332:
! 333: /* Assumes that g=gamma(a,prec). Unchecked */
! 334: GEN
! 335: incgam4(GEN a, GEN x, GEN g, long prec)
! 336: {
! 337: GEN p1, z = cgetr(prec);
! 338: long av = avma;
! 339:
! 340: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
! 341: if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0)
! 342: p1 = incgam2(a,x,prec);
! 343: else
! 344: p1 = gsub(g, incgam3(a,x,prec));
! 345: affrr(p1, z);
! 346: avma = av; return z;
! 347: }
! 348:
! 349: GEN
! 350: incgam0(GEN a, GEN x, GEN z,long prec)
! 351: {
! 352: return z? incgam4(a,x,z,prec): incgam(a,x,prec);
! 353: }
! 354:
! 355: GEN
! 356: eint1(GEN x, long prec)
! 357: {
! 358: long av = avma,tetpil,l,n,i;
! 359: GEN p1,p2,p3,p4,run,y;
! 360:
! 361: if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;}
! 362: if (signe(x) >= 0)
! 363: {
! 364: if (expo(x) >= 4)
! 365: return gerepileupto(av, incgam2_0(x));
! 366:
! 367: l = lg(x); consteuler(l);
! 368: n = -bit_accuracy(l)-1;
! 369:
! 370: run = realun(l);
! 371: p4 = p3 = p2 = p1 = run;
! 372: for (i=2; expo(p2)>=n; i++)
! 373: {
! 374: p1 = addrr(p1,divrs(run,i)); /* p1 = sum_{i=1} 1/i */
! 375: p4 = divrs(mulrr(x,p4),i); /* p4 = sum_{i=1} x^(i-1)/i */
! 376: p2 = mulrr(p4,p1);
! 377: p3 = addrr(p2,p3);
! 378: }
! 379: p3 = mulrr(x,mulrr(mpexp(negr(x)),p3));
! 380: p1 = addrr(mplog(x),geuler);
! 381: return gerepileupto(av, subrr(p3,p1));
! 382: }
! 383: else
! 384: { /* written by Manfred Radimersky */
! 385: l = lg(x);
! 386: n = bit_accuracy(l);
! 387: /* IS: line split to avoid a Workshop cc-5.0 bug (Sun BugID #4254959) */
! 388: n = 3 * n / 4;
! 389: y = negr(x);
! 390: if(gcmpgs(y, n) < 0) {
! 391: p1 = p2 = p3 = y;
! 392: p4 = gzero;
! 393: i = 2;
! 394: while(gcmp(p3, p4)) {
! 395: p4 = p3;
! 396: p1 = gmul(p1, gdivgs(y, i));
! 397: p2 = gdivgs(p1, i);
! 398: p3 = gadd(p3, p2);
! 399: i++;
! 400: }
! 401: consteuler(l);
! 402: p1 = gadd(mplog(y), geuler);
! 403: y = gadd(p3, p1);
! 404: } else {
! 405: p1 = gdivsg(1, y);
! 406: p2 = realun(l);
! 407: p3 = p2;
! 408: p4 = gzero;
! 409: i = 1;
! 410: while(gcmp(p3, p4)) {
! 411: p4 = p3;
! 412: p2 = gmulgs(p2, i);
! 413: p2 = gmul(p2, p1);
! 414: p3 = gadd(p3, p2);
! 415: i++;
! 416: }
! 417: p1 = gdiv(mpexp(y), y);
! 418: y = gmul(p3, p1);
! 419: }
! 420: tetpil = avma;
! 421: y = gerepile(av, tetpil, negr(y));
! 422: }
! 423: return y;
! 424: }
! 425:
! 426: GEN
! 427: veceint1(GEN C, GEN nmax, long prec)
! 428: {
! 429: long av,av1,k,n,nstop,i,cd,nmin,G,a;
! 430: GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit;
! 431:
! 432: if (!nmax) return eint1(C,prec);
! 433:
! 434: if (signe(nmax)<=0) return cgetg(1,t_VEC);
! 435: if (DEBUGLEVEL>1) fprintferr("Entering veceint1:\n");
! 436: if (typ(C) != t_REAL || lg(C) > prec)
! 437: { y = cgetr(prec); gaffect(C,y); C = y; }
! 438: if (signe(C) <= 0) err(talker,"negative or zero constant in veceint1");
! 439:
! 440: G = -bit_accuracy(prec);
! 441: n=itos(nmax); y=cgetg(n+1,t_VEC);
! 442: for(i=1; i<=n; i++) y[i]=lgetr(prec);
! 443: av=avma;
! 444:
! 445: nstop = itos(gceil(divsr(4,C)));
! 446: if (nstop<1) nstop=1;
! 447: if (nstop>n) nstop=n;
! 448: nmin=n-10; if (nmin<nstop) nmin=nstop;
! 449: if(DEBUGLEVEL>1) fprintferr("nstop = %ld\n",nstop);
! 450:
! 451: e1=mpexp(negr(mulsr(n,C)));
! 452: e2=mpexp(mulsr(10,C));
! 453: unr = realun(prec);
! 454: zeror=realzero(prec); deninit=negr(unr);
! 455: f=cgetg(3,t_COL); M2=cgetg(3,t_VEC); av1=avma;
! 456:
! 457: F0=(GEN)y[n];
! 458: affrr(eint1(mulsr(n,C),prec), F0);
! 459: do
! 460: {
! 461: if (DEBUGLEVEL>1) fprintferr("%ld ",n);
! 462: minvn=divrs(unr,-n); mcn=mulrr(C,minvn);
! 463: M2[1] = (long)zeror; M2[2] = lsubrr(minvn,C);
! 464: f[1]=(long)zeror; f[2]=ldivrs(e1,-n);
! 465: affrr(mulrr(e1,e2), e1);
! 466: vdiff=cgetg(2,t_VEC); vdiff[1]=f[2];
! 467: for (cd=a=1,n--; n>=nmin; n--,a++)
! 468: {
! 469: F = F0;
! 470: ap = stoi(a); den = deninit;
! 471: for (k=1;;)
! 472: {
! 473: if (k>cd)
! 474: {
! 475: cd++; p1 = (GEN)f[2];
! 476: f[2] = lmul(M2,f);
! 477: f[1] = (long)p1;
! 478: M2[1] = laddrr((GEN)M2[1],mcn);
! 479: M2[2] = laddrr((GEN)M2[2],minvn);
! 480: vdiff = concatsp(vdiff,(GEN)f[2]);
! 481: }
! 482: p1 = mulrr(mulri(den,ap),(GEN)vdiff[k]);
! 483: if (expo(p1) < G) { affrr(F,(GEN)y[n]); break; }
! 484: F = addrr(F,p1); ap = mulis(ap,a);
! 485: k++; den = divrs(den,-k);
! 486: }
! 487: }
! 488: avma=av1; n++; F0=(GEN)y[n]; nmin -= 10;
! 489: if (nmin < nstop) nmin=nstop;
! 490: }
! 491: while(n>nstop);
! 492: for(i=1; i<=nstop; i++)
! 493: affrr(eint1(mulsr(i,C),prec), (GEN)y[i]);
! 494: if (DEBUGLEVEL>1) fprintferr("\n");
! 495: avma=av; return y;
! 496: }
! 497:
! 498: GEN
! 499: gerfc(GEN x, long prec)
! 500: {
! 501: long av=avma;
! 502: GEN p1,p2;
! 503:
! 504: if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; }
! 505: av = avma; p1 = incgam(ghalf,gsqr(x),prec);
! 506: p2 = mpsqrt(mppi(lg(x)));
! 507: p1 = divrr(p1,p2);
! 508: if (signe(x) < 0) p1 = subsr(2,p1);
! 509: return gerepileupto(av,p1);
! 510: }
! 511:
! 512: /***********************************************************************/
! 513: /** **/
! 514: /** FONCTION ZETA DE RIEMANN **/
! 515: /** **/
! 516: /***********************************************************************/
! 517:
! 518: static GEN
! 519: czeta(GEN s, long prec)
! 520: {
! 521: long av,n,p,n1,l,flag1,flag2,flag3,i,i2;
! 522: double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf;
! 523: GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp;
! 524:
! 525: l=precision(s);
! 526: if (typ(s)==t_COMPLEX)
! 527: {
! 528: if (!l) l=prec;
! 529: res=cgetg(3,t_COMPLEX); res[1]=lgetr(l); res[2]=lgetr(l);
! 530: av=avma;
! 531: p1=cgetg(3,t_COMPLEX); p1[1]=lgetr(l+1); p1[2]=lgetr(l+1);
! 532: gaffect(s,p1); s=p1; sig=(GEN)s[1];
! 533: }
! 534: else
! 535: {
! 536: res = cgetr(l); av=avma;
! 537: p1=cgetr(l+1); affrr(s,p1); sig=s=p1;
! 538: }
! 539:
! 540: if (signe(sig)>0 && expo(sig)>-2) flag1 = 0;
! 541: else
! 542: {
! 543: if (gcmp0(gimag(s)) && gcmp0(gfrac(gmul2n(sig,-1))))
! 544: {
! 545: if (gcmp0(sig)) gaffect(gneg_i(ghalf),res); else gaffsg(0,res);
! 546: avma=av; return res;
! 547: }
! 548: flag1=1; s=gsub(gun,s); sig=greal(s);
! 549: }
! 550: ssig=rtodbl(sig); st=fabs(rtodbl(gimag(s))); maxbeta = pow(3.0,-2.5);
! 551: if (st)
! 552: {
! 553: ns = ssig*ssig + st*st;
! 554: alpha=pariC2*(prec-2)-0.39-0.5*(ssig-1.0)*log(ns)-log(ssig)+ssig*2*pariC1;
! 555: beta=(alpha+ssig)/st-atan(ssig/st);
! 556: if (beta<=0)
! 557: {
! 558: if (ssig>=1.0)
! 559: {
! 560: p=0; sn=sqrt(ns);
! 561: n=(long)(ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig)));
! 562: }
! 563: else
! 564: {
! 565: p=1; sn=ssig+1; n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
! 566: }
! 567: }
! 568: else
! 569: {
! 570: if (beta<maxbeta) xinf=beta+pow(3*beta,1.0/3.0);
! 571: else
! 572: {
! 573: double eps=0.0087, x00 = beta+PI/2.0, y00,x11;
! 574: for(;;)
! 575: {
! 576: y00=x00*x00; x11=(beta+atan(x00))*(1+y00)/y00-1/x00;
! 577: if (x00-x11 < eps) break;
! 578: x00 = x11;
! 579: }
! 580: xinf=x11;
! 581: }
! 582: sp=1.0-ssig+st*xinf;
! 583: if (sp>0)
! 584: {
! 585: p=(long)ceil(sp/2.0); sn=ssig+2*p-1;
! 586: n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
! 587: }
! 588: else
! 589: {
! 590: p=0; sn=sqrt(ns);
! 591: n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
! 592: }
! 593: }
! 594: }
! 595: else
! 596: {
! 597: beta=pariC2*(prec-2)+0.61+ssig*2*pariC1-ssig*log(ssig);
! 598: if (beta>0)
! 599: {
! 600: p=(long)ceil(beta/2.0); sn=ssig+2*p-1;
! 601: n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
! 602: }
! 603: else
! 604: {
! 605: p=0; sn=sqrt(ssig*ssig+st*st);
! 606: n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
! 607: }
! 608: }
! 609: if (n < 46340) { flag2=1; n1=n*n; } else flag2=0;
! 610: y=gun; ms=gneg_i(s); p1=cgetr(prec+1);
! 611: for (i=2; i<=n; i++)
! 612: {
! 613: affsr(i,p1); p2 = gexp(gmul(ms,mplog(p1)), prec+1);
! 614: y = gadd(y,p2);
! 615: }
! 616: flag3 = (2*p < 46340);
! 617: mpbern(p,prec+1); p31=cgetr(prec+1); z=gzero;
! 618: for (i=p; i>=1; i--)
! 619: {
! 620: i2=i<<1;
! 621: p1=gmul(gaddsg(i2-1,s),gaddsg(i2,s));
! 622: p1=flag3? gdivgs(p1,i2*(i2+1)): gdivgs(gdivgs(p1,i2),i2+1);
! 623: p1=flag2? gdivgs(p1,n1): gdivgs(gdivgs(p1,n),n);
! 624: p3 = bern(i);
! 625: if (bernzone[2]>prec+1) { affrr(p3,p31); p3=p31; }
! 626: z=gadd(divrs(p3,i),gmul(p1,z));
! 627: }
! 628: p1=gsub(gdivsg(n,gsubgs(s,1)),ghalf);
! 629: z=gmul(gadd(p1,gmul(s,gdivgs(z,n<<1))),p2);
! 630: y = gadd(y,z);
! 631: if (flag1)
! 632: {
! 633: pitemp=mppi(prec+1); setexpo(pitemp,2);
! 634: y=gmul(gmul(y,ggamma(s,prec+1)),gpui(pitemp,ms,prec+1));
! 635: setexpo(pitemp,0);
! 636: y=gmul2n(gmul(y,gcos(gmul(pitemp,s),prec+1)),1);
! 637: }
! 638: gaffect(y,res); avma=av; return res;
! 639: }
! 640:
! 641: /* y = binomial(n,k-2). Return binomial(n,k) */
! 642: static GEN
! 643: next_bin(GEN y, long n, long k)
! 644: {
! 645: y = divrs(mulrs(y, n-k+2), k-1);
! 646: return divrs(mulrs(y, n-k+1), k);
! 647: }
! 648:
! 649: static GEN
! 650: izeta(long k, long prec)
! 651: {
! 652: long av=avma,av2,tetpil,kk,n,li, limit;
! 653: GEN y,p1,pitemp,qn,z,q,binom;
! 654:
! 655: if (!k) return gneg(ghalf);
! 656: if (k<0)
! 657: {
! 658: if ((k&1) == 0) return gzero;
! 659: y=bernreal(1-k,prec); tetpil=avma;
! 660: return gerepile(av,tetpil,divrs(y,k-1));
! 661: }
! 662: if (k > bit_accuracy(prec)+1) return realun(prec);
! 663: pitemp=mppi(prec); setexpo(pitemp,2);
! 664: if ((k&1) == 0)
! 665: {
! 666: p1 = mulrr(gpuigs(pitemp,k),absr(bernreal(k,prec)));
! 667: y = mpfactr(k,prec); tetpil=avma;
! 668: y=divrr(p1,y); setexpo(y,expo(y)-1);
! 669: return gerepile(av,tetpil,y);
! 670: }
! 671: binom = realun(prec+1);
! 672: q = mpexp(pitemp); kk = k+1;
! 673: li = -(1+bit_accuracy(prec));
! 674: if ((k&3)==3)
! 675: {
! 676: for (n=0; n <= kk>>1; n+=2)
! 677: {
! 678: p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
! 679: if (n) { binom = next_bin(binom,kk,n); setlg(binom,prec+1); }
! 680: p1 = mulrr(binom,p1);
! 681: if (n == kk>>1) setexpo(p1, expo(p1)-1);
! 682: if ((n>>1)&1) setsigne(p1,-signe(p1));
! 683: y = n? addrr(y,p1): p1;
! 684: }
! 685: y=mulrr(divrr(gpuigs(pitemp,k),mpfactr(kk,prec)),y);
! 686:
! 687: av2 = avma; limit = stack_lim(av2,1);
! 688: qn=gsqr(q); z=divsr(1,addrs(q,-1));
! 689: for (n=2; ; n++)
! 690: {
! 691: p1=divsr(1,mulir(gpuigs(stoi(n),k),addrs(qn,-1)));
! 692:
! 693: z=addrr(z,p1); if (expo(p1)< li) break;
! 694: qn=mulrr(qn,q);
! 695: if (low_stack(limit,stack_lim(av2,1)))
! 696: {
! 697: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;
! 698: if (DEBUGMEM>1) err(warnmem,"izeta");
! 699: gerepilemany(av2,gptr,2);
! 700: }
! 701: }
! 702: setexpo(z,expo(z)+1); tetpil = avma;
! 703: y = addrr(y,z); setsigne(y,-signe(y));
! 704: }
! 705: else
! 706: {
! 707: GEN p2 = divrs(pitemp, k-1);
! 708: for (n=0; n <= k>>1; n+=2)
! 709: {
! 710: p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
! 711: if (n) binom = next_bin(binom,kk,n);
! 712: p1 = mulrr(binom,p1);
! 713: p1 = mulsr(kk-(n<<1),p1);
! 714: if ((n>>1)&1) setsigne(p1,-signe(p1));
! 715: y = n? addrr(y,p1): p1;
! 716: }
! 717: y = mulrr(divrr(gpuigs(pitemp,k),mpfactr(kk,prec)),y);
! 718: y = divrs(y,k-1);
! 719: av2 = avma; limit = stack_lim(av2,1);
! 720: qn = q; z=gzero;
! 721: for (n=1; ; n++)
! 722: {
! 723: p1=mulir(gpuigs(stoi(n),k),gsqr(addrs(qn,-1)));
! 724: p1=divrr(addrs(mulrr(qn,addsr(1,mulsr(n<<1,p2))),-1),p1);
! 725:
! 726: z=addrr(z,p1); if (expo(p1) < li) break;
! 727: qn=mulrr(qn,q);
! 728: if (low_stack(limit,stack_lim(av2,1)))
! 729: {
! 730: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;
! 731: if (DEBUGMEM>1) err(warnmem,"izeta");
! 732: gerepilemany(av2,gptr,2);
! 733: }
! 734: }
! 735: setexpo(z,expo(z)+1); tetpil = avma;
! 736: y = subrr(y,z);
! 737: }
! 738: return gerepile(av,tetpil,y);
! 739: }
! 740:
! 741: GEN
! 742: gzeta(GEN x, long prec)
! 743: {
! 744: if (gcmp1(x)) err(talker,"argument equal to one in zeta");
! 745: switch(typ(x))
! 746: {
! 747: case t_INT:
! 748: return izeta(itos(x),prec);
! 749:
! 750: case t_REAL: case t_COMPLEX:
! 751: return czeta(x,prec);
! 752:
! 753: case t_INTMOD: case t_PADIC:
! 754: err(typeer,"gzeta");
! 755: case t_SER:
! 756: err(impl,"zeta of power series");
! 757: }
! 758: return transc(gzeta,x,prec);
! 759: }
! 760:
! 761: void
! 762: gzetaz(GEN x, GEN y)
! 763: {
! 764: long av=avma, prec = precision(y);
! 765:
! 766: if (!prec) err(infprecer,"gzetaz");
! 767: gaffect(gzeta(x,prec),y); avma=av;
! 768: }
! 769:
! 770: /***********************************************************************/
! 771: /** **/
! 772: /** FONCTIONS POLYLOGARITHME **/
! 773: /** **/
! 774: /***********************************************************************/
! 775:
! 776: /* validity domain contains .005 < |x| < 230 */
! 777: static GEN
! 778: cxpolylog(long m, GEN x, long prec)
! 779: {
! 780: long av=avma,li,i,n,bern_upto;
! 781: GEN p1,z,h,q,s;
! 782:
! 783: if (gcmp1(x)) return izeta(m,prec);
! 784: z=glog(x,prec); h=gneg_i(glog(gneg_i(z),prec));
! 785: for (i=1; i<m; i++) h = gadd(h,ginv(stoi(i)));
! 786:
! 787: bern_upto=m+50; mpbern(bern_upto,prec);
! 788: q=gun; s=izeta(m,prec);
! 789: for (n=1; n<=m+1; n++)
! 790: {
! 791: q=gdivgs(gmul(q,z),n);
! 792: s=gadd(s,gmul((n==(m-1))? h: izeta(m-n,prec),q));
! 793: }
! 794:
! 795: n = m+3; z=gsqr(z); li = -(bit_accuracy(prec)+1);
! 796:
! 797: for(;;)
! 798: {
! 799: q = gdivgs(gmul(q,z),(n-1)*n);
! 800: p1 = gmul(izeta(m-n,prec),q);
! 801: s = gadd(s,p1);
! 802: if (gexpo(p1) < li) break;
! 803: n+=2;
! 804: if (n>=bern_upto) { bern_upto += 50; mpbern(bern_upto,prec); }
! 805: }
! 806: return gerepileupto(av,s);
! 807: }
! 808:
! 809: GEN
! 810: polylog(long m, GEN x, long prec)
! 811: {
! 812: long av,av1,limpile,tetpil,l,e,sx,i;
! 813: GEN p1,p2,n,y,logx,logx2;
! 814: GEN *gptr[4];
! 815:
! 816: if (m<0) err(talker,"negative index in polylog");
! 817: if (!m) return gneg(ghalf);
! 818: if (gcmp0(x)) return gcopy(x);
! 819: av=avma;
! 820: if (m==1)
! 821: {
! 822: p1=glog(gsub(gun,x),prec); tetpil=avma;
! 823: return gerepile(av,tetpil,gneg(p1));
! 824: }
! 825: l=precision(x);
! 826: if (!l) { l=prec; x=gmul(x, realun(l)); }
! 827: e=gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
! 828: if (e>0)
! 829: {
! 830: logx=glog(x,l); sx=gsigne(gimag(x));
! 831: if (!sx)
! 832: {
! 833: if (m&1)
! 834: sx = gsigne(greal(gsub(gun,x)));
! 835: else
! 836: sx = -gsigne(greal(x));
! 837: }
! 838: x=ginv(x);
! 839: }
! 840: av1=avma; limpile=stack_lim(av1,1);
! 841: y=x; n=gun; p1=x;
! 842: do
! 843: {
! 844: n=addsi(1,n); p1=gmul(x,p1); p2=gdiv(p1,gpuigs(n,m));
! 845: y=gadd(y,p2);
! 846: if (low_stack(limpile, stack_lim(av1,1)))
! 847: {
! 848: if(DEBUGMEM>1) err(warnmem,"polylog");
! 849: gptr[0]=&y; gptr[1]=&n; gptr[2]=&p1; gptr[3]=&p2;
! 850: gerepilemany(av1,gptr,4);
! 851: }
! 852: }
! 853: while (gexpo(p2) > -bit_accuracy(l));
! 854: tetpil=avma;
! 855: if (e<=0) return gerepile(av,tetpil,gcopy(y));
! 856: logx2=gsqr(logx); p1=gneg_i(ghalf);
! 857: for (i=m-2; i>=0; i-=2)
! 858: {
! 859: p1 = gadd(izeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
! 860: }
! 861: if (m&1) p1 = gmul(logx,p1);
! 862: p2=cgetg(3,t_COMPLEX); p2[1]=zero;
! 863: p2[2]=ldiv(gmulsg(sx,mppi(l)),mpfact(m-1));
! 864: p1=gadd(gmul2n(p1,1), gmul(p2,gpuigs(logx,m-1)));
! 865: if ((m&1) == 0) y = gneg_i(y);
! 866: tetpil=avma; return gerepile(av,tetpil,gadd(p1,y));
! 867: }
! 868:
! 869: GEN
! 870: dilog(GEN x, long prec)
! 871: {
! 872: GEN p1,p2,p3,y;
! 873: long av=avma,tetpil;
! 874:
! 875: if (gcmpgs(gnorm(x),1)<1)
! 876: {
! 877: tetpil=avma; return gerepile(av,tetpil,polylog(2,x,prec));
! 878: }
! 879: y=gneg_i(polylog(2,ginv(x),prec)); p3=mppi(prec); p2=gdivgs(gsqr(p3),6);
! 880: p1=glog(gneg_i(x),prec); p1=gadd(p2,gmul2n(gsqr(p1),-1));
! 881: p1 = gneg_i(p1); tetpil=avma;
! 882: return gerepile(av,tetpil,gadd(y,p1));
! 883: }
! 884:
! 885: GEN
! 886: polylogd0(long m, GEN x, long flag, long prec)
! 887: {
! 888: long k,l,av,tetpil,fl,m2;
! 889: GEN p1,p2,p3,y;
! 890:
! 891: m2=m&1; av=avma;
! 892: if (gcmp0(x)) return gcopy(x);
! 893: if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
! 894: l=precision(x);
! 895: if (!l) { l=prec; x=gmul(x,realun(l)); }
! 896: p1=gabs(x,prec); fl=0;
! 897: if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
! 898:
! 899: p1=gneg_i(glog(p1,prec)); p2=gun;
! 900: y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
! 901: for (k=1; k<m; k++)
! 902: {
! 903: p2=gdivgs(gmul(p2,p1),k);
! 904: p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
! 905: tetpil=avma; y=gadd(y,gmul(p2,p3));
! 906: }
! 907: if (m2)
! 908: {
! 909: if (flag)
! 910: p2 = gdivgs(gmul(p2,p1),-2*m);
! 911: else
! 912: p2 = gdivgs(gmul(glog(gnorm(gsub(gun,x)),prec),p2),2*m);
! 913: tetpil=avma; y=gadd(y,p2);
! 914: }
! 915: if (fl) { tetpil=avma; return gerepile(av,tetpil,gneg(y)); }
! 916: return gerepile(av,tetpil,y);
! 917: }
! 918:
! 919: GEN
! 920: polylogd(long m, GEN x, long prec)
! 921: {
! 922: return polylogd0(m,x,0,prec);
! 923: }
! 924:
! 925: GEN
! 926: polylogdold(long m, GEN x, long prec)
! 927: {
! 928: return polylogd0(m,x,1,prec);
! 929: }
! 930:
! 931: GEN
! 932: polylogp(long m, GEN x, long prec)
! 933: {
! 934: long k,l,av,tetpil,fl,m2;
! 935: GEN p1,y;
! 936:
! 937: m2=m&1; av=avma;
! 938: if (gcmp0(x)) return gcopy(x);
! 939: if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
! 940: l=precision(x);
! 941: if (!l) { l=prec; x=gmul(x,realun(l)); }
! 942: p1=gabs(x,prec); fl=0;
! 943: if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
! 944:
! 945: p1=gmul2n(glog(p1,prec),1); mpbern(m>>1,prec);
! 946: y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
! 947:
! 948: if (m==1)
! 949: {
! 950: p1=gmul2n(p1,-2); tetpil=avma; y=gadd(y,p1);
! 951: }
! 952: else
! 953: {
! 954: GEN p2=gun, p3, p4, p5, p51=cgetr(prec);
! 955:
! 956: for (k=1; k<m; k++)
! 957: {
! 958: p2=gdivgs(gmul(p2,p1),k);
! 959: if (!(k&1) || k==1)
! 960: {
! 961: if (k!=1)
! 962: {
! 963: p5=(GEN)bern(k>>1);
! 964: if (bernzone[2]>prec) { affrr(p5,p51); p5=p51; }
! 965: p4=gmul(p2,p5);
! 966: }
! 967: else p4=gneg_i(gmul2n(p2,-1));
! 968: p3=polylog(m-k,x,prec); p3=m2?greal(p3):gimag(p3);
! 969: tetpil=avma; y=gadd(y,gmul(p4,p3));
! 970: }
! 971: }
! 972: }
! 973: if (fl) { tetpil=avma; return gerepile(av,tetpil,gneg(y)); }
! 974: return gerepile(av,tetpil,y);
! 975: }
! 976:
! 977: GEN
! 978: gpolylog(long m, GEN x, long prec)
! 979: {
! 980: long i,lx,av=avma,tetpil,v,n;
! 981: GEN y,p1,p2;
! 982:
! 983: if (m<=0)
! 984: {
! 985: p1=polx[0]; p2=gsub(gun,p1);
! 986: for (i=1; i<=(-m); i++)
! 987: p1=gmul(polx[0],gadd(gmul(p2,derivpol(p1)),gmulsg(i,p1)));
! 988: p1=gdiv(p1,gpuigs(p2,1-m)); tetpil=avma;
! 989: return gerepile(av,tetpil,poleval(p1,x));
! 990: }
! 991:
! 992: switch(typ(x))
! 993: {
! 994: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
! 995: case t_COMPLEX: case t_QUAD:
! 996: return polylog(m,x,prec);
! 997:
! 998: case t_INTMOD: case t_PADIC: case t_POLMOD:
! 999: p1=roots((GEN)x[1],prec); lx=lg(p1); p2=cgetg(lx,t_COL);
! 1000: for (i=1; i<lx; i++) p2[i]=lpoleval((GEN)x[2],(GEN)p1[i]);
! 1001: tetpil=avma; y=cgetg(lx,t_COL);
! 1002: for (i=1; i<lx; i++) y[i]=(long)polylog(m,(GEN)p2[i],prec);
! 1003: return gerepile(av,tetpil,y);
! 1004:
! 1005: case t_POL: case t_RFRAC: case t_RFRACN:
! 1006: p1=tayl(x,gvar(x),precdl); tetpil=avma;
! 1007: return gerepile(av,tetpil,gpolylog(m,p1,prec));
! 1008:
! 1009: case t_SER:
! 1010: if (!m) return gneg(ghalf);
! 1011: if (m==1)
! 1012: {
! 1013: p1=glog(gsub(gun,x),prec); tetpil=avma;
! 1014: return gerepile(av,tetpil,gneg(p1));
! 1015: }
! 1016: if (valp(x)<=0) err(impl,"polylog around a!=0");
! 1017: v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2);
! 1018: for (i=n; i>=1; i--)
! 1019: {
! 1020: p1=gadd(gpuigs(stoi(i),-m),y); tetpil=avma;
! 1021: y=gmul(x,p1);
! 1022: }
! 1023: return gerepile(av,tetpil,y);
! 1024:
! 1025: case t_VEC: case t_COL: case t_MAT:
! 1026: lx=lg(x); y=cgetg(lx,typ(x));
! 1027: for (i=1; i<lx; i++)
! 1028: y[i]=(long)gpolylog(m,(GEN)x[i],prec);
! 1029: return y;
! 1030: }
! 1031: err(typeer,"gpolylog");
! 1032: return NULL; /* not reached */
! 1033: }
! 1034:
! 1035: void
! 1036: gpolylogz(long m, GEN x, GEN y)
! 1037: {
! 1038: long av=avma, prec = precision(y);
! 1039:
! 1040: if (!prec) err(infprecer,"gpolylogz");
! 1041: gaffect(gpolylog(m,x,prec),y); avma=av;
! 1042: }
! 1043:
! 1044: GEN
! 1045: polylog0(long m, GEN x, long flag, long prec)
! 1046: {
! 1047: switch(flag)
! 1048: {
! 1049: case 0: return gpolylog(m,x,prec);
! 1050: case 1: return polylogd(m,x,prec);
! 1051: case 2: return polylogdold(m,x,prec);
! 1052: case 3: return polylogp(m,x,prec);
! 1053: default: err(flagerr,"polylog");
! 1054: }
! 1055: return NULL; /* not reached */
! 1056: }
! 1057:
! 1058: static GEN
! 1059: qq(GEN x, long prec)
! 1060: {
! 1061: long tx=typ(x);
! 1062:
! 1063: if (tx==t_PADIC) return x;
! 1064: if (is_scalar_t(tx))
! 1065: {
! 1066: long l=precision(x);
! 1067: GEN p1,p2,q;
! 1068:
! 1069: if (tx != t_COMPLEX || gcmp((GEN)x[2],gzero)<=0)
! 1070: err(talker,"argument must belong to upper half-plane");
! 1071:
! 1072: if (!l) l=prec; p1=mppi(l); setexpo(p1,2);
! 1073: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */
! 1074: q=gmul(x,p2); return gexp(q,prec);
! 1075: }
! 1076: if (tx != t_POL && tx != t_SER)
! 1077: err(talker,"bad argument for modular function");
! 1078:
! 1079: return tayl(x,gvar(x),precdl);
! 1080: }
! 1081:
! 1082: static GEN
! 1083: inteta(GEN q)
! 1084: {
! 1085: long tx=typ(q);
! 1086: GEN p1,ps,qn,y;
! 1087:
! 1088: y=gun; qn=gun; ps=gun;
! 1089: if (tx==t_PADIC)
! 1090: {
! 1091: for(;;)
! 1092: {
! 1093: p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
! 1094: y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
! 1095: p1=y; y=gadd(y,ps); if (gegal(p1,y)) return y;
! 1096: }
! 1097: }
! 1098: else
! 1099: {
! 1100: long l,v, av = avma, lim = stack_lim(av,3);
! 1101:
! 1102: if (is_scalar_t(tx)) l = -bit_accuracy(precision(q));
! 1103: else
! 1104: { v=gvar(q); tx = 0; }
! 1105: for(;;)
! 1106: {
! 1107: p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
! 1108: y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
! 1109: y=gadd(y,ps);
! 1110: if (tx)
! 1111: { if (gexpo(ps)-gexpo(y) < l) return y; }
! 1112: else
! 1113: { if (gval(ps,v)>=precdl) return y; }
! 1114: if (low_stack(lim, stack_lim(av,3)))
! 1115: {
! 1116: GEN *gptr[3];
! 1117: if(DEBUGMEM>1) err(warnmem,"inteta");
! 1118: gptr[0]=&y; gptr[1]=&qn; gptr[2]=&ps;
! 1119: gerepilemany(av,gptr,3);
! 1120: }
! 1121: }
! 1122: }
! 1123: }
! 1124:
! 1125: GEN
! 1126: eta(GEN x, long prec)
! 1127: {
! 1128: long av = avma;
! 1129: GEN q = qq(x,prec);
! 1130: return gerepileupto(av,inteta(q));
! 1131: }
! 1132:
! 1133: /* returns the true value of eta(x) for Im(x) > 0, using reduction */
! 1134: GEN
! 1135: trueeta(GEN x, long prec)
! 1136: {
! 1137: long tx=typ(x), av=avma, tetpil,l;
! 1138: GEN p1,p2,q,q24,n,z,m,unapprox;
! 1139:
! 1140: if (!is_scalar_t(tx)) err(typeer,"trueeta");
! 1141: if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0)
! 1142: err(talker,"argument must belong to upper half-plane");
! 1143: l=precision(x); if (l) prec=l;
! 1144: p1=mppi(prec); setexpo(p1,2);
! 1145: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */
! 1146: z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */
! 1147: unapprox=gsub(gun,gpuigs(stoi(10),-8));
! 1148: m=gun;
! 1149: for(;;)
! 1150: {
! 1151: n=ground(greal(x));
! 1152: if (signe(n)) {x=gsub(x,n); m=gmul(m,powgi(z,n));}
! 1153: if (gcmp(gnorm(x), unapprox)>=0) break;
! 1154: m=gmul(m,gsqrt(gdiv(gi,x),prec)); x=gdivsg(-1,x);
! 1155: }
! 1156: q=gmul(p2,x);
! 1157: q24=gexp(gdivgs(q,24),prec); q=gpuigs(q24,24);
! 1158: p1=gmul(m,q24); q = inteta(q); tetpil=avma;
! 1159: return gerepile(av,tetpil,gmul(p1,q));
! 1160: }
! 1161:
! 1162: GEN
! 1163: eta0(GEN x, long flag,long prec)
! 1164: {
! 1165: return flag? trueeta(x,prec): eta(x,prec);
! 1166: }
! 1167:
! 1168: GEN
! 1169: jell(GEN x, long prec)
! 1170: {
! 1171: long av=avma,tetpil,tx = typ(x);
! 1172: GEN p1;
! 1173:
! 1174: if (!is_scalar_t(tx))
! 1175: {
! 1176: GEN p2,q = qq(x,prec);
! 1177: p1=gdiv(inteta(gsqr(q)), inteta(q));
! 1178: p1=gmul2n(gsqr(p1),1);
! 1179: p1=gmul(q,gpuigs(p1,12));
! 1180: p2=gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
! 1181: p1=gmulsg(48,p1); tetpil=avma;
! 1182: return gerepile(av,tetpil,gadd(p2,p1));
! 1183: }
! 1184: p1=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
! 1185: p1=gsqr(gsqr(gsqr(p1)));
! 1186: p1=gadd(gmul2n(gsqr(p1),8), ginv(p1)); tetpil=avma;
! 1187: return gerepile(av,tetpil,gpuigs(p1,3));
! 1188: }
! 1189:
! 1190: GEN
! 1191: wf2(GEN x, long prec)
! 1192: {
! 1193: long av=avma,tetpil;
! 1194: GEN p1,p2;
! 1195:
! 1196: p1=gsqrt(gdeux,prec);
! 1197: p2=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
! 1198: tetpil=avma;
! 1199: return gerepile(av,tetpil,gmul(p1,p2));
! 1200: }
! 1201:
! 1202: GEN
! 1203: wf1(GEN x, long prec)
! 1204: {
! 1205: long av=avma,tetpil;
! 1206: GEN p1,p2;
! 1207:
! 1208: p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec);
! 1209: tetpil=avma;
! 1210: return gerepile(av,tetpil,gdiv(p1,p2));
! 1211: }
! 1212:
! 1213: GEN
! 1214: wf(GEN x, long prec)
! 1215: {
! 1216: long av=avma,tetpil;
! 1217: GEN p1,p2;
! 1218:
! 1219: p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec));
! 1220: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=ldivrs(mppi(prec),-24);
! 1221: p2=gexp(p2,prec); tetpil=avma;
! 1222: return gerepile(av,tetpil,gmul(p1,p2));
! 1223: }
! 1224:
! 1225: GEN
! 1226: weber0(GEN x, long flag,long prec)
! 1227: {
! 1228: switch(flag)
! 1229: {
! 1230: case 0: return wf(x,prec);
! 1231: case 1: return wf1(x,prec);
! 1232: case 2: return wf2(x,prec);
! 1233: default: err(flagerr,"weber");
! 1234: }
! 1235: return NULL; /* not reached */
! 1236: }
! 1237:
! 1238: static GEN
! 1239: sagm(GEN x, long prec)
! 1240: {
! 1241: GEN p1,a,b,a1,b1,y;
! 1242: long av,tetpil,l,ep;
! 1243:
! 1244: if (gcmp0(x)) return gcopy(x);
! 1245: switch(typ(x))
! 1246: {
! 1247: case t_REAL:
! 1248: l = precision(x); y = cgetr(l); av=avma;
! 1249: a1 = x; b1 = realun(l);
! 1250: l = 5-bit_accuracy(prec);
! 1251: do
! 1252: {
! 1253: a=a1; b=b1; a1 = addrr(a,b);
! 1254: setexpo(a1,expo(a1)-1);
! 1255: b1=mpsqrt(mulrr(a,b));
! 1256: }
! 1257: while (expo(subrr(b1,a1))-expo(b1) >= l);
! 1258: affrr(a1,y); avma=av; return y;
! 1259:
! 1260: case t_COMPLEX:
! 1261: if (gcmp0((GEN)x[2]))
! 1262: return transc(sagm,(GEN)x[1],prec);
! 1263: av=avma; l=precision(x); if (l) prec=l;
! 1264: a1=x; b1=gun; l = 5-bit_accuracy(prec);
! 1265: do
! 1266: {
! 1267: a=a1; b=b1;
! 1268: a1=gmul2n(gadd(a,b),-1);
! 1269: b1=gsqrt(gmul(a,b),prec);
! 1270: }
! 1271: while (gexpo(gsub(b1,a1))-gexpo(b1) >= l);
! 1272: tetpil=avma; return gerepile(av,tetpil,gcopy(a1));
! 1273:
! 1274: case t_PADIC:
! 1275: av=avma; a1=x; b1=gun; l=precp(x);
! 1276: do
! 1277: {
! 1278: a=a1; b=b1;
! 1279: a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),0);
! 1280: p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
! 1281: if (ep<=0) { b1=gneg_i(b1); p1=gsub(b1,a1); ep=valp(p1)-valp(b1); }
! 1282: }
! 1283: while (ep<l && !gcmp0(p1));
! 1284: tetpil=avma; return gerepile(av,tetpil,gcopy(a1));
! 1285:
! 1286: case t_SER:
! 1287: av=avma; a1=x; b1=gun; l=lg(x)-2;
! 1288: do
! 1289: {
! 1290: a=a1; b=b1;
! 1291: a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),0);
! 1292: p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
! 1293: }
! 1294: while (ep<l && !gcmp0(p1));
! 1295: tetpil=avma; return gerepile(av,tetpil,gcopy(a1));
! 1296:
! 1297: case t_INTMOD:
! 1298: err(impl,"agm of mod");
! 1299: }
! 1300: return transc(sagm,x,prec);
! 1301: }
! 1302:
! 1303: GEN
! 1304: agm(GEN x, GEN y, long prec)
! 1305: {
! 1306: long av,tetpil, ty=typ(y);
! 1307: GEN z;
! 1308:
! 1309: if (is_matvec_t(ty))
! 1310: {
! 1311: ty=typ(x);
! 1312: if (is_matvec_t(ty)) err(talker,"agm of two vector/matrices");
! 1313: z=x; x=y; y=z;
! 1314: }
! 1315: if (gcmp0(y)) return gcopy(y);
! 1316: av=avma; z=sagm(gdiv(x,y),prec); tetpil=avma;
! 1317: return gerepile(av,tetpil,gmul(y,z));
! 1318: }
! 1319:
! 1320: GEN
! 1321: logagm(GEN q)
! 1322: {
! 1323: long av=avma,prec=lg(q),tetpil,s,n,lim;
! 1324: GEN y,q4,q1;
! 1325:
! 1326: if (typ(q)!=t_REAL) err(typeer,"logagm");
! 1327: if (signe(q)<=0) err(talker,"non positive argument in logagm");
! 1328: if (expo(q)<0) s=1; else { q=ginv(q); s=0; }
! 1329: lim = - (bit_accuracy(prec) >> 1);
! 1330: for (n=0; expo(q)>=lim ; n++) { q1=q; q=gsqr(q); }
! 1331: q4=gmul2n(q,2);
! 1332: if (!n) q1=gsqrt(q,prec);
! 1333: y=divrr(mppi(prec), agm(addsr(1,q4),gmul2n(q1,2),prec));
! 1334: tetpil=avma; y=gmul2n(y,-n); if (s) setsigne(y,-1);
! 1335: return gerepile(av,tetpil,y);
! 1336: }
! 1337:
! 1338: GEN
! 1339: glogagm(GEN x, long prec)
! 1340: {
! 1341: long av,tetpil;
! 1342: GEN y,p1,p2;
! 1343:
! 1344: switch(typ(x))
! 1345: {
! 1346: case t_REAL:
! 1347: if (signe(x)>=0) return logagm(x);
! 1348:
! 1349: y=cgetg(3,t_COMPLEX); y[2]=lmppi(lg(x));
! 1350: setsigne(x,1); y[1]=(long)logagm(x);
! 1351: setsigne(x,-1); return y;
! 1352:
! 1353: case t_COMPLEX:
! 1354: y=cgetg(3,t_COMPLEX); y[2]=larg(x,prec);
! 1355: av=avma; p1=glogagm(gnorm(x),prec); tetpil=avma;
! 1356: y[1]=lpile(av,tetpil,gmul2n(p1,-1));
! 1357: return y;
! 1358:
! 1359: case t_PADIC:
! 1360: return palog(x);
! 1361:
! 1362: case t_SER:
! 1363: av=avma; if (valp(x)) err(negexper,"glogagm");
! 1364: p1=gdiv(derivser(x),x);
! 1365: p1=integ(p1,varn(x));
! 1366: if (gcmp1((GEN)x[2])) return gerepileupto(av,p1);
! 1367: p2=glogagm((GEN)x[2],prec); tetpil=avma;
! 1368: return gerepile(av,tetpil,gadd(p1,p2));
! 1369:
! 1370: case t_INTMOD:
! 1371: err(typeer,"glogagm");
! 1372: }
! 1373: return transc(glogagm,x,prec);
! 1374: }
! 1375:
! 1376: GEN
! 1377: theta(GEN q, GEN z, long prec)
! 1378: {
! 1379: long av=avma,tetpil,l,n;
! 1380: GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold;
! 1381:
! 1382: if (gexpo(q)>=0) err(thetaer1);
! 1383: l=precision(q); if (l) prec=l;
! 1384: p1=realun(prec); z=gmul(p1,z);
! 1385: if (!l) q=gmul(p1,q);
! 1386: zy = gimag(z);
! 1387: if (gcmp0(zy)) k=gzero;
! 1388: else
! 1389: {
! 1390: lq=glog(q,prec); k=ground(gdiv(zy,greal(lq)));
! 1391: if (!gcmp0(k)) { zold=z; z=gadd(z,gdiv(gmul(lq,k),gi)); }
! 1392: }
! 1393: y=gsin(z,prec); n=0; qn=gun;
! 1394: ps2=gsqr(q); ps=gneg_i(ps2);
! 1395: do
! 1396: {
! 1397: n++; p1=gsin(gmulsg(2*n+1,z),prec);
! 1398: qnold=qn; qn=gmul(qn,ps);
! 1399: ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
! 1400: }
! 1401: while (gexpo(qnold) >= -bit_accuracy(prec));
! 1402: if (signe(k))
! 1403: {
! 1404: y=gmul(y,gmul(gpui(q,gsqr(k),prec),
! 1405: gexp(gmul2n(gmul(gmul(gi,zold),k),1),prec)));
! 1406: if (mod2(k)) y=gneg_i(y);
! 1407: }
! 1408: p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
! 1409: return gerepile(av,tetpil,gmul(p1,y));
! 1410: }
! 1411:
! 1412: GEN
! 1413: thetanullk(GEN q, long k, long prec)
! 1414: {
! 1415: long av=avma,tetpil,l,n;
! 1416: GEN p1,ps,qn,y,ps2;
! 1417:
! 1418: if (gexpo(q)>=0) err(thetaer1);
! 1419: if (!(k&1)) return gzero;
! 1420: n=0; qn=gun; ps2=gsqr(q); ps=gneg_i(ps2);
! 1421: y=gun; l=precision(q);
! 1422: if (!l) { l=prec; q=gmul(q,realun(l)); }
! 1423:
! 1424: do
! 1425: {
! 1426: n++; p1=gpuigs(stoi(n+n+1),k); qn=gmul(qn,ps);
! 1427: ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
! 1428: }
! 1429: while (gexpo(p1) >= -bit_accuracy(l));
! 1430: p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
! 1431: if (k&2) { p1=gneg_i(p1); tetpil=avma; }
! 1432: return gerepile(av,tetpil,gmul(p1,y));
! 1433: }
! 1434:
! 1435: GEN
! 1436: jbesselh(GEN n, GEN z, long prec)
! 1437: {
! 1438: long av,tetpil,k,l,i,lz;
! 1439: GEN y,p1,p2,zinv,p0,s,c;
! 1440:
! 1441: if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh");
! 1442: k=itos(n); if (k<0) err(impl,"ybesselh");
! 1443:
! 1444: switch(typ(z))
! 1445: {
! 1446: case t_REAL: case t_COMPLEX:
! 1447: if (gcmp0(z)) return gzero;
! 1448: av=avma; zinv=ginv(z);
! 1449: l=precision(z); if (l>prec) prec=l;
! 1450: gsincos(z,&s,&c,prec);
! 1451: p1=gmul(zinv,s);
! 1452: if (k)
! 1453: {
! 1454: p0=p1; p1=gmul(zinv,gsub(p0,c));
! 1455: for (i=2; i<=k; i++)
! 1456: {
! 1457: p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0);
! 1458: p0=p1; p1=p2;
! 1459: }
! 1460: }
! 1461: p2=gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec);
! 1462: tetpil=avma; return gerepile(av,tetpil,gmul(p2,p1));
! 1463:
! 1464: case t_VEC: case t_COL: case t_MAT:
! 1465: lz=lg(z); y=cgetg(lz,typ(z));
! 1466: for (i=1; i<lz; i++)
! 1467: y[i]=(long)jbesselh(n,(GEN)z[i],prec);
! 1468: return y;
! 1469:
! 1470: case t_INT: case t_FRAC: case t_FRACN:
! 1471: av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma;
! 1472: return gerepile(av,tetpil,jbesselh(n,p1,prec));
! 1473:
! 1474: case t_QUAD:
! 1475: av=avma; p1=gmul(z,realun(prec)); tetpil=avma;
! 1476: return gerepile(av,tetpil,jbesselh(n,p1,prec));
! 1477:
! 1478: case t_POL: case t_RFRAC: case t_RFRACN:
! 1479: av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma;
! 1480: return gerepile(av,tetpil,jbesselh(n,p1,prec));
! 1481:
! 1482: case t_POLMOD:
! 1483: av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL);
! 1484: for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]);
! 1485: tetpil=avma; y=cgetg(lz,t_COL);
! 1486: for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)p2[i],prec);
! 1487: return gerepile(av,tetpil,y);
! 1488:
! 1489: case t_PADIC:
! 1490: err(impl,"p-adic jbessel function");
! 1491: case t_SER:
! 1492: err(impl,"jbessel of power series");
! 1493: }
! 1494: err(typeer,"jbesselh");
! 1495: return NULL; /* not reached */
! 1496: }
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 / src / basemath