source: trunk/CrypPlugins/QuadraticSieve/msieve/common/driver.c @ 1373

Last change on this file since 1373 was 1373, checked in by Sven Rech, 12 years ago

some msieve modifications

File size: 15.0 KB
Line 
1/*--------------------------------------------------------------------
2This source distribution is placed in the public domain by its author,
3Jason Papadopoulos. You may use it for any purpose, free of charge,
4without having to notify anyone. I disclaim any responsibility for any
5errors.
6
7Optionally, please be nice and tell me if you find this source to be
8useful. Again optionally, if you add to the functionality present here
9please consider making those additions public too, so that others may
10benefit from your work.
11
12$Id: driver.c 23 2009-07-20 02:59:07Z jasonp_sf $
13--------------------------------------------------------------------*/
14
15#include <common.h>
16#include <wrapper.h>
17
18/*--------------------------------------------------------------------*/
19msieve_obj * msieve_obj_new(char *input_integer, uint32 flags,
20                            char *savefile_name, char *logfile_name,
21                            char *nfs_fbfile_name,
22                            uint32 seed1, uint32 seed2, uint32 max_relations,
23                            uint64 nfs_lower, uint64 nfs_upper,
24                            enum cpu_type cpu,
25                            uint32 cache_size1, uint32 cache_size2,
26                            uint32 num_threads, uint32 mem_mb) {
27
28        msieve_obj *obj = (msieve_obj *)xcalloc((size_t)1, sizeof(msieve_obj));
29
30        if (obj == NULL)
31                return obj;
32       
33        obj->input = input_integer;
34        obj->flags = flags;
35        obj->seed1 = seed1;
36        obj->seed2 = seed2;
37        obj->max_relations = max_relations;
38        obj->nfs_lower = nfs_lower;
39        obj->nfs_upper = nfs_upper;
40        obj->cpu = cpu;
41        obj->cache_size1 = cache_size1;
42        obj->cache_size2 = cache_size2;
43        obj->num_threads = num_threads;
44        obj->mem_mb = mem_mb;
45        obj->logfile_name = MSIEVE_DEFAULT_LOGFILE;
46        if (logfile_name)
47                obj->logfile_name = logfile_name;
48        obj->nfs_fbfile_name = MSIEVE_DEFAULT_NFS_FBFILE;
49        if (nfs_fbfile_name)
50                obj->nfs_fbfile_name = nfs_fbfile_name;
51        savefile_init(&obj->savefile, savefile_name);
52       
53        return obj;
54}
55
56/*--------------------------------------------------------------------*/
57msieve_obj * msieve_obj_free(msieve_obj *obj) {
58
59        msieve_factor *curr_factor;
60
61        curr_factor = obj->factors;
62        while (curr_factor != NULL) {
63                msieve_factor *next_factor = curr_factor->next;
64                free(curr_factor->number);
65                free(curr_factor);
66                curr_factor = next_factor;
67        }
68
69        savefile_free(&obj->savefile);
70        free(obj);
71        return NULL;
72}
73
74/*--------------------------------------------------------------------*/
75static uint32 msieve_run_core(msieve_obj *obj, mp_t *n, 
76                                factor_list_t *factor_list) {
77
78        uint32 i;
79        uint32 bits;
80
81        /* detect if n is a perfect power. Try extracting any root
82           whose value would exceed the trial factoring bound */
83
84        i = 2;
85        while (1) {
86                mp_t n2;
87                if (mp_iroot(n, i, &n2) == 0) {
88                        factor_list_add(obj, factor_list, &n2);
89                        return 1;
90                }
91                if (n2.nwords == 1 && n2.val[0] < PRECOMPUTED_PRIME_BOUND)
92                        break;
93                i++;
94        }
95
96        /* If n is small enough, use custom routines */
97
98        bits = mp_bits(n);
99        if (bits <= SMALL_COMPOSITE_CUTOFF_BITS) {
100                mp_t n1, n2;
101                mp_clear(&n1);
102                mp_clear(&n2);
103
104                if (bits <= 60) {
105                        i = squfof(n);
106                        if (i > 1) {
107                                n1.nwords = 1;
108                                n1.val[0] = i;
109                                mp_divrem_1(n, i, &n2);
110                        }
111                        else {
112                                tinyqs(n, &n1, &n2);
113                        }
114                }
115                else {
116                        tinyqs(n, &n1, &n2);
117                }
118
119                if (!mp_is_zero(&n1) && !mp_is_zero(&n2) &&
120                    !mp_is_one(&n1) && !mp_is_one(&n2)) {
121                        factor_list_add(obj, factor_list, &n1);
122                        factor_list_add(obj, factor_list, &n2);
123                        return 1;
124                }
125                else {
126                        /* if the tiny factoring routines failed, then
127                           the heavy artillery routines would also fail
128                           (the input is too small). Just give up */
129
130                        printf("error: tiny factoring failed\n");
131                        return 0;
132                }
133        }
134
135        /* Beyond this point we use the heavy artillery. */
136        //return factor_mpqs(obj, n, factor_list);
137        get_trivial_factorlist(factor_list, obj);
138}
139
140/*--------------------------------------------------------------------*/
141void msieve_run(msieve_obj *obj) {
142
143        char *n_string;
144        int32 status;
145        uint32 i;
146        mp_t n, reduced_n;
147        factor_list_t factor_list;
148        time_t start_time;
149
150        /* convert the input number to an mp_t */
151
152        status = evaluate_expression(obj->input, &n);
153        if (status < 0 || mp_is_zero(&n)) {
154                printf("error %d converting '%s'\n", status, obj->input);
155                obj->flags |= MSIEVE_FLAG_FACTORIZATION_DONE;
156                return;
157        }
158        n_string = mp_sprintf(&n, 10, obj->mp_sprintf_buf);
159
160        /* print startup banner */
161
162        logprintf(obj, "\n");
163        logprintf(obj, "\n");
164        logprintf(obj, "Msieve v. %d.%02d\n", MSIEVE_MAJOR_VERSION, 
165                                        MSIEVE_MINOR_VERSION);
166        start_time = time(NULL);
167        if (obj->flags & MSIEVE_FLAG_LOG_TO_STDOUT) {
168                printf("%s", ctime(&start_time));
169        }
170
171        logprintf(obj, "random seeds: %08x %08x\n", obj->seed1, obj->seed2);
172        logprintf(obj, "factoring %s (%d digits)\n", 
173                                n_string, strlen(n_string));
174
175        /* handle trivial inputs */
176
177        if (mp_is_zero(&n) || mp_is_one(&n)) {
178                add_next_factor(obj, &n, MSIEVE_PRIME);
179                obj->flags |= MSIEVE_FLAG_FACTORIZATION_DONE;
180                return;
181        }
182
183        /* perform trial division */
184
185        factor_list_init(&factor_list);
186        trial_factor(obj, &n, &reduced_n, &factor_list);
187        if (mp_is_one(&reduced_n))
188                goto clean_up;
189
190        /* save the remaining cofactor of n; if composite,
191           run Pollard Rho unconditionally */
192
193        if (factor_list_add(obj, &factor_list, &reduced_n) > 0)
194                rho(obj, &reduced_n, &reduced_n, &factor_list);
195
196        /* if still composite, and not quite small,
197           run P+-1 and ECM. These could take quite a while
198           to run, so allow interruptions to exit gracefully */
199
200        if (factor_list_max_composite(&factor_list) > 
201                                SMALL_COMPOSITE_CUTOFF_BITS) {
202
203                obj->flags |= MSIEVE_FLAG_SIEVING_IN_PROGRESS;
204                ecm_pp1_pm1(obj, &reduced_n, &reduced_n, &factor_list);
205                obj->flags &= ~MSIEVE_FLAG_SIEVING_IN_PROGRESS;
206                if (obj->flags & MSIEVE_FLAG_STOP_SIEVING)
207                        goto clean_up;
208        }
209
210        /* while forward progress is still being made */
211
212        while (1) {
213                uint32 num_factors = factor_list.num_factors;
214                status = 0;
215
216                /* process the next composite factor of n. Only
217                   attempt one factorization at a time, since
218                   the underlying list of factors could change */
219
220                for (i = 0; i < num_factors; i++) {
221                        final_factor_t *f = factor_list.final_factors[i];
222
223                        if (f->type == MSIEVE_COMPOSITE) {
224                                mp_t new_n;
225                                mp_copy(&f->factor, &new_n);
226                                status = msieve_run_core(obj, &new_n,
227                                                        &factor_list);
228                                break;
229                        }
230                }
231                if (status == 0 || (obj->flags & MSIEVE_FLAG_STOP_SIEVING))
232                        break;
233        }
234
235clean_up:
236        factor_list_free(&n, &factor_list, obj);
237        if (!(obj->flags & MSIEVE_FLAG_STOP_SIEVING))
238                obj->flags |= MSIEVE_FLAG_FACTORIZATION_DONE;
239        i = (uint32)(time(NULL) - start_time);
240        logprintf(obj, "elapsed time %02d:%02d:%02d\n", i / 3600,
241                                (i % 3600) / 60, i % 60);
242}
243
244/*--------------------------------------------------------------------*/
245void add_next_factor(msieve_obj *obj, mp_t *n, 
246                        enum msieve_factor_type factor_type) {
247
248        msieve_factor *new_factor = (msieve_factor *)xmalloc(
249                                        sizeof(msieve_factor));
250        char *type_string;
251        char *tmp;
252        size_t len;
253
254        if (factor_type == MSIEVE_PRIME)
255                type_string = "p";
256        else if (factor_type == MSIEVE_COMPOSITE)
257                type_string = "c";
258        else
259                type_string = "prp";
260
261        /* Copy n. We could use strdup(), but that causes
262           warnings for gcc on AMD64 */
263
264        tmp = mp_sprintf(n, 10, obj->mp_sprintf_buf);
265        len = strlen(tmp) + 1;
266        new_factor->number = (char *)xmalloc((size_t)len);
267        memcpy(new_factor->number, tmp, len);
268
269        new_factor->factor_type = factor_type;
270        new_factor->next = NULL;
271
272        if (obj->factors != NULL) {
273                msieve_factor *curr_factor = obj->factors;
274                while (curr_factor->next != NULL)
275                        curr_factor = curr_factor->next;
276                curr_factor->next = new_factor;
277        }
278        else {
279                obj->factors = new_factor;
280        }
281
282        logprintf(obj, "%s%d factor: %s\n", type_string, 
283                                (int32)(len - 1),
284                                new_factor->number);
285}
286
287/*--------------------------------------------------------------------*/
288void logprintf(msieve_obj *obj, char *fmt, ...) {
289
290        va_list ap;
291
292        /* do *not* initialize 'ap' and use it twice; this
293           causes crashes on AMD64 */
294
295        if (obj->flags & MSIEVE_FLAG_USE_LOGFILE) {
296                time_t t = time(NULL);
297                char buf[64];
298                FILE *logfile = fopen(obj->logfile_name, "a");
299
300                if (logfile == NULL) {
301                        fprintf(stderr, "cannot open logfile\n");
302                        throwException("cannot open logfile\n");
303                }
304
305                va_start(ap, fmt);
306                buf[0] = 0;
307                strcpy(buf, ctime(&t));
308                *(strchr(buf, '\n')) = 0;
309                fprintf(logfile, "%s  ", buf);
310                vfprintf(logfile, fmt, ap);
311                fclose(logfile);
312                va_end(ap);
313        }
314        if (obj->flags & MSIEVE_FLAG_LOG_TO_STDOUT) {
315                va_start(ap, fmt);
316                vfprintf(stdout, fmt, ap);
317                va_end(ap);
318        }
319}
320
321/*--------------------------------------------------------------------*/
322void fill_prime_list(prime_list_t *prime_list,
323                        uint32 max_list_size,
324                        uint32 max_prime) {
325        uint32 i, j;
326        uint32 num_alloc;
327        prime_sieve_t prime_sieve;
328        uint32 *list;
329
330        num_alloc = 1500;
331        list = prime_list->list = (uint32 *)xmalloc(
332                                        num_alloc * sizeof(uint32));
333
334        init_prime_sieve(&prime_sieve, 0, max_prime);
335        for (i = j = 0; i < max_list_size; i++) {
336                uint32 prime = get_next_prime(&prime_sieve);
337
338                if (prime > max_prime)
339                        break;
340                if (j == num_alloc) {
341                        num_alloc *= 2;
342                        list = prime_list->list = (uint32 *)xrealloc(
343                                        prime_list->list,
344                                        num_alloc * sizeof(uint32));
345                }
346                list[j++] = prime;
347        }
348        prime_list->num_primes = j;
349        free_prime_sieve(&prime_sieve);
350}
351
352/*--------------------------------------------------------------------*/
353void factor_list_init(factor_list_t *list) {
354
355        memset(list, 0, sizeof(factor_list_t));
356}
357
358/*--------------------------------------------------------------------*/
359uint32 factor_list_max_composite(factor_list_t *list) {
360
361        uint32 i, bits;
362
363        /* Find the number of bits in the largest composite factor,
364           and return that (to give calling code an estimate of
365           how much work would be left if it stopped trying to
366           find new factors now) */
367
368        for (i = bits = 0; i < list->num_factors; i++) {
369                final_factor_t *curr_factor = list->final_factors[i];
370
371                if (curr_factor->type == MSIEVE_COMPOSITE) {
372                        uint32 curr_bits = mp_bits(&curr_factor->factor);
373                        bits = MAX(bits, curr_bits);
374                }
375        }
376
377        return bits;
378}
379
380/*--------------------------------------------------------------------*/
381static int compare_factors (const void *x, const void *y) {
382        final_factor_t **xx = (final_factor_t **)x;
383        final_factor_t **yy = (final_factor_t **)y;
384        return mp_cmp(&((*xx)->factor), &((*yy)->factor));
385}
386
387void factor_list_free(mp_t *n, factor_list_t *list, msieve_obj *obj) {
388
389        uint32 i;
390        mp_t q, r, tmp;
391
392        /* if there's only one factor and it's listed as composite,
393           don't save anything (it would just confuse people) */
394
395        if (list->num_factors == 1 && 
396            list->final_factors[0]->type == MSIEVE_COMPOSITE) {
397                free(list->final_factors[0]);
398                return;
399        }       
400
401        /* sort the factors in order of increasing size */
402
403        if (list->num_factors > 1)
404                qsort(list->final_factors, (size_t)list->num_factors, 
405                                sizeof(final_factor_t *), compare_factors);
406
407        /* report each factor every time it appears in n */
408
409        mp_copy(n, &tmp);
410        for (i = 0; i < list->num_factors; i++) {
411                final_factor_t *curr_factor = list->final_factors[i];
412
413                while (1) {
414                        mp_divrem(&tmp, &curr_factor->factor, &q, &r);
415                        if (mp_is_zero(&q) || !mp_is_zero(&r))
416                                break;
417
418                        mp_copy(&q, &tmp);
419                        add_next_factor(obj, &curr_factor->factor, 
420                                        curr_factor->type);
421                }
422
423                free(curr_factor);
424        }
425}
426
427/*--------------------------------------------------------------------*/
428static void factor_list_add_core(msieve_obj *obj, 
429                                factor_list_t *list, 
430                                mp_t *new_factor) {
431
432        /* recursive routine to do the actual adding of factors.
433           Upon exit, the factors in 'list' will be mutually coprime */
434
435        uint32 i;
436        mp_t tmp1, tmp2, common, q, r;
437        uint32 num_factors = list->num_factors;
438
439        mp_clear(&tmp1); tmp1.nwords = tmp1.val[0] = 1;
440        mp_clear(&tmp2); tmp2.nwords = tmp2.val[0] = 1;
441        mp_clear(&common);
442
443        /* compare new_factor to all the factors
444           already in the list */
445
446        for (i = 0; i < num_factors; i++) {
447                final_factor_t *curr_factor = list->final_factors[i];
448
449                /* skip new_factor if element i of the current
450                   list would duplicate it */
451
452                if (mp_cmp(new_factor, &curr_factor->factor) == 0)
453                        return;
454
455                /* if new_factor has a factor C in common with element
456                   i in the list, remove all instances of C, remove
457                   factor i, compress the list and postprocess */
458
459                mp_gcd(new_factor, &curr_factor->factor, &common);
460                if (!mp_is_one(&common)) {
461
462                        mp_copy(new_factor, &tmp1);
463                        while (1) {
464                                mp_divrem(&tmp1, &common, &q, &r);
465                                if (mp_is_zero(&q) || !mp_is_zero(&r))
466                                        break;
467                                mp_copy(&q, &tmp1);
468                        }
469
470                        mp_copy(&curr_factor->factor, &tmp2);
471                        while (1) {
472                                mp_divrem(&tmp2, &common, &q, &r);
473                                if (mp_is_zero(&q) || !mp_is_zero(&r))
474                                        break;
475                                mp_copy(&q, &tmp2);
476                        }
477
478                        free(list->final_factors[i]);
479                        list->final_factors[i] = 
480                                list->final_factors[--list->num_factors];
481                        break;
482                }
483        }
484
485        if (i < num_factors) {
486
487                /* there is overlap between new_factor and one
488                   of the factors previously found. In the worst
489                   case there are three new factors to deal with */
490
491                if (!mp_is_one(&tmp1))
492                        factor_list_add_core(obj, list, &tmp1);
493                if (!mp_is_one(&tmp2))
494                        factor_list_add_core(obj, list, &tmp2);
495                if (!mp_is_one(&common))
496                        factor_list_add_core(obj, list, &common);
497        }
498        else {
499                /* list doesn't need to be modified, except
500                   to append new_factor. We go to some trouble to
501                   avoid unnecessary primality tests when new_factor
502                   is small */
503
504                i = list->num_factors++;
505                list->final_factors[i] = (final_factor_t *)xmalloc(
506                                                sizeof(final_factor_t));
507                if (new_factor->nwords <= 2 &&
508                    ((uint64)new_factor->val[1] << 32 | 
509                                        new_factor->val[0]) <
510                    ((uint64)PRECOMPUTED_PRIME_BOUND * 
511                                        PRECOMPUTED_PRIME_BOUND)) {
512                        list->final_factors[i]->type = MSIEVE_PRIME;
513                }
514                else {
515                        list->final_factors[i]->type = (mp_is_prime(
516                                                new_factor, 
517                                                &obj->seed1, &obj->seed2)) ?
518                                                MSIEVE_PROBABLE_PRIME : 
519                                                MSIEVE_COMPOSITE;
520                }
521                mp_copy(new_factor, &(list->final_factors[i]->factor));
522        }
523}
524
525/*--------------------------------------------------------------------*/
526uint32 factor_list_add(msieve_obj *obj, factor_list_t *list, 
527                                mp_t *new_factor) {
528
529        if (!mp_is_zero(new_factor) && !mp_is_one(new_factor))
530                factor_list_add_core(obj, list, new_factor);
531       
532        return factor_list_max_composite(list);
533}
534
535/*--------------------------------------------------------------------*/
536void free_cycle_list(la_col_t *cycle_list, uint32 num_cycles) {
537
538        uint32 i;
539
540        for (i = 0; i < num_cycles; i++)
541                free(cycle_list[i].cycle.list);
542        free(cycle_list);
543}
544
545/*------------------------------------------------------------------*/
546uint32 merge_relations(uint32 *merge_array,
547                  uint32 *src1, uint32 n1,
548                  uint32 *src2, uint32 n2) {
549
550        /* Given two sorted lists of integers, merge
551           the lists into a single sorted list. We assume
552           each list contains no duplicate entries.
553           If a particular entry occurs in both lists,
554           don't add it to the final list at all. Returns
555           the number of elements in the resulting list */
556
557        uint32 i1, i2;
558        uint32 num_merge;
559
560        i1 = i2 = 0;
561        num_merge = 0;
562
563        while (i1 < n1 && i2 < n2) {
564                uint32 val1 = src1[i1];
565                uint32 val2 = src2[i2];
566
567                if (val1 < val2) {
568                        merge_array[num_merge++] = val1;
569                        i1++;
570                }
571                else if (val1 > val2) {
572                        merge_array[num_merge++] = val2;
573                        i2++;
574                }
575                else {
576                        i1++;
577                        i2++;
578                }
579        }
580
581        while (i1 < n1)
582                merge_array[num_merge++] = src1[i1++];
583        while (i2 < n2)
584                merge_array[num_merge++] = src2[i2++];
585
586        return num_merge;
587}
588
Note: See TracBrowser for help on using the repository browser.