You are not logged in.
Pages: 1
{
\\ Jacobi theta functions
theta1=((q,z)->theta(q,z));
theta2=((q,z)->2*q^(1/4)*suminf(n=0,q^(n*n+n)*cos((2*n+1)*z)));
theta3=((q,z)->1+2*suminf(n=1,q^n^2*cos(2*n*z)));
theta4=((q,z)->theta3(-q,z));
\\ Lattice reduction
ellreducez=((w,z)->my(t=ellperiods(w),a,b);b=Mat([real(t),imag(t)]~);a=round(Mat([real(w),imag(w)]~)^(-1)*b);b=round(b^(-1)*[real(z),imag(z)]~);[z-t*b,a*b]);
ellperiodratio=((w)->w=ellperiods(w);w[1]/w[2]);
\\ Elliptic nome
ellnome=((k)->if(k,exp(-Pi*agm(1,sqrt(1-k^2))/agm(1,k)),(k^2)/2));
invellnome=((q)->my(s=4*q^(1/2),t);if(s,t=log(q)/(Pi*I);s*=(eta(t/2)*eta(2*t)^2/eta(t)^3)^4);s);
nometomodularangle=((q)->q=sqrt(q);4*suminf(n=0,my(t=2*n+1);(-1)^n*q^t/(t*(1+q^(2*t)))));
\\ j-invariant
invellj=((j)->my(prec=bitprecision(j));if(exponent(j)>getlocalbitprec()/2+10,return(I*log(1.*j-744)/(2*Pi)));if(prec<oo,j=bitprecision(j,prec+64));ellperiodratio(ellinit(ellfromj(j))));
\\ Logarithm of Dedekind eta function
logeta=((t)->
if(imag(t)<=0,error("domain error in modular function: Im(argument) <= 0"));
my(w=[t,1],a,b,c,d,M);t=ellperiods(w);M=round(Mat([real(t),imag(t)]~)^(-1)*Mat([real(w),imag(w)]~));c=M[1,2];t=t[1]/t[2];if(c,M*=sign(c);a=M[1,1];b=M[2,1];c=M[1,2];d=M[2,2];I*Pi*((a+d)/(12*c)+sumdedekind(-d,c))+log(-I*(c*t+d))/2,I*Pi*M[2,1]/12)+Pi*I*t/12+log(eta(t))
);
\\ Neville theta functions
thetaS=((k,z)->if(k,my(m=k^2,a=agm(1,sqrt(1-m)),q=exp(-Pi*a/agm(1,k)));theta1(q,z*a)*sqrt(a)*m^(-1/4)*(1-m)^(-1/4),sin(z)));
thetaC=((k,z)->if(k,my(m=k^2,a=agm(1,sqrt(1-m)),q=exp(-Pi*a/agm(1,k)));theta2(q,z*a)*sqrt(a)*m^(-1/4),cos(z)));
thetaD=((k,z)->if(k,my(m=k^2,a=agm(1,sqrt(1-m)),q=exp(-Pi*a/agm(1,k)));theta3(q,z*a)*sqrt(a),1.));
thetaN=((k,z)->if(k,my(m=k^2,a=agm(1,sqrt(1-m)),q=exp(-Pi*a/agm(1,k)));theta4(q,z*a)*sqrt(a)*(1-m)^(-1/4),1.));
\\ Jacobi elliptic functions
jacobiSN=((k,z)->thetaS(k,z)/thetaN(k,z));
jacobiCN=((k,z)->thetaC(k,z)/thetaN(k,z));
jacobiTN=((k,z)->thetaS(k,z)/thetaC(k,z));
jacobiDN=((k,z)->thetaD(k,z)/thetaN(k,z));
jacobiZN=((k,z)->((z)->thetaN(bitprecision(k,getlocalbitprec()),z))'(z)/thetaN(k,z));
jacobiEpsilon=((k,z)->jacobiZN(k,z)+z*ellE(k)/ellK(k));
jacobiAM=((k,z)->if(k,my(c=thetaC(k,z),s=thetaS(k,z),a=c+I*s,b=c-I*s,t=norm(a)-norm(b));z=if(t<0,I,-I)*log(if(t<0,b,a)/thetaN(k,z));if(t,z,real(z)),z));
\\ Alternative implementations
jacobiSN=((k,z)->my(a,b,w,p,t,m=k^2);if(m==0,return(sin(z)),m==1,return(tanh(z)));a=ellK(k);b=I*ellK(sqrt(1-m));[z,t]=ellreducez(2*[a,b],z);w=[4*a,2*b];t=(-1)^t[1];z*=1.;my(lp=getlocalbitprec(),m1=1+m,e=exponent(z));if(!z||max(0,exponent(m1))+2*e<-lp,return(t*z));w=ellperiods(w);e-=exponent(w[2]);if(e<0,[w,z]=bitprecision([w,z],lp-e));p=iferr(ellwp(w,z,1),E,return(t*z),errname(E)=="e_DOMAIN");-t*2.*(p[1]-m1/6)/p[2]);
jacobiCN=((k,z)->my(a,b,w,p,t,m=k^2);if(m==0,return(cos(z)),m==1,return(1/cosh(z)));a=ellK(k);b=I*ellK(sqrt(1-m));[z,t]=ellreducez(2*[a,b],z);w=2*[a+b,a-b];t=(-1)^(t[1]+t[2]);p=iferr(ellwp(w,z),E,return(t*1.),errname(E)=="e_DOMAIN");t*(1-6/(1+4*m+12*p)));
jacobiTN=((k,z)->my(a,b,w,p,t,m=k^2);if(m==0,return(tan(z)),m==1,return(sinh(z)));a=ellK(k);b=I*ellK(sqrt(1-m));[z,t]=ellreducez(2*[a,b],z);w=[2*a,4*b];t=(-1)^t[2];z*=1.;my(lp=getlocalbitprec(),m1=2-m,e=exponent(z));if(!z||max(0,exponent(m1))+2*e<-lp,return(t*z));w=ellperiods(w);e-=exponent(w[2]);if(e<0,[w,z]=bitprecision([w,z],lp-e));p=iferr(ellwp(w,z,1),E,return(t*z),errname(E)=="e_DOMAIN");-t*2.*(p[1]+m1/6)/p[2]);
jacobiDN=((k,z)->my(a,b,w,p,t,m=k^2);if(m==0,return(1.),m==1,return(1/cosh(z)));a=ellK(k);b=I*ellK(sqrt(1-m));[z,t]=ellreducez(2*[a,b],z);w=[2*a,4*b];t=(-1)^t[2];p=iferr(ellwp(w,z),E,return(t*1.),errname(E)=="e_DOMAIN");t*(1-6*m/(4+m+12*p)));
jacobiZN=((k,z)->if(k==0,return(0.),k^2==1,return(tanh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-k^2)),w=[2*a,2*b]);ellzeta(w,z+b)-ellzeta(w,b)-ellzeta(w,a)/a*z);
jacobiEpsilon=((k,z)->if(k==0,return(z),k^2==1,return(tanh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-k^2)),w=[2*a,2*b]);ellzeta(w,z+b)-ellzeta(w,b)+z*(2-k^2)/3);
jacobiAM=((k,z)->if(k,my(a=agm(1,k),b=agm(1,sqrt(1-k^2)),c=a*z/2,t);if(b,b=a*Pi/(2*b);t=real(c)\/real(b);if(t,c-=t*b;t*=Pi));t+2*(atan(tanh(c))+if(b,suminf(i=1,atan(tanh(c+i*b))+atan(tanh(c-i*b))),0)),z));
\\ Elliptic exponential and logarithm
ellexpnum=((E,z)->my(P=ellztopoint(E,z));if(#P<2,0.,-P[1]/P[2]));
elllognum=((E,z)->if(!z,return(z));my(R=polroots(E.a6*z^2+(E.a3*z+E.a4*z^2)*x+(-1+E.a1*z+E.a2*z^2)*x^2+z^2*x^3),r,x);foreach(R,t,my(a=abs(t));if(a>r,r=a;x=t));ellreducez(E.omega,ellpointtoz(E,[x,-x/z]))[1]);
\\ Incomplete elliptic integrals
incellF=((k,phi)->my(m=k^2,t,s,c,d,E);if(m==1,2*atanh(tan(phi/2)),m,t=round(real(phi)/Pi);phi-=t*Pi;if(!phi,return(1.*phi+2*t*ellK(k)));s=sin(phi);c=s^2/(1+cos(phi));d=1+sqrt(1-m*s^2);E=ellinit(-[(1-m+m^2)/3,(2-3*m-3*m^2+2*m^3)/27]);2*(ellreducez(E.omega,ellpointtoz(E,[d/c-(1+m)/3,s*(m*c-d)/c^2]))[1]+t*ellK(k)),phi*1.));
incellE=((k,phi)->jacobiEpsilon(k,incellF(k,phi)));
incellD=((k,phi)->my(u=incellF(k,phi));(u-jacobiEpsilon(k,u))/k^2);
jacobiZeta=((k,phi)->jacobiZN(k,incellF(k,phi)));
invellwp=((w,z)->w=ellperiods(w);my(g2=elleisnum(w,4,1),g3=elleisnum(w,6,1),E=ellinit(-[g2,g3]/4),x,y,l);z=Vec(z);x=z[1];l=#z<2;y=if(l,sqrt(4*x^3-x*g2-g3),z[2]);z=ellreducez(w,ellpointtoz(E,[x,-y/2]))[1];if(l&&real(z)>=0,z,-z));
invjacobiSN=((k,v)->incellF(k,asin(v)));
invjacobiCN=((k,v)->incellF(k,acos(v)));
invjacobiTN=((k,v)->incellF(k,atan(v)));
invjacobiDN=((k,v)->incellF(1/k,acos(v))/k);
\\ Carlson elliptic integrals
my(ellR_=((x,y,z)->my(t=(x+y+z)/3,[e1,e2,e3]=[x-t,y-t,z-t],E=ellinit([e1*e2+e1*e3+e2*e3,e1*e2*e3]),c);c=ellreducez(E.omega,ellpointtoz(E,[t,-sqrt(x)*sqrt(y)*sqrt(z)]))[1];[c,t,E]));
ellRC=((x,y)->if(x==y,1/sqrt(x),my(sx=sqrt(x),sy=sqrt(y));acos(sx/sy)/(sqrt(1-x/y)*sy)));
ellRD=((x,y,z)->my([u,P,E]=ellR_(x,y,z));3*(u*(z-P)-ellzeta(E,u)+sqrt(x)*sqrt(y)/sqrt(z))/((x-z)*(y-z)));
ellRE=((x,y)->my(sx=sqrt(x),sy=sqrt(y),t=sx*sy);2*((sx+sy)*ellE((sx-sy)/(sx+sy))/Pi-if(t,t/(2*agm(sx,sy)),0)));
ellRF=((x,y,z)->ellR_(x,y,z)[1]);
ellRG=((x,y,z)->my([u,P,E]=ellR_(x,y,z));(u*P+ellzeta(E,u))/2);
ellRJ=((x,y,z,p)->my(c=(!x+!y+!z)/2+!p);3/2*intnum(t=[0,-c],[oo,-5/2],1/(sqrt(t+x)*sqrt(t+y)*sqrt(t+z)*(t+p))));
ellRK=((x,y)->1/agm(sqrt(x),sqrt(y)));
ellRM=((x,y,p)->4/(3*Pi)*ellRJ(0,x,y,p));
\\ Complete elliptic integrals
ellD=((k)->(ellK(k)-ellE(k))/k^2);
\\ Lemniscatic elliptic functions
sinlemn=((z)->my(w=2*ellK(I),q=round(z/w),s=(-1)^(real(q)+imag(q)),p,r);z=1.*(z-q*w);if(z&&-2*(q=exponent(z))<(r=bitprecision(z)),[w,z]=bitprecision([w*[1+I,1-I],z],r-q);iferr(p=ellwp(w,z,1),E,,errname(E)=="e_DOMAIN"));s*if(p,-bitprecision(2*p[1]/p[2],r),z));
coslemn=((z)->my(p=2*ellK(I),q=round(z/p),s=(-1)^(real(q)+imag(q)));z-=q*p;iferr(p=ellwp(p*[1+I,1-I],z),E,return(s*1.),errname(E)=="e_DOMAIN");s*(1-2/(2*p+1)));
asinlemn=((z)->z*hypergeom([1/4,1/2],5/4,z^4));
acoslemn=((z)->ellK(I)-asinlemn(z));
\\ Associated Weierstrass sigma functions
my(ellsigma_=((w,z,i)->my([t,e]=ellperiods(w,1),a=Mat([real(t),imag(t)]~),b=round(a^(-1)*Mat([real(w),imag(w)]~))*Col(i)%2,c);if(real(t[1]/t[2])>0,b[2]*=-1);[a,e]=[t*b,e*b]/2;c=real(z*conj(a));if(iferr(c>0,E,c=t;c[2-b[1]]=a;return(exp(-ellwp(t,a)*z^2/2)*ellsigma(c,z)/ellsigma(t,z)),errname(E)=="e_TYPE2"),z=-z);ellsigma(w,z+a)/ellsigma(w,a)*exp(-z*e)));
ellsigma1=((w,z)->ellsigma_(w,z,[1,0]));
ellsigma2=((w,z)->ellsigma_(w,z,[0,1]));
ellsigma3=((w,z)->ellsigma_(w,z,[1,1]));
\\ Dixon elliptic functions
DixonSM=((z)->my(a,b);iferr([a,b]=ellwp(ellinit([0,-1/108]),z,1),E,return(.),errname(E)=="e_DOMAIN");-6*a/(3*b-1));
DixonCM=((z)->my(a,b);iferr([a,b]=ellwp(ellinit([0,-1/108]),z,1),E,return(1.),errname(E)=="e_DOMAIN");(3*b+1)/(3*b-1));
DixonThetaS=((z)->ellsigma(ellinit([0,1/4]),z));
DixonThetaC=((z)->my(t=gamma(1/3)^3/(sqrt(3)*(2*Pi)),E=ellinit([0,1/4]));exp(Pi*(2*z/t-1)/sqrt(27))*ellsigma(E,t-z));
DixonThetaM=((z)->DixonThetaC(-z));
ArcDixonSM=((z)->z*hypergeom([1/3,2/3],4/3,z^3));
ArcDixonCM=((z)->gamma(1/3)^3/(sqrt(3)*(2*Pi))-ArcDixonSM(z));
\\ q-Pochhammer
qpoch=((a,q)->my(p=getlocalbitprec(),r=1,t=1);[a,q]=precision([-a,q],p+10);localbitprec(p+10);bitprecision((1+a)*(1+suminf(i=1,r*=q;t*=a*r/(1-r);t)),p));
\\ Genus 2 arithmetic geometric mean
genus2agm=((args[..])->my(p=getlocalbitprec());localbitprec(p+64);bitprecision(if(
#args==4,my([a,b,c,d]=bitprecision(args,p+64));while(exponent(abs(b-a)+abs(c-a)+abs(d-a))-exponent(a)>=-p,my(sa=sqrt(a),sb=sqrt(b),sc=sqrt(c),sd=sqrt(d));[a,b,c,d]=[(a+b+c+d)/4,(sa*sb+sc*sd)/2,(sa*sc+sb*sd)/2,(sa*sd+sb*sc)/2]);a,
#args==6,my([a,b,c,d,e,f]=bitprecision(args,p+64));while(exponent(abs(a-b)+abs(c-d)+abs(e-f))-exponent([a,b,c,d,e,f])>=-p,my(ab=a*b,cd=c*d,ef=e*f,ac=a-c,ad=a-d,ae=a-e,af=a-f,bc=b-c,bd=b-d,be=b-e,bf=b-f,ce=c-e,cf=c-f,de=d-e,df=d-f,n1=ab-cd,n2=ab-ef,n3=cd-ef,d1=ac+bd,d2=ae+bf,d3=ce+df,s1=sqrt(ac)*sqrt(ad)*sqrt(bc)*sqrt(bd),s2=sqrt(ae)*sqrt(af)*sqrt(be)*sqrt(bf),s3=sqrt(ce)*sqrt(cf)*sqrt(de)*sqrt(df),r1=[n1+s1,n1-s1]/d1,r2=[n2+s2,n2-s2]/d2,r3=[n3+s3,n3-s3]/d3,E=oo);for(i=1,2,for(j=1,2,for(k=1,2,my(t=abs(r2[i]-r1[j])+abs(r1[3-j]-r3[k])+abs(r3[3-j]-r2[3-k]));if(E>t,E=t;a=r2[i];b=r1[j];c=r1[3-j];d=r3[i];e=r3[3-i];f=r2[3-i];break)))));[a,c,e],
error("genus2agm: must be 4 or 6 arguments")
),p));
}Last edited by lanxiyu (2023-12-04 20:07:21)
Offline
Pages: 1