Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 2

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)
\\ 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
\\ 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)
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,-I*log((thetaC(k,z)+I*thetaS(k,z))/thetaN(k,z)),z)
\\ Alternative implementations
jacobiSN(k,z)=if(k==0,return(sin(z)),k^2==1,return(tanh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-k^2)),w=[4*a,2*b],p);p=iferr(ellwp(w,z,1),E,return(0),errname(E)=="e_DOMAIN");-2*(p[1]-ellwp(w,2*a))/p[2]
jacobiCN(k,z)=if(k==0,return(cos(z)),k^2==1,return(1/cosh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-k^2)),w=[4*a,2*(a+b)],p);p=iferr(ellwp(w,z),E,return(1),errname(E)=="e_DOMAIN");(p-ellwp(w,a))/(p-ellwp(w,b))
jacobiDN(k,z)=if(k==0,return(1.),k^2==1,return(1/cosh(z)));my(a=ellK(k),b=I*ellK(sqrt(1-k^2)),w=[2*a,4*b],p);p=iferr(ellwp(w,z),E,return(1),errname(E)=="e_DOMAIN");(p-ellwp(w,a+b))/(p-ellwp(w,b))
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
\\ Incomplete elliptic integrals
incellF(k,phi)=my(s=sin(phi),c=1-cos(phi),d=1+sqrt(1-(k*s)^2),E=ellinit(-[(1-k^2+k^4)/3,(2-3*k^2-3*k^4+2*k^6)/27]));2*ellpointtoz(E,[d/c-(1+k^2)/3,s*(k^2*c-d)/c^2])
incellE(k,phi)=jacobiEpsilon(k,incellF(k,phi))
\\ Lemniscatic elliptic functions
sinlemn(z)=my(p);iferr(p=ellwp(ellinit([1/4,0]),z,1),E,return(0.),errname(E)=="e_DOMAIN");-2*p[1]/p[2]
coslemn(z)=my(p);iferr(p=ellwp(ellinit([1/4,0]),z),E,return(1.),errname(E)=="e_DOMAIN");(2*p-1)/(2*p+1)
asinlemn(z)=z*hypergeom([1/4,1/2],5/4,z^4)
acoslemn(z)=ellK(I)-asinlemn(z)
\\ q-Pochhammer
qpoch(a,q)=my(p=getlocalbitprec(),r=1,t=1);a=precision(-a,p+10);localbitprec(p+10);bitprecision((1+a)*(1+suminf(i=1,r*=q;t*=a*r/(1-r);t)),p)

Last edited by lanxiyu (2023-03-13 03:43:55)