1 /*
2 Copyright (C) 2010 Fredrik Johansson
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 #include <gmp.h>
13 #include "flint.h"
14 #include "fmpz.h"
15 #include "fmpz_vec.h"
16 #include "mpn_extras.h"
17 #include "ulong_extras.h"
18
19 static slong trial_cutoff[15] = {4, 4, 4, 6, 11, 18, 31, 54, 97, 172, 309, 564, 1028, 1900, 3512};
20
21 /*
22 Tuning values to give a roughly 1/3 chance of finding a factor of the given
23 number of bits. Parameters are {bits, B1, curves}. B2 should be taken to be
24 100*B1. Tuning values are a bit rough from 62 bits on.
25 */
26 static slong ecm_tuning[][3] =
27 {
28 {0, 0, 0}, {2, 1, 1}, {4, 3, 1}, {6, 5, 1}, {8, 7, 1},
29 {10, 9, 2}, {12, 11, 2}, {14, 13, 2}, {16, 15, 2}, {18, 17, 2},
30 {20, 19, 5}, {22, 25, 8}, {24, 32, 10}, {26, 47, 11}, {28, 67, 13},
31 {30, 102, 13}, {32, 126, 13}, {34, 207, 15}, {36, 293, 16}, {38, 415, 17},
32 {40, 610, 18}, {42, 920, 18}, {44, 1270, 20}, {46, 1800, 20}, {48, 2650, 20},
33 {50, 3850, 21}, {52, 5300, 22}, {54, 8500, 22}, {56, 10000, 26}, {58, 12000, 33},
34 {60, 14000, 42}, {62, 15000, 57}, {64, 16500, 72}, {66, 18000, 87}, {68, 22000, 102},
35 {70, 26000, 117}, {72, 30000, 131}, {74, 40000, 146}, {76, 50000, 161}, {78, 60000, 175},
36 {80, 70000, 190}, {82, 80000, 205}, {84, 100000, 220}, {86, 140000, 240}, {88, 190000, 255},
37 {90, 240000, 291}, {92, 280000, 318}, {94, 320000, 345}, {96, 370000, 372}, {98, 420000, 400},
38 {100, 470000, 430}
39 };
40
_is_prime(const fmpz_t n,int proved)41 int _is_prime(const fmpz_t n, int proved)
42 {
43 if (proved)
44 return fmpz_is_prime(n);
45 else
46 return fmpz_is_probabprime(n);
47 }
48
remove_found_factors(fmpz_factor_t factor,fmpz_t n,fmpz_t f)49 void remove_found_factors(fmpz_factor_t factor, fmpz_t n, fmpz_t f)
50 {
51 slong i;
52 fmpz_factor_t fac;
53
54 fmpz_tdiv_q(n, n, f);
55
56 fmpz_factor_init(fac);
57 fmpz_factor_no_trial(fac, f);
58
59 for (i = 0; i < fac->num; i++)
60 fac->exp[i] += fmpz_remove(n, n, fac->p + i);
61
62 _fmpz_factor_concat(factor, fac, 1);
63
64 fmpz_factor_clear(fac);
65 }
66
fmpz_factor_smooth(fmpz_factor_t factor,const fmpz_t n,slong bits,int proved)67 int fmpz_factor_smooth(fmpz_factor_t factor, const fmpz_t n,
68 slong bits, int proved)
69 {
70 ulong exp;
71 mp_limb_t p;
72 __mpz_struct * xsrc;
73 mp_ptr xd;
74 mp_size_t xsize;
75 slong found;
76 slong trial_stop;
77 slong * idx;
78 slong i, b, bits2, istride;
79 const mp_limb_t * primes;
80 int ret = 0;
81
82 TMP_INIT;
83
84 if (!COEFF_IS_MPZ(*n))
85 {
86 fmpz_factor_si(factor, *n);
87 return 1;
88 }
89
90 _fmpz_factor_set_length(factor, 0);
91
92 /* Get sign and size */
93 xsrc = COEFF_TO_PTR(*n);
94 if (xsrc->_mp_size < 0)
95 {
96 xsize = -(xsrc->_mp_size);
97 factor->sign = -1;
98 }
99 else
100 {
101 xsize = xsrc->_mp_size;
102 factor->sign = 1;
103 }
104
105 /* Just a single limb */
106 if (xsize == 1)
107 {
108 _fmpz_factor_extend_factor_ui(factor, xsrc->_mp_d[0]);
109 return 1;
110 }
111
112 /* Create a temporary copy to be mutated */
113 TMP_START;
114 xd = TMP_ALLOC(xsize * sizeof(mp_limb_t));
115 flint_mpn_copyi(xd, xsrc->_mp_d, xsize);
116
117 /* Factor out powers of two */
118 xsize = flint_mpn_remove_2exp(xd, xsize, &exp);
119 if (exp != 0)
120 _fmpz_factor_append_ui(factor, UWORD(2), exp);
121
122 if (bits <= 0)
123 {
124 flint_printf("(fmpz_factor_smooth) Number of bits must be at least 1\n");
125 flint_abort();
126 }
127
128 if (bits <= 15)
129 trial_stop = trial_cutoff[bits - 1];
130 else
131 trial_stop = 3512;
132
133 b = fmpz_sizeinbase(n, 2) - exp;
134 idx = (slong *) flint_malloc((5 + b/4)*sizeof(slong));
135
136 found = flint_mpn_factor_trial_tree(idx, xd, xsize, trial_stop);
137
138 if (found)
139 {
140 primes = n_primes_arr_readonly(trial_stop);
141
142 for (i = 0; i < found; i++)
143 {
144 p = primes[idx[i]];
145
146 exp = 1;
147 xsize = flint_mpn_divexact_1(xd, xsize, p);
148
149 /* Check if p^2 divides n */
150 if (flint_mpn_divisible_1_p(xd, xsize, p))
151 {
152 xsize = flint_mpn_divexact_1(xd, xsize, p);
153 exp = 2;
154 }
155
156 /* If we're up to cubes, then maybe there are higher powers */
157 if (exp == 2 && flint_mpn_divisible_1_p(xd, xsize, p))
158 {
159 xsize = flint_mpn_divexact_1(xd, xsize, p);
160 xsize = flint_mpn_remove_power_ascending(xd, xsize, &p, 1, &exp);
161 exp += 3;
162 }
163
164 _fmpz_factor_append_ui(factor, p, exp);
165 }
166 }
167
168 if (xsize == 1)
169 {
170 /* Any single-limb factor left? */
171 if (xd[0] != 1)
172 _fmpz_factor_extend_factor_ui(factor, xd[0]);
173
174 ret = 1;
175 }
176 else
177 {
178 fmpz_t n2, f;
179 __mpz_struct * data;
180
181 fmpz_init2(n2, xsize);
182
183 data = _fmpz_promote(n2);
184 flint_mpn_copyi(data->_mp_d, xd, xsize);
185 data->_mp_size = xsize;
186
187 if (proved != -1 && _is_prime(n2, proved))
188 {
189 _fmpz_factor_append(factor, n2, 1);
190 ret = 1;
191 }
192 else
193 {
194 fmpz_t root;
195
196 fmpz_init(root);
197
198 exp = fmpz_is_perfect_power(root, n2);
199
200 if (exp != 0)
201 {
202 fmpz_factor_t fac;
203
204 fmpz_factor_init(fac);
205
206 ret = fmpz_factor_smooth(fac, root, bits, proved);
207 fmpz_set_ui(n2, 1);
208
209 _fmpz_factor_concat(factor, fac, exp);
210
211 fmpz_factor_clear(fac);
212 }
213 else if (bits >= 16) /* trial factored already up to 15 bits */
214 {
215 int found;
216 flint_rand_t state;
217
218 fmpz_init(f);
219 flint_randinit(state);
220
221 /* currently only tuning values up to factors of 100 bits */
222 bits = FLINT_MIN(bits, 100);
223 bits2 = (bits + 1)/2;
224
225 /* tuning is in increments of 2 bits */
226 istride = 3;
227 /* start with 18-22 bits, advance by 6 bits at a time */
228 for (i = 9 + (bits2 % 3); i <= bits2; i += istride)
229 {
230 found = fmpz_factor_ecm(f, ecm_tuning[i][2],
231 ecm_tuning[i][1], ecm_tuning[i][1]*100, state, n2);
232
233 if (found != 0)
234 {
235 /* make sure all prime divisors in factor are removed from n2 */
236 remove_found_factors(factor, n2, f);
237
238 if (fmpz_is_one(n2))
239 {
240 ret = 1;
241
242 break;
243 }
244
245 /* if what remains is below the bound, just factor it */
246 if (fmpz_sizeinbase(n2, 2) < bits)
247 {
248 fmpz_factor_no_trial(factor, n2);
249
250 ret = 1;
251
252 break;
253 }
254
255 if (_is_prime(n2, proved))
256 {
257 _fmpz_factor_append(factor, n2, 1);
258
259 ret = 1;
260
261 break;
262 }
263
264 i -= istride; /* redo with the same parameters if factor found */
265 }
266 }
267
268 flint_randclear(state);
269 fmpz_clear(f);
270 }
271 }
272
273 if (ret != 1 && !fmpz_is_one(n2))
274 _fmpz_factor_append(factor, n2, 1); /* place cofactor in factor struct */
275 else
276 ret = 1;
277
278 fmpz_clear(n2);
279 }
280
281 flint_free(idx);
282
283 TMP_END;
284 return ret;
285 }
286
287