# Computing Eisenstein series and modular discriminant

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);
}
}
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);
}
_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.

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).$

• Yesterday I was trying to compute such things on pari/gp (what I really need is $\sum_{n\ge 1, (n,k)=1} a(n) e^{2i \pi n z}$ with $a_n$ the coefs of some modular forms, in particular $a(n^2) = \chi(n)$ for arbitrary Dirichlet character $\chi$). What software do you use, and do you have some codes ? Oct 17, 2017 at 17:32
• @reuns I use PARI/GP and Mathematica depending on what is more convenient. Of course I have codes, Some of it is on my website and lots in OEIS sequence entries. Oct 17, 2017 at 18:00
• $q = e^{2i \pi z}$. The $q$-series for modular forms works well for $|q| \le r$ but not for $|q| \to 1$, unfortunately what I need is to evaluate modular forms near $\Im(z) = 0, |q|=1$ (this is where the Riemann hypothesis is encoded, where infinitely many cusps meet). Thus I need to find $ad-bc=1$ maximizing $\Im( \frac{az+b}{cz+d}) = \frac{\Im(z)}{|cz+d|^2}$ Oct 17, 2017 at 18:22
• @reuns Ask another question. Not here. Oct 17, 2017 at 18:24
• The answer to my question is to use the LLL algorithm to minimize $|cz+d|^2$ ie. maximize $\Im(\gamma z)$ so that $\gamma z$ is in the $|\tau|> 1, \Re(\tau) \in [-1,1]$ fundamental domain (a matlab implemention of LLL there). Do you have some sort of "modular forms toolbox" in pari/gp or anything else ? Oct 17, 2017 at 19:19