1 /*
2 Copyright (C) 2015 Kushagra Singh
3
4 This file is part of FLINT.
5
6 FLINT 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 <https://www.gnu.org/licenses/>.
10 */
11
12 /* Outer wrapper for ECM
13 makes calls to stage I and stage II (one) */
14
15 #include <gmp.h>
16 #include "flint.h"
17 #include "fmpz.h"
18 #include "mpn_extras.h"
19
20 static
21 ulong n_ecm_primorial[] =
22 {
23 #ifdef FLINT64
24
25 UWORD(2), UWORD(6), UWORD(30), UWORD(210), UWORD(2310), UWORD(30030),
26 UWORD(510510), UWORD(9699690), UWORD(223092870), UWORD(6469693230),
27 UWORD(200560490130), UWORD(7420738134810), UWORD(304250263527210),
28 UWORD(13082761331670030), UWORD(614889782588491410)
29 /* 15 values */
30
31 #else
32
33 UWORD(2), UWORD(6), UWORD(30), UWORD(210), UWORD(2310), UWORD(30030),
34 UWORD(510510), UWORD(9699690)
35 /* 9 values */
36
37 #endif
38 };
39
40 #ifdef FLINT64
41 #define num_n_ecm_primorials 15
42 #else
43 #define num_n_ecm_primorials 9
44 #endif
45
46 int
fmpz_factor_ecm(fmpz_t f,mp_limb_t curves,mp_limb_t B1,mp_limb_t B2,flint_rand_t state,const fmpz_t n_in)47 fmpz_factor_ecm(fmpz_t f, mp_limb_t curves, mp_limb_t B1, mp_limb_t B2,
48 flint_rand_t state, const fmpz_t n_in)
49 {
50 fmpz_t sig, nm8;
51 mp_limb_t P, num, maxP, mmin, mmax, mdiff, prod, maxj, n_size, cy;
52 int i, j, ret;
53 ecm_t ecm_inf;
54 __mpz_struct *fac, *mpz_ptr;
55 mp_ptr n, mpsig;
56
57 TMP_INIT;
58
59 const mp_limb_t *prime_array;
60 n_size = fmpz_size(n_in);
61
62 if (n_size == 1)
63 {
64 ret = n_factor_ecm(&P, curves, B1, B2, state, fmpz_get_ui(n_in));
65 fmpz_set_ui(f, P);
66 return ret;
67 }
68
69 fmpz_factor_ecm_init(ecm_inf, n_size);
70
71 TMP_START;
72
73 n = TMP_ALLOC(n_size * sizeof(mp_limb_t));
74 mpsig = TMP_ALLOC(n_size * sizeof(mp_limb_t));
75
76 if ((!COEFF_IS_MPZ(* n_in)))
77 {
78 count_leading_zeros(ecm_inf->normbits, fmpz_get_ui(n_in));
79 n[0] = fmpz_get_ui(n_in);
80 n[0] <<= ecm_inf->normbits;
81 }
82 else
83 {
84 mpz_ptr = COEFF_TO_PTR(* n_in);
85 count_leading_zeros(ecm_inf->normbits, mpz_ptr->_mp_d[n_size - 1]);
86 if (ecm_inf->normbits)
87 mpn_lshift(n, mpz_ptr->_mp_d, n_size, ecm_inf->normbits);
88 else
89 flint_mpn_copyi(n, mpz_ptr->_mp_d, n_size);
90 }
91
92 flint_mpn_preinvn(ecm_inf->ninv, n, n_size);
93 ecm_inf->one[0] = UWORD(1) << ecm_inf->normbits;
94
95 fmpz_init(sig);
96 fmpz_init(nm8);
97 fmpz_sub_ui(nm8, n_in, 8);
98
99 ret = 0;
100 fac = _fmpz_promote(f);
101 mpz_realloc(fac, fmpz_size(n_in));
102
103 /************************ STAGE I PRECOMPUTATIONS ************************/
104
105 num = n_prime_pi(B1); /* number of primes under B1 */
106
107 /* compute list of primes under B1 for stage I */
108 prime_array = n_primes_arr_readonly(num);
109
110 /************************ STAGE II PRECOMPUTATIONS ***********************/
111
112 maxP = n_sqrt(B2);
113
114 /* Selecting primorial */
115
116 j = 1;
117 while ((j < num_n_ecm_primorials) && (n_ecm_primorial[j] < maxP))
118 j += 1;
119
120 P = n_ecm_primorial[j - 1];
121
122 mmin = (B1 + (P/2)) / P;
123 mmax = ((B2 - P/2) + P - 1)/P; /* ceil */
124 if (mmax < mmin)
125 {
126 flint_printf("Exception (ecm). B1 > B2 encountered.\n");
127 flint_abort();
128 }
129 maxj = (P + 1)/2;
130 mdiff = mmax - mmin + 1;
131
132 /* compute GCD_table */
133
134 ecm_inf->GCD_table = flint_malloc(maxj + 1);
135
136 for (j = 1; j <= maxj; j += 2)
137 {
138 if ((j%2) && n_gcd(j, P) == 1)
139 ecm_inf->GCD_table[j] = 1;
140 else
141 ecm_inf->GCD_table[j] = 0;
142 }
143
144 /* compute prime table */
145
146 ecm_inf->prime_table = flint_malloc(mdiff * sizeof(unsigned char*));
147
148 for (i = 0; i < mdiff; i++)
149 ecm_inf->prime_table[i] = flint_malloc((maxj + 1) * sizeof(unsigned char));
150
151 for (i = 0; i < mdiff; i++)
152 {
153 for (j = 1; j <= maxj; j += 2)
154 {
155 ecm_inf->prime_table[i][j] = 0;
156
157 /* if (i + mmin)*P + j
158 is prime, mark 1. Can be possibly prime
159 only if gcd(j, P) = 1 */
160
161 if (ecm_inf->GCD_table[j] == 1)
162 {
163 prod = (i + mmin)*P + j;
164 if (n_is_prime(prod))
165 ecm_inf->prime_table[i][j] = 1;
166
167 prod = (i + mmin)*P - j;
168 if (n_is_prime(prod))
169 ecm_inf->prime_table[i][j] = 1;
170 }
171 }
172 }
173
174
175 /****************************** TRY "CURVES" *****************************/
176
177 for (j = 0; j < curves; j++)
178 {
179 fmpz_randm(sig, state, nm8);
180 fmpz_add_ui(sig, sig, 7);
181
182 mpn_zero(mpsig, ecm_inf->n_size);
183
184 if ((!COEFF_IS_MPZ(*sig)))
185 {
186 mpsig[0] = fmpz_get_ui(sig);
187 if (ecm_inf->normbits)
188 {
189 cy = mpn_lshift(mpsig, mpsig, 1, ecm_inf->normbits);
190 if (cy)
191 mpsig[1] = cy;
192 }
193 }
194 else
195 {
196 mpz_ptr = COEFF_TO_PTR(*sig);
197
198 if (ecm_inf->normbits)
199 {
200 cy = mpn_lshift(mpsig, mpz_ptr->_mp_d, mpz_ptr->_mp_size, ecm_inf->normbits);
201 if (cy)
202 mpsig[mpz_ptr->_mp_size] = cy;
203 } else
204 {
205 flint_mpn_copyi(mpsig, mpz_ptr->_mp_d, mpz_ptr->_mp_size);
206 }
207 }
208
209 /************************ SELECT CURVE ************************/
210
211 ret = fmpz_factor_ecm_select_curve(fac->_mp_d, mpsig, n, ecm_inf);
212
213 if (ret)
214 {
215 /* Found factor while selecting curve,
216 very very lucky :) */
217
218 if (ecm_inf->normbits)
219 mpn_rshift(fac->_mp_d, fac->_mp_d, ret, ecm_inf->normbits);
220 MPN_NORM(fac->_mp_d, ret);
221
222 fac->_mp_size = ret;
223 _fmpz_demote_val(f);
224 ret = -1;
225 goto cleanup;
226 }
227
228 if (ret != -1)
229 {
230
231 /************************** STAGE I ***************************/
232
233 ret = fmpz_factor_ecm_stage_I(fac->_mp_d, prime_array, num, B1, n, ecm_inf);
234
235 if (ret)
236 {
237 /* Found factor after stage I */
238 if (ecm_inf->normbits)
239 mpn_rshift(fac->_mp_d, fac->_mp_d, ret, ecm_inf->normbits);
240 MPN_NORM(fac->_mp_d, ret);
241
242 fac->_mp_size = ret;
243 _fmpz_demote_val(f);
244 ret = 1;
245 goto cleanup;
246 }
247 /************************** STAGE II ***************************/
248
249 ret = fmpz_factor_ecm_stage_II(fac->_mp_d, B1, B2, P, n, ecm_inf);
250
251 if (ret)
252 {
253 /* Found factor after stage II */
254 if (ecm_inf->normbits)
255 mpn_rshift(fac->_mp_d, fac->_mp_d, ret, ecm_inf->normbits);
256 MPN_NORM(fac->_mp_d, ret);
257 fac->_mp_size = ret;
258 _fmpz_demote_val(f);
259 ret = 2;
260 goto cleanup;
261 }
262 }
263 }
264
265 cleanup:
266
267
268 flint_free(ecm_inf->GCD_table);
269 for (i = 0; i < mdiff; i++)
270 flint_free(ecm_inf->prime_table[i]);
271 flint_free(ecm_inf->prime_table);
272
273 fmpz_factor_ecm_clear(ecm_inf);
274
275 fmpz_clear(nm8);
276 fmpz_clear(sig);
277
278 TMP_END;
279
280 return ret;
281 }
282