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

fix log table for large moduli + add order 2 & 3

Browse source
This commit is contained in:
Pascal 2016-03-15 15:39:26 +01:00
parent 387c7f19e3
commit 32a1c966b9
12 changed files with 214 additions and 87 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
dlog.h
View file

@ -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);

6
dlog/bsgs.c
View file

@ -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();
}

7
dlog/bsgs_init.c
View file

@ -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

12
dlog/crt_init.c
View file

@ -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

3
dlog/modpe.c
View file

@ -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

32
dlog/order23.c Normal file
View file

@ -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;
}

39
dlog/order23_init.c Normal file
View file

@ -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;
}

5
dlog/power_init.c
View file

@ -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

3
dlog/precomp.c
View file

@ -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;

74
dlog/precomp_init.c
View file

@ -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;

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

@ -27,7 +27,10 @@
#include "profiler.h"
#include <math.h>
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);

3
dlog/table_init.c
View file

@ -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