1710.sagews

Pollard $p-1$

%auto
def lcm_upto(B):
    return prod([p^int(math.log(B)/math.log(p)) for p in prime_range(B+1)])

lcm_upto(10^2)

69720375229712477164533808935312303556800

LCM([1..10^2])

69720375229712477164533808935312303556800

time v = lcm_upto(10^5)

Time: CPU 0.06 s, Wall: 0.06 s

time w = LCM([1..10^5])

Time: CPU 0.53 s, Wall: 1.33 s

v == w

True

len(str(v))

43452

%auto
def pollard(N, B=10^5, stop=10):
    m = prod([p^int(math.log(B)/math.log(p)) 
             for p in prime_range(B+1)])
    for a in [2..stop]:
        x = (Mod(a,N)^m - 1).lift()
        if x == 0: continue
        g = gcd(x, N)
        if g != 1 or g != N: return g
    return 1

pollard(779167, B = 5)

1

pollard(779167, B=15)

2003

It does not matter if $N$ is big.  All that matters is that $N$ has some factor $p$ with $p-1$ smooth.  Pretty amazing, eh?

N = 2003*next_prime(ZZ.random_element(10^100)); N

12200812832238143616017533240222183238927105188567926260963719543150415852204779120537357146219613867529

pollard(N, B=15)

2003

Elliptic Curve Factorization Method

def ecm(N, B=10^3, trials=10):
    m = lcm_upto(B)
    R = Integers(N)
    # Sneaky: make Sage think that R is a field:
    def f(): return True
    R.is_field = f
    for i in range(trials):
        while True:
            a = R.random_element()
            if gcd(4*a.lift()^3 + 27, N) == 1: break
        try:
            m * EllipticCurve([a, 1])([0,1])
        except ZeroDivisionError, msg:
            print msg
            # msg: "Inverse of <int> does not exist"
            return gcd(Integer(str(msg).split()[2]), N)
    return 1

N = 2003*next_prime(10^20); N

200300000000000000078117

time ecm(N, B=20)

Inverse of 65550682212630542772312 does not exist 2003 Time: CPU 0.03 s, Wall: 0.09 s

factor(2002)

2 * 7 * 11 * 13

Next, find a factor $p$ where $p-1$ is not smooth (so Pollard rho doesn't work).

factor(389-1)

2^2 * 97

N = 389 * next_prime(10^30); N

389000000000000000000000000022173

pollard(N, B=90)

1

ecm(N,B=90)

389

pollard(N, B=100)

389

for a in [-7..7]:
    print factor(389-a)

2^2 * 3^2 * 11 5 * 79 2 * 197 3 * 131 2^3 * 7^2 17 * 23 2 * 3 * 5 * 13 389 2^2 * 97 3^2 * 43 2 * 193 5 * 7 * 11 2^7 * 3 383 2 * 191

N = next_prime(10^18) * next_prime(10^30); N

1000000000000000003000000000057000000000000000171

time ecm(N,B=10^3)

1 Time: CPU 1.63 s, Wall: 1.68 s

time ecm(N,B=10^4)

1 Time: CPU 10.83 s, Wall: 17.63 s

time ecm(N,B=10^4)

1 Time: CPU 10.45 s, Wall: 18.15 s

Sage includes GMP-ECM, which is the world's best public implementation of the ECM algorithm.

See http://www.loria.fr/~zimmerma/records/ecmnet.html

A = sage.interfaces.ecm.ECM(1, watch=True)

A.factor(N)

[1000000000000000003, 1000000000000000000000000000057]

sage.interfaces.ecm.ECM?

File: /home/wstein/sage/local/lib/python2.6/site-packages/sage/interfaces/ecm.py

Type: <type 'classobj'>

Definition: sage.interfaces.ecm.ECM(n, watch=False)

Docstring:

        Create an interface to the GMP-ECM elliptic curve method
        factorization program.

        See http://ecm.gforge.inria.fr/.

        AUTHORS:
            These people wrote GMP-ECM:
              Pierrick Gaudry, Jim Fougeron,
              Laurent Fousse, Alexander Kruppa,
              Dave Newman, Paul Zimmermann
            
        William Stein and Robert Bradshaw -- wrote the Sage interface to GMP-ECM

        INPUT:
            B1 -- stage 1 bound
            B2 -- stage 2 bound (or interval B2min-B2max)

            x0 -- x        use x as initial point
            sigma -- s     use s as curve generator [ecm]
            A -- a         use a as curve parameter [ecm]
            k -- n         perform >= n steps in stage 2
            power -- n     use x^n for Brent-Suyama's extension
            dickson -- n   use n-th Dickson's polynomial for Brent-Suyama's extension
            c -- n         perform n runs for each input
            pm1 --         perform P-1 instead of ECM
            pp1 --         perform P+1 instead of ECM
            q --           quiet mode
            v --           verbose mode
            timestamp --   print a time stamp with each number
            mpzmod --      use GMP's mpz_mod for mod reduction
            modmuln --     use Montgomery's MODMULN for mod reduction
            redc --        use Montgomery's REDC for mod reduction
            nobase2 --     disable special base-2 code
            base2 -- n     force base 2 mode with 2^n+1 (n>0) or 2^n-1 (n<0)
            save -- file   save residues at end of stage 1 to file
            savea -- file  like -save, appends to existing files
            resume -- file resume residues from file, reads from
                           stdin if file is "-"
            primetest --   perform a primality test on input
            treefile -- f  store product tree of F in files f.0 f.1 ...
            i -- n         increment B1 by this constant on each run
            I -- f         auto-calculated increment for B1 multiplied by 'f' scale factor
            inp -- file    Use file as input (instead of redirecting stdin)
            b --           Use breadth-first mode of file processing
            d --           Use depth-first mode of file processing (default)
            one --         Stop processing a candidate if a factor is found (looping mode)
            n --           run ecm in 'nice' mode (below normal priority)
            nn --          run ecm in 'very nice' mode (idle priority)
            t -- n         Trial divide candidates before P-1, P+1 or ECM up to n
            ve -- n        Verbosely show short (< n character) expressions on each loop
            cofdec --      Force cofactor output in decimal (even if expressions are used)
            B2scale -- f   Multiplies the default B2 value by f
            go -- val      Preload with group order val, which can be a simple expression,
                           or can use N as a placeholder for the number being factored.
            prp -- cmd     use shell command cmd to do large primality tests
            prplen -- n    only candidates longer than this number of digits are 'large'
            prpval -- n    value>=0 which indicates the prp command foundnumber to be PRP.
            prptmp -- file outputs n value to temp file prior to running (NB. gets deleted)
            prplog -- file otherwise get PRP results from this file (NB. gets deleted)
            prpyes -- str  literal string found in prplog file when number is PRP
            prpno -- str   literal string found in prplog file when number is composite