Nom

Number Notes — Some notes on the Number toolbox.

A difficult case for Fermat's method

This is a numerical experiment with a difficult number to factor with Fermat's method. Consider

n = 33554431 
= 1801 * 18631
= 1801 * ( 31  *  601 )
   

This number comes from

2^50 - 1 = 3.11.31.251.601.1801.4051
   

The algorithm wants to factor n = u * v with

u = 1801 < v = 18631
   

The algorithm begins with

a = sqrt(n) = 5793, b = a^2 - n = 67
   

and ends with

a=10216, b=8415
   

There are 10216 - 5793 = 4423 iterations required.

The following factor_fermatdo function applies the Fermat method recursively.

function result = factor_fermatdo ( n )
    a = ceil(sqrt(n))
    while ( %t )
        b2 = a * a - n
        // Computes if b is a square
        b = ceil(sqrt(b2))
        mprintf ( "a=%d, b=%d, n=%d, a^2 - b^2=%d, |n-(a^2-b^2)|=%d\n", ..
          a, b, n, a^2 - b^2 , abs(n-a^2+b^2))
        if ( b^2 == b2 ) then
            // b2 is a square so that b2 = b^2 and 
            // n = a^2 - b^2
            break
        end
        a = a + 1
    end
    // n = (a - b) * (a + b)
    p1 = a - b
    p2 = a + b
    isprime = ( ( p1==1 ) & ( p2==n ) )
    // Factor p1 and p2 recursively, if required
    result = []
    if ( isprime ) then
        // If n is prime, there is no need for recursive calls
        // since n = 1 x n (where the 1 factor is discarded)
        result ( 1 , $ + 1 ) = p2
    else
        // Factorize p1
        otherfactors1 = factor_fermatdo ( p1 )
        // Be p1 prime or composite, add its factors.
        result = [ result otherfactors1 ]
        // Factorize p2
        otherfactors2 = factor_fermatdo ( p2 )
        // Be p2 prime or composite, add its factors.
        result = [ result otherfactors2 ]
    end
endfunction
   

We first decompose 33554431 into primes, which requires a large number of iterations.

-->factor_fermatdo ( 33554431 )
a=5793, b=67, n=33554431, a^2 - b^2=33554360, |n-(a^2-b^2)|=71
a=5794, b=127, n=33554431, a^2 - b^2=33554307, |n-(a^2-b^2)|=124
a=5795, b=167, n=33554431, a^2 - b^2=33554136, |n-(a^2-b^2)|=295
a=5796, b=198, n=33554431, a^2 - b^2=33554412, |n-(a^2-b^2)|=19
a=5797, b=226, n=33554431, a^2 - b^2=33554133, |n-(a^2-b^2)|=298
a=5798, b=250, n=33554431, a^2 - b^2=33554304, |n-(a^2-b^2)|=127
[...]
a=6069, b=1811, n=33554431, a^2 - b^2=33553040, |n-(a^2-b^2)|=1391
a=6070, b=1814, n=33554431, a^2 - b^2=33554304, |n-(a^2-b^2)|=127
a=6071, b=1818, n=33554431, a^2 - b^2=33551917, |n-(a^2-b^2)|=2514
a=6072, b=1821, n=33554431, a^2 - b^2=33553143, |n-(a^2-b^2)|=1288
a=6073, b=1824, n=33554431, a^2 - b^2=33554353, |n-(a^2-b^2)|=78
a=6074, b=1828, n=33554431, a^2 - b^2=33551892, |n-(a^2-b^2)|=2539
[...]
[...]
a=7000, b=3931, n=33554431, a^2 - b^2=33547239, |n-(a^2-b^2)|=7192
a=7001, b=3932, n=33554431, a^2 - b^2=33553377, |n-(a^2-b^2)|=1054
a=7002, b=3934, n=33554431, a^2 - b^2=33551648, |n-(a^2-b^2)|=2783
a=7003, b=3936, n=33554431, a^2 - b^2=33549913, |n-(a^2-b^2)|=4518
a=7004, b=3938, n=33554431, a^2 - b^2=33548172, |n-(a^2-b^2)|=6259
a=7005, b=3939, n=33554431, a^2 - b^2=33554304, |n-(a^2-b^2)|=127
[...]
a=8000, b=5518, n=33554431, a^2 - b^2=33551676, |n-(a^2-b^2)|=2755
a=8001, b=5520, n=33554431, a^2 - b^2=33545601, |n-(a^2-b^2)|=8830
a=8002, b=5521, n=33554431, a^2 - b^2=33550563, |n-(a^2-b^2)|=3868
a=8003, b=5523, n=33554431, a^2 - b^2=33544480, |n-(a^2-b^2)|=9951
a=8004, b=5524, n=33554431, a^2 - b^2=33549440, |n-(a^2-b^2)|=4991
a=8005, b=5525, n=33554431, a^2 - b^2=33554400, |n-(a^2-b^2)|=31
[...]
a=9000, b=6889, n=33554431, a^2 - b^2=33541679, |n-(a^2-b^2)|=12752
a=9001, b=6890, n=33554431, a^2 - b^2=33545901, |n-(a^2-b^2)|=8530
a=9002, b=6891, n=33554431, a^2 - b^2=33550123, |n-(a^2-b^2)|=4308
a=9003, b=6892, n=33554431, a^2 - b^2=33554345, |n-(a^2-b^2)|=86
a=9004, b=6894, n=33554431, a^2 - b^2=33544780, |n-(a^2-b^2)|=9651
a=9005, b=6895, n=33554431, a^2 - b^2=33549000, |n-(a^2-b^2)|=5431
[...]
a=10000, b=8152, n=33554431, a^2 - b^2=33544896, |n-(a^2-b^2)|=9535
a=10001, b=8153, n=33554431, a^2 - b^2=33548592, |n-(a^2-b^2)|=5839
a=10002, b=8154, n=33554431, a^2 - b^2=33552288, |n-(a^2-b^2)|=2143
a=10003, b=8156, n=33554431, a^2 - b^2=33539673, |n-(a^2-b^2)|=14758
a=10004, b=8157, n=33554431, a^2 - b^2=33543367, |n-(a^2-b^2)|=11064
a=10005, b=8158, n=33554431, a^2 - b^2=33547061, |n-(a^2-b^2)|=7370
[...]
a=10210, b=8408, n=33554431, a^2 - b^2=33549636, |n-(a^2-b^2)|=4795
a=10211, b=8409, n=33554431, a^2 - b^2=33553240, |n-(a^2-b^2)|=1191
a=10212, b=8411, n=33554431, a^2 - b^2=33540023, |n-(a^2-b^2)|=14408
a=10213, b=8412, n=33554431, a^2 - b^2=33543625, |n-(a^2-b^2)|=10806
a=10214, b=8413, n=33554431, a^2 - b^2=33547227, |n-(a^2-b^2)|=7204
a=10215, b=8414, n=33554431, a^2 - b^2=33550829, |n-(a^2-b^2)|=3602
a=10216, b=8415, n=33554431, a^2 - b^2=33554431, |n-(a^2-b^2)|=0
   

This proves that 33554431 = p1 * p2, with:

p1 = a - b = 1801
p2 = a + b = 18631
   

This proves that 33554431 = 1801 * 18631. But 1801 is prime. We must now decompose 18631 into primes.

-->factor_fermatdo ( 18631 )
a=137, b=12, n=18631, a^2 - b^2=18625, |n-(a^2-b^2)|=6
a=138, b=21, n=18631, a^2 - b^2=18603, |n-(a^2-b^2)|=28
a=139, b=27, n=18631, a^2 - b^2=18592, |n-(a^2-b^2)|=39
a=140, b=32, n=18631, a^2 - b^2=18576, |n-(a^2-b^2)|=55
a=141, b=36, n=18631, a^2 - b^2=18585, |n-(a^2-b^2)|=46
a=142, b=40, n=18631, a^2 - b^2=18564, |n-(a^2-b^2)|=67
a=143, b=43, n=18631, a^2 - b^2=18600, |n-(a^2-b^2)|=31
a=144, b=46, n=18631, a^2 - b^2=18620, |n-(a^2-b^2)|=11
a=145, b=49, n=18631, a^2 - b^2=18624, |n-(a^2-b^2)|=7
a=146, b=52, n=18631, a^2 - b^2=18612, |n-(a^2-b^2)|=19
[...]
a=6, b=3, n=31, a^2 - b^2=27, |n-(a^2-b^2)|=4
a=7, b=5, n=31, a^2 - b^2=24, |n-(a^2-b^2)|=7
a=8, b=6, n=31, a^2 - b^2=28, |n-(a^2-b^2)|=3
a=9, b=8, n=31, a^2 - b^2=17, |n-(a^2-b^2)|=14
a=10, b=9, n=31, a^2 - b^2=19, |n-(a^2-b^2)|=12
a=11, b=10, n=31, a^2 - b^2=21, |n-(a^2-b^2)|=10
a=12, b=11, n=31, a^2 - b^2=23, |n-(a^2-b^2)|=8
a=13, b=12, n=31, a^2 - b^2=25, |n-(a^2-b^2)|=6
a=14, b=13, n=31, a^2 - b^2=27, |n-(a^2-b^2)|=4
a=15, b=14, n=31, a^2 - b^2=29, |n-(a^2-b^2)|=2
a=16, b=15, n=31, a^2 - b^2=31, |n-(a^2-b^2)|=0
a=25, b=5, n=601, a^2 - b^2=600, |n-(a^2-b^2)|=1
a=26, b=9, n=601, a^2 - b^2=595, |n-(a^2-b^2)|=6
a=27, b=12, n=601, a^2 - b^2=585, |n-(a^2-b^2)|=16
a=28, b=14, n=601, a^2 - b^2=588, |n-(a^2-b^2)|=13
a=29, b=16, n=601, a^2 - b^2=585, |n-(a^2-b^2)|=16
a=30, b=18, n=601, a^2 - b^2=576, |n-(a^2-b^2)|=25
a=31, b=19, n=601, a^2 - b^2=600, |n-(a^2-b^2)|=1
a=32, b=21, n=601, a^2 - b^2=583, |n-(a^2-b^2)|=18
a=33, b=23, n=601, a^2 - b^2=560, |n-(a^2-b^2)|=41
a=34, b=24, n=601, a^2 - b^2=580, |n-(a^2-b^2)|=21
a=35, b=25, n=601, a^2 - b^2=600, |n-(a^2-b^2)|=1
a=36, b=27, n=601, a^2 - b^2=567, |n-(a^2-b^2)|=34
a=37, b=28, n=601, a^2 - b^2=585, |n-(a^2-b^2)|=16
[...]
a=294, b=293, n=601, a^2 - b^2=587, |n-(a^2-b^2)|=14
a=295, b=294, n=601, a^2 - b^2=589, |n-(a^2-b^2)|=12
a=296, b=295, n=601, a^2 - b^2=591, |n-(a^2-b^2)|=10
a=297, b=296, n=601, a^2 - b^2=593, |n-(a^2-b^2)|=8
a=298, b=297, n=601, a^2 - b^2=595, |n-(a^2-b^2)|=6
a=299, b=298, n=601, a^2 - b^2=597, |n-(a^2-b^2)|=4
a=300, b=299, n=601, a^2 - b^2=599, |n-(a^2-b^2)|=2
a=301, b=300, n=601, a^2 - b^2=601, |n-(a^2-b^2)|=0
 ans  =
    31.    601.  
   

Then 18631 = 31 * 601

Factoring into primes

This is a discussion about the algorithm behind the number_factor function.

The Fermat algorithm has some problems when the factor a is far from sqrt(n). The following number is an example :

33554431
= 1801 * 18631
= 1801 * ( 31  *  601 )
   

The "fasttrialdivision" algorithm requires a lot of memory. It fails to factor the Mersenne number M42 = 2^42 - 1 (but successfully factors M41 = 2^41 - 1).

The memorylesstd algorithm does not require any memory. But it is not able to factor the following numbers M53 = 2^53 - 1

Although the memorylesstd algorithm is able to factor M51 = 2^51 - 1, it declares that the status is failure, because it had no time to explore the factors before the maximum number of loops is reached. The same is true for M_47 = 2^47 - 1 correctly factored as 2351.4513.13264529 but declared as failure.

The "pollard-rho" algorithm successively factors the following primes :

==================
Factoring n=2^48-1 
> Success
2^48 - 1 = 3.3.5.7.13.17.97.241.257.673
==================
Factoring n=2^49-1 
> Failure !!!
2^49 - 1 = 127 * 4432676798593 > failure
But 4432676798593 is a prime number.
==================
Factoring n=2^50-1 
2^50 - 1 = 3.11.31.251.601.1801.4051
==================
Factoring n=2^51-1 
2^51 - 1 = 7.103.2143.11119.131071 after 3 attemps
> Factor 1227133513 is not prime
==================
Factoring n=2^52-1 
> Success
2^52 - 1 = 3.5.53.157.1613.2731.8191
=====================
   

After all, this is a success for the pollard-rho method.

Primality test

The primality test for n = 4432676798593 which is a prime, is very long. This number comes as a factor of M_49 = 2^49 - 1 = 127.4432676798593 Only the "210k" algorithm returns an answer.

The probableprime test fails in this case, because intermediate terms involved in the power of a^u modulo n are too large to be stored in 53 bits.

Author

Michael Baudin

Bibliography

“The Art of Computer Programming”, Donald Knuth, Volume 2, Seminumerical Algorithms.

“Introduction to algorithms”, Cormen, Leiserson, Rivest, Stein, 2001, McGraw-Hill.