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