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