fix bug on char_next + add log precomputations

acb_dirichlet.h

@@ -23,6 +23,7 @@ extern "C" {
 typedef struct
 {
     ulong q;                /* modulus */
+    nmod_t mod;             /* modulus with precomputed inverse */
     ulong q_even;           /* even part of modulus */
     ulong q_odd;            /* odd part of modulus */
     ulong phi_q;            /* phi(q) = group size */
@@ -35,6 +36,7 @@ typedef struct
     ulong * generators;     /* generator for each prime p[k] lifted mod q */
     ulong * phi;            /* phi(k) = phi(p[k]^e[k])       */
     ulong * PHI;            /* PHI(k) = expo / phi(k)        */
+    dlog_precomp_t * dlog;  /* precomputed data for discrete log mod p^e */
 }
 acb_dirichlet_group_struct;
 
@@ -101,6 +103,9 @@ void acb_dirichlet_char_conrey(acb_dirichlet_char_t chi, const acb_dirichlet_gro
 void acb_dirichlet_char_normalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
 void acb_dirichlet_char_denormalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
 
+ulong acb_dirichlet_char_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi);
+int acb_dirichlet_char_is_odd(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi);
+
 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);
 
@@ -115,6 +120,10 @@ void acb_dirichlet_chi_vec(ulong *v, ulong nv, const acb_dirichlet_group_t G, co
 
 void acb_dirichlet_char_vec(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec);
 
+
+void acb_dirichlet_zeta(acb_t res, ulong order, slong prec);
+void acb_dirichlet_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

acb_dirichlet/arb_quadratic_powers.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_arb_quadratic_powers(arb_ptr v, slong nv, const arb_t x, slong prec)
+{
+    slong i;
+    arb_t dx, x2;
+    arb_init(dx);
+    arb_init(x2);
+    arb_set(dx, x);
+    arb_mul(x2, x, x, prec);
+    for (i = 0; i < nv; i++)
+    {
+        if (i == 0)
+            arb_one(v + i);
+        else if (i == 1)
+            arb_set_round(v + i, x, prec);
+        else
+        {
+            arb_mul(dx, dx, x2, prec);
+            arb_mul(v + i, v + i - 1, dx, prec);
+        }
+    }
+    arb_clear(dx);
+    arb_clear(x2);
+}

acb_dirichlet/char_conductor.c

@@ -0,0 +1,56 @@
+/*=============================================================================
+
+    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_char_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
+{
+    int k, f;
+    ulong cond = 1;
+    if (G->neven >= 1 && chi->expo[0] == 1)
+        cond = 4;
+    if (G->neven == 2 && chi->expo[1])
+    {
+        ulong l2 = chi->expo[1];
+        f = n_remove(&l2, 2);
+        cond = G->primepowers[1] >> f;
+    }
+    for (k = G->neven; k < G->neven; k++)
+    {
+        if (chi->expo[k])
+        {
+            ulong p, lp;
+            p = G->primes[k];
+            lp = chi->expo[k];
+            f = n_remove(&lp, p);
+            if (f)
+                cond *= n_pow(p, G->exponents[k] - f);
+            else
+                cond *= G->primepowers[k];
+        }
+    }
+    return cond;
+}

acb_dirichlet/char_is_odd.c

@@ -0,0 +1,40 @@
+/*=============================================================================
+
+    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
+acb_dirichlet_char_is_odd(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
+{
+    slong k, odd = 0;
+    for (k = 0; k < G->num; k++)
+    {
+        if (k == 1 && G->neven == 2)
+            continue;
+        if (chi->expo[k] % 2 == 1)
+            odd = 1 - odd;
+    }
+    return odd;
+}

acb_dirichlet/char_next.c

@@ -35,8 +35,7 @@ acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
   /* update index */
   for (k=0; k < G->num ; k++)
   {
-    /* chi->n = n_mulmod(chi->n, G->generators[k], G->q); */
-    chi->n = chi->n * G->generators[k] % G->q;
+    chi->n = nmod_mul(chi->n, G->generators[k], G->mod);
     chi->expo[k] += G->PHI[k];
     if (chi->expo[k] < G->expo)
       break;

acb_dirichlet/char_next_primitive.c

@@ -36,8 +36,7 @@ acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_
   k = 0;
   if (G->neven == 2)
   {
-      /* chi->n = n_mulmod(chi->n, G->generators[0], G->q); */
-      chi->n = chi->n * G->generators[0] % G->q;
+      chi->n = nmod_mul(chi->n, G->generators[0], G->mod);
       if (++chi->expo[0] < G->expo)
           return 0;
       chi->expo[0] = 0;
@@ -45,13 +44,11 @@ acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_
   }
   for (; k < G->num ; k++)
   {
-      /* chi->n = n_mulmod(chi->n, G->generators[k], G->q); */
-      chi->n = chi->n * G->generators[k] % G->q;
+      chi->n = nmod_mul(chi->n, G->generators[k], G->mod);
       chi->expo[k] += G->PHI[k];
       if (chi->expo[k] % G->primes[k] == 0)
       {
-          /* chi->n = n_mulmod(chi->n, G->generators[k], G->q); */
-          chi->n = chi->n * G->generators[k] % G->q;
+          chi->n = nmod_mul(chi->n, G->generators[k], G->mod);
           chi->expo[k] += G->PHI[k];
       }
       if (chi->expo[k] < G->expo)

acb_dirichlet/chi_vec_loop.c

@@ -29,24 +29,28 @@
 void
 acb_dirichlet_chi_vec_loop(ulong *v, ulong nv, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
 {
-  int j;
-  ulong t, k;
-  acb_dirichlet_conrey_t x;
-  acb_dirichlet_conrey_init(x, G);
-  acb_dirichlet_conrey_one(x, G);
-  t = v[1] = 0;
-  while ( (j = acb_dirichlet_conrey_next(x, G)) < G->num )
-  {
-    /* exponents were modified up to j */
-    for (k = 0; k < j; k++)
-      t = (t + chi->expo[k] * x->log[k]) % chi->order;
-    if (x->n < nv)
-      v[x->n] = t;
-  }
-  /* fix result outside primes */
-  acb_dirichlet_vec_set_null(v, nv, G);
-  /* copy outside modulus */
-  for (k = G->q + 1; k < nv ; k++ )
-    v[k] = v[k - G->q];
-  acb_dirichlet_conrey_clear(x);
+    int j;
+    ulong t, k;
+    acb_dirichlet_conrey_t x;
+    acb_dirichlet_conrey_init(x, G);
+    acb_dirichlet_conrey_one(x, G);
+
+    for (k = 0; k < nv; k++)
+        v[k] = CHI_NULL;
+
+    t = v[1] = 0;
+    while ( (j = acb_dirichlet_conrey_next(x, G)) < G->num )
+    {
+        /* exponents were modified up to j */
+        for (k = 0; k <= j; k++)
+            t = (t + chi->expo[k]) % chi->order;
+        if (x->n < nv)
+            v[x->n] = t;
+    }
+    /* fix result outside primes */
+    /*acb_dirichlet_vec_set_null(v, nv, G);*/
+    /* copy outside modulus */
+    for (k = G->q; k < nv ; k++ )
+        v[k] = v[k - G->q];
+    acb_dirichlet_conrey_clear(x);
 }

acb_dirichlet/conrey_exp.c

@@ -30,8 +30,7 @@ acb_dirichlet_conrey_exp(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G
 {
     ulong k, n = 1;
     for (k = G->neven; k < G->num; k++)
-        /* n = n_mulmod(n, n_powmod(G->generators[k], x->log[k], G->q), G->q); */
-        n = n * n_powmod(G->generators[k], x->log[k], G->q) % G->q;
+        n = nmod_mul(n, nmod_pow_ui(G->generators[k], x->log[k], G->mod), G->mod);
     x->n = n;
     return n;
 }

acb_dirichlet/conrey_first_primitive.c

@@ -26,15 +26,23 @@
 #include "acb_dirichlet.h"
 
 void
-acb_dirichlet_conrey_first_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) {
+acb_dirichlet_conrey_first_primitive(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
+{
     ulong k;
     if (G->q % 4 == 2)
     {
-      flint_printf("Exception (acb_dirichlet_conrey_first_primitive). no primitive element mod %wu.\n",G->q);
-      abort();
+        flint_printf("Exception (acb_dirichlet_conrey_first_primitive). No primitive element mod %wu.\n",G->q);
+        abort();
     }
+    x->n = 1;
     for (k = 0; k < G->num ; k++)
-      x->log[k] = 1;
-    if (G->neven == 2)
-      x->log[0] = 0;
+    {
+        if (k == 0 && G->neven == 2)
+            x->log[k] = 0;
+        else
+        {
+            x->n = nmod_mul(x->n, G->generators[k], G->mod);
+            x->log[k] = 1;
+        }
+    }
 }

acb_dirichlet/conrey_log.c

@@ -35,20 +35,31 @@ acb_dirichlet_conrey_log(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G
     ulong k, pk, gk;
     /* even part */
     if (G->neven >= 1)
-      x->log[0] = (m % 4 == 3);
-    if (G->neven == 2)
     {
-        ulong q_even = G->q_even;
-        ulong g2 = 5;
-        ulong m2 = (m % 4 == 3) ? -m % q_even : m % q_even;
-        x->log[1] = n_discrete_log_bsgs(m2, g2, q_even);
+      x->log[0] = (m % 4 == 3);
+      if (G->neven == 2)
+      {
+        ulong m2 = (x->log[0]) ? -m % G->q_even : m % G->q_even;
+        if (G->dlog == NULL)
+            x->log[1] = n_discrete_log_bsgs(m2, 5, G->q_even);
+        else
+            x->log[1] = dlog_precomp(G->dlog[1], m2);
+      }
     }
     /* odd part */
-    for (k = G->neven; k < G->num; k++)
+    if (G->dlog == NULL)
     {
-        pk = G->primepowers[k];
-        gk = G->generators[k] % pk;
-        x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk);
+        for (k = G->neven; k < G->num; k++)
+        {
+            pk = G->primepowers[k];
+            gk = G->generators[k] % pk;
+            x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk);
+        }
+    }
+    else
+    {
+        for (k = G->neven; k < G->num; k++)
+            x->log[k] = dlog_precomp(G->dlog[k], m % G->primepowers[k]);
     }
     /* keep value m */
     x->n = m;

acb_dirichlet/conrey_mul.c

@@ -30,6 +30,6 @@ acb_dirichlet_conrey_mul(acb_dirichlet_conrey_t c, const acb_dirichlet_group_t G
 {
     ulong k;
     for (k = 0; k < G->num ; k++)
-      c->log[k] = a->log[k] + b->log[k] % G->phi[k];
-    c->n = a->n * b->n % G->q;
+      c->log[k] = (a->log[k] + b->log[k]) % G->phi[k];
+    c->n = nmod_mul(a->n, b->n, G->mod);
 }

acb_dirichlet/conrey_next.c

@@ -28,16 +28,16 @@
 int
 acb_dirichlet_conrey_next(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
 {
-  /* update index */
-  ulong k;
-  for (k=0; k < G->num ; k++)
-  {
-    /* x->n = n_mulmod(x->n, G->generators[k], G->q); */
-    x->n = x->n * G->generators[k] % G->q;
-    if (++x->log[k] < G->phi[k])
-      break;
-    x->log[k] = 0;  
-  }
-  /* return last index modified */
-  return k;
+    /* update index */
+    int k;
+    for (k=0; k < G->num; k++)
+    {
+        x->n = nmod_mul(x->n, G->generators[k], G->mod);
+        x->log[k] += 1;
+        if (x->log[k] < G->phi[k])
+            break;
+        x->log[k] = 0;
+    }
+    /* return last index modified */
+    return k;
 }

acb_dirichlet/conrey_next_primitive.c

@@ -33,8 +33,7 @@ acb_dirichlet_conrey_next_primitive(acb_dirichlet_conrey_t x, const acb_dirichle
   ulong k = 0;
   if (G->neven == 2)
   {
-      /* x->n = n_mulmod(x->n, G->generators[0], G->q); */
-      x->n = x->n * G->generators[0] % G->q;
+      x->n = nmod_mul(x->n, G->generators[0], G->mod);
       if (++x->log[0] == 1)
           return 0;
       x->log[0] = 0;
@@ -42,12 +41,10 @@ acb_dirichlet_conrey_next_primitive(acb_dirichlet_conrey_t x, const acb_dirichle
   }
   for (; k < G->num ; k++)
   {
-    /* x->n = n_mulmod(x->n, G->generators[k], G->q); */
-    x->n = x->n * G->generators[k] % G->q;
+    x->n = nmod_mul(x->n, G->generators[k], G->mod);
     if (++x->log[k] % G->primes[k] == 0)
     {
-        /* x->n = n_mulmod(x->n, G->generators[k], G->q); */
-        x->n = x->n * G->generators[k] % G->q;
+        x->n = nmod_mul(x->n, G->generators[k], G->mod);
         ++x->log[k];
     }
     if (x->log[k] < G->phi[k])

acb_dirichlet/conrey_one.c

@@ -26,9 +26,10 @@
 #include "acb_dirichlet.h"
 
 void
-acb_dirichlet_conrey_one(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G) {
+acb_dirichlet_conrey_one(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G)
+{
     ulong k;
     for (k = 0; k < G->num ; k++)
-      x->log[k] = 0;
+        x->log[k] = 0;
     x->n = 1;
 }

acb_dirichlet/conrey_print.c

@@ -29,8 +29,9 @@ void
 acb_dirichlet_conrey_print(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
 {
     int k;
-    flint_printf("[%wu",x->log[0]);
+    if (G->num)
+        flint_printf("[%wu", x->log[0]);
     for (k = 1; k < G->num; k++)
-        flint_printf(", %wu",x->log[k]);
+        flint_printf(", %wu", x->log[k]);
     flint_printf("]\n");
 }

acb_dirichlet/group_dlog_precompute.c

@@ -0,0 +1,42 @@
+/*=============================================================================
+
+    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_group_dlog_precompute(acb_dirichlet_group_t G, ulong num)
+{
+    slong k;
+    G->dlog = flint_malloc(G->num * sizeof(dlog_precomp_t));
+    for (k = G->neven; k < G->num; k++)
+    {
+        ulong p, e, pe, a;
+        p = G->primes[k];
+        e = G->exponents[k];
+        pe = G->primepowers[k];
+        a = G->generators[k] % pe;
+        dlog_precomp_modpe_init(G->dlog[k], a, p, e, pe, num);
+    }
+}

acb_dirichlet/group_init.c

@@ -41,6 +41,8 @@ acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
     n_factor_t fac;
 
     G->q = q;
+    nmod_init(&G->mod, q);
+
     G->q_odd = q;
     G->q_even = 1;
 
@@ -63,6 +65,7 @@ acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
     G->generators = flint_malloc(G->num * sizeof(ulong));
     G->phi = flint_malloc(G->num * sizeof(ulong));
     G->PHI = flint_malloc(G->num * sizeof(ulong));
+    G->dlog = NULL;
 
     /* even part */
     G->expo = G->phi_q = 1;

acb_dirichlet/test/t-thetanull.c

@@ -0,0 +1,143 @@
+/*=============================================================================
+
+    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 <math.h>
+
+#define PI   3.14159265358
+#define LOG2 0.69314718055
+
+int main()
+{
+    slong prec = 80;
+    ulong q;
+
+    flint_printf("theta null....");
+    fflush(stdout);
+
+    /* check the only theta functions
+     *   theta(chi) = sum chi(k)* k^odd * exp(-Pi * k^2 / q)
+     * vanishing at 1 correspond to two specific
+     * characters of moduli 300 and 600 + conjugates
+     */
+
+    for (q = 3; q < 800; q ++)
+    {
+        acb_dirichlet_group_t G;
+        acb_dirichlet_conrey_t x;
+        acb_dirichlet_char_t chi;
+        ulong * v, nv, k;
+
+        acb_t zeta, sum;
+        acb_ptr z;
+
+        arb_t eq;
+        arb_ptr t, kt, tt;
+
+        if (q % 4 == 2)
+            continue;
+
+        acb_dirichlet_group_init(G, q);
+        acb_dirichlet_conrey_init(x, G);
+        acb_dirichlet_char_init(chi, G);
+
+        acb_init(zeta);
+        acb_dirichlet_zeta(zeta, G->expo, prec);
+
+        z = _acb_vec_init(G->expo);
+        _acb_vec_set_powers(z, zeta, G->expo, prec);
+
+        nv = ceil(sqrt((double)q * prec * LOG2 / PI)) + 2;
+        v = flint_malloc(nv * sizeof(ulong));
+
+        arb_init(eq);
+        arb_const_pi(eq, prec);
+        arb_div_ui(eq, eq, q, prec);
+        arb_neg(eq, eq);
+        arb_exp(eq, eq, prec);
+
+        t = _arb_vec_init(nv);
+        acb_dirichlet_arb_quadratic_powers(t, nv, eq, prec);
+
+        kt = _arb_vec_init(nv);
+        for (k = 1; k < nv; k++)
+            arb_mul_ui(kt + k, t + k, k, prec);
+
+        /* theta function on primitive characters */
+        acb_init(sum);
+        acb_dirichlet_conrey_first_primitive(x, G);
+
+        while (1) {
+            ulong m;
+            acb_zero(sum);
+
+            acb_dirichlet_char_conrey(chi, G, x);
+            acb_dirichlet_chi_vec_loop(v, nv, G, chi);
+
+            m = G->expo / chi->order;
+            /*
+            flint_printf("Theta(chi_%wu(%wu)) (m=%wu)\n", q, chi->n, m);
+            */
+            tt = acb_dirichlet_char_is_odd(G, chi) ? kt : t;
+            for (k = 1; k < nv; k++)
+                if (v[k] != CHI_NULL)
+                    acb_addmul_arb(sum, z + (v[k] * m), tt + k, prec);
+
+            if ((q == 300 && (chi->n == 271 || chi->n == 131))
+                    || (q == 600 && (chi->n == 11 || chi->n == 91)))
+            {
+                if (!acb_contains_zero(sum))
+                {
+                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n);
+                acb_printd(sum, 10);
+                flint_printf("\n");
+                abort();
+                }
+            }
+            else if (acb_contains_zero(sum))
+            {
+                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n);
+                acb_printd(sum, 10);
+                flint_printf("\n");
+                abort();
+            }
+
+            if (acb_dirichlet_conrey_next_primitive(x, G) == G->num)
+                break;
+        }
+        _acb_vec_clear(z, G->expo);
+        _arb_vec_clear(t, nv);
+        acb_clear(zeta);
+        acb_clear(sum);
+        arb_clear(eq);
+        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/zeta.c

@@ -0,0 +1,36 @@
+/*=============================================================================
+
+    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_zeta(acb_t res, ulong order, slong prec)
+{
+    fmpq_t t;
+    fmpq_init(t);
+    fmpq_set_si(t, 2, order);
+    arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+    fmpq_clear(t);
+}

dlog.h

@@ -191,6 +191,7 @@ ulong dlog_bsgs(const dlog_bsgs_t t, ulong b);
 ulong dlog_rho(const dlog_rho_t t, ulong b);
 /*#define dlog_bsgs(t, b) n_discrete_log_bsgs_table(t, b)*/
 
+void dlog_precomp_modpe_init(dlog_precomp_t pre, ulong a, ulong p, ulong e, ulong pe, ulong num);
 void dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num);
 void dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num);
 void dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, ulong pe, ulong num);