Line data Source code
1 : /* Choice of best parameters for stage 2.
2 :
3 : Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2010 Paul Zimmermann,
4 : Alexander Kruppa, Dave Newman.
5 :
6 : This file is part of the ECM Library.
7 :
8 : The ECM Library is free software; you can redistribute it and/or modify
9 : it under the terms of the GNU Lesser General Public License as published by
10 : the Free Software Foundation; either version 3 of the License, or (at your
11 : option) any later version.
12 :
13 : The ECM Library is distributed in the hope that it will be useful, but
14 : WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 : or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 : License for more details.
17 :
18 : You should have received a copy of the GNU Lesser General Public License
19 : along with the ECM Library; see the file COPYING.LIB. If not, see
20 : http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 : 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22 :
23 : #include <stdio.h>
24 : #include <gmp.h>
25 : #include "ecm-impl.h"
26 :
27 : /*
28 : Compute (d, d2, k) such that:
29 : (0) k >= k0
30 : (1) d is a multiple of 6
31 : (2) k * d * (eulerphi(d)/2) * d2 / eulerphi(d2) >= B2 - B2min
32 : (3) gcd(d, d2) == 1
33 : (4) k is minimal, subject to previous conditions
34 : (5) if parameter po2 is != 0, rounds dF up to a power of 2
35 : Return non-zero iff an error occurred (too large step 2 bound).
36 : */
37 :
38 : /* How we test whether given d,d2,dF,k,i0 parameters cover the desired
39 : B2min-B2 range:
40 :
41 : In stage 2 we generate all values p = f(i * d) +- f(j * d2) with
42 :
43 : 1. gcd (i, d2) == 1,
44 : 2. gcd (j, d) == 1,
45 : 3. j == 1 (mod 6),
46 : 4. 6|d
47 : 5. 1 <= j <= d - 5, (it's -5, not just -1, because of 3. and 4.)
48 : 6. i0 <= i <= i1
49 : 7. gcd (d, d2) == 1
50 :
51 : where f(x) is x^S or the S-th Dickson polynomial g_{S,-1}(x). Extra
52 : factors included by S>1 are not considered in this analysis, we assume
53 : S=1, f(x)=x so that p = i * d +- j * d2.
54 :
55 : (Note: i values greater than stated in 3. may be generated if we have
56 : to round up dF, for example to a power of 2. However, the root generation
57 : code can put anything it likes in those extra roots, so we make no
58 : assumption here that this will extend the range of the i values.)
59 :
60 : Hence the values at the high end of the stage 2 range that are not
61 : generated are
62 :
63 : p = (i1 + n) * d +- j * d2, n > 0
64 :
65 : and the smallest one of those is
66 :
67 : p = (i1 + 1) * d - (d - 5) * d2
68 : = d * (i1 - d2 + 1) + 5 * d2
69 :
70 : At the low end of stage 2, values not generated are
71 :
72 : p = (i0 - n) * d +- j * d2, n > 0
73 :
74 : the largest one being
75 :
76 : p = (i0 - 1) * d + (d - 5) * d2
77 : = d * (i0 + d2 - 1) - 5*d2
78 :
79 : Thus, values p that are coprime do d*d2 and
80 : d * (i0 + d2 - 1) - 5*d2 + 1 <= p <= d * (i1 - d2 + 1) + 5 * d2 - 1
81 : are included in stage 2.
82 :
83 :
84 : The number of roots of G we compute is k * dF. For d2 == 1, this means
85 : i1 = i0 + k * dF - 1 (-1 because both i0 and i1 are included).
86 :
87 : For d2 > 1, values j not coprime to d2 are skipped (see condition 1).
88 : The number of values in [1, i0] that are not coprime to d2 (with d2 prime)
89 : is floor (i0 / d2); in [1, i1] it is floor (i1 / d2).
90 : So we require that
91 : k * dF >= i1 - i0 + 1 - (floor (i1 / d2) - floor (i0 / d2))
92 :
93 :
94 : */
95 :
96 : int
97 1566 : bestD (root_params_t *root_params, unsigned long *finalk,
98 : unsigned long *finaldF, mpz_t B2min, mpz_t B2, int po2, int use_ntt,
99 : double maxmem, int treefile, mpmod_t modulus)
100 : {
101 : /* the following list contains successive values of b with
102 : increasing values of eulerphi(b). It was generated by the following
103 : Maple program:
104 : l := [[1,1]]:
105 : for b from 12 by 6 do
106 : d:=numtheory[phi](b)/2;
107 : while d <= l[nops(l)][2] do l:=subsop(nops(l)=NULL, l) od;
108 : n := nops(l);
109 : if b>1.1*l[n][1] then l := [op(l), [b,d]]; lprint(l) fi;
110 : od:
111 : */
112 : #define N 109
113 : static unsigned int l[N] = {12, 18, 30, 42, 60, 90, 120, 150, 210, 240, 270, 330, 420, 510, 630, 840, 1050, 1260, 1470, 1680, 1890, 2310, 2730, 3150, 3570, 3990, 4620, 5460, 6090, 6930, 8190, 9240, 10920, 12180, 13860, 16170, 18480, 20790, 23100, 30030, 34650, 39270, 43890, 48510, 60060, 66990, 78540, 90090, 99330, 120120, 133980, 150150, 180180, 210210, 240240, 270270, 300300, 334950, 371280, 420420, 510510, 570570, 600600, 630630, 746130, 870870, 1021020, 1141140, 1291290, 1531530, 1711710, 1891890, 2081310, 2312310, 2552550, 2852850, 3183180, 3573570, 3993990, 4594590, 5105100, 5705700, 6322470, 7147140, 7987980, 8978970, 10210200, 11741730, 13123110, 14804790, 16546530, 19399380, 21411390, 23993970, 26816790, 29609580, 33093060, 36606570, 40330290, 44414370, 49639590, 54624570, 60090030, 67897830, 77597520, 87297210, 96996900, 107056950, 118107990};
114 : #define Npo2 23
115 : static unsigned int lpo2[Npo2] = {12, 30, 60, 120, 240, 510, 1020, 2310,
116 : 4620, 9240, 19110, 39270, 79170, 158340,
117 : 324870, 690690, 1345890, 2852850, 5705700,
118 : 11741730, 23130030, 48498450, 96996900};
119 :
120 1566 : unsigned long i, d1 = 0, d2 = 0, dF = 0, phid, k, maxN;
121 : mpz_t j, t, i0, i1;
122 1566 : int r = 0;
123 :
124 1566 : if (mpz_cmp (B2, B2min) < 0)
125 : {
126 : /* No stage 2. Set relevant parameters to 0. Leave B2, B2min the same */
127 100 : *finalk = 0;
128 100 : *finaldF = 0;
129 100 : return 0;
130 : }
131 :
132 1466 : mpz_init (i0);
133 1466 : mpz_init (i1);
134 1466 : mpz_init (j);
135 1466 : mpz_init (t);
136 1466 : k = *finalk; /* User specified k value passed in via finalk */
137 :
138 : /* Look for largest dF we can use while satisfying the maxmem parameter */
139 1466 : maxN = (po2) ? Npo2 : N;
140 1466 : if (maxmem != 0.)
141 : {
142 283 : for (i = 0; i < maxN; i++)
143 : {
144 283 : int lg_dF, sp_num = 0;
145 : double memory;
146 :
147 283 : d1 = (po2) ? lpo2[i] : l[i];
148 283 : phid = eulerphi (d1) / 2;
149 283 : dF = (po2) ? 1U << ceil_log2 (phid) : phid;
150 283 : lg_dF = ceil_log2 (dF);
151 :
152 283 : if (use_ntt)
153 200 : sp_num = (2 * mpz_sizeinbase (modulus->orig_modulus, 2) + lg_dF) /
154 200 : SP_NUMB_BITS + 4;
155 :
156 283 : memory = memory_use (dF, sp_num, (treefile) ? 0 : lg_dF, modulus);
157 283 : outputf (OUTPUT_DEVVERBOSE,
158 : "Estimated mem for dF = %.0d, sp_num = %d: %.0f\n",
159 : dF, sp_num, memory);
160 283 : if (memory > maxmem)
161 21 : break;
162 : }
163 21 : maxN = i;
164 : }
165 :
166 13254 : for (i = 0; i < maxN; i++)
167 : {
168 13233 : d1 = (po2) ? lpo2[i] : l[i];
169 13233 : phid = eulerphi (d1) / 2;
170 13233 : dF = (po2) ? 1U << ceil_log2 (phid) : phid;
171 : /* Look for smallest prime < 25 that does not divide d1 */
172 : /* The caller can force d2 = 1 by setting root_params->d2 != 0 */
173 13233 : d2 = 1;
174 13233 : if (root_params->d2 == 0)
175 34539 : for (d2 = 5; d2 < 25; d2 += 2)
176 : {
177 34539 : if (d2 % 3 == 0)
178 3848 : continue;
179 30691 : if (d1 % d2 > 0)
180 13233 : break;
181 : }
182 :
183 13233 : if (d2 >= 25 || d2 - 1 > dF)
184 3007 : d2 = 1;
185 :
186 : #if 0
187 : /* The code to init roots of G can handle negative i0 now. */
188 : if (d2 > 1 && mpz_cmp_ui (B2min, (d1 - 1) * d2 - d1) <= 0)
189 : d2 = 1; /* Would make i0 < 0 */
190 : #endif
191 :
192 13233 : mpz_set_ui (i0, d1 - 1);
193 13233 : mpz_mul_ui (i0, i0, d2);
194 13233 : mpz_set (j, B2);
195 13233 : mpz_add (i1, j, i0); /* i1 = B2 + (d1 - 1) * d2 */
196 13233 : mpz_set (j, B2min);
197 13233 : mpz_sub (i0, j, i0); /* i0 = B2min - (d1 - 1) * d2 */
198 13233 : mpz_cdiv_q_ui (i0, i0, d1); /* i0 = ceil ((B2min - (d1 - 1) * d2) / d1) */
199 13233 : mpz_fdiv_q_ui (i1, i1, d1); /* i1 = floor ((B2 + (d1 - 1) * d2) / d1) */
200 :
201 : /* How many roots of G will we need ? */
202 13233 : mpz_sub (j, i1, i0);
203 13233 : mpz_add_ui (j, j, 1);
204 :
205 : /* Integer multiples of d2 are skipped (if d2 > 1) */
206 13233 : if (d2 > 1)
207 : {
208 10226 : mpz_fdiv_q_ui (t, i1, d2);
209 10226 : mpz_sub (j, j, t);
210 10226 : mpz_fdiv_q_ui (t, i0, d2);
211 10226 : mpz_add (j, j, t); /* j -= floor (i1 / d2) - floor (i0 / d2) */
212 : }
213 :
214 : /* How many blocks will we need ? Divide lines by dF, rounding up */
215 13233 : mpz_cdiv_q_ui (j, j, dF);
216 :
217 13233 : if ((k != ECM_DEFAULT_K && mpz_cmp_ui (j, k) <= 0) ||
218 11830 : (k == ECM_DEFAULT_K && mpz_cmp_ui (j, (po2) ? 6 : 2) <= 0))
219 : break;
220 : }
221 :
222 1466 : if (i == maxN)
223 : {
224 21 : if (k != ECM_DEFAULT_K)
225 : {
226 : /* The user asked for a specific k and we couldn't satisfy the
227 : condition. Nothing we can do ... */
228 10 : outputf (OUTPUT_ERROR, "Error: too large step 2 bound, increase -k\n");
229 10 : r = ECM_ERROR;
230 10 : goto clear_and_exit;
231 : }
232 11 : else if (!mpz_fits_ulong_p (j))
233 : {
234 : /* Can't fit the number of blocks in an unsigned long. Nothing we
235 : can do ... */
236 10 : outputf (OUTPUT_ERROR, "Error: stage 2 interval too large, cannot "
237 : "generate suitable parameters.\nTry a smaller B2 value.\n");
238 10 : r = ECM_ERROR;
239 10 : goto clear_and_exit;
240 : }
241 1 : if (maxN == 0)
242 : {
243 : /* We can't do a stage 2 at all with the memory the user allowed.
244 : Nothing we can do ... */
245 0 : outputf (OUTPUT_ERROR, "Error: stage 2 not possible with memory "
246 : "allowed by -maxmem.\n");
247 0 : r = ECM_ERROR;
248 0 : goto clear_and_exit;
249 : }
250 : /* else: We can fit the number of blocks into an unsigned int, so we use
251 : it. This may be a very large value for huge B2-B2min, the user
252 : is going to notice sooner or later */
253 : }
254 :
255 : /* If the user specified a number of blocks, we'll use that no matter what.
256 : Since j may be smaller than k, this may increase the B2 limit */
257 1446 : if (k == ECM_DEFAULT_K)
258 1325 : k = mpz_get_ui (j);
259 :
260 : /* Now that we have the number of blocks, compute real i1. There will be
261 : k * dF roots of G computed, starting at i0, skipping all that are not
262 : coprime to d2. While d2 is prime, that means: are not multiples of d2.
263 : Hence we want i1 so that
264 : i1 - floor(i1 / d2) - i0 + ceil((i0 / d2) == k * dF
265 : i1 - floor(i1 / d2) == k * dF + i0 - ceil((i0 / d2)
266 : */
267 :
268 1446 : mpz_set_ui (j, k);
269 1446 : mpz_mul_ui (j, j, dF);
270 1446 : if (d2 == 1)
271 : {
272 103 : mpz_add (i1, i0, j);
273 103 : mpz_sub_ui (i1, i1, 1);
274 : }
275 : else
276 : {
277 1343 : mpz_add (j, j, i0);
278 1343 : mpz_cdiv_q_ui (t, i0, d2);
279 1343 : mpz_sub (j, j, t); /* j = k * dF + i0 - ceil((i0 / d2) */
280 1343 : mpz_fdiv_qr_ui (j, t, j, d2 - 1);
281 1343 : mpz_mul_ui (j, j, d2);
282 1343 : mpz_add (i1, j, t);
283 : }
284 :
285 1446 : root_params->d1 = d1;
286 1446 : root_params->d2 = d2;
287 1446 : mpz_set (root_params->i0, i0);
288 1446 : *finaldF = dF;
289 1446 : *finalk = k;
290 :
291 : /* We want B2' the largest integer that satisfies
292 : i1 = floor ((B2' + (d1 - 1) * d2) / d1)
293 : = floor ((B2'-d2)/d1) + d2
294 : i1 - d2 = floor ((B2'-d2)/d1)
295 : (B2'-d2)/d1 < i1-d2+1
296 : B2'-d2 < (i1-d2+1) * d1
297 : B2' < (i1-d2+1) * d1 + d2
298 : B2' = (i1-d2+1) * d1 + d2 - 1
299 : */
300 :
301 1446 : mpz_sub_ui (i1, i1, d2 - 1);
302 1446 : mpz_mul_ui (B2, i1, d1);
303 1446 : mpz_add_ui (B2, B2, d2 - 1);
304 :
305 1466 : clear_and_exit:
306 1466 : mpz_clear (t);
307 1466 : mpz_clear (j);
308 1466 : mpz_clear (i1);
309 1466 : mpz_clear (i0);
310 :
311 1466 : return r;
312 : }
|
