1 /*=============================================================================
2
3 This file is part of Antic.
4
5 Antic is free software: you can redistribute it and/or modify it under
6 the terms of the GNU Lesser General Public License (LGPL) as published
7 by the Free Software Foundation; either version 2.1 of the License, or
8 (at your option) any later version. See <http://www.gnu.org/licenses/>.
9
10 =============================================================================*/
11 /******************************************************************************
12
13 Copyright (C) 2012 William Hart
14
15 ******************************************************************************/
16
17 #include <stdlib.h>
18 #include <gmp.h>
19 #include "qfb.h"
20
21 /*
22 find which power of the base is the exponent of f
23 */
find_power(qfb_t f,fmpz_t n,ulong base)24 ulong find_power(qfb_t f, fmpz_t n, ulong base)
25 {
26 ulong s = 1;
27
28 do
29 {
30 qfb_pow_ui(f, f, n, base);
31 s *= base;
32 } while (!qfb_is_principal_form(f, n));
33
34 return s;
35 }
36
qfb_exponent_element_stage2(qfb_t f,fmpz_t n,ulong B2_sqrt)37 ulong qfb_exponent_element_stage2(qfb_t f, fmpz_t n, ulong B2_sqrt)
38 {
39 qfb_t pow, pow2, f2;
40 fmpz_t L, r;
41 slong i, i2, ret = 0;
42 slong depth = FLINT_BIT_COUNT(B2_sqrt) + 1;
43 qfb_hash_t * qhash = qfb_hash_init(depth);
44
45 fmpz_init(L);
46 fmpz_init(r);
47 fmpz_abs(L, n);
48 fmpz_root(L, L, 4);
49
50 qfb_init(f2);
51 qfb_init(pow);
52 qfb_init(pow2);
53 qfb_hash_insert(qhash, f, NULL, 1, depth);
54
55 qfb_nucomp(f2, f, f, n, L); /* large primes are odd */
56 qfb_reduce(f2, f2, n);
57
58 qfb_set(pow, f);
59
60 for (i = 1; i < B2_sqrt - 1; i += 2) /* baby steps */
61 {
62 qfb_nucomp(pow, pow, f2, n, L);
63 qfb_reduce(pow, pow, n);
64
65 qfb_hash_insert(qhash, pow, NULL, i + 2, depth);
66 }
67
68 qfb_nucomp(pow, pow, f, n, L); /* compute f^B2_sqrt */
69 qfb_reduce(pow, pow, n);
70
71 qfb_nucomp(pow, pow, pow, n, L); /* we hash for a form or its inverse,
72 so we can jump by f^(2x B2_sqrt) */
73 qfb_reduce(pow, pow, n);
74 qfb_set(pow2, pow);
75
76 for(i = 2; i <= B2_sqrt; i += 2) /* giant steps */
77 {
78 i2 = qfb_hash_find(qhash, pow2, depth);
79 if (i2 != -1) /* found collision */
80 {
81 fmpz_set_ui(r, B2_sqrt);
82 fmpz_mul_ui(r, r, i);
83 if (fmpz_sgn(qhash[i2].q->b) == fmpz_sgn(pow2->b))
84 fmpz_sub_ui(r, r, qhash[i2].iter);
85 else
86 fmpz_add_ui(r, r, qhash[i2].iter);
87 ret = (fmpz_size(r) > 1 ? 0 : fmpz_get_ui(r)); /* we probably should be more aggressive here */
88 break;
89 }
90
91 qfb_nucomp(pow2, pow2, pow, n, L);
92 qfb_reduce(pow2, pow2, n);
93 }
94
95 fmpz_clear(r);
96 fmpz_clear(L);
97 qfb_clear(f2);
98 qfb_clear(pow);
99 qfb_clear(pow2);
100 qfb_hash_clear(qhash, depth);
101
102 return ret;
103 }
104
105 typedef struct
106 {
107 ulong pr;
108 qfb_t pow;
109 slong i;
110 } qfb_restart_t;
111
112 #define go_restart \
113 do { \
114 need_restarts = 0; \
115 if (j == 0) \
116 { \
117 qfb_pow(f2, f2, n, exponent); \
118 goto do_restart1; \
119 } \
120 j--; \
121 qfb_set(pow, restart[j].pow); \
122 qfb_pow(pow, pow, n, exponent); \
123 i = restart[j].i; \
124 n_primes_jump_after(iter, restart[j].pr); \
125 pr = restart[j].pr; \
126 goto do_restart; \
127 } while (0)
128
qfb_exponent_element(fmpz_t exponent,qfb_t f,fmpz_t n,ulong B1,ulong B2_sqrt)129 int qfb_exponent_element(fmpz_t exponent, qfb_t f, fmpz_t n, ulong B1, ulong B2_sqrt)
130 {
131 slong i, j, iters = 1024, restart_inc;
132 qfb_t pow, oldpow, f2;
133 ulong pr, oldpr, s2, sqrt, exp;
134 fmpz_t prod, L, pow2;
135 int ret = 1, num_restarts = 0, need_restarts = 1, clean2 = 0;
136 qfb_restart_t restart[128];
137 n_primes_t iter;
138 ulong hi, lo;
139 double quot;
140 mp_bitcnt_t bits0;
141
142 n_primes_init(iter);
143
144 sqrt = n_sqrt(B1);
145 bits0 = FLINT_BIT_COUNT(B1);
146
147 fmpz_set_ui(exponent, 1);
148
149 qfb_init(f2);
150 qfb_init(pow);
151 qfb_init(oldpow);
152
153 fmpz_init(pow2);
154 fmpz_init(prod);
155 fmpz_init(L);
156
157 fmpz_abs(L, n);
158 fmpz_root(L, L, 4);
159
160 qfb_set(f2, f);
161
162 do_restart1:
163
164 if (qfb_is_principal_form(f2, n))
165 goto cleanup2;
166
167 /* raise to appropriate power of 2 */
168 qfb_set(oldpow, f2);
169 fmpz_set_ui(pow2, 2);
170 fmpz_pow_ui(pow2, pow2, bits0);
171 qfb_pow(pow, oldpow, n, pow2);
172 if (qfb_is_principal_form(pow, n))
173 {
174 ulong s = find_power(oldpow, n, 2);
175 fmpz_mul_ui(exponent, exponent, s);
176
177 goto cleanup2;
178 }
179
180 for (i = 0; i < 128; i++)
181 qfb_init(restart[i].pow);
182
183 clean2 = 1;
184
185 n_prime_pi_bounds(&lo, &hi, B1);
186 restart_inc = (((hi/1024 + 127)/128))*1024;
187 quot = (double) B2_sqrt/(double) lo;
188
189 pr = 2;
190 n_primes_jump_after(iter, 2);
191 i = 0;
192
193 do_restart:
194
195 /* keep raising iters by a factor of 2 until pr exceeds B1 or exponent found */
196 do
197 {
198 /* raise to prime powers until the identity is found */
199 for ( ; i < iters; )
200 {
201 if ((i % restart_inc) == 0 && need_restarts)
202 {
203 qfb_set(restart[num_restarts].pow, pow);
204 restart[num_restarts].i = i;
205 restart[num_restarts].pr = pr;
206 num_restarts++;
207 }
208
209 j = FLINT_MIN(i + 1024, iters);
210 qfb_set(oldpow, pow);
211 oldpr = pr;
212
213 for ( ; i < j; i+=4)
214 {
215 pr = n_primes_next(iter);
216 fmpz_set_ui(prod, pr);
217 if (pr < sqrt)
218 {
219 exp = bits0/FLINT_BIT_COUNT(pr);
220 fmpz_pow_ui(prod, prod, exp);
221 }
222 pr = n_primes_next(iter);
223 if (pr < sqrt)
224 {
225 exp = bits0/FLINT_BIT_COUNT(pr);
226 fmpz_set_ui(pow2, pr);
227 fmpz_pow_ui(pow2, pow2, exp);
228 fmpz_mul(prod, prod, pow2);
229 } else
230 fmpz_mul_ui(prod, prod, pr);
231 pr = n_primes_next(iter);
232 if (pr < sqrt)
233 {
234 exp = bits0/FLINT_BIT_COUNT(pr);
235 fmpz_set_ui(pow2, pr);
236 fmpz_pow_ui(pow2, pow2, exp);
237 fmpz_mul(prod, prod, pow2);
238 } else
239 fmpz_mul_ui(prod, prod, pr);
240 pr = n_primes_next(iter);
241 if (pr < sqrt)
242 {
243 exp = bits0/FLINT_BIT_COUNT(pr);
244 fmpz_set_ui(pow2, pr);
245 fmpz_pow_ui(pow2, pow2, exp);
246 fmpz_mul(prod, prod, pow2);
247 } else
248 fmpz_mul_ui(prod, prod, pr);
249 qfb_pow_with_root(pow, pow, n, prod, L);
250 }
251
252 /* identity is found, compute exponent recursively */
253 if (qfb_is_principal_form(pow, n))
254 {
255 qfb_set(pow, oldpow);
256 n_primes_jump_after(iter, oldpr);
257 pr = oldpr;
258 while (1)
259 {
260 pr = n_primes_next(iter);
261
262 if (pr < sqrt)
263 {
264 ulong k;
265 exp = bits0/FLINT_BIT_COUNT(pr);
266
267 for (k = 1; k <= exp; k++)
268 {
269 qfb_pow_ui(pow, pow, n, pr);
270 if (qfb_is_principal_form(pow, n))
271 {
272 fmpz_set_ui(pow2, pr);
273 fmpz_pow_ui(pow2, pow2, k);
274 fmpz_mul(exponent, exponent, pow2);
275 for (j = 0; j < num_restarts; j++)
276 {
277 qfb_set(pow, restart[j].pow);
278 qfb_pow(pow, pow, n, exponent);
279 if (qfb_is_principal_form(pow, n))
280 break;
281 }
282 go_restart;
283 }
284 }
285 } else
286 {
287 qfb_pow_ui(pow, pow, n, pr);
288 if (qfb_is_principal_form(pow, n))
289 {
290 fmpz_mul_ui(exponent, exponent, pr);
291 for (j = 0; j < num_restarts; j++)
292 {
293 qfb_set(pow, restart[j].pow);
294 qfb_pow(pow, pow, n, exponent);
295 if (qfb_is_principal_form(pow, n))
296 break;
297 }
298 go_restart;
299 }
300 }
301 }
302 }
303 }
304
305 /* stage 2 */
306 s2 = qfb_exponent_element_stage2(pow, n, (ulong) ((double) iters * quot));
307 if (s2 && n_is_prime(s2)) /* we probably should be more aggressive here */
308 {
309 fmpz_mul_ui(exponent, exponent, s2);
310 for (j = 0; j < num_restarts; j++)
311 {
312 qfb_set(pow, restart[j].pow);
313 qfb_pow(pow, pow, n, exponent);
314 if (qfb_is_principal_form(pow, n))
315 break;
316 }
317 go_restart;
318 }
319
320 iters = FLINT_MIN(2*iters, hi);
321 } while (pr <= B1);
322
323 ret = 0;
324
325 cleanup2:
326
327 if (clean2)
328 for (i = 0; i < 128; i++)
329 qfb_clear(restart[i].pow);
330
331 qfb_clear(f2);
332 qfb_clear(pow);
333 qfb_clear(oldpow);
334 fmpz_clear(pow2);
335 fmpz_clear(prod);
336 fmpz_clear(L);
337 n_primes_clear(iter);
338
339 return ret;
340 }
341