Last active
December 20, 2023 01:56
-
-
Save stla/bc49500c72a6f9bc38e166ed03f86b98 to your computer and use it in GitHub Desktop.
Weierstrass p-function and its first three derivatives with JavaScript.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
Weierstrass p-function and its first three derivatives. | |
Author: Stéphane Laurent. | |
Depends on math.js <https://mathjs.org/> | |
and jacobi.js <https://gist.github.com/stla/955af9c1c2713a47883bc927a759bdc7>. | |
*/ | |
var weierstrass = (function () { | |
const jtheta1 = jacobi.jtheta1; | |
const jtheta2 = jacobi.jtheta2; | |
const jtheta3 = jacobi.jtheta3; | |
const jtheta4 = jacobi.jtheta4; | |
const g2_from_omega1_and_tau = function(omega1, tau) { | |
const j2square = math.square(jtheta2(0, tau)); | |
const j3square = math.square(jtheta3(0, tau)); | |
const j2fourth = math.square(j2square); | |
const j3fourth = math.square(j3square); | |
const omega1square = math.square(omega1); | |
const pisquare = Math.PI * Math.PI; | |
return math.multiply( | |
math.divide(pisquare*pisquare/12, math.square(omega1square)), | |
math.subtract( | |
math.add(math.square(j2fourth), math.square(j3fourth)), | |
math.multiply(j2fourth, j3fourth) | |
) | |
); | |
}; | |
const wp_from_tau = function(z, tau) { // wp(z, omega1 = 1/2, omega2 = tau/2) | |
const j2square = math.square(jtheta2(0, tau)); | |
const j3square = math.square(jtheta3(0, tau)); | |
const j4overj1 = math.divide(jtheta4(z, tau), jtheta1(z, tau)); | |
const pisquare = Math.PI * Math.PI; | |
return math.subtract( | |
math.multiply( | |
pisquare, | |
math.multiply( | |
j2square, | |
math.multiply( | |
j3square, | |
math.square(j4overj1) | |
) | |
) | |
), | |
math.multiply( | |
pisquare / 3, | |
math.add( | |
math.square(j2square), math.square(j3square) | |
) | |
) | |
); | |
}; | |
const wp_from_omega1_and_tau = function(z, omega1, tau) { | |
const double_omega1 = math.multiply(2, omega1); | |
return math.divide( | |
wp_from_tau( | |
math.divide(z, double_omega1), tau | |
), | |
math.square(double_omega1) | |
); | |
}; | |
const weierDerivative = function(z, omega1, tau) { | |
const double_omega1 = math.multiply(2, omega1); | |
const w1 = math.divide(double_omega1, Math.PI); | |
const z1 = math.unaryMinus(math.divide(z, double_omega1)); | |
const j1 = jtheta1(z1, tau); | |
const j2 = jtheta2(z1, tau); | |
const j3 = jtheta3(z1, tau); | |
const j4 = jtheta4(z1, tau); | |
const j1prime0 = jacobi.jtheta1Prime0(tau); | |
const f = math.divide( | |
math.multiply(j1prime0, math.square(j1prime0)), | |
math.multiply( | |
math.multiply(j1, math.square(j1)), | |
math.multiply( | |
jtheta2(0, tau), | |
math.multiply(jtheta3(0, tau), jtheta4(0, tau)) | |
) | |
) | |
); | |
return math.divide( | |
math.multiply( | |
math.multiply(2, f), math.multiply(j2, math.multiply(j3, j4)) | |
), | |
math.multiply(w1, math.square(w1)) | |
); | |
}; | |
const wp = function(z, omega1, omega2, derivative = 0) { | |
if([0, 1, 2, 3].indexOf(derivative) === -1) { | |
throw "The value of `derivative` must be an integer between 0 and 3."; | |
} | |
const tau = math.divide(omega2, omega1); | |
if(math.im(tau) <= 0) { | |
throw "The ratio `omega2/omega1` must have a nonnegative imaginary part."; | |
} | |
let weier; | |
if(derivative !== 1) { | |
weier = wp_from_omega1_and_tau(z, omega1, tau); | |
if(derivative === 0) { | |
return weier; | |
} | |
if(derivative === 2) { | |
const g2 = g2_from_omega1_and_tau(omega1, tau); | |
return math.subtract( | |
math.multiply(6, math.square(weier)), | |
math.divide(g2, 2) | |
); | |
} | |
} | |
const weierPrime = weierDerivative(z, omega1, tau); | |
if(derivative === 1) { | |
return weierPrime; | |
} // else derivative = 3 | |
return math.multiply(12, math.multiply(weier, weierPrime)); | |
}; | |
const distance = function(a, b) { | |
return math.abs(math.subtract(a, b)); | |
}; | |
const agm = function(a, b) { | |
const eps = Number.EPSILON; | |
while(distance(a, b) >= eps) { | |
const a1 = math.divide(math.add(a, b), 2); | |
const b1 = math.sqrt(math.multiply(a, b)); | |
if(distance(a, a1) < eps && distance(b, b1) < eps){ | |
break; | |
} | |
a = a1; | |
b = b1; | |
} | |
return math.divide(math.add(a, b), 2); | |
}; | |
const eisensteinG4 = function(tau) { | |
const pisquare = Math.PI * Math.PI; | |
const j2 = math.square(math.square(math.square(jtheta2(0, tau)))); | |
const j3 = math.square(math.square(math.square(jtheta3(0, tau)))); | |
const j4 = math.square(math.square(math.square(jtheta4(0, tau)))); | |
return math.multiply( | |
pisquare * pisquare / 90, | |
math.add(j2, math.add(j3, j4)) | |
); | |
}; | |
const eisensteinG6 = function(tau) { | |
const pithird = Math.PI * Math.PI * Math.PI; | |
const j2fourth = math.square(math.square(jtheta2(0, tau))); | |
const j3fourth = math.square(math.square(jtheta3(0, tau))); | |
const j4fourth = math.square(math.square(jtheta4(0, tau))); | |
return math.multiply( | |
pithird * pithird / 945, | |
math.subtract( | |
math.add( | |
math.multiply(j3fourth, math.square(j3fourth)), | |
math.multiply(j4fourth, math.square(j4fourth)) | |
), | |
math.multiply( | |
math.multiply(3, math.square(j2fourth)), | |
math.add(j3fourth, j4fourth) | |
) | |
) | |
); | |
}; | |
const kleinjinv = function(j) { | |
if(math.equal(j, 0)) { | |
return math.complex(0.5, Math.sqrt(3)/2); | |
} | |
const j2 = math.square(j); | |
const j3 = math.multiply(j, j2); | |
const t = math.subtract( | |
math.add( | |
math.multiply(2304, j2), | |
math.multiply( | |
12288, | |
math.sqrt( | |
math.multiply( | |
3, | |
math.subtract( | |
math.multiply(1728, j2), j3 | |
) | |
) | |
) | |
) | |
), | |
math.add(j3, math.multiply(884736, j)) | |
); | |
const u = math.pow(t, 1/3); | |
const x = math.subtract( | |
math.add( | |
math.divide(u, 768), | |
math.subtract(1, math.divide(j, 768)) | |
), | |
math.divide( | |
math.subtract(math.multiply(1536, j), j2), | |
math.multiply(768, u) | |
) | |
); | |
const lbd = math.multiply( | |
-0.5, | |
math.subtract( | |
-1, math.sqrt(math.subtract(1, math.multiply(4, x))) | |
) | |
); | |
return math.multiply( | |
math.complex(0, 1), | |
math.divide( | |
agm(1, math.sqrt(math.subtract(1, lbd))), | |
agm(1, math.sqrt(lbd)) | |
) | |
); | |
}; | |
const omega1_and_tau = function(g2, g3) { | |
if(math.equal(g2, 0)) { | |
const omega1 = math.divide( | |
1.52995403705719335008, | |
math.pow(g3, 1/6) | |
); | |
const tau = math.complex(0.5, Math.sqrt(3)/2); | |
return [omega1, tau]; | |
} | |
const g2cube = math.multiply(g2, math.square(g2)); | |
const j = math.divide( | |
math.multiply(1728, g2cube), | |
math.subtract(g2cube, math.multiply(27, math.square(g3))) | |
); | |
const tau = kleinjinv(j); | |
let omega1; | |
if(math.equal(g3, 0)) { | |
omega1 = math.multiply( | |
math.complex(0, 1), | |
math.sqrt( | |
math.sqrt( | |
math.multiply( | |
math.divide(15/4, g2), | |
eisensteinG4(tau) | |
) | |
) | |
) | |
); | |
} else { | |
omega1 = math.sqrt( | |
math.divide( | |
math.multiply( | |
7, math.multiply(eisensteinG6(tau), g2) | |
), | |
math.multiply( | |
12, math.multiply(eisensteinG4(tau), g3) | |
) | |
) | |
); | |
} | |
return [omega1, tau]; | |
}; | |
const halfPeriods = function(g2, g3) { | |
const omega1_tau = omega1_and_tau(g2, g3); | |
return { | |
omega1: omega1_tau[0], | |
omega2: math.multiply(omega1_tau[0], omega1_tau[1]) | |
} | |
}; | |
return { | |
agm: agm, | |
halfPeriods: halfPeriods, | |
wp: wp | |
}; | |
})(); |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment