Computing Eisenstein series and modular discriminant

Question

I begin to self learn elliptic functions and modular forms by coding (from scratch) the special functions such as $E_6(\tau), E_6(\tau)$ (Eisenstein series), $g_2(\tau),g_3(\tau)$ (invariants in Weierstrass elliptic functions), and the discriminant $\Delta$.

It seems that my codes for $g_2(\tau),g_3(\tau)$ are correct, the Java codes for $g_3(\tau)$ read

public static double[] g3(double[] tao, int max) { // Eisenstein series
    double[] sum = new double[2];
    for (int i = -max; i <= max; i++) {
        for (int j = -max; j <= max; j++) {
            if (0 == i && 0 == j)
                continue;
            double[] omega = M.add(new double[] { j, 0 }, M.scale(tao, i));
            double[] c2 = M.square(omega);
            double[] c4 = M.square(c2);
            double[] c6 = M.mul(c2, c4);
            double[] v = M.inverse(c6);
            _add(sum, v);
        }
    }
    return scale(sum, 140);
}

which implements $$g_3(\tau)= 140\sum_{ (m,n) \neq 0} (m+n\tau)^{-6}$$

The codes succeeded in generating an image similar to https://en.wikipedia.org/wiki/Weierstrass%27s_elliptic_functions#/media/File:Gee_three_real.jpeg

However, my codes for $\Delta= g_2^3- 27g_3^2$

    double[] g2=g2(tao, far) ;
    double[] g3=g3(tao, far) ;
    double[] g2_3 = cubic(g2);
    double[] g3_2 = square(g3);
    return sub(g2_3, scale(g3_2, 27));

seem problematic as the resultant image differs from https://en.wikipedia.org/wiki/Weierstrass%27s_elliptic_functions#/media/File:Discriminant_real_part.jpeg.

I try to double-check whether my calculation of $g_2, g_3$ is correct by calculating $$E_4(\tau)= 1+240 \sum_{n=1} \frac{n^3 q^n}{1-q^n} $$ $$E_6(\tau)= 1-504 \sum_{n=1} \frac{n^5 q^n}{1-q^n} $$ where $q= e^{2 \pi i \tau}$, since $$g_2= \frac{4}{3} \pi^4 E_4$$ $$g_3= \frac{8}{27} \pi^6 E_6$$

Codes for $E_6$ are

public static double[] E6(double[] q, int far) { // Eisenstein series by q
    double[] sum6 = new double[2];
    for (int i = 1; i < far; i++) {
        double[] qn = power(q, i);
        double[] fm = { 1 - qn[0], -qn[1] };
        double[] fz = scale(qn, power(i, 5)); // n^5 q^n
        double[] vb = divide(fz, fm);
        _add(sum6, vb);
    }
    _scale(sum6, -504);
    double[] E6 = { 1 + sum6[0], sum6[1] };
    return E6;
}

However the calculated values of $E_6$ is far from consistent with that of $g_3$. Something must be wrong. I am not ready to blame computer's inaccuracy of float number calculation (rounding errors for example), because the inconsistency is significant.

ps. I am equipped with a set of elementary methods such as

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 void _add(double[] a,  double[] b) {
    for(int i=0;i<a.length;i++)
      a[i]+=b[i];
}
public static double[] divide(double[] c1, double[] c2) {
    double a = c1[0];
    double b = c1[1];
    double c = c2[0];
    double d = c2[1];
    double x = (a * c + b * d) / (c * c + d * d);
    double y = (b * c - a * d) / (c * c + d * d);
    return new double[] { x, y };
}
public static double[] mul(double[] c1, double[] c2) {
    double a = c1[0];
    double b = c1[1];
    double c = c2[0];
    double d = c2[1];
    return new double[] { a * c - b * d, b * c + a * d };
}

which are unlikely to be wrong.

Answer 1

A big problem with the calculation of sums over the nonzero periods is that they need a lot of terms, so if you did not use enough terms you will get relatively inaccurate result. On the other hand, sum of q-series only needs a handful of terms to get accurate results. For example, with $q=.1^6$ using $5$ terms of the $q$-series gives $\;g_2=129.909959,\; g_3=284.712485,\;$ but using $-100\leq n,m \leq 100$ gives $\;g_2=129.909727,\; g_3=284.712485.\;$ Notice t hat $g_3$ is good in both, but not for $g_2$. In case you are interested, my code using PARI/GP is:

default(realprecision, 9); /* decimal place precision */
q = .1^6; N = 5;
E4 = 1 + 240*sum(n=1, N, n^3*q^n/(1-q^n));
E6 = 1 - 504*sum(n=1, N, n^5*q^n/(1-q^n));
g2E = 4/3*Pi^4*E4; g3E = 8/27*Pi^6*E6; print([g2E, g3E]);
tau = log(q)/(2*Pi*I); M = 100;
S4 = sum(n=-M, M, sum(m=-M, M, if(n==0&&m==0, 0, (m+n*tau)^-4 )));
S6 = sum(n=-M, M, sum(m=-M, M, if(n==0&&m==0, 0, (m+n*tau)^-6 )));
g2S = 60*S4; g3S = 140*S6; print([g2S, g3S]);

The calculation of the discriminant $\Delta$ is a more delicate matter. Computing the Eisenstein series requires many more terms if $|q|$ is not small. The formula $\Delta= g_2^3- 27g_3^2$ is subject to catastrophic cancellation, if $|q|$ is not small. There are alternative formulas that can be used. For example, $\Delta = (2\pi)^{12}q\prod_{n=1}^\infty (1-q^n)^{24}.$ Also, because $E_4,E_6,\Delta$ are modular forms, you can use transformation formulas to make $|q|$ small. For example, $\Delta(-1/\tau)=(\tau/i)^{12}\Delta(\tau).$