\p 67; /* default precision 38 decimal digits */ /* for higher precision, use values of \p 133 or \p 105 */ /* L, B, scnt, strmat, rtrmat, strm, rtrm, sfuncprec, Period idelta are global variables set by init(base), rtrmat set by riemup */ centerat = 0; /* optional centerat value for sexp */ renormup = 1; /* used in the recenterup routine, for experiments, try renormup=0, =2 */ curprecision=0; lastprecision=-2; slowmode = 0; /* old slowmode, idelta changes. new "fastmode" uses idelta=0.12*I, sexpup uses itself */ fastdelta = 12*I/100; /* default idelta for use when (slowmode==0); must be a rational number, or else there is precision loss for \p 134 */ debugloop = 0; /* debug loop plot control */ B = exp(1); etaB = exp(1/exp(1)); L = 0.318 + 1.337*I; idelta = 0.1*I; /* I*imag(exp(Pi*I*(1/(strm*2)))), set by init(b) initialization routine */ imult = 1; /* also set by init(b), corresponds to e^2pi*idelta */ recenter = 0; /* updated by recenterup called by init(b) */ renorm = 1; /* set by init(b) */ gennewsexp = 0; /* use the riemB function for bases0), lastabs = curabs; lastl=L; L=log(L)*rlnB; curabs = abs(lastl-L); s--; ); Period = 2*Pi*I/(L*log(B)+log(log(B))); if (B0), lastabs = curabs; lastl=L2; L2=exp(L2*lnB); curabs = abs(lastl-L2); s--; ); Period2 = -2*Pi*I/(L2*log(B)+log(log(B))); ); /* precision goal bits for loop(t) */ z=1.0; precis=precision(z); y=1.0;while(z<>(z+y),y=y/2); goalbits = abs(log(sqrt(y)*50)/log(2))-10; /* subtract 10 for the last iteration */ print (" base " B); print (" fixed point " L); s=10^-((precis/2)-1.5); setsfuncprec(s); /* set sfuncprec; calculate value for scnt. This is the number iterations required for the superfunction, isuperf routines */ LxlnB = L*lnB; rlogLxlnB = 1/log(LxlnB); /*idelta=0.031*I; */ /* works better for initialization renormup */ idelta=fastdelta; imult = exp(-idelta*I*2*Pi); xterminc = 0; /* xterminc used to stabilize isuperf for bases<1.7 */ if ((B<1.7), xterminc=1); if ((B<1.55), xterminc=2); if ((B<1.49), xterminc=3); if ((B<1.47), xterminc=0); /* limited convergence range for isuperf(z), superf/riemaprx are not effected */ centerat=0; if ((B>20), centerat=-0.5); if ((B>5000), centerat=-0.65); strmmult=3/2; curprecision=0; sieinit; rtrm=3; for (s=1, rtrm, rie[s]=0); print (" Pseudo Period " Period); default(format, "g0.32"); if (BetaB, print ("superf2 only valid for B1, if (imag(z)>0, print ("!!! WARNING, riemaprx(z) much better than sexp(z) for imag(z)>I !!!"), print ("!!! WARNING, splicesexp(z) much better than sexp(z) for imag(z)<-I !!!") ); ); z=z+recenter-centerat; t=0; while ((real(z)<-0.5), z=z+1; t--; ); /* real z converges over a very large range by mapping back to the unit approximation range */ while ((real(z)> 0.5), z=z-1; t++; ); y = z; s = sie[1]; for (i=2,strmat, s=s+sie[i]*y; y=y*z; ); while ((t<0), s=log(s)*rlnB; t++; ); while ((t>0), s=exp(s*lnB); t--; ); return(s); } shortsexp(z) = { /* only used for recenterup generation, needs sexp(z) without recenter */ local (i,s,w,y); y = z; s = sie[1]; for (i=2,strmat, s=s+sie[i]*y; y=y*z; ); return(s); } riemaprx(w) = { local(s,w1,y,y1,z,z1,tot); if (imag(w)<0, print ("!!! WARNING, splicesexp(z) much better than riemaprx(z) for imag(z)<-I !!!")); if ((B>etaB) || gennewsexp, z = exp(w*2*Pi*I); tot = w+rie[1]; y = z; for (s=2,rtrmat, tot=tot+rie[s]*y; y=y*z; ); z = superf(tot); return(z); ); if ((B=ideltai), rieest[t]=riemaprx(scrc[t]+centerat), rieest[t]=precision(sexp(scrc[t]+centerat),precis); /* make sure sexp has full precision before 100's of isuperf log iterations, riemaprx naturally has precision */ ); ); strm = floor(st*strmmult); sie = vector (strm,i,0); strmat=strm; for (s=1,strm, tot=0; for (t=1,st, tot=tot+real(rieest[t]); rieest[t]=rieest[t]*conj(scrc[t]); ); sie[s]=tot/st; if (noboundchk==0, if ((s%2==0), /* odd Taylor series terms should always be positive */ if ((sie[s]<0), strmat=s-2;s=strm); ); if (((s%2) && (s>1)), /* even Taylor series, sie[s-2]<0 => sie[s]<0 */ if (((sie[s]>0) && (sie[s-2]<0)), strmat=s-2;s=strm); ); ); ); } erroravg(w) = { local(tot,x); x = sexp(-0.5); print ("sexp(-0.5)= " x); tot=0; for (s=1,150, z=-0.5-1/300+(s/150)+idelta; tot=tot+abs(sexp(z)-riemaprx(z)); ); tot=-log(tot/150)/log(2); default(format, "g0.10"); print (tot " Riemann/sexp binary precision bits I=" idelta); print (loopcount "=loopcount " recenter " recenter/renorm " renorm); if (precis>67, default(format, "g0.64"), default(format, "g0.32")); return(tot); } ideltaupd(c) = { local(y); y=fastdelta*1.0; if (slowmode, y=exp(Pi*I*(1/(2*c)))); idelta =I*imag(y); imult = exp(-idelta*I*2*Pi); } setsfuncprec(newprec) = { /* set sfuncprec; calculate value for scnt. */ /* This is the number iterations required for the superfunction, isuperf routines */ local(y,s); sfuncprec=newprec; y=0.5; s=0; while ((abs(y-L)>sfuncprec), y=log(y)*rlnB; s++; ); /* sfuncprec defaults to 1E-8 */ scnt=s; if (Bsfuncprec), y=exp(y*lnB); s++; ); /* sfuncprec defaults to 1E-8 */ scnt2=s; ); } riemcontour(z) = { local(y); y=sexp(z); y = isuperf(y)-z; return(y); } slog(z,est) = { local (x,y,s,slop,lastyz,curyz,lest,ly,jt); /* inverse sexp(z), using "est" as initial estimate */ /* est is an optional parameter, important near the singularity */ /* the 0.01 radius used for the slope means won't converge within 0.01 of singularity */ jt=0; if ((imag(z)==0) && (est==0), while (z>B,z=rlnB*log(z);jt++); if (z<0,z=exp(z*lnB);jt=-1); ); lastyz=100; y=splicesexp(est); curyz=abs(y-z); lest=est+curyz*0.05; ly=splicesexp(lest); s=1; while ((curyz>sfuncprec) && ((curyz0.1, print (curyz " bad result, need better initial est, slog(z,est)")); return(est+jt); } /* split the imag(z) into three cases, for optimal precision in the complex plane */ splicesexp(z) = { if (imag(z)>=0.8, riemaprx(z), if (imag(z)<=-0.8, conj(riemaprx(conj(z))), sexp(z))); } sexptaylor( w,r) = { local(rinv,s,t,x,y,z,tot,t_est,tcrc); /* outputs tseries taylor series, the complex taylor series for sexp, */ /* centered at "w", radius "r", use stitchsexp for the complex plane */ /* if the taylor series for slog is desired, replace stitchsexp with slog */ t_est = vector (240,i,0); tcrc = vector (240,i,0); if (r==0,r=1); rinv = 1/r; for(s=1, 240, x=-1/(120*2)+(s/120); tcrc[s]=exp(Pi*I*x); ); /* uses the splicesexp to allow for centering anywhere in the complex plane */ /* if the taylor series for slog is desired, replace stitchsexp with slog here */ for (t=1,240, t_est[t] = splicesexp(w+r*tcrc[t]); ); for (s=1,200, tot=0; for (t=1,240, tot=tot+t_est[t]; t_est[t]=t_est[t]*conj(tcrc[t]); ); tseries[s]=tot/240; ); for (s=2,200, tseries[s]=tseries[s]*(rinv)^(s-1); ); print ("sexp taylor series; first 60 terms of tseries[1..200] centered at " w); if ((imag(w)==0) && (real(w)>-2), tseries=real(tseries), default(format,"g0.22")); for (s=1,61,z=tseries[s];if (s>10, print1("a" s-1 "= "), print1("a" s-1 "= "));if (real(z)<0, print(z), print(" " z))); } slogtaylor(w,r,est) = { local(rinv,s,t,x,y,z,tot,t_est,tcrc); /* outputs tseries taylor series, the complex taylor series for sexp, */ /* centered at "w", radius "r", use stitchsexp for the complex plane */ /* if the taylor series for slog is desired, replace stitchsexp with slog */ /* est is an optional initial estimate for slog(w+r) */ t_est = vector (240,i,0); tcrc = vector (240,i,0); if (r==0,r=1); rinv = 1/r; for(s=1, 240, x=-1/(120*2)+(s/120); tcrc[s]=exp(Pi*I*x); ); /* uses the splicesexp to allow for centering anywhere in the complex plane */ /* if the taylor series for slog is desired, replace stitchsexp with slog here */ if (est==0, est=slog(w+r)); for (t=1,240, est=slog(w+r*tcrc[t],est);t_est[t]=est); for (s=1,200, tot=0; for (t=1,240, tot=tot+t_est[t]; t_est[t]=t_est[t]*conj(tcrc[t]); ); tseries[s]=tot/240; ); for (s=2,200, tseries[s]=tseries[s]*(rinv)^(s-1); ); print ("slog taylor series; first 60 terms of tseries[1..200] centered at " w); if (imag(w)==0, tseries=real(tseries), default(format,"g0.22")); for (s=1,61,z=tseries[s];if (s>10, print1("a" s-1 "= "), print1("a" s-1 "= "));if (real(z)<0, print(z), print(" " z))); } taprx(z) = { /* sexp approximation from taylor series generated with the sexptaylor routine */ local (i,s,y); y = z; s = tseries[1]; for (i=2,200, s=s+tseries[i]*y; y=y*z; ); return(s); } /* Cheta, sexp_eta stuff added to kneser.gp */ /* inicheta, cheta, sexpeta, genpoly, */ /* sexp_e because its the complex base change function, too much fun! */ initcheta(w) = { /* stuff added for base cheta, sexpeta, invcheta, invsexpeta, sexp_e */ /* could be a little more agressive here, aiming for 75% precision, option for usematrix=0 */ chterms = 2*(floor(precis/2)-8)+1; chdelta = (chterms-1)*2; fcheta = genpoly(chterms,chdelta); invchetr = imag(cheta(-chdelta+(chterms-1)*I/2)); invprecis = 10000.*(10^-chterms); /* invprecis = 100*abs(chetaerr(I+0.5,chdelta)); */ /* invprecis = 1E-47; */ chetadlt=cheta(-chdelta); fsexpeta = genpoly(chterms,chdelta,1); sxpetadlt=sexpeta(chdelta); return(0); } cheta(w) = { local(i,x,y,e1); e1=exp(-1); x=w; i=chdelta; while (real(x)>+1/2,x--;i++); while (real(x)<-1/2,x++;i--); y = eval(fcheta); if (i>0, for (s=1,i, y=exp(y*e1))); if (i<0, for (s=1,-i,y=log(y)*exp(1))); return(y); } sexpeta(w) = { local(i,x,y,e,e1); e=exp(1); e1=exp(-1); x=w; i=-chdelta; while (real(x)>+1/2,x--;i++); while (real(x)<-1/2,x++;i--); y = eval(fsexpeta); if (i>0, for (s=1,i, y=exp(y*e1))); if (i<0, for (s=1,-i,y=log(y)*e)); return(y); } genpoly(mterms,mdelta,forsexp) = { local(i,ht,xr,yr,initc,e,e1,fl,m25,s25,r25); xr = vector (mterms,i,0); yr = vector (mterms,i,0); if (mdelta==0,mdelta=chdelta); ht = (mterms-1)/2; e=exp(1); e1=1/e; if (forsexp==0, initc=2*e; for (s=1,mdelta-ht,initc=log(initc)*e); for (i=1,mterms, yr[mterms+1-i]=initc;initc=log(initc)*e); ); if (forsexp, initc=1; for (s=1,mdelta-ht,initc=exp(initc*e1)); for (i=1,mterms, yr[i]=initc;initc=exp(initc*e1)); ); if (usematrix==1, m25=matrix(ht,ht,i,j,i^(2*j)); s25=matrix(ht,1,i,j,(yr[ht+1+i]+yr[ht+1-i])/2-yr[ht+1]); r25=matsolve(m25,s25); fl=0; for (s=1,ht,fl=fl+r25[s,1]*x^(s*2)); m25=matrix(ht,ht,i,j,i^(2*j-1)); s25=matrix(ht,1,i,j,(yr[ht+1+i]-yr[ht+1-i])/2); r25=matsolve(m25,s25); for (s=1,ht,fl=fl+r25[s,1]*x^(s*2-1)); fl = fl + yr[ht+1]; ); if (usematrix==2, m25=matrix(mterms,mterms,i,j,if ((i==1+ht) && (j==1), 1, (i-1-ht)^(j-1))); s25=matrix(mterms,1,i,j,yr[i]); r25=matsolve(m25,s25); fl=0; for (s=1,mterms,fl=fl+r25[s,1]*x^(s-1)); ); if (usematrix==0, for(i=1,mterms,xr[i]=i-1-ht); fl = polinterpolate(xr,yr); ); return(fl); } invcheta(z,est) = { local(i,t,x,mz,e,lastyz,curyz,lest,ly,y); e=exp(1); mz=z; t=0; if (est==0, x=0; mz=z; while ((abs(mz-chetadlt)>invchetr) && (t<1000), t++; mz=log(mz)*e; ); if ((t==1000), print ("too close to e " z " " mz)); est = 2/(chetadlt/e-1) - chdelta - 2/(mz/e-1); ); /* inverse sexp(z), using "est" as initial estimate */ /* est is an optional parameter, important near the singularity */ /* the 0.05 radius used for the slope means won't converge within 0.05 of singularity */ lastyz=100; y=cheta(est); curyz=abs(y-mz); lest=est+curyz*0.05; ly=cheta(lest); i=40; while ((curyz>invprecis) && ((curyzinvchetr) && (t<1000), t++; mz=exp(mz*e1); ); if ((t==1000), print ("too close to e " z " " mz)); est = chdelta - 2/(1-sxpetadlt*e1) + 2/(1-mz*e1); ); /* inverse sexp(z), using "est" as initial estimate */ /* est is an optional parameter, important near the singularity */ /* the 0.05 radius used for the slope means won't converge within 0.05 of singularity */ lastyz=100; y=sexpeta(est); curyz=abs(y-mz); lest=est+curyz*0.05; ly=sexpeta(lest); i=40; while ((curyz>invprecis) && ((curyz4.584,w--); */ /* cheta(z) overflow protection */ x=cheta(z+bias_e+w)/exp(1)-1; if ((w>0) && (imag(z)==0), for (s=1,w,x=log(x))); if ((w>0) && (imag(z)<>0), y=cheta(z+bias_e+w-1)/exp(1)-1; x=log(x); /* while (imag(y)>(imag(x)+Pi), x=x+2*Pi*I); */ /* while (imag(y)<(imag(x)-Pi), x=x-2*Pi*I); */ if (imag(y)>(imag(x)+Pi), x=x+2*Pi*I*floor((imag(y)-imag(x)+Pi)/(2*Pi))); if (imag(y)<(imag(x)-Pi), x=x-2*Pi*I*floor((imag(x)-imag(y)+Pi)/(2*Pi))); for (s=2,w,x=log(x)); ); return(x); } /* gfunc/gie/gtaylor/gaprx is a free extra taylor series routine for use by anything you want */ /* gfunc(z) = { cheta(invcheta(z)+0.5); } */ /* gfunc(z) = { sexp_e(z,3); } */ /* gfunc(z) = { invcheta(z) } */ /* gfunc(z) = { sexpeta(z) } */ /* gfunc(z) = { cheta(z) } */ gfunc(z)=cheta(z); gtaylor( w,r,samples) = { local(rinv,s,t,x,y,z,tot,t_est,tcrc,halfsamples); /* outputs gie taylor series, the complex taylor series for sexp, */ /* centered at "w", radius "r", use stitchsexp for the complex plane */ /* if the taylor series for slog is desired, replace stitchsexp with slog */ /* default halfsamples=240 */ if (samples==0, samples=240); /* no matter how many sample points, the default gie series size is 200 halfsamples */ halfsamples=samples/2; t_est = vector (samples,i,0); tcrc = vector (samples,i,0); gcent = w; if (r==0,r=1); rinv = 1/r; for(s=1, samples, x=-1/(samples)+(s/halfsamples); tcrc[s]=exp(Pi*I*x); ); /* uses the gfunc to allow for centering anywhere in the complex plane */ /* if the taylor series for slog is desired, replace stitchsexp with slog here */ for (t=1,samples, t_est[t] = gfunc(w+r*tcrc[t]); ); for (s=1,200, tot=0; for (t=1,samples, tot=tot+t_est[t]; t_est[t]=t_est[t]*conj(tcrc[t]); ); gie[s]=tot/(samples); ); for (s=2,200, gie[s]=gie[s]*(rinv)^(s-1); ); print ("sexp taylor series; first 60 halfsamples of gie[1..200] centered at " w); default(format,"g0.22"); for (s=1,81,z=gie[s];if (s>10, print1("a" s-1 "= "), print1("a" s-1 "= "));if (real(z)<0, print(z), print(" " z))); } gaprx(z) = { /* inverse ffunc approximation from unit circle Taylor series. Converges nicely for -1sfuncprec) && ((curyz0.1, print (curyz " bad result, need better initial est, slog(z,est)")); return(est); } loop(t) = { local (rt,st,s,xmlt,y,precbits,lastloop,slowprecision); if (loopcount>0, t=t+loopcount); lastloop=0; if ((loopcount==0), if (t==1, lastloop=1); /* t=1, only one loop */ loopcount++; /* this is the initialization loop, use a reasonable set of initial values */ ideltaupd(12); setsfuncprec(1E-8); riemup(18); /* calculate a 18 term rie array from the initial sexp approximation */ sexpup(12);/* calculate a 9 term sie array from the riemaprx */ recenterup; /* recenter and renormzlize the sie array */ lastprecision=0; curprecision = erroravg; /* calculate the current error term, to be used in the next iteration */ ); precbits=0; while ((lastloop==0), if ((curprecisiongoalbits, who cares? */ loopcount=-1; print ("UNEXPECTED LOSS: curprecisiongoalbits, /* if-then-else */ lastloop=1; /* IF: last iteration */ xmlt=lastradius, /* IF: use 0.65 for last iteration, and don't update precision */ setsfuncprec(2^(-precbits)); /* ELSE: this is not the last iteration, update the precision, and the scnt count used for superf/isuperf */ ); /* st will become the new number of terms and sample points in the sexp taylor series array */ y = precbits*xmlt + 1.5 + log(abs(sie[strmat]))/log(2); st = strmat + floor(y); if (st<12, st=12); ideltaupd(st); /* update idelta to match the value in st */ /* determine the value to update for the number of points in the riemaprx taylor series */ y = floor((log(2)/log(imult))*(precbits - 2.5 + log(abs(rie[rtrmat]))/log(2))); riemup(y); sexpup(st); recenterup; /* recenter and renormzlize the sie array */ lastprecision=curprecision; curprecision = erroravg; /* calculate the current error term, to be used in the next iteration */ if (debugloop, if (loopcount>=debugloop, /*ploth(t=-0.6,0.6,x=t+idelta*1.1;y=sexp(x)-riemaprx(x);[real(y),imag(y)]);*/ ploth(t=-0.6,0.6,x=t+idelta*1.1;y=sexp(x)-riemaprx(x);z=isuperf(sexp(x))-isuperf(riemaprx(x));[real(y),imag(y),real(z),imag(z)]); /*ploth(t=-0.6,0.6,x=t+idelta*1.1;z=isuperf(sexp(x))-isuperf(riemaprx(x));[real(z),imag(z)]);*/ ); ); ); if (precis>67, default(format, "g0.64")); return(1); } base_e_2_10(w) = { init(exp(1)); loop; dumparray; init(2); loop(0); dumparray; init(10); loop(0); dumparray; } dumparray(w) = { local (r,s,t); default(format, "g0.10"); print ("sie sexp Taylor series approximation terms"); if (centerat<>0, print ("sexp Tayler series centerat " centerat)); default(format, "g0.64"); write ("kneser.txt", " "); write ("kneser.txt", "base " B); write ("kneser.txt", "fixed point " L); write ("kneser.txt", "Pseudo Period " Period); if (centerat<>0, write ("kneser.txt", "sexp Tayler series centerat " centerat)); write ("kneser.txt", " "); write ("kneser.txt", "iterations " scnt " used for superf/isuperf"); write ("kneser.txt", strmat " strm(s) out of " strm " sexp(z) generates " 2*rtrm " Riemann samples"); write ("kneser.txt", rtrmat " rtrm(s) out of " rtrm " riemaprx(z) generates " strm " sexp samples"); write ("kneser.txt", curprecision " Riemann/sexp binary precision bits"); write ("kneser.txt", loopcount "=loopcount " recenter " recenter/renorm " renorm); write ("kneser.txt", " "); write ("kneser.txt", "sie sexp Taylor series approximation terms"); if (loopcount==-1, write ("kneser.txt", "UNEXPECTED PRECISION LOSS. EXITED EARLY, curprecision=idelta, 0.012*I"); print ("splicesexp(z) sexp/riemaprx(z) valid everywhere in the complex plane"); print ("superf(z) superfunction of z, base B"); print ("isuperf(z) inverse superfunction of z, base B"); print ("dumparray dumps sexp/Riemann Taylor series arrays to kneser.txt file"); print ("cheta(z) base eta functions: cheta(z) sexpeta(z) invcheta(z) invsexpeta(z)"); print ("help print this menu, 'morestuff' for more special functions"); print ("loop or loop(0), loop until optimal precision, ~104 bits, takes ~30 seconds"); print (""); } morestuff(w) = { print (""); print ("morestuff morestuff menu for Kneser.gp program"); print ("functions"); print ("shortsexp(z) used for recenterup generation, sexp(z) without recenter, centerat"); print ("sexpup; updates the sexp(z) function sie array from riemaprx(z) function"); print ("riemup; updates the riemaprx(z) function rie array from sexp(z) function"); print ("recenterup; updates the recenter/renorm values, based on renormup"); print ("new functions"); print ("slog(z,est) inverse sexp(z), using optional "est" as initial estimate"); print ("sexptaylor(w,r) generates tseries sexp taylor series, centered at w, radius r"); print ("slogtaylor(w,r) generates tseries slog taylor series, centered at w, radius r"); print ("taprx(z) evaluates tseries sexptaylor/slogtaylor series at z"); print ("gfunc(z) set to whatever you want, gtaylor(z,r) genseries gaprx(z)"); print ("superf2(z) only for B