[dirichlet] add conrey type to handle logs + char type

- try to handle even and odd components the same way in the dirichlet group

- switch from phi_q_odd to smaller expo = exponent of the group

  all character orders divide this number, and a character of that order exists

- use conrey logarithm to reuse log and to loop efficiently over the group

  (see the diff on l.c, only 1 log in computed instead of 2 * q)

- NOT TESTED, for the moment it just compiles, I know some errors
  (e.g. the FIXME in group_init.c : the generators have to be lifted mod q)
  this commit is just a proof of concept.

acb_dirichlet.h

@@ -24,18 +24,30 @@ typedef struct
     ulong q;                /* modulus */
     ulong q_even;           /* even part of modulus */
     ulong q_odd;            /* odd part of modulus */
-    ulong phi_q;            /* phi(q) */
-    ulong phi_q_odd;        /* phi(q_odd) */
-    slong num;              /* number of odd prime factors */
-    ulong * primes;         /* odd primes p[k] */
+    ulong phi_q;            /* phi(q) = group size */
+    ulong expo;             /* group exponent = lcm(phi(q_even), phi(p[k]^e[k]) ) */
+    slong neven;            /* number of even components */
+    slong num;              /* number of prime components (even + odd) */
+    ulong * primes;         /* primes p[k] */
     ulong * exponents;      /* exponents e[k] */
-    ulong * generators;     /* generator for each odd prime p[k] */
-    ulong * PHI;            /* PHI(k) = phi(q_odd) / phi(p[k]^e[k])  */
+    ulong * generators;     /* generator for each prime p[k] */
+    ulong * phi;            /* phi(k) = phi(p[k]^e[k])       */
+    ulong * PHI;            /* PHI(k) = expo / phi(k)        */
 }
 acb_dirichlet_group_struct;
 
 typedef acb_dirichlet_group_struct acb_dirichlet_group_t[1];
 
+/* elements of the group, keep both number and log */
+typedef struct
+{
+  ulong n;           /* number */
+  ulong * log;       /* s.t. prod generators[k]^log[k] = number */
+}
+acb_conrey_struct;
+
+typedef acb_conrey_struct acb_conrey_t[1];
+
 void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s,
     const signed char * chi, int mod, int reciprocal, slong prec);
 
@@ -45,8 +57,52 @@ void acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q);
 
 void acb_dirichlet_group_clear(acb_dirichlet_group_t G);
 
+void acb_conrey_init(acb_conrey_t x, const acb_dirichlet_group_t G);
+
+void acb_conrey_clear(acb_conrey_t x);
+
+void acb_conrey_one(acb_conrey_t x, const acb_dirichlet_group_t G);
+
+void acb_conrey_log(acb_conrey_t x, const acb_dirichlet_group_t G, ulong m);
+
+int  acb_conrey_next(acb_conrey_t x, const acb_dirichlet_group_t G);
+
+long n_dirichlet_chi_conrey(const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b);
+
+long n_dirichlet_chi(const acb_dirichlet_group_t G, ulong m, ulong n);
+
+void acb_dirichlet_chi_conrey(acb_t res, const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b, slong prec);
+
 void acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec);
 
+/* introducing character type */
+
+/* character = reduced exponents, keep its order and number */
+typedef struct
+{
+  ulong q;      /* modulus */
+  ulong n;      /* number */
+  ulong order;  /* order */
+  ulong * expo; /* reduced exponents ( order * log[k] / gcd( ) ) */
+}
+acb_dirichlet_char_struct;
+
+typedef acb_dirichlet_char_struct acb_dirichlet_char_t[1];
+
+void acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
+
+void acb_dirichlet_char_clear(acb_dirichlet_char_t chi);
+
+void acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n);
+
+long n_dirichlet_char_eval(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n);
+
+void acb_dirichlet_char_eval(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec);
+
+void n_dirichlet_char_vec(long *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong nv);
+
+void acb_dirichlet_char_vec(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec);
+
 #ifdef __cplusplus
 }
 #endif

acb_dirichlet/char.c

@@ -0,0 +1,66 @@
+/*=============================================================================
+
+    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) 2015 Jonathan Bober
+    Copyright (C) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G) {
+  chi->expo = flint_malloc(G->num * sizeof(ulong));
+}
+
+void
+acb_dirichlet_char_clear(acb_dirichlet_char_t chi) {
+    flint_free(chi->expo);
+}
+
+void
+acb_dirichlet_char_log(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, const acb_conrey_t x)
+{
+  ulong k, g;
+  g = G->expo;
+  /* modify exponents */
+  for (k = 0; k < G->num; k++)
+  {
+    chi->expo[k] = (x->log[k] * G->PHI[k]) % G->expo;
+    g = n_gcd(g, chi->expo[k]);
+  }
+  for (k = 0; k < G->num; k++)
+    chi->expo[k] = chi->expo[k] / g;
+
+  chi->q = G->q;
+  chi->order = G->expo / g;
+  chi->n = x->n;
+}
+
+void
+acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n)
+{
+  acb_conrey_t x;
+  x->log = chi->expo;
+  acb_conrey_log(x, G, n);
+  acb_dirichlet_char_log(chi, G, x);
+}
+

acb_dirichlet/chareval.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) 2015 Jonathan Bober
+    Copyright (C) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+long
+n_dirichlet_char_eval(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n)
+{
+  ulong v = 0, k;
+  acb_conrey_t x;
+  acb_conrey_init(x, G);
+  acb_conrey_log(x, G, n);
+  for (k = 0; k < G->num; k++)
+    v = (v + chi->expo[k] * x->log[k]) % chi->order;
+  acb_conrey_clear(x);
+  return v;
+}
+
+void
+fmpq_dirichlet_char_eval(fmpq_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n)
+{
+    fmpq_set_si(res, n_dirichlet_char_eval(G, chi, n), chi->order);
+}
+
+void
+acb_dirichlet_char_eval(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong n, slong prec)
+{
+    fmpq_t t;
+    fmpq_init(t);
+    fmpq_dirichlet_char_eval(t, G, chi, n);
+    arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+    fmpq_clear(t);
+}

acb_dirichlet/charevalvec.c

@@ -0,0 +1,67 @@
+/*=============================================================================
+
+    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) 2015 Jonathan Bober
+    Copyright (C) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+static void
+set_non_invertible_values(long *v, const acb_dirichlet_group_t G, ulong nv)
+{
+  ulong k, l;
+  if (G->q_even > 1)
+  {
+    for (k = 2; k < nv; k += 2)
+      v[k] = -1;
+  }
+  for (l = 0; l < G->num; l++)
+  {
+    ulong p = G->primes[k];
+    for (k = p; k < nv; k += p)
+      v[k] = -1;
+  }
+}
+
+void
+n_dirichlet_char_vec(long *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, ulong nv)
+{
+  ulong t, k, j;
+  acb_conrey_t x;
+  acb_conrey_init(x, G);
+  acb_conrey_one(x, G);
+  t = v[1] = 0;
+  while( (j = acb_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 */
+  set_non_invertible_values(v, G, nv);
+  /* copy outside modulus */
+  for (k = G->q + 1; k < nv ; k++ )
+    v[k] = v[k-G->q];
+  acb_conrey_clear(x);
+}

acb_dirichlet/chi.c

@@ -14,114 +14,50 @@
 
 /* todo: modular arithmetic */
 
-static ulong
-chi_odd_exponent(const acb_dirichlet_group_t G, ulong m, ulong n)
+long
+n_dirichlet_chi_conrey(const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b)
 {
-    ulong x, k, pk, gk, logm, logn;
-
-    x = 0;
-
-    for (k = 0; k < G->num; k++)
-    {
-        pk = n_pow(G->primes[k], G->exponents[k]);
-        gk = G->generators[k] % pk;
-
-        logm = n_discrete_log_bsgs(m % pk, gk, pk);
-        logn = n_discrete_log_bsgs(n % pk, gk, pk);
-
-        x = (x + G->PHI[k] * logm * logn) % G->phi_q_odd;
-    }
-
-    return x;
+  ulong x, k;
+  x = 0;
+  for (k = 0; k < G->num; k++)
+    x = (x + G->PHI[k] * a->log[k] * b->log[k]) % G->expo;
+  return x;
 }
 
-static ulong
-chi_even_exponent(const acb_dirichlet_group_t G, ulong m, ulong n)
+long
+n_dirichlet_chi(const acb_dirichlet_group_t G, ulong m, ulong n)
 {
-    ulong x;
-    ulong q_even = G->q_even;
+  ulong x;
+  acb_conrey_t a, b;
+  acb_conrey_init(a, G);
+  acb_conrey_init(b, G);
 
-    if (q_even <= 2)
-        return 0;
+  acb_conrey_log(a, G, m);
+  acb_conrey_log(b, G, n);
+  x = n_dirichlet_chi_conrey(G, a, b);
 
-    x = 0;
+  acb_conrey_clear(a);
+  acb_conrey_clear(b);
 
-    if ((m % 4 == 3) && (n % 4 == 3))
-        x = q_even / 8;
+  return x;
+}
 
-    if (q_even > 4)
-    {
-        ulong g2, logm, logn;
-
-        g2 = 5;
-
-        if (m % 4 == 3)
-        {
-            m = n_negmod(m, q_even);
-        }
-
-        if (n % 4 == 3)
-        {
-            n = n_negmod(n, q_even);
-        }
-
-        logm = n_discrete_log_bsgs(m % q_even, g2, q_even);
-        logn = n_discrete_log_bsgs(n % q_even, g2, q_even);
-        x += logm * logn;
-    }
-
-    return x % (q_even / 4);
+void
+acb_dirichlet_chi_conrey(acb_t res, const acb_dirichlet_group_t G, const acb_conrey_t a, const acb_conrey_t b, slong prec)
+{
+  fmpq_t t;
+  fmpq_init(t);
+  fmpq_set_si(t, n_dirichlet_chi_conrey(G, a, b), G->expo);
+  arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+  fmpq_clear(t);
 }
 
 void
 acb_dirichlet_chi(acb_t res, const acb_dirichlet_group_t G, ulong m, ulong n, slong prec)
 {
-    fmpq_t t, u;
-
-    ulong odd_part, even_part;
-    ulong q_even = G->q_even;
-    ulong q_odd = G->q_odd;
-
-    if ((q_even > 1 && (n % 2 == 0)) || (q_odd > 1 && (n_gcd(q_odd, n) != 1)))
-    {
-        acb_zero(res);
-        return;
-    }
-
-    if (q_even > 2)
-    {
-        if (q_even == 4)
-        {
-            if (m % 4 == 3 && n % 4 == 3)
-                even_part = q_even / 2;  /* -1 */
-            else
-                even_part = 0;           /* 1 */
-        }
-        else
-        {
-            even_part = 4 * chi_even_exponent(G, m % q_even, n % q_even);
-        }
-    }
-    else
-    {
-        even_part = 0;
-    }
-
-    if (q_odd > 1)
-        odd_part = chi_odd_exponent(G, m % q_odd, n % q_odd);
-    else
-        odd_part = 0;
-
-    fmpq_init(t);
-    fmpq_init(u);
-
-    fmpq_set_si(t, 2 * even_part, q_even);
-    fmpq_set_si(u, 2 * odd_part, G->phi_q_odd);
-    fmpq_add(t, t, u);
-
-    arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
-
-    fmpq_clear(t);
-    fmpq_clear(u);
+  fmpq_t t;
+  fmpq_init(t);
+  fmpq_set_si(t, n_dirichlet_chi(G, m, n), G->expo);
+  arb_sin_cos_pi_fmpq(acb_imagref(res), acb_realref(res), t, prec);
+  fmpq_clear(t);
 }
-

acb_dirichlet/conrey_log.c

@@ -0,0 +1,91 @@
+/*=============================================================================
+
+    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) 2015 Jonathan Bober
+    Copyright (C) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_dirichlet.h"
+
+/* todo: modular arithmetic
+   discrete log can be computed along exponents or using p-adic log
+*/
+
+void
+acb_conrey_init(acb_conrey_t x, const acb_dirichlet_group_t G) {
+    x->log = flint_malloc(G->num * sizeof(ulong));
+}
+
+void
+acb_conrey_one(acb_conrey_t x, const acb_dirichlet_group_t G) {
+    ulong k;
+    for (k = 0; k < G->num ; k++)
+      x->log[k] = 0;
+    x->n = 1;
+}
+
+void
+acb_conrey_clear(acb_conrey_t x) {
+    flint_free(x->log);
+}
+
+/* TODO: use precomputations in G if present */
+void
+acb_conrey_log(acb_conrey_t x, const acb_dirichlet_group_t G, ulong m)
+{
+    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) ? n_negmod(m, q_even) : m % q_even;
+        x->log[1] = n_discrete_log_bsgs(m2, g2, q_even);
+    }
+    /* odd part */
+    for (k = G->neven; k < G->num; k++)
+    {
+        pk = n_pow(G->primes[k], G->exponents[k]);
+        gk = G->generators[k] % pk;
+        x->log[k] = n_discrete_log_bsgs(m % pk, gk, pk);
+    }
+    /* keep value m */
+    x->n = m;
+}
+
+int
+acb_conrey_next(acb_conrey_t x, const acb_dirichlet_group_t G)
+{
+  /* update index */
+  ulong k;
+  for (k=0; k < G->num ; k++)
+  {
+    x->n = x->n * G->generators[k];
+    if(x->log[k]++ < G->phi[k])
+      break;
+    x->log[k] = 0;  
+  }
+  /* return last index modified */
+  return k;
+}

acb_dirichlet/group_init.c

@@ -36,6 +36,7 @@ void
 acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
 {
     slong k;
+    ulong e2 = 0;
     n_factor_t fac;
 
     G->q = q;
@@ -46,39 +47,56 @@ acb_dirichlet_group_init(acb_dirichlet_group_t G, ulong q)
     {
         G->q_odd /= 2;
         G->q_even *= 2;
+        e2++;
     }
 
     n_factor_init(&fac);
     n_factor(&fac, G->q_odd, 1);
 
-    G->num = fac.num;
+    /* number of components at p=2 */
+    G->neven = (e2 >= 3) ? 2 : (e2 == 2) ? 1 : 0;
+    G->num = G->neven + fac.num;
     G->primes = flint_malloc(G->num * sizeof(ulong));
     G->exponents = flint_malloc(G->num * sizeof(ulong));
     G->generators = flint_malloc(G->num * sizeof(ulong));
+    G->phi = flint_malloc(G->num * sizeof(ulong));
     G->PHI = flint_malloc(G->num * sizeof(ulong));
 
-    for (k = 0; k < G->num; k++)
+    /* even part */
+    if (G->q_even <= 2)
     {
-        G->primes[k] = fac.p[k];
-        G->exponents[k] = fac.exp[k];
+      G->expo = G->phi_q = 1;
+    }
+    else
+    {
+      G->phi_q = G->q_even / 2;
+      G->expo = G->phi_q / 2;
     }
 
-    G->phi_q_odd = 1;
-    for (k = 0; k < G->num; k++)
-        G->phi_q_odd *= (G->primes[k] - 1) * n_pow(G->primes[k], G->exponents[k]-1);
+    for (k = 0; k < G->neven; k++)
+    {
+        G->primes[k] = 2;
+        G->exponents[k] = (k==0) ? 1 : e2-2;
+        G->generators[k] = (k==0) ? -1 : 5;
+        G->phi[k] = (k==0) ? 1 : G->expo;
+    }
 
-    if (G->q_even == 1)
-        G->phi_q = G->phi_q_odd;
-    else
-        G->phi_q = G->phi_q_odd * (G->q_even / 2);
+    for (k = G->neven; k < G->num; k++)
+    {
+        ulong phik, p1;
+        G->primes[k] = fac.p[k - G->neven];
+        G->exponents[k] = fac.exp[k - G->neven];
+        G->generators[k] = primitive_root_p_and_p2(G->primes[k]);
+        p1 = G->primes[k] - 1;
+        phik = p1 * n_pow(G->primes[k], G->exponents[k]-1);
+        G->expo *= phik / n_gcd(G->expo, p1);
+        G->phi_q *= phik;
+        G->phi[k] = phik;
+    }
 
     for (k = 0; k < G->num; k++)
     {
-        ulong phi;
-
-        G->generators[k] = primitive_root_p_and_p2(G->primes[k]);
-        phi = n_pow(G->primes[k], G->exponents[k] - 1) * (G->primes[k] - 1);
-        G->PHI[k] = G->phi_q_odd / phi;
+        G->PHI[k] = G->expo / G->phi[k];
+        /* FIXME: generators[k] should be lifted mod q! */
     }
 }
-

acb_dirichlet/l.c

@@ -16,26 +16,30 @@ acb_dirichlet_l(acb_t res, const acb_t s,
     const acb_dirichlet_group_t G, ulong m, slong prec)
 {
     acb_t chi, t, u, a;
-    ulong k;
+    acb_conrey_t cm, cn;
 
     acb_init(chi);
+    acb_conrey_init(cm, G);
+    acb_conrey_init(cn, G);
     acb_init(t);
     acb_init(u);
     acb_init(a);
 
+    acb_conrey_log(cm, G, m);
+    acb_conrey_one(cn, G);
     acb_zero(t);
 
-    for (k = 1; k <= G->q; k++)
-    {
-        acb_dirichlet_chi(chi, G, m, k, prec);
+    while (1) {
+        /* todo: use n_dirichlet_chi and precomputed roots instead */
+        acb_dirichlet_chi_conrey(chi, G, cm, cn, prec);
 
-        if (!acb_is_zero(chi))
-        {
-            acb_set_ui(a, k);
-            acb_div_ui(a, a, G->q, prec);
-            acb_hurwitz_zeta(u, s, a, prec);
-            acb_addmul(t, chi, u, prec);
-        }
+        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);
+
+        if (acb_conrey_next(cn, G) == G->num)
+          break;
     }
 
     acb_set_ui(u, G->q);
@@ -43,6 +47,9 @@ acb_dirichlet_l(acb_t res, const acb_t s,
     acb_pow(u, u, a, prec);
     acb_mul(res, t, u, prec);
 
+    acb_conrey_clear(cm);
+    acb_conrey_clear(cn);
+
     acb_clear(chi);
     acb_clear(t);
     acb_clear(u);