Invariant $g_2$ in Weierstrass's elliptic function I programmed the invariant
$$
g_2 = 60\sum_{(m,n)\neq{0,0}} \left(m+ni  \right)^{-4}
$$
(periods $\omega_1=1,\omega_1=i$, https://en.wikipedia.org/wiki/Weierstrass%27s_elliptic_functions)
by (java codes)
public static double[] g2() {
    double[] sum = { 0, 0 };
    int iteration =100; //the result in not sensitive to the number of iterations
    for (int i = 0; i < iteration; i++) {
        for (int j = 0; j < iteration; j++) {
            if (0 == i && 0 == j)
                continue;
            double[] grid = {i, j }; // grid m,n
            double[] s4 = M.square(M.square(grid));
            sum = M.add(sum, M.inverse(s4));
        }
    }
    return M.scale(sum, 60);
}

The program reports $g2=(112.207,0)$ while the correct value is 189.072 (Ferguson p4). Wiki also tells $g_2=2(e_1^2+e_2^2)$, so I am convinced that $g_2=4e_1^2=189.072$ as Hoffman(p113) and other resources say $e_1=6.875$. 
I cannot understand this difference since the formula is quite simple and I have no bugs in this very short program.
Reference:
H.Ferguson，A.Gray，S. Markvorsen, Costa's minimal surface via Mathematica, 1995.
D.A. Hoffman，W. Meeks III, A complete embedded minimal surface in $ \bf {R}^3$ with genus one and three ends, 1985
The complete program:
public class TestConstants3 { // D.A. Hoffman, 1985
    static final int iteration = 100;
    public static void main(String[] args) {
        double[] g2 = g2();
        System.out.println("g2: "+g2[0]);
    }
public static double[] g2() {
    double[] sum = { 0, 0 };
    for (int i = 0; i < iteration; i++) {
        for (int j = 0; j < iteration; j++) {
            if (0 == i && 0 == j)
                continue;
            double[] grid = { i, j }; // grid m,n
            double[] s4 = square(square(grid));
            sum = add(sum, inverse(s4));
        }
    }
    return scale(sum, 60);
}

public static double[] add(double[] a, double[] b) {
    double[] re = new double[a.length];
    for (int i = 0; i < a.length; i++)
        re[i] = a[i] + b[i];
    return re;
}

public static double[] square(double[] c) {
    double x = c[0];
    double y = c[1];
    return new double[] { x * x - y * y, 2 * x * y };
}

public static double[] inverse(double[] c) {
    double x = c[0];
    double y = c[1];
    double d = x * x + y * y;
    return new double[] { x / d, -y / d };
}

public static double[] scale(double[] a, double s) {
    double[] re = new double[a.length];
    for (int i = 0; i < a.length; i++)
        re[i] = a[i] * s;
    return re;
}

}
 A: I will outline an explicit evaluation of such constant through the following preliminary lemmas:
$$ \sum_{n\geq 1}\frac{\coth(\pi n)}{n^3}=\frac{7\pi^3}{180},\qquad \sum_{n\geq 1}\frac{1}{\sinh^2(\pi n)}=\frac{1}{6}-\frac{1}{2\pi} \tag{1}$$
$$ \sum_{n\geq 1}\frac{1}{\sinh^4(\pi n)}=-\frac{11}{90}+\frac{1}{3 \pi }+\frac{\Gamma\left(\frac{1}{4}\right)^8}{1920 \pi ^6}\tag{2}$$
As a reference, we may consider Zucker, The Summation of Series of Hyperbolic Functions.
We have:
$$ \frac{g_2}{60}= 4\zeta(4)+4\sum_{m,n\geq 1}\frac{(m^2+n^2)^2-8m^2 n^2}{\left(m^2+n^2\right)^4}=4K\zeta(2)-32\sum_{m,n\geq 1}\frac{m^2 n^2}{(m^2+n^2)^4}\tag{3}$$
where $K$ is Catalan's constant. On the other hand:
$$\small \sum_{m,n\geq 1}\frac{m^2 n^2}{(m^2+n^2)^4}=\sum_{n\geq 1}\left(\frac{\pi  \coth(\pi n)}{32 n^3}+\frac{\pi^2}{32 n^2 \sinh(\pi n)^2}-\frac{\pi^4}{24\sinh^2(\pi n)}-\frac{\pi^4}{16\sinh^4(\pi n)}\right)\tag{4} $$
hence:
$$\small 32\sum_{m,n\geq 1}\frac{m^2 n^2}{(m^2+n^2)^4}=\frac{176 \pi ^6-3 \Gamma\left(\frac{1}{4}\right)^8}{2880 \pi ^2}+\pi^2\sum_{n\geq 1}\frac{1}{n^2\sinh^2(\pi n)} \tag{5}$$
where:
$$ \sum_{n\geq 1}\frac{\pi^2}{n^2\sinh^2(\pi n)}=\sum_{n\geq 1}\left[\frac{1}{n^4}-\sum_{m\geq 1}\left(\frac{2}{n^2(n^2+m^2)}-\frac{4}{(m^2+n^2)^2}\right)\right]\tag{6}$$
equals:
$$ \zeta(4)+4(K \zeta(2)-\zeta(4))-\sum_{n\geq 1}\frac{-1+n \pi\coth(\pi n)}{n^4}=\frac{2 K \pi^2}{3}-\frac{11 \pi^4}{180}\tag{7}$$
and we get:

$$ g_2 = \color{red}{\frac{\Gamma\left(\frac{1}{4}\right)^8}{16\pi^2}}=\frac{4\pi^4}{\text{AGM}(1,\sqrt{2})^4}\tag{8} $$
  With just one step of the AGM mean we have that $\frac{1024 \pi ^4}{\left(1+2^{1/4}\right)^8}=189.0621\ldots$ is already a very accurate approximation of $g_2$.

Here it is the pseudocode for an improved evaluation:


*

*Set $a\leftarrow 1$ and $b\leftarrow \sqrt{2}$

*While $(b-a)$ is greater than the working precision do $c\leftarrow\frac{a+b}{2}, b\leftarrow\sqrt{ab}, a\leftarrow c$;

*Return $\frac{4\pi^4}{a^4}$.

A: According to Jack D'Aurizio
public static double g2() {
    double a= 1;
    double b=Math.sqrt(2);
    for (int i = 0; i <3; i++) {
         double  c= 0.5*(a+b);
           b= Math.sqrt(a*b);
           a=c;
    }
    double v= Math.PI/a;
    double v2=v*v;
    return 4*v2*v2;
}

gives 189.0727201.
Now I am looking for a corresponding algorithm for evaluating Weierstrass's elliptic function $\wp (z)$.
