Line data Source code
1 : /* Dickman's rho function (to compute probability of success of ecm).
2 :
3 : Copyright 2004, 2005, 2006, 2008, 2009, 2010, 2011, 2012, 2013
4 : Alexander Kruppa, Paul Zimmermann.
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 : /* define TESTDRIVE to compile rho as a stand-alone program, in which case
24 : you need to have libgsl installed */
25 :
26 : #include "config.h"
27 : #if defined(TESTDRIVE)
28 : #define _ISOC99_SOURCE 1
29 : #endif
30 : #if defined(DEBUG_NUMINTEGRATE) || defined(TESTDRIVE)
31 : # include <stdio.h>
32 : #endif
33 : #include <stdlib.h>
34 : #include <math.h>
35 : #if defined(TESTDRIVE)
36 : #include <string.h>
37 : #include <primesieve.h>
38 : #endif
39 : #if defined(TESTDRIVE)
40 : #include <gsl/gsl_math.h>
41 : #include <gsl/gsl_sf_expint.h>
42 : #include <gsl/gsl_integration.h>
43 : #endif
44 : #include "ecm-impl.h"
45 :
46 : /* For Suyama's curves, we have a known torsion factor of 12 = 2^2*3^1, and
47 : an average extra exponent of 1/2 for 2, and 1/3 for 3 due to the probability
48 : that the group order divided by 12 is divisible by 2 or 3, thus on average
49 : we should have 2^2.5*3^1.333 ~ 24.5, however experimentally we have
50 : 2^3.323*3^1.687 ~ 63.9 (see Alexander Kruppa's thesis, Table 5.1 page 96,
51 : row sigma=2, http://tel.archives-ouvertes.fr/tel-00477005/en/).
52 : The exp(ECM_EXTRA_SMOOTHNESS) value takes into account the extra
53 : smoothness with respect to a random number. */
54 : #ifndef ECM_EXTRA_SMOOTHNESS
55 : #define ECM_EXTRA_SMOOTHNESS 3.134
56 : #endif
57 :
58 : #define M_PI_SQR 9.869604401089358619 /* Pi^2 */
59 : #define M_PI_SQR_6 1.644934066848226436 /* Pi^2/6 */
60 : /* gsl_math.h defines M_EULER */
61 : #ifndef M_EULER
62 : #define M_EULER 0.577215664901532861
63 : #endif
64 : #define M_EULER_1 0.422784335098467139 /* 1 - Euler */
65 :
66 : #ifndef MAX
67 : #define MAX(x,y) ((x) > (y) ? (x) : (y))
68 : #endif
69 : #ifndef MIN
70 : #define MIN(x,y) ((x) < (y) ? (x) : (y))
71 : #endif
72 :
73 : void rhoinit (int, int); /* used in stage2.c */
74 :
75 : static double *rhotable = NULL;
76 : static int invh = 0;
77 : static double h = 0.;
78 : static int tablemax = 0;
79 : #if defined(TESTDRIVE)
80 : #define PRIME_PI_MAX 10000
81 : #define PRIME_PI_MAP(x) (((x)+1)/2)
82 : /* The number of primes up to i. Use prime_pi[PRIME_PI_MAP(i)].
83 : Only correct for i >= 2. */
84 : static unsigned int prime_pi[PRIME_PI_MAP(PRIME_PI_MAX)+1];
85 : #endif
86 :
87 : /* Fixme: need prime generating funcion without static state variables */
88 : const unsigned char primemap[667] = {
89 : 254, 223, 239, 126, 182, 219, 61, 249, 213, 79, 30, 243, 234, 166, 237, 158,
90 : 230, 12, 211, 211, 59, 221, 89, 165, 106, 103, 146, 189, 120, 30, 166, 86,
91 : 86, 227, 173, 45, 222, 42, 76, 85, 217, 163, 240, 159, 3, 84, 161, 248, 46,
92 : 253, 68, 233, 102, 246, 19, 58, 184, 76, 43, 58, 69, 17, 191, 84, 140, 193,
93 : 122, 179, 200, 188, 140, 79, 33, 88, 113, 113, 155, 193, 23, 239, 84, 150,
94 : 26, 8, 229, 131, 140, 70, 114, 251, 174, 101, 146, 143, 88, 135, 210, 146,
95 : 216, 129, 101, 38, 227, 160, 17, 56, 199, 38, 60, 129, 235, 153, 141, 81,
96 : 136, 62, 36, 243, 51, 77, 90, 139, 28, 167, 42, 180, 88, 76, 78, 38, 246,
97 : 25, 130, 220, 131, 195, 44, 241, 56, 2, 181, 205, 205, 2, 178, 74, 148, 12,
98 : 87, 76, 122, 48, 67, 11, 241, 203, 68, 108, 36, 248, 25, 1, 149, 168, 92,
99 : 115, 234, 141, 36, 150, 43, 80, 166, 34, 30, 196, 209, 72, 6, 212, 58, 47,
100 : 116, 156, 7, 106, 5, 136, 191, 104, 21, 46, 96, 85, 227, 183, 81, 152, 8,
101 : 20, 134, 90, 170, 69, 77, 73, 112, 39, 210, 147, 213, 202, 171, 2, 131, 97,
102 : 5, 36, 206, 135, 34, 194, 169, 173, 24, 140, 77, 120, 209, 137, 22, 176, 87,
103 : 199, 98, 162, 192, 52, 36, 82, 174, 90, 64, 50, 141, 33, 8, 67, 52, 182,
104 : 210, 182, 217, 25, 225, 96, 103, 26, 57, 96, 208, 68, 122, 148, 154, 9, 136,
105 : 131, 168, 116, 85, 16, 39, 161, 93, 104, 30, 35, 200, 50, 224, 25, 3, 68,
106 : 115, 72, 177, 56, 195, 230, 42, 87, 97, 152, 181, 28, 10, 104, 197, 129,
107 : 143, 172, 2, 41, 26, 71, 227, 148, 17, 78, 100, 46, 20, 203, 61, 220, 20,
108 : 197, 6, 16, 233, 41, 177, 130, 233, 48, 71, 227, 52, 25, 195, 37, 10, 48,
109 : 48, 180, 108, 193, 229, 70, 68, 216, 142, 76, 93, 34, 36, 112, 120, 146,
110 : 137, 129, 130, 86, 38, 27, 134, 233, 8, 165, 0, 211, 195, 41, 176, 194, 74,
111 : 16, 178, 89, 56, 161, 29, 66, 96, 199, 34, 39, 140, 200, 68, 26, 198, 139,
112 : 130, 129, 26, 70, 16, 166, 49, 9, 240, 84, 47, 24, 210, 216, 169, 21, 6, 46,
113 : 12, 246, 192, 14, 80, 145, 205, 38, 193, 24, 56, 101, 25, 195, 86, 147, 139,
114 : 42, 45, 214, 132, 74, 97, 10, 165, 44, 9, 224, 118, 196, 106, 60, 216, 8,
115 : 232, 20, 102, 27, 176, 164, 2, 99, 54, 16, 49, 7, 213, 146, 72, 66, 18, 195,
116 : 138, 160, 159, 45, 116, 164, 130, 133, 120, 92, 13, 24, 176, 97, 20, 29, 2,
117 : 232, 24, 18, 193, 1, 73, 28, 131, 48, 103, 51, 161, 136, 216, 15, 12, 244,
118 : 152, 136, 88, 215, 102, 66, 71, 177, 22, 168, 150, 8, 24, 65, 89, 21, 181,
119 : 68, 42, 82, 225, 179, 170, 161, 89, 69, 98, 85, 24, 17, 165, 12, 163, 60,
120 : 103, 0, 190, 84, 214, 10, 32, 54, 107, 130, 12, 21, 8, 126, 86, 145, 1, 120,
121 : 208, 97, 10, 132, 168, 44, 1, 87, 14, 86, 160, 80, 11, 152, 140, 71, 108,
122 : 32, 99, 16, 196, 9, 228, 12, 87, 136, 11, 117, 11, 194, 82, 130, 194, 57,
123 : 36, 2, 44, 86, 37, 122, 49, 41, 214, 163, 32, 225, 177, 24, 176, 12, 138,
124 : 50, 193, 17, 50, 9, 197, 173, 48, 55, 8, 188, 145, 130, 207, 32, 37, 107,
125 : 156, 48, 143, 68, 38, 70, 106, 7, 73, 142, 9, 88, 16, 2, 37, 197, 196, 66,
126 : 90, 128, 160, 128, 60, 144, 40, 100, 20, 225, 3, 132, 81, 12, 46, 163, 138,
127 : 164, 8, 192, 71, 126, 211, 43, 3, 205, 84, 42, 0, 4, 179, 146, 108, 66, 41,
128 : 76, 131, 193, 146, 204, 28};
129 :
130 : #ifdef TESTDRIVE
131 : unsigned long
132 : gcd (unsigned long a, unsigned long b)
133 : {
134 : unsigned long t;
135 :
136 : while (b != 0)
137 : {
138 : t = a % b;
139 : a = b;
140 : b = t;
141 : }
142 :
143 : return a;
144 : }
145 :
146 : unsigned long
147 : eulerphi (unsigned long n)
148 : {
149 : unsigned long phi = 1, p;
150 :
151 : for (p = 2; p * p <= n; p += 2)
152 : {
153 : if (n % p == 0)
154 : {
155 : phi *= p - 1;
156 : n /= p;
157 : while (n % p == 0)
158 : {
159 : phi *= p;
160 : n /= p;
161 : }
162 : }
163 :
164 : if (p == 2)
165 : p--;
166 : }
167 :
168 : /* now n is prime */
169 :
170 : return (n == 1) ? phi : phi * (n - 1);
171 : }
172 :
173 :
174 : /* The number of positive integers up to x that have no prime factor <= y,
175 : for x >= y >= 2. Uses Buchstab's identity */
176 : unsigned long
177 : Buchstab_Phi(unsigned long x, unsigned long y)
178 : {
179 : unsigned long p, s;
180 : primesieve_iterator pg[1];
181 :
182 : if (x < 1)
183 : return 0;
184 : if (x <= y)
185 : return 1;
186 : #if 0
187 : if (x < y^2)
188 : return(1 + primepi(x) - primepi (y)));
189 : #endif
190 :
191 : s = 1;
192 : primesieve_init (pg);
193 : primesieve_skipto (pg, y, x+1);
194 : for (p = primesieve_next_prime (pg); p <= x; p = primesieve_next_prime (pg))
195 : s += Buchstab_Phi(x / p, p - 1);
196 : return (s);
197 : }
198 :
199 :
200 : /* The number of positive integers up to x that have no prime factor
201 : greter than y, for x >= y >= 2. Uses Buchstab's identity */
202 : unsigned long
203 : Buchstab_Psi(const unsigned long x, const unsigned long y)
204 : {
205 : unsigned long r, p;
206 : primesieve_iterator pg[1];
207 :
208 : if (x <= y)
209 : return (x);
210 :
211 : if (y == 1UL)
212 : return (1);
213 :
214 : /* If y^2 > x, then
215 : Psi(x,y) = x - \sum_{y < p < x, p prime} floor(x/p)
216 :
217 : We separate the sum into ranges where floor(x/p) = k,
218 : which is x/(k+1) < p <= x/k.
219 : We also need to satisfy y < p, so we need k < x/y - 1,
220 : or k_max = ceil (x/y) - 2.
221 : The primes y < p <= x/(k_max + 1) are summed separately. */
222 : if (x <= PRIME_PI_MAX && x < y * y)
223 : {
224 : unsigned long kmax = x / y - 1;
225 : unsigned long s1, s2, k;
226 :
227 : s1 = (kmax + 1) * (prime_pi [PRIME_PI_MAP(x / (kmax + 1))] -
228 : prime_pi [PRIME_PI_MAP(y)]);
229 : s2 = 0;
230 : for (k = 1; k <= kmax; k++)
231 : s2 += prime_pi[PRIME_PI_MAP(x / k)];
232 : s2 -= kmax * prime_pi [PRIME_PI_MAP(x / (kmax+1))];
233 : return (x - s1 - s2);
234 : }
235 :
236 : r = 1;
237 : primesieve_init (pg);
238 : for (p = primesieve_next_prime (pg); p <= y; p = primesieve_next_prime (pg))
239 : r += Buchstab_Psi (x / p, p);
240 : return (r);
241 : }
242 :
243 : #endif /* TESTDRIVE */
244 :
245 :
246 : #if defined(TESTDRIVE)
247 : static double
248 : Li (const double x)
249 : {
250 : return (- gsl_sf_expint_E1 (- log(x)));
251 : }
252 : #endif
253 :
254 : /*
255 : Evaluate dilogarithm via the sum
256 : \Li_{2}(z)=\sum_{k=1}^{\infty} \frac{z^k}{k^2},
257 : see http://mathworld.wolfram.com/Dilogarithm.html
258 : Assumes |z| <= 0.5, for which the sum converges quickly.
259 : */
260 :
261 : static double
262 36620 : dilog_series (const double z)
263 : {
264 36620 : double r = 0.0, zk; /* zk = z^k */
265 : int k, k2; /* k2 = k^2 */
266 : /* Doubles have 53 bits in significand, with |z| <= 0.5 the k+1-st term
267 : is <= 1/(2^k k^2) of the result, so 44 terms should do */
268 1647900 : for (k = 1, k2 = 1, zk = z; k <= 44; k2 += 2 * k + 1, k++, zk *= z)
269 1611280 : r += zk / (double) k2;
270 :
271 36620 : return r;
272 : }
273 :
274 : static double
275 36620 : dilog (double x)
276 : {
277 : ASSERT(x <= -1.0); /* dilog(1-x) is called from rhoexact for 2 < x <= 3 */
278 :
279 36620 : if (x <= -2.0)
280 0 : return -dilog_series (1./x) - M_PI_SQR_6 - 0.5 * log(-1./x) * log(-1./x);
281 : else /* x <= -1.0 */
282 : {
283 : /* L2(z) = -L2(1 - z) + 1/6 * Pi^2 - ln(1 - z)*ln(z)
284 : L2(z) = -L2(1/z) - 1/6 * Pi^2 - 0.5*ln^2(-1/z)
285 : ->
286 : L2(z) = -(-L2(1/(1-z)) - 1/6 * Pi^2 - 0.5*ln^2(-1/(1-z))) + 1/6 * Pi^2 - ln(1 - z)*ln(z)
287 : = L2(1/(1-z)) - 1/6 * Pi^2 + 0.5*ln(1 - z)^2 - ln(1 - z)*ln(-z)
288 : z in [-1, -2) -> 1/(1-z) in [1/2, 1/3)
289 : */
290 36620 : double log1x = log (1. - x);
291 36620 : return dilog_series (1. / (1. - x))
292 36620 : - M_PI_SQR_6 + log1x * (0.5 * log1x - log (-x));
293 : }
294 : }
295 :
296 : #if 0
297 : static double
298 : L2 (double x)
299 : {
300 : return log (x) * (1 - log (x-1)) + M_PI_SQR_6 - dilog (1 - x);
301 : }
302 : #endif
303 :
304 : static double
305 99321 : rhoexact (double x)
306 : {
307 : ASSERT(x <= 3.);
308 99321 : if (x <= 0.)
309 122 : return 0.;
310 99199 : else if (x <= 1.)
311 31237 : return 1.;
312 67962 : else if (x <= 2.)
313 31342 : return 1. - log (x);
314 : else /* 2 < x <= 3 thus -2 <= 1-x < -1 */
315 36620 : return 1. - log (x) * (1. - log (x - 1.)) + dilog (1. - x) + 0.5 * M_PI_SQR_6;
316 : }
317 :
318 :
319 : #if defined(TESTDRIVE)
320 :
321 : /* The Buchstab omega(x) function, exact for x <= 4 where it can be
322 : evaluated without numerical integration, and approximated by
323 : exp(gamma) for larger x. */
324 :
325 : static double
326 : Buchstab_omega (const double x)
327 : {
328 : /* magic = dilog(-1) + 1 = Pi^2/12 + 1 */
329 : const double magic = 1.82246703342411321824;
330 :
331 : if (x < 1.) return (0.);
332 : if (x <= 2.) return (1. / x);
333 : if (x <= 3.) return ((log (x - 1.) + 1.) / x);
334 : if (x <= 4.)
335 : return ((dilog(2. - x) + (1. + log(x - 2.)) * log(x - 1.) + magic) / x);
336 :
337 : /* If argument is out of range, return the limiting value for
338 : $x->\infty$: e^-gamma.
339 : For x only a little larger than 4., this has relative error 2.2e-6,
340 : for larger x the error rapidly drops further */
341 :
342 : return 0.56145948356688516982;
343 : }
344 :
345 : #endif
346 :
347 : void
348 214 : rhoinit (int parm_invh, int parm_tablemax)
349 : {
350 : int i;
351 :
352 214 : if (parm_invh == invh && parm_tablemax == tablemax)
353 0 : return;
354 :
355 214 : if (rhotable != NULL)
356 : {
357 92 : free (rhotable);
358 92 : rhotable = NULL;
359 92 : invh = 0;
360 92 : h = 0.;
361 92 : tablemax = 0;
362 : }
363 :
364 : /* The integration below expects 3 * invh > 4 */
365 214 : if (parm_tablemax == 0 || parm_invh < 2)
366 92 : return;
367 :
368 122 : invh = parm_invh;
369 122 : h = 1. / (double) invh;
370 122 : tablemax = parm_tablemax;
371 :
372 122 : rhotable = (double *) malloc (parm_invh * parm_tablemax * sizeof (double));
373 122 : ASSERT_ALWAYS(rhotable != NULL);
374 :
375 93818 : for (i = 0; i < (3 < parm_tablemax ? 3 : parm_tablemax) * invh; i++)
376 93696 : rhotable[i] = rhoexact (i * h);
377 :
378 218746 : for (i = 3 * invh; i < parm_tablemax * invh; i++)
379 : {
380 : /* rho(i*h) = 1 - \int_{1}^{i*h} rho(x-1)/x dx
381 : = rho((i-4)*h) - \int_{(i-4)*h}^{i*h} rho(x-1)/x dx */
382 :
383 218624 : rhotable[i] = rhotable[i - 4] - 2. / 45. * (
384 218624 : 7. * rhotable[i - invh - 4] / (double)(i - 4)
385 218624 : + 32. * rhotable[i - invh - 3] / (double)(i - 3)
386 218624 : + 12. * rhotable[i - invh - 2] / (double)(i - 2)
387 218624 : + 32. * rhotable[i - invh - 1] / (double)(i - 1)
388 218624 : + 7. * rhotable[i - invh] / (double)i );
389 218624 : if (rhotable[i] < 0.)
390 : {
391 : #ifndef DEBUG_NUMINTEGRATE
392 2928 : rhotable[i] = 0.;
393 : #else
394 : printf (stderr, "rhoinit: rhotable[%d] = %.16f\n", i,
395 : rhotable[i]);
396 : exit (EXIT_FAILURE);
397 : #endif
398 : }
399 : }
400 : }
401 :
402 : /* assumes alpha < tablemax */
403 : static double
404 207846 : dickmanrho (double alpha)
405 : {
406 : ASSERT(alpha < tablemax);
407 :
408 207846 : if (alpha <= 3.)
409 5625 : return rhoexact (alpha);
410 : {
411 202221 : int a = floor (alpha * invh);
412 202221 : double rho1 = rhotable[a];
413 202221 : double rho2 = (a + 1) < tablemax * invh ? rhotable[a + 1] : 0;
414 202221 : return rho1 + (rho2 - rho1) * (alpha * invh - (double) a);
415 : }
416 : }
417 :
418 : #if 0
419 : static double
420 : dickmanrhosigma (double alpha, double x)
421 : {
422 : if (alpha <= 0.)
423 : return 0.;
424 : if (alpha <= 1.)
425 : return 1.;
426 : if (alpha < tablemax)
427 : return dickmanrho (alpha) + M_EULER_1 * dickmanrho (alpha - 1.) / log (x);
428 :
429 : return 0.;
430 : }
431 :
432 : static double
433 : dickmanrhosigma_i (int ai, double x)
434 : {
435 : if (ai <= 0)
436 : return 0.;
437 : if (ai <= invh)
438 : return 1.;
439 : if (ai < tablemax * invh)
440 : return rhotable[ai] - M_EULER * rhotable[ai - invh] / log(x);
441 :
442 : return 0.;
443 : }
444 : #endif
445 :
446 : /* return the value of the "local" Dickman rho function, for numbers near x
447 : (as opposed to numbers <= x for the original Dickman rho function).
448 : Reference: PhD thesis of Alexander Kruppa,
449 : http://docnum.univ-lorraine.fr/public/SCD_T_2010_0054_KRUPPA.pdf,
450 : equation (5.6) page 100 */
451 : static double
452 1138430 : dickmanlocal (double alpha, double x)
453 : {
454 1138430 : if (alpha <= 1.)
455 0 : return rhoexact (alpha);
456 1138430 : if (alpha < tablemax)
457 103923 : return dickmanrho (alpha) - M_EULER * dickmanrho (alpha - 1.) / log (x);
458 1034507 : return 0.;
459 : }
460 :
461 : static double
462 665319 : dickmanlocal_i (int ai, double x)
463 : {
464 665319 : if (ai <= 0)
465 0 : return 0.;
466 665319 : if (ai <= invh)
467 76800 : return 1.;
468 588519 : if (ai <= 2 * invh && ai < tablemax * invh)
469 77090 : return rhotable[ai] - M_EULER / log (x);
470 511429 : if (ai < tablemax * invh)
471 : {
472 508665 : double logx = log (x);
473 508665 : return rhotable[ai] - (M_EULER * rhotable[ai - invh]
474 508665 : + M_EULER_1 * rhotable[ai - 2 * invh] / logx) / logx;
475 : }
476 :
477 2764 : return 0.;
478 : }
479 :
480 : static int
481 10208050 : isprime(unsigned long n)
482 : {
483 : unsigned int r;
484 :
485 10208050 : if (n % 2 == 0)
486 5104050 : return (n == 2);
487 5104000 : if (n % 3 == 0)
488 1701200 : return (n == 3);
489 3402800 : if (n % 5 == 0)
490 680400 : return (n == 5);
491 :
492 2722400 : if (n / 30 >= sizeof (primemap))
493 0 : abort();
494 :
495 2722400 : r = n % 30; /* 8 possible values: 1,7,11,13,17,19,23,29 */
496 2722400 : r = (r * 16 + r) / 64; /* maps the 8 values onto 0, ..., 7 */
497 :
498 2722400 : return ((primemap[n / 30] & (1 << r)) != 0);
499 : }
500 :
501 : /* return the sum in Equation (5.10) page 102 of Alexander Kruppa's
502 : PhD thesis */
503 : static double
504 630 : dickmanmu_sum (const unsigned long B1, const unsigned long B2,
505 : const double x)
506 : {
507 630 : double s = 0.;
508 630 : const double inv_logB1 = 1. / log(B1);
509 630 : const double logx = log(x);
510 : unsigned long p;
511 :
512 10208680 : for (p = B1 + 1; p <= B2; p++)
513 10208050 : if (isprime(p))
514 1134260 : s += dickmanlocal ((logx - log(p)) * inv_logB1, x / p) / p;
515 :
516 630 : return s;
517 : }
518 :
519 : /* return the probability that a number < x has its 2nd largest prime factor
520 : less than x^(1/alpha) and its largest prime factor less than x^(beta/alpha)
521 : */
522 : static double
523 1290 : dickmanmu (double alpha, double beta, double x)
524 : {
525 : double a, b, sum;
526 : int ai, bi, i;
527 1290 : ai = ceil ((alpha - beta) * invh);
528 1290 : if (ai > tablemax * invh)
529 648 : ai = tablemax * invh;
530 1290 : a = (double) ai * h;
531 1290 : bi = floor ((alpha - 1.) * invh);
532 1290 : if (bi > tablemax * invh)
533 674 : bi = tablemax * invh;
534 1290 : b = (double) bi * h;
535 1290 : sum = 0.;
536 77543 : for (i = ai + 1; i < bi; i++)
537 76253 : sum += dickmanlocal_i (i, x) / (alpha - i * h);
538 1290 : sum += 0.5 * dickmanlocal_i (ai, x) / (alpha - a);
539 1290 : sum += 0.5 * dickmanlocal_i (bi, x) / (alpha - b);
540 1290 : sum *= h;
541 1290 : sum += (a - alpha + beta) * 0.5 * (dickmanlocal_i (ai, x) / (alpha - a) + dickmanlocal (alpha - beta, x) / beta);
542 1290 : sum += (alpha - 1. - b) * 0.5 * (dickmanlocal (alpha - 1., x) + dickmanlocal_i (bi, x) / (alpha - b));
543 :
544 1290 : return sum;
545 : }
546 :
547 : static double
548 300 : brentsuyama (double B1, double B2, double N, double nr)
549 : {
550 : double a, alpha, beta, sum;
551 : int ai, i;
552 300 : alpha = log (N) / log (B1);
553 300 : beta = log (B2) / log (B1);
554 300 : ai = floor ((alpha - beta) * invh);
555 300 : if (ai > tablemax * invh)
556 60 : ai = tablemax * invh;
557 300 : a = (double) ai * h;
558 300 : sum = 0.;
559 583606 : for (i = 1; i < ai; i++)
560 583306 : sum += dickmanlocal_i (i, N) / (alpha - i * h) * (1 - exp (-nr * pow (B1, (-alpha + i * h))));
561 300 : sum += 0.5 * (1 - exp(-nr / pow (B1, alpha)));
562 300 : sum += 0.5 * dickmanlocal_i (ai, N) / (alpha - a) * (1 - exp(-nr * pow (B1, (-alpha + a))));
563 300 : sum *= h;
564 300 : sum += 0.5 * (alpha - beta - a) * (dickmanlocal_i (ai, N) / (alpha - a) + dickmanlocal (alpha - beta, N) / beta);
565 :
566 300 : return sum;
567 : }
568 :
569 : static double
570 160 : brsudickson (double B1, double B2, double N, double nr, int S)
571 : {
572 : int i, f;
573 : double sum;
574 160 : sum = 0;
575 160 : f = eulerphi (S) / 2;
576 1060 : for (i = 1; i <= S / 2; i++)
577 900 : if (gcd (i, S) == 1)
578 300 : sum += brentsuyama (B1, B2, N, nr * (gcd (i - 1, S) + gcd (i + 1, S) - 4) / 2);
579 :
580 160 : return sum / (double)f;
581 : }
582 :
583 : static double
584 0 : brsupower (double B1, double B2, double N, double nr, int S)
585 : {
586 : int i, f;
587 : double sum;
588 0 : sum = 0;
589 0 : f = eulerphi (S);
590 0 : for (i = 1; i < S; i++)
591 0 : if (gcd (i, S) == 1)
592 0 : sum += brentsuyama (B1, B2, N, nr * (gcd (i - 1, S) - 2));
593 :
594 0 : return sum / (double)f;
595 : }
596 :
597 : /* Assume N is as likely smooth as a number around N/exp(delta) */
598 :
599 : static double
600 1310 : prob (double B1, double B2, double N, double nr, int S, double delta)
601 : {
602 1310 : const double sumthresh = 20000.;
603 : double alpha, beta, stage1, stage2, brsu;
604 1310 : const double effN = N / exp (delta);
605 :
606 : ASSERT(rhotable != NULL);
607 :
608 : /* What to do if rhotable is not initialised and asserting is not enabled?
609 : For now, bail out with 0. result. Not really pretty, either */
610 1310 : if (rhotable == NULL)
611 0 : return 0.;
612 :
613 1310 : if (B1 < 2. || N <= 1.)
614 20 : return 0.;
615 :
616 1290 : if (effN <= B1)
617 0 : return 1.;
618 :
619 : #ifdef TESTDRIVE
620 : printf ("B1 = %f, B2 = %f, N = %.0f, nr = %f, S = %d\n", B1, B2, N, nr, S);
621 : #endif
622 :
623 1290 : alpha = log (effN) / log (B1);
624 1290 : stage1 = dickmanlocal (alpha, effN);
625 1290 : stage2 = 0.;
626 1290 : if (B2 > B1)
627 : {
628 1290 : if (B1 < sumthresh)
629 : {
630 630 : stage2 += dickmanmu_sum (B1, MIN(B2, sumthresh), effN);
631 630 : beta = log (B2) / log (MIN(B2, sumthresh));
632 : }
633 : else
634 660 : beta = log (B2) / log (B1);
635 :
636 1290 : if (beta > 1.)
637 1290 : stage2 += dickmanmu (alpha, beta, effN);
638 : }
639 1290 : brsu = 0.;
640 1290 : if (S < -1)
641 160 : brsu = brsudickson (B1, B2, effN, nr, -S * 2);
642 1290 : if (S > 1)
643 0 : brsu = brsupower (B1, B2, effN, nr, S * 2);
644 :
645 : #ifdef TESTDRIVE
646 : printf ("stage 1 : %f, stage 2 : %f, Brent-Suyama : %f\n", stage1, stage2, brsu);
647 : #endif
648 :
649 1290 : return (stage1 + stage2 + brsu) > 0. ? (stage1 + stage2 + brsu) : 0.;
650 : }
651 :
652 : double
653 660 : ecmprob (double B1, double B2, double N, double nr, int S)
654 : {
655 660 : return prob (B1, B2, N, nr, S, ECM_EXTRA_SMOOTHNESS);
656 : }
657 :
658 : /* see Willemien Ekkelkamp's Phd thesis:
659 : https://openaccess.leidenuniv.nl/bitstream/handle/1887/14567/proefschrift_041109.pdf?sequence=2:
660 : Bach-Peralta formula page 12
661 : Corollary 4 page 18 with a 2nd-order term */
662 : double
663 650 : pm1prob (double B1, double B2, double N, double nr, int S, const mpz_t go)
664 : {
665 : mpz_t cof;
666 : /* A prime power q^k divides p-1, p prime, with probability 1/(q^k-q^(k-1))
667 : not with probability 1/q^k as for random numbers. This is taken into
668 : account by the "smoothness" value here; a prime p-1 is about as likely
669 : smooth as a random number around (p-1)/exp(smoothness).
670 : smoothness = \sum_{q in Primes} log(q)/(q-1)^2 */
671 : /* Note that this routine is also called for P+1, where we assume the same
672 : behaviour as with P-1. However, if x0=6/5, Kruppa writes in his PhD
673 : thesis that we get smoothness = 1.92012; with x0=2/7, we get
674 : smoothness = 2.05093. */
675 650 : double smoothness = 1.2269688;
676 : unsigned long i;
677 :
678 650 : if (go != NULL && mpz_cmp_ui (go, 1UL) > 0)
679 : {
680 : double res;
681 50 : mpz_init (cof);
682 50 : mpz_set (cof, go);
683 4950 : for (i = 2; i < 100; i++)
684 4900 : if (mpz_divisible_ui_p (cof, i))
685 : {
686 : /* If we know that q divides p-1 with probability 1, we need to
687 : adjust the smoothness parameter */
688 50 : smoothness -= log ((double) i) / (double) ((i-1)*(i-1));
689 : /* printf ("pm1prob: Dividing out %lu\n", i); */
690 100 : while (mpz_divisible_ui_p (cof, i))
691 50 : mpz_tdiv_q_ui (cof, cof, i);
692 : }
693 : /* printf ("pm1prob: smoothness after dividing out go primes < 100: %f\n",
694 : smoothness); */
695 50 : res = prob (B1, B2, N, nr, S, smoothness + log(mpz_get_d (cof)));
696 50 : mpz_clear (cof);
697 50 : return res;
698 : }
699 :
700 600 : return prob (B1, B2, N, nr, S, smoothness);
701 : }
702 :
703 : #if defined(TESTDRIVE)
704 :
705 : /* Compute probability for primes p == r (mod m) */
706 :
707 : static double
708 : pm1prob_rm (double B1, double B2, double N, double nr, int S, unsigned long r,
709 : unsigned long m)
710 : {
711 : unsigned long cof;
712 : double smoothness = 1.2269688;
713 : unsigned long p;
714 :
715 : cof = m;
716 :
717 : for (p = 2UL; p < 100UL; p++)
718 : if (cof % p == 0UL) /* For each prime in m */
719 : {
720 : unsigned long cof_r, k, i;
721 : /* Divisibility by i is determined by r and m. We need to
722 : adjust the smoothness parameter. In P-1, we had estimated the
723 : expected value for the exponent of p as p/(p-1)^2. Undo that. */
724 : smoothness -= (double)p / ((p-1)*(p-1)) * log ((double) p);
725 : /* The expected value for the exponent of this prime is k s.t.
726 : p^k || r, plus 1/(p-1) if p^k || m as well */
727 : cof_r = gcd (r - 1UL, m);
728 : for (k = 0UL; cof_r % p == 0UL; k++)
729 : cof_r /= p;
730 : smoothness += k * log ((double) p);
731 :
732 : cof_r = m;
733 : for (i = 0UL; cof_r % p == 0UL; i++)
734 : cof_r /= p;
735 :
736 : if (i == k)
737 : smoothness += (1./(p - 1.) * log ((double) p));
738 :
739 : while (cof % p == 0UL)
740 : cof /= p;
741 : printf ("pm1prob_rm: p = %lu, k = %lu, i = %lu, new smoothness = %f\n",
742 : p, i, k, smoothness);
743 : }
744 :
745 : return prob (B1, B2, N, nr, S, smoothness);
746 : }
747 :
748 :
749 : /* The \Phi(x,y) function gives the number of natural numbers <= x
750 : that have no prime factor <= y, see Tenenbaum,
751 : "Introduction the analytical and probabilistic number theory", III.6.
752 : This function estimates the \Phi(x,y) function via eq. (48) of the 1st
753 : edition resp. equation (6.49) of the 3rd edition of Tenenbaum's book. */
754 :
755 : static double
756 : integrand1 (double x, double *y)
757 : {
758 : return pow (*y, x) / x * log(x-1.);
759 : }
760 :
761 :
762 : static double
763 : integrand2 (double v, double *y)
764 : {
765 : return Buchstab_omega (v) * pow (*y, v);
766 : }
767 :
768 :
769 : /* Return approximate number of integers n with x1 < n <= x2
770 : that have no prime factor <= y */
771 :
772 : double
773 : no_small_prime (double x1, double x2, double y)
774 : {
775 : double u1, u2;
776 : ASSERT (x1 >= 2.);
777 : ASSERT (x2 >= x1);
778 : ASSERT (y >= 2.);
779 : if (x1 == x2 || x2 <= y)
780 : return 0.;
781 : if (x1 < y)
782 : x1 = y;
783 :
784 : u1 = log(x1)/log(y);
785 : u2 = log(x2)/log(y);
786 :
787 : /* If no prime factors <= sqrt(x2), numbers must be a primes > y */
788 : if (x2 <= y*y)
789 : return (Li(x2) - Li(x1));
790 :
791 : if (u2 <= 3)
792 : {
793 : double r, abserr;
794 : size_t neval;
795 : gsl_function f;
796 :
797 : f.function = (double (*) (double, void *)) &integrand1;
798 : f.params = &y;
799 :
800 : /* intnum(v=1,u,buchstab(v)*y^v) */
801 :
802 : /* First part: intnum(v=u1, u, y^v/v*log(v-1.)) */
803 : gsl_integration_qng (&f, MAX(u1, 2.) , u2, 0., 0.001, &r, &abserr, &neval);
804 :
805 : /* Second part: intnum(v=u1, u2, y^v/v) = Li(x2) - Li(x1) */
806 : r += Li (x2) - Li (x1);
807 :
808 : return r;
809 : }
810 :
811 : {
812 : double r, abserr;
813 : size_t neval;
814 : gsl_function f;
815 :
816 : f.function = (double (*) (double, void *)) &integrand2;
817 : f.params = &y;
818 :
819 : gsl_integration_qng (&f, u1, u2, 0., 0.001, &r, &abserr, &neval);
820 : return r;
821 : }
822 : }
823 :
824 :
825 : static double
826 : integrand3 (double p, double *param)
827 : {
828 : const double x1 = param[0];
829 : const double x2 = param[1];
830 : const double y = param[2];
831 :
832 : return no_small_prime (x1 / p, x2 / p, y) / log(p);
833 : }
834 :
835 :
836 : double
837 : no_small_prime_factor (const double x1, const double x2, const double y,
838 : const double z1, const double z2)
839 : {
840 : double r, abserr, param[3];
841 : size_t neval;
842 : gsl_function f;
843 :
844 : param[0] = x1;
845 : param[1] = x2;
846 : param[2] = y;
847 : f.function = (double (*) (double, void *)) &integrand3;
848 : f.params = ¶m;
849 :
850 : gsl_integration_qng (&f, z1, z2, 0., 0.01, &r, &abserr, &neval);
851 :
852 : return r;
853 : }
854 :
855 : #endif
856 :
857 :
858 : #ifdef TESTDRIVE
859 : int
860 : main (int argc, char **argv)
861 : {
862 : double B1, B2, N, nr, r, m;
863 : int S;
864 : unsigned long p, i, pi;
865 : primesieve_iterator pg[1];
866 :
867 : primesieve_init (pg);
868 : i = pi = 0;
869 : for (p = primesieve_next_prime (pg); p <= PRIME_PI_MAX; p = primesieve_next_prime (pg))
870 : {
871 : for ( ; i < p; i++)
872 : prime_pi[PRIME_PI_MAP(i)] = pi;
873 : pi++;
874 : }
875 : for ( ; i < p; i++)
876 : prime_pi[PRIME_PI_MAP(i)] = pi;
877 :
878 :
879 : if (argc < 2)
880 : {
881 : printf ("Usage: rho <B1> <B2> <N> <nr> <S> [<r> <m>]\n\n\n");
882 : printf (" Calculate the probability of ECM/P-1 finding a factor near N\n"
883 : " with B1/B2, evaluating nr random distinct points in stage 2,\n"
884 : " with a degree -S Dickson polynomial (if S < 0) or\n"
885 : " S'th power as the Brent-Suyama function\n\n");
886 : printf (" <B1> B1 limit.\n");
887 : printf (" <B2> B2 limit.\n");
888 : printf (" <N> N of similiar size, or number of bits in factor (if < 50).\n");
889 : printf (" <nr> Number of random points evaluated in stage 2.\n");
890 : printf (" <S> Degree of Brent-Suyama polynomial in stage 2.\n");
891 : printf (" [<r> <m>] Limit P-1 to primes p == r (mod m).\n");
892 : return 1;
893 : }
894 :
895 : if (strcmp (argv[1], "-Buchstab_Phi") == 0)
896 : {
897 : unsigned long x, y, r;
898 : if (argc < 4)
899 : {
900 : printf ("-Buchstab_Phi needs x and y paramters\n");
901 : exit (EXIT_FAILURE);
902 : }
903 : x = strtoul (argv[2], NULL, 10);
904 : y = strtoul (argv[3], NULL, 10);
905 : r = Buchstab_Phi (x, y);
906 : printf ("Buchstab_Phi (%lu, %lu) = %lu\n", x, y, r);
907 : exit (EXIT_SUCCESS);
908 : }
909 : else if (strcmp (argv[1], "-Buchstab_Psi") == 0)
910 : {
911 : unsigned long x, y, r;
912 : if (argc < 4)
913 : {
914 : printf ("-Buchstab_Psi needs x and y paramters\n");
915 : exit (EXIT_FAILURE);
916 : }
917 : x = strtoul (argv[2], NULL, 10);
918 : y = strtoul (argv[3], NULL, 10);
919 : r = Buchstab_Psi (x, y);
920 : printf ("Buchstab_Psi (%lu, %lu) = %lu\n", x, y, r);
921 : exit (EXIT_SUCCESS);
922 : }
923 : else if (strcmp (argv[1], "-nsp") == 0)
924 : {
925 : double x1, x2, y, r;
926 :
927 : if (argc < 5)
928 : {
929 : printf ("-nsp needs x1, x2, and y paramters\n");
930 : exit (EXIT_FAILURE);
931 : }
932 : x1 = atof (argv[2]);
933 : x2 = atof (argv[3]);
934 : y = atof (argv[4]);
935 : r = no_small_prime (x1, x2, y);
936 : printf ("no_small_prime(%f, %f, %f) = %f\n", x1, x2, y, r);
937 : exit (EXIT_SUCCESS);
938 : }
939 : else if (strcmp (argv[1], "-nspf") == 0)
940 : {
941 : double x1, x2, y, z1, z2, r;
942 :
943 : if (argc < 7)
944 : {
945 : printf ("-nspf needs x1, x2, y, z1, and z2 paramters\n");
946 : exit (EXIT_FAILURE);
947 : }
948 : x1 = atof (argv[2]);
949 : x2 = atof (argv[3]);
950 : y = atof (argv[4]);
951 : z1 = atof (argv[5]);
952 : z2 = atof (argv[6]);
953 : r = no_small_prime_factor (x1, x2, y, z1, z2);
954 : printf ("no_small_prime(%f, %f, %f, %f, %f) = %f\n", x1, x2, y, z1, z2, r);
955 : exit (EXIT_SUCCESS);
956 : }
957 :
958 :
959 : if (argc < 6)
960 : {
961 : printf ("Need 5 or 7 arguments: B1 B2 N nr S [r m]\n");
962 : exit (EXIT_FAILURE);
963 : }
964 :
965 : B1 = atof (argv[1]);
966 : B2 = atof (argv[2]);
967 : N = atof (argv[3]);
968 : nr = atof (argv[4]);
969 : S = atoi (argv[5]);
970 : r = 0; m = 1;
971 : if (argc > 7)
972 : {
973 : r = atoi (argv[6]);
974 : m = atoi (argv[7]);
975 : }
976 :
977 : rhoinit (256, 10);
978 : if (N < 50.)
979 : {
980 : double sum;
981 : sum = ecmprob(B1, B2, exp2 (N), nr, S);
982 : sum += 4. * ecmprob(B1, B2, 3./2. * exp2 (N), nr, S);
983 : sum += ecmprob(B1, B2, 2. * exp2 (N), nr, S);
984 : sum *= 1./6.;
985 : printf ("ECM: %.16f\n", sum);
986 :
987 : sum = pm1prob_rm (B1, B2, exp2 (N), nr, S, r, m);
988 : sum += 4. * pm1prob_rm (B1, B2, 3./2. * exp2 (N), nr, S, r, m);
989 : sum += pm1prob_rm (B1, B2, 2. * exp2 (N), nr, S, r, m);
990 : sum *= 1./6.;
991 : printf ("P-1: %.16f\n", sum);
992 : }
993 : else
994 : {
995 : printf ("ECM: %.16f\n", ecmprob(B1, B2, N, nr, S));
996 : printf ("P-1: %.16f\n", pm1prob_rm (B1, B2, N, nr, S, r, m));
997 : }
998 : rhoinit (0, 0);
999 : return 0;
1000 : }
1001 : #endif
|
