add theta functions

acb_dirichlet.h

@@ -132,6 +132,7 @@ ulong acb_dirichlet_char_conductor(const acb_dirichlet_group_t G, const acb_diri
 void acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
 void acb_dirichlet_char_first_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
 
+ulong acb_dirichlet_ui_chi_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const acb_dirichlet_conrey_t x);
 ulong acb_dirichlet_ui_chi(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n);
 
 void acb_dirichlet_ui_vec_set_null(ulong *v, const acb_dirichlet_group_t G, slong nv);
@@ -140,12 +141,37 @@ void acb_dirichlet_ui_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G,
 void acb_dirichlet_ui_chi_vec_sieve(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
 void acb_dirichlet_ui_chi_vec(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv);
 
+/* precompute powers of a root of unity */
+typedef struct
+{
+    ulong order;
+    acb_ptr z;
+    ulong m;
+    ulong M;
+    acb_ptr Z;
+}
+acb_dirichlet_powers_struct;
+
+typedef acb_dirichlet_powers_struct acb_dirichlet_powers_t[1];
+
+void acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec);
+void acb_dirichlet_powers_clear(acb_dirichlet_powers_t t);
+void acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec);
+
 void acb_dirichlet_nth_root(acb_t res, ulong order, slong prec);
 
 void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec);
+void acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv, slong prec);
 
 void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec);
-ulong acb_dirichlet_theta_length(ulong q, double x, slong prec);
+
+ulong acb_dirichlet_theta_length_d(ulong q, double x, slong prec);
+ulong acb_dirichlet_theta_length(ulong q, const arb_t x, slong prec);
+void acb_dirichlet_chi_theta(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t x, slong prec);
+
+void acb_dirichlet_gauss_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
+void acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
+void acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec);
 
 #ifdef __cplusplus
 }

acb_dirichlet/arb_theta.c

@@ -0,0 +1,149 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+#include "acb_poly.h"
+
+/* x = Pi / q * t^2 */
+static void
+acb_dirichlet_arb_theta_argt(arb_t x, ulong q, const arb_t t, slong prec)
+{
+    arb_const_pi(x, prec);
+    arb_div_ui(x, x, q, prec);
+    arb_mul(x, x, t, prec);
+    arb_mul(x, x, t, prec);
+    arb_neg(x, x);
+    arb_exp(x, x, prec);
+}
+
+/* small order, multiply by chi at the end */
+void
+acb_dirichlet_arb_theta_smallorder(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
+{
+    slong k;
+    ulong order = z->order;
+    arb_t xk2, kxk2, dx, x2;
+    acb_ptr t;
+    arb_init(xk2);
+    arb_init(dx);
+    arb_init(x2);
+    arb_init(kxk2);
+
+    arb_set(dx, x);
+    arb_set(xk2, x);
+    arb_mul(x2, x, x, prec);
+
+    t = _acb_vec_init(order);
+    _acb_vec_zero(t, order);
+
+    arb_set(acb_realref(t + 0), xk2);
+
+    /* TODO: reduce precision at each step */
+    for (k = 2; k < len; k++)
+    {
+        arb_mul(dx, dx, x2, prec);
+        arb_mul(xk2, xk2, dx, prec);
+        if (a[k] != ACB_DIRICHLET_CHI_NULL)
+        {
+           if (parity)
+           {
+               arb_mul_si(kxk2, xk2, k, prec);
+               arb_add(acb_realref(t + a[k]), acb_realref(t + a[k]), kxk2, prec);
+           }
+           else
+           {
+               arb_add(acb_realref(t + a[k]), acb_realref(t + a[k]), xk2, prec);
+           }
+        }
+    }
+
+    /* now Hörner */
+    _acb_poly_evaluate(res, t, order, z->z + 1, prec);
+
+    _acb_vec_clear(t, order);
+    arb_clear(xk2);
+    arb_clear(x2);
+    arb_clear(dx);
+}
+
+void
+acb_dirichlet_arb_theta_naive(acb_t res, const arb_t x, int parity, const ulong *a, const acb_dirichlet_powers_t z, slong len, slong prec)
+{
+    slong k;
+    arb_t xk2, dx, x2;
+    acb_t zk;
+    arb_init(xk2);
+    arb_init(dx);
+    arb_init(x2);
+    acb_init(zk);
+
+    arb_set(dx, x);
+    arb_set(xk2, dx);
+    arb_mul(x2, dx, dx, prec);
+
+    acb_set_arb(res, xk2);
+
+    /* TODO: reduce prec */
+    for (k = 2; k < len; k++)
+    {
+        arb_mul(dx, dx, x2, prec);
+        arb_mul(xk2, xk2, dx, prec);
+        if (a[k] != ACB_DIRICHLET_CHI_NULL)
+        {
+            acb_dirichlet_power(zk, z, a[k], prec); 
+            if (parity)
+                acb_mul_si(zk, zk, k, prec);
+            acb_addmul_arb(res, zk, xk2, prec);
+        }
+    }
+
+    arb_clear(xk2);
+    arb_clear(x2);
+    arb_clear(dx);
+    acb_clear(zk);
+}
+
+void
+acb_dirichlet_chi_theta(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const arb_t t, slong prec)
+{
+    slong len;
+    ulong * a;
+    arb_t x;
+    acb_dirichlet_powers_t z;
+
+    len = acb_dirichlet_theta_length(G->q, t, prec);
+
+    a = flint_malloc(len * sizeof(ulong));
+    acb_dirichlet_ui_chi_vec(a, G, chi, len);
+    acb_dirichlet_powers_init(z, chi->order, len, prec);
+
+    arb_init(x);
+    acb_dirichlet_arb_theta_argt(x, G->q, t, prec);
+    acb_dirichlet_arb_theta_naive(res, x, chi->parity, a, z, len, prec);
+
+    arb_clear(x);
+    flint_free(a);
+    acb_dirichlet_powers_clear(z);
+}

acb_dirichlet/chi_vec.c

@@ -0,0 +1,51 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_chi_vec(acb_ptr v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv, slong prec)
+{
+    slong k;
+    ulong * a;
+    acb_dirichlet_powers_t t;
+
+    a = flint_malloc(nv * sizeof(ulong));
+    acb_dirichlet_ui_chi_vec(a, G, chi, nv);
+
+    acb_dirichlet_powers_init(t, chi->order, nv, prec);
+
+    acb_zero(v + 0);
+    for (k = 0; k < nv; k++)
+    {
+        if (a[k] != ACB_DIRICHLET_CHI_NULL)
+            acb_dirichlet_power(v + k, t, a[k], prec);
+        else
+            *(v + k) = *(v + 0);
+    }
+
+    acb_dirichlet_powers_clear(t);
+    flint_free(a);
+}

acb_dirichlet/gauss_sum.c

@@ -0,0 +1,51 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
+{
+    if (chi->order <= 2)
+    {
+        if (chi->parity)
+        {
+            arb_sqrt_ui(acb_imagref(res), G->q, prec);
+            arb_zero(acb_realref(res));
+        }
+        else
+        {
+            arb_sqrt_ui(acb_realref(res), G->q, prec);
+            arb_zero(acb_imagref(res));
+        }
+    }
+    else
+    {
+        if (acb_dirichlet_theta_length_d(G->q, 1, prec) > G->q)
+            acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
+        else
+            acb_dirichlet_gauss_sum_theta(res, G, chi, prec);
+    }
+}

acb_dirichlet/gauss_sum_naive.c

@@ -0,0 +1,45 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_gauss_sum_naive(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
+{
+    acb_t z;
+    acb_ptr v;
+
+    v = _acb_vec_init(G->q);
+    acb_dirichlet_chi_vec(v, G, chi, G->q, prec);
+
+    acb_init(z);
+    acb_dirichlet_nth_root(z, G->q, prec);
+
+    _acb_poly_evaluate(res, v, G->q, z, prec);
+
+    acb_clear(z);
+    _acb_vec_clear(v, G->q);
+}

acb_dirichlet/gauss_sum_theta.c

@@ -0,0 +1,65 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
+{
+    arb_t x;
+    acb_t eps;
+
+    arb_init(x);
+
+    if ((G->q == 300 && (chi->n == 271 || chi->n == 131))
+            || (G->q == 600 && (chi->n == 11 || chi->n == 91)))
+    {
+        /* or could use l'Hopital rule */
+        acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
+        return;
+    }
+
+    arb_one(x);
+    acb_dirichlet_chi_theta(res, G, chi, x, prec);
+    acb_init(eps);
+    acb_conj(eps, res);
+    acb_div(eps, res, eps, prec);
+
+    if (chi->parity)
+    {
+        arb_zero(acb_realref(res));
+        arb_sqrt_ui(acb_imagref(res), G->q, prec);
+    }
+    else
+    {
+        arb_zero(acb_imagref(res));
+        arb_sqrt_ui(acb_realref(res), G->q, prec);
+    }
+
+    acb_mul(res, res, eps, prec);
+
+    arb_clear(x);
+    acb_clear(eps);
+}

acb_dirichlet/l.c

@@ -12,31 +12,35 @@
 #include "acb_dirichlet.h"
 
 void
-acb_dirichlet_l(acb_t res, const acb_t s,
-    const acb_dirichlet_group_t G, ulong m, slong prec)
+acb_dirichlet_l_hurwitz(acb_t res, const acb_t s,
+    const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
 {
-    acb_t chi, t, u, a;
-    acb_dirichlet_conrey_t cm, cn;
+    ulong chin;
+    acb_t z, t, u, a;
+    acb_ptr xz;
+    acb_dirichlet_conrey_t cn;
 
-    acb_init(chi);
-    acb_dirichlet_conrey_init(cm, G);
     acb_dirichlet_conrey_init(cn, G);
+    acb_init(z);
     acb_init(t);
     acb_init(u);
     acb_init(a);
 
-    acb_dirichlet_conrey_log(cm, G, m);
     acb_dirichlet_conrey_one(cn, G);
     acb_zero(t);
 
+    acb_dirichlet_nth_root(z, chi->order, prec);
+    xz = _acb_vec_init(chi->order);
+    _acb_vec_set_powers(xz, z, chi->order, prec);
+
     while (1) {
-        /* todo: use n_dirichlet_chi and precomputed roots instead */
-        acb_dirichlet_pairing_conrey(chi, G, cm, cn, prec);
+        chin = acb_dirichlet_ui_chi_conrey(G, chi, cn);
 
         acb_set_ui(a, cn->n);
         acb_div_ui(a, a, G->q, prec);
         acb_hurwitz_zeta(u, s, a, prec);
-        acb_addmul(t, chi, u, prec);
+
+        acb_addmul(t, xz + chin, u, prec);
 
         if (acb_dirichlet_conrey_next(cn, G) == G->num)
           break;
@@ -47,10 +51,10 @@ acb_dirichlet_l(acb_t res, const acb_t s,
     acb_pow(u, u, a, prec);
     acb_mul(res, t, u, prec);
 
-    acb_dirichlet_conrey_clear(cm);
     acb_dirichlet_conrey_clear(cn);
 
-    acb_clear(chi);
+    _acb_vec_clear(xz, chi->order);
+    acb_clear(z);
     acb_clear(t);
     acb_clear(u);
     acb_clear(a);

acb_dirichlet/nth_root.c

@@ -25,8 +25,8 @@
 
 #include "acb_dirichlet.h"
 
-void
-acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
+static void
+_acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
 {
     fmpq_t t;
     fmpq_init(t);
@@ -34,3 +34,23 @@ acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
     arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
     fmpq_clear(t);
 }
+
+void
+acb_dirichlet_nth_root(acb_t res, ulong order, slong prec)
+{
+    switch (order)
+    {
+       case 1:
+           acb_one(res);
+           break;
+       case 2:
+           acb_set_si(res, -1);
+           break;
+       case 4:
+           acb_onei(res);
+           break;
+       default:
+           _acb_dirichlet_nth_root(res, order, prec);
+           break;
+    }
+}

acb_dirichlet/power.c

@@ -0,0 +1,51 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+#include "acb_poly.h"
+
+/* assume won't be modified */
+void
+acb_dirichlet_power(acb_t z, const acb_dirichlet_powers_t t, ulong n, slong prec)
+{
+    if (n < t->m)
+    {
+        /* acb_set(z, t->z + n); */
+        *z = *(t->z + n);
+    }
+    else
+    {
+        ulong q, r;
+        q = n / t->m;
+        r = n % t->m;
+        if (q >= t->M)
+        {
+            flint_printf("acb_dirichlet_power: power %wu not available "
+                    "in table of size %wu * %wu.", n, t->m, t->M);
+            abort();
+        }
+        acb_mul(z, t->Z + q, t->z + r, prec);
+    }
+}

acb_dirichlet/powers_clear.c

@@ -0,0 +1,35 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_powers_clear(acb_dirichlet_powers_t t)
+{
+    _acb_vec_clear(t->z, t->m);
+    if (t->M)
+        _acb_vec_clear(t->Z, t->M);
+}

acb_dirichlet/powers_init.c

@@ -0,0 +1,57 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_powers_init(acb_dirichlet_powers_t t, ulong order, slong num, slong prec)
+{
+    ulong m;
+    acb_t zeta;
+    t->order = order;
+    m = (num == 1) ? 1 : num * (prec / 64 + 1);
+    if (m > order)
+        m = order;
+    t->m = m;
+
+    acb_init(zeta);
+    acb_dirichlet_nth_root(zeta, order, prec);
+    t->z = _acb_vec_init(m);
+    _acb_vec_set_powers(t->z, zeta, m, prec);
+
+    if (order > m)
+    {
+        t->M = (order / m) + 1;
+        t->Z = _acb_vec_init(t->M);
+        acb_pow_ui(zeta, zeta, m, prec);
+        _acb_vec_set_powers(t->Z, zeta, t->M, prec);
+    }
+    else
+    {
+        t->M = 0;
+        t->Z = NULL;
+    }
+}

acb_dirichlet/test/t-gauss.c

@@ -0,0 +1,95 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+int main()
+{
+    slong prec = 128;
+    ulong q;
+
+    flint_printf("gauss....");
+    fflush(stdout);
+
+    /* check Gauss sums */
+
+    for (q = 3; q < 200; q ++)
+    {
+        acb_dirichlet_group_t G;
+        acb_dirichlet_conrey_t x;
+        acb_dirichlet_char_t chi;
+
+        acb_t s1, s2, s3;
+
+        if (q % 4 == 2)
+            /* no primitive character mod q */
+            continue;
+
+        acb_dirichlet_group_init(G, q);
+        acb_dirichlet_conrey_init(x, G);
+        acb_dirichlet_char_init(chi, G);
+
+        acb_init(s1);
+        acb_init(s2);
+        acb_init(s3);
+        acb_dirichlet_conrey_first_primitive(x, G);
+
+        while (1) {
+
+            acb_dirichlet_char_conrey(chi, G, x);
+
+            acb_dirichlet_gauss_sum_naive(s1, G, chi, prec);
+            acb_dirichlet_gauss_sum_theta(s2, G, chi, prec);
+            acb_dirichlet_gauss_sum(s3, G, chi, prec);
+
+            if (!acb_overlaps(s1, s2)
+                    || !acb_overlaps(s1, s3))
+            {
+                flint_printf("FAIL: G(chi_%wu(%wu))\n\n", q, chi->n);
+                flint_printf("\nnaive ", q, x->n);
+                acb_printd(s1, 25);
+                flint_printf("\ntheta ", q, x->n);
+                acb_printd(s2, 25);
+                flint_printf("\ndefault ", q, x->n);
+                acb_printd(s3, 25);
+                abort();
+            }
+
+            if (acb_dirichlet_conrey_next_primitive(x, G) == G->num)
+                break;
+        }
+        acb_clear(s1);
+        acb_clear(s2);
+        acb_clear(s3);
+
+        acb_dirichlet_group_clear(G);
+        acb_dirichlet_char_clear(chi);
+        acb_dirichlet_conrey_clear(x);
+    }
+
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/test/t-thetanull.c

@@ -66,7 +66,7 @@ int main()
         z = _acb_vec_init(G->expo);
         _acb_vec_set_powers(z, zeta, G->expo, prec);
 
-        nv = acb_dirichlet_theta_length(q, 1, prec);
+        nv = acb_dirichlet_theta_length_d(q, 1, prec);
         v = flint_malloc(nv * sizeof(ulong));
 
         arb_init(eq);

acb_dirichlet/theta_length.c

@@ -29,8 +29,22 @@
 #define LOG2 0.69314718055
 
 ulong
-acb_dirichlet_theta_length(ulong q, double x, slong prec)
+acb_dirichlet_theta_length_d(ulong q, double x, slong prec)
 {
     double a = PI / (double)q * x * x;
     return ceil(sqrt(((double)prec * LOG2 - log(2 * a)) / a));
 }
+
+ulong
+acb_dirichlet_theta_length(ulong q, const arb_t x, slong prec)
+{
+    double dx;
+    ulong len;
+    arf_t ax;
+    arf_init(ax);
+    arb_get_lbound_arf(ax, x, 53);
+    dx = arf_get_d(ax, ARF_RND_DOWN);
+    len = acb_dirichlet_theta_length_d(q, dx, prec);
+    arf_clear(ax);
+    return len;
+}

acb_dirichlet/ui_chi.c

@@ -34,13 +34,13 @@ acb_dirichlet_ui_chi(const acb_dirichlet_group_t G, const acb_dirichlet_char_t c
     }
     else
     {
-        ulong v = 0, k;
+        ulong v;
         acb_dirichlet_conrey_t x;
         acb_dirichlet_conrey_init(x, G);
         acb_dirichlet_conrey_log(x, G, n);
-        for (k = 0; k < G->num; k++)
-            v = (v + chi->expo[k] * x->log[k]) % chi->order;
-        acb_dirichlet_conrey_clear(x);
+
+        v = acb_dirichlet_ui_chi_conrey(G, chi, x);
+
         return v;
     }
 }

acb_dirichlet/ui_chi_conrey.c

@@ -0,0 +1,37 @@
+/*=============================================================================
+
+    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) 2016 Pascal Molin
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+ulong
+acb_dirichlet_ui_chi_conrey(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, const acb_dirichlet_conrey_t x)
+{
+        ulong v = 0, k;
+
+        for (k = 0; k < G->num; k++)
+            v = (v + chi->expo[k] * x->log[k]) % chi->order;
+
+        return v;
+}