1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
/*
Copyright (C) 2016 Pascal Molin
This file is part of FLINT.
FLINT 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 3 of the License, or
(at your option) any later version. See <https://www.gnu.org/licenses/>.
*/
#include <math.h>
#include "nmod.h"
#include "dlog.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)
{
#if vbs
ulong smooth = 0, sievecount = 0, logcount = 0, missed = 0;
#endif
#if 0
ulong limcount;
#endif
ulong logcost;
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_ui_mul_ui(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_ui_mul_ui(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);
#if vbs
sievecount++;
#endif
/* 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 vbs
if (logcost < cost)
sievecount--;
else
missed++;
logcount++;
#endif
vp = nmod_ui_mul_ui(dlog_precomp(pre, p), va, order);
}
for (k = p, m = 1; k < nv; k += p, m++)
{
if (v[m] == DLOG_NOT_FOUND)
continue;
#if vbs
smooth++;
#endif
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];
}