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