You are not logged in.
Pages: 1
{
\\ Header
my(jacobi_,ellpointtoz_,ellR_,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));
\\ 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)->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)->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],t=ellzeta(w,z+b));if(imag(k)||imag(z),t-ellzeta(w,b),real(t))-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,b],t=ellzeta(w,z+b));if(imag(k)||imag(z),t-ellzeta(w,b),real(t))+z*(2-k^2)/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,n,x,x2,a,b,w1,w2,E);if(m==1,2*atanh(tan(phi/2)),m,n=round(real(phi)/Pi);w1=2*ellK(k);phi-=n*Pi;if(!phi,return(1.*phi+n*w1));x=tan(phi/2);x2=x^2;a=1-2*m;b=a*x2;b=(1+b+sqrt(1+2*b+x2^2))/(2*x2);E=ellinit(-[0,a,0,m*(1-m),0]);a=2*ellpointtoz(E,[b,-b/x]);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,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)/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)->(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],z],r-q);iferr(p=sqrt(ellwp(w,z));z=bitprecision(sign(real(z*p))/p,r),E,,errname(E)=="e_DOMAIN"));s*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");if(real(z)*imag(z)*imag(p)>=0,p=real(p));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
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));
\\ 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));
\\ 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 (2025-03-25 18:49:35)
Offline
Pages: 1