Number Notes — Some notes on the Number toolbox.
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
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.
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.