Annotation of OpenXM_contrib/pari/src/basemath/trans2.c, Revision 1.1
1.1 ! maekawa 1: /********************************************************************/
! 2: /** **/
! 3: /** TRANSCENDENTAL FONCTIONS **/
! 4: /** (part 2) **/
! 5: /** **/
! 6: /********************************************************************/
! 7: /* $Id: trans2.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
! 8: #include "pari.h"
! 9:
! 10: GEN mpsin(GEN x);
! 11: static GEN mpach(GEN x);
! 12:
! 13: /********************************************************************/
! 14: /** **/
! 15: /** FONCTION ARCTG **/
! 16: /** **/
! 17: /********************************************************************/
! 18:
! 19: static GEN
! 20: mpatan(GEN x)
! 21: {
! 22: long l,l1,l2,n,m,u,i,av0,av,lp,e,sx,s;
! 23: double alpha,beta,gama=1.0,delta,fi;
! 24: GEN y,p1,p2,p3,p4,p5,unr;
! 25:
! 26: sx=signe(x);
! 27: if (!sx)
! 28: {
! 29: y=cgetr(3); y[1]=x[1]; y[2]=0; return y;
! 30: }
! 31: l = lp = lg(x);
! 32: if (sx<0) setsigne(x,1);
! 33: u = cmprs(x,1);
! 34: if (!u)
! 35: {
! 36: y=mppi(l+1); setexpo(y,-1);
! 37: if (sx<0)
! 38: {
! 39: setsigne(x,-1);
! 40: setsigne(y,-1);
! 41: }
! 42: return y;
! 43: }
! 44:
! 45: e = expo(x);
! 46: if (e>0) lp += (e>>TWOPOTBITS_IN_LONG);
! 47:
! 48: y=cgetr(lp); av0=avma;
! 49: p1=cgetr(l+1); affrr(x,p1); setsigne(x,sx);
! 50: if (u==1) divsrz(1,p1,p1);
! 51: e = expo(p1);
! 52: if (e<-100)
! 53: alpha = log(PI)-e*LOG2;
! 54: else
! 55: {
! 56: alpha = rtodbl(p1);
! 57: alpha = log(PI/atan(alpha));
! 58: }
! 59: beta = (bit_accuracy(l)>>1) * LOG2;
! 60: delta=LOG2+beta-alpha/2;
! 61: if (delta<=0) { n=1; m=0; }
! 62: else
! 63: {
! 64: fi=alpha-2*LOG2;
! 65: if (delta>=gama*fi*fi/LOG2)
! 66: {
! 67: n=(long)(1+sqrt(gama*delta/LOG2));
! 68: m=(long)(1+sqrt(delta/(gama*LOG2))-fi/LOG2);
! 69: }
! 70: else
! 71: {
! 72: n=(long)(1+beta/fi); m=0;
! 73: }
! 74: }
! 75: l2=l+1+(m>>TWOPOTBITS_IN_LONG);
! 76: p2=cgetr(l2); p3=cgetr(l2);
! 77: affrr(p1,p2); av=avma;
! 78: for (i=1; i<=m; i++)
! 79: {
! 80: p5 = mulrr(p2,p2); setlg(p5,l2);
! 81: p5 = addsr(1,p5); setlg(p5,l2);
! 82: p5 = mpsqrt(p5);
! 83: p5 = addsr(1,p5); setlg(p5,l2);
! 84: divrrz(p2,p5,p2); avma=av;
! 85: }
! 86: mulrrz(p2,p2,p3); l1=4;
! 87: unr=realun(l2); setlg(unr,4);
! 88: p4=cgetr(l2); setlg(p4,4);
! 89: divrsz(unr,2*n+1,p4);
! 90:
! 91: s=0; e=expo(p3); av=avma;
! 92: for (i=n; i>=1; i--)
! 93: {
! 94: setlg(p3,l1); p5 = mulrr(p4,p3);
! 95: s -= e; l1 += (s>>TWOPOTBITS_IN_LONG);
! 96: if (l1>l2) l1=l2;
! 97: s %= BITS_IN_LONG;
! 98: setlg(unr,l1); p5 = subrr(divrs(unr,2*i-1), p5);
! 99: setlg(p4,l1); affrr(p5,p4); avma=av;
! 100: }
! 101: setlg(p4,l2); p4 = mulrr(p2,p4);
! 102: setexpo(p4, expo(p4)+m);
! 103: if (u==1)
! 104: {
! 105: p5 = mppi(lp+1); setexpo(p5,0);
! 106: p4 = subrr(p5,p4);
! 107: }
! 108: if (sx == -1) setsigne(p4,-signe(p4));
! 109: affrr(p4,y); avma=av0; return y;
! 110: }
! 111:
! 112: GEN
! 113: gatan(GEN x, long prec)
! 114: {
! 115: long av,tetpil;
! 116: GEN y,p1;
! 117:
! 118: switch(typ(x))
! 119: {
! 120: case t_REAL:
! 121: return mpatan(x);
! 122:
! 123: case t_COMPLEX:
! 124: av=avma; p1=cgetg(3,t_COMPLEX);
! 125: p1[1]=lneg((GEN)x[2]);
! 126: p1[2]=x[1]; tetpil=avma;
! 127: y=gerepile(av,tetpil,gath(p1,prec));
! 128: p1=(GEN)y[1]; y[1]=y[2]; y[2]=(long)p1;
! 129: setsigne(p1,-signe(p1)); return y;
! 130:
! 131: case t_SER:
! 132: av=avma; if (valp(x)<0) err(negexper,"gatan");
! 133: p1=gdiv(derivser(x), gaddsg(1,gsqr(x)));
! 134: y=integ(p1,varn(x)); if (valp(x)) return gerepileupto(av,y);
! 135:
! 136: p1=gatan((GEN)x[2],prec); tetpil=avma;
! 137: return gerepile(av,tetpil,gadd(p1,y));
! 138:
! 139: case t_INTMOD: case t_PADIC:
! 140: err(typeer,"gatan");
! 141: }
! 142: return transc(gatan,x,prec);
! 143: }
! 144:
! 145: void
! 146: gatanz(GEN x, GEN y)
! 147: {
! 148: long av=avma,prec = precision(y);
! 149:
! 150: if (!prec) err(infprecer,"gatanz");
! 151: gaffect(gatan(x,prec),y); avma=av;
! 152: }
! 153:
! 154: /********************************************************************/
! 155: /** **/
! 156: /** FONCTION ARCSINUS **/
! 157: /** **/
! 158: /********************************************************************/
! 159:
! 160: /* x is non zero |x| <= 1 */
! 161: static GEN
! 162: mpasin(GEN x)
! 163: {
! 164: long l,u,v,av;
! 165: GEN y,p1;
! 166:
! 167: u=cmprs(x,1); v=cmpsr(-1,x);
! 168: if (!u || !v)
! 169: {
! 170: y=mppi(lg(x)); setexpo(y,0);
! 171: if (signe(x)<0) setsigne(y,-1);
! 172: return y;
! 173: }
! 174: l = lg(x); y=cgetr(l); av=avma;
! 175:
! 176: p1 = cgetr(l+1); mulrrz(x,x,p1);
! 177: subsrz(1,p1,p1);
! 178: divrrz(x,mpsqrt(p1),p1);
! 179: affrr(mpatan(p1),y); if (signe(x)<0) setsigne(y,-1);
! 180: avma=av; return y;
! 181: }
! 182:
! 183: GEN
! 184: gasin(GEN x, long prec)
! 185: {
! 186: long av,tetpil,l,sx;
! 187: GEN y,p1;
! 188:
! 189: switch(typ(x))
! 190: {
! 191: case t_REAL: sx=signe(x);
! 192: if (!sx) { y=cgetr(3); y[1]=x[1]; y[2]=0; return y; }
! 193: if (sx<0) setsigne(x,1);
! 194: if (cmpsr(1,x)>=0) { setsigne(x,sx); return mpasin(x); }
! 195:
! 196: y=cgetg(3,t_COMPLEX);
! 197: y[1]=lmppi(lg(x)); setexpo(y[1],0);
! 198: y[2]=lmpach(x);
! 199: if (sx<0)
! 200: {
! 201: setsigne(y[1],-signe(y[1]));
! 202: setsigne(y[2],-signe(y[2]));
! 203: setsigne(x,sx);
! 204: }
! 205: return y;
! 206:
! 207: case t_COMPLEX:
! 208: av=avma; p1=cgetg(3,t_COMPLEX);
! 209: p1[1]=lneg((GEN)x[2]);
! 210: p1[2]=x[1]; tetpil=avma;
! 211: y=gerepile(av,tetpil,gash(p1,prec));
! 212: l=y[1]; y[1]=y[2];
! 213: y[2]=l; gnegz((GEN)l,(GEN)l);
! 214: return y;
! 215:
! 216: case t_SER:
! 217: if (gcmp0(x)) return gcopy(x);
! 218:
! 219: av=avma; if (valp(x)<0) err(negexper,"gasin");
! 220: p1=gdiv(derivser(x), gsqrt(gsubsg(1,gsqr(x)),prec));
! 221: y=integ(p1,varn(x)); if (valp(x)) return gerepileupto(av,y);
! 222:
! 223: p1=gasin((GEN)x[2],prec); tetpil=avma;
! 224: return gerepile(av,tetpil,gadd(p1,y));
! 225:
! 226: case t_INTMOD: case t_PADIC:
! 227: err(typeer,"gasin");
! 228: }
! 229: return transc(gasin,x,prec);
! 230: }
! 231:
! 232: void
! 233: gasinz(GEN x, GEN y)
! 234: {
! 235: long av=avma, prec = precision(y);
! 236:
! 237: if (!prec) err(infprecer,"gasinz");
! 238: gaffect(gasin(x,prec),y); avma=av;
! 239: }
! 240:
! 241: /********************************************************************/
! 242: /** **/
! 243: /** FONCTION ARCCOSINUS **/
! 244: /** **/
! 245: /********************************************************************/
! 246:
! 247: /* |x|<=1 */
! 248: static GEN
! 249: mpacos(GEN x)
! 250: {
! 251: long l,u,v,sx,av;
! 252: GEN y,p1,p2;
! 253:
! 254: u=cmprs(x,1); v=cmpsr(-1,x); sx = signe(x);
! 255: if (!sx)
! 256: {
! 257: l=expo(x)>>TWOPOTBITS_IN_LONG; if (l>=0) l = -1;
! 258: y=mppi(2-l); setexpo(y,0); return y;
! 259: }
! 260: l=lg(x);
! 261: if (!u)
! 262: {
! 263: y = cgetr(3);
! 264: y[1] = evalexpo(-(bit_accuracy(l)>>1));
! 265: y[2] = 0; return y;
! 266: }
! 267: if (!v) return mppi(l);
! 268: y=cgetr(l); av=avma;
! 269:
! 270: p1=cgetr(l+1);
! 271: if (expo(x)<0)
! 272: {
! 273: mulrrz(x,x,p1);
! 274: subsrz(1,p1,p1);
! 275: p1 = mpsqrt(p1); divrrz(x,p1,p1);
! 276: p1 = mpatan(p1);
! 277: p2 = mppi(l); setexpo(p2,0);
! 278: p1 = subrr(p2,p1);
! 279: }
! 280: else
! 281: {
! 282: p2=cgetr(l+1);
! 283: if (sx>0) addsrz(1,x,p1); else subsrz(1,x,p1);
! 284: subsrz(2,p1,p2);
! 285: mulrrz(p1,p2,p1);
! 286: p1 = mpsqrt(p1); divrrz(p1,x,p1);
! 287: p1 = mpatan(p1);
! 288: if (sx<0) { setlg(p1,l); p1 = addrr(mppi(l),p1); }
! 289: }
! 290: affrr(p1,y); avma=av; return y;
! 291: }
! 292:
! 293: GEN
! 294: gacos(GEN x, long prec)
! 295: {
! 296: long av,tetpil,l,sx;
! 297: GEN y,p1;
! 298:
! 299: switch(typ(x))
! 300: {
! 301: case t_REAL: sx=signe(x);
! 302: if (sx<0) setsigne(x,1);
! 303: if (cmprs(x,1)<=0) { setsigne(x,sx); return mpacos(x); }
! 304:
! 305: y=cgetg(3,t_COMPLEX); y[2]=lmpach(x);
! 306: if (sx<0) y[1]=lmppi(lg(x));
! 307: else
! 308: {
! 309: y[1]=zero; setsigne(y[2],-signe(y[2]));
! 310: }
! 311: setsigne(x,sx); return y;
! 312:
! 313: case t_COMPLEX:
! 314: y=gach(x,prec);
! 315: l=y[1]; y[1]=y[2]; y[2]=l;
! 316: setsigne(y[2],-signe(y[2])); return y;
! 317:
! 318: case t_SER: av=avma;
! 319: if (valp(x)<0) err(negexper,"gacos");
! 320: p1=integ(gdiv(derivser(x), gsqrt(gsubsg(1,gsqr(x)),prec)),varn(x));
! 321: if (gcmp1((GEN)x[2]) && !valp(x))
! 322: {
! 323: tetpil=avma; return gerepile(av,tetpil,gneg(p1));
! 324: }
! 325: if (valp(x)) { y=mppi(prec); setexpo(y,0); }
! 326: else y=gacos((GEN)x[2],prec);
! 327: tetpil=avma; return gerepile(av,tetpil,gsub(y,p1));
! 328:
! 329: case t_INTMOD: case t_PADIC:
! 330: err(typeer,"gacos");
! 331: }
! 332: return transc(gacos,x,prec);
! 333: }
! 334:
! 335: void
! 336: gacosz(GEN x, GEN y)
! 337: {
! 338: long av=avma, prec = precision(y);
! 339:
! 340: if (!prec) err(infprecer,"gacosz");
! 341: gaffect(gacos(x,prec),y); avma=av;
! 342: }
! 343:
! 344: /********************************************************************/
! 345: /** **/
! 346: /** FONCTION ARGUMENT **/
! 347: /** **/
! 348: /********************************************************************/
! 349:
! 350: /* we know that x and y are not both 0 */
! 351: static GEN
! 352: mparg(GEN x, GEN y)
! 353: {
! 354: long prec,sx,sy;
! 355: GEN theta,pitemp;
! 356:
! 357: sx=signe(x); sy=signe(y);
! 358: if (!sy)
! 359: {
! 360: if (sx>0)
! 361: {
! 362: theta=cgetr(3); theta[1]=y[1]-expo(x);
! 363: theta[2]=0; return theta;
! 364: }
! 365: return mppi(lg(x));
! 366: }
! 367: prec = lg(y); if (prec<lg(x)) prec = lg(x);
! 368: if (!sx)
! 369: {
! 370: theta=mppi(prec); setexpo(theta,0);
! 371: if (sy<0) setsigne(theta,-1);
! 372: return theta;
! 373: }
! 374:
! 375: if (expo(x)-expo(y) > -2)
! 376: {
! 377: theta = mpatan(divrr(y,x));
! 378: if (sx>0) return theta;
! 379: pitemp=mppi(prec);
! 380: if (sy>0) return addrr(pitemp,theta);
! 381: return subrr(theta,pitemp);
! 382: }
! 383: theta = mpatan(divrr(x,y));
! 384: pitemp=mppi(prec); setexpo(pitemp,0);
! 385: if (sy>0) return subrr(pitemp,theta);
! 386:
! 387: theta = addrr(pitemp,theta);
! 388: setsigne(theta,-signe(theta)); return theta;
! 389: }
! 390:
! 391: static GEN
! 392: rfix(GEN x,long prec)
! 393: {
! 394: GEN p1;
! 395: switch(typ(x))
! 396: {
! 397: case t_INT: case t_FRAC: case t_FRACN:
! 398: p1=cgetr(prec); gaffect(x,p1); return p1;
! 399: }
! 400: return x;
! 401: }
! 402:
! 403: static GEN
! 404: sarg(GEN x, GEN y, long prec)
! 405: {
! 406: long av=avma;
! 407: x = rfix(x,prec); y = rfix(y,prec);
! 408: return gerepileupto(av,mparg(x,y));
! 409: }
! 410:
! 411: GEN
! 412: garg(GEN x, long prec)
! 413: {
! 414: GEN p1;
! 415: long av,tx=typ(x),tetpil;
! 416:
! 417: if (gcmp0(x)) err(talker,"zero argument in garg");
! 418: switch(tx)
! 419: {
! 420: case t_REAL:
! 421: prec=lg(x); /* fall through */
! 422: case t_INT: case t_FRAC: case t_FRACN:
! 423: return (gsigne(x)>0)? realzero(prec): mppi(prec);
! 424:
! 425: case t_QUAD:
! 426: av=avma; gaffsg(1,p1=cgetr(prec)); p1=gmul(p1,x);
! 427: tetpil=avma; return gerepile(av,tetpil,garg(p1,prec));
! 428:
! 429: case t_COMPLEX:
! 430: return sarg((GEN)x[1],(GEN)x[2],prec);
! 431:
! 432: case t_VEC: case t_COL: case t_MAT:
! 433: return transc(garg,x,prec);
! 434: }
! 435: err(typeer,"garg");
! 436: return NULL; /* not reached */
! 437: }
! 438:
! 439: /********************************************************************/
! 440: /** **/
! 441: /** FONCTION COSINUS HYPERBOLIQUE **/
! 442: /** **/
! 443: /********************************************************************/
! 444:
! 445: static GEN
! 446: mpch(GEN x)
! 447: {
! 448: long l,av;
! 449: GEN y,p1;
! 450:
! 451: if (gcmp0(x)) return gaddsg(1,x);
! 452:
! 453: l=lg(x); y=cgetr(l); av=avma;
! 454: p1=mpexp(x); p1 = addrr(p1, divsr(1,p1));
! 455: setexpo(p1, expo(p1)-1);
! 456: affrr(p1,y); avma=av; return y;
! 457: }
! 458:
! 459: GEN
! 460: gch(GEN x, long prec)
! 461: {
! 462: long av,tetpil;
! 463: GEN p1;
! 464:
! 465: switch(typ(x))
! 466: {
! 467: case t_REAL:
! 468: return mpch(x);
! 469:
! 470: case t_COMPLEX:
! 471: av=avma; p1=gexp(x,prec);
! 472: p1=gadd(p1,ginv(p1)); tetpil=avma;
! 473: return gerepile(av,tetpil,gmul2n(p1,-1));
! 474:
! 475: case t_SER:
! 476: av=avma; p1=gexp(x,prec); p1=gadd(p1,gdivsg(1,p1));
! 477: tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-1));
! 478:
! 479: case t_INTMOD: case t_PADIC:
! 480: err(typeer,"gch");
! 481: }
! 482: return transc(gch,x,prec);
! 483: }
! 484:
! 485: void
! 486: gchz(GEN x, GEN y)
! 487: {
! 488: long av=avma,prec = precision(y);
! 489:
! 490: if (!prec) err(infprecer,"gchz");
! 491: gaffect(gch(x,prec),y); avma=av;
! 492: }
! 493:
! 494: /********************************************************************/
! 495: /** **/
! 496: /** FONCTION SINUS HYPERBOLIQUE **/
! 497: /** **/
! 498: /********************************************************************/
! 499:
! 500: static GEN
! 501: mpsh(GEN x)
! 502: {
! 503: long l,av;
! 504: GEN y,p1;
! 505:
! 506: if (!signe(x))
! 507: {
! 508: y=cgetr(3); y[1]=x[1]; y[2]=0; return y;
! 509: }
! 510: l=lg(x); y=cgetr(l); av=avma;
! 511: p1=mpexp(x); p1 = addrr(p1, divsr(-1,p1));
! 512: setexpo(p1, expo(p1)-1);
! 513: affrr(p1,y); avma=av; return y;
! 514: }
! 515:
! 516: GEN
! 517: gsh(GEN x, long prec)
! 518: {
! 519: long av,tetpil;
! 520: GEN p1;
! 521:
! 522: switch(typ(x))
! 523: {
! 524: case t_REAL:
! 525: return mpsh(x);
! 526:
! 527: case t_COMPLEX:
! 528: av=avma; p1=gexp(x,prec);
! 529: p1=gsub(p1,ginv(p1)); tetpil=avma;
! 530: return gerepile(av,tetpil,gmul2n(p1,-1));
! 531:
! 532: case t_SER:
! 533: if (gcmp0(x)) return gcopy(x);
! 534:
! 535: av=avma; p1=gexp(x,prec); p1=gsub(p1,gdivsg(1,p1));
! 536: tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-1));
! 537:
! 538: case t_INTMOD: case t_PADIC:
! 539: err(typeer,"gsh");
! 540: }
! 541: return transc(gsh,x,prec);
! 542: }
! 543:
! 544: void
! 545: gshz(GEN x, GEN y)
! 546: {
! 547: long av=avma, prec = precision(y);
! 548:
! 549: if (!prec) err(infprecer,"gshz");
! 550: gaffect(gsh(x,prec),y); avma=av;
! 551: }
! 552:
! 553: /********************************************************************/
! 554: /** **/
! 555: /** FONCTION TANGENTE HYPERBOLIQUE **/
! 556: /** **/
! 557: /********************************************************************/
! 558:
! 559: static GEN
! 560: mpth(GEN x)
! 561: {
! 562: long l,av;
! 563: GEN y,p1,p2;
! 564:
! 565: if (!signe(x))
! 566: {
! 567: y=cgetr(3); y[1]=x[1]; y[2]=0;
! 568: return y;
! 569: }
! 570: l=lg(x); y=cgetr(l); av=avma;
! 571:
! 572: p1=cgetr(l+1); affrr(x,p1);
! 573: setexpo(p1,expo(p1)+1);
! 574: p1 = mpexp1(p1);
! 575: p2 = addsr(2,p1); setlg(p2,l+1);
! 576: p1 = divrr(p1,p2);
! 577: affrr(p1,y); avma=av; return y;
! 578: }
! 579:
! 580: GEN
! 581: gth(GEN x, long prec)
! 582: {
! 583: long av,tetpil;
! 584: GEN p1,p2,p3;
! 585:
! 586: switch(typ(x))
! 587: {
! 588: case t_REAL:
! 589: return mpth(x);
! 590:
! 591: case t_COMPLEX:
! 592: av=avma; p1=gexp(gmul2n(x,1),prec);
! 593: p1=gdivsg(-2,gaddgs(p1,1)); tetpil=avma;
! 594: return gerepile(av,tetpil,gaddsg(1,p1));
! 595:
! 596: case t_SER:
! 597: if (gcmp0(x)) return gcopy(x);
! 598:
! 599: av=avma; p1=gexp(gmul2n(x ,1),prec);
! 600: p2=gsubgs(p1,1); p3=gaddgs(p1,1);
! 601: tetpil=avma; return gerepile(av,tetpil,gdiv(p2,p3));
! 602:
! 603: case t_INTMOD: case t_PADIC:
! 604: err(typeer,"gth");
! 605: }
! 606: return transc(gth,x,prec);
! 607: }
! 608:
! 609: void
! 610: gthz(GEN x, GEN y)
! 611: {
! 612: long av=avma, prec = precision(y);
! 613:
! 614: if (!prec) err(infprecer,"gthz");
! 615: gaffect(gth(x,prec),y); avma=av;
! 616: }
! 617:
! 618: /********************************************************************/
! 619: /** **/
! 620: /** FONCTION ARGUMENT SINUS HYPERBOLIQUE **/
! 621: /** **/
! 622: /********************************************************************/
! 623:
! 624: /* x is non-zero */
! 625: static GEN
! 626: mpash(GEN x)
! 627: {
! 628: long s=signe(x),av;
! 629: GEN y,p1;
! 630:
! 631: y=cgetr(lg(x)); av=avma;
! 632: p1 = (s<0)? negr(x): x;
! 633: p1 = addrr(p1,mpsqrt(addsr(1,mulrr(p1,p1))));
! 634: p1 = mplog(p1); if (s<0) setsigne(p1,-signe(p1));
! 635: affrr(p1,y); avma=av; return y;
! 636: }
! 637:
! 638: GEN
! 639: gash(GEN x, long prec)
! 640: {
! 641: long av,tetpil,sx,sy,sz;
! 642: GEN y,p1;
! 643:
! 644: if (gcmp0(x)) return gcopy(x);
! 645: switch(typ(x))
! 646: {
! 647: case t_REAL:
! 648: return mpash(x);
! 649:
! 650: case t_COMPLEX:
! 651: av=avma; p1=gaddsg(1,gsqr(x));
! 652: p1=gadd(x,gsqrt(p1,prec));
! 653: tetpil=avma; y=glog(p1,prec);
! 654: sz=gsigne((GEN)y[1]);
! 655: sx=gsigne((GEN)p1[1]);
! 656: sy=gsigne((GEN)p1[2]);
! 657: if (sx>0 || (!sx && sy*sz<=0)) return gerepile(av,tetpil,y);
! 658:
! 659: y=gneg_i(y); p1=cgetg(3,t_COMPLEX); p1[1]=zero;
! 660: p1[2]=lmppi(prec); if (sy<0) setsigne(p1[2],-1);
! 661: tetpil=avma;
! 662: return gerepile(av,tetpil,gadd(y,p1));
! 663:
! 664: case t_SER:
! 665: if (gcmp0(x)) return gcopy(x);
! 666: if (valp(x)<0) err(negexper,"gash");
! 667:
! 668: av=avma; p1=gdiv(derivser(x),gsqrt(gaddsg(1,gsqr(x)),0));
! 669: y = integ(p1,varn(x)); if (valp(x)) return gerepileupto(av,y);
! 670: p1=gash((GEN)x[2],prec);
! 671: tetpil=avma; return gerepile(av,tetpil,gadd(p1,y));
! 672:
! 673: case t_INTMOD: case t_PADIC:
! 674: err(typeer,"gash");
! 675: }
! 676: return transc(gash,x,prec);
! 677: }
! 678:
! 679: void
! 680: gashz(GEN x, GEN y)
! 681: {
! 682: long av=avma, prec = precision(y);
! 683:
! 684: if (!prec) err(infprecer,"gashz");
! 685: gaffect(gash(x,prec),y); avma=av;
! 686: }
! 687:
! 688: /********************************************************************/
! 689: /** **/
! 690: /** FONCTION ARGUMENT COSINUS HYPERBOLIQUE **/
! 691: /** **/
! 692: /********************************************************************/
! 693:
! 694: static GEN
! 695: mpach(GEN x)
! 696: {
! 697: long l,av;
! 698: GEN y,p1;
! 699:
! 700: l=lg(x); y=cgetr(l); av=avma;
! 701:
! 702: p1=cgetr(l+1); affrr(x,p1);
! 703: p1 = mulrr(p1,p1);
! 704: subrsz(p1,1,p1);
! 705: p1 = mpsqrt(p1);
! 706: p1 = mplog(addrr(x,p1));
! 707: affrr(p1,y); avma=av; return y;
! 708: }
! 709:
! 710: GEN
! 711: gach(GEN x, long prec)
! 712: {
! 713: long av,tetpil;
! 714: GEN y,p1;
! 715:
! 716: switch(typ(x))
! 717: {
! 718: case t_REAL:
! 719: if (gcmpgs(x,1)>=0) return mpach(x);
! 720:
! 721: y=cgetg(3,t_COMPLEX);
! 722: if (gcmpgs(x,-1)>=0)
! 723: {
! 724: y[2]=lmpacos(x); y[1]=zero;
! 725: return y;
! 726: }
! 727: av=avma; p1=mpach(gneg_i(x)); tetpil=avma;
! 728: y[1]=lpile(av,tetpil,gneg(p1));
! 729: y[2]=lmppi(lg(x));
! 730: return y;
! 731:
! 732: case t_COMPLEX:
! 733: av=avma; p1=gaddsg(-1,gsqr(x));
! 734: p1=gadd(x,gsqrt(p1,prec)); tetpil=avma;
! 735: y=glog(p1,prec);
! 736: if (signe(y[2])<0) { tetpil=avma; y=gneg(y); }
! 737: return gerepile(av,tetpil,y);
! 738:
! 739: case t_SER:
! 740: av=avma; if (valp(x)<0) err(negexper,"gach");
! 741: p1=gdiv(derivser(x),gsqrt(gsubgs(gsqr(x),1),prec));
! 742: y=integ(p1,varn(x));
! 743: if (!valp(x) && gcmp1((GEN)x[2])) return gerepileupto(av,y);
! 744: if (valp(x))
! 745: {
! 746: p1=cgetg(3,t_COMPLEX); p1[1]=zero; p1[2]=lmppi(prec);
! 747: setexpo(p1[2],0);
! 748: }
! 749: else p1=gach((GEN)x[2],prec);
! 750: tetpil=avma; return gerepile(av,tetpil,gadd(p1,y));
! 751:
! 752: case t_INTMOD: case t_PADIC:
! 753: err(typeer,"gach");
! 754: }
! 755: return transc(gach,x,prec);
! 756: }
! 757:
! 758: void
! 759: gachz(GEN x, GEN y)
! 760: {
! 761: long av=avma,prec = precision(y);
! 762:
! 763: if (!prec) err(infprecer,"gachz");
! 764: gaffect(gach(x,prec),y); avma=av;
! 765: }
! 766:
! 767: /********************************************************************/
! 768: /** **/
! 769: /** FONCTION ARGUMENT TANGENTE HYPERBOLIQUE **/
! 770: /** **/
! 771: /********************************************************************/
! 772:
! 773: static GEN
! 774: mpath(GEN x)
! 775: {
! 776: long av;
! 777: GEN y,p1;
! 778:
! 779: if (!signe(x))
! 780: {
! 781: y=cgetr(3); y[1]=x[1]; y[2]=0;
! 782: return y;
! 783: }
! 784: y=cgetr(lg(x)); av=avma;
! 785: p1 = addrs(divsr(2,subsr(1,x)),-1);
! 786: affrr(mplog(p1), y); avma=av;
! 787: setexpo(y,expo(y)-1); return y;
! 788: }
! 789:
! 790: GEN
! 791: gath(GEN x, long prec)
! 792: {
! 793: long av,tetpil;
! 794: GEN y,p1;
! 795:
! 796: switch(typ(x))
! 797: {
! 798: case t_REAL:
! 799: if (expo(x)<0) return mpath(x);
! 800:
! 801: av=avma; p1=addrs(divsr(2,addsr(-1,x)),1);
! 802: tetpil=avma; y=cgetg(3,t_COMPLEX);
! 803: p1=mplog(p1); setexpo(p1,expo(p1)-1);
! 804: y[1]=(long)p1;
! 805: y[2]=lmppi(lg(x)); setexpo(y[2],0);
! 806: return gerepile(av,tetpil,y);
! 807:
! 808: case t_COMPLEX:
! 809: av=avma;
! 810: p1=gaddgs(gdivsg(2,gsubsg(1,x)),-1);
! 811: p1=glog(p1,prec); tetpil=avma;
! 812: return gerepile(av,tetpil,gmul2n(p1,-1));
! 813:
! 814: case t_SER:
! 815: av=avma; if (valp(x)<0) err(negexper,"gath");
! 816: p1=gdiv(derivser(x), gsubsg(1,gsqr(x)));
! 817: y = integ(p1,varn(x)); if (valp(x)) return gerepileupto(av,y);
! 818: p1=gath((GEN)x[2],prec); tetpil=avma;
! 819: return gerepile(av,tetpil,gadd(p1,y));
! 820:
! 821: case t_INTMOD: case t_PADIC:
! 822: err(typeer,"gath");
! 823: }
! 824: return transc(gath,x,prec);
! 825: }
! 826:
! 827: void
! 828: gathz(GEN x, GEN y)
! 829: {
! 830: long av=avma, prec = precision(y);
! 831:
! 832: if (!prec) err(infprecer,"gathz");
! 833: gaffect(gath(x,prec),y); avma=av;
! 834: }
! 835:
! 836: /********************************************************************/
! 837: /** **/
! 838: /** FONCTION TABLEAU DES NOMBRES DE BERNOULLI **/
! 839: /** **/
! 840: /********************************************************************/
! 841:
! 842: /* pour calculer B_0,B_2,...,B_2*nb */
! 843: void
! 844: mpbern(long nb, long prec)
! 845: {
! 846: long n,m,i,j,d,d1,d2,l,av,av2,code0;
! 847: GEN p1,p2;
! 848:
! 849: if (nb<0) nb=0;
! 850: if (bernzone)
! 851: {
! 852: if (bernzone[1]>=nb && bernzone[2]>=prec) return;
! 853: gunclone(bernzone);
! 854: }
! 855: d = 3 + prec*(nb+1); bernzone=newbloc(d);
! 856: bernzone[0]=evallg(d);
! 857: bernzone[1]=nb;
! 858: bernzone[2]=prec;
! 859: av=avma; l = prec+1; p1=realun(l);
! 860:
! 861: code0 = evaltyp(t_REAL) | evallg(prec);
! 862: *(bern(0)) = code0; affsr(1,bern(0));
! 863: p2 = p1; av2=avma;
! 864: for (i=1; i<=nb; i++)
! 865: {
! 866: if (i>1)
! 867: {
! 868: n=8; m=5; d = d1 = i-1; d2 = 2*i-3;
! 869: for (j=d; j>0; j--)
! 870: {
! 871: p2 = bern(j);
! 872: if (j<d) p2 = addrr(p2,p1);
! 873: else
! 874: {
! 875: affrr(p2,p1); p2=p1;
! 876: }
! 877: p2 = mulsr(n*m,p2); setlg(p2,l);
! 878: p2 = divrs(p2, d1*d2); affrr(p2,p1);
! 879: n+=4; m+=2; d1--; d2-=2;
! 880: }
! 881: p2 = addsr(1,p1); setlg(p2,l);
! 882: }
! 883: p2 = divrs(p2,2*i+1); p2 = subsr(1,p2);
! 884: setexpo(p2, expo(p2) - 2*i);
! 885: *(bern(i)) = code0; affrr(p2,bern(i)); avma=av2;
! 886: }
! 887: avma=av;
! 888: }
! 889:
! 890: GEN
! 891: bernreal(long n, long prec)
! 892: {
! 893: GEN B;
! 894:
! 895: if (n==1) { B=cgetr(prec); affsr(-1,B); setexpo(B,-1); return B; }
! 896: if (n<0 || n&1) return gzero;
! 897: n >>= 1; mpbern(n+1,prec); B=cgetr(prec);
! 898: affrr(bern(n),B); return B;
! 899: }
! 900:
! 901: /* k > 0 */
! 902: static GEN
! 903: bernfracspec(long k)
! 904: {
! 905: long n,av,lim, K = k+1;
! 906: GEN s,c,N,h;
! 907:
! 908: c = N = stoi(K); s = gun; h = gzero;
! 909: av = avma; lim = stack_lim(av,2);
! 910: for (n=2; ; n++)
! 911: {
! 912: c = gdivgs(gmulgs(c,n-k-2),n);
! 913: h = gadd(h, gdivgs(gmul(c,s), n));
! 914: if (n==K) return gerepileupto(av,h);
! 915: N[2] = n; s = addii(s, gpuigs(N,k));
! 916: if (low_stack(lim, stack_lim(av,2)))
! 917: {
! 918: GEN *gptr[3]; gptr[0]=&c; gptr[1]=&h; gptr[2]=&s;
! 919: if (DEBUGMEM>1) err(warnmem,"bernfrac");
! 920: gerepilemany(av,gptr,3);
! 921: }
! 922: }
! 923: }
! 924:
! 925: GEN
! 926: bernfrac(long k)
! 927: {
! 928: if (!k) return gun;
! 929: if (k==1) return gneg(ghalf);
! 930: if (k<0 || k&1) return gzero;
! 931: return bernfracspec(k);
! 932: }
! 933:
! 934: static GEN
! 935: bernvec2(long k)
! 936: {
! 937: GEN B = cgetg(k+2,t_VEC);
! 938: ulong i;
! 939:
! 940: for (i=1; i<=k; i++)
! 941: B[i+1]=(long)bernfracspec(i<<1);
! 942: B[1]=un; return B;
! 943: }
! 944:
! 945: /* mpbern as exact fractions */
! 946: GEN
! 947: bernvec(long nb)
! 948: {
! 949: long n,m,i,j,d1,d2,av,tetpil;
! 950: GEN p1,y;
! 951:
! 952: if (nb < 300) return bernvec2(nb);
! 953:
! 954: y=cgetg(nb+2,t_VEC); y[1]=un;
! 955: for (i=1; i<=nb; i++)
! 956: {
! 957: av=avma; p1=gzero;
! 958: n=8; m=5; d1=i-1; d2=2*i-3;
! 959: for (j=d1; j>0; j--)
! 960: {
! 961: p1 = gmulsg(n*m, gadd(p1,(GEN)y[j+1]));
! 962: p1 = gdivgs(p1, d1*d2);
! 963: n+=4; m+=2; d1--; d2-=2;
! 964: }
! 965: p1 = gsubsg(1,gdivgs(gaddsg(1,p1),2*i+1));
! 966: tetpil=avma; p1 = gmul2n(p1,-2*i);
! 967: y[i+1] = lpile(av,tetpil,p1);
! 968: }
! 969: return y;
! 970: }
! 971:
! 972: /********************************************************************/
! 973: /** **/
! 974: /** FONCTION GAMMA **/
! 975: /** **/
! 976: /********************************************************************/
! 977:
! 978: static GEN
! 979: mpgamma(GEN x)
! 980: {
! 981: GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;
! 982: long l,l1,l2,u,i,k,e,s,sx,n,p,av,av1;
! 983: double alpha,beta,dk;
! 984:
! 985: sx=signe(x); if (!sx) err(gamer2);
! 986: l=lg(x); y=cgetr(l); av=avma;
! 987:
! 988: l2=l+1; p1=cgetr(l2);
! 989: u = (expo(x)<-1 || sx<0);
! 990: if (!u) p2 = x;
! 991: else
! 992: {
! 993: p2=gfrac(x); if (gcmp0(p2)) err(gamer2);
! 994: p2 = subsr(1,x);
! 995: }
! 996: affrr(p2,p1);
! 997: alpha=rtodbl(p1); beta=((bit_accuracy(l)>>1)*LOG2/PI) - alpha;
! 998: if (beta>=0) n=(long)(1 + pariK2*beta); else n=0;
! 999: if (n)
! 1000: {
! 1001: p = (long)(1 + PI*(alpha+n));
! 1002: l2 += n>>TWOPOTBITS_IN_LONG;
! 1003: p2 = cgetr(l2); addsrz(n,p1,p2);
! 1004: }
! 1005: else
! 1006: {
! 1007: dk = pariK4*alpha/(l-2); beta = log(dk)/LOG2;
! 1008: if (beta>1.) beta += log(beta)/LOG2;
! 1009: p = (long)((bit_accuracy(l)>>1)/beta + 1);
! 1010: p2 = p1;
! 1011: }
! 1012: mpbern(p,l2); p3=mplog(p2);
! 1013:
! 1014: p4=realun(l2); setexpo(p4,-1);
! 1015: p6 = subrr(p2,p4); p6 = mulrr(p6,p3);
! 1016: p6 = subrr(p6,p2);
! 1017: pitemp = mppi(l2); setexpo(pitemp,2);
! 1018: p7 = mplog(pitemp); setexpo(pitemp,1);
! 1019: setexpo(p7,-1); addrrz(p6,p7,p4);
! 1020:
! 1021: affrr(ginv(gsqr(p2)), p3); e=expo(p3);
! 1022:
! 1023: p5=cgetr(l2); setlg(p5,4);
! 1024: p71=cgetr(l2); p7 = bern(p);
! 1025: if (bernzone[2]>4) { setlg(p71,4); affrr(p7,p71); p7=p71; }
! 1026: p7 = divrs(p7, 2*p*(2*p-1)); affrr(p7,p5);
! 1027:
! 1028: s=0; l1=4; av1=avma;
! 1029: for (k=p-1; k>0; k--)
! 1030: {
! 1031: setlg(p3,l1); p6 = mulrr(p3,p5); p7 = bern(k);
! 1032: if (bernzone[2]>l1) { setlg(p71,l1); affrr(p7,p71); p7=p71; }
! 1033: p7 = divrs(p7, (2*k)*(2*k-1));
! 1034: s -= e; l1 += (s>>TWOPOTBITS_IN_LONG); if (l1>l2) l1=l2;
! 1035: s &= (BITS_IN_LONG-1); p7 = addrr(p7,p6);
! 1036: setlg(p5,l1); affrr(p7,p5); avma=av1;
! 1037: }
! 1038: setlg(p5,l2); p6 = divrr(p5,p2);
! 1039: p4 = addrr(p4,p6); setlg(p4,l2); p4 = mpexp(p4);
! 1040: for (i=1; i<=n; i++)
! 1041: {
! 1042: addsrz(-1,p2,p2); p4 = divrr(p4,p2);
! 1043: }
! 1044: if (u)
! 1045: {
! 1046: setlg(pitemp,l+1); p1 = mulrr(pitemp,x);
! 1047: p4 = mulrr(mpsin(p1), p4); p4 = divrr(pitemp,p4);
! 1048: }
! 1049: affrr(p4,y); avma=av; return y;
! 1050: }
! 1051:
! 1052: static GEN
! 1053: cxgamma(GEN x, long prec)
! 1054: {
! 1055: GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;
! 1056: long l,l1,l2,u,i,k,e,s,n,p,av,av1;
! 1057: double alpha,beta,dk;
! 1058:
! 1059: l = precision(x); if (!l) l = prec;
! 1060: l2 = l+1; y=cgetg(3,t_COMPLEX);
! 1061: y[1]=lgetr(l); y[2]=lgetr(l); av=avma;
! 1062:
! 1063: p1=cgetg(3,t_COMPLEX); p1[1]=lgetr(l2); p1[2]=lgetr(l2);
! 1064: u = (typ(x[1])!=t_REAL && gsigne((GEN)x[1])<=0)? 1: (gexpo((GEN)x[1]) < -1);
! 1065: p2 = u? gsub(gun,x): x;
! 1066: gaffect(p2,p1);
! 1067:
! 1068: alpha=rtodbl(gabs(p1,DEFAULTPREC));
! 1069: beta = (bit_accuracy(l)>>1) * LOG2 / PI - alpha;
! 1070: if (beta>=0) n=(long)(1 + pariK2*beta); else n=0;
! 1071: if (n)
! 1072: {
! 1073: p = (long)(1+PI*(alpha+n));
! 1074: l2 += n>>TWOPOTBITS_IN_LONG;
! 1075: p2=cgetg(3,t_COMPLEX); p2[1]=lgetr(l2); p2[2]=lgetr(l2);
! 1076: addsrz(n,(GEN)p1[1],(GEN)p2[1]);
! 1077: affrr((GEN)p1[2],(GEN)p2[2]);
! 1078: }
! 1079: else
! 1080: {
! 1081: dk = pariK4*alpha/(l-2); beta = log(dk)/LOG2;
! 1082: if (beta>1.) beta += log(beta)/LOG2;
! 1083: p = (long)((bit_accuracy(l)>>1)/beta + 1);
! 1084: p2 = p1;
! 1085: }
! 1086: mpbern(p,l2); p3 = glog(p2,l2);
! 1087:
! 1088: p4=cgetg(3,t_COMPLEX);
! 1089: p4[1] = (long)realun(l2); setexpo(p4[1],-1);
! 1090: p4[2] = (long)rcopy((GEN)p2[2]);
! 1091: subrrz((GEN)p2[1],(GEN)p4[1],(GEN)p4[1]);
! 1092: gmulz(p4,p3,p4); gsubz(p4,p2,p4);
! 1093:
! 1094: pitemp = mppi(l2); setexpo(pitemp,2);
! 1095: p7 = mplog(pitemp); setexpo(pitemp,1);
! 1096: setexpo(p7,-1); addrrz((GEN)p4[1],p7, (GEN)p4[1]);
! 1097:
! 1098: gaffect(ginv(gsqr(p2)), p3); e=gexpo(p3);
! 1099:
! 1100: p5=cgetg(3,t_COMPLEX);
! 1101: p5[1]=lgetr(l2); setlg(p5[1],4);
! 1102: p5[2]=lgetr(l2); setlg(p5[2],4);
! 1103: p71=cgetr(l2); p7 = bern(p);
! 1104: if (bernzone[2]>4) { setlg(p71,4); affrr(p7,p71); p7=p71; }
! 1105: p7 = divrs(p7, 2*p*(2*p-1)); gaffect(p7,p5);
! 1106:
! 1107: s=0; l1=4; av1=avma;
! 1108: for (k=p-1; k>0; k--)
! 1109: {
! 1110: setlg(p3[1],l1); setlg(p3[2],l1);
! 1111: p6 = gmul(p3,p5); p7 = bern(k);
! 1112: if (bernzone[2]>l1) { setlg(p71,l1); affrr(p7,p71); p7=p71; }
! 1113: p7 = divrs(p7, (2*k)*(2*k-1));
! 1114: s -= e; l1 += (s>>TWOPOTBITS_IN_LONG); if (l1>l2) l1=l2;
! 1115: s &= (BITS_IN_LONG-1); p7 = addrr(p7, (GEN)p6[1]);
! 1116: setlg(p5[1],l1); affrr(p7, (GEN)p5[1]); p7 = (GEN)p6[2];
! 1117: setlg(p5[2],l1); affrr(p7, (GEN)p5[2]); avma=av1;
! 1118: }
! 1119: setlg(p5[1],l2); setlg(p5[2],l2);
! 1120: p6 = gdiv(p5,p2); setlg(p6[1],l2); setlg(p6[2],l2);
! 1121: p4 = gadd(p4,p6); setlg(p4[1],l2); setlg(p4[2],l2);
! 1122: p4 = gexp(p4,l2);
! 1123: for (i=1; i<=n; i++)
! 1124: {
! 1125: addsrz(-1,(GEN)p2[1],(GEN)p2[1]); p4 = gdiv(p4,p2);
! 1126: }
! 1127: if (u)
! 1128: {
! 1129: setlg(pitemp,l+1); p1 = gmul(pitemp,x);
! 1130: p4 = gmul(gsin(p1,l+1), p4); p4 = gdiv(pitemp,p4);
! 1131: }
! 1132: gaffect(p4,y); avma=av; return y;
! 1133: }
! 1134:
! 1135: GEN
! 1136: ggamma(GEN x, long prec)
! 1137: {
! 1138: switch(typ(x))
! 1139: {
! 1140: case t_INT:
! 1141: if (signe(x)<=0) err(gamer2);
! 1142: return transc(ggamma,x,prec);
! 1143:
! 1144: case t_REAL:
! 1145: return mpgamma(x);
! 1146:
! 1147: case t_COMPLEX:
! 1148: return gcmp0((GEN)x[2])? ggamma((GEN)x[1],prec): cxgamma(x,prec);
! 1149:
! 1150: case t_SER:
! 1151: return gexp(glngamma(x,prec),prec);
! 1152:
! 1153: case t_PADIC:
! 1154: err(impl,"p-adic gamma function");
! 1155: case t_INTMOD:
! 1156: err(typeer,"ggamma");
! 1157: }
! 1158: return transc(ggamma,x,prec);
! 1159: }
! 1160:
! 1161: void
! 1162: ggammaz(GEN x, GEN y)
! 1163: {
! 1164: long av=avma, prec = precision(y);
! 1165:
! 1166: if (!prec) err(infprecer,"ggammaz");
! 1167: gaffect(ggamma(x,prec),y); avma=av;
! 1168: }
! 1169:
! 1170: static GEN
! 1171: mplngamma(GEN x)
! 1172: {
! 1173: GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;
! 1174: long l,l1,l2,u,i,k,e,f,s,sx,n,p,av,av1;
! 1175: double alpha,beta,dk;
! 1176:
! 1177: sx=signe(x); if (!sx) err(talker,"zero argument in mplngamma");
! 1178: cgetg(3, t_COMPLEX); l=lg(x); y=cgetr(l); av=avma;
! 1179:
! 1180: l2=l+1; p1=cgetr(l2);
! 1181: u = (expo(x)<-1 || sx<0);
! 1182: if (!u) p2 = x;
! 1183: else
! 1184: {
! 1185: p2=gfrac(x); if (gcmp0(p2)) err(gamer2);
! 1186: p2 = subsr(1,x);
! 1187: }
! 1188: affrr(p2,p1);
! 1189: if (expo(p1)>1000)
! 1190: {
! 1191: n=0; beta = log(pariK4/(l-2))/LOG2+expo(p1);
! 1192: beta += log(beta)/LOG2;
! 1193: p = (long)((bit_accuracy(l)>>1)/beta + 1);
! 1194: p2 = p1;
! 1195: }
! 1196: else
! 1197: {
! 1198: alpha=rtodbl(p1); beta = ((bit_accuracy(l)>>1) * LOG2 / PI) - alpha;
! 1199: if (beta>=0) n=(long)(1 + pariK2*beta); else n=0;
! 1200: if (n)
! 1201: {
! 1202: p=(long)(1+PI*(alpha+n));
! 1203: l2 += n>>TWOPOTBITS_IN_LONG;
! 1204: p2 = cgetr(l2); addsrz(n,p1,p2);
! 1205: }
! 1206: else
! 1207: {
! 1208: dk = pariK4*alpha/(l-2); beta = log(dk)/LOG2;
! 1209: if (beta>1.) beta += (log(beta)/LOG2);
! 1210: p = (long)((bit_accuracy(l)>>1)/beta + 1);
! 1211: p2 = p1;
! 1212: }
! 1213: }
! 1214: mpbern(p,l2); p3=mplog(p2);
! 1215:
! 1216: p4=realun(l2); setexpo(p4,-1);
! 1217: p6 = subrr(p2,p4); p6 = mulrr(p6,p3);
! 1218: p6 = subrr(p6,p2);
! 1219: pitemp = mppi(l2); setexpo(pitemp,2);
! 1220: p7 = mplog(pitemp); setexpo(pitemp,1);
! 1221: setexpo(p7,-1); addrrz(p6,p7,p4);
! 1222:
! 1223: affrr(ginv(gsqr(p2)), p3); e=expo(p3);
! 1224:
! 1225: p5=cgetr(l2); setlg(p5,4);
! 1226: p71=cgetr(l2); p7 = bern(p);
! 1227: if (bernzone[2]>4) { setlg(p71,4); affrr(p7,p71); p7=p71; }
! 1228: p7 = divrs(p7, 2*p*(2*p-1)); affrr(p7,p5);
! 1229:
! 1230: s=0; l1=4; av1=avma;
! 1231: for (k=p-1; k>0; k--)
! 1232: {
! 1233: setlg(p3,l1); p6 = mulrr(p3,p5); p7 = bern(k);
! 1234: if (bernzone[2]>l1) { setlg(p71,l1); affrr(p7,p71); p7=p71; }
! 1235: p7 = divrs(p7, (2*k)*(2*k-1));
! 1236: s -= e; l1 += (s>>TWOPOTBITS_IN_LONG); if (l1>l2) l1=l2;
! 1237: s &= (BITS_IN_LONG-1); p7 = addrr(p7,p6);
! 1238: setlg(p5,l1); affrr(p7,p5); avma=av1;
! 1239: }
! 1240: setlg(p5,l2); p6 = divrr(p5,p2);
! 1241: p4 = addrr(p4,p6); setlg(p4,l2);
! 1242: if (n)
! 1243: {
! 1244: for (i=1; i<=n; i++)
! 1245: {
! 1246: addsrz(-1,p2,p2); p7 = (i==1)? rcopy(p2): mulrr(p7,p2);
! 1247: }
! 1248: f=signe(p7); if (f<0) setsigne(p7,1);
! 1249: subrrz(p4,mplog(p7),p4);
! 1250: }
! 1251: else f=1;
! 1252: if (u)
! 1253: {
! 1254: setlg(pitemp,l+1); p1 = mulrr(pitemp,x);
! 1255: p1 = divrr(pitemp,mpsin(p1));
! 1256: if (signe(p1) < 0) { setsigne(p1,1); f = -f; }
! 1257: p4 = subrr(mplog(p1),p4);
! 1258: }
! 1259: if (f<0) /* t_COMPLEX result */
! 1260: {
! 1261: p2 = y - 3;
! 1262: p2[1] = (long)y; affrr(p4,y); avma = av;
! 1263: p2[2] = (long)mppi(l); return p2;
! 1264: }
! 1265: /* t_REAL result */
! 1266: y[3] = y[0]; y += 3;
! 1267: affrr(p4,y); avma = (long)y; return y;
! 1268: }
! 1269:
! 1270: static GEN
! 1271: cxlngamma(GEN x, long prec)
! 1272: {
! 1273: GEN y,p1,p2,p3,p4,p5,p6,p7,p71,pitemp;
! 1274: long l,l1,l2,u,i,k,e,s,n,p,av,av1;
! 1275: double alpha,beta,dk;
! 1276:
! 1277: l = precision(x); if (!l) l = prec;
! 1278: l2 = l+1; y=cgetg(3,t_COMPLEX);
! 1279: y[1]=lgetr(l); y[2]=lgetr(l); av=avma;
! 1280:
! 1281: p1=cgetg(3,t_COMPLEX); p1[1]=lgetr(l2); p1[2]=lgetr(l2);
! 1282: u = (typ(x[1])!=t_REAL && gsigne((GEN)x[1])<=0)? 1: (gexpo((GEN)x[1]) < -1);
! 1283: if (u && (gcmp0((GEN)x[2]) || gexpo((GEN)x[2])>16)) u = 0;
! 1284: p2 = u? gsub(gun,x): x;
! 1285: gaffect(p2,p1);
! 1286:
! 1287: p2=gabs(p1,DEFAULTPREC);
! 1288: if (expo(p2)>1000)
! 1289: {
! 1290: n=0; beta = log(pariK4/(l-2)) / LOG2 + expo(p1);
! 1291: beta += log(beta)/LOG2;
! 1292: p = (long)((bit_accuracy(l)>>1)/beta + 1);
! 1293: p2 = p1;
! 1294: }
! 1295: else
! 1296: {
! 1297: alpha=rtodbl(p2); beta = ((bit_accuracy(l)>>1) * LOG2 / PI) - alpha;
! 1298: if (beta>=0) n=(long)(1+pariK2*beta); else n=0;
! 1299: if (n)
! 1300: {
! 1301: p = (long)(1+PI*(alpha+n));
! 1302: l2 += n>>TWOPOTBITS_IN_LONG;
! 1303: p2=cgetg(3,t_COMPLEX); p2[1]=lgetr(l2); p2[2]=lgetr(l2);
! 1304: addsrz(n,(GEN)p1[1],(GEN)p2[1]);
! 1305: affrr((GEN)p1[2],(GEN)p2[2]);
! 1306: }
! 1307: else
! 1308: {
! 1309: dk = pariK4*alpha/(l-2); beta = log(dk)/LOG2;
! 1310: if (beta>1.) beta += log(beta)/LOG2;
! 1311: p = (long)((bit_accuracy(l)>>1)/beta + 1);
! 1312: p2 = p1;
! 1313: }
! 1314: }
! 1315: mpbern(p,l2); p3 = glog(p2,l2);
! 1316:
! 1317: p4=cgetg(3,t_COMPLEX);
! 1318: p4[1] = (long)realun(l2); setexpo(p4[1],-1);
! 1319: subrrz((GEN)p2[1],(GEN)p4[1],(GEN)p4[1]);
! 1320: p4[2] = (long)rcopy((GEN)p2[2]);
! 1321: gmulz(p4,p3,p4); gsubz(p4,p2,p4);
! 1322:
! 1323: pitemp = mppi(l2); setexpo(pitemp,2);
! 1324: p7 = mplog(pitemp); setexpo(pitemp,1);
! 1325: setexpo(p7,-1); addrrz((GEN)p4[1],p7, (GEN)p4[1]);
! 1326:
! 1327: gaffect(ginv(gsqr(p2)),p3); e=gexpo(p3);
! 1328:
! 1329: p5=cgetg(3,t_COMPLEX);
! 1330: p5[1]=lgetr(l2); setlg(p5[1],4);
! 1331: p5[2]=lgetr(l2); setlg(p5[2],4);
! 1332: p71=cgetr(l2); p7 = bern(p);
! 1333: if (bernzone[2]>4) { setlg(p71,4); affrr(p7,p71); p7=p71; }
! 1334: p7 = divrs(p7, 2*p*(2*p-1)); gaffect(p7,p5);
! 1335:
! 1336: s=0; l1=4; av1=avma;
! 1337: for (k=p-1; k>0; k--)
! 1338: {
! 1339: setlg(p3[1],l1); setlg(p3[2],l1);
! 1340: p6 = gmul(p3,p5); p7 = bern(k);
! 1341: if (bernzone[2]>l1) { setlg(p71,l1); affrr(p7,p71); p7=p71; }
! 1342: p7 = divrs(p7, (2*k)*(2*k-1));
! 1343: s -= e; l1 += (s>>TWOPOTBITS_IN_LONG); if (l1>l2) l1=l2;
! 1344: s &= (BITS_IN_LONG-1); p7 = addrr(p7, (GEN)p6[1]);
! 1345: setlg(p5[1],l1); affrr(p7, (GEN)p5[1]); p7 = (GEN)p6[2];
! 1346: setlg(p5[2],l1); affrr(p7, (GEN)p5[2]); avma=av1;
! 1347: }
! 1348: setlg(p5[1],l2); setlg(p5[2],l2);
! 1349: p6 = gdiv(p5,p2); setlg(p6[1],l2); setlg(p6[2],l2);
! 1350: p4 = gadd(p4,p6); setlg(p4[1],l2); setlg(p4[2],l2);
! 1351:
! 1352: if (n)
! 1353: {
! 1354: p7 = cgetg(3,t_COMPLEX); p7[2] = p2[2];
! 1355: for (i=1; i<=n; i++)
! 1356: {
! 1357: addsrz(-1,(GEN)p2[1], (GEN)p2[1]);
! 1358: if (i==1) p7[1] = lrcopy((GEN)p2[1]); else p7 = gmul(p7,p2);
! 1359: }
! 1360: gsubz(p4,glog(p7,l+1), p4);
! 1361: }
! 1362: if (u)
! 1363: {
! 1364: setlg(pitemp,l+1); p1 = gmul(pitemp,x);
! 1365: p1 = gdiv(pitemp,gsin(p1,l+1));
! 1366: p4 = gsub(glog(p1,l+1),p4);
! 1367: }
! 1368: affrr((GEN)p4[1], (GEN)y[1]); setlg(p4[2],l+1);
! 1369:
! 1370: p1 = subrr(pitemp, (GEN)p4[2]); setexpo(pitemp,2);
! 1371: p1 = gfloor(divrr(p1,pitemp));
! 1372: p1 = addrr(mulir(p1,pitemp), (GEN)p4[2]);
! 1373: affrr(p1, (GEN)y[2]); avma=av; return y;
! 1374: }
! 1375:
! 1376: GEN
! 1377: glngamma(GEN x, long prec)
! 1378: {
! 1379: long i,av,tetpil,n;
! 1380: GEN y,p1;
! 1381:
! 1382: switch(typ(x))
! 1383: {
! 1384: case t_INT:
! 1385: if (signe(x)<=0) err(gamer2);
! 1386: return transc(glngamma,x,prec);
! 1387:
! 1388: case t_REAL:
! 1389: return mplngamma(x);
! 1390:
! 1391: case t_COMPLEX:
! 1392: return gcmp0((GEN)x[2])? glngamma((GEN)x[1],prec): cxlngamma(x,prec);
! 1393:
! 1394: case t_SER:
! 1395: av=avma; if (valp(x)) err(negexper,"glngamma");
! 1396: if (!gcmp1((GEN)x[2])) err(impl,"lngamma around a!=1");
! 1397: p1=gsubsg(1,x); n=(lg(x)-3)/valp(p1);
! 1398: y=ggrando(polx[varn(x)],lg(x)-2);
! 1399: for (i=n; i>=2; i--)
! 1400: y = gmul(p1,gadd(gdivgs(gzeta(stoi(i),prec),i),y));
! 1401: y = gadd(mpeuler(prec),y); tetpil=avma;
! 1402: return gerepile(av,tetpil,gmul(p1,y));
! 1403:
! 1404: case t_PADIC:
! 1405: err(impl,"p-adic lngamma function");
! 1406: case t_INTMOD:
! 1407: err(typeer,"glngamma");
! 1408: }
! 1409: return transc(glngamma,x,prec);
! 1410: }
! 1411:
! 1412: void
! 1413: glngammaz(GEN x, GEN y)
! 1414: {
! 1415: long av=avma, prec = precision(y);
! 1416:
! 1417: if (!prec) err(infprecer,"glngammaz");
! 1418: gaffect(glngamma(x,prec),y); avma=av;
! 1419: }
! 1420:
! 1421: /********************************************************************/
! 1422: /** **/
! 1423: /** FONCTION GAMMA DES DEMI-ENTIERS **/
! 1424: /** **/
! 1425: /********************************************************************/
! 1426:
! 1427: static GEN
! 1428: mpgamd(long x, long prec)
! 1429: {
! 1430: long i,j,a,l,av;
! 1431: GEN y,p1,p2;
! 1432:
! 1433: a = labs(x);
! 1434: l = prec+1+(a>>TWOPOTBITS_IN_LONG);
! 1435: if (l > LGBITS>>1) err(talker,"argument too large in ggamd");
! 1436: y=cgetr(prec); av=avma;
! 1437:
! 1438: p1 = mpsqrt(mppi(l));
! 1439: j = 1; p2 = realun(l);
! 1440: for (i=1; i<a; i++)
! 1441: {
! 1442: j += 2;
! 1443: p2 = mulsr(j,p2); setlg(p2,l);
! 1444: }
! 1445: if (x>=0) p1 = mulrr(p1,p2);
! 1446: else
! 1447: {
! 1448: p1 = divrr(p1,p2); if (a&1) setsigne(p1,-1);
! 1449: }
! 1450: setexpo(p1, expo(p1)-x);
! 1451: affrr(p1,y); avma=av; return y;
! 1452: }
! 1453:
! 1454: void
! 1455: mpgamdz(long s, GEN y)
! 1456: {
! 1457: long av=avma;
! 1458:
! 1459: affrr(mpgamd(s,lg(y)),y); avma=av;
! 1460: }
! 1461:
! 1462: GEN
! 1463: ggamd(GEN x, long prec)
! 1464: {
! 1465: long av,tetpil;
! 1466:
! 1467: switch(typ(x))
! 1468: {
! 1469: case t_INT:
! 1470: return mpgamd(itos(x),prec);
! 1471:
! 1472: case t_REAL: case t_FRAC: case t_FRACN: case t_COMPLEX: case t_QUAD:
! 1473: av=avma; x = gadd(x,ghalf); tetpil=avma;
! 1474: return gerepile(av,tetpil,ggamma(x,prec));
! 1475:
! 1476: case t_INTMOD: case t_PADIC:
! 1477: err(typeer,"ggamd");
! 1478: case t_SER:
! 1479: err(impl,"gamd of a power series");
! 1480: }
! 1481: return transc(ggamd,x,prec);
! 1482: }
! 1483:
! 1484: void
! 1485: ggamdz(GEN x, GEN y)
! 1486: {
! 1487: long av=avma, prec = precision(y);
! 1488:
! 1489: if (!prec) err(infprecer,"ggamdz");
! 1490: gaffect(ggamd(x,prec),y); avma=av;
! 1491: }
! 1492:
! 1493: /********************************************************************/
! 1494: /** **/
! 1495: /** FONCTION PSI **/
! 1496: /** **/
! 1497: /********************************************************************/
! 1498:
! 1499: #if 1
! 1500: static GEN
! 1501: mppsi(GEN z) /* version p=2 */
! 1502: {
! 1503: long l,n,k,x,xx,av,av1,tetpil;
! 1504: GEN zk,u,v,a,b;
! 1505:
! 1506: av=avma; l=lg(z);
! 1507: x=(long)(1 + (bit_accuracy(l)>>1)*LOG2 + 1.58*rtodbl(absr(z)));
! 1508: if (expo(z)>=15 || x>46340) err(impl,"psi(x) for x>=29000");
! 1509: xx=x*x; n=(long)(1+3.591*x);
! 1510: affsr(x,a=cgetr(l));
! 1511: a = mplog(a);
! 1512: gaffect(a,u=cgetr(l));
! 1513: gaffsg(1,b=cgetr(l));
! 1514: gaffsg(1,v=cgetr(l)); av1=avma;
! 1515: for (k=1; k<=n; k++)
! 1516: {
! 1517: zk = (k>1)? addsr(k-1,z): z;
! 1518: divrrz(mulsr(xx,b),gsqr(zk),b);
! 1519: divrrz(subrr(divrr(mulsr(xx,a),zk),b),zk,a);
! 1520: addrrz(u,a,u); addrrz(v,b,v); avma=av1;
! 1521: }
! 1522: tetpil=avma; return gerepile(av,tetpil,divrr(u,v));
! 1523: }
! 1524: #else
! 1525: static GEN /* by Manfred Radimersky */
! 1526: mppsi(GEN z)
! 1527: {
! 1528: long head, tail;
! 1529: long len, num, k;
! 1530: GEN a, b, f, s, x;
! 1531:
! 1532: head = avma;
! 1533: len = lg(z);
! 1534: num = bit_accuracy(len);
! 1535:
! 1536: if(signe(z) < 0) {
! 1537: x = subsr(1, z);
! 1538: s = mppsi(x);
! 1539: f = mulrr(mppi(len), z);
! 1540: mpsincos(f, &a, &b);
! 1541: f = divrr(b, a);
! 1542: a = mulrr(mppi(len), f);
! 1543: tail = avma;
! 1544: gerepile(head, tail, subrr(s, a));
! 1545: }
! 1546:
! 1547: a = cgetr(len);
! 1548: affsr(0, a);
! 1549: x = cgetr(len);
! 1550: affrr(z, x);
! 1551: tail = avma;
! 1552: while(cmprs(x, num >> 2) < 0) {
! 1553: gaddz(a, divsr(1, x), a);
! 1554: gaddgsz(x, 1, x);
! 1555: avma = tail;
! 1556: }
! 1557:
! 1558: s = mplog(x);
! 1559: gsubz(s, a, s);
! 1560: b = gmul2n(x, 1);
! 1561: gsubz(s, divsr(1, b), s);
! 1562:
! 1563: mpbern(num, len);
! 1564:
! 1565: affsr(-1, a);
! 1566: gdivgsz(a, 2, a);
! 1567: f = mulrr(x, x);
! 1568: f = divsr(1, f);
! 1569: k = 1;
! 1570: do {
! 1571: gmulz(a, f, a);
! 1572: gdivgsz(a, k, b);
! 1573: gmulz(b, bern(k), b);
! 1574: gaddz(s, b, s);
! 1575: k++;
! 1576: } while(expo(s) - expo(b) < num);
! 1577:
! 1578: tail = avma;
! 1579: return gerepile(head, tail, gcopy(s));
! 1580: }
! 1581: #endif
! 1582:
! 1583: #if 1
! 1584: static GEN
! 1585: cxpsi(GEN z, long prec) /* version p=2 */
! 1586: {
! 1587: long l,n,k,x,xx,av,av1,tetpil;
! 1588: GEN zk,u,v,a,b,p1;
! 1589:
! 1590: l = precision(z); if (!l) l = prec; av=avma;
! 1591: x=(long)(1 + (bit_accuracy(l)>>1)*LOG2 + 1.58*rtodbl(gabs(z,DEFAULTPREC)));
! 1592: xx=x*x; n=(long)(1+3.591*x);
! 1593: a=cgetg(3,t_COMPLEX); a[1]=lgetr(l); a[2]=lgetr(l); gaffsg(x,a);
! 1594: b=cgetg(3,t_COMPLEX); b[1]=lgetr(l); b[2]=lgetr(l); gaffsg(1,b);
! 1595: u=cgetg(3,t_COMPLEX); u[1]=lgetr(l); u[2]=lgetr(l);
! 1596: v=cgetg(3,t_COMPLEX); v[1]=lgetr(l); v[2]=lgetr(l); gaffsg(1,v);
! 1597: p1=glog(a,l); gaffect(p1,a); gaffect(p1,u); av1=avma;
! 1598: for (k=1; k<=n; k++)
! 1599: {
! 1600: zk=(k>1) ? gaddsg(k-1,z) : z;
! 1601: gdivz(gmulsg(xx,b),gsqr(zk),b);
! 1602: gdivz(gsub(gdiv(gmulsg(xx,a),zk),b),zk,a);
! 1603: gaddz(u,a,u); gaddz(v,b,v); avma=av1;
! 1604: }
! 1605: tetpil=avma; return gerepile(av,tetpil,gdiv(u,v));
! 1606: }
! 1607: #else
! 1608: GEN
! 1609: cxpsi(GEN z, long prec) /* by Manfred Radimersky */
! 1610: {
! 1611: long head, tail;
! 1612: long bord, nubi, k;
! 1613: GEN a, b, f, s, w, wn;
! 1614:
! 1615: head = avma;
! 1616: nubi = bit_accuracy(prec);
! 1617: bord = (nubi * nubi) >> 4;
! 1618:
! 1619: if(signe((GEN) z[1]) < 0) {
! 1620: w = gsubsg(1, z);
! 1621: s = cxpsi(w, prec);
! 1622: f = gmul(mppi(prec), z);
! 1623: gsincos(f, &a, &b, prec);
! 1624: f = gdiv(b, a);
! 1625: a = gmul(mppi(prec), f);
! 1626: tail = avma;
! 1627: gerepile(head, tail, gsub(s, a));
! 1628: }
! 1629:
! 1630: a = cgetg(3, t_COMPLEX);
! 1631: a[1] = (long) cgetr(prec);
! 1632: a[2] = (long) cgetr(prec);
! 1633: gaffsg(0, a);
! 1634: w = cgetg(3, t_COMPLEX);
! 1635: w[1] = (long) cgetr(prec);
! 1636: w[2] = (long) cgetr(prec);
! 1637: gaffect(z, w);
! 1638: wn = gnorm(w);
! 1639: tail = avma;
! 1640: while(cmprs(wn, bord) < 0) {
! 1641: gaddz(a, gdivsg(1, w), a);
! 1642: gaddgsz((GEN) w[1], 1, (GEN) w[1]);
! 1643: gaffect(gnorm(w), wn);
! 1644: avma = tail;
! 1645: }
! 1646:
! 1647: s = glog(w, prec);
! 1648: gsubz(s, a, s);
! 1649: b = gmul2n(w, 1);
! 1650: gsubz(s, gdivsg(1, b), s);
! 1651:
! 1652: mpbern(nubi, prec);
! 1653:
! 1654: gaffsg(-1, a);
! 1655: gdivgsz(a, 2, a);
! 1656: f = gmul(w, w);
! 1657: f = gdivsg(1, f);
! 1658: k = 1;
! 1659: tail = avma;
! 1660: do {
! 1661: gmulz(a, f, a);
! 1662: gdivgsz(a, k, b);
! 1663: gmulz(b, bern(k), b);
! 1664: gaddz(s, b, s);
! 1665: bord = expo(gnorm(s)) - expo(gnorm(b));
! 1666: k++;
! 1667: avma = tail;
! 1668: } while(bord < nubi << 1);
! 1669:
! 1670: tail = avma;
! 1671: return gerepile(head, tail, gcopy(s));
! 1672: }
! 1673: #endif
! 1674:
! 1675: GEN
! 1676: gpsi(GEN x, long prec)
! 1677: {
! 1678: switch(typ(x))
! 1679: {
! 1680: case t_REAL:
! 1681: return mppsi(x);
! 1682:
! 1683: case t_COMPLEX:
! 1684: return cxpsi(x,prec);
! 1685:
! 1686: case t_INTMOD: case t_PADIC:
! 1687: err(typeer,"gpsi");
! 1688: case t_SER:
! 1689: err(impl,"psi of power series");
! 1690: }
! 1691: return transc(gpsi,x,prec);
! 1692: }
! 1693:
! 1694: void
! 1695: gpsiz(GEN x, GEN y)
! 1696: {
! 1697: long av=avma, prec = precision(y);
! 1698:
! 1699: if (!prec) err(infprecer,"gpsiz");
! 1700: gaffect(gpsi(x,prec),y); avma=av;
! 1701: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to trans2.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