Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 48

PARI/GP code for Elliptic functions

{
\\ 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(w,p,t,m=k^2);if(m==0,return(sin(z)),m==1,return(tanh(z)));w=2*[ellK(k),I*ellK(sqrt(1-m))];[z,t]=ellreducez(w,z);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/2,1),E,return(t*z),errname(E)=="e_DOMAIN");t*p[2]/(m-(p[1]+m1/3)^2));
jacobiCN=((k,z)->my(w,p,t,m=k^2);if(m==0,return(cos(z)),m==1,return(1/cosh(z)));w=2*[ellK(k),I*ellK(sqrt(1-m))];[z,t]=ellreducez(w,z);t=(-1)^(t[1]+t[2]);p=iferr(ellwp(w,z/2),E,return(t*1.),errname(E)=="e_DOMAIN");t*(1+2*(p+(1-2*m)/3)/(m-(p+(1+m)/3)^2)));
jacobiTN=((k,z)->my(w,p,t,m=k^2);if(m==0,return(tan(z)),m==1,return(sinh(z)));w=2*[ellK(k),I*ellK(sqrt(1-m))];[z,t]=ellreducez(w,z);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/2,1),E,return(t*z),errname(E)=="e_DOMAIN");t*p[2]/(1-m-(p[1]-m1/3)^2));
jacobiDN=((k,z)->my(w,p,t,m=k^2);if(m==0,return(1.),m==1,return(1/cosh(z)));w=2*[ellK(k),I*ellK(sqrt(1-m))];[z,t]=ellreducez(w,z);t=(-1)^t[2];p=iferr(ellwp(w,z/2),E,return(t*1.),errname(E)=="e_DOMAIN");t*(1+2*m*((m-2)/3+p)/(m-(p+(1+m)/3)^2)));
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,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(m=k^2,a=agm(1,if(real(k)<0,-k,k)),b=agm(1,sqrt(1-m)),c,d,t,rb);if(real(m)<0.5,z*=2*I*b;b*=Pi/a;rb=real(b);t=real(z)\/(2*rb);z-=2*t*b;d=bitprecision(z);c=sum(n=0,ceil(log(2)*(if(d==oo,getlocalbitprec(),d)/2+abs(real(z)))/rb)+1,my(s=(2*n+1)*b);(-1)^n*(log(-expm1(z-s))-log(-expm1(-z-s))))-z/2;c=if(t%2,c=((imag(z)-real(z)*imag(b)/rb)\(2*Pi)*2+1)*Pi-I*c,I*c);if(imag(m)||real(z),c,real(c)),c=a*z/2;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
DixonLambda=((a)->if(a==-1,return([2*Pi/sqrt(27)]));my(s=sqrt(3),g=gamma(1/3)^3,m=g*hypergeom([1/3,1/3],2/3,-a^3)/Pi,n=Pi^2*a*hypergeom([2/3,2/3],4/3,-a^3)/(3*g));[s*m/6+4*n,-(s-3*I)*m/12-2*(1+s*I)*n,-(s+3*I)*m/12-2*(1-s*I)*n]);
DixonKappa=((a)->if(a==-1,return([2*Pi/sqrt(27)-1]));my(s=sqrt(3),g=gamma(1/3)^3,m=Pi^2*hypergeom([-1/3,2/3],1/3,-a^3)/(3*g),n=a^2*g*hypergeom([1/3,4/3],5/3,-a^3)/Pi);[a+4*m+s*n/12,a-2*(1+s*I)*m-(s-3*I)*n/24,a-2*(1-s*I)*m-(s+3*I)*n/24]);
DixonSM=((a,z)->if(a==-1,z*=sqrt(3)/2;return(sin(z)/sin(Pi/3+z)));my(x,y,b=a^3);iferr([x,y]=ellwp(ellinit([a*(8-b)/48,(b*(20+b)-8)/864]),z,1),E,return(.),errname(E)=="e_DOMAIN");-(a^2/2+2*x)/(y+a*x-(a^3+4)/12));
DixonCM=((a,z)->if(a==-1,z*=sqrt(3)/2;return(sin(Pi/3-z)/sin(Pi/3+z)));my(x,y,b=a^3);iferr([x,y]=ellwp(ellinit([a*(8-b)/48,(b*(20+b)-8)/864]),z,1),E,return(1.),errname(E)=="e_DOMAIN");b=(a^3+4)/12;(y-a*x+b)/(y+a*x-b));
DixonF=((a,z)->if(a==-1,my(t=z*sqrt(3)/2);return(z-sin(t)/sin(Pi/3+t)));my(E=ellinit([-a*(24+a^3*27)/16,(8+a^3*(36+a^3*27))/32]),l=DixonLambda(a)[1],k=DixonKappa(a)[1]);3/4*a*(2+a*(z+l))+ellzeta(E,z+l)-k);
DixonThetaS=((a,z)->if(a==-1,my(t=sqrt(3)/2);return(exp(z*(z-1)/2)/t*sin(z*t)));my(E=ellinit([-a*(24+a^3*27)/16,(8+a^3*(36+a^3*27))/32]));exp((4*a*z+3*a^2*z^2)/8)*ellsigma(E,z));
DixonThetaC=((a,z)->if(a==-1,my(t=sqrt(3)/2);return(exp(z*(z-1)/2)/t*sin(Pi/3-z*t)));my(E=ellinit([-a*(24+a^3*27)/16,(8+a^3*(36+a^3*27))/32]),l=DixonLambda(a)[1],k=DixonKappa(a)[1]);exp((4*a*z+3*a^2*(l-z)^2-4*(k-a)*(l-2*z))/8)*ellsigma(E,l-z));
DixonThetaM=((a,z)->if(a==-1,my(t=sqrt(3)/2);return(exp(z*(z-1)/2)/t*sin(Pi/3+z*t)));my(E=ellinit([-a*(24+a^3*27)/16,(8+a^3*(36+a^3*27))/32]),l=DixonLambda(a)[1],k=DixonKappa(a)[1]);exp((4*a*z+3*a^2*(l+z)^2-4*(k-a)*(l+2*z))/8)*ellsigma(E,l+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 (2024-04-14 00:34:08)