diff --git a/acb_dirichlet/profile/p-gauss.c b/acb_dirichlet/profile/p-gauss.c
index 2eb63063..d1ab2ee0 100644
--- a/acb_dirichlet/profile/p-gauss.c
+++ b/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);
diff --git a/acb_dirichlet/profile/p-vec.c b/acb_dirichlet/profile/p-vec.c
index 3ce17a02..098bfdd3 100644
--- a/acb_dirichlet/profile/p-vec.c
+++ b/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);
}
diff --git a/acb_dirichlet/ui_chi_vec_primeloop.c b/acb_dirichlet/ui_chi_vec_primeloop.c
index f05cc408..7e1311da 100644
--- a/acb_dirichlet/ui_chi_vec_primeloop.c
+++ b/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);
diff --git a/dlog.h b/dlog.h
index 8a4805ec..3a475044 100644
--- a/dlog.h
+++ b/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
diff --git a/dlog/vec_add_precomp.c b/dlog/vec_add_precomp.c
new file mode 100644
index 00000000..9bcfeeac
--- /dev/null
+++ b/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 .
+*/
+
+#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);
+}
diff --git a/dlog/vec_init.c b/dlog/vec_fill.c
similarity index 100%
rename from dlog/vec_init.c
rename to dlog/vec_fill.c
diff --git a/dlog/vec_sieve.c b/dlog/vec_sieve.c
index eb877c1f..693adcdc 100644
--- a/dlog/vec_sieve.c
+++ b/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];
}
diff --git a/dlog/vec_sieve_add_precomp.c b/dlog/vec_sieve_add_precomp.c
new file mode 100644
index 00000000..b753a0b7
--- /dev/null
+++ b/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 .
+*/
+
+#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);
+}
diff --git a/dlog/vec_sieve_precomp.c b/dlog/vec_sieve_precomp.c
new file mode 100644
index 00000000..71b405b6
--- /dev/null
+++ b/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 .
+*/
+
+#include "dlog.h"
+#include
+
+#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];
+
+}