hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

split files

Browse source
This commit is contained in:
Pascal 2016-03-17 17:27:28 +01:00
parent 5ecc8e1f25
commit 4ece4cbacb
11 changed files with 399 additions and 363 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

12
dlog.h
View file

@ -197,17 +197,23 @@ void dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong
void dlog_precomp_clear(dlog_precomp_t pre);
ulong dlog_precomp(const dlog_precomp_t pre, ulong b);
#define NOT_FOUND UWORD_MAX
#define LOOP_MAX_FACTOR 6
#define G_SMALL 0
#define G_BIG 1
void dlog_n_factor_group(n_factor_t * fac, ulong bound);
#define NOT_FOUND UWORD_MAX
typedef struct
{
ulong m, logm;
} log_pair_t;
ulong dlog_vec_pindex_factorgcd(log_pair_t * v, ulong nv, ulong p, nmod_t mod, ulong a, ulong na, ulong loga, ulong logm1, nmod_t order, int maxtry);
void dlog_vec_loop(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_eratos_ph(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_eratos_subgroup(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_eratos(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_sieve_ph(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_sieve_subgroup(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_sieve(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order);
void dlog_vec_crt(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);

63
dlog/factor_group.c Normal file
View file

@ -0,0 +1,63 @@
/*=============================================================================
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 "dlog.h"
/* group components up to bound and return cofactor */
void
dlog_n_factor_group(n_factor_t * fac, ulong bound)
{
int i, j, k;
ulong m[FLINT_MAX_FACTORS_IN_LIMB];
ulong c = 1;
i = 0;
for (k = fac->num - 1; k >= 0; k--)
{
ulong qe = n_pow(fac->p[k], fac->exp[k]);
if (qe > bound)
c *= qe;
else
{
/* try to insert somewhere in m */
for (j = 0; j < i && (m[j] * qe > bound); j++);
if (j == i)
m[i++] = qe;
else
m[j] *= qe;
}
}
for (j = 0; j < i; j++)
{
fac->p[j] = m[j];
fac->exp[j] = G_SMALL;
}
if (c > 1)
{
fac->p[i] = c;
fac->exp[i] = G_BIG;
i++;
}
fac->num = i;
}

24
dlog/profile/p-precomp.c
View file

@ -28,7 +28,10 @@
#include <math.h>
#define NUMPRIMES 400
#define CSV 0
#define LOG 0
#define CSV 1
#define JSON 2
#define OUT JSON
typedef void (*log_f) (ulong p, ulong a, ulong num);
@ -133,7 +136,7 @@ int main()
for (i = 0; i < nl; i++)
{
int f;
#if LOG
#if OUT == LOG
flint_printf("%d * logs mod primes of size %d.....\n", l[i], nbits);
#endif
@ -147,17 +150,26 @@ int main()
continue;
if (f == 1 && nbits >= 30 && l[i] > 10)
continue;
#if LOG
#if OUT == LOG
flint_printf("%-20s... ",n[f]);
fflush(stdout);
#elif OUT == CSV
flint_printf("%-8s, %2d, %4d, %3d, ",n[f],nbits,l[i],np);
#elif OUT == JSON
flint_printf("{ \"name\": \"%s\", \"bits\": %d, \"nlogs\": %d, \"nprimes\": %d, \"time\": ",
n[f],nbits,l[i],np);
#endif
TIMEIT_ONCE_START
for (j = 0; j < np; j ++)
(func[f])(p[j], a[j], l[i]);
#if LOG == 0
flint_printf("%-8s, %2d, %4d, %3d, ",n[f],nbits,l[i],np);
#endif
TIMEIT_ONCE_STOP
#if OUT == JSON
flint_printf("}\n");
#else
flint_printf("\n");
#endif
}
}
}

22
dlog/profile/p-vec.c
View file

@ -27,6 +27,11 @@
#include "profiler.h"
#define NPRIMES 640
#define LOG 0
#define CSV 1
#define JSON 2
#define OUT JSON
typedef void (*vec_f) (ulong *v, ulong nv, ulong a, ulong va, const nmod_t mod, ulong na, const nmod_t order);
@ -60,7 +65,6 @@ int main()
p = flint_malloc(np * sizeof(nmod_t));
a = flint_malloc(np * sizeof(ulong));
flint_randinit(state);
for (i = 0; i < ni; i++)
@ -76,8 +80,10 @@ int main()
v = flint_malloc(nv[j] * sizeof(ulong));
#if OUT == LOG
flint_printf("log(1..%wu) mod %d primes of size %d bits....\n", nv[j], np, bits[i]);
fflush(stdout);
#endif
for (l = 0; l < nf; l++)
{
@ -86,8 +92,16 @@ int main()
if (l == 2 && i > 5)
continue;
#if OUT == LOG
flint_printf("%-20s... ",n[l]);
fflush(stdout);
#elif OUT == CSV
flint_printf("%-8s, %2d, %4d, %3d, ",n[l],bits[i],nv[j],np);
#elif OUT == JSON
flint_printf("{ \"name\": \"%s\", \"bits\": %d, \"nv\": %d, \"nprimes\": %d, \"time\": ",
n[l],bits[i],nv[j],np);
#endif
TIMEIT_ONCE_START
for (k = 0; k < np; k++)
{
@ -97,6 +111,12 @@ int main()
(func[l])(v, nv[j], a[k], 1, p[k], p[k].n - 1, order);
}
TIMEIT_ONCE_STOP
#if OUT == JSON
flint_printf("}\n");
#else
flint_printf("\n");
#endif
}
flint_free(v);
}

46
dlog/vec_crt.c
View file

@ -27,46 +27,6 @@
#define SIEVE_START 100
/* group components up to bound and return cofactor */
#define G_SMALL 0
#define G_BIG 1
static void
n_factor_group(n_factor_t * fac, ulong bound)
{
int i, j, k;
ulong m[FLINT_MAX_FACTORS_IN_LIMB];
ulong c = 1;
i = 0;
for (k = fac->num - 1; k >= 0; k--)
{
ulong qe = n_pow(fac->p[k], fac->exp[k]);
if (qe > bound)
c *= qe;
else
{
/* try to insert somewhere in m */
for (j = 0; j < i && (m[j] * qe > bound); j++);
if (j == i)
m[i++] = qe;
else
m[j] *= qe;
}
}
for (j = 0; j < i; j++)
{
fac->p[j] = m[j];
fac->exp[j] = G_SMALL;
}
if (c > 1)
{
fac->p[i] = c;
fac->exp[i] = G_BIG;
i++;
}
fac->num = i;
}
/* assume v[k] = -1 for bad primes? */
/* loop on small components and if needed keep one subgroup for DLOG + sieve */
void
@ -79,7 +39,7 @@ dlog_vec_crt(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t
maxloop = LOOP_MAX_FACTOR * nv;
n_factor_init(&fac);
n_factor(&fac, na, 1);
n_factor_group(&fac, maxloop);
dlog_n_factor_group(&fac, maxloop);
for (k = 0; k < fac.num; k++)
{
ulong m, M, aM, uM, vaM;
@ -93,9 +53,9 @@ dlog_vec_crt(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t
else
{
if (nv <= SIEVE_START)
dlog_vec_eratos_ph(v, nv, aM, vaM, M, mod, m, order);
dlog_vec_eratos_subgroup(v, nv, aM, vaM, M, mod, m, order);
else
dlog_vec_sieve_ph(v, nv, aM, vaM, M, mod, m, order);
dlog_vec_sieve_subgroup(v, nv, aM, vaM, M, mod, m, order);
}
}
}

28
dlog/vec_eratos.c
View file

@ -25,34 +25,8 @@
#include "dlog.h"
void
dlog_vec_eratos_ph(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order)
{
ulong p, pmax;
dlog_precomp_t pre;
n_primes_t iter;
/* discrete log on primes */
pmax = (nv < mod.n) ? nv : mod.n;
dlog_precomp_n_init(pre, a, mod.n, na, n_prime_pi(nv));
n_primes_init(iter);
while ((p = n_primes_next(iter)) < pmax)
{
ulong k, pM, wp;
if (mod.n % p == 0)
continue; /* won't be attained another time */
pM = (M) ? nmod_pow_ui(p, M, mod) : p;
wp = nmod_mul(dlog_precomp(pre, pM), va, order);
for (pM = p; pM < nv; pM *= p)
for (k = pM; k < nv; k += pM)
v[k] = nmod_add(v[k], wp, order);
}
n_primes_clear(iter);
}
void
dlog_vec_eratos(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order)
{
dlog_vec_eratos_ph(v, nv, a, va, 1, mod, na, order);
dlog_vec_eratos_subgroup(v, nv, a, va, 1, mod, na, order);
}

52
dlog/vec_eratos_subgroup.c Normal file
View file

@ -0,0 +1,52 @@
/*=============================================================================
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 "dlog.h"
void
dlog_vec_eratos_subgroup(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order)
{
ulong p, pmax;
dlog_precomp_t pre;
n_primes_t iter;
/* discrete log on primes */
pmax = (nv < mod.n) ? nv : mod.n;
dlog_precomp_n_init(pre, a, mod.n, na, n_prime_pi(nv));
n_primes_init(iter);
while ((p = n_primes_next(iter)) < pmax)
{
ulong k, pM, wp;
if (mod.n % p == 0)
continue; /* won't be attained another time */
pM = (M) ? nmod_pow_ui(p, M, mod) : p;
wp = nmod_mul(dlog_precomp(pre, pM), va, order);
for (pM = p; pM < nv; pM *= p)
for (k = pM; k < nv; k += pM)
v[k] = nmod_add(v[k], wp, order);
}
n_primes_clear(iter);
}

6
dlog/vec_loop.c
View file

@ -33,8 +33,8 @@ dlog_vec_loop(ulong * v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod
ulong vx = 0;
for (x = a; x != 1; x = nmod_mul(x, a, mod))
{
vx = nmod_add(vx, va, order);
for(xp = x; xp < nv; xp+=mod.n)
v[xp] = nmod_add(v[xp], vx, order);
vx = nmod_add(vx, va, order);
for(xp = x; xp < nv; xp+=mod.n)
v[xp] = nmod_add(v[xp], vx, order);
}
}

114
dlog/vec_pindex_factorgcd.c Normal file
View file

@ -0,0 +1,114 @@
/*=============================================================================
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 "dlog.h"
#include <math.h>
#define vbs 0
#define FACTOR_RATIO 4
static int
factor_until(ulong * n, ulong nlim, const ulong * p, ulong pmax, ulong * fp, int * fe)
{
int i, j;
for (i = 0, j = 0; *n >= nlim && p[j] < pmax; j++)
{
int e = n_remove(n, p[j]);
if (e)
{
fp[i] = p[j];
fe[i] = e;
i++;
}
}
return i;
}
ulong
dlog_vec_pindex_factorgcd(log_pair_t * v, ulong nv, ulong p, nmod_t mod, ulong a, ulong na, ulong loga, ulong logm1, nmod_t order, int maxtry)
{
int nm = 0, ng = 0;
ulong pm, logm, pmax;
ulong u[2], r[2], t;
ulong up[15], rp[15];
int ue[15], re[15];
const ulong * prime;
prime = n_primes_arr_readonly(p);
pmax = p / FACTOR_RATIO;
pm = p;
logm = 0;
while (nm++ < maxtry)
{
int i, j, iu, ir;
ulong logr;
pm = nmod_mul(pm, a, mod);
logm = nmod_add(logm, loga, order);
/*
if (2 * pm > mod.n)
{
pm = nmod_neg(pm, mod);
logm = nmod_add(logm, logm1, order);
}
*/
/* half gcd u * pm + v * mod = r, ignore v */
u[0] = 0; r[0] = mod.n;
u[1] = 1; r[1] = pm;
i = 1; j = 0; /* flip flap */
while (r[i] > u[i])
{
ng++;
if (r[i] < nv && v[r[i]].m == r[i] && u[i] < nv && v[u[i]].m == u[i])
{
/* early smooth detection: occurs for primes < 30 bits */
ulong x;
/* chi(-1)^j*chi(u)*chi(p)*chi(m)=chi(r) */
x = nmod_sub(v[r[i]].logm, nmod_add(v[u[i]].logm, logm, order), order);
if (j)
x = nmod_add(x, logm1, order);
return x;
}
j = i; i = 1 - i; /* switch */
t = r[i] / r[j];
r[i] = r[i] % r[j];
u[i] = u[i] + t * u[j]; /* (-1)^j */
};
/* try to factor both r[i] and u[i] */
iu = factor_until(&u[i], nv, prime, pmax, up, ue);
if (u[i] >= nv || v[u[i]].m < u[i])
continue;
ir = factor_until(&r[i], nv, prime, pmax, rp, re);
if (r[i] >= nv || v[r[i]].m < r[i])
continue;
/* log(u)+log(p)+log(m)=log(r) */
logr = 0;
for (i=0; i < ir; i++)
logr = nmod_add(logr, (re[i] * v[rp[i]].logm) % order.n, order);
for (i=0; i < iu; i++)
logm = nmod_add(logm, (ue[i] * v[up[i]].logm) % order.n, order);
return nmod_sub(logr, logm, order);
}
return NOT_FOUND;
}

281
dlog/vec_sieve.c
View file

@ -24,288 +24,9 @@
******************************************************************************/
#include "dlog.h"
#include <math.h>
#define vbs 0
#define FACTOR_RATIO 4
static ulong
logp_sieve_gcd(log_pair_t * v, ulong nv, ulong p, nmod_t mod, ulong a, ulong na, ulong loga, ulong logm1, nmod_t order, int maxtry)
{
int i, j, nm = 0, ng = 0;
ulong pm, logm;
ulong u[2], r[2], t;
#if vbs > 1
flint_printf("\nEnter logp_sieve p=%wu mod %wu...\n", p, mod);
#endif
pm = p;
logm = 0;
while (nm++ < maxtry)
{
pm = nmod_mul(pm, a, mod);
logm = nmod_add(logm, loga, order);
if (2 * pm > mod.n)
{
pm = nmod_neg(pm, mod);
logm = nmod_add(logm, logm1, order);
}
/* half gcd u * pm + v * mod = r, ignore v */
u[0] = 0; r[0] = mod.n;
u[1] = 1; r[1] = pm;
i = 1; j = 0; /* flip flap */
do {
ng++;
#if vbs > 1
if (r[i] < nv && u[i] < nv)
flint_printf("[r=%wu, v[r]=%wu, u=%wu, v[u]=%wu -- nv=%wu]\n",
r[i], v[r[i]].m, u[i], v[u[i]].m, nv);
#endif
if (r[i] < nv && v[r[i]].m == r[i] && u[i] < nv && v[u[i]].m == u[i])
{
ulong x;
/* chi(-1)^j*chi(u)*chi(p)*chi(m)=chi(r) */
x = nmod_sub(v[r[i]].logm, nmod_add(v[u[i]].logm, logm, order), order);
if (j)
x = nmod_add(x, logm1, order);
return x;
}
j = i; i = 1 - i; /* switch */
t = r[i] / r[j];
r[i] = r[i] % r[j];
u[i] = u[i] + t * u[j]; /* (-1)^j */
} while (r[i] > 0 && u[i] < nv);
}
return NOT_FOUND;
}
static ulong
logp_sieve_factor(log_pair_t * v, ulong nv, ulong p, nmod_t mod, ulong a, ulong na, ulong loga, ulong logm1, nmod_t order, int maxtry)
{
int nm = 0;
ulong pm, logm;
const ulong * prime;
prime = n_primes_arr_readonly(p);
pm = p;
logm = 0;
while (pm < mod.n)
{
/*nm++;*/ /* init ignored */
pm *= a;
logm = nmod_add(logm, loga, order);
}
pm = pm % mod.n;
do {
int i, j, ind[15], exp[15];
/* find multiplier m */
if (2 * pm > mod.n)
{
pm = nmod_neg(pm, mod);
logm = nmod_add(logm, logm1, order);
}
for (i = 0, j = 0; j < p && pm >= nv && 4 * prime[j] < p; j++)
{
int e = n_remove(&pm, prime[j]);
if (e)
{
ind[i] = j;
exp[i] = e;
i++;
}
}
if (pm < nv && v[pm].m == pm)
{
/* goal! */
ulong x = v[pm].logm;
/* chi(m)*chi(p)=chi(pm)*prod */
for(j = 0; j < i; j++)
x = nmod_add(x, nmod_mul(exp[j], v[prime[ind[j]]].logm, order), order);
x = nmod_sub(logm, x, order);
/*flint_printf("managed %d mults / %d for p=%wu, pm=%wu\n",nm,maxtry,p,pm);*/
return x;
}
nm++;
pm = nmod_mul(pm, a, mod);
logm = nmod_add(logm, loga, order);
} while (nm < maxtry);
return NOT_FOUND;
}
static int
factor_until(ulong * n, ulong nlim, const ulong * p, ulong pmax, ulong * fp, int * fe)
{
int i, j;
for (i = 0, j = 0; *n >= nlim && p[j] < pmax; j++)
{
int e = n_remove(n, p[j]);
if (e)
{
fp[i] = p[j];
fe[i] = e;
i++;
}
}
return i;
}
static ulong
logp_sieve_factorgcd(log_pair_t * v, ulong nv, ulong p, nmod_t mod, ulong a, ulong na, ulong loga, ulong logm1, nmod_t order, int maxtry)
{
int nm = 0, ng = 0;
ulong pm, logm, pmax;
ulong u[2], r[2], t;
ulong up[15], rp[15];
int ue[15], re[15];
const ulong * prime;
prime = n_primes_arr_readonly(p);
pmax = p / FACTOR_RATIO;
pm = p;
logm = 0;
while (nm++ < maxtry)
{
int i, j, iu, ir;
ulong logr;
pm = nmod_mul(pm, a, mod);
logm = nmod_add(logm, loga, order);
/*
if (2 * pm > mod.n)
{
pm = nmod_neg(pm, mod);
logm = nmod_add(logm, logm1, order);
}
*/
/* half gcd u * pm + v * mod = r, ignore v */
u[0] = 0; r[0] = mod.n;
u[1] = 1; r[1] = pm;
i = 1; j = 0; /* flip flap */
while (r[i] > u[i])
{
ng++;
if (r[i] < nv && v[r[i]].m == r[i] && u[i] < nv && v[u[i]].m == u[i])
{
/* early smooth detection: occurs for primes < 30 bits */
ulong x;
/* chi(-1)^j*chi(u)*chi(p)*chi(m)=chi(r) */
x = nmod_sub(v[r[i]].logm, nmod_add(v[u[i]].logm, logm, order), order);
if (j)
x = nmod_add(x, logm1, order);
return x;
}
j = i; i = 1 - i; /* switch */
t = r[i] / r[j];
r[i] = r[i] % r[j];
u[i] = u[i] + t * u[j]; /* (-1)^j */
};
/* try to factor both r[i] and u[i] */
iu = factor_until(&u[i], nv, prime, pmax, up, ue);
if (u[i] >= nv || v[u[i]].m < u[i])
continue;
ir = factor_until(&r[i], nv, prime, pmax, rp, re);
if (r[i] >= nv || v[r[i]].m < r[i])
continue;
/* log(u)+log(p)+log(m)=log(r) */
logr = 0;
for (i=0; i < ir; i++)
logr = nmod_add(logr, (re[i] * v[rp[i]].logm) % order.n, order);
for (i=0; i < iu; i++)
logm = nmod_add(logm, (ue[i] * v[up[i]].logm) % order.n, order);
return nmod_sub(logr, logm, order);
}
return NOT_FOUND;
}
void
dlog_vec_sieve_ph(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order)
{
ulong smooth = 0, sievecount = 0, logcount = 0, missed = 0;
ulong logcost, limcount;
ulong k, p, p1, pmax, logm1;
log_pair_t * w;
dlog_precomp_t pre;
n_primes_t iter;
ulong X, aX, vaX;
/* store size */
w = flint_malloc( nv * sizeof(log_pair_t));
for (k = 0; k < nv; k++)
{
w[k].m = 1;
w[k].logm = 0; /* could be v[k]... */
}
/* 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);
/*flint_printf("## single log cost: %wu\n", pre->cost);*/
logcost = pre->cost;
if (logcost < 15)
{
/* p1 = pmax; */
limcount = mod.n;
}
else
{
limcount = ceil(pow((double)mod.n,1./2.3) * 40 / logcost);
logm1 = (mod.n % 2) ? 0 : dlog_precomp(pre, mod.n - 1);
}
/* find big power of gen */
X = n_nextprime(na / 2 + 10, 0);
X = (na % 257) ? 257 % na : 1031 % na ; /* FIXME! */
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 = log(mod.n)/log(p);
ulong m, wp;
if (mod.n % p == 0) /* FIXME: those primes could be known... */
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 || (wp = logp_sieve_factorgcd(w, nv, p, mod, aX, na, vaX, logm1, order, cost)) == NOT_FOUND)
{
if (logcost < cost)
sievecount--;
else
missed++;
logcount++;
wp = nmod_mul(dlog_precomp(pre, nmod_pow_ui(p, M, mod)), va, order);
}
for (k = p, m = 1; k < nv; k += p, m++)
{
w[k].m = w[m].m * p;
w[k].logm = nmod_add(w[m].logm, wp, order);
if (w[k].m == k)
smooth++;
}
}
/* write in v */
for (k = 0; k < nv; k++)
if (v[k] != NOT_FOUND)
v[k] = nmod_add(v[k], w[k].logm, 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);
flint_free(w);
}
void
dlog_vec_sieve(ulong *v, ulong nv, ulong a, ulong va, nmod_t mod, ulong na, nmod_t order)
{
dlog_vec_sieve_ph(v, nv, a, va, 1, mod, na, order);
dlog_vec_sieve_subgroup(v, nv, a, va, 1, mod, na, order);
}

114
dlog/vec_sieve_subgroup.c Normal file
View file

@ -0,0 +1,114 @@
/*=============================================================================
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 "dlog.h"
#include <math.h>
#define vbs 0
void
dlog_vec_sieve_subgroup(ulong *v, ulong nv, ulong a, ulong va, ulong M, nmod_t mod, ulong na, nmod_t order)
{
ulong smooth = 0, sievecount = 0, logcount = 0, missed = 0;
ulong logcost, limcount;
ulong k, p, p1, pmax, logm1;
log_pair_t * w;
dlog_precomp_t pre;
n_primes_t iter;
ulong X, aX, vaX;
/* store size */
w = flint_malloc( nv * sizeof(log_pair_t));
for (k = 0; k < nv; k++)
{
w[k].m = 1;
w[k].logm = 0; /* could be v[k]... */
}
/* 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);
/*flint_printf("## single log cost: %wu\n", pre->cost);*/
logcost = pre->cost;
if (logcost < 15)
{
/* p1 = pmax; */
limcount = mod.n;
}
else
{
limcount = ceil(pow((double)mod.n,1./2.3) * 40 / logcost);
logm1 = (mod.n % 2) ? 0 : dlog_precomp(pre, mod.n - 1);
}
/* find big power of gen */
X = n_nextprime(na / 2 + 10, 0);
X = (na % 257) ? 257 % na : 1031 % na ; /* FIXME! */
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 = log(mod.n)/log(p);
ulong m, wp;
if (mod.n % p == 0) /* FIXME: those primes could be known... */
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 || (wp = dlog_vec_pindex_factorgcd(w, nv, p, mod, aX, na, vaX, logm1, order, cost)) == NOT_FOUND)
{
if (logcost < cost)
sievecount--;
else
missed++;
logcount++;
wp = nmod_mul(dlog_precomp(pre, nmod_pow_ui(p, M, mod)), va, order);
}
for (k = p, m = 1; k < nv; k += p, m++)
{
w[k].m = w[m].m * p;
w[k].logm = nmod_add(w[m].logm, wp, order);
if (w[k].m == k)
smooth++;
}
}
/* write in v */
for (k = 0; k < nv; k++)
if (v[k] != NOT_FOUND)
v[k] = nmod_add(v[k], w[k].logm, 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);
flint_free(w);
}