1 #define _DEFAULT_SOURCE // for random()
2 #include <time.h>
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include <math.h>
6 #include <stdbool.h>
7 #include <limits.h>
8 #include <sys/time.h>
9 #include <gmp.h>
10 #include "ecm.h"
11
12 #include "smooth_detect.h"
13
14 #ifndef MIN
15 #define MIN(l,o) ((l) < (o) ? (l) : (o))
16 #endif
17 #ifndef MAX
18 #define MAX(h,i) ((h) > (i) ? (h) : (i))
19 #endif
20
21 // For a bug in ecm ?
22 static double default_B1done;
23
24 #define EA_THRESHOLD 0.8 // Heuristic for early-abort: low = keep many
25
26 // Expected effort to extract of prime of p bits.
27 // This has been computed for B1min = 200.0, but this is rather stable
28 // for interesting values of p.
29 #define MAX_PBIT 100
30 static const double expected_effort[MAX_PBIT] = {
31 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, // 0 - 9
32 200, 200, 200, 200, 200, 200, 200, 211, 219, 243, // 10 - 19
33 296, 280, 351, 380, 482, 588, 653, 719, 834, 985, // 20 - 29
34 1280, 1489, 2352, 2324, 2714, 3813, 5141, 5891, 7986, 8816, // 30 - 39
35 10087, 12264, 14191, 18827, 25633, 25491, 37803, 39392, 44290, 51925, // 40-49
36 87203, 79943, 110007, 121644, 147602, // 50 - 54
37 174995, 199245, 257190, 279228, 345960, // 55 - 59
38 351682, 530703, 640140, 759310, 775311, // 60 - 64
39 960249, 1267879, 1174122, 1272107, 1589907, // 65 - 69
40 2258437, 2004235, 2903959, 3002629, 3888904, // 70 - 74
41 4373729, 4899345, 5218152, 6269843, 7063446, // 75 - 79
42 9553542, 9891138, 10623352, 13795248, 17574109, // 80 - 84
43 18790448, 23529670, 24757303, 30897420, 31188626, // 85 - 89
44 36647830, 41007546, 45692079, 53881307, 69709984, // 90 - 94
45 81857570, 84194044, 107900749, 94609433, 136660173 // 95 - 99
46 };
47
48
get_time()49 double get_time() {
50 return (double)(clock()) / CLOCKS_PER_SEC;
51 }
52
53
54 ///////////////////////////////////////////////////////////////////////////
55 // Candidate: structure to store a pair of numbers currently being
56 // factored.
57 ///////////////////////////////////////////////////////////////////////////
58
cand_init(cand_t c)59 void cand_init(cand_t c) {
60 mpz_init(c->u0);
61 mpz_init(c->v0);
62 mpz_init(c->u);
63 mpz_init(c->v);
64 }
65
cand_set(cand_t c,const cand_t c0)66 void cand_set(cand_t c, const cand_t c0) {
67 mpz_set(c->u0, c0->u0);
68 mpz_set(c->v0, c0->v0);
69 mpz_set(c->u, c0->u);
70 mpz_set(c->v, c0->v);
71 c->lpu = c0->lpu;
72 c->lpv = c0->lpv;
73 c->effort = c0->effort;
74 c->id = c0->id;
75 }
76
cand_init_copy(cand_t c,const cand_t c0)77 void cand_init_copy(cand_t c, const cand_t c0) {
78 cand_init(c);
79 cand_set(c, c0);
80 }
81
cand_clear(cand_t c)82 void cand_clear(cand_t c) {
83 mpz_clear(c->u0);
84 mpz_clear(c->v0);
85 mpz_clear(c->u);
86 mpz_clear(c->v);
87 }
88
cand_swap(cand_t c,cand_t d)89 void cand_swap(cand_t c, cand_t d) {
90 mpz_swap(c->u0, d->u0);
91 mpz_swap(c->v0, d->v0);
92 mpz_swap(c->u, d->u);
93 mpz_swap(c->v, d->v);
94 unsigned int t;
95 t = c->lpu; c->lpu = d->lpu; d->lpu = t;
96 t = c->lpv; c->lpv = d->lpv; d->lpv = t;
97 double tt;
98 tt = c->effort; c->effort = d->effort; d->effort = tt;
99 unsigned long id;
100 id = c->id; c->id = d->id; d->id = id;
101 }
102
cand_print(const cand_t c)103 void cand_print(const cand_t c) {
104 printf("Candidate id = %lu\n", c->id);
105 gmp_printf("u0=%Zd\n", c->u0);
106 gmp_printf("v0=%Zd\n", c->v0);
107 gmp_printf("u=%Zd (%lu bits) largest prime so far of %lu bits\n",
108 c->u, mpz_sizeinbase (c->u, 2), c->lpu);
109 gmp_printf("v=%Zd (%lu bits) largest prime so far of %lu bits\n",
110 c->v, mpz_sizeinbase (c->v, 2), c->lpv);
111 printf("effort=%.0f\n", c->effort);
112 }
113
cand_is_factored(const cand_t c)114 bool cand_is_factored(const cand_t c) {
115 return (mpz_cmp_ui(c->u, 1) == 0) && (mpz_cmp_ui(c->v, 1) == 0);
116 }
117
118 // if effort is such that all primes up to b-bits have been removed,
119 // and if an unfactored part has less than 3*b bits, then there are only
120 // two factors. If size is more than 2*bound, it can not be smooth.
121 // Heuristic: we consider that if the effort is twice the average for
122 // detecting p bits, then all p bits prime have been removed.
cand_is_probably_not_smooth(const cand_t c,unsigned int bound)123 bool cand_is_probably_not_smooth(const cand_t c, unsigned int bound) {
124 double eff = c->effort;
125 unsigned int bits = 0;
126 while ((bits < MAX_PBIT) && (2*expected_effort[bits] < eff)) {
127 bits++;
128 }
129 if (bits == MAX_PBIT) {
130 return false; // can not conclude; no data for this effort
131 }
132 // bits -= 3; // take some margin
133
134 for (unsigned int k = 2; k < 4; ++k) {
135 unsigned long bu = mpz_sizeinbase(c->u, 2);
136 if ((bu < (k+1)*bits) && (bu > k*bound)) {
137 // printf("Probably not smooth, level %d: %lu bits!\n", k, bu);
138 return true;
139 }
140 unsigned long bv = mpz_sizeinbase(c->v, 2);
141 if ((bv < (k+1)*bits) && (bv > k*bound)) {
142 // printf("Probably not smooth, level %d: %lu bits!\n", k, bv);
143 return true;
144 }
145 }
146 return false;
147 }
148
cand_update_check_prime(cand_t c)149 void cand_update_check_prime(cand_t c) {
150 if (mpz_probab_prime_p(c->u, 1)) {
151 c->lpu = MAX(c->lpu, mpz_sizeinbase(c->u, 2));
152 mpz_set_ui(c->u, 1);
153 }
154 if (mpz_probab_prime_p(c->v, 1)) {
155 c->lpv = MAX(c->lpv, mpz_sizeinbase(c->v, 2));
156 mpz_set_ui(c->v, 1);
157 }
158 }
159
cand_set_original_values(cand_t c,const mpz_t u0,const mpz_t v0,unsigned long id)160 void cand_set_original_values(cand_t c, const mpz_t u0, const mpz_t v0,
161 unsigned long id) {
162 mpz_set(c->u0, u0);
163 mpz_set(c->v0, v0);
164 mpz_set(c->u, u0);
165 mpz_set(c->v, v0);
166 c->lpu = 0;
167 c->lpv = 0;
168 c->effort = 0.0;
169 c->id = id;
170 cand_update_check_prime(c);
171 }
172
cand_set_presieved_values(cand_t c,const mpz_t u0,const mpz_t v0,const mpz_t u,const mpz_t v,unsigned int lpu,unsigned int lpv,unsigned long id)173 void cand_set_presieved_values(cand_t c, const mpz_t u0, const mpz_t v0,
174 const mpz_t u, const mpz_t v,
175 unsigned int lpu, unsigned int lpv,
176 unsigned long id) {
177 mpz_set(c->u0, u0);
178 mpz_set(c->v0, v0);
179 mpz_set(c->u, u);
180 mpz_set(c->v, v);
181 c->lpu = lpu;
182 c->lpv = lpv;
183 c->effort = 0.0;
184 c->id = id;
185 cand_update_check_prime(c);
186 }
187
188 // Cost is the bitsize of the largest unfactored part
189 // For factored numbers, this is INT_MAX
cost(const cand_t c)190 int cost(const cand_t c) {
191 if (cand_is_factored(c))
192 return INT_MAX;
193 else
194 return MAX(mpz_sizeinbase(c->u, 2), mpz_sizeinbase(c->v, 2));
195 }
196
197 // For use in qsort
cand_cmp(const void * c1,const void * c2)198 int cand_cmp(const void *c1, const void *c2) {
199 cand_t *C1 = (cand_t *)c1;
200 cand_t *C2 = (cand_t *)c2;
201 if (cost(*C1) > cost(*C2))
202 return 1;
203 if (cost(*C1) < cost(*C2))
204 return -1;
205 return 0;
206 }
207
208
209 ///////////////////////////////////////////////////////////////////////////
210 // Pool: contains a list of candidates that are still interesting to
211 // try to factor
212 /////////////////////////////////////////////////////////////////////////////
213
214 // The list of candidates must always be sorted by increasing cost.
215 typedef struct {
216 cand_t *list;
217 unsigned int n;
218 } pool_s;
219 typedef pool_s pool_t[1];
220
pool_init(pool_t p)221 void pool_init(pool_t p) {
222 p->list = NULL;
223 p->n = 0;
224 }
225
pool_clear(pool_t p)226 void pool_clear(pool_t p) {
227 for(unsigned int i = 0; i < p->n; i++)
228 cand_clear(p->list[i]);
229 free(p->list);
230 }
231
232 // if pool has been modified, sort it again
pool_sort(pool_t p)233 void pool_sort(pool_t p) {
234 qsort(p->list, p->n, sizeof(cand_t), cand_cmp);
235 }
236
237 // Insert candidate in pool, keeping it sorted according to cost
pool_insert(pool_t p,const cand_t c)238 void pool_insert(pool_t p, const cand_t c) {
239 p->list = (cand_t *)realloc(p->list, (p->n + 1)*sizeof(cand_t));
240 cand_init_copy(p->list[p->n], c);
241 for(unsigned int i = p->n;
242 i > 0 && cost(p->list[i-1]) > cost(p->list[i]);
243 i--) {
244 cand_swap(p->list[i-1], p->list[i]);
245 }
246 p->n += 1;
247 }
248
249 // Remove candidates that have a prime more than lmax, and those that
250 // are fully factored. If there are still more than max_size elements,
251 // keep those with smallest cost.
pool_purge(pool_t p,unsigned int max_size,unsigned int lmax)252 void pool_purge(pool_t p, unsigned int max_size, unsigned int lmax) {
253 unsigned int i, j;
254 j = 0;
255 for (i = 0; i < p->n; ++i) {
256 if (!cand_is_factored(p->list[i]) &&
257 !cand_is_probably_not_smooth(p->list[i], lmax) &&
258 (MAX(p->list[i]->lpu, p->list[i]->lpv) <= lmax)) {
259 cand_swap(p->list[j], p->list[i]);
260 j++;
261 }
262 }
263 unsigned int newsize = MIN(j, max_size);
264 for (i = newsize; i < p->n; ++i)
265 cand_clear(p->list[i]);
266 p->n = newsize;
267 p->list = (cand_t *)realloc(p->list, (p->n)*sizeof(cand_t));
268 }
269
pool_print(const pool_t p)270 void pool_print(const pool_t p) {
271 for (unsigned int i = 0; i < p->n; ++i) {
272 printf("%u: %lu %lu (%u)\n", i,
273 mpz_sizeinbase(p->list[i]->u, 2),
274 mpz_sizeinbase(p->list[i]->v, 2),
275 MAX(p->list[i]->lpu, p->list[i]->lpv));
276 }
277 }
278
279 ///////////////////////////////////////////////////////////////////////////
280 // Stats: keep track of statistics about the exepected number of bits
281 // obtained after running x curves.
282 ///////////////////////////////////////////////////////////////////////////////
283
284 #define MAX_CPT 2048
285 typedef struct {
286 double aver_gain[MAX_CPT]; // average number of bits removed after i curves
287 unsigned long nb_test[MAX_CPT]; // size of the sample on which this avearge was done
288 } stats_s;
289 typedef stats_s stats_t[1];
290
stats_init(stats_t S)291 void stats_init(stats_t S) {
292 for(unsigned int i = 0; i < MAX_CPT; ++i) {
293 S->aver_gain[i] = 0.0;
294 S->nb_test[i] = 0;
295 }
296 }
297
298 #if 0
299 void stats_clear(stats_t S) { }
300 #endif
301
stats_update(stats_t S,double gain,unsigned int i)302 void stats_update(stats_t S, double gain, unsigned int i) {
303 if (i >= MAX_CPT) {
304 abort();
305 }
306 double newav = S->aver_gain[i]*S->nb_test[i] + gain;
307 S->nb_test[i]++;
308 newav /= S->nb_test[i];
309 S->aver_gain[i] = newav;
310 }
311
stats_is_below_average(stats_t S,double gain,unsigned int i)312 int stats_is_below_average(stats_t S, double gain, unsigned int i) {
313 if (S->nb_test[i] <= 100)
314 return 0; // not enough data to conclude
315 else
316 return gain < EA_THRESHOLD*S->aver_gain[i];
317 }
318
stats_print(stats_t S)319 void stats_print(stats_t S) {
320 printf("[");
321 for (int i = 1; i < 500; ++i) {
322 if (S->nb_test[i] <= 100)
323 break;
324 printf(" %.0f", S->aver_gain[i]);
325 }
326 printf(" ]\n");
327 }
328
329 ///////////////////////////////////////////////////////////////////////////
330 // General data structure for an instance of the search of smooth numbers
331 ///////////////////////////////////////////////////////////////////////////
332
333 typedef struct {
334 pool_t pool;
335 stats_t stats;
336 void *param_next_cand;
337 int (*next_cand)(cand_t, void *); // new candidate put in first arg.
338 unsigned long target; // smoothness bound (in bits)
339 double current_effort; // current effort per candidate.
340 double max_effort;
341 unsigned long max_pool_size;
342 double minB1;
343 } context_s;
344 typedef context_s context_t[1];
345
remove_small_factors(mpz_t z)346 double remove_small_factors(mpz_t z) {
347 double gain = 0.0;
348 while(mpz_divisible_ui_p(z, 2)) {
349 mpz_divexact_ui(z, z, 2);
350 gain += 1.0;
351 }
352 while(mpz_divisible_ui_p(z, 3)) {
353 mpz_divexact_ui(z, z, 3);
354 gain += log2(3.0);
355 }
356 return gain;
357 }
358
359 // get a B1, so that we can quickly cover the target effort
get_B1_from_effort(double effort,double minB1)360 double get_B1_from_effort(double effort, double minB1) {
361 double B1 = minB1;
362 double S = B1;
363 while (S < effort) {
364 B1 += sqrt(B1);
365 S += B1;
366 }
367 return B1;
368 }
369
increase_effort(context_t ctx)370 void increase_effort(context_t ctx) {
371 ctx->current_effort += sqrt(ctx->current_effort);
372 ctx->current_effort = MIN(ctx->current_effort, ctx->max_effort);
373 }
374
my_ecm_factor(mpz_t f,mpz_t z,double B1)375 void my_ecm_factor(mpz_t f, mpz_t z, double B1) {
376 ecm_params ecm_par;
377 ecm_init(ecm_par);
378 long sig = random();
379 mpz_set_ui(ecm_par->sigma, sig);
380 ecm_par->B1done = default_B1done; /* issue with ECM 6.4.x */
381 ecm_factor(f, z, B1, ecm_par);
382 ecm_clear(ecm_par);
383 }
384
385 // One step of smoothness detection: get a new candidate, run a bunch of
386 // ECMs, update the pool, and the stats.
387 // Return value:
388 // 1: smooth
389 // 0: non-smooth
390 // -1: early stop, no more candidate to test.
smooth_detect_one_step(cand_t winner,context_t ctx)391 int smooth_detect_one_step(cand_t winner, context_t ctx) {
392 cand_t C;
393 cand_init(C);
394 // Get a new candidate, not obviously non-smooth
395 double gain_u = mpz_sizeinbase(C->u0, 2) - mpz_sizeinbase(C->u, 2);
396 double gain_v = mpz_sizeinbase(C->v0, 2) - mpz_sizeinbase(C->v, 2);
397 do {
398 int ret = ctx->next_cand(C, ctx->param_next_cand);
399 if (ret == 0) {
400 cand_clear(C);
401 return -1; // early stop, no more candidate.
402 }
403 gain_u += remove_small_factors(C->u);
404 gain_v += remove_small_factors(C->v);
405 cand_update_check_prime(C);
406 } while (C->lpu >= ctx->target || C->lpv >= ctx->target);
407
408 // Start a loop of ECM
409 double B1 = ctx->minB1; // initial B1
410 int cpt = 0; // number of curves tried on this number
411 mpz_t f; // output of ecm
412 mpz_init(f);
413 while (!cand_is_probably_not_smooth(C, ctx->target)
414 && !cand_is_factored(C)
415 && C->effort < ctx->current_effort) {
416 cpt++;
417 // u-side
418 if (mpz_cmp_ui(C->u, 1) != 0) {
419 my_ecm_factor(f, C->u, B1);
420 gain_u += log2(mpz_get_d(f));
421 stats_update(ctx->stats, gain_u, cpt);
422 if (mpz_cmp_ui(f, 1) > 0) {
423 mpz_divexact(C->u, C->u, f);
424 C->lpu = MAX(C->lpu, mpz_sizeinbase(f, 2));
425 cand_update_check_prime(C);
426 if (C->lpu >= ctx->target)
427 break;
428 }
429 }
430 // v-side
431 if (mpz_cmp_ui(C->v, 1) != 0) {
432 my_ecm_factor(f, C->v, B1);
433 gain_v += log2(mpz_get_d(f));
434 stats_update(ctx->stats, gain_v, cpt);
435 if (mpz_cmp_ui(f, 1) > 0) {
436 mpz_divexact(C->v, C->v, f);
437 C->lpv = MAX(C->lpv, mpz_sizeinbase(f, 2));
438 cand_update_check_prime(C);
439 if (C->lpv >= ctx->target)
440 break;
441 }
442 }
443 // if both gain are below average, then abort this candidate
444 if (stats_is_below_average(ctx->stats, gain_u, cpt) &&
445 stats_is_below_average(ctx->stats, gain_v, cpt)) {
446 // mark it as unsmooth, by increasing lpu
447 C->lpu = UINT_MAX;
448 break;
449 }
450 // remember current effort for this number, and increase B1.
451 C->effort += B1;
452 B1 += sqrt(B1);
453 }
454 mpz_clear(f);
455
456 // If we had a prime larger than expected, then abort.
457 unsigned int l = MAX(C->lpu, C->lpv);
458 if (l > ctx->target) {
459 cand_clear(C);
460 increase_effort(ctx);
461 return 0;
462 }
463
464 // Did we factor the candidate completely?
465 // If so, print result, otherwise, insert in pool
466 if (cand_is_factored(C)) {
467 cand_print(C);
468 cand_set(winner, C);
469 cand_clear(C);
470 return 1;
471 } else {
472 pool_insert(ctx->pool, C);
473 }
474 cand_clear(C);
475
476 B1 = get_B1_from_effort(ctx->current_effort, ctx->minB1);
477
478 // Run ECM on the numbers in pool, so that they all have received more
479 // or less the same effort.
480 mpz_init(f);
481 for (unsigned int i = 0; i < ctx->pool->n; ++i) {
482 cand_s * c = &ctx->pool->list[i][0];
483 double effort = ctx->current_effort;
484 // more effort for the most promising candidates!
485 // the correcting factors below are completely heuristic...
486 if (ctx->pool->n > 5) {
487 if (i == 0) { effort *= 2; }
488 if (i == 1) { effort *= 1.6; }
489 if (i == 2) { effort *= 1.3; }
490 if (i == 3) { effort *= 1.1; }
491 }
492 while (!cand_is_factored(c) && c->effort < effort) {
493 if (mpz_cmp_ui(c->u, 1) != 0) {
494 my_ecm_factor(f, c->u, B1);
495 if (mpz_cmp_ui(f, 1) > 0) {
496 mpz_divexact(c->u, c->u, f);
497 c->lpu = MAX(c->lpu, mpz_sizeinbase(f, 2));
498 }
499 }
500 if (mpz_cmp_ui(c->v, 1) != 0) {
501 my_ecm_factor(f, c->v, B1);
502 if (mpz_cmp_ui(f, 1) > 0) {
503 mpz_divexact(c->v, c->v, f);
504 c->lpv = MAX(c->lpv, mpz_sizeinbase(f, 2));
505 }
506 }
507 cand_update_check_prime(c);
508 c->effort += B1;
509 if (MAX(c->lpu, c->lpv) > ctx->target)
510 break;
511 }
512 }
513 mpz_clear(f);
514
515 // Go through pool: detect winners, purge bad candidates, keep best.
516 for (unsigned int i = 0; i < ctx->pool->n; ++i) {
517 cand_s * c = &ctx->pool->list[i][0];
518 unsigned int l = MAX(c->lpu, c->lpv);
519 if (cand_is_factored(c) && l <= ctx->target) {
520 cand_print(c);
521 cand_set(winner, c);
522 return 1;
523 }
524 }
525 pool_sort(ctx->pool);
526 pool_purge(ctx->pool, ctx->max_pool_size, ctx->target);
527
528 increase_effort(ctx);
529 return 0;
530 }
531
smooth_detect(cand_t C,int (* next_cand)(cand_t,void *),void * param_next_cand,unsigned long target,const smooth_detect_param_s * param)532 void smooth_detect(cand_t C, int (*next_cand)(cand_t, void *),
533 void *param_next_cand, unsigned long target,
534 const smooth_detect_param_s* param) {
535 /* fix issue with ECM 6.4.x */
536 {
537 ecm_params params;
538 ecm_init(params);
539 default_B1done = params->B1done;
540 ecm_clear(params);
541 }
542
543 const smooth_detect_param_s default_param = {2000.0, 1e20, 10, 1, 100.0};
544 if (param == NULL) {
545 param = &default_param;
546 }
547
548 // Create a context.
549 context_t ctx;
550 pool_init(ctx->pool);
551 stats_init(ctx->stats);
552 ctx->param_next_cand = param_next_cand;
553 ctx->next_cand = next_cand;
554 ctx->target = target;
555 ctx->current_effort = param->min_effort;
556 ctx->max_effort = param->max_effort;
557 ctx->max_pool_size = param->max_pool_size;
558 ctx->minB1 = param->minB1;
559
560 double tm = get_time();
561 int cpt = 0;
562 int found = 0;
563 while (found == 0) {
564 found = smooth_detect_one_step(C, ctx);
565 cpt++;
566 if (param->verbose && (cpt % 20 == 0)) {
567 printf("***** Pool status after %d candidates in %.1fs\n", cpt,
568 get_time()-tm);
569 printf("current_effort = %.0f\n", ctx->current_effort);
570 printf("current max B1 = %.0f\n",
571 get_B1_from_effort(ctx->current_effort, ctx->minB1));
572 printf("current stats:\n");
573 stats_print(ctx->stats);
574 pool_print(ctx->pool);
575 }
576 }
577 }
578
579