NTIC A Taste of Modernity

Section 12.6 A Taste of Modernity

Now, these methods are the beginnings of how people really factor big numbers. Typically, one does trial division up to a certain size (maybe the first few hundred or thousand primes), then perhaps some modification of Fermat to make sure that there aren't any factors close to the square root if you are attacking something like RSA where that would otherwise be advantageous.

🔗

Then what?

🔗

There are many answers to this question, some of which involve cool things called continued fractions or number fields. See Exercises 12.7.13–12.7.15 to investigate these, starting with a simpler (but related) algorithm to Lenstra's in Exercise 12.7.13. See [E.5.2, Chapter 9] for an elementary approach to Exercise 12.7.15.

🔗

Another important modern technique that is beginning to show up in introductory textbooks is Lenstra's elliptic curve method; see once again either [E.4.19] or [E.2.10, Chapter 18.7] for details at the level of this text.

🔗

We won't touch more on those topics, but Algorithm 12.4.7 brought up a concept important in factoring, not just finding primes. Namely, we could come up with some probabilistic/random methods. That's right, we are going to try to find a factor randomly!

🔗

Subsection 12.6.1 The Pollard Rho algorithm

🔗

Here is the essence of this random (or ‘Monte Carlo’) approach; it is highly recursive, like many good algorithms.

🔗

Algorithm 12.6.1. Generic routine for “random” factoring.

Follow these steps.

Pick some polynomial that will be easy to compute mod $(n) .$
Plug in an essentially random seed value. (Often the seed is $2$ or $3 .$ )
Compute the polynomial's value at the seed.
If that has a non-trivial gcd with $n,$ we have a factor. Otherwise, plug the new value back into the polynomial, and repeat (and hope it eventually succeeds).

🔗

Below is code for the method we'll discuss in this section. It has a modification to the generic algorithm which I will discuss below.

xxxxxxxxxx
 
def PollardRhoFactor(n,kstop=50,seed=2):
    d=1
    a,b=seed,seed
    k=1
    def f(x):
        return (x^2+1)%n
    while (d==1 or d==n):
        a = f(a)
        b = f(f(b))
        d=gcd(a-b,n)
        k=k+1
        if k>kstop:
            print("Pollard Rho breaking off after %s rounds"%k)
            break
    if d>1:
        print("Pollard Rho took %s rounds"%k)
        print("The number it tried in the last round was %s, which shares factor %s"%(a-b,d))
        print("And %s is a factor of %s since %s * %s=%s"%(d,n,d,n/d,d*(n/d)))

🔗

The essence of the method is that by plugging in the values of the polynomial modulo

n,

we are generating a ‘pseudo-random’ sequence of numbers.

$x_{0} \equiv 2$ (mod $n$ )
$x_{1} \equiv f (x_{0})$ (mod $n$ )
$x_{2} \equiv f (x_{1})$ (mod $n$ )
$x_{i + 1} \equiv f (x_{i}),$ all (mod $n$ ).

🔗

Such a ‘pseudo-random’ sequence might be better than the sequences we used for trial division or Fermat factorization, precisely because it will likely hit some small(ish) factors and some fairly large factors, a good mix. It might also be good that it could give us numbers which, although not a factor of

n,

might at least share a factor with

n .

🔗

A first choice of seed and polynomial might be

x_{0} = 2

and

f (x) = x^{2} + 1 .

These choices could be different, but they are typical; John Pollard's original paper used

f (x) = x^{2} - 1,

for instance.

🔗

Example 12.6.2.

Let's try computing what we get for some specific numbers. Picking $n = 8051$ as in [E.2.4, Example 3.25], we get results as in the following interact.

xxxxxxxxxx
 
var('x')
@interact
def _(seed=2,n=8051,poly=x^2+1,trials=(10,[2..50])):
    f(x)=poly
    for i in range(trials):
        pretty_print(html("$x_{%s}=%s$"%(i,seed)))
        seed = (ZZ(f(seed)) % ZZ(n))
    pretty_print(html("$x_{%s}=%s$"%(i+1,seed)))

Notice that for $n = 8051,$ the term $x_{4} \equiv x_{19}$ (mod $n$ ), so the sequence, while seeming fairly random, will not come up with every possible divisor. With seed $3,$ $x_{13} \equiv x_{28}$ is the first repeat; if instead we change the modulus to $8053,$ there is no repeat through $x_{50} .$ So you can see the output will be hard to predict.

🔗

Although the outputs of

x_{i + 1} = f (x_{i})

might already seem to be fairly random, we will not actually try to find common divisors of these numbers with

n .

Instead, we will try to see if all the differences

x_{i} - x_{j}

share a common factor with

n,

using the (highly efficient to compute) greatest common divisor. That gives a lot more opportunities to find a common factor than just comparing with

x_{i}!

And hopefully it's just as (or more) ‘random’, and just as effective at finding factors.

🔗

However, having all possible differences to check might slow things down too much. So instead there is a final modification to the algorithm.

🔗

First, since there are finitely many possible outcomes modulo any particular modulus

d,

the sequence of results will eventually repeat not just modulo

n,

but for any

d .

In particular, suppose

d

is a divisor of

n

such that

gcd (x_{i} - x_{j}, n) = d

for a specific pair

x_{i}

and

x_{j}

with

j > i .

🔗

Now consider the values of the sequence of

x_{ℓ}

modulo

d .

Because polynomials are well-defined in modular arithmetic, we have that

x_{i} \equiv x_{j} (mod d)

🔗

implies that, for any

m

we have

x_{m + i} \equiv f^{m} (x_{i}) \equiv f^{m} (x_{j}) \equiv x_{m + j} (mod d)

🔗

as well. When we let

m = j - i,

this becomes

x_{j} \equiv x_{2 j - i} (mod d)

🔗

which means the sequence (modulo

d,

the common divisor, not

n

) repeats itself every

j - i

terms (after

i,

of course, if

j - i < i

🔗

Now let

s

be an integer such that

s (j - i) \geq i .

Then

x_{s (j - i)}

appears after the periodic behavior (modulo

d

) begins, so

x_{s (j - i)} \equiv x_{s (j - i) + (j - i)} \equiv \dots \equiv x_{s (j - i) + s (j - i)} (mod d) .

🔗

If we now let

k = s (j - i)

this congruence means

x_{k} \equiv x_{2 k} (mod d)

🔗

d

is a divisor of

x_{2 k} - x_{k}

specifically, not just

x_{i} - x_{j} .

🔗

Finally, this means instead of checking all possible differences for a common divisor, we only have to check

gcd (x_{2 k} - x_{k}, n)

for all

k

in the algorithm.

🔗

Algorithm 12.6.3. Pollard Rho factoring algorithm.

Follow these steps.

Pick some polynomial $f (x)$ that will be easy to compute (mod $n$ ) (such as $x^{2} + 1,$ though other quadratics might be used).
Plug in an essentially random seed value $x_{0} .$ (Often the seed is $2$ or $3 .$ )
Compute the polynomial's value at the seed, $f (x_{0}) = x_{1} .$
Continue plugging in $f (x_{i}) = x_{i + 1},$ modulo $n .$
For each $k$ we check whether

$1 < gcd (x_{2 k} - x_{k}, n) = d < n .$

🔗

Since the algorithm doesn't always find a factor for any given combination of seed, number, and polynomial, there is nothing to prove per se. However, probabilistically (just like with Miller-Rabin) it should succeed for

k

in the neighborhood of the size of the square root of the smallest factor of

n .

(This is also the justification for introducing the algorithm in the original papers introducing this method and its variants.) So if

n

has a big, but not too big, divisor, this test should help us find that divisor.

🔗

Example 12.6.4.

Let's try this with $n = 9991 .$ Keeping $x^{2} + 1$ and seed $2,$ the numbers we would get for the first six rounds are

x_{0} = 2, x_{1} = 5, x_{2} = 26, x_{3} = 677, x_{4} = 8735, x_{5} = 8950

x_{6} = 4654, x_{7} = 9220, x_{8} = 4973, x_{9} = 3005, x_{10} = 8153 .

This gives differences as follows:

$x_{2} - x_{1} = 26 - 5 = 21$
$x_{4} - x_{2} = 8735 - 26 = 8709$
$x_{6} - x_{3} = 4654 - 677 = 3977$
$x_{8} - x_{4} = 4973 - 8735 = - 3762$
$x_{10} - x_{5} = 8153 - 8950 = - 797$

These factor as follows:

xxxxxxxxxx
 
factor(21), factor(8709), factor(3977), factor(3672), factor(797)

That is an impressive list of eight different prime factors that could potentially be shared with $9991$ in just five iterations. These differences have the following gcds with $9991 :$

xxxxxxxxxx
 
gcd(9991,21), gcd(9991,8709), gcd(9991,3977), gcd(9991,3672), gcd(9991,797)

Indeed the third one already caught a common divisor with $9991 .$

xxxxxxxxxx
 
PollardRhoFactor(8051)

🔗

Remark 12.6.5.

This method is called the Pollard rho method because (apparently) a very imaginative eye can interpret the $x_{i}$ eventually repeating (mod $d$ ) (in the example, $d = 97$ ) as a tail and then a loop, i.e. a Greek $ρ .$ John Pollard actually has another method named after him, the $p - 1$ method, which we will not explore; however, it is related to some of the more advanced methods we mentioned in the introduction to this section.

🔗

Example 12.6.6.

Sometimes the rho method doesn't come up with an answer quickly, or at all.

xxxxxxxxxx
 
PollardRhoFactor(991*997)

Here we took 50 rounds without success, using the seed $2 .$ Because of the $ρ$ repetition, it will never succeed. So what do you do then – bring out the really advanced methods? Not at all – just as with primality testing, you just change your starting point to try again!

xxxxxxxxxx
 
PollardRhoFactor(991*997,seed=3)

🔗

Subsection 12.6.2 More factorization

🔗

In general, there are other such probabilistic algorithms, and they are quite successful with factoring numbers which might have reasonably sized but not gigantic factors.

🔗

Historical remark 12.6.7. Factoring Fermat.

The big success of this algorithm was the 1980 factorization of $F_{8}$ (the eighth Fermat number) by Pollard and Richard Brent (see Brent's website). They used one of several variants of the method, due to Brent, to found the previously unknown prime factor⁶ $1238926361552897 .$ Finding this factor of $F_{8}$ took, in total, 2 hours on a UNIVAC 1100/42 (which for the time period was very fast, indeed).

If you want to memorize this historic number, Brent provides the phrase “I am now entirely persuaded to employ the method, a handy trick, on gigantic composite numbers” to help you remember it.

Interestingly, the other (much larger) factor was not proven to be prime until significantly later; and as of this writing (January 2021) even $F_{12}$ has not been fully factored! See Subsection 17.5.2 for even more information.

🔗

Things don't automatically work quickly even with today's far more powerful hardware.

xxxxxxxxxx
 
PollardRhoFactor(2^(2^8)+1,1000000) # one million rounds!

🔗

Hmm, what now? Let's change the seed.

xxxxxxxxxx
 
PollardRhoFactor(2^(2^8)+1,1000000,seed=3)

🔗

No one method will factor every number quickly. Luckily, we have bigger guns at our disposal in Sage (especially in the component program Pari), that polish thing off rather more quickly.

xxxxxxxxxx
 
factor(2^(2^8)+1)

🔗

That is a little better than two hours on a mainframe, or even on your computer, I hope you'll agree.

🔗

Real factorization algorithms use several different methods to attack different types of factors. We can try to simulate this in a basic way by creating a Sage interact. Evaluate the first cell to define things (don't forget to keep the rho method defined); then you can evaluate the second cell, which is the interact.

xxxxxxxxxx
 
def TrialDivFactor(n):
    p = next_prime(1)
    top = ceil(math.sqrt(n))
    while p < top:
        if mod(n,p)==0:
            break
        p=next_prime(p)
    if n==1:
        print("1 is not prime")
    elif p==n:
        print(n,"is prime")
    elif mod(n,p)==0:
        print(n,"factors as",p,"times",n/p)
    else:
        print(n,"is prime")
​
def FermatFactor(n,verbose=False):
    if n%2==0:
        raise TypeError("Input must be odd!")
    s=ceil(math.sqrt(n))
    top=(n+1)/2
    while is_square(s^2-n)==0:
        if verbose:
            print(s,"squared minus",n,"is",s^2-n,"which is not a perfect square")
        s=s+1
    t=sqrt(s^2-n)
    print("Fermat found that",s,"squared minus",t,"squared equals",n)
    if s^2==n:
        print("So",n,"was already a perfect square,",s,"times",s)
    elif s<top:
        print("So",s+t,"times",s-t,"equals",(s-t)*(s+t),"which is",n)
    elif s==top:
        print("So Fermat did not find a factor, which means",n,"is prime!")

xxxxxxxxxx
 
@interact
def _(n=991*997,method=['trial','Fermat','Pollard Rho']):
    if method=='trial':
        TrialDivFactor(n)
    if method=='Fermat':
        FermatFactor(n)
    if method=='Pollard Rho':
        PollardRhoFactor(n)

🔗

Sage note 12.6.8. Building interacts.

An interact is just a Sage/Python function, except with @interact before it. There are many different input widgets you can use; this one demonstrates using a list and an input box which takes any input. See the interact documentation or Quickstart for many examples and more details.

🔗

If you think this sort of thing is cool, the Cunningham Project is a place to explore. I particularly like their Most Wanted lists. The idea is this:

The Cunningham Project seeks to factor the numbers $b^{n} \pm 1$ for $b = 2, 3, 5, 6, 7, 10, 11, 12,$ up to high powers $n .$

🔗

Another interesting resource is Sage developer Paul Zimmermann's Integer Factoring Records. Finally, Wagstaff's The joy of factoring [E.4.12] has tons of awesome examples and procedures – far too many, really, as well as an excellent discussion of how to speed up trial division, etc.