From 32a1c966b90cec3d9e3b22c20b0c5739a231509f Mon Sep 17 00:00:00 2001 From: Pascal Date: Tue, 15 Mar 2016 15:39:26 +0100 Subject: [PATCH] fix log table for large moduli + add order 2 & 3 --- dlog.h | 22 +++++++--- dlog/bsgs.c | 6 +-- dlog/bsgs_init.c | 7 +-- dlog/crt_init.c | 12 ++++- dlog/modpe.c | 3 +- dlog/order23.c | 32 ++++++++++++++ dlog/order23_init.c | 39 +++++++++++++++++ dlog/power_init.c | 5 ++- dlog/precomp.c | 3 ++ dlog/precomp_init.c | 74 +++++++++++++++++++++---------- dlog/profile/p-precomp.c | 95 ++++++++++++++++++++-------------------- dlog/table_init.c | 3 +- 12 files changed, 214 insertions(+), 87 deletions(-) create mode 100644 dlog/order23.c create mode 100644 dlog/order23_init.c diff --git a/dlog.h b/dlog.h index ce63e62d..a4adbaf2 100644 --- a/dlog.h +++ b/dlog.h @@ -32,7 +32,7 @@ enum { - DLOG_MODPE, DLOG_CRT, DLOG_POWER, DLOG_BSGS, DLOG_TABLE + DLOG_MODPE, DLOG_CRT, DLOG_POWER, DLOG_BSGS, DLOG_TABLE, DLOG_23 }; typedef struct dlog_precomp_struct dlog_precomp_struct; @@ -139,12 +139,15 @@ dlog_power_struct; typedef dlog_power_struct dlog_power_t[1]; +typedef ulong dlog_order23_t[1]; + /* generic decomposition */ /*typedef */ struct dlog_precomp_struct { int type; + ulong cost; union { dlog_table_t table; @@ -152,20 +155,23 @@ dlog_precomp_struct dlog_crt_t crt; dlog_power_t power; dlog_modpe_t modpe; + dlog_order23_t order23; } t; }; typedef dlog_precomp_struct dlog_precomp_t[1]; -void dlog_table_init(dlog_table_t t, ulong a, ulong mod); +ulong dlog_order23_init(dlog_order23_t t, ulong a); +ulong dlog_table_init(dlog_table_t t, ulong a, ulong mod); +ulong dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num); +ulong dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num); +ulong dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num); +ulong dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m); void dlog_1modpe_init(dlog_1modpe_t t, ulong a1, ulong p, ulong e); -void dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num); -void dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num); -void dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num); -void dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m); void dlog_rho_init(dlog_rho_t t, ulong a, ulong mod, ulong n); /*#define dlog_bsgs_init(t, a, n, m) bsgs_table_init(t, a, n, m)*/ +void dlog_order23_clear(dlog_order23_t t); void dlog_table_clear(dlog_table_t t); void dlog_1modpe_clear(dlog_1modpe_t t); void dlog_crt_clear(dlog_crt_t t); @@ -175,6 +181,7 @@ void dlog_bsgs_clear(dlog_bsgs_t t); void dlog_rho_clear(dlog_rho_t t); /*#define dlog_bsgs_clear(t) bsgs_table_clear(t)*/ +ulong dlog_order23(const dlog_order23_t t, ulong b); ulong dlog_table(const dlog_table_t t, ulong b); ulong dlog_crt(const dlog_crt_t t, ulong b); ulong dlog_power(const dlog_power_t t, ulong b); @@ -198,6 +205,9 @@ typedef struct } log_pair_t; 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(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(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/bsgs.c b/dlog/bsgs.c index b966814e..0f85c679 100644 --- a/dlog/bsgs.c +++ b/dlog/bsgs.c @@ -33,7 +33,7 @@ dlog_bsgs(const dlog_bsgs_t t, ulong b) apow_t c, * x; c.ak = b; - for (i = 0; i <= t->g; i++) + for (i = 0; i < t->g; i++) { x = bsearch(&c, t->table, t->m, sizeof(apow_t), (int(*)(const void*,const void*))apow_cmp); @@ -42,7 +42,7 @@ dlog_bsgs(const dlog_bsgs_t t, ulong b) c.ak = nmod_mul(c.ak, t->am, t->mod); } flint_printf("Exception (n_discrete_log_bsgs). discrete log not found.\n"); - flint_printf(" table size %wu, cosize %wu mod %wu. %wu not found\n", - t->m, t->g, t->mod.n, b); + flint_printf(" table size %wu, cosize %wu mod %wu. %wu not found (a^-m=%wu)\n", + t->m, t->g, t->mod.n, b, t->am); abort(); } diff --git a/dlog/bsgs_init.c b/dlog/bsgs_init.c index 2a0c68a1..42fb2d19 100644 --- a/dlog/bsgs_init.c +++ b/dlog/bsgs_init.c @@ -34,16 +34,16 @@ apow_cmp(const apow_t * x, const apow_t * y) } /* set size of table m=sqrt(nk) to compute k logs in a group of size n */ -void +ulong dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m) { ulong k, ak; + if (m >= n) m = n + 1; t->table = (apow_t *)flint_malloc(m * sizeof(apow_t)); - if (m > n) m = n; nmod_init(&t->mod, mod); t->m = m; - t->g = n / m; + t->g = n / m + 1; for (k = 0, ak = 1; k < m; k++) { @@ -54,6 +54,7 @@ dlog_bsgs_init(dlog_bsgs_t t, ulong a, ulong mod, ulong n, ulong m) t->am = nmod_inv(ak, t->mod); qsort(t->table, m, sizeof(apow_t), (int(*)(const void*,const void*))apow_cmp); + return t->g; } void diff --git a/dlog/crt_init.c b/dlog/crt_init.c index 17d41da8..381e247f 100644 --- a/dlog/crt_init.c +++ b/dlog/crt_init.c @@ -25,12 +25,13 @@ #include "dlog.h" -void +ulong dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num) { int k; n_factor_t fac; ulong * M, * u; + ulong cost = 0; n_factor_init(&fac); n_factor(&fac, n, 1); @@ -54,8 +55,17 @@ dlog_crt_init(dlog_crt_t t, ulong a, ulong mod, ulong n, ulong num) M[k] = n / mk; u[k] = nmod_mul(M[k], n_invmod(M[k] % mk, mk), t->n); /* depends on the power */ +#if 0 + flint_printf("[sub-crt -- init for size %wu mod %wu]\n", mk, mod); +#endif dlog_precomp_pe_init(t->pre[k], nmod_pow_ui(a, M[k], t->mod), mod, p, e, mk, num); + cost += t->pre[k]->cost; } +#if 0 + if (cost > 500) + flint_printf("[crt init for size %wu mod %wu -> cost %wu]\n", n,mod,cost); +#endif + return cost; } void diff --git a/dlog/modpe.c b/dlog/modpe.c index a41095f2..1941d44b 100644 --- a/dlog/modpe.c +++ b/dlog/modpe.c @@ -26,7 +26,7 @@ #include "dlog.h" #include "padic.h" -void +ulong dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num) { t->p = p; @@ -34,6 +34,7 @@ dlog_modpe_init(dlog_modpe_t t, ulong a, ulong p, ulong e, ulong pe, ulong num) t->pe = pe; dlog_precomp_n_init(t->modp, a, p, p-1, num); dlog_1modpe_init(t->modpe, n_powmod(a, p - 1, pe), p, e); + return t->modp->cost + e; } void diff --git a/dlog/order23.c b/dlog/order23.c new file mode 100644 index 00000000..ae4d14d3 --- /dev/null +++ b/dlog/order23.c @@ -0,0 +1,32 @@ +/*============================================================================= + + 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" + +ulong +dlog_order23(const dlog_order23_t t, ulong b) +{ + return (b == *t) ? 1 : (b == 1) ? 0 : 2; +} diff --git a/dlog/order23_init.c b/dlog/order23_init.c new file mode 100644 index 00000000..412818de --- /dev/null +++ b/dlog/order23_init.c @@ -0,0 +1,39 @@ +/*============================================================================= + + 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" + +ulong +dlog_order23_init(dlog_order23_t t, ulong a) +{ + * t = a; + return 0; +} + +void +dlog_order23_clear(dlog_order23_t t) +{ + return; +} diff --git a/dlog/power_init.c b/dlog/power_init.c index 5ee59af8..61409db3 100644 --- a/dlog/power_init.c +++ b/dlog/power_init.c @@ -25,7 +25,7 @@ #include "dlog.h" -void +ulong dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num) { int k; @@ -40,7 +40,8 @@ dlog_power_init(dlog_power_t t, ulong a, ulong mod, ulong p, ulong e, ulong num) for (k = 1; k < e; k++) t->apk[k] = nmod_pow_ui(t->apk[k-1], p, t->mod); - dlog_precomp_p_init(t->pre, nmod_inv(t->apk[e-1], t->mod), mod, p, num); + dlog_precomp_p_init(t->pre, nmod_inv(t->apk[e-1], t->mod), mod, p, e * num); + return e * t->pre->cost; } void diff --git a/dlog/precomp.c b/dlog/precomp.c index 830b3975..8b03ff6b 100644 --- a/dlog/precomp.c +++ b/dlog/precomp.c @@ -47,6 +47,9 @@ dlog_precomp(const dlog_precomp_t pre, ulong b) case DLOG_BSGS: return dlog_bsgs(pre->t.bsgs, b); break; + case DLOG_23: + return dlog_order23(pre->t.order23, b); + break; default: abort(); break; diff --git a/dlog/precomp_init.c b/dlog/precomp_init.c index af7f531d..b2229294 100644 --- a/dlog/precomp_init.c +++ b/dlog/precomp_init.c @@ -26,30 +26,56 @@ #include "math.h" #include "dlog.h" +#define SMALL_LIM 50 +#define TABLE_LIM 50 #define TABLE_P_LIM 50 +#define TABLE_MODPE_LIM 50 #define TABLE_PE_LIM 50 #define TABLE_N_LIM 50 #define BSGS_LIM 500 +void +dlog_precomp_small_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num) +{ + if (n <= 3) + { + pre->type = DLOG_23; + pre->cost = dlog_order23_init(pre->t.order23, a); + } + else + { + if ( mod < TABLE_LIM ) + { + pre->type = DLOG_TABLE; + pre->cost = dlog_table_init(pre->t.table, a, mod); + } + else + { + pre->type = DLOG_BSGS; + pre->cost = dlog_bsgs_init(pre->t.bsgs, a, mod, n, n); + } + } +} /* log mod p^e */ void dlog_precomp_modpe_init(dlog_precomp_t pre, ulong a, ulong p, ulong e, ulong pe, ulong num) { - if ( pe < TABLE_PE_LIM ) + if ( pe < TABLE_MODPE_LIM ) { - pre->type = DLOG_TABLE; - dlog_table_init(pre->t.table, a, pe); + dlog_precomp_small_init(pre, a, pe, pe - pe / p, num); return; } - if (e > 1) - { - pre->type = DLOG_MODPE; - dlog_modpe_init(pre->t.modpe, a, p, e, pe, num); - } - else - { - dlog_precomp_n_init(pre, a, p, p - 1, num); + else { + if (e > 1) + { + pre->type = DLOG_MODPE; + pre->cost = dlog_modpe_init(pre->t.modpe, a, p, e, pe, num); + } + else + { + dlog_precomp_n_init(pre, a, p, p - 1, num); + } } } @@ -60,20 +86,21 @@ dlog_precomp_n_init(dlog_precomp_t pre, ulong a, ulong mod, ulong n, ulong num) if (n%2 && n_is_probabprime(n)) dlog_precomp_p_init(pre, a, mod, n, num); else { - if ( mod < TABLE_N_LIM ) + if (n < TABLE_N_LIM) + { + dlog_precomp_small_init(pre, a, mod, n, num); + } + else { - pre->type = DLOG_TABLE; - dlog_table_init(pre->t.table, a, mod); - } else { if (n < BSGS_LIM) { ulong m; m = (2 * num < n) ? ceil(sqrt((double) n * num)) : n; pre->type = DLOG_BSGS; - dlog_bsgs_init(pre->t.bsgs, a, mod, n, m); + pre->cost = dlog_bsgs_init(pre->t.bsgs, a, mod, n, m); } else { pre->type = DLOG_CRT; - dlog_crt_init(pre->t.crt, a, mod, n, num); + pre->cost = dlog_crt_init(pre->t.crt, a, mod, n, num); } } } @@ -85,15 +112,14 @@ dlog_precomp_p_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong num) { if ( p < TABLE_P_LIM ) { - pre->type = DLOG_TABLE; - dlog_table_init(pre->t.table, a, mod); + dlog_precomp_small_init(pre, a, mod, p, num); } else { ulong m; m = (2 * num < p) ? ceil(sqrt((double) p * num)) : p; pre->type = DLOG_BSGS; - dlog_bsgs_init(pre->t.bsgs, a, mod, p, m); + pre->cost = dlog_bsgs_init(pre->t.bsgs, a, mod, p, m); } } @@ -102,8 +128,7 @@ dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, u { if ( pe < TABLE_PE_LIM ) { - pre->type = DLOG_TABLE; - dlog_table_init(pre->t.table, a, mod); + dlog_precomp_small_init(pre, a, mod, pe, num); } else { @@ -114,7 +139,7 @@ dlog_precomp_pe_init(dlog_precomp_t pre, ulong a, ulong mod, ulong p, ulong e, u else { pre->type = DLOG_POWER; - dlog_power_init(pre->t.power, a, mod, p, e, num); + pre->cost = dlog_power_init(pre->t.power, a, mod, p, e, num); } } } @@ -139,6 +164,9 @@ dlog_precomp_clear(dlog_precomp_t pre) case DLOG_BSGS: dlog_bsgs_clear(pre->t.bsgs); break; + case DLOG_23: + dlog_order23_clear(pre->t.order23); + break; default: abort(); break; diff --git a/dlog/profile/p-precomp.c b/dlog/profile/p-precomp.c index 3b0638d9..c4aa9c1b 100644 --- a/dlog/profile/p-precomp.c +++ b/dlog/profile/p-precomp.c @@ -27,7 +27,10 @@ #include "profiler.h" #include -typedef void (*multilog_f) (ulong p, ulong a, ulong num); +#define NUMPRIMES 400 +#define CSV 0 + +typedef void (*log_f) (ulong p, ulong a, ulong num); void flog_table(ulong p, ulong a, ulong num) @@ -37,7 +40,7 @@ flog_table(ulong p, ulong a, ulong num) dlog_table_init(t, a, p); for (k = 1; k < num; k++) { - if (k == p) continue; + if (k % p == 0) continue; dlog_table(t, k % p); } dlog_table_clear(t); @@ -50,7 +53,7 @@ flog_bsgs(ulong p, ulong a, ulong num) dlog_bsgs_init(t, a, p, p-1, ceil( sqrt((double)num * p))); for (k = 1; k < num; k++) { - if (k == p) continue; + if (k % p == 0) continue; dlog_bsgs(t, k % p); } dlog_bsgs_clear(t); @@ -63,7 +66,7 @@ flog_rho(ulong p, ulong a, ulong num) dlog_rho_init(t, a, p, p-1); for (k = 1; k < num; k++) { - if (k == p) continue; + if (k % p == 0) continue; dlog_rho(t, k % p); } dlog_rho_clear(t); @@ -76,7 +79,7 @@ flog_crt(ulong p, ulong a, ulong num) dlog_crt_init(t, a, p, p-1, num); for (k = 1; k < num; k++) { - if (k == p) continue; + if (k % p == 0) continue; dlog_crt(t, k % p); } dlog_crt_clear(t); @@ -89,37 +92,38 @@ flog_gen(ulong p, ulong a, ulong num) dlog_precomp_n_init(t, a, p, p-1, num); for (k = 1; k < num; k++) { - if (k == p) continue; + if (k % p == 0) continue; dlog_precomp(t, k % p); } dlog_precomp_clear(t); } -void -dlog_bench(multilog_f flog, const ulong * p, const ulong * a, int np, int num) -{ - int i; - ulong q; - TIMEIT_ONCE_START - for (i = 0; i < np; i ++) - flog(p[i], a[i], num); - TIMEIT_ONCE_STOP -} int main() { slong iter, k, nv, nref, r, nr; ulong minq, maxq; ulong * rand; - int nbits, nl = 4; - int vl[5] = { 1, 10, 100, 1000 , 5000}; + int nbits, nl = 5; + int l[5] = { 1, 10, 100, 1000 , 5000}; + + int nf = 4; + log_f func[4] = { flog_table, flog_bsgs, flog_crt, flog_gen }; + char * n[4] = { "table", "bsgs", "crt", "generic" }; + + int np = NUMPRIMES; + flint_rand_t state; flint_randinit(state); - for (nbits = 10; nbits < 52; nbits += 10) + for (nbits = 10; nbits < 50; nbits += 5) { - int i, np = 100; - ulong p[100], a[100]; + int i; + ulong p[NUMPRIMES], a[NUMPRIMES]; + + if (nbits >= 25) + np /= 2; + for (i=0; i < np; i++) { p[i] = n_randprime(state, nbits, 0); @@ -128,36 +132,33 @@ int main() for (i = 0; i < nl; i++) { - flint_printf("%wu * logs mod primes of size %wu.....\n", vl[i], nbits); + int f; +#if LOG + flint_printf("%d * logs mod primes of size %d.....\n", l[i], nbits); +#endif - if (nbits <= 20) - { - flint_printf("table.......... "); - fflush(stdout); - dlog_bench(flog_table, p, a, np, vl[i]); - } - - if (nbits <= 40 || vl[i] <= 10) + for (f = 0; f < nf; f++) { - flint_printf("bsgs........... "); + int j; + + + /* skip useless */ + if (f == 0 && nbits >= 20) + continue; + if (f == 1 && nbits >= 30 && l[i] > 10) + continue; +#if LOG + flint_printf("%-20s... ",n[f]); fflush(stdout); - dlog_bench(flog_bsgs, p, a, np, vl[i]); +#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 (nbits <= 20 || vl[i] == 1) - { - flint_printf("rho............ "); - fflush(stdout); - dlog_bench(flog_rho, p, a, np, vl[i]); - } - - flint_printf("crt............ "); - fflush(stdout); - dlog_bench(flog_crt, p, a, np, vl[i]); - - flint_printf("generic........ "); - fflush(stdout); - dlog_bench(flog_gen, p, a, np, vl[i]); } } flint_randclear(state); diff --git a/dlog/table_init.c b/dlog/table_init.c index 7a19cf68..5ed9f4f1 100644 --- a/dlog/table_init.c +++ b/dlog/table_init.c @@ -26,7 +26,7 @@ #include "dlog.h" /* assume mod is small so no overflow */ -void +ulong dlog_table_init(dlog_table_t t, ulong a, ulong mod) { int k; @@ -41,6 +41,7 @@ dlog_table_init(dlog_table_t t, ulong a, ulong mod) ak = (ak * a) % mod; } while (ak != 1); + return 1; } void