Annotation of OpenXM_contrib/pari-2.2/src/basemath/trans2.c, Revision 1.1.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