add theta sum and j-invariant evaluation

acb_modular.h

@@ -133,6 +133,12 @@ void acb_modular_addseq_theta(long * exponents, long * aindex, long * bindex, lo
 
 void acb_modular_addseq_eta(long * exponents, long * aindex, long * bindex, long num);
 
+void acb_modular_theta_1234_sum(acb_t theta1, acb_t theta2,
+        acb_t theta3, acb_t theta4,
+    const acb_t w, int w_is_unit, const acb_t q, long prec);
+
+void acb_modular_j(acb_t z, const acb_t tau, long prec);
+
 #ifdef __cplusplus
 }
 #endif

acb_modular/j.c

@@ -0,0 +1,89 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_modular.h"
+
+void
+acb_modular_j(acb_t z, const acb_t tau, long prec)
+{
+    psl2z_t g;
+    arf_t one_minus_eps;
+    acb_t tau_prime, t1, t2, t3, t4, w, q;
+
+    psl2z_init(g);
+    arf_init(one_minus_eps);
+    acb_init(tau_prime);
+    acb_init(t1);
+    acb_init(t2);
+    acb_init(t3);
+    acb_init(t4);
+    acb_init(w);
+    acb_init(q);
+
+    arf_set_ui_2exp_si(one_minus_eps, 63, -6);
+    acb_modular_fundamental_domain_approx(tau_prime, g, tau,
+        one_minus_eps, prec);
+
+    acb_one(w);
+    acb_exp_pi_i(q, tau_prime, prec);
+    acb_modular_theta_1234_sum(t1, t2, t3, t4, w, 1, q, prec);
+
+    /* theta2 ^ 8 */
+    acb_mul(t2, t2, t2, prec);
+    acb_mul(t2, t2, t2, prec);
+    acb_mul(t2, t2, q, prec);
+    acb_mul(t2, t2, t2, prec);
+
+    /* theta3 ^ 8 */
+    acb_mul(t3, t3, t3, prec);
+    acb_mul(t3, t3, t3, prec);
+    acb_mul(t3, t3, t3, prec);
+
+    /* theta4 ^ 8 */
+    acb_mul(t4, t4, t4, prec);
+    acb_mul(t4, t4, t4, prec);
+    acb_mul(t4, t4, t4, prec);
+
+    acb_mul(z, t2, t3, prec);
+    acb_mul(z, z, t4, prec);
+
+    acb_add(t2, t2, t3, prec);
+    acb_add(t2, t2, t4, prec);
+    acb_cube(t2, t2, prec);
+
+    acb_div(z, t2, z, prec);
+    acb_mul_2exp_si(z, z, 5);
+
+    psl2z_clear(g);
+    arf_clear(one_minus_eps);
+    acb_clear(tau_prime);
+    acb_clear(t1);
+    acb_clear(t2);
+    acb_clear(t3);
+    acb_clear(t4);
+    acb_clear(w);
+    acb_clear(q);
+}
+

acb_modular/test/t-j.c

@@ -0,0 +1,145 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_modular.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("j....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* Test SL2Z invariance */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        acb_t tau1, tau2, z1, z2;
+        long e0, prec0, prec1, prec2;
+        psl2z_t g;
+
+        psl2z_init(g);
+        acb_init(tau1);
+        acb_init(tau2);
+        acb_init(z1);
+        acb_init(z2);
+
+        e0 = 1 + n_randint(state, 100);
+        prec0 = 2 + n_randint(state, 2000);
+        prec1 = 2 + n_randint(state, 2000);
+        prec2 = 2 + n_randint(state, 2000);
+
+        acb_randtest(tau1, state, prec0, e0);
+        acb_randtest(tau2, state, prec0, e0);
+        acb_randtest(z1, state, prec0, e0);
+        acb_randtest(z2, state, prec0, e0);
+
+        psl2z_randtest(g, state, 1 + n_randint(state, 200));
+
+        acb_modular_transform(tau2, g, tau1, prec0);
+
+        acb_modular_j(z1, tau1, prec1);
+        acb_modular_j(z2, tau2, prec2);
+
+        if (!acb_overlaps(z1, z2))
+        {
+            printf("FAIL (overlap)\n");
+            printf("tau1 = "); acb_print(tau1); printf("\n\n");
+            printf("tau2 = "); acb_print(tau2); printf("\n\n");
+            printf("z1 = "); acb_print(z1); printf("\n\n");
+            printf("z2 = "); acb_print(z2); printf("\n\n");
+            abort();
+        }
+
+        acb_modular_j(tau1, tau1, prec2);
+
+        if (!acb_overlaps(z1, tau1))
+        {
+            printf("FAIL (aliasing)\n");
+            printf("tau1 = "); acb_print(tau1); printf("\n\n");
+            printf("tau2 = "); acb_print(tau2); printf("\n\n");
+            printf("z1 = "); acb_print(z1); printf("\n\n");
+            printf("z2 = "); acb_print(z2); printf("\n\n");
+            abort();
+        }
+
+        acb_clear(tau1);
+        acb_clear(tau2);
+        acb_clear(z1);
+        acb_clear(z2);
+        psl2z_clear(g);
+    }
+
+    /* Test special values */
+    for (iter = 0; iter < 100; iter++)
+    {
+        acb_t tau, z;
+        long prec;
+
+        acb_init(tau);
+        acb_init(z);
+
+        prec = 2 + n_randint(state, 2000);
+
+        acb_randtest(z, state, prec, 10);
+
+        acb_onei(tau);
+        acb_modular_j(z, tau, prec);
+        acb_sub_ui(z, z, 1728, prec);
+
+        if (!acb_contains_zero(z))
+        {
+            printf("FAIL (value 1)\n");
+            printf("tau = "); acb_print(tau); printf("\n\n");
+            printf("z = "); acb_print(z); printf("\n\n");
+            abort();
+        }
+
+        acb_set_ui(tau, 2);
+        acb_div_ui(tau, tau, 3, prec);
+        acb_exp_pi_i(tau, tau, prec);
+
+        acb_modular_j(z, tau, prec);
+
+        if (!acb_contains_zero(z))
+        {
+            printf("FAIL (value 2)\n");
+            printf("tau = "); acb_print(tau); printf("\n\n");
+            printf("z = "); acb_print(z); printf("\n\n");
+            abort();
+        }
+
+        acb_clear(tau);
+        acb_clear(z);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_modular/test/t-theta_1234_sum.c

@@ -0,0 +1,123 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_modular.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("theta_1234_sum....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* Very weak test, just testing the error bounds and not
+       that we compute the right functions */
+    for (iter = 0; iter < 10000; iter++)
+    {
+        acb_t t1a, t1b, t2a, t2b, t3a, t3b, t4a, t4b, w, q;
+        int w_is_unit;
+        long prec0, e0, prec1, prec2;
+
+        acb_init(t1a);
+        acb_init(t1b);
+        acb_init(t2a);
+        acb_init(t2b);
+        acb_init(t3a);
+        acb_init(t3b);
+        acb_init(t4a);
+        acb_init(t4b);
+        acb_init(w);
+        acb_init(q);
+
+        e0 = 1 + n_randint(state, 100);
+        prec0 = 2 + n_randint(state, 3000);
+        prec1 = 2 + n_randint(state, 3000);
+        prec2 = 2 + n_randint(state, 3000);
+
+        if (n_randint(state, 2))
+        {
+            arb_randtest(acb_realref(q), state, prec0, e0);
+            arb_zero(acb_imagref(q));
+            acb_exp_pi_i(w, q, prec0);
+            w_is_unit = n_randint(state, 2);
+        }
+        else
+        {
+            acb_randtest(w, state, prec0, e0);
+            w_is_unit = 0;
+        }
+
+        acb_randtest(q, state, prec0, e0);
+
+        acb_randtest(t1a, state, prec0, e0);
+        acb_randtest(t1b, state, prec0, e0);
+        acb_randtest(t2a, state, prec0, e0);
+        acb_randtest(t2b, state, prec0, e0);
+        acb_randtest(t3a, state, prec0, e0);
+        acb_randtest(t3b, state, prec0, e0);
+        acb_randtest(t4a, state, prec0, e0);
+        acb_randtest(t4b, state, prec0, e0);
+
+        acb_modular_theta_1234_sum(t1a, t2a, t3a, t4a, w, w_is_unit, q, prec1);
+        acb_modular_theta_1234_sum(t1b, t2b, t3b, t4b, w, w_is_unit & n_randint(state, 2), q, prec2);
+
+        if (!acb_overlaps(t1a, t1b) || !acb_overlaps(t2a, t2b)
+            || !acb_overlaps(t3a, t3b) || !acb_overlaps(t4a, t4b))
+        {
+            printf("FAIL (overlap)\n");
+            printf("q = "); acb_print(q); printf("\n\n");
+            printf("w = "); acb_print(w); printf("\n\n");
+            printf("t1a = "); acb_print(t1a); printf("\n\n");
+            printf("t1b = "); acb_print(t1b); printf("\n\n");
+            printf("t2a = "); acb_print(t2a); printf("\n\n");
+            printf("t2b = "); acb_print(t2b); printf("\n\n");
+            printf("t3a = "); acb_print(t3a); printf("\n\n");
+            printf("t3b = "); acb_print(t3b); printf("\n\n");
+            printf("t4a = "); acb_print(t4a); printf("\n\n");
+            printf("t4b = "); acb_print(t4b); printf("\n\n");
+            abort();
+        }
+
+        acb_clear(t1a);
+        acb_clear(t1b);
+        acb_clear(t2a);
+        acb_clear(t2b);
+        acb_clear(t3a);
+        acb_clear(t3b);
+        acb_clear(t4a);
+        acb_clear(t4b);
+        acb_clear(w);
+        acb_clear(q);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_modular/theta_1234_sum.c

@@ -0,0 +1,358 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_modular.h"
+
+static void
+acb_mul_approx(acb_t z, acb_t tmp1, acb_t tmp2, const acb_t x, const acb_t y, long wprec, long prec)
+{
+    if (prec <= 1024)
+    {
+        acb_mul(z, x, y, wprec);
+    }
+    else if (x == y)
+    {
+        acb_set_round(tmp1, x, wprec);
+        acb_mul(z, tmp1, tmp1, wprec);
+    }
+    else
+    {
+        acb_set_round(tmp1, x, wprec);
+        acb_set_round(tmp2, y, wprec);
+        acb_mul(z, tmp1, tmp2, wprec);
+    }
+}
+
+void mag_sub_lower(mag_t z, const mag_t x, const mag_t y);
+
+double
+mag_get_log2_d_approx(const mag_t x)
+{
+    if (mag_is_zero(x))
+    {
+        return COEFF_MIN;
+    }
+    else if (mag_is_inf(x))
+    {
+        return COEFF_MAX;
+    }
+    else if (COEFF_IS_MPZ(MAG_EXP(x)))
+    {
+        if (fmpz_sgn(MAG_EXPREF(x)) < 0)
+            return COEFF_MIN;
+        else
+            return COEFF_MAX;
+    }
+    else
+    {
+        long e = MAG_EXP(x);
+
+        if (e < -20 || e > 20)
+            return e;
+        else
+            return e + 1.4426950408889634074 *
+                mag_d_log_upper_bound(MAG_MAN(x) * ldexp(1.0, -MAG_BITS));
+    }
+}
+
+void
+acb_modular_theta_1234_sum(acb_t theta1, acb_t theta2,
+        acb_t theta3, acb_t theta4,
+    const acb_t w, int w_is_unit, const acb_t q, long prec)
+{
+    mag_t err, qmag, wmag, vmag;
+    double log2q_approx, log2w_approx, log2term_approx;
+    long e, e1, e2, k, k1, k2, N, WN, term_prec;
+    long *exponents, *aindex, *bindex;
+    acb_ptr qpow, wpow, vpow;
+    acb_t tmp1, tmp2, v;
+    int q_is_real, w_is_one;
+
+    q_is_real = arb_is_zero(acb_imagref(q));
+    w_is_one = acb_is_one(w);
+
+    acb_init(tmp1);
+    acb_init(tmp2);
+    acb_init(v);
+    mag_init(err);
+
+    mag_init(qmag);
+    mag_init(wmag);
+    mag_init(vmag);
+
+    if (w_is_one)
+        acb_one(v);
+    else
+        acb_inv(v, w, prec);
+
+    acb_get_mag(qmag, q);
+    log2q_approx = mag_get_log2_d_approx(qmag);
+
+    if (w_is_unit)
+    {
+        mag_one(wmag);
+        mag_one(vmag);
+        log2w_approx = 0.0;
+    }
+    else
+    {
+        acb_get_mag(wmag, w);
+        acb_get_mag(vmag, v);
+        mag_max(wmag, wmag, vmag);
+        log2w_approx = mag_get_log2_d_approx(wmag);
+    }
+
+    if (log2q_approx >= 0.0)
+    {
+        N = 1;
+        mag_inf(err);
+    }
+    else  /* Pick N and compute error bound */
+    {
+        mag_t den;
+        mag_init(den);
+
+        N = 1;
+        while (0.05 * N * N < prec)
+        {
+            log2term_approx = log2q_approx * ((N+2)*(N+2)/4) + (N+2)*log2w_approx;
+            if (log2term_approx < -prec - 2)
+                break;
+            N++;
+        }
+
+        if (w_is_unit)
+        {
+            mag_one(den);
+            mag_sub_lower(den, den, qmag);  /* 1 - |q| is good enough */
+        }
+        else  /* denominator: 1 - |q|^(floor((N+1)/2)+1) * max(|w|,1/|w|) */
+        {
+            mag_pow_ui(err, qmag, (N + 1) / 2 + 1);
+            mag_mul(err, err, wmag);
+            mag_one(den);
+            mag_sub_lower(den, den, err);
+        }
+
+        /* no convergence */
+        if (mag_is_zero(den))
+        {
+            N = 1;
+            mag_inf(err);
+        }
+        else if (w_is_unit)
+        {
+            mag_pow_ui(err, qmag, ((N + 2) * (N + 2)) / 4);
+            mag_div(err, err, den);
+            mag_mul_2exp_si(err, err, 1);
+        }
+        else
+        {
+            mag_pow_ui(err, qmag, ((N + 2) * (N + 2)) / 4);
+            mag_pow_ui(vmag, wmag, N + 2);
+            mag_mul(err, err, vmag);
+            mag_div(err, err, den);
+            mag_mul_2exp_si(err, err, 1);
+        }
+
+        mag_clear(den);
+    }
+
+    exponents = flint_malloc(sizeof(long) * 3 * N);
+    aindex = exponents + N;
+    bindex = aindex + N;
+
+    qpow = _acb_vec_init(N);
+
+    acb_modular_addseq_theta(exponents, aindex, bindex, N);
+    acb_set_round(qpow + 0, q, prec);
+
+    acb_zero(theta1);
+    acb_zero(theta2);
+    acb_zero(theta3);
+    acb_zero(theta4);
+
+    WN = (N + 3) / 2;
+
+    /* compute powers of w^2 and v = 1/w^2 */
+    if (!w_is_one)
+    {
+        wpow = _acb_vec_init(WN);
+        vpow = _acb_vec_init(WN + 1);
+
+        acb_mul(tmp1, w, w, prec);
+        acb_mul(tmp2, v, v, prec);
+
+        _acb_vec_set_powers(wpow, tmp1, WN, prec);
+        _acb_vec_set_powers(vpow, tmp2, WN + 1, prec);
+    }
+    else
+    {
+        wpow = vpow = NULL;
+    }
+
+    for (k = 0; k < N; k++)
+    {
+        e = exponents[k];
+
+        log2term_approx = e * log2q_approx + (k+2) * log2w_approx;
+        term_prec = FLINT_MIN(FLINT_MAX(prec + log2term_approx + 16.0, 16.0), prec);
+
+        if (k > 0)
+        {
+            k1 = aindex[k];
+            k2 = bindex[k];
+            e1 = exponents[k1];
+            e2 = exponents[k2];
+
+            if (e == e1 + e2)
+            {
+                acb_mul_approx(qpow + k, tmp1, tmp2, qpow + k1, qpow + k2, term_prec, prec);
+            }
+            else if (e == 2 * e1 + e2)
+            {
+                acb_mul_approx(qpow + k, tmp1, tmp2, qpow + k1, qpow + k1, term_prec, prec);
+                acb_mul_approx(qpow + k, tmp1, tmp2, qpow + k, qpow + k2, term_prec, prec);
+            }
+            else
+            {
+                printf("exponent not in addition sequence!\n");
+                abort();
+            }
+        }
+
+        if (w_is_one)
+        {
+            if (k % 2 == 0)
+            {
+                acb_add(theta3, theta3, qpow + k, prec);
+
+                if (k % 4 == 0)
+                    acb_sub(theta4, theta4, qpow + k, prec);
+                else
+                    acb_add(theta4, theta4, qpow + k, prec);
+            }
+            else
+            {
+                acb_add(theta2, theta2, qpow + k, prec);
+            }
+        }
+        else
+        {
+            if (k % 2 == 0)
+            {
+                acb_add(tmp1, wpow + k / 2 + 1, vpow + k / 2 + 1, term_prec);
+                acb_mul(tmp1, qpow + k, tmp1, term_prec);
+
+                acb_add(theta3, theta3, tmp1, prec);
+
+                if (k % 4 == 0)
+                    acb_sub(theta4, theta4, tmp1, prec);
+                else
+                    acb_add(theta4, theta4, tmp1, prec);
+            }
+            else
+            {
+                if (k / 2 + 1 > WN - 1)
+                    abort();
+                if (k / 2 + 2 > WN + 1 - 1)
+                    abort();
+
+                acb_add(tmp1, wpow + k / 2 + 1, vpow + k / 2 + 2, term_prec);
+                acb_mul(tmp1, qpow + k, tmp1, term_prec);
+                acb_add(theta2, theta2, tmp1, prec);
+
+                acb_sub(tmp1, wpow + k / 2 + 1, vpow + k / 2 + 2, term_prec);
+                acb_mul(tmp1, qpow + k, tmp1, term_prec);
+                if (k % 4 == 1)
+                    acb_sub(theta1, theta1, tmp1, prec);
+                else
+                    acb_add(theta1, theta1, tmp1, prec);
+            }
+        }
+    }
+
+    if (w_is_one)
+    {
+        acb_mul_2exp_si(theta2, theta2, 1);
+        acb_mul_2exp_si(theta3, theta3, 1);
+        acb_mul_2exp_si(theta4, theta4, 1);
+
+        acb_add_ui(theta2, theta2, 2, prec);
+        acb_add_ui(theta3, theta3, 1, prec);
+        acb_add_ui(theta4, theta4, 1, prec);
+    }
+    else
+    {
+        /* w * [(1 - w^-2) + series] */
+        acb_sub(theta1, theta1, vpow + 1, prec);
+        acb_mul(theta1, theta1, w, prec);
+        acb_add(theta1, theta1, w, prec);
+
+        /* multiply by -i */
+        acb_mul_onei(theta1, theta1);
+        acb_neg(theta1, theta1);
+
+        /* w * [(1 + w^-2) + series] */
+        acb_add(theta2, theta2, vpow + 1, prec);
+        acb_mul(theta2, theta2, w, prec);
+        acb_add(theta2, theta2, w, prec);
+
+        acb_add_ui(theta3, theta3, 1, prec);
+        acb_add_ui(theta4, theta4, 1, prec);
+    }
+
+    if (q_is_real && w_is_unit)    /* result must be real */
+    {
+        arb_add_error_mag(acb_realref(theta1), err);
+        arb_add_error_mag(acb_realref(theta2), err);
+        arb_add_error_mag(acb_realref(theta3), err);
+        arb_add_error_mag(acb_realref(theta4), err);
+    }
+    else
+    {
+        acb_add_error_mag(theta1, err);
+        acb_add_error_mag(theta2, err);
+        acb_add_error_mag(theta3, err);
+        acb_add_error_mag(theta4, err);
+    }
+
+    if (!w_is_one)
+    {
+        _acb_vec_clear(wpow, WN);
+        _acb_vec_clear(vpow, WN + 1);
+    }
+
+    flint_free(exponents);
+    _acb_vec_clear(qpow, N);
+    acb_clear(tmp1);
+    acb_clear(tmp2);
+    acb_clear(v);
+    mag_clear(err);
+    mag_clear(qmag);
+    mag_clear(wmag);
+    mag_clear(vmag);
+}
+

arb/sin_cos_pi.c

@@ -32,6 +32,12 @@ arb_sin_pi(arb_t y, const arb_t x, long prec)
     arb_t u;
     fmpz_t v;
 
+    if (!arb_is_finite(x))
+    {
+        arb_indeterminate(y);
+        return;
+    }
+
     if (arf_cmpabs_2exp_si(arb_midref(x), FLINT_MAX(65536, (4*prec))) > 0)
     {
         arf_zero(arb_midref(y));
@@ -81,6 +87,12 @@ arb_cos_pi(arb_t y, const arb_t x, long prec)
     arb_t u;
     fmpz_t v;
 
+    if (!arb_is_finite(x))
+    {
+        arb_indeterminate(y);
+        return;
+    }
+
     if (arf_cmpabs_2exp_si(arb_midref(x), FLINT_MAX(65536, (4*prec))) > 0)
     {
         arf_zero(arb_midref(y));
@@ -130,6 +142,13 @@ arb_sin_cos_pi(arb_t s, arb_t c, const arb_t x, long prec)
     arb_t u;
     fmpz_t v;
 
+    if (!arb_is_finite(x))
+    {
+        arb_indeterminate(s);
+        arb_indeterminate(c);
+        return;
+    }
+
     if (arf_cmpabs_2exp_si(arb_midref(x), FLINT_MAX(65536, (4*prec))) > 0)
     {
         arf_zero(arb_midref(s));

doc/source/acb_modular.rst

@@ -130,29 +130,89 @@ Modular transformations
     `|\operatorname{Re}(z)| \le 1/2 + \varepsilon` where `\varepsilon` is
     specified by *tol*. Returns zero if this is false or cannot be determined.
 
-
-The Dedekind eta function
--------------------------------------------------------------------------------
-
-To be done
-
-.. function:: void acb_modular_addseq_eta(long * exponents, long * aindex, long * bindex, long num)
-
-    Constructs an addition sequence for the first *num* generalized pentagonal
-    numbers (excluding zero), i.e. 1, 2, 5, 7, 12, 15, 22, 26, 35, 40 etc.
-
 Jacobi theta functions
 -------------------------------------------------------------------------------
 
-To be done
-
 .. function:: void acb_modular_addseq_theta(long * exponents, long * aindex, long * bindex, long num)
 
     Constructs an addition sequence for the first *num* squares and triangular
     numbers interleaved (excluding zero), i.e. 1, 2, 4, 6, 9, 12, 16, 20, 25, 30 etc.
 
-Eisenstein series
+.. function:: void acb_modular_theta_1234_sum(acb_t theta1, acb_t theta2, acb_t theta3, acb_t theta4, const acb_t w, int w_is_unit, const acb_t q, long prec)
+
+    Simultaneously evaluates
+
+    .. math ::
+
+        q^{-1/4} \theta_1 = 2 \sum_{n=0}^{\infty} (-1)^n q^{n(n+1)} \sin((2n+1) \pi z)
+
+        q^{-1/4} \theta_2 = 2 \sum_{n=0}^{\infty} q^{n(n+1)} \cos((2n+1) \pi z)
+
+        \theta_3 = 1 + 2 \sum_{n=1}^{\infty} q^{n^2} \cos(2n \pi z)
+
+        \theta_4 = 1 + 2 \sum_{n=1}^{\infty} (-1)^n q^{n^2} \cos(2n \pi z)
+
+    given `w = \exp(\pi i z)` and `q = \exp(\pi i \tau)`.
+    If *w_is_unit* is nonzero, *w* is assumed to lie on the unit circle,
+    i.e. *z* is assumed to be real.
+
+    Note that the factor `q^{1/4}` is removed from `\theta_1` and `\theta_2`.
+    To get the true theta function values, the user has to multiply
+    this factor back. This convention avoids an unnecessary root extraction,
+    since the user can compute `q^{1/4} = \exp(\pi i \tau / 4)` followed by
+    `q = (q^{1/4})^4`, and in many cases when computing products or quotients
+    of theta functions, the factor `q^{1/4}` can be eliminated entirely.
+
+    This function is intended for `|q| \ll 1`. It can be called with any
+    `q`, but will return useless intervals if convergence is not rapid.
+    For general evaluation of theta functions, the user should only call
+    this function after applying a suitable modular transformation.
+
+    We consider the sums together, alternatingly updating `(\theta_1, \theta_2)`
+    or `(\theta_3, \theta_4)`. For `k = 0, 1, 2, \ldots`, the powers of `q`
+    are `\lfloor (k+2)^2 / 4 \rfloor = 1, 2, 4, 6, 9` etc. and the powers of `w` are
+    `\pm (k+2) = \pm 2, \pm 3, \pm 4, \ldots` etc.
+    For some integer `N \ge 1`, the summation is stopped just before term
+    `k = N`. The error can then be bounded as
+
+    .. math ::
+
+        \frac{2 |q|^E \max(|w|,|w^{-1}|)^{N+2}}{1 - |q|^{\lfloor (N+1)/2 \rfloor + 1} \max(|w|,|w^{-1}|)}
+
+    where `E = \lfloor (N+2)^2 / 4 \rfloor`, assuming that the denominator
+    is positive.
+    This is simply the bound for a geometric series, with the leading
+    factor 2 coming from the fact that we sum both negative and positive
+    powers of `w`.
+
+    To actually evaluate the series, when `w \ne 1`, we write the even
+    cosine terms as `w^{2n} + w^{-2n}`, the odd cosine terms as
+    `w (w^{2n} + w^{-2n-2})`, and the sine terms as `w (w^{2n} - w^{-2n-2})`.
+    This way we only need even powers of `w` and `w^{-1}`.
+    The implementation is not yet optimized for real `z`, in which case
+    further work can be saved.
+
+    This function does not permit aliasing between input and output
+    arguments.
+
+The Dedekind eta function
 -------------------------------------------------------------------------------
 
-To be done
+.. function:: void acb_modular_addseq_eta(long * exponents, long * aindex, long * bindex, long num)
+
+    Constructs an addition sequence for the first *num* generalized pentagonal
+    numbers (excluding zero), i.e. 1, 2, 5, 7, 12, 15, 22, 26, 35, 40 etc.
+
+Modular forms
+-------------------------------------------------------------------------------
+
+.. function:: void acb_modular_j(acb_t r, const acb_t tau, long prec)
+
+    Computes Klein's j-invariant `j(\tau)` given `\tau` in the upper
+    half-plane. The function is normalized so that `j(i) = 1728`.
+    We first move `\tau` to the fundamental domain, which does not change
+    the value of the function. Then we use the formula
+    `j(\tau) = 32 (\theta_2^8 + \theta_3^8 + \theta_4^8)^3 / (\theta_2 \theta_3 \theta_4)^8`
+    where `\theta_k` is the respective theta constant evaluated at `\tau`.
+