use dlog_precomp in chi_vec

acb_dirichlet/profile/p-gauss.c

@@ -33,7 +33,7 @@ int main(int argc, char *argv[])
         acb_dirichlet_gauss_sum_theta,
         acb_dirichlet_gauss_sum,
         acb_dirichlet_gauss_sum };
-    char * name[5] = { "naive", "factor", "theta", "default", "default+log" };
+    char * name[5] = { "naive", "factor", "theta", "default", "precomp" };
 
     int i, ni = 5;
     ulong qmin[5] =  {  3,  50, 500, 1000, 10000 };
@@ -96,7 +96,7 @@ int main(int argc, char *argv[])
 
                     acb_dirichlet_group_init(G, q);
                     if (l == 4)
-                        acb_dirichlet_group_dlog_precompute(G, q);
+                        acb_dirichlet_group_dlog_precompute(G, 1);
                     acb_dirichlet_char_init(chi, G);
                     acb_dirichlet_char_first_primitive(chi, G);
                     acb_init(res);

acb_dirichlet/profile/p-vec.c

@@ -42,7 +42,7 @@ do_dlog_primeloop(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_c
 }
 
 static void
-do_sieve(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
+do_eratos(ulong *v, const acb_dirichlet_group_t G, const acb_dirichlet_char_t chi, slong nv)
 {
 	slong k, p, pmax;
 	n_primes_t iter;
@@ -81,12 +81,16 @@ int main(int argc, char *argv[])
     flint_rand_t state;
     slong r, nr;
 
-    int l, nf = 6;
-    do_f func[6] = { do_empty, acb_dirichlet_ui_chi_vec_loop, do_dlog_primeloop,
-        acb_dirichlet_ui_chi_vec_primeloop, do_sieve,
+    int l, nf = 9;
+    do_f func[9] = { do_empty, acb_dirichlet_ui_chi_vec_loop, do_dlog_primeloop,
+        acb_dirichlet_ui_chi_vec_primeloop, do_eratos,
+        acb_dirichlet_ui_chi_vec,
+        acb_dirichlet_ui_chi_vec,
+        acb_dirichlet_ui_chi_vec,
         acb_dirichlet_ui_chi_vec };
-    char * name[6] = { "char only", "big loop", "prime loops",
-        "prime dlog_vec", "manual sieve", "default" };
+    char * name[9] = { "char only", "big loop", "prime loops",
+        "prime dlog_vec", "manual eratos", "default",
+        "precomp 1", "precomp 20", "precomp 100" };
 
     int i, ni = 5;
     ulong qmin[5] =  {   2, 1000, 3000, 10000, 100000 };
@@ -131,12 +135,16 @@ int main(int argc, char *argv[])
         {
 
             if (out == LOG)
-                flint_printf("%wu * chi(rand, 1..%wu) for all %wu <= q <= %wu....\n", nr, nv[j], qmin[i], qmax[i]);
+                flint_printf("%wu * ui_chi(rand, 1..%wu) for all %wu <= q <= %wu....\n", nr, nv[j], qmin[i], qmax[i]);
 
             for (l = 0; l < nf; l++)
             {
                 ulong q;
 
+                /* eratos too slow */
+                if (l == 4 && i > 2)
+                    continue;
+
                 if (out == LOG)
                     flint_printf("%-14s ...  ", name[l]);
                 else if (out == CSV)
@@ -155,12 +163,18 @@ int main(int argc, char *argv[])
                     acb_dirichlet_group_init(G, q);
                     acb_dirichlet_char_init(chi, G);
 
+                    if (l >= 6)
+                        acb_dirichlet_group_dlog_precompute(G, (l == 6) ? 1 : (l==7) ? 20 : 100);
+
                     for (r = 0; r < nr; r++)
                     {
                         acb_dirichlet_char(chi, G, rand[r] % q);
                         func[l](v, G, chi, nv[j]);
                     }
 
+                    if (l >= 6)
+                        acb_dirichlet_group_dlog_clear(G);
+
                     acb_dirichlet_char_clear(chi);
                     acb_dirichlet_group_clear(G);
                 }

acb_dirichlet/ui_chi_vec_primeloop.c

@@ -62,7 +62,10 @@ acb_dirichlet_ui_chi_vec_primeloop(ulong *v, const acb_dirichlet_group_t G, cons
         acb_dirichlet_prime_group_struct P = G->P[l];
 
         /* FIXME: there may be some precomputed dlog in P if needed */
-        dlog_vec_add(v, nv, P.g, chi->expo[l], P.pe, P.phi, chi->order);
+        if (P.dlog == NULL)
+            dlog_vec_add(v, nv, P.g, chi->expo[l], P.pe, P.phi, chi->order);
+        else
+            dlog_vec_add_precomp(v, nv, P.dlog, P.g, chi->expo[l], P.pe, P.phi, chi->order);
 
     }
     acb_dirichlet_ui_vec_set_null(v, G, nv);

dlog.h

@@ -259,4 +259,9 @@ void dlog_vec_sieve(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na
 void dlog_vec_add(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
 void dlog_vec(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
 
+
+void dlog_vec_sieve_precomp(ulong *v, ulong nv, dlog_precomp_t pre,  ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
+void dlog_vec_sieve_add_precomp(ulong *v, ulong nv, dlog_precomp_t pre, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
+void dlog_vec_add_precomp(ulong *v, ulong nv, dlog_precomp_t pre, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
+
 #endif

dlog/vec_add_precomp.c

@@ -0,0 +1,23 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "dlog.h"
+
+void
+dlog_vec_add_precomp(ulong *v, ulong nv, dlog_precomp_t pre, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order)
+{
+    if (va == 0)
+        return;
+    if (na * DLOG_LOOP_MAX_FACTOR < nv)
+        dlog_vec_loop_add(v, nv, a, va, mod, na, order);
+    else
+        dlog_vec_sieve_add_precomp(v, nv, pre, a, va, mod, na, order);
+}

dlog/vec_sieve.c

@@ -18,80 +18,10 @@
 void
 dlog_vec_sieve(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order)
 {
-    ulong smooth = 0, sievecount = 0, logcount = 0, missed = 0;
-    ulong logcost;
-#if 0
-    ulong limcount;
-#endif
-    ulong k, p, p1, pmax, logm1;
+    ulong p1 = 50; /* FIXME: tune this limit! */
     dlog_precomp_t pre;
-    n_primes_t iter;
-    ulong X, aX, vaX;
 
-    dlog_vec_fill(v, nv, DLOG_NOT_FOUND);
-    v[1] = 0;
-
-    logm1 = (na % 2) ? 0 : nmod_mul(na / 2, va, order);
-
-    /* discrete log on first primes, then sieve */
-    pmax = (nv < mod.n) ? nv : mod.n;
-    p1 = 50; /* FIXME: tune this limit! */
     dlog_precomp_n_init(pre, a, mod.n, na, p1);
-    logcost = pre->cost;
-
-#if 0
-    if (logcost < 15)
-    {
-        /* p1 = pmax; */
-        limcount = mod.n;
-    }
-    else
-    {
-        limcount = ceil(pow((double)mod.n,1./2.3) * 40 / logcost);
-    }
-#endif
-
-    /* take big power of gen */
-    X = n_nextprime(3 * na / 2, 0) % na;
-    aX = nmod_pow_ui(a, X, mod);
-    vaX = nmod_mul(va, X % order.n, order);
-
-    n_primes_init(iter);
-    while ((p = n_primes_next(iter)) < pmax)
-    {
-        double cost;
-        ulong m, vp;
-        if (mod.n % p == 0)
-            continue; /* won't be attained another time */
-        cost = log(mod.n)/log(p);
-        cost = pow(cost,cost);
-        sievecount++;
-        /* if (p < p1 || (wp = logp_sieve(w, nv, p, mod.n, logm1, order, logcost)) == NOT_FOUND) */
-        /* if (smooth < limcount || (wp = logp_sieve_factor(w, nv, p, mod.n, a, na, va, logm1, order, logcost)) == NOT_FOUND)*/
-        if (logcost < cost || (vp = dlog_vec_pindex_factorgcd(v, nv, p, mod, aX, na, vaX, logm1, order, cost)) == DLOG_NOT_FOUND)
-        {
-            if (logcost < cost)
-                sievecount--;
-            else
-                missed++;
-            logcount++;
-            vp = nmod_mul(dlog_precomp(pre, p), va, order);
-        }
-        for (k = p, m = 1; k < nv; k += p, m++)
-        {
-            if (v[m] == DLOG_NOT_FOUND)
-                continue;
-            smooth++;
-            v[k] = nmod_add(v[m],  vp, order);
-        }
-    }
-#if vbs
-    if (missed)
-        flint_printf("[sieve: got %wu / %wu, n = %wu, cost %wu, logs %wu, sieve %wu missed %wu]\n",
-                smooth, limcount, mod.n, logcost, logcount, sievecount, missed);
-#endif
-    n_primes_clear(iter);
+    dlog_vec_sieve_precomp(v, nv, pre, a, va, mod, na, order);
     dlog_precomp_clear(pre);
-    for (k = mod.n + 1; k < nv; k++)
-        v[k] = v[k - mod.n];
 }

dlog/vec_sieve_add_precomp.c

@@ -0,0 +1,26 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "dlog.h"
+
+void
+dlog_vec_sieve_add_precomp(ulong *v, ulong nv, dlog_precomp_t pre, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order)
+{
+    ulong * w, k;
+    /* store size */
+    w = flint_malloc(nv * sizeof(ulong));
+    dlog_vec_sieve_precomp(w, nv, pre, a, va, mod, na, order);
+    /* write in v */
+    for (k = 0; k < nv; k++)
+        if (v[k] != DLOG_NOT_FOUND)
+            v[k] = nmod_add(v[k], w[k], order);
+    flint_free(w);
+}

dlog/vec_sieve_precomp.c

@@ -0,0 +1,95 @@
+/*
+    Copyright (C) 2016 Pascal Molin
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "dlog.h"
+#include <math.h>
+
+#define vbs 0
+
+/* TODO: tune the limit dlog -> index calculus */
+void
+dlog_vec_sieve_precomp(ulong *v, ulong nv, dlog_precomp_t pre,  ulong a, ulong va, nmod_t mod, ulong na, nmod_t order)
+{
+    ulong smooth = 0, sievecount = 0, logcount = 0, missed = 0;
+    ulong logcost;
+#if 0
+    ulong limcount;
+#endif
+    ulong k, p, pmax, logm1;
+    n_primes_t iter;
+    ulong X, aX, vaX;
+
+    dlog_vec_fill(v, nv, DLOG_NOT_FOUND);
+    v[1] = 0;
+
+    logm1 = (na % 2) ? 0 : nmod_mul(na / 2, va, order);
+
+    /* discrete log on first primes, then sieve */
+    pmax = (nv < mod.n) ? nv : mod.n;
+    logcost = pre->cost;
+
+#if 0
+    if (logcost < 15)
+    {
+        /* p1 = pmax; */
+        limcount = mod.n;
+    }
+    else
+    {
+        limcount = ceil(pow((double)mod.n,1./2.3) * 40 / logcost);
+    }
+#endif
+
+    /* take big power of gen */
+    X = n_nextprime(3 * na / 2, 0) % na;
+    aX = nmod_pow_ui(a, X, mod);
+    vaX = nmod_mul(va, X % order.n, order);
+
+    n_primes_init(iter);
+    while ((p = n_primes_next(iter)) < pmax)
+    {
+        double cost;
+        ulong m, vp;
+        if (mod.n % p == 0)
+            continue; /* won't be attained another time */
+        cost = log(mod.n)/log(p);
+        cost = pow(cost,cost);
+        sievecount++;
+        /* if (p < p1 || (wp = logp_sieve(w, nv, p, mod.n, logm1, order, logcost)) == NOT_FOUND) */
+        /* if (smooth < limcount || (wp = logp_sieve_factor(w, nv, p, mod.n, a, na, va, logm1, order, logcost)) == NOT_FOUND)*/
+        if (logcost < cost || (vp = dlog_vec_pindex_factorgcd(v, nv, p, mod, aX, na, vaX, logm1, order, cost)) == DLOG_NOT_FOUND)
+        {
+            if (logcost < cost)
+                sievecount--;
+            else
+                missed++;
+            logcount++;
+            vp = nmod_mul(dlog_precomp(pre, p), va, order);
+        }
+        for (k = p, m = 1; k < nv; k += p, m++)
+        {
+            if (v[m] == DLOG_NOT_FOUND)
+                continue;
+            smooth++;
+            v[k] = nmod_add(v[m],  vp, order);
+        }
+    }
+#if vbs
+    if (missed)
+        flint_printf("[sieve: got %wu / %wu, n = %wu, cost %wu, logs %wu, sieve %wu missed %wu]\n",
+                smooth, limcount, mod.n, logcost, logcount, sievecount, missed);
+#endif
+    n_primes_clear(iter);
+
+    for (k = mod.n + 1; k < nv; k++)
+        v[k] = v[k - mod.n];
+
+}