Re: misc/43192: games/factor takes forever on some large inputs

[Matthias Drochner <M.Drochner%fz-juelich.de@localhost>]

The following reply was made to PR misc/43192; it has been noted by GNATS.

From: Matthias Drochner <M.Drochner%fz-juelich.de@localhost>
To: <dholland-bugs%netbsd.org@localhost>, <lhf%impa.br@localhost>
Cc: <gnats-bugs%netbsd.org@localhost>
Subject: Re: misc/43192: games/factor takes forever on some large inputs
Date: Fri, 23 Apr 2010 19:59:19 +0200

 --==_Exmh_6376895038440
 Content-Type: text/plain; charset="us-ascii"
 Content-Transfer-Encoding: quoted-printable
 
 
 David Holland wrote:
 > It looks to me like the code we have is a mixture of the Pollard p-1
 > algorithm and the Pollard rho algorithm
 
 For me it looks like p-1, only that the exponent ("i" in the code,
 "m" in the reference you cited) is not calculated by the Least-
 Common-Multiple algorithm but just tried from 2 with increments
 of 1. This way, useful exponents are hit eventually as well, but
 it is extremely wasteful. In particular, since prime factors
 up to ~65000 are already ruled out at this point, it doesn't make
 much sense to start with such small numbers. (Also, since we
 have a table of prime numbers in the code already, it would be
 easy to reuse it for the L.C.M. calculation.)
 Unfortunately, the L.C.M. exponent will grow very quickly for
 useful Power Smoothness values ("B" in the reference). For
 a B of 100000 which is in the order of magnitude of the prime
 factors already considered, log(m) is 100051.6. So it would
 be a number with tenthousends of digits.
 
 Luiz Henrique de Figueiredo wrote:
 > Pollard p-1 seems to have a hard time with
 > 1234567890123456789012345678901...
 
 Indeed. The factors are 7742394596501 and 159455563099482401.
 The p-1 values factor into
 2^2*5^3*15484789193 and 2^5*5^2*29*283*24286518079.
 If I understood that stuff correctly, this means that one
 would need an exponent calculated for a "B" of >1.5e10.
 That's certainly too much to be practical.
 There is of course the statistical possibility that a factor
 is found with smaller exponents, but I still got the impression
 that p-1 is not the right algorithm here.
 
 As a test, I've just modified the pollard_pminus1() function
 within factor.c to do the rho algorithm. The code isn't complete,
 it just uses 2^x+1 as a pseudorandom generator polynom, with a
 fixed start value of 1. It factors Luiz's number within 70s
 on a 2GHz CPU. I'll append the code, if you want to play with it.
 (just replace the function, should just plug in)
 
 best regards
 Matthias
 
 
 ---------------------------------------------------------------------------=
 ---------------------
 ---------------------------------------------------------------------------=
 ---------------------
 Forschungszentrum Juelich GmbH
 52425 Juelich
 Sitz der Gesellschaft: Juelich
 Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
 Vorsitzende des Aufsichtsrats: MinDir'in Baerbel Brumme-Bothe
 Geschaeftsfuehrung: Prof. Dr. Achim Bachem (Vorsitzender),
 Dr. Ulrich Krafft (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt,
 Prof. Dr. Sebastian M. Schmidt
 ---------------------------------------------------------------------------=
 ---------------------
 ---------------------------------------------------------------------------=
 ---------------------
 
 --==_Exmh_6376895038440
 Content-Type: text/plain; name="rho.c"; charset="us-ascii"
 Content-Description: rho.c
 Content-Disposition: attachment; filename="rho.c"
 
 static void
 pollard_pminus1(BIGNUM *val)
 {
        unsigned int steps_taken, steps_limit, a;
        BIGNUM *x, *y, *tmp, *diff, *num;
 
        x = BN_new();
        y = BN_new();
        tmp = BN_new();
        diff = BN_new();
        num = BN_new();
 
        if (!BN_mod_word(val, 2)) {
                printf(" 2");
                BN_div_word(val, 2);
                if (!BN_is_one(val))
                        pollard_pminus1(val);
                return;
        }
 
 startover:
        steps_taken = 0;
        steps_limit = 2;
        a = 1; /* XXX random */
        BN_set_word(x, 1); /* XXX random */
        BN_copy(y, x);
 
        for (;;) {
                BN_sqr(x, x, ctx);
                BN_add_word(x, a);
                BN_mod(tmp, x, val, ctx);
                BN_copy(x, tmp);
                if (BN_cmp(x, y) > 0)
                        BN_sub(diff, x, y);
                else
                        BN_sub(diff, y, x);
                if (BN_is_zero(diff)) {
                        printf("loop\n"); /* XXX retry with other a and x0 */
                        return;
                }
 
                BN_gcd(tmp, diff, val, ctx);
 
                if (!BN_is_one(tmp)) {
                        if (BN_is_prime(tmp, PRIME_CHECKS, NULL, NULL,
                            NULL) == 1) {
                                putchar(' ');
                                BN_print_dec_fp(stdout, tmp);
                        } else
                                pollard_pminus1(tmp);
                        fflush(stdout);
 
                        BN_div(num, NULL, val, tmp, ctx);
                        if (BN_is_one(num))
                                return;
                        if (BN_is_prime(num, PRIME_CHECKS, NULL, NULL,
                            NULL) == 1) {
                                putchar(' ');
                                BN_print_dec_fp(stdout, num);
                                fflush(stdout);
                                return;
                        }
                        BN_copy(val, num);
                        goto startover;
                }
 
                steps_taken++;
                if (steps_taken == steps_limit) {
                        BN_copy(y, x);
                        steps_taken = 0;
                        steps_limit *= 2;
                }
        }
 }
 
 --==_Exmh_6376895038440--