Annotation of OpenXM_contrib/pari-2.2/src/basemath/trans3.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: trans3.c,v 1.30 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 3) **/
20: /** **/
21: /********************************************************************/
22: #include "pari.h"
23:
24: GEN
25: cgetc(long l)
26: {
27: GEN u = cgetg(3,t_COMPLEX); u[1]=lgetr(l); u[2]=lgetr(l);
28: return u;
29: }
30:
31: static GEN
32: fix(GEN x, long l)
33: {
34: GEN y;
35: if (typ(x) == t_REAL) return x;
36: y = cgetr(l); gaffect(x, y); return y;
37: }
38:
39: /***********************************************************************/
40: /** **/
41: /** K BESSEL FUNCTION **/
42: /** **/
43: /***********************************************************************/
44:
45: GEN
46: kbessel0(GEN nu, GEN gx, long flag, long prec)
47: {
48: switch(flag)
49: {
50: case 0: return kbessel(nu,gx,prec);
51: case 1: return kbessel2(nu,gx,prec);
52: default: err(flagerr,"besselk");
53: }
54: return NULL; /* not reached */
55: }
56:
57: GEN
58: kbessel(GEN nu, GEN gx, long prec)
59: {
60: GEN x,y,p1,p2,zf,zz,s,t,q,r,u,v,e,f,c,d,ak,pitemp,nu2,w;
61: long l,lbin,av,av1,k,k2,l1,n2,n;
62:
63: if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec);
64: l = (typ(gx)==t_REAL)? lg(gx): prec;
65: y=cgetr(l); l1=l+1;
66: av=avma; x = fix(gx, l);
67: u=cgetr(l1); v=cgetr(l1); c=cgetr(l1); d=cgetr(l1);
68: e=cgetr(l1); f=cgetr(l1);
69: nu2=gmulgs(gsqr(nu),-4);
70: n = (long) (bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gnorm(nu)))) / 2;
71: n2=(n<<1); pitemp=mppi(l1);
72: /* this 10 should really be a 5, but then kbessel(15.99) enters oo loop */
73: lbin = 10 - bit_accuracy(l); av1=avma;
74: if (gcmpgs(x,n)<0)
75: {
76: zf=gsqrt(gdivgs(pitemp,n2),prec);
77: zz=cgetr(l1); gaffect(ginv(stoi(n2<<2)), zz);
78: s=gun; t=gzero;
79: for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
80: {
81: if (k2 & ~(MAXHALFULONG>>1))
82: p1 = gadd(mulss(k2,k2),nu2);
83: else
84: p1 = gaddsg(k2*k2,nu2);
85: ak=gdivgs(gmul(p1,zz),-k);
86: s=gadd(gun,gmul(ak,s));
87: t=gaddsg(k2,gmul(ak,t));
88: }
89: gmulz(s,zf,u); t=gmul2n(t,-1);
90: gdivgsz(gadd(gmul(t,zf),gmul(u,nu)),-n2,v); avma=av1;
91: affsr(n2,q=cgetr(l1));
92: r=gmul2n(x,1); av1=avma;
93: for(;;)
94: {
95: p1=divsr(5,q); if (expo(p1) >= -1) p1=ghalf;
96: p2=subsr(1,divrr(r,q));
97: if (gcmp(p1,p2)>0) p1=p2;
98: gnegz(p1,c); gaffsg(1,d); affrr(u,e); affrr(v,f);
99: for (k=1; ; k++)
100: {
101: w=gadd(gmul(gsubsg(k,ghalf),u), gmul(subrs(q,k),v));
102: w=gadd(w, gmul(nu,gsub(u,gmul2n(v,1))));
103: gdivgsz(gmul(q,v),k,u);
104: gdivgsz(w,k,v);
105: gmulz(d,c,d);
106: gaddz(e,gmul(d,u),e); p1=gmul(d,v);
107: gaddz(f,p1,f);
108: if (gexpo(p1)-gexpo(f) <= lbin) break;
109: avma=av1;
110: }
111: p1=u; u=e; e=p1;
112: p1=v; v=f; f=p1;
113: gmulz(q,gaddsg(1,c),q);
114: if (expo(subrr(q,r)) <= lbin) break;
115: }
116: gmulz(u,gpui(gdivgs(x,n),nu,prec),y);
117: }
118: else
119: {
120: p2=gmul2n(x,1);
121: zf=gsqrt(gdiv(pitemp,p2),prec);
122: zz=ginv(gmul2n(p2,2)); s=gun;
123: for (k=n2,k2=2*n2-1; k > 0; k--,k2-=2)
124: {
125: if (k2 & ~(MAXHALFULONG>>1))
126: p1=gadd(mulss(k2,k2),nu2);
127: else
128: p1=gaddsg(k2*k2,nu2);
129: ak=gdivgs(gmul(p1,zz),k);
130: s=gsub(gun,gmul(ak,s));
131: }
132: gmulz(s,zf,y);
133: }
134: gdivz(y,mpexp(x),y); avma=av; return y;
135: }
136:
137: /***********************************************************************/
138: /* **/
139: /** FONCTION U(a,b,z) GENERALE **/
140: /** ET CAS PARTICULIERS **/
141: /** **/
142: /***********************************************************************/
143: /* Assume gx > 0 and a,b complex */
144: /* This might one day be extended to handle complex gx */
145: /* see Temme, N. M. "The numerical computation of the confluent */
146: /* hypergeometric function U(a,b,z)" in Numer. Math. 41 (1983), */
147: /* no. 1, 63--82. */
148: GEN
149: hyperu(GEN a, GEN b, GEN gx, long prec)
150: {
151: GEN x,y,p1,p2,p3,zf,zz,s,t,q,r,u,v,e,f,c,d,w,a1,gn;
152: long l,lbin,av,av1,av2,k,l1,n,ex;
153:
154: if(gsigne(gx) <= 0) err(talker,"hyperu's third argument must be positive");
155: ex = (iscomplex(a) || iscomplex(b));
156:
157: l = (typ(gx)==t_REAL)? lg(gx): prec;
158: if (ex) y=cgetc(l); else y=cgetr(l);
159: l1=l+1; av=avma; x = fix(gx, l);
160: a1=gaddsg(1,gsub(a,b));
161: if (ex)
162: {
163: u=cgetc(l1); v=cgetc(l1); c=cgetc(l1);
164: d=cgetc(l1); e=cgetc(l1); f=cgetc(l1);
165: }
166: else
167: {
168: u=cgetr(l1); v=cgetr(l1); c=cgetr(l1);
169: d=cgetr(l1); e=cgetr(l1); f=cgetr(l1);
170: }
171: n=(long)(bit_accuracy(l)*LOG2 + PI*sqrt(gtodouble(gabs(gmul(a,a1),l1))));
172: lbin = 10-bit_accuracy(l); av1=avma;
173: if (gcmpgs(x,n)<0)
174: {
175: gn=stoi(n); zf=gpui(gn,gneg_i(a),l1);
176: zz=gdivsg(-1,gn); s=gun; t=gzero;
177: for (k=n-1; k>=0; k--)
178: {
179: p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
180: s=gaddsg(1,gmul(p1,s));
181: t=gadd(gaddsg(k,a),gmul(p1,t));
182: }
183: gmulz(s,zf,u); t=gmul(t,zz); gmulz(t,zf,v); avma=av1;
184: q=cgetr(l1); affsr(n,q); r=x; av1=avma;
185: do
186: {
187: p1=divsr(5,q); if (expo(p1)>= -1) p1=ghalf;
188: p2=subsr(1,divrr(r,q)); if (gcmp(p1,p2)>0) p1=p2;
189: gnegz(p1,c); avma=av1;
190: k=0; gaffsg(1,d);
191: gaffect(u,e); gaffect(v,f);
192: p3=gsub(q,b); av2=avma;
193: for(;;)
194: {
195: w=gadd(gmul(gaddsg(k,a),u),gmul(gaddsg(-k,p3),v));
196: k++; gdivgsz(gmul(q,v),k,u);
197: gdivgsz(w,k,v);
198: gmulz(d,c,d);
199: gaddz(e,gmul(d,u),e); p1=gmul(d,v);
200: gaddz(f,p1,f);
201: if (gexpo(p1)-gexpo(f) <= lbin) break;
202: avma=av2;
203: }
204: p1=u; u=e; e=p1;
205: p1=v; v=f; f=p1;
206: gmulz(q,gadd(gun,c),q);
207: p1=subrr(q,r); ex=expo(p1); avma=av1;
208: }
209: while (ex>lbin);
210: }
211: else
212: {
213: zf=gpui(x,gneg_i(a),l1);
214: zz=gdivsg(-1,x); s=gun;
215: for (k=n-1; k>=0; k--)
216: {
217: p1=gdivgs(gmul(gmul(gaddgs(a,k),gaddgs(a1,k)),zz),k+1);
218: s=gadd(gun,gmul(p1,s));
219: }
220: u = gmul(s,zf);
221: }
222: gaffect(u,y); avma=av; return y;
223: }
224:
225: GEN
226: kbessel2(GEN nu, GEN x, long prec)
227: {
228: long av = avma,tetpil;
229: GEN p1,p2,x2,a,pitemp;
230:
231: if (typ(x)==t_REAL) prec = lg(x);
232: x2=gshift(x,1); pitemp=mppi(prec);
233: if (gcmp0(gimag(nu))) a=cgetr(prec); else a=cgetc(prec);
234: gaddz(gun,gshift(nu,1),a);
235: p1=hyperu(gshift(a,-1),a,x2,prec);
236: p1=gmul(gmul(p1,gpui(x2,nu,prec)),mpsqrt(pitemp));
237: p2=gexp(x,prec); tetpil=avma;
238: return gerepile(av,tetpil,gdiv(p1,p2));
239: }
240:
241: GEN
242: incgam(GEN s, GEN x, long prec)
243: {
244: GEN p1,z = cgetr(prec);
245: long av = avma;
246:
247: if (gcmp0(x)) return ggamma(s,prec);
248: if (typ(x)!=t_REAL) { gaffect(x,z); x=z; }
249: if (gcmp(subrs(x,1),s) > 0 || gsigne(greal(s)) <= 0)
250: p1 = incgam2(s,x,prec);
251: else
252: p1 = gsub(ggamma(s,prec), incgam3(s,x,prec));
253: affrr(p1,z); avma = av; return z;
254: }
255:
256: /* = incgam2(0, x, prec). typ(x) = t_REAL. Optimized for eint1 */
257: static GEN
258: incgam2_0(GEN x)
259: {
260: long l,n,i;
261: GEN p1;
262: double m,mx;
263:
264: l = lg(x); mx = rtodbl(x);
265: m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
266: p1 = divsr(-n, addsr(n<<1,x));
267: for (i=n-1; i >= 1; i--)
268: p1 = divsr(-i, addrr(addsr(i<<1,x), mulsr(i,p1)));
269: return mulrr(divrr(mpexp(negr(x)), x), addrr(realun(l),p1));
270: }
271:
272: GEN
273: incgam1(GEN a, GEN x, long prec)
274: {
275: GEN p2,p3,y, z = cgetr(prec);
276: long av=avma, av1, l,n,i;
277: double m,mx;
278:
279: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
280: l=lg(x); mx=rtodbl(x);
281: m=(long) bit_accuracy(l)*LOG2; n=(long)(m/(log(m)-(1+log(mx))));
282: p2 = cgetr(l); affrr(addir(gun,gsub(x,a)), p2);
283: p3 = subrs(p2, n+1); av1 = avma;
284: for (i=n; i>=1; i--)
285: {
286: affrr(addrr(subrs(p2,i), divrr(mulsr(i,x),p3)), p3);
287: avma = av1;
288: }
289: y = mulrr(mpexp(negr(x)),gpui(x,a,prec));
290: affrr(divrr(y,p3), z);
291: avma = av; return z;
292: }
293:
294: /* assume x != 0 */
295: GEN
296: incgam2(GEN a, GEN x, long prec)
297: {
298: GEN b,p1,p2,p3,y, z = cgetr(prec);
299: long av = avma,av1,l,n,i;
300: double m,mx;
301:
302: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
303: if (gcmp0(a)) { affrr(incgam2_0(x), z); avma = av; return z; }
304: l=lg(x); mx=rtodbl(x);
305: m = (bit_accuracy(l)*LOG2 + mx)/4; n=(long)(1+m*m/mx);
306: i = typ(a);
307: if (i == t_REAL) b = addsr(-1,a);
308: else
309: { /* keep b integral : final powering more efficient */
310: gaffect(a,p1=cgetr(prec));
311: b = (i == t_INT)? addsi(-1,a): addsr(-1,p1);
312: a = p1;
313: }
314: p2 = cgetr(l); gaffect(subrr(x,a),p2);
315: p3 = divrr(addsr(-n,a), addrs(p2,n<<1));
316: av1 = avma;
317: for (i=n-1; i>=1; i--)
318: {
319: affrr(divrr(addsr(-i,a), addrr(addrs(p2,i<<1),mulsr(i,p3))), p3);
320: avma = av1;
321: }
322: y = gmul(mpexp(negr(x)), gpui(x,b,prec));
323: affrr(mulrr(y, addsr(1,p3)), z);
324: avma = av; return z;
325: }
326:
327: GEN
328: incgam3(GEN s, GEN x, long prec)
329: {
330: GEN b,p1,p2,p3,y, z = cgetr(prec);
331: long av = avma, av1,l,n,i;
332:
333: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
334: l=lg(x); n = -bit_accuracy(l)-1;
335: p3 = realun(l);
336: p2 = rcopy(p3);
337: i = typ(s);
338: if (i == t_REAL) b = s;
339: else
340: {
341: gaffect(s,p1=cgetr(prec));
342: b = (i == t_INT)? s: p1;
343: s = p1;
344: }
345: if (signe(s) <= 0)
346: {
347: p1 = gcvtoi(s, &i);
348: if (i < 5 - bit_accuracy(prec))
349: err(talker,"negative argument too close to an integer in incgamc");
350: }
351: av1 = avma;
352: for (i=1; expo(p2) >= n; i++)
353: {
354: affrr(divrr(mulrr(x,p2), addsr(i,s)), p2);
355: affrr(addrr(p2,p3), p3); avma = av1;
356: }
357: y = gdiv(mulrr(mpexp(negr(x)), gpui(x,b,prec)), s);
358: affrr(mulrr(y,p3), z);
359: avma = av; return z;
360: }
361:
362: /* Assumes that g=gamma(a,prec). Unchecked */
363: GEN
364: incgam4(GEN a, GEN x, GEN g, long prec)
365: {
366: GEN p1, z = cgetr(prec);
367: long av = avma;
368:
369: if (typ(x) != t_REAL) { gaffect(x,z); x=z; }
370: if (gcmp(subrs(x,1),a) > 0 || gsigne(greal(a)) <= 0)
371: p1 = incgam2(a,x,prec);
372: else
373: p1 = gsub(g, incgam3(a,x,prec));
374: affrr(p1, z);
375: avma = av; return z;
376: }
377:
378: GEN
379: incgam0(GEN a, GEN x, GEN z,long prec)
380: {
381: return z? incgam4(a,x,z,prec): incgam(a,x,prec);
382: }
383:
384: GEN
385: eint1(GEN x, long prec)
386: {
387: long av = avma,tetpil,l,n,i;
388: GEN p1,p2,p3,p4,run,y;
389:
390: if (typ(x) != t_REAL) { gaffect(x,p1=cgetr(prec)); x = p1;}
391: if (signe(x) >= 0)
392: {
393: if (expo(x) >= 4)
394: return gerepileupto(av, incgam2_0(x));
395:
396: l = lg(x);
397: n = -bit_accuracy(l)-1;
398:
399: run = realun(l);
400: p4 = p3 = p2 = p1 = run;
401: for (i=2; expo(p2)>=n; i++)
402: {
403: p1 = addrr(p1,divrs(run,i)); /* p1 = sum_{i=1} 1/i */
404: p4 = divrs(mulrr(x,p4),i); /* p4 = sum_{i=1} x^(i-1)/i */
405: p2 = mulrr(p4,p1);
406: p3 = addrr(p2,p3);
407: }
408: p3 = mulrr(x,mulrr(mpexp(negr(x)),p3));
409: p1 = addrr(mplog(x), mpeuler(l));
410: return gerepileupto(av, subrr(p3,p1));
411: }
412: else
413: { /* written by Manfred Radimersky */
414: l = lg(x);
415: n = bit_accuracy(l);
416: /* IS: line split to avoid a Workshop cc-5.0 bug (Sun BugID #4254959) */
417: n = 3 * n / 4;
418: y = negr(x);
419: if(gcmpgs(y, n) < 0) {
420: p1 = p2 = p3 = y;
421: p4 = gzero;
422: i = 2;
423: while(gcmp(p3, p4)) {
424: p4 = p3;
425: p1 = gmul(p1, gdivgs(y, i));
426: p2 = gdivgs(p1, i);
427: p3 = gadd(p3, p2);
428: i++;
429: }
430: p1 = gadd(mplog(y), mpeuler(l));
431: y = gadd(p3, p1);
432: } else {
433: p1 = gdivsg(1, y);
434: p2 = realun(l);
435: p3 = p2;
436: p4 = gzero;
437: i = 1;
438: while(gcmp(p3, p4)) {
439: p4 = p3;
440: p2 = gmulgs(p2, i);
441: p2 = gmul(p2, p1);
442: p3 = gadd(p3, p2);
443: i++;
444: }
445: p1 = gdiv(mpexp(y), y);
446: y = gmul(p3, p1);
447: }
448: tetpil = avma;
449: y = gerepile(av, tetpil, negr(y));
450: }
451: return y;
452: }
453:
454: GEN
455: veceint1(GEN C, GEN nmax, long prec)
456: {
457: long av,av1,k,n,nstop,i,cd,nmin,G,a, chkpoint;
458: GEN y,e1,e2,F0,F,M2,f,den,minvn,mcn,p1,vdiff,ap,unr,zeror,deninit;
459:
460: if (!nmax) return eint1(C,prec);
461:
462: if (signe(nmax)<=0) return cgetg(1,t_VEC);
463: if (DEBUGLEVEL>1) fprintferr("Entering veceint1:\n");
464: if (typ(C) != t_REAL || lg(C) > prec)
465: { y = cgetr(prec); gaffect(C,y); C = y; }
466: if (signe(C) <= 0) err(talker,"negative or zero constant in veceint1");
467:
468: G = -bit_accuracy(prec);
469: n=itos(nmax); y=cgetg(n+1,t_VEC);
470: for(i=1; i<=n; i++) y[i]=lgetr(prec);
471: av=avma;
472:
473: nstop = itos(gceil(divsr(4,C)));
474: if (nstop<1) nstop=1;
475: if (nstop>n) nstop=n;
476: nmin=n-10; if (nmin<nstop) nmin=nstop;
477: if(DEBUGLEVEL>1) fprintferr("nstop = %ld\n",nstop);
478:
479: e1=mpexp(negr(mulsr(n,C)));
480: e2=mpexp(mulsr(10,C));
481: unr = realun(prec);
482: zeror=realzero(prec); deninit=negr(unr);
483: f=cgetg(3,t_COL); M2=cgetg(3,t_VEC); av1=avma;
484:
485: F0=(GEN)y[n]; chkpoint = n;
486: affrr(eint1(mulsr(n,C),prec), F0);
487: do
488: {
489: if (DEBUGLEVEL>1 && n < chkpoint)
490: { fprintferr("%ld ",n) ; chkpoint -= (itos(nmax) / 20); }
491: minvn=divrs(unr,-n); mcn=mulrr(C,minvn);
492: M2[1] = (long)zeror; M2[2] = lsubrr(minvn,C);
493: f[1]=(long)zeror; f[2]=ldivrs(e1,-n);
494: affrr(mulrr(e1,e2), e1);
495: vdiff=cgetg(2,t_VEC); vdiff[1]=f[2];
496: for (cd=a=1,n--; n>=nmin; n--,a++)
497: {
498: F = F0;
499: ap = stoi(a); den = deninit;
500: for (k=1;;)
501: {
502: if (k>cd)
503: {
504: cd++; p1 = (GEN)f[2];
505: f[2] = lmul(M2,f);
506: f[1] = (long)p1;
507: M2[1] = laddrr((GEN)M2[1],mcn);
508: M2[2] = laddrr((GEN)M2[2],minvn);
509: vdiff = concatsp(vdiff,(GEN)f[2]);
510: }
511: p1 = mulrr(mulri(den,ap),(GEN)vdiff[k]);
512: if (expo(p1) < G) { affrr(F,(GEN)y[n]); break; }
513: F = addrr(F,p1); ap = mulis(ap,a);
514: k++; den = divrs(den,-k);
515: }
516: }
517: avma=av1; n++; F0=(GEN)y[n]; nmin -= 10;
518: if (nmin < nstop) nmin=nstop;
519: }
520: while(n>nstop);
521: for(i=1; i<=nstop; i++)
522: affrr(eint1(mulsr(i,C),prec), (GEN)y[i]);
523: if (DEBUGLEVEL>1) fprintferr("\n");
524: avma=av; return y;
525: }
526:
527: GEN
528: gerfc(GEN x, long prec)
529: {
530: long av=avma;
531: GEN p1,p2;
532:
533: if (typ(x)!=t_REAL) { p1=cgetr(prec); gaffect(x,p1); x=p1; }
534: av = avma; p1 = incgam(ghalf,gsqr(x),prec);
535: p2 = mpsqrt(mppi(lg(x)));
536: p1 = divrr(p1,p2);
537: if (signe(x) < 0) p1 = subsr(2,p1);
538: return gerepileupto(av,p1);
539: }
540:
541: /***********************************************************************/
542: /** **/
543: /** FONCTION ZETA DE RIEMANN **/
544: /** **/
545: /***********************************************************************/
546:
547: #if 0
548: static GEN
549: czeta(GEN s, long prec)
550: {
551: long av,n,p,n1,flag1,flag2,i,i2;
552: double st,sp,sn,ssig,ns,alpha,beta,maxbeta,xinf;
553: GEN y,z,res,sig,ms,p1,p2,p3,p31,pitemp;
554:
555: i=precision(s); if (i) prec = i;
556: if (typ(s)==t_COMPLEX)
557: {
558: res = cgetc(prec); av=avma;
559: p1 = cgetc(prec+1);
560: gaffect(s,p1); s=p1; sig=(GEN)s[1];
561: }
562: else
563: {
564: res = cgetr(prec); av=avma;
565: p1 = cgetr(prec+1); affrr(s,p1); sig = s = p1;
566: }
567:
568: if (signe(sig)>0 && expo(sig)>-2) flag1 = 0;
569: else
570: {
571: if (gcmp0(gimag(s)) && gcmp0(gfrac(gmul2n(sig,-1))))
572: {
573: if (gcmp0(sig)) gaffect(gneg_i(ghalf),res); else gaffsg(0,res);
574: avma=av; return res;
575: }
576: flag1 = 1; s = gsub(gun,s); sig = greal(s);
577: }
578: ssig=rtodbl(sig); st=fabs(rtodbl(gimag(s))); maxbeta = pow(3.0,-2.5);
579: if (st)
580: {
581: ns = ssig*ssig + st*st;
582: alpha=pariC2*(prec-2)-0.39-0.5*(ssig-1.0)*log(ns)-log(ssig)+ssig*2*pariC1;
583: beta=(alpha+ssig)/st-atan(ssig/st);
584: if (beta<=0)
585: {
586: if (ssig>=1.0)
587: {
588: p=0; sn=sqrt(ns);
589: n=(long)(ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig)));
590: }
591: else
592: {
593: p=1; sn=ssig+1; n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
594: }
595: }
596: else
597: {
598: if (beta<maxbeta) xinf=beta+pow(3*beta,1.0/3.0);
599: else
600: {
601: double eps=0.0087, x00 = beta+PI/2.0, y00,x11;
602: for(;;)
603: {
604: y00=x00*x00; x11=(beta+atan(x00))*(1+y00)/y00-1/x00;
605: if (x00-x11 < eps) break;
606: x00 = x11;
607: }
608: xinf=x11;
609: }
610: sp=1.0-ssig+st*xinf;
611: if (sp>0)
612: {
613: p=(long)ceil(sp/2.0); sn=ssig+2*p-1;
614: n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
615: }
616: else
617: {
618: p=0; sn=sqrt(ns);
619: n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
620: }
621: }
622: }
623: else
624: {
625: beta=pariC2*(prec-2)+0.61+ssig*2*pariC1-ssig*log(ssig);
626: if (beta>0)
627: {
628: p=(long)ceil(beta/2.0); sn=ssig+2*p-1;
629: n=(long)ceil(sqrt(sn*sn+st*st)/(2*PI));
630: }
631: else
632: {
633: p=0; sn=sqrt(ssig*ssig+st*st);
634: n=(long)ceil(exp(pariC2*(prec-2)/ssig)*pow(sn/(2*ssig),1.0/ssig));
635: }
636: }
637: if (n < 46340) n1=n*n; else n1=0;
638: y=gun; ms=gneg_i(s); p1=cgetr(prec+1); p2=gun;
639: for (i=2; i<=n; i++)
640: {
641: affsr(i,p1); p2 = gexp(gmul(ms,mplog(p1)), prec+1);
642: y = gadd(y,p2);
643: }
644: flag2 = (2*p < 46340);
645: mpbern(p,prec+1); p31=cgetr(prec+1); z=gzero;
646: for (i=p; i>=1; i--)
647: {
648: i2=i<<1;
649: p1=gmul(gaddsg(i2-1,s),gaddsg(i2,s));
650: p1=flag2? gdivgs(p1,i2*(i2+1)): gdivgs(gdivgs(p1,i2),i2+1);
651: p1=n1? gdivgs(p1,n1): gdivgs(gdivgs(p1,n),n);
652: p3 = bern(i);
653: if (bernzone[2]>prec+1) { affrr(p3,p31); p3=p31; }
654: z=gadd(divrs(p3,i),gmul(p1,z));
655: }
656: p1=gsub(gdivsg(n,gsubgs(s,1)),ghalf);
657: z=gmul(gadd(p1,gmul(s,gdivgs(z,n<<1))),p2);
658: y = gadd(y,z);
659: if (flag1)
660: {
661: pitemp=mppi(prec+1); setexpo(pitemp,2);
662: y=gmul(gmul(y,ggamma(s,prec+1)),gpui(pitemp,ms,prec+1));
663: setexpo(pitemp,0);
664: y=gmul2n(gmul(y,gcos(gmul(pitemp,s),prec+1)),1);
665: }
666: gaffect(y,res); avma = av; return res;
667: }
668:
669: /* y = binomial(n,k-2). Return binomial(n,k) */
670: static GEN
671: next_bin(GEN y, long n, long k)
672: {
673: y = divrs(mulrs(y, n-k+2), k-1);
674: return divrs(mulrs(y, n-k+1), k);
675: }
676:
677: static GEN
678: izeta(long k, long prec)
679: {
680: long av=avma,av2,kk,n,li, limit;
681: GEN y,p1,pi2,qn,z,q,binom;
682:
683: /* treat trivial cases */
684: if (!k) return gneg(ghalf);
685: if (k < 0)
686: {
687: if ((k&1) == 0) return gzero;
688: y = bernreal(1-k,prec);
689: return gerepileuptoleaf(av, divrs(y,k-1));
690: }
691: if (k > bit_accuracy(prec)+1) return realun(prec);
692: pi2 = mppi(prec); setexpo(pi2,2); /* 2Pi */
693: if ((k&1) == 0)
694: {
695: p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec)));
696: y = divrr(p1, mpfactr(k,prec)); setexpo(y,expo(y)-1);
697: return gerepileuptoleaf(av, y);
698: }
699: /* k > 1 odd */
700: binom = realun(prec+1);
701: q = mpexp(pi2); kk = k+1; /* >= 4 */
702: li = -(1+bit_accuracy(prec));
703: y = NULL; /* gcc -Wall */
704: if ((k&3)==3)
705: {
706: for (n=0; n <= kk>>1; n+=2)
707: {
708: p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
709: if (n) { binom = next_bin(binom,kk,n); setlg(binom,prec+1); }
710: p1 = mulrr(binom,p1);
711: if (n == kk>>1) setexpo(p1, expo(p1)-1);
712: if ((n>>1)&1) setsigne(p1,-signe(p1));
713: y = n? addrr(y,p1): p1;
714: }
715: y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y);
716:
717: av2 = avma; limit = stack_lim(av2,1);
718: qn = gsqr(q); z = divsr(1,addrs(q,-1));
719: for (n=2; ; n++)
720: {
721: p1 = divsr(1,mulir(gpuigs(stoi(n),k),addrs(qn,-1)));
722:
723: z = addrr(z,p1); if (expo(p1)< li) break;
724: qn = mulrr(qn,q);
725: if (low_stack(limit,stack_lim(av2,1)))
726: {
727: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;
728: if (DEBUGMEM>1) err(warnmem,"izeta");
729: gerepilemany(av2,gptr,2);
730: }
731: }
732: setexpo(z,expo(z)+1);
733: y = addrr(y,z); setsigne(y,-signe(y));
734: }
735: else
736: {
737: GEN p2 = divrs(pi2, k-1);
738: for (n=0; n <= k>>1; n+=2)
739: {
740: p1 = mulrr(bernreal(kk-n,prec),bernreal(n,prec));
741: if (n) binom = next_bin(binom,kk,n);
742: p1 = mulrr(binom,p1);
743: p1 = mulsr(kk-(n<<1),p1);
744: if ((n>>1)&1) setsigne(p1,-signe(p1));
745: y = n? addrr(y,p1): p1;
746: }
747: y = mulrr(divrr(gpuigs(pi2,k),mpfactr(kk,prec)),y);
748: y = divrs(y,k-1);
749: av2 = avma; limit = stack_lim(av2,1);
750: qn = q; z=gzero;
751: for (n=1; ; n++)
752: {
753: p1=mulir(gpuigs(stoi(n),k),gsqr(addrs(qn,-1)));
754: p1=divrr(addrs(mulrr(qn,addsr(1,mulsr(n<<1,p2))),-1),p1);
755:
756: z=addrr(z,p1); if (expo(p1) < li) break;
757: qn=mulrr(qn,q);
758: if (low_stack(limit,stack_lim(av2,1)))
759: {
760: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&qn;
761: if (DEBUGMEM>1) err(warnmem,"izeta");
762: gerepilemany(av2,gptr,2);
763: }
764: }
765: setexpo(z,expo(z)+1);
766: y = subrr(y,z);
767: }
768: return gerepileuptoleaf(av, y);
769: }
770: #else
771:
772: /* return x^n, assume n > 0 */
773: static long
774: pows(long x, long n)
775: {
776: long i, y = x;
777: for (i=1; i<n; i++) y *= x;
778: return y;
779: }
780:
781: /* return n^-s, n > 1 odd. tab[q] := q^-s, q prime power */
782: static GEN
783: n_s(long n, GEN *tab)
784: {
785: byteptr d = diffptr + 2;
786: GEN x = NULL;
787: long p, e;
788:
789: for (p = 3; n > 1; p += *d++)
790: {
791: e = svaluation(n, p, &n);
792: if (e)
793: {
794: GEN y = tab[pows(p,e)];
795: if (!x) x = y; else x = gmul(x,y);
796: }
797: }
798: return x;
799: }
800:
801: GEN czeta(GEN s0, long prec);
802:
803: /* assume k != 1 */
804: static GEN
805: izeta(long k, long prec)
806: {
807: long av = avma;
808: GEN y,p1,pi2;
809:
810: /* treat trivial cases */
811: if (!k) { y = realun(prec); setexpo(y,-1); setsigne(y,-1); return y; }
812: if (k < 0)
813: {
814: if ((k&1) == 0) return gzero;
815: y = bernreal(1-k,prec);
816: return gerepileuptoleaf(av, divrs(y,k-1));
817: }
818: if (k > bit_accuracy(prec)+1) return realun(prec);
819: if ((k&1) == 0)
820: {
821: pi2 = mppi(prec); setexpo(pi2,2); /* 2Pi */
822: p1 = mulrr(gpuigs(pi2,k),absr(bernreal(k,prec)));
823: y = divrr(p1, mpfactr(k,prec)); setexpo(y,expo(y)-1);
824: return gerepileuptoleaf(av, y);
825: }
826: /* k > 1 odd */
827: return czeta(stoi(k), prec);
828: }
829:
830: extern GEN rpowsi(ulong a, GEN n, long prec);
831: extern GEN divrs2_safe(GEN x, long i);
832: extern void dcxlog(double s, double t, double *a, double *b);
833: extern double dnorm(double s, double t);
834: extern GEN trans_fix_arg(long *prec, GEN *s0, GEN *sig, long *av, GEN *res);
835:
836: /* s0 a t_INT, t_REAL or t_COMPLEX.
837: * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd)
838: * */
839: GEN
840: czeta(GEN s0, long prec)
841: {
842: GEN s, u, a, y, res, tes, sig, invn2, p1, unr;
843: GEN sim, ms, s1, s2, s3, s4, s5, *tab, tabn;
844: long p, i, sqn, nn, lim, lim2, ct, av, av2 = avma, avlim;
845: int funeq = 0;
846: byteptr d;
847:
848: if (DEBUGLEVEL) timer2();
849: s = trans_fix_arg(&prec,&s0,&sig,&av,&res);
850: if (gcmp0(s)) { y = gneg(ghalf); goto END; }
851: if (signe(sig) <= 0 || expo(sig) < -1)
852: { /* s <--> 1-s */
853: if (typ(s0) == t_INT)
854: {
855: p = itos(s0); avma = av2;
856: return izeta(p,prec);
857: }
858: funeq = 1; s = gsub(gun, s); sig = greal(s);
859: }
860: if (gcmp(sig, stoi(bit_accuracy(prec) + 1)) > 0) { y = gun; goto END; }
861:
862: { /* find "optimal" parameters [lim, nn] */
863: double ssig = rtodbl(sig);
864: double st = rtodbl(gimag(s));
865: double ns = dnorm(ssig,st), l,l2;
866: long la = 1;
867:
868: if (typ(s0) == t_INT)
869: {
870: long ss = itos(s0);
871: switch(ss)
872: { /* should depend on prec ? */
873: case 3: la = 6; break;
874: default: la = 3; break;
875: }
876: }
877: if (dnorm(ssig-1,st) < 0.1) /* |s - 1| < 0.1 */
878: l2 = -(ssig - 0.5);
879: else
880: { /* l2 = Re( (s - 1/2) log (s-1) ) */
881: double rlog, ilog; /* log(s-1) */
882: dcxlog(ssig-1,st, &rlog,&ilog);
883: l2 = (ssig - 0.5)*rlog - st*ilog;
884: }
885: l = (pariC2*(prec-2) - l2 + ssig*2*pariC1) / (2. * (1.+ log((double)la)));
886: l2 = sqrt(ns)/2;
887: if (l < l2) l = l2;
888: lim = (long) ceil(l); if (lim < 2) lim = 2;
889: l2 = (lim+ssig/2.-.25);
890: nn = 1 + (long)ceil( sqrt(l2*l2 + st*st/4) * la / PI );
891: if (DEBUGLEVEL) fprintferr("lim, nn: [%ld, %ld]\n",lim,nn);
892: if ((ulong)nn >= maxprime()) err(primer1);
893: }
894: prec++; unr = realun(prec); /* one extra word of precision */
895:
896: tab = (GEN*)cgetg(nn, t_VEC); /* table of q^(-s), q = p^e */
897: d = diffptr + 1;
898: if (typ(s0) == t_INT)
899: { /* no explog for 1/p^s */
900: for (p=2; p < nn; p += *d++)
901: tab[p] = divrr(unr, rpowsi(p, s0, prec));
902: a = divrr(unr, rpowsi(nn, s0, prec));
903: }
904: else
905: { /* general case */
906: ms = gneg(s); p1 = cgetr(prec);
907: for (p=2; p < nn; p += *d++)
908: {
909: affsr(p, p1);
910: tab[p] = gexp(gmul(ms, mplog(p1)), prec);
911: }
912: affsr(nn,p1);
913: a = gexp(gmul(ms, mplog(p1)), prec);
914: }
915: sqn = (long)sqrt(nn-1.);
916: d = diffptr + 2; /* fill in odd prime powers */
917: for (p=3; p <= sqn; p += *d++)
918: {
919: ulong oldq = p, q = p*p;
920: while (q<(ulong)nn) { tab[q] = gmul(tab[p], tab[oldq]); oldq = q; q *= p; }
921: }
922: if (DEBUGLEVEL) msgtimer("tab[q^-s] from 1 to N-1");
923:
924: tabn = cgetg(nn, t_VECSMALL); ct = 0;
925: for (i = nn-1; i; i>>=1) tabn[++ct] = (i-1)>>1;
926: sim = y = unr;
927: for (i=ct; i > 1; i--)
928: {
929: long j, av2 = avma;
930: for (j=tabn[i]+1; j<=tabn[i-1]; j++)
931: sim = gadd(sim, n_s(2*j+1, tab));
932: sim = gerepileupto(av2, sim);
933: y = gadd(sim, gmul(tab[2],y));
934: }
935: y = gadd(y, gmul2n(a,-1));
936: if (DEBUGLEVEL) msgtimer("sum from 1 to N-1");
937:
938: invn2 = divrs(unr, nn*nn); lim2 = lim<<1;
939: tes = bernreal(lim2, prec);
940: if (typ(s0) == t_INT)
941: {
942: av2 = avma; avlim = stack_lim(av2,3);
943: for (i=lim2-2; i>=2; i-=2)
944: { /* using single prec (when (s0 + i) < 2^31) not faster (even at \p28) */
945: u = mulri(mulrr(tes,invn2), mulii(addsi(i,s0), addsi(i-1,s0)));
946: tes = addrr(bernreal(i,prec), divrs2_safe(u, i+1)); /* u / (i+1)(i+2) */
947: if (low_stack(avlim,stack_lim(av2,3)))
948: {
949: if(DEBUGMEM>1) err(warnmem,"czeta");
950: tes = gerepileuptoleaf(av2, tes);
951: }
952: }
953: u = gmul(gmul(tes,invn2), gmul2n(mulii(s0, addsi(-1,s0)), -1));
954: tes = gmulsg(nn, gaddsg(1, u));
955: }
956: else /* typ(s0) != t_INT */
957: {
958: s1 = gsub(gmul2n(s,1), unr);
959: s2 = gmul(s, gsub(s,unr));
960: s3 = gmul2n(invn2,3);
961: av2 = avma; avlim = stack_lim(av2,3);
962: s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
963: s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
964: for (i = lim2-2; i>=2; i -= 2)
965: {
966: s5 = gsub(s5, s4);
967: s4 = gsub(s4, s3);
968: tes = gadd(bernreal(i,prec), gdivgs(gmul(s5,tes), (i+1)*(i+2)));
969: if (low_stack(avlim,stack_lim(av2,3)))
970: {
971: GEN *gptr[3]; gptr[0]=&tes; gptr[1]=&s5; gptr[2]=&s4;
972: if(DEBUGMEM>1) err(warnmem,"czeta");
973: gerepilemany(av2,gptr,3);
974: }
975: }
976: u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
977: tes = gmulsg(nn, gaddsg(1, u));
978: }
979: if (DEBUGLEVEL) msgtimer("Bernoulli sum");
980: /* y += tes a / (s-1) */
981: y = gadd(y, gmul(tes, gdiv(a, gsub(s, unr))));
982:
983: END:
984: if (funeq)
985: {
986: GEN pitemp = mppi(prec); setexpo(pitemp,2); /* 2Pi */
987: y = gmul(gmul(y, ggamma(s,prec)), gpui(pitemp,gneg(s),prec));
988: setexpo(pitemp,0); /* Pi/2 */
989: y = gmul2n(gmul(y, gcos(gmul(pitemp,s),prec)), 1);
990: }
991: gaffect(y,res); avma = av; return res;
992: }
993: #endif
994:
995: GEN
996: gzeta(GEN x, long prec)
997: {
998: if (gcmp1(x)) err(talker,"argument equal to one in zeta");
999: switch(typ(x))
1000: {
1001: case t_INT:
1002: if (is_bigint(x))
1003: {
1004: if (signe(x) > 0) return realun(prec);
1005: if (signe(x) < 0 && mod2(x) == 0) return realzero(prec);
1006: }
1007: return izeta(itos(x),prec);
1008:
1009: case t_REAL: case t_COMPLEX:
1010: return czeta(x,prec);
1011:
1012: case t_INTMOD: case t_PADIC:
1013: err(typeer,"gzeta");
1014: case t_SER:
1015: err(impl,"zeta of power series");
1016: }
1017: return transc(gzeta,x,prec);
1018: }
1019:
1020: void
1021: gzetaz(GEN x, GEN y)
1022: {
1023: long av=avma, prec = precision(y);
1024:
1025: if (!prec) err(infprecer,"gzetaz");
1026: gaffect(gzeta(x,prec),y); avma=av;
1027: }
1028:
1029: /***********************************************************************/
1030: /** **/
1031: /** FONCTIONS POLYLOGARITHME **/
1032: /** **/
1033: /***********************************************************************/
1034:
1035: /* validity domain contains .005 < |x| < 230
1036: * Li_m(x = e^z) = sum_n=0 zeta(m-n) z^n / n!
1037: * with zeta(1) := H_m - log(-z) */
1038: static GEN
1039: cxpolylog(long m, GEN x, long prec)
1040: {
1041: long av=avma,li,i,n,bern_upto;
1042: GEN p1,z,h,q,s;
1043:
1044: if (gcmp1(x)) return izeta(m,prec);
1045: z=glog(x,prec); h=gneg_i(glog(gneg_i(z),prec));
1046: for (i=1; i<m; i++) h = gadd(h,ginv(stoi(i)));
1047:
1048: bern_upto=m+50; mpbern(bern_upto,prec);
1049: q=gun; s=izeta(m,prec);
1050: for (n=1; n<=m+1; n++)
1051: {
1052: q=gdivgs(gmul(q,z),n);
1053: s=gadd(s,gmul((n==(m-1))? h: izeta(m-n,prec),q));
1054: }
1055:
1056: n = m+3; z=gsqr(z); li = -(bit_accuracy(prec)+1);
1057:
1058: for(;;)
1059: {
1060: q = gdivgs(gmul(q,z),(n-1)*n);
1061: p1 = gmul(izeta(m-n,prec),q);
1062: s = gadd(s,p1);
1063: if (gexpo(p1) < li) break;
1064: n+=2;
1065: if (n>=bern_upto) { bern_upto += 50; mpbern(bern_upto,prec); }
1066: }
1067: return gerepileupto(av,s);
1068: }
1069:
1070: GEN
1071: polylog(long m, GEN x, long prec)
1072: {
1073: long av,av1,limpile,l,e,i,G,sx;
1074: GEN X,z,p1,p2,n,y,logx;
1075:
1076: if (m<0) err(talker,"negative index in polylog");
1077: if (!m) return gneg(ghalf);
1078: if (gcmp0(x)) return gcopy(x);
1079: av = avma;
1080: if (m==1)
1081: return gerepileupto(av, gneg(glog(gsub(gun,x), prec)));
1082:
1083: l = precision(x);
1084: if (!l) { l=prec; x=gmul(x, realun(l)); }
1085: e = gexpo(gnorm(x)); if (!e || e== -1) return cxpolylog(m,x,prec);
1086: X = (e > 0)? ginv(x): x;
1087: G = -bit_accuracy(l);
1088: n = icopy(gun);
1089: av1=avma; limpile=stack_lim(av1,1);
1090: y = p1 = X;
1091: for (i=2; ; i++)
1092: {
1093: n[2] = i; p1 = gmul(X,p1); p2 = gdiv(p1,gpuigs(n,m));
1094: y = gadd(y,p2);
1095: if (gexpo(p2) <= G) break;
1096:
1097: if (low_stack(limpile, stack_lim(av1,1)))
1098: { GEN *gptr[2]; gptr[0]=&y; gptr[1]=&p1;
1099: if(DEBUGMEM>1) err(warnmem,"polylog");
1100: gerepilemany(av1,gptr,2);
1101: }
1102: }
1103: if (e < 0) return gerepileupto(av, y);
1104:
1105: sx = gsigne(gimag(x));
1106: if (!sx)
1107: {
1108: if (m&1) sx = gsigne(gsub(gun,greal(x)));
1109: else sx = - gsigne(greal(x));
1110: }
1111: z = cgetg(3,t_COMPLEX);
1112: z[1] = zero;
1113: z[2] = ldivri(mppi(l), mpfact(m-1));
1114: if (sx < 0) z[2] = lnegr((GEN)z[2]);
1115: logx = glog(x,l);
1116:
1117: if (m == 2)
1118: { /* same formula as below, written more efficiently */
1119: y = gneg_i(y);
1120: p1 = gmul2n(gsqr(gsub(logx, z)), -1); /* = (log(-x))^2 / 2 */
1121: p1 = gadd(divrs(gsqr(mppi(l)), 6), p1);
1122: p1 = gneg_i(p1);
1123: }
1124: else
1125: {
1126: GEN logx2 = gsqr(logx); p1 = gneg_i(ghalf);
1127: for (i=m-2; i>=0; i-=2)
1128: p1 = gadd(izeta(m-i,l), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
1129: if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
1130: p1 = gadd(gmul2n(p1,1), gmul(z,gpuigs(logx,m-1)));
1131: }
1132:
1133: return gerepileupto(av, gadd(y,p1));
1134: }
1135:
1136: GEN
1137: dilog(GEN x, long prec)
1138: {
1139: return gpolylog(2, x, prec);
1140: }
1141:
1142: GEN
1143: polylogd0(long m, GEN x, long flag, long prec)
1144: {
1145: long k,l,av,fl,m2;
1146: GEN p1,p2,p3,y;
1147:
1148: m2=m&1; av=avma;
1149: if (gcmp0(x)) return gcopy(x);
1150: if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
1151: l=precision(x);
1152: if (!l) { l=prec; x=gmul(x,realun(l)); }
1153: p1=gabs(x,prec); fl=0;
1154: if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
1155:
1156: p1=gneg_i(glog(p1,prec)); p2=gun;
1157: y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
1158: for (k=1; k<m; k++)
1159: {
1160: p2=gdivgs(gmul(p2,p1),k);
1161: p3=m2?greal(polylog(m-k,x,prec)):gimag(polylog(m-k,x,prec));
1162: y=gadd(y,gmul(p2,p3));
1163: }
1164: if (m2)
1165: {
1166: if (flag)
1167: p2 = gdivgs(gmul(p2,p1),-2*m);
1168: else
1169: p2 = gdivgs(gmul(glog(gnorm(gsub(gun,x)),prec),p2),2*m);
1170: y=gadd(y,p2);
1171: }
1172: if (fl) y = gneg(y);
1173: return gerepileupto(av, y);
1174: }
1175:
1176: GEN
1177: polylogd(long m, GEN x, long prec)
1178: {
1179: return polylogd0(m,x,0,prec);
1180: }
1181:
1182: GEN
1183: polylogdold(long m, GEN x, long prec)
1184: {
1185: return polylogd0(m,x,1,prec);
1186: }
1187:
1188: GEN
1189: polylogp(long m, GEN x, long prec)
1190: {
1191: long k,l,av,fl,m2;
1192: GEN p1,y;
1193:
1194: m2=m&1; av=avma;
1195: if (gcmp0(x)) return gcopy(x);
1196: if (gcmp1(x) && m>=2) return m2?izeta(m,prec):gzero;
1197: l=precision(x);
1198: if (!l) { l=prec; x=gmul(x,realun(l)); }
1199: p1=gabs(x,prec); fl=0;
1200: if (gcmpgs(p1,1)>0) { x=ginv(x); p1=gabs(x,prec); fl=!m2; }
1201:
1202: p1=gmul2n(glog(p1,prec),1); mpbern(m>>1,prec);
1203: y=polylog(m,x,prec); y=m2?greal(y):gimag(y);
1204:
1205: if (m==1)
1206: {
1207: p1=gmul2n(p1,-2); y=gadd(y,p1);
1208: }
1209: else
1210: {
1211: GEN p2=gun, p3, p4, p5, p51=cgetr(prec);
1212:
1213: for (k=1; k<m; k++)
1214: {
1215: p2=gdivgs(gmul(p2,p1),k);
1216: if (!(k&1) || k==1)
1217: {
1218: if (k!=1)
1219: {
1220: p5=(GEN)bern(k>>1);
1221: if (bernzone[2]>prec) { affrr(p5,p51); p5=p51; }
1222: p4=gmul(p2,p5);
1223: }
1224: else p4=gneg_i(gmul2n(p2,-1));
1225: p3=polylog(m-k,x,prec); p3=m2?greal(p3):gimag(p3);
1226: y=gadd(y,gmul(p4,p3));
1227: }
1228: }
1229: }
1230: if (fl) y = gneg(y);
1231: return gerepileupto(av, y);
1232: }
1233:
1234: GEN
1235: gpolylog(long m, GEN x, long prec)
1236: {
1237: long i,lx,av=avma,v,n;
1238: GEN y,p1,p2;
1239:
1240: if (m<=0)
1241: {
1242: p1=polx[0]; p2=gsub(gun,p1);
1243: for (i=1; i<=(-m); i++)
1244: p1=gmul(polx[0],gadd(gmul(p2,derivpol(p1)),gmulsg(i,p1)));
1245: p1=gdiv(p1,gpuigs(p2,1-m));
1246: return gerepileupto(av, poleval(p1,x));
1247: }
1248:
1249: switch(typ(x))
1250: {
1251: case t_INT: case t_REAL: case t_FRAC: case t_FRACN:
1252: case t_COMPLEX: case t_QUAD:
1253: return polylog(m,x,prec);
1254:
1255: case t_INTMOD: case t_PADIC:
1256: err(impl, "padic polylogarithm");
1257: case t_POLMOD:
1258: p1=roots((GEN)x[1],prec); lx=lg(p1); p2=cgetg(lx,t_COL);
1259: for (i=1; i<lx; i++) p2[i]=lpoleval((GEN)x[2],(GEN)p1[i]);
1260: y=cgetg(lx,t_COL);
1261: for (i=1; i<lx; i++) y[i]=(long)polylog(m,(GEN)p2[i],prec);
1262: return gerepileupto(av, y);
1263:
1264: case t_POL: case t_RFRAC: case t_RFRACN:
1265: p1=tayl(x,gvar(x),precdl);
1266: return gerepileupto(av, gpolylog(m,p1,prec));
1267:
1268: case t_SER:
1269: if (!m) return gneg(ghalf);
1270: if (m==1)
1271: {
1272: p1=glog(gsub(gun,x),prec);
1273: return gerepileupto(av, gneg(p1));
1274: }
1275: if (valp(x)<=0) err(impl,"polylog around a!=0");
1276: v=varn(x); n=(lg(x)-2)/valp(x); y=ggrando(polx[v],lg(x)-2);
1277: for (i=n; i>=1; i--)
1278: {
1279: p1=gadd(gpuigs(stoi(i),-m),y);
1280: y=gmul(x,p1);
1281: }
1282: return gerepileupto(av, y);
1283:
1284: case t_VEC: case t_COL: case t_MAT:
1285: lx=lg(x); y=cgetg(lx,typ(x));
1286: for (i=1; i<lx; i++)
1287: y[i]=(long)gpolylog(m,(GEN)x[i],prec);
1288: return y;
1289: }
1290: err(typeer,"gpolylog");
1291: return NULL; /* not reached */
1292: }
1293:
1294: void
1295: gpolylogz(long m, GEN x, GEN y)
1296: {
1297: long av=avma, prec = precision(y);
1298:
1299: if (!prec) err(infprecer,"gpolylogz");
1300: gaffect(gpolylog(m,x,prec),y); avma=av;
1301: }
1302:
1303: GEN
1304: polylog0(long m, GEN x, long flag, long prec)
1305: {
1306: switch(flag)
1307: {
1308: case 0: return gpolylog(m,x,prec);
1309: case 1: return polylogd(m,x,prec);
1310: case 2: return polylogdold(m,x,prec);
1311: case 3: return polylogp(m,x,prec);
1312: default: err(flagerr,"polylog");
1313: }
1314: return NULL; /* not reached */
1315: }
1316:
1317: static GEN
1318: qq(GEN x, long prec)
1319: {
1320: long tx=typ(x);
1321:
1322: if (tx==t_PADIC) return x;
1323: if (is_scalar_t(tx))
1324: {
1325: long l=precision(x);
1326: GEN p1,p2,q;
1327:
1328: if (tx != t_COMPLEX || gcmp((GEN)x[2],gzero)<=0)
1329: err(talker,"argument must belong to upper half-plane");
1330:
1331: if (!l) l=prec; p1=mppi(l); setexpo(p1,2);
1332: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */
1333: q=gmul(x,p2); return gexp(q,prec);
1334: }
1335: if (tx != t_POL && tx != t_SER)
1336: err(talker,"bad argument for modular function");
1337:
1338: return tayl(x,gvar(x),precdl);
1339: }
1340:
1341: static GEN
1342: inteta(GEN q)
1343: {
1344: long tx=typ(q);
1345: GEN p1,ps,qn,y;
1346:
1347: y=gun; qn=gun; ps=gun;
1348: if (tx==t_PADIC)
1349: {
1350: if (valp(q) <= 0)
1351: err(talker,"non-positive valuation in inteta");
1352: for(;;)
1353: {
1354: p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
1355: y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
1356: p1=y; y=gadd(y,ps); if (gegal(p1,y)) return y;
1357: }
1358: }
1359: else
1360: {
1361: long l = 0, v = 0; /* gcc -Wall */
1362: long av = avma, lim = stack_lim(av,3);
1363:
1364: if (is_scalar_t(tx)) l = -bit_accuracy(precision(q));
1365: else
1366: {
1367: v=gvar(q); tx = 0;
1368: if (valp(q) <= 0)
1369: err(talker,"non-positive valuation in inteta");
1370: }
1371: for(;;)
1372: {
1373: p1=gneg_i(gmul(ps,gmul(q,gsqr(qn))));
1374: y=gadd(y,p1); qn=gmul(qn,q); ps=gmul(p1,qn);
1375: y=gadd(y,ps);
1376: if (tx)
1377: { if (gexpo(ps)-gexpo(y) < l) return y; }
1378: else
1379: { if (gval(ps,v)>=precdl) return y; }
1380: if (low_stack(lim, stack_lim(av,3)))
1381: {
1382: GEN *gptr[3];
1383: if(DEBUGMEM>1) err(warnmem,"inteta");
1384: gptr[0]=&y; gptr[1]=&qn; gptr[2]=&ps;
1385: gerepilemany(av,gptr,3);
1386: }
1387: }
1388: }
1389: }
1390:
1391: GEN
1392: eta(GEN x, long prec)
1393: {
1394: long av = avma;
1395: GEN q = qq(x,prec);
1396: return gerepileupto(av,inteta(q));
1397: }
1398:
1399: /* returns the true value of eta(x) for Im(x) > 0, using reduction */
1400: GEN
1401: trueeta(GEN x, long prec)
1402: {
1403: long tx=typ(x), av=avma, tetpil,l;
1404: GEN p1,p2,q,q24,n,z,m,unapprox;
1405:
1406: if (!is_scalar_t(tx)) err(typeer,"trueeta");
1407: if (tx != t_COMPLEX || gsigne((GEN)x[2])<=0)
1408: err(talker,"argument must belong to upper half-plane");
1409: l=precision(x); if (l) prec=l;
1410: p1=mppi(prec); setexpo(p1,2);
1411: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=(long)p1; /* 2*I*Pi */
1412: z=gexp(gdivgs(p2,24),prec); /* exp(2*I*Pi/24) */
1413: unapprox=gsub(gun,gpuigs(stoi(10),-8));
1414: m=gun;
1415: for(;;)
1416: {
1417: n=ground(greal(x));
1418: if (signe(n)) {x=gsub(x,n); m=gmul(m,powgi(z,n));}
1419: if (gcmp(gnorm(x), unapprox)>=0) break;
1420: m=gmul(m,gsqrt(gdiv(gi,x),prec)); x=gdivsg(-1,x);
1421: }
1422: q=gmul(p2,x);
1423: q24=gexp(gdivgs(q,24),prec); q=gpuigs(q24,24);
1424: p1=gmul(m,q24); q = inteta(q); tetpil=avma;
1425: return gerepile(av,tetpil,gmul(p1,q));
1426: }
1427:
1428: GEN
1429: eta0(GEN x, long flag,long prec)
1430: {
1431: return flag? trueeta(x,prec): eta(x,prec);
1432: }
1433:
1434: GEN
1435: jell(GEN x, long prec)
1436: {
1437: long av=avma,tetpil,tx = typ(x);
1438: GEN p1;
1439:
1440: if (!is_scalar_t(tx) || tx == t_PADIC)
1441: {
1442: GEN p2,q = qq(x,prec);
1443: p1=gdiv(inteta(gsqr(q)), inteta(q));
1444: p1=gmul2n(gsqr(p1),1);
1445: p1=gmul(q,gpuigs(p1,12));
1446: p2=gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
1447: p1=gmulsg(48,p1); tetpil=avma;
1448: return gerepile(av,tetpil,gadd(p2,p1));
1449: }
1450: p1=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
1451: p1=gsqr(gsqr(gsqr(p1)));
1452: p1=gadd(gmul2n(gsqr(p1),8), ginv(p1)); tetpil=avma;
1453: return gerepile(av,tetpil,gpuigs(p1,3));
1454: }
1455:
1456: GEN
1457: wf2(GEN x, long prec)
1458: {
1459: long av=avma,tetpil;
1460: GEN p1,p2;
1461:
1462: p1=gsqrt(gdeux,prec);
1463: p2=gdiv(trueeta(gmul2n(x,1),prec),trueeta(x,prec));
1464: tetpil=avma;
1465: return gerepile(av,tetpil,gmul(p1,p2));
1466: }
1467:
1468: GEN
1469: wf1(GEN x, long prec)
1470: {
1471: long av=avma,tetpil;
1472: GEN p1,p2;
1473:
1474: p1=trueeta(gmul2n(x,-1),prec); p2=trueeta(x,prec);
1475: tetpil=avma;
1476: return gerepile(av,tetpil,gdiv(p1,p2));
1477: }
1478:
1479: GEN
1480: wf(GEN x, long prec)
1481: {
1482: long av=avma,tetpil;
1483: GEN p1,p2;
1484:
1485: p1=gdiv(trueeta(gmul2n(gaddgs(x,1),-1),prec),trueeta(x,prec));
1486: p2=cgetg(3,t_COMPLEX); p2[1]=zero; p2[2]=ldivrs(mppi(prec),-24);
1487: p2=gexp(p2,prec); tetpil=avma;
1488: return gerepile(av,tetpil,gmul(p1,p2));
1489: }
1490:
1491: GEN
1492: weber0(GEN x, long flag,long prec)
1493: {
1494: switch(flag)
1495: {
1496: case 0: return wf(x,prec);
1497: case 1: return wf1(x,prec);
1498: case 2: return wf2(x,prec);
1499: default: err(flagerr,"weber");
1500: }
1501: return NULL; /* not reached */
1502: }
1503:
1504: static GEN
1505: sagm(GEN x, long prec)
1506: {
1507: GEN p1,a,b,a1,b1,y;
1508: long l,ep;
1509: ulong av;
1510:
1511: if (gcmp0(x)) return gcopy(x);
1512: switch(typ(x))
1513: {
1514: case t_REAL:
1515: l = precision(x); y = cgetr(l); av=avma;
1516: a1 = x; b1 = realun(l);
1517: l = 5-bit_accuracy(prec);
1518: do
1519: {
1520: a=a1; b=b1; a1 = addrr(a,b);
1521: setexpo(a1,expo(a1)-1);
1522: b1=mpsqrt(mulrr(a,b));
1523: }
1524: while (expo(subrr(b1,a1))-expo(b1) >= l);
1525: affrr(a1,y); avma=av; return y;
1526:
1527: case t_COMPLEX:
1528: if (gcmp0((GEN)x[2]))
1529: return transc(sagm,(GEN)x[1],prec);
1530: av=avma; l=precision(x); if (l) prec=l;
1531: a1=x; b1=gun; l = 5-bit_accuracy(prec);
1532: do
1533: {
1534: a=a1; b=b1;
1535: a1=gmul2n(gadd(a,b),-1);
1536: b1=gsqrt(gmul(a,b),prec);
1537: }
1538: while (gexpo(gsub(b1,a1))-gexpo(b1) >= l);
1539: return gerepilecopy(av,a1);
1540:
1541: case t_PADIC:
1542: av=avma; a1=x; b1=gun; l=precp(x);
1543: do
1544: {
1545: a=a1; b=b1;
1546: a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),0);
1547: p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
1548: if (ep<=0) { b1=gneg_i(b1); p1=gsub(b1,a1); ep=valp(p1)-valp(b1); }
1549: }
1550: while (ep<l && !gcmp0(p1));
1551: return gerepilecopy(av,a1);
1552:
1553: case t_SER:
1554: av=avma; a1=x; b1=gun; l=lg(x)-2;
1555: do
1556: {
1557: a=a1; b=b1;
1558: a1=gmul2n(gadd(a,b),-1); b1=gsqrt(gmul(a,b),prec);
1559: p1=gsub(b1,a1); ep=valp(p1)-valp(b1);
1560: }
1561: while (ep<l && !gcmp0(p1));
1562: return gerepilecopy(av,a1);
1563:
1564: case t_INTMOD:
1565: err(impl,"agm of mod");
1566: }
1567: return transc(sagm,x,prec);
1568: }
1569:
1570: GEN
1571: agm(GEN x, GEN y, long prec)
1572: {
1573: long av,tetpil, ty=typ(y);
1574: GEN z;
1575:
1576: if (is_matvec_t(ty))
1577: {
1578: ty=typ(x);
1579: if (is_matvec_t(ty)) err(talker,"agm of two vector/matrices");
1580: z=x; x=y; y=z;
1581: }
1582: if (gcmp0(y)) return gcopy(y);
1583: av=avma; z=sagm(gdiv(x,y),prec); tetpil=avma;
1584: return gerepile(av,tetpil,gmul(y,z));
1585: }
1586:
1587: GEN
1588: logagm(GEN q)
1589: {
1590: long av=avma,prec=lg(q),tetpil,s,n,lim;
1591: GEN y,q4,q1;
1592:
1593: if (typ(q)!=t_REAL) err(typeer,"logagm");
1594: if (signe(q)<=0) err(talker,"non positive argument in logagm");
1595: if (expo(q)<0) s=1; else { q=ginv(q); s=0; }
1596: lim = - (bit_accuracy(prec) >> 1);
1597: q1 = NULL; /* gcc -Wall */
1598: for (n=0; expo(q)>=lim; n++) { q1=q; q=gsqr(q); }
1599: q4=gmul2n(q,2);
1600: if (!n) q1=gsqrt(q,prec);
1601: y=divrr(mppi(prec), agm(addsr(1,q4),gmul2n(q1,2),prec));
1602: tetpil=avma; y=gmul2n(y,-n); if (s) setsigne(y,-1);
1603: return gerepile(av,tetpil,y);
1604: }
1605:
1606: GEN
1607: glogagm(GEN x, long prec)
1608: {
1609: long av,tetpil;
1610: GEN y,p1,p2;
1611:
1612: switch(typ(x))
1613: {
1614: case t_REAL:
1615: if (signe(x)>=0) return logagm(x);
1616:
1617: y=cgetg(3,t_COMPLEX); y[2]=lmppi(lg(x));
1618: setsigne(x,1); y[1]=(long)logagm(x);
1619: setsigne(x,-1); return y;
1620:
1621: case t_COMPLEX:
1622: y=cgetg(3,t_COMPLEX); y[2]=larg(x,prec);
1623: av=avma; p1=glogagm(gnorm(x),prec); tetpil=avma;
1624: y[1]=lpile(av,tetpil,gmul2n(p1,-1));
1625: return y;
1626:
1627: case t_PADIC:
1628: return palog(x);
1629:
1630: case t_SER:
1631: av=avma; if (valp(x)) err(negexper,"glogagm");
1632: p1=gdiv(derivser(x),x);
1633: p1=integ(p1,varn(x));
1634: if (gcmp1((GEN)x[2])) return gerepileupto(av,p1);
1635: p2=glogagm((GEN)x[2],prec); tetpil=avma;
1636: return gerepile(av,tetpil,gadd(p1,p2));
1637:
1638: case t_INTMOD:
1639: err(typeer,"glogagm");
1640: }
1641: return transc(glogagm,x,prec);
1642: }
1643:
1644: GEN
1645: theta(GEN q, GEN z, long prec)
1646: {
1647: long av=avma,tetpil,l,n;
1648: GEN ps,qn,qnold,y,zy,lq,ps2,p1,k,zold;
1649:
1650: l=precision(q); if (l) prec=l;
1651: p1=realun(prec); z=gmul(p1,z);
1652: if (!l) q=gmul(p1,q);
1653: if (gexpo(q)>=0) err(thetaer1);
1654: zy = gimag(z);
1655: zold = NULL; /* gcc -Wall */
1656: if (gcmp0(zy)) k=gzero;
1657: else
1658: {
1659: lq=glog(q,prec); k=ground(gdiv(zy,greal(lq)));
1660: if (!gcmp0(k)) { zold=z; z=gadd(z,gdiv(gmul(lq,k),gi)); }
1661: }
1662: y=gsin(z,prec); n=0; qn=gun;
1663: ps2=gsqr(q); ps=gneg_i(ps2);
1664: do
1665: {
1666: n++; p1=gsin(gmulsg(2*n+1,z),prec);
1667: qnold=qn; qn=gmul(qn,ps);
1668: ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
1669: }
1670: while (gexpo(qnold) >= -bit_accuracy(prec));
1671: if (signe(k))
1672: {
1673: y=gmul(y,gmul(gpui(q,gsqr(k),prec),
1674: gexp(gmul2n(gmul(gmul(gi,zold),k),1),prec)));
1675: if (mod2(k)) y=gneg_i(y);
1676: }
1677: p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
1678: return gerepile(av,tetpil,gmul(p1,y));
1679: }
1680:
1681: GEN
1682: thetanullk(GEN q, long k, long prec)
1683: {
1684: long av=avma,tetpil,l,n;
1685: GEN p1,ps,qn,y,ps2;
1686:
1687: l=precision(q);
1688: if (!l) { l=prec; q=gmul(q,realun(l)); }
1689: if (gexpo(q)>=0) err(thetaer1);
1690:
1691: if (!(k&1)) { avma = av; return gzero; }
1692: ps2=gsqr(q); ps=gneg_i(ps2);
1693: qn = gun; y = gun; n = 0;
1694:
1695: do
1696: {
1697: n++; p1=gpuigs(stoi(n+n+1),k); qn=gmul(qn,ps);
1698: ps=gmul(ps,ps2); p1=gmul(p1,qn); y=gadd(y,p1);
1699: }
1700: while (gexpo(p1) >= -bit_accuracy(l));
1701: p1=gmul2n(gsqrt(gsqrt(q,prec),prec),1); tetpil=avma;
1702: if (k&2) { p1=gneg_i(p1); tetpil=avma; }
1703: return gerepile(av,tetpil,gmul(p1,y));
1704: }
1705:
1706: GEN
1707: jbesselh(GEN n, GEN z, long prec)
1708: {
1709: long av,tetpil,k,l,i,lz;
1710: GEN y,p1,p2,zinv,p0,s,c;
1711:
1712: if (typ(n)!=t_INT) err(talker,"not an integer index in jbesselh");
1713: k=itos(n); if (k<0) err(impl,"ybesselh");
1714:
1715: switch(typ(z))
1716: {
1717: case t_REAL: case t_COMPLEX:
1718: if (gcmp0(z)) return gzero;
1719: av=avma; zinv=ginv(z);
1720: l=precision(z); if (l>prec) prec=l;
1721: gsincos(z,&s,&c,prec);
1722: p1=gmul(zinv,s);
1723: if (k)
1724: {
1725: p0=p1; p1=gmul(zinv,gsub(p0,c));
1726: for (i=2; i<=k; i++)
1727: {
1728: p2=gsub(gmul(gmulsg(2*i-1,zinv),p1),p0);
1729: p0=p1; p1=p2;
1730: }
1731: }
1732: p2=gsqrt(gdiv(gmul2n(z,1),mppi(prec)),prec);
1733: tetpil=avma; return gerepile(av,tetpil,gmul(p2,p1));
1734:
1735: case t_VEC: case t_COL: case t_MAT:
1736: lz=lg(z); y=cgetg(lz,typ(z));
1737: for (i=1; i<lz; i++)
1738: y[i]=(long)jbesselh(n,(GEN)z[i],prec);
1739: return y;
1740:
1741: case t_INT: case t_FRAC: case t_FRACN:
1742: av=avma; p1=cgetr(prec); gaffect(z,p1); tetpil=avma;
1743: return gerepile(av,tetpil,jbesselh(n,p1,prec));
1744:
1745: case t_QUAD:
1746: av=avma; p1=gmul(z,realun(prec)); tetpil=avma;
1747: return gerepile(av,tetpil,jbesselh(n,p1,prec));
1748:
1749: case t_POL: case t_RFRAC: case t_RFRACN:
1750: av=avma; p1=tayl(z,gvar(z),precdl); tetpil=avma;
1751: return gerepile(av,tetpil,jbesselh(n,p1,prec));
1752:
1753: case t_POLMOD:
1754: av=avma; p1=roots((GEN)z[1],prec); lz=lg(p1); p2=cgetg(lz,t_COL);
1755: for (i=1; i<lz; i++) p2[i]=lpoleval((GEN)z[2],(GEN)p1[i]);
1756: tetpil=avma; y=cgetg(lz,t_COL);
1757: for (i=1; i<lz; i++) y[i]=(long)jbesselh(n,(GEN)p2[i],prec);
1758: return gerepile(av,tetpil,y);
1759:
1760: case t_PADIC:
1761: err(impl,"p-adic jbessel function");
1762: case t_SER:
1763: err(impl,"jbessel of power series");
1764: }
1765: err(typeer,"jbesselh");
1766: return NULL; /* not reached */
1767: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to trans3.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