Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 54

PARI/GP test cases for Weierstrass elliptic function

{
my(
f=default(format),
p=default(realprecision),
bp=default(realbitprecision),
w,w1,w2,w3,g2,g3,z1,z2,z3,z4,p1,p2,p3,p4,P
);
default(realbitprecision,bp+64);
default(format,f);
setrand(getwalltime());
f=strprintf("%%.%dg\n",p);
w1=tan(random(1.)*Pi/2);
w2=I*tan(random(1.)*Pi/2);
w3=-w1-w2;
w=2*[w1,w2];
g2=elleisnum(w,4,1);
g3=elleisnum(w,6,1);
P=((z,n)->my(a,b,x,y,t);[x,y]=ellwp(w,z,1);t=vector(n);a='x;t[1]=x;for(i=2,n,t[i]=if(i%2,a=(4*'x^3-g2*'x-g3)*deriv(b,'x)+(6*'x^2-g2/2)*b;subst(a,'x,x),b=deriv(a,'x);y*subst(b,'x,x)));return(t));
z1=random(1.)*w1;p1=P(z1,6);
z2=random(1.)*w1;p2=P(z2,4);
z3=random(1.)*w1;p3=P(z3,4);
z4=random(1.)*w1;p4=P(z4,4);

printf(f, matdet([
	1,p1[1];
	1,p2[1]
]));
printf(f, ellsigma(w,z1-z2)*ellsigma(w,z1+z2)/(ellsigma(w,z1)*ellsigma(w,z2))^2);

printf(f, matdet([
	1,p1[1],p1[2];
	1,p2[1],p2[2];
	1,p3[1],p3[2]
]));
printf(f, -2!*ellsigma(w,z1-z2)*ellsigma(w,z1-z3)*ellsigma(w,z2-z3)*ellsigma(w,z1+z2+z3)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z3))^3);

printf(f, matdet([
	1,p1[1],p1[2],p1[3];
	1,p2[1],p2[2],p2[3];
	1,p3[1],p3[2],p3[3];
	1,p4[1],p4[2],p4[3]
]));
printf(f, -2!*3!*ellsigma(w,z1-z2)*ellsigma(w,z1-z3)*ellsigma(w,z1-z4)*ellsigma(w,z2-z3)*ellsigma(w,z2-z4)*ellsigma(w,z3-z4)*ellsigma(w,z1+z2+z3+z4)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z3)*ellsigma(w,z4))^4);

printf(f, matdet([
	1,p1[1];
	0,p1[2]
]));
printf(f, -ellsigma(w,2*z1)/ellsigma(w,z1)^4);

printf(f, matdet([
	1,p1[1],p1[2];
	0,p1[2],p1[3];
	0,p1[3],p1[4]
]));
printf(f, 2!^2*ellsigma(w,3*z1)/ellsigma(w,z1)^9);

printf(f, matdet([
	1,p1[1],p1[2],p1[3];
	0,p1[2],p1[3],p1[4];
	0,p1[3],p1[4],p1[5];
	0,p1[4],p1[5],p1[6]
]));
printf(f, -(2!*3!)^2*ellsigma(w,4*z1)/ellsigma(w,z1)^16);

printf(f, matdet([
	1,p1[1],p1[2];
	1,p2[1],p2[2];
	1,p3[1],p3[2]
])/matdet([
	1,p1[1];
	1,p2[1]
]));
printf(f, -2!*ellsigma(w,z1-z3)*ellsigma(w,z2-z3)*ellsigma(w,z1+z2+z3)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z1+z2)*ellsigma(w,z3)^3));

printf(f, matdet([
	1,p1[1],p1[2],p1[3];
	1,p2[1],p2[2],p2[3];
	1,p3[1],p3[2],p3[3];
	1,p4[1],p4[2],p4[3]
])/matdet([
	1,p1[1],p1[2];
	1,p2[1],p2[2];
	1,p3[1],p3[2]
]));
printf(f, 3!*ellsigma(w,z1-z4)*ellsigma(w,z2-z4)*ellsigma(w,z3-z4)*ellsigma(w,z1+z2+z3+z4)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z3)*ellsigma(w,z1+z2+z3)*ellsigma(w,z4)^4));

printf(f, matdet([
	1,p1[1],p1[2];
	0,p1[2],p1[3];
	1,p2[1],p2[2]
])/matdet([
	1,p1[1];
	0,p1[2]
]));
printf(f, -2!*ellsigma(w,z1-z2)^2*ellsigma(w,2*z1+z2)/(ellsigma(w,z1)^2*ellsigma(w,2*z1)*ellsigma(w,z2)^3));

printf(f, matdet([
	1,p1[1],p1[2],p1[3];
	0,p1[2],p1[3],p1[4];
	0,p1[3],p1[4],p1[5];
	1,p2[1],p2[2],p2[3]
])/matdet([
	1,p1[1],p1[2];
	0,p1[2],p1[3];
	0,p1[3],p1[4]
]));
printf(f, 3!*ellsigma(w,z1-z2)^3*ellsigma(w,3*z1+z2)/(ellsigma(w,z1)^3*ellsigma(w,3*z1)*ellsigma(w,z2)^4));

printf(f, (1/2)*matdet([
	1,p1[2];
	1,p2[2]
])/matdet([
	1,p1[1];
	1,p2[1]
]));
printf(f, ellzeta(w,z1+z2)-ellzeta(w,z1)-ellzeta(w,z2));

printf(f, (1/3)*matdet([
	1,p1[1],p1[3];
	1,p2[1],p2[3];
	1,p3[1],p3[3]
])/matdet([
	1,p1[1],p1[2];
	1,p2[1],p2[2];
	1,p3[1],p3[2]
]));
printf(f, ellzeta(w,z1+z2+z3)-ellzeta(w,z1)-ellzeta(w,z2)-ellzeta(w,z3));

printf(f, (1/4)*matdet([
	1,p1[1],p1[2],p1[4];
	1,p2[1],p2[2],p2[4];
	1,p3[1],p3[2],p3[4];
	1,p4[1],p4[2],p4[4]
])/matdet([
	1,p1[1],p1[2],p1[3];
	1,p2[1],p2[2],p2[3];
	1,p3[1],p3[2],p3[3];
	1,p4[1],p4[2],p4[3]
]));
printf(f, ellzeta(w,z1+z2+z3+z4)-ellzeta(w,z1)-ellzeta(w,z2)-ellzeta(w,z3)-ellzeta(w,z4));

default(realbitprecision,bp);
}

Last edited by lanxiyu (2024-11-28 12:24:59)