fix and test chi_vec_primeloop

acb_dirichlet.h

@@ -96,6 +96,7 @@ void acb_dirichlet_acb_pairing(acb_t res, const acb_dirichlet_group_t G, ulong m
 typedef struct
 {
     ulong q;           /* modulus */
+    /*? acb_dirichet_conrey_struct n; */
     ulong n;           /* number */
     ulong order;       /* order */
     ulong * expo;      /* reduced exponents ( order * log[k] / gcd( ) ) */
@@ -105,7 +106,7 @@ acb_dirichlet_char_struct;
 
 typedef acb_dirichlet_char_struct acb_dirichlet_char_t[1];
 
-ACB_DIRICHLET_INLINE int
+ACB_DIRICHLET_INLINE ulong
 acb_dirichlet_char_order(const acb_dirichlet_char_t chi)
 {
     return chi->order;
@@ -127,7 +128,6 @@ void acb_dirichlet_char_normalize(acb_dirichlet_char_t chi, const acb_dirichlet_
 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);
@@ -146,6 +146,7 @@ void acb_dirichlet_char_vec(acb_t res, const acb_dirichlet_group_t G, const acb_
 
 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);
+ulong acb_dirichlet_theta_length(ulong q, double x, slong prec);
 
 #ifdef __cplusplus
 }

acb_dirichlet/chi_vec_primeloop.c

@@ -25,6 +25,39 @@
 
 #include "acb_dirichlet.h"
 
+void
+chi_vec_evenpart(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
+{
+
+    ulong c3, c4, x;
+
+    if (G->neven >= 1 && (c3 = chi->expo[0]))
+    {
+        for (x = 3; x < nv; x += 4)
+            v[x] = c3;
+    }
+    if (G->neven == 2 && (c4 = chi->expo[1]))
+    {
+        ulong g, pe, vx, xp;
+        vx = c4;
+        pe = G->primepowers[1];
+        g = G->generators[1] % pe;
+
+        for (x = g; x > 1;)
+        {
+
+            for (xp = x; xp < nv; xp += pe)
+                v[xp] = (v[xp] + vx) % chi->order;
+
+            for (xp = pe - x; xp < nv; xp += pe)
+                v[xp] = (v[xp] + vx) % chi->order;
+
+            x = (x * g) % pe;
+            vx = (vx + c4) % chi->order;
+        }
+    }
+}
+
 /* loop over primary components */
 void
 acb_dirichlet_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
@@ -34,25 +67,30 @@ acb_dirichlet_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G, const a
     for (k = 1; k < nv; k++)
         v[k] = 0;
 
-    for (l = 1; l < G->num; l++)
+    if (G->neven)
+        chi_vec_evenpart(v, G, chi, nv);
+
+    for (l = G->neven; l < G->num; l++)
     {
-        long p, pe, g, x, vp, xp;
-        long j, vj;
+        ulong p, pe, g, x, vx, vp, xp;
         p = G->primes[l];
         pe = G->primepowers[l];
         g = G->generators[l] % pe;
-        vj = vp = chi->expo[l];
+        vx = vp = chi->expo[l];
+
+        if (vp == 0)
+            continue;
         /* for each x = g^j mod p^e,
          * set a[x] += j*vp
          * and use periodicity */
-
-        for (j = 1, x = g; x > 1; j++)
+        for (x = g; x > 1;)
         {
-            for (xp = x; xp < nv; xp += pe)
-                v[xp] = (v[xp] + vj) % chi->order;
 
-            x = (x*g) % pe;
-            vj = (vj + vp) % chi->order;
+            for (xp = x; xp < nv; xp += pe)
+                v[xp] = (v[xp] + vx) % chi->order;
+
+            x = (x * g) % pe;
+            vx = (vx + vp) % chi->order;
         }
     }
     acb_dirichlet_vec_set_null(v, G, nv);

acb_dirichlet/test/t-chars.c

@@ -26,6 +26,7 @@
 #include "acb_dirichlet.h"
 
 
+/* order of an element knowing the factorization of a multiple */
 static ulong
 nmod_order_precomp(ulong a, nmod_t mod, ulong expo, n_factor_t fac)
 {

acb_dirichlet/test/t-conrey.c

@@ -34,7 +34,7 @@ int main()
     fflush(stdout);
     flint_randinit(state);
 
-    for (iter = 0; iter < 10000; iter++)
+    for (iter = 0; iter < 3000; iter++)
     {
         acb_dirichlet_group_t G;
         acb_dirichlet_conrey_t x;
@@ -109,7 +109,7 @@ int main()
                if (G->exponents[k] == 1)
                    ref *= p - 2;
                else
-                   ref *= (p * (p -2) + 1) * n_pow(p, G->exponents[k] - 2);
+                   ref *= (p * (p - 2) + 1) * n_pow(p, G->exponents[k] - 2);
            }
         }
 

acb_dirichlet/test/t-thetanull.c

@@ -24,14 +24,10 @@
 ******************************************************************************/
 
 #include "acb_dirichlet.h"
-#include <math.h>
-
-#define PI   3.14159265358
-#define LOG2 0.69314718055
 
 int main()
 {
-    slong prec = 80;
+    slong prec = 64;
     ulong q;
 
     flint_printf("thetanull....");
@@ -43,7 +39,7 @@ int main()
      * characters of moduli 300 and 600 + conjugates
      */
 
-    for (q = 3; q < 800; q ++)
+    for (q = 3; q < 1000; q ++)
     {
         acb_dirichlet_group_t G;
         acb_dirichlet_conrey_t x;
@@ -57,6 +53,7 @@ int main()
         arb_ptr t, kt, tt;
 
         if (q % 4 == 2)
+            /* no primitive character mod q */
             continue;
 
         acb_dirichlet_group_init(G, q);
@@ -69,7 +66,7 @@ int main()
         z = _acb_vec_init(G->expo);
         _acb_vec_set_powers(z, zeta, G->expo, prec);
 
-        nv = ceil(sqrt((double)q * prec * LOG2 / PI)) + 2;
+        nv = acb_dirichlet_theta_length(q, 1, prec);
         v = flint_malloc(nv * sizeof(ulong));
 
         arb_init(eq);
@@ -94,12 +91,9 @@ int main()
             acb_zero(sum);
 
             acb_dirichlet_char_conrey(chi, G, x);
-            acb_dirichlet_chi_vec_loop(v, G, chi, nv);
+            acb_dirichlet_chi_vec(v, G, chi, nv);
 
             m = G->expo / chi->order;
-            /*
-            flint_printf("Theta(chi_%wu(%wu)) (m=%wu)\n", q, chi->n, m);
-            */
             tt = acb_dirichlet_char_parity(chi) ? kt : t;
 
             for (k = 1; k < nv; k++)
@@ -133,6 +127,7 @@ int main()
         acb_clear(zeta);
         acb_clear(sum);
         arb_clear(eq);
+        flint_free(v);
         acb_dirichlet_group_clear(G);
         acb_dirichlet_char_clear(chi);
         acb_dirichlet_conrey_clear(x);

acb_dirichlet/test/t-vec.c

@@ -0,0 +1,90 @@
+/*=============================================================================
+
+    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"
+
+static ulong
+vec_diff(ulong * v, ulong * ref, ulong nv)
+{
+    ulong k;
+    for (k = 1; k < nv; k++)
+        if (ref[k] != v[k])
+            return k;
+    return 0;
+}
+
+int main()
+{
+    ulong q;
+
+    flint_printf("vec....");
+    fflush(stdout);
+
+    for (q = 2; q < 600; q ++)
+    {
+        acb_dirichlet_group_t G;
+        acb_dirichlet_conrey_t x;
+        acb_dirichlet_char_t chi;
+        ulong * v1, * v2, nv, k;
+
+        acb_dirichlet_group_init(G, q);
+        acb_dirichlet_conrey_init(x, G);
+        acb_dirichlet_char_init(chi, G);
+
+        nv = 100;
+        v1 = flint_malloc(nv * sizeof(ulong));
+        v2 = flint_malloc(nv * sizeof(ulong));
+
+        acb_dirichlet_conrey_one(x, G);
+
+        while (1) {
+
+            acb_dirichlet_char_conrey(chi, G, x);
+
+            acb_dirichlet_chi_vec_loop(v1, G, chi, nv);
+            acb_dirichlet_chi_vec_primeloop(v2, G, chi, nv);
+
+            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);
+                abort();
+            }
+
+            if (acb_dirichlet_conrey_next(x, G) == G->num)
+                break;
+        }
+
+        flint_free(v1);
+        flint_free(v2);
+        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/theta_length.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"
+#include <math.h>
+#define PI   3.14159265358
+#define LOG2 0.69314718055
+
+ulong
+acb_dirichlet_theta_length(ulong q, double x, slong prec)
+{
+    double a = PI / (double)q * x * x;
+    return ceil(sqrt(((double)prec * LOG2 - log(2 * a)) / a));
+}