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