1 /*
2 Copyright (C) 2016 Pascal Molin
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "dlog.h"
13 #include <math.h>
14
15 static ulong
dlog_single(ulong b,ulong a,const nmod_t mod,ulong n)16 dlog_single(ulong b, ulong a, const nmod_t mod, ulong n)
17 {
18 if (n < 50)
19 {
20 int k;
21 ulong ak = 1;
22
23 for (k=0; k < n; k++)
24 {
25 if (ak == b)
26 return k;
27 ak = nmod_mul(ak, a, mod);
28 }
29
30 flint_printf("FAIL[dlog single]: log(%wu,%wu) mod %wu not found (size %wu)\n",
31 b, a, mod.n, n);
32 flint_abort();
33 return 0; /* dummy return because flint_abort() is not declared noreturn */
34 }
35 else
36 {
37 dlog_rho_t t;
38 dlog_rho_init(t, a, mod.n, n);
39 return dlog_rho(t, b);
40 }
41 }
42
43 /* solve log knowing equation e = f * log(b) [n] */
44 static ulong
dlog_quotient(const dlog_rho_t t,ulong e,ulong f,ulong g,ulong b)45 dlog_quotient(const dlog_rho_t t, ulong e, ulong f, ulong g, ulong b)
46 {
47 ulong r, b_ar, an;
48 nmod_t n = t->n;
49
50 if (g == n.n)
51 {
52 flint_printf("FAIL[dlog quotient]: trivial relation e = %wu, f = %wu mod %wu\n",
53 e, f, n.n);
54 flint_abort();
55 }
56
57 nmod_init(&n, n.n / g);
58 e = e / g;
59 f = f / g;
60 r = nmod_div(e, f, n);
61 an = nmod_pow_ui(t->a, n.n, t->mod);
62 b_ar = nmod_div(b, nmod_pow_ui(t->a, r, t->mod), t->mod);
63
64 return r + n.n * dlog_single(b_ar, an, t->mod, g);
65 }
66
67 #define RWALK 20
68 ulong
dlog_rho(const dlog_rho_t t,ulong b)69 dlog_rho(const dlog_rho_t t, ulong b)
70 {
71 int j, k, l;
72 ulong m[RWALK], n[RWALK], ab[RWALK];
73 ulong x[2], e[2], f[2], g;
74 flint_rand_t state;
75
76 flint_randinit(state);
77
78 do {
79
80 for (k = 0; k < RWALK; k++)
81 {
82 m[k] = 1 + n_randint(state, t->n.n - 1);
83 n[k] = 1 + n_randint(state, t->n.n - 1);
84 ab[k] = nmod_mul(nmod_pow_ui(t->a, m[k], t->mod), nmod_pow_ui(b, n[k], t->mod), t->mod);
85 }
86
87 /* x[l] = a^e[l] * b^f[l] */
88 x[0] = x[1] = 1;
89 e[0] = e[1] = 0;
90 f[0] = f[1] = 0;
91
92 do {
93
94 for(j = 0; j < 3; j++)
95 {
96 l = (j > 0);
97 k = floor( (double) RWALK * x[l] / t->mod.n );
98 x[l] = nmod_mul(x[l], ab[k], t->mod);
99 e[l] = nmod_add(e[l], m[k], t->n);
100 f[l] = nmod_add(f[l], n[k], t->n);
101 }
102
103 } while (x[0] != x[1]);
104
105 } while (e[0] == e[1] && f[0] == f[1]);
106
107 flint_randclear(state);
108
109 /* e = f * log(b) */
110 e[0] = nmod_sub(e[0], e[1], t->n);
111 f[0] = nmod_sub(f[1], f[0], t->n);
112
113 if (!t->nisprime && (g = n_gcd(f[0], t->n.n)) > 1)
114 return dlog_quotient(t, e[0], f[0], g, b);
115 else
116 return nmod_div(e[0], f[0], t->n);
117 }
118