Weierstrass p-function and its first three derivatives with JavaScript.

weierstrass.js

/*
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
  };
  
})();