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