Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 64

PARI/GP code for Elliptic functions

{
\\ Header
my(thetaall_,jacobi_,ellpointtoz_,ellR_,elllemn_,ellsigma_);
\\ Jacobi theta functions
theta1=((q,z)->theta(q,z));
theta2=((q,z)->my(s,t);s=real(z);if(s>0,z=-z);t=theta(q,Pi/2+z);if(!imag(q)&&real(q)>0&&(!s||!real(z)),real(t),t));
theta3=((q,z)->theta4(-q,z));
theta4=((q,z)->my(s,t);s=imag(z);if(s>0,z=-z);t=-I*exp(I*z)*q^(1/4)*theta(q,z-I*log(q)/2);if(!imag(q)&&(!s||!real(z)),real(t),t));
thetaall_=if(arity(theta)>2,
	((z,tau)->call(theta,[z,tau,0])),
	((z,tau)->my(q=exp(I*Pi*tau));z*=Pi;[theta3(q,z),theta4(q,z),theta2(q,z),-theta1(q,z)])
);
\\ Elliptic nome
ellnome=((k)->my(m=k^2);if(k,exp(-Pi*agm(1,sqrt(1-m))/agm(1,k)),m/16));
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)->my(r,e,p,t);if(r=!imag(q),q=real(q));e=exponent(1.*q);p=bitprecision(q);if(p==oo,p=getlocalbitprec());if(!q||-e>p,return(4*sqrt(q)));t=ceil(p-e/2+10);if(r&&q<=1&&2*t*log(2)*(1-q)<Pi^2,localbitprec(p);return(Pi/2));localbitprec(t);t=-I*log(bitprecision(q,t))/(4*Pi);t=2*I*(logeta(t+1/4)-logeta(t-1/4))+Pi/12;bitprecision(if(r,abs(t)*sqrt(sign(q)),t),p));
\\ j-invariant
invellj=((j)->if(!j,return((-1)^(2/3)));my(a=tan(asin(12^(3/2)/sqrt(j))/3));a*=2/(a+sqrt(3));I*agm(1,sqrt(1-a))/agm(1,sqrt(a)));
\\ 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=1.*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+if(9*imag(t)<bitprecision(t),log(eta(t)),0)
);
\\ Neville theta functions
thetaS=((k,z)->my(m=k^2,a,q);if(m==0,sin(z),m==1,sinh(z),a=agm(1,sqrt(1-m));q=exp(-Pi*a/agm(1,k));theta1(q,z*a)*sqrt(a)*q^(-1/4)*(q/m)^(1/4)*(1-m)^(-1/4)));
thetaC=((k,z)->my(m=k^2,a,q);if(m==0,cos(z),m==1,1.,a=agm(1,sqrt(1-m));q=exp(-Pi*a/agm(1,k));theta2(q,z*a)*sqrt(a)*q^(-1/4)*(q/m)^(1/4)));
thetaD=((k,z)->my(m=k^2,a,q);if(m==0,1.,m==1,1.,a=agm(1,sqrt(1-m));q=exp(-Pi*a/agm(1,k));theta3(q,z*a)*sqrt(a)));
thetaN=((k,z)->my(m=k^2,a,q);if(m==0,1.,m==1,cosh(z),a=agm(1,sqrt(1-m));q=exp(-Pi*a/agm(1,k));theta4(q,z*a)*sqrt(a)*(1-m)^(-1/4)));
\\ Jacobi elliptic functions
jacobi_=((k,z)->if(real(k)<0,k=-k);my(m=k^2,m1=1-m,w1=2*ellK(k),w2=I*Pi/agm(1,k),w,v1,v2,v,a,b);if(real(m)>real(m1),if(norm(m)>1,b=sign(imag(w1));w1-=b*w2),norm(m1)>1,a=sign(real(w2));w2-=a*w1);w=[w1,w2];[v1,v2]=v=round(Mat([real(w),imag(w)]~)^(-1)*[real(z),imag(z)]~);[w,z-w*v,v1-a*v2,v2-b*v1,a,b]);
jacobiSN=((k,z)->my(w,p,t,m=k^2,v1,v2);if(m==0,return(sin(z)),m==1,return(tanh(z)));[w,z,v1,v2]=jacobi_(k,z);t=(-1)^v1;z*=1.;my(lp=getlocalbitprec(),m1=1+m,e=exponent(z));if(!z||max(0,exponent(m1))+2*e<-lp,return(t*z));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,v1,v2);if(m==0,return(cos(z)),m==1,return(1/cosh(z)));[w,z,v1,v2]=jacobi_(k,z);t=(-1)^(v1+v2);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,v1,v2);if(m==0,return(tan(z)),m==1,return(sinh(z)));[w,z,v1,v2]=jacobi_(k,z);t=(-1)^v2;z*=1.;my(lp=getlocalbitprec(),m1=2-m,e=exponent(z));if(!z||max(0,exponent(m1))+2*e<-lp,return(t*z));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,v1,v2);if(m==0,return(1.),m==1,return(1/cosh(z)));[w,z,v1,v2]=jacobi_(k,z);t=(-1)^v2;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)->my(m=k^2);if(m==0,return(0.),m==1,return(tanh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-m)),w=2*[a,-b],t=ellzeta(w,z+b),e=elleta(w)/2);if(imag(k)||imag(z),t+e[2],real(t))-e[1]/a*z);
jacobiEpsilon=((k,z)->my(m=k^2);if(m==0,return(z),m==1,return(tanh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-m)),w=2*[a,-b],t=ellzeta(w,z+b));if(imag(k)||imag(z),t+elleta(w)[2]/2,real(t))+z*(2-m)/3);
jacobiAM=((k,z)->my(w,p,m=k^2,v1,v2,a,b,n,q);if(m==0,return(z),m==1,return(2*atan(tanh(z/2))));[w,z,v1,v2,a,b]=jacobi_(k,z);z*=1.;my(lp=getlocalbitprec(),e=exponent(z));q=z;until(1,if(!z||exponent(m)+2*e<-lp,break);e-=exponent(w[2]);if(e<0,[w,z]=bitprecision([w,z],lp-e));p=iferr(ellwp(w,z/2,1),E,break,errname(E)=="e_DOMAIN");q=2*atan(2*((2*m-1)/3-p[1])/p[2]);n=2*(real(z/w[1]-q/Pi)\/2));if(v2%2,my(w1=w[1],w2=w[2]+a*w1);n=2*(imag(z/w2)\imag(w1/w2))-1-n;q=-q);(v1+n)*Pi+q);
\\ Elliptic exponential and logarithm
ellpointtoz_=((E,P,w=ellperiods(E))->my(z=ellpointtoz(E,P));z-w*round(Mat([real(w),imag(w)]~)^(-1)*[real(z),imag(z)]~));
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(Pol([z^2,-1+E.a1*z+E.a2*z^2,E.a3*z+E.a4*z^2,E.a6*z^2])),r,x);foreach(R,t,my(a=abs(t));if(a>r,r=a;x=t));ellpointtoz_(E,[x,-x/z]));
\\ Incomplete elliptic integrals
incellF=((k,phi)->my(m=k^2,a=2*m-1,b=m*(m-1),E=ellinit([0,a,0,b,0]));phi*=1.;if(E===[],if(real(a)<0,phi,2*atanh(tan(phi/2))),my(n,x,c,d,w1,w2);n=real(phi)\/Pi;w1=2*ellK(k);phi-=n*Pi;x=k*phi;if(!x||bitprecision(x)<-2*exponent(x),return(phi+n*w1));x=cotan(phi/2);c=x^2-a;d=sqrt(c^2-4*b);if(real(c*conj(d))<0,d=-d);c=(c+d)/2;a=2*ellpointtoz(E,[c,-c*x]);if(real(k)<0,k=-k);w2=I*Pi/agm(1,k);b=imag(a/w1)\/imag(w2/w1);a-=b*w2;if(b%2,a=w1-a);n+=2*(real(phi/Pi-a/w1)\/2);a+n*w1));
incellE=((k,phi)->jacobiEpsilon(k,incellF(k,phi)));
incellD=((k,phi)->my(m=k^2,u,p,e);if(!m,return(phi/2-sin(2*phi)/4));u=incellF(k,phi);if(!u,return(u^3/3));p=bitprecision(u);e=exponent(m);if(e<0,p-=e);localbitprec(p);[k,u]=bitprecision([k,u],p);(u-jacobiEpsilon(k,u))/m);
jacobiZeta=((k,phi)->jacobiZN(k,incellF(k,phi)));
invellwp=((w,z)->w=ellperiods(w);my(g2=elleisnum(w,4)/12,g3=-elleisnum(w,6)/216,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=ellpointtoz_(E,[x,-y/2],w);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)->my(m=k^2,m1=1-m,kp=sqrt(m1),w=ellK(kp),s=norm(v)<abs(kp),u=I*(incellF(kp,asin(if(s,v/kp,1/v)))-w));if(real(u/w)<0,u=-u);if(s,u=ellK(k)-u);if(!imag(k)&&!imag(v)&&(v=real(v))<=1&&v>=real(kp),real(u),u));
\\ Carlson elliptic integrals
ellR_=((x,y,z)->my(t=x+y+z,E=ellinit([0,t,0,x*y+x*z+y*z,x*y*z]),c);c=ellpointtoz_(E,[0,-sqrt(x)*sqrt(y)*sqrt(z)]);if(real(c)<0,c=-c);[c,t/3,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)->my(m=k^2,lp=getlocalbitprec(),e=exponent(m));if(e<-lp,return(Pi/4));if(e<0,lp-=e);localbitprec(lp);k=bitprecision(k,lp);(ellK(k)-ellE(k))/m);
\\ Lemniscatic elliptic functions
elllemn_=((z,n)->my(lp=getlocalbitprec(),b,L);localbitprec(b=lp+64);z=bitprecision(1.*z,b);bitprecision(if(z&&-4*exponent(z)<=lp,z/=ellK(I);z-=I*(b=round(imag(z)));b%=4;L=thetaall_(z/2,I);if(b%2,n=n+2);if(n==1,b=2-b);I^b*L[5-n]/L[n],if(n==1,z,z*=z;1.-z/(1.+z/2))),lp));
sinlemn=((z)->elllemn_(z,1));
coslemn=((z)->elllemn_(z,2));
asinlemn=((z)->z*hypergeom([1/4,1/2],5/4,z^4));
acoslemn=((z)->ellK(I)-asinlemn(z));
\\ Associated Weierstrass sigma functions
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='x)->ellsigma_(w,z,[1,0]));
ellsigma2=((w,z='x)->ellsigma_(w,z,[0,1]));
ellsigma3=((w,z='x)->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,g,m,n);if(!imag(a)&&abs(a)>1,g=2*Pi/a^2;s=g*hypergeom([2/3,4/3],1,1+1/a^3)/sqrt(27);if(a>0,m=real(s);n=I*abs(imag(s));[a-2*m,a+m-n,a+m+n],n=g*I*hypergeom([2/3,4/3],2,-1/a^3)/9;[a+s,a+s+n,a-2*s-n]),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));
invDixonSM=((a,z)->my(b=(1+a^3)/3,E=ellinit([-a,-a^2,b,a*b,-b^2/3]));elllognum(E,z));
\\ q-Pochhammer
qpoch=((a,q)->my(p=getlocalbitprec(),r=1,t=1);[a,q]=bitprecision([-a,q],p+10);localbitprec(p+10);bitprecision((1+a)*(1+suminf(i=1,r*=q;t*=a*r/(1-r);t)),p));
qfactorial=((n,q)->my(r=1,t=1);while(n-->0,t=1+q*t;r*=t);r);
qbinomial=((n,m,q)->m=min(m,n-m);if(m<0,return(0));prod(i=0,m-1,(1-q^(n-i))/(1-q^(i+1))));
qgamma=((z,q)->my(t);if(norm(q)>1,t=1/q;return(self()(z,t)*t^(-(z-1)*(z-2)/2)));t=1-q;if(!t,return(gamma(z)));t^(1-z)*qpoch(q,q)/qpoch(q^z,q));
qpsi=((z,q)->my(t);if(norm(q)>1,t=1/q;return((3/2-z)*log(t)+self()(z,t)));t=1-q;if(!t,return(psi(z)));z=q^(z-1);-log(t)+log(q)*suminf(n=1,z^n*q^(n^2)*(1-z*q^(2*n))/((1-q^n)*(1-z*q^n))));
\\ 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));
\\ Inverse of incomplete elliptic integrals
invincellE=((k,y)->my(t=y*Pi/(2*ellE(k)));solve(x=t-1/3,t+1/3,incellE(k,x)-y));
}

Last edited by lanxiyu (2025-06-17 12:42:43)