global(eps,tabx0,tabw0,tabxp,tabxm,tabwp,tabwm,nt,NB,AM); global(MYDEBUG); MYDEBUG=0; /* The init functions have the following purposes: 1) They fill the value tabx0 = phi(0) and arrays of abcissas tabxp[] = phi(k/2^m) (positive x) and also of tabxm[] = phi(-k/2^m) (negative x) unless the phi function is odd, in which case this is useless. 2) They fill the corresponding arrays of weights tabw0 = phi'(0) and tabwp[] = phi'(k/2^m) (and possibly also of tabwm[] = phi'(-k/2^m)). 3) They set the global variable eps to the desired accuracy (depending on the GP default). 4) They compute the global variable nt which says that the weights tabwp[k] and tabwm[k] are negligible with respect to eps if k > nt. In particular the tabxx[] arrays are indexed from 1 to nt+1. */ /* The variable AM is not very important, it should be in [1,2]. According to the authors, the optimal value is Pi/2. */ AM=1; \\ Evaluate the function f at u, and increment global evaluation counter. \\ Heavy in GP, would be better in C. E(f_,u) = NB++; eval(Str(f_,"(",u,")")) \\ phi(t)=tanh(sinh(t)) : from -1 to 1, hence from a to b compact interval. inittanhsinh(m)= { local(h,lim,et,ct,st,ext,ex); eps=10.^(-precision(1.)); h=2.^(-m); lim=ceil(20*2^m); tabxp=vector(lim); tabwp=vector(lim); tabw0=AM; nt=-1; ex=exp(h); et=1.; for (k=1,lim, et*=ex; ct=(et+1/et)/2; st=et-ct; ext=2/(exp(2*AM*st)+1); tabxp[k]=1-ext; tabwp[k]=AM*ct*ext*(1+tabxp[k]); if (exteps,nt=max(k-1,0);break) ); } \\ User-defined change of variable phi(t)=ph(t), where t always goes from \\ -\infty to +\infty and phi(t) goes either from 0 to \infty if fl=0, from \\ -\infty to \infty if fl=1, from -1 to 1 if fl=-1. Default 0. initphigen(ph,m,fl=0)= { local(h,lim,NBI,php,te); php=concat(ph,"'"); eps=10.^(-precision(1.)); h=2.^(-m); lim=20*2^m; NBI=NB; nt=-1; tabxp=vector(lim); tabwp=vector(lim); if (fl>=0, tabxm=vector(lim); tabwm=vector(lim)); tabx0=E(ph,0); tabw0=E(php,0); for (k=1,lim, tabxp[k]=E(ph,k*h); tabwp[k]=E(php,k*h); if (fl>=0, tabxm[k]=E(ph,-k*h); tabwm[k]=E(php,-k*h); te=abs(tabxm[k]); if (fl==1, te=exp(-te)), te=1-tabxp[k] ); if (te0, print("warning: integral from infty to infty or from -infty to -infty"); return (0), return (b[1]*intninfinf(f_,m,fl)) ) ) ) } \\ (2i\pi)^{-1}\int_{|z-a|=R}f(z)dz = \sum_{|z_0-a|