From 4ece4cbacb6664ea3a4a007d0f68caa73ef41fa3 Mon Sep 17 00:00:00 2001 From: Pascal Date: Thu, 17 Mar 2016 17:27:28 +0100 Subject: [PATCH] split files --- dlog.h | 12 +- dlog/factor_group.c | 63 ++++++++ dlog/profile/p-precomp.c | 24 ++- dlog/profile/p-vec.c | 22 ++- dlog/vec_crt.c | 46 +----- dlog/vec_eratos.c | 28 +--- dlog/vec_eratos_subgroup.c | 52 +++++++ dlog/vec_loop.c | 6 +- dlog/vec_pindex_factorgcd.c | 114 +++++++++++++++ dlog/vec_sieve.c | 281 +----------------------------------- dlog/vec_sieve_subgroup.c | 114 +++++++++++++++ 11 files changed, 399 insertions(+), 363 deletions(-) create mode 100644 dlog/factor_group.c create mode 100644 dlog/vec_eratos_subgroup.c create mode 100644 dlog/vec_pindex_factorgcd.c create mode 100644 dlog/vec_sieve_subgroup.c diff --git a/dlog.h b/dlog.h index a4adbaf2..4de69b7c 100644 --- a/dlog.h +++ b/dlog.h @@ -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); diff --git a/dlog/factor_group.c b/dlog/factor_group.c new file mode 100644 index 00000000..65bb43e8 --- /dev/null +++ b/dlog/factor_group.c @@ -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; +} diff --git a/dlog/profile/p-precomp.c b/dlog/profile/p-precomp.c index c4aa9c1b..56e9eba9 100644 --- a/dlog/profile/p-precomp.c +++ b/dlog/profile/p-precomp.c @@ -28,7 +28,10 @@ #include #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 } } } diff --git a/dlog/profile/p-vec.c b/dlog/profile/p-vec.c index bfedf73a..cdb9fc58 100644 --- a/dlog/profile/p-vec.c +++ b/dlog/profile/p-vec.c @@ -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); } diff --git a/dlog/vec_crt.c b/dlog/vec_crt.c index 09bb4928..918573bb 100644 --- a/dlog/vec_crt.c +++ b/dlog/vec_crt.c @@ -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); } } } diff --git a/dlog/vec_eratos.c b/dlog/vec_eratos.c index fd33374e..7abbf7b6 100644 --- a/dlog/vec_eratos.c +++ b/dlog/vec_eratos.c @@ -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); } diff --git a/dlog/vec_eratos_subgroup.c b/dlog/vec_eratos_subgroup.c new file mode 100644 index 00000000..c17e4b70 --- /dev/null +++ b/dlog/vec_eratos_subgroup.c @@ -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); +} diff --git a/dlog/vec_loop.c b/dlog/vec_loop.c index 66640ea0..3bfcaab3 100644 --- a/dlog/vec_loop.c +++ b/dlog/vec_loop.c @@ -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); } } diff --git a/dlog/vec_pindex_factorgcd.c b/dlog/vec_pindex_factorgcd.c new file mode 100644 index 00000000..0458ab65 --- /dev/null +++ b/dlog/vec_pindex_factorgcd.c @@ -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 + +#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; +} diff --git a/dlog/vec_sieve.c b/dlog/vec_sieve.c index d2263c39..98ee3c3f 100644 --- a/dlog/vec_sieve.c +++ b/dlog/vec_sieve.c @@ -24,288 +24,9 @@ ******************************************************************************/ #include "dlog.h" -#include - -#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); } diff --git a/dlog/vec_sieve_subgroup.c b/dlog/vec_sieve_subgroup.c new file mode 100644 index 00000000..ba0b0b21 --- /dev/null +++ b/dlog/vec_sieve_subgroup.c @@ -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 + +#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); +}