add conrey index inside char

acb_dirichlet.h

@@ -85,6 +85,15 @@ void acb_dirichlet_conrey_init(acb_dirichlet_conrey_t x, const acb_dirichlet_gro
 void acb_dirichlet_conrey_clear(acb_dirichlet_conrey_t x);
 void acb_dirichlet_conrey_print(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
 
+ACB_DIRICHLET_INLINE void
+acb_dirichlet_conrey_copy(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t y)
+{
+    slong k;
+    x->n = y->n;
+    for (k = 0; k < G->num; k++)
+        x->log[k] = y->log[k];
+}
+
 int acb_dirichlet_conrey_parity(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
 ulong acb_dirichlet_conrey_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
 void acb_dirichlet_conrey_log(acb_dirichlet_conrey_t x, const acb_dirichlet_group_t G, ulong m);
@@ -112,10 +121,9 @@ void acb_dirichlet_pairing(acb_t res, const acb_dirichlet_group_t G, ulong m, ul
 typedef struct
 {
     ulong q;           /* modulus */
-    /*? acb_dirichet_conrey_struct n; */
-    ulong n;           /* number */
     ulong order;       /* order */
-    ulong * expo;      /* reduced exponents ( order * log[k] / gcd( ) ) */
+    acb_dirichlet_conrey_t x;
+    ulong * expo;      /* reduced exponents ( log[k] * PHI[k] / gcd( ) ) */
     int parity;        /* 0 for even char, 1 for odd */
     ulong conductor;
 }
@@ -129,6 +137,12 @@ acb_dirichlet_char_order(const acb_dirichlet_char_t chi)
     return chi->order;
 }
 
+ACB_DIRICHLET_INLINE ulong
+acb_dirichlet_char_conductor(const acb_dirichlet_char_t chi)
+{
+    return chi->conductor;
+}
+
 ACB_DIRICHLET_INLINE int
 acb_dirichlet_char_parity(const acb_dirichlet_char_t chi)
 {
@@ -141,14 +155,20 @@ void acb_dirichlet_char_print(const acb_dirichlet_group_t G, const acb_dirichlet
 
 void acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n);
 void acb_dirichlet_char_conrey(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x);
+void acb_dirichlet_char_set_expo(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
 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);
+void acb_dirichlet_char_mul(acb_dirichlet_char_t chi12, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2);
+
+/*ulong acb_dirichlet_char_conductor(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);
 
+int acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G);
+int acb_dirichlet_char_next_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);
 

acb_dirichlet/char.c

@@ -28,8 +28,6 @@
 void
 acb_dirichlet_char(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, ulong n)
 {
-    acb_dirichlet_conrey_t x;
-    x->log = chi->expo;
-    acb_dirichlet_conrey_log(x, G, n);
-    acb_dirichlet_char_conrey(chi, G, x);
+    acb_dirichlet_conrey_log(chi->x, G, n);
+    acb_dirichlet_char_conrey(chi, G, chi->x);
 }

acb_dirichlet/char_clear.c

@@ -28,5 +28,6 @@
 void
 acb_dirichlet_char_clear(acb_dirichlet_char_t chi)
 {
+    acb_dirichlet_conrey_clear(chi->x);
     flint_free(chi->expo);
 }

acb_dirichlet/char_conrey.c

@@ -30,15 +30,17 @@
 void
 acb_dirichlet_char_conrey(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G, const acb_dirichlet_conrey_t x)
 {
-    ulong k;
+    /* assume chi->x already set if x == NULL */
+    if (x == NULL)
+        x = chi->x;
+    else
+        acb_dirichlet_conrey_copy(chi->x, G, x);
+
     chi->q = G->q;
-    chi->n = x->n;
     chi->parity = acb_dirichlet_conrey_parity(G, x);
     chi->conductor = acb_dirichlet_conrey_conductor(G, x);
 
-    for (k = 0; k < G->num; k++)
-        chi->expo[k] = (x->log[k] * G->PHI[k]) % G->expo;
-
+    acb_dirichlet_char_set_expo(chi, G);
     /* optional: divide by gcd to obtain true order */
     acb_dirichlet_char_normalize(chi, G);
 }

acb_dirichlet/char_first_primitive.c

@@ -28,11 +28,7 @@
 void
 acb_dirichlet_char_first_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
-    acb_dirichlet_conrey_t x;
-    chi->q = G->q;
-    x->log = chi->expo;
-    acb_dirichlet_conrey_first_primitive(x, G);
-    chi->n = x->n;
-    chi->conductor = chi->q;
+    acb_dirichlet_conrey_first_primitive(chi->x, G);
+    acb_dirichlet_char_conrey(chi, G, NULL);
     acb_dirichlet_char_normalize(chi, G);
 }

acb_dirichlet/char_init.c

@@ -28,5 +28,6 @@
 void
 acb_dirichlet_char_init(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
+    acb_dirichlet_conrey_init(chi->x, G);
     chi->expo = flint_malloc(G->num * sizeof(ulong));
 }

acb_dirichlet/char_mul.c

@@ -25,8 +25,9 @@
 
 #include "acb_dirichlet.h"
 
-ulong
-acb_dirichlet_char_conductor(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
+void
+acb_dirichlet_char_mul(acb_dirichlet_char_t chi12, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2)
 {
-    return chi->conductor;
+    acb_dirichlet_conrey_mul(chi12->x, G, chi1->x, chi2->x);
+    acb_dirichlet_char_conrey(chi12, G, chi12->x);
 }

acb_dirichlet/char_next.c

@@ -28,24 +28,9 @@
 int
 acb_dirichlet_char_next(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
-    ulong k;
-
-    acb_dirichlet_char_denormalize(chi, G);
-
-    /* update index */
-    for (k=0; k < G->num ; k++)
-    {
-        chi->n = nmod_mul(chi->n, G->generators[k], G->mod);
-        chi->expo[k] += G->PHI[k];
-        if (chi->expo[k] < G->expo)
-            break;
-        chi->expo[k] = 0;  
-    }
-
-    /* should do it with log, but log is not kept yet in the struct */
-    chi->conductor = acb_dirichlet_ui_conductor(G, chi->n);
-    acb_dirichlet_char_normalize(chi, G);
-
+    int k;
+    k = acb_dirichlet_conrey_next(chi->x, G);
+    acb_dirichlet_char_conrey(chi, G, NULL);
     /* return last index modified */
     return k;
 }

acb_dirichlet/char_next_primitive.c

@@ -28,38 +28,9 @@
 int
 acb_dirichlet_char_next_primitive(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
-    ulong k;
-
-    acb_dirichlet_char_denormalize(chi, G);
-
-    /* update index */
-    k = 0;
-    if (G->neven == 2)
-    {
-        chi->n = nmod_mul(chi->n, G->generators[0], G->mod);
-        chi->expo[0]++;
-        if (chi->expo[0] < G->expo)
-            return 0;
-        chi->expo[0] = 0;
-        k = 1;
-    }
-
-    for (; k < G->num ; k++)
-    {
-        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 = nmod_mul(chi->n, G->generators[k], G->mod);
-            chi->expo[k] += G->PHI[k];
-        }
-        if (chi->expo[k] < G->expo)
-            break;
-        chi->expo[k] = G->PHI[k];  
-    }
-
-    acb_dirichlet_char_normalize(chi, G);
-
+    int k;
+    k = acb_dirichlet_conrey_next(chi->x, G);
+    acb_dirichlet_char_conrey(chi, G, NULL);
     /* return last index modified */
     return k;
 }

acb_dirichlet/char_normalize.c

@@ -25,6 +25,14 @@
 
 #include "acb_dirichlet.h"
 
+void
+acb_dirichlet_char_set_expo(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
+{
+    slong k;
+    for (k = 0; k < G->num; k++)
+        chi->expo[k] = (chi->x->log[k] * G->PHI[k]) % G->expo;
+}
+
 void
 acb_dirichlet_char_normalize(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {

acb_dirichlet/char_one.c

@@ -28,13 +28,10 @@
 void
 acb_dirichlet_char_one(acb_dirichlet_char_t chi, const acb_dirichlet_group_t G)
 {
-    ulong k;
+    acb_dirichlet_conrey_one(chi->x, G);
     chi->q = G->q;
-    chi->n = 1;
-
-    for (k = 0; k < G->num; k++)
-        chi->expo[k] = 0;
-
-    chi->order = 1;
     chi->conductor = 1;
+    chi->parity = 0;
+    acb_dirichlet_char_set_expo(chi, G);
+    acb_dirichlet_char_normalize(chi, G);
 }

acb_dirichlet/char_print.c

@@ -29,7 +29,9 @@ void
 acb_dirichlet_char_print(const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi)
 {
     acb_dirichlet_conrey_t x;
+    flint_printf("chi_%wu(%wu,.) of order %wu, index ", G->q, chi->x->n, chi->order);
+    acb_dirichlet_conrey_print(G, chi->x);
+    flint_printf(" and exponents ");
     x->log = chi->expo;
-    flint_printf("chi_%wu(%wu,.) of order %wu and index ", G->q, chi->n, chi->order);
     acb_dirichlet_conrey_print(G, x);
 }

acb_dirichlet/chi_theta_arb.c

@@ -39,7 +39,7 @@ acb_dirichlet_arb_theta_argt(arb_t x, ulong q, const arb_t t, 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 t, slong prec)
+acb_dirichlet_chi_theta_arb(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;

acb_dirichlet/conrey_print.c

@@ -33,5 +33,5 @@ acb_dirichlet_conrey_print(const acb_dirichlet_group_t G, const acb_dirichlet_co
         flint_printf("[%wu", x->log[0]);
     for (k = 1; k < G->num; k++)
         flint_printf(", %wu", x->log[k]);
-    flint_printf("]\n");
+    flint_printf("]");
 }

acb_dirichlet/gauss_sum.c

@@ -28,7 +28,9 @@
 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->x->n == 1)
+        acb_set_si(res, -1);
+    else if (chi->order <= 2)
     {
         if (chi->parity)
         {

acb_dirichlet/gauss_sum_theta.c

@@ -33,8 +33,8 @@ acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const ac
 
     arb_init(x);
 
-    if ((G->q == 300 && (chi->n == 271 || chi->n == 131))
-            || (G->q == 600 && (chi->n == 11 || chi->n == 91)))
+    if ((G->q == 300 && (chi->x->n == 271 || chi->x->n == 131))
+            || (G->q == 600 && (chi->x->n == 11 || chi->x->n == 91)))
     {
         /* or could use l'Hopital rule */
         acb_dirichlet_gauss_sum_naive(res, G, chi, prec);
@@ -42,7 +42,7 @@ acb_dirichlet_gauss_sum_theta(acb_t res, const acb_dirichlet_group_t G, const ac
     }
 
     arb_one(x);
-    acb_dirichlet_chi_theta(res, G, chi, x, prec);
+    acb_dirichlet_chi_theta_arb(res, G, chi, x, prec);
     acb_init(eps);
     acb_conj(eps, res);
     acb_div(eps, res, eps, prec);

acb_dirichlet/jacobi_sum.c

@@ -0,0 +1,64 @@
+/*=============================================================================
+
+    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_jacobi_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1, const acb_dirichlet_char_t chi2, slong prec)
+{
+    if (chi1->x->n == 1 || chi2->x->n == 1)
+    {
+        if (chi1->x->n == 1 && chi2->x->n == 1)
+            acb_set_si(res, G->q - 2);
+        else
+            acb_set_si(res, - 1);
+    }
+    else if (nmod_mul(chi1->x->n, chi2->x->n, G->mod) == 1)
+    {
+        if (chi1->parity)
+            acb_one(res);
+        else
+            acb_set_si(res, -1);
+    }
+    else
+    {
+        acb_dirichlet_char_t chi12;
+        acb_t tmp;
+
+        acb_dirichlet_char_init(chi12, G);
+        acb_init(tmp);
+
+        acb_dirichlet_gauss_sum(res, G, chi1, prec);
+        acb_dirichlet_gauss_sum(tmp, G, chi2, prec);
+        acb_mul(res, res, tmp, prec);
+        acb_dirichlet_char_mul(chi12, G, chi1, chi2);
+        acb_dirichlet_gauss_sum(tmp, G, chi12, prec);
+        acb_div(res, res, tmp, prec);
+
+        acb_dirichlet_char_clear(chi12);
+        acb_clear(tmp);
+    }
+
+}

acb_dirichlet/profile/p-conrey.c

@@ -46,7 +46,7 @@ int main()
     TIMEIT_ONCE_STOP
     nref = n;
 
-    flint_printf("conrey elements.... ");
+    flint_printf("conrey.... ");
     TIMEIT_ONCE_START
     for (n = 0, q = 2; q <= maxq; q++)
     {
@@ -70,6 +70,31 @@ int main()
             abort();
     }
 
+    flint_printf("chars.... ");
+    TIMEIT_ONCE_START
+    for (n = 0, q = 2; q <= maxq; q++)
+    {
+        acb_dirichlet_group_t G;
+        acb_dirichlet_char_t chi;
+
+        acb_dirichlet_group_init(G, q);
+        acb_dirichlet_char_init(chi, G);
+
+        acb_dirichlet_char_one(chi, G);
+        n++;
+
+        for (; acb_dirichlet_char_next(chi, G) < G->num; n++);
+        acb_dirichlet_char_clear(chi);
+        acb_dirichlet_group_clear(G);
+    }
+    TIMEIT_ONCE_STOP
+    if (n != nref)
+    {
+            flint_printf("FAIL: wrong number of elements %wu != %wu\n\n",n, nref);
+            abort();
+    }
+
+
     flint_cleanup();
     flint_printf("PASS\n");
     return EXIT_SUCCESS;

acb_dirichlet/profile/p-vec.c

@@ -51,7 +51,7 @@ vecloop(dir_f dir, ulong minq, ulong maxq, ulong * rand, ulong nr, ulong * v, ul
         for (r = 0; r < nr; r++)
         {
             acb_dirichlet_char(chi, G, rand[r] % q);
-            dir(v, nv, G, chi);
+            dir(v, G, chi, nv);
         }
 
         acb_dirichlet_char_clear(chi);
@@ -97,15 +97,15 @@ int main()
 
         flint_printf("big loop................ ");
         fflush(stdout);
-        vecloop(acb_dirichlet_chi_vec_loop, minq, maxq, rand, nr, v, nv);
+        vecloop(acb_dirichlet_ui_chi_vec_loop, minq, maxq, rand, nr, v, nv);
 
         flint_printf("med loop................ ");
         fflush(stdout);
-        vecloop(acb_dirichlet_chi_vec_primeloop, minq, maxq, rand, nr, v, nv);
+        vecloop(acb_dirichlet_ui_chi_vec_primeloop, minq, maxq, rand, nr, v, nv);
 
         flint_printf("sieve................... ");
         fflush(stdout);
-        vecloop(acb_dirichlet_chi_vec_sieve, minq, maxq, rand, nr, v, nv);
+        vecloop(acb_dirichlet_ui_chi_vec_sieve, minq, maxq, rand, nr, v, nv);
 
         /*
         flint_printf("generic........ ");

acb_dirichlet/test/t-chars.c

@@ -76,13 +76,13 @@ int main()
             }
 
             cond = acb_dirichlet_ui_conductor(G, m);
-            if (cond != acb_dirichlet_char_conductor(G, chi))
+            if (cond != chi->conductor)
             {
                 flint_printf("FAIL: conductor\n\n");
                 flint_printf("q = %wu\n\n", q);
                 flint_printf("m = %wu\n\n", m);
                 flint_printf("conductor(m) = %wu\n\n", cond);
-                flint_printf("chi->conductor = %wu\n\n", acb_dirichlet_char_conductor(G, chi));
+                flint_printf("chi->conductor = %wu\n\n", chi->conductor);
                 abort();
             }
 

acb_dirichlet/test/t-gauss.c

@@ -67,7 +67,7 @@ int main()
             if (!acb_overlaps(s1, s2)
                     || !acb_overlaps(s1, s3))
             {
-                flint_printf("FAIL: G(chi_%wu(%wu))\n\n", q, chi->n);
+                flint_printf("FAIL: G(chi_%wu(%wu))\n\n", q, chi->x->n);
                 flint_printf("\nnaive ", q, x->n);
                 acb_printd(s1, 25);
                 flint_printf("\ntheta ", q, x->n);

acb_dirichlet/test/t-thetanull.c

@@ -100,22 +100,26 @@ int main()
                 if (v[k] != ACB_DIRICHLET_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 ((q == 300 && (chi->x->n == 271 || chi->x->n == 131))
+                    || (q == 600 && (chi->x->n == 11 || chi->x->n == 91)))
             {
                 if (!acb_contains_zero(sum))
                 {
-                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n);
+                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->x->n);
                 acb_printd(sum, 10);
                 flint_printf("\n");
+                acb_dirichlet_char_print(G, chi);
+                flint_printf("\n");
                 abort();
                 }
             }
             else if (acb_contains_zero(sum))
             {
-                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->n);
+                flint_printf("FAIL: Theta(chi_%wu(%wu))=", q, chi->x->n);
                 acb_printd(sum, 10);
                 flint_printf("\n");
+                acb_dirichlet_char_print(G, chi);
+                flint_printf("\n");
                 abort();
             }
 

acb_dirichlet/test/t-vec.c

@@ -69,7 +69,7 @@ int main()
             if ((k = vec_diff(v1, v2, nv)))
             {
                 flint_printf("FAIL: chi(%wu,%wu) = %wu != chi(%wu,%wu) = %wu mod %wu (modulus = %wu)\n",
-                        chi->n, k, v1[k], chi->n, k, v2[k], chi->order, q);
+                        chi->x->n, k, v1[k], chi->x->n, k, v2[k], chi->order, q);
                 abort();
             }
 

doc/source/acb_dirichlet.rst

@@ -213,6 +213,13 @@ unity.
 
     There are no restrictions on *n*.
 
+Gauss and Jacobi sums
+-------------------------------------------------------------------------------
+
+.. function:: void acb_dirichlet_gauss_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong prec)
+
+.. function:: void acb_dirichlet_jacobi_sum(acb_t res, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi1,  const acb_dirichlet_char_t chi2, slong prec)
+
 Euler products
 -------------------------------------------------------------------------------
 
@@ -250,3 +257,18 @@ Simple functions
     Note that the alternating character `\{1,-1\}` is not itself
     a Dirichlet character.
 
+
+Implementation notes
+-------------------------------------------------------------------------------
+
+The current implementation introduces a *char* type which contains a *conrey*
+index plus additional information which
+
+- make evaluation of a single character a bit faster
+
+- have some initialization cost.
+
+Even if it is straiforward to convert a *conrey* index to the
+corresponding *char*, looping is faster at the
+level of conrey representation. Things can be improved on this aspect
+but it makes code more intricate.