factorisation -- Sage

# A quick test to ensure that we are dealing with a composite N_max = 100 # upper bound for testing def prime_exclude(n): R = Integers(n) a = R(2) if a != R(0) and a^(n-1) != R(1): return true return false def MR_step( n, a): #print 'Performing a MR step for n, a =', n, a r = 0 d = n-1 while d%2 == 0: d = d/2 r = r+1 # we have 2^r*d = n-1 b = Integers(n)(a) b = b^d #print b if b == Integers(n)(1) : return true s = 0 while b != Integers(n)(-1) and s < r: s = s+1 b = b^2 # print b if s == r : return false if s < r and b == Integers(n)(-1): return true raise ArithmeticError def prime_test_MR( n, prob = 10^(-80) ): if n == 2: return true if n != 2 and n%2 == 0: return false if n == 1: return false #if prime_exclude(n) : # return false prob_fail = 1 while prob_fail >= prob: a = randint(1,n-1) if gcd(a, n) > 1 : return false if MR_step(n, a) == false : return false prob_fail = prob_fail/4 return true def prime_exclude_MR( n ): return not(prime_test_MR(n)) def test_prime_exclude( f ): print 'Checking if the test', f, 'works:' print 'Checking how times test fails up to', N_max, '(must be 0):' , sum([ f(p) for p in range(N_max) if p in Primes()]) N_tot = sum([ 1 for p in range(N_max) if p not in Primes()]) N_excl = sum([ f(p) for p in range(N_max) if p not in Primes()]) print 'Checking the ratio of the number of accurate predictions to the number of all composites:', N_excl, ' : ', N_tot, ' = ', RR(N_excl/N_tot).n(digits = 4) time test_prime_exclude( prime_exclude ) print ' ' print ' ' time test_prime_exclude( prime_exclude_MR )

Checking if the test <function prime_exclude at 0x1709df7c> works:
Checking how times test fails up to 100 (must be 0): 0
Checking the ratio of the number of accurate predictions to the number
of all composites: 74  :  75  =  0.9867
Time: CPU 0.02 s, Wall: 0.19 s
 
 
Checking if the test <function prime_exclude_MR at 0xabd133c>
works:
Checking how times test fails up to 100 (must be 0): 0
Checking the ratio of the number of accurate predictions to the number
of all composites: 75  :  75  =  1.000
Time: CPU 0.53 s, Wall: 0.53 s

Checking if the test <function prime_exclude at 0x1709df7c> works:
Checking how times test fails up to 100 (must be 0): 0
Checking the ratio of the number of accurate predictions to the number of all composites: 74  :  75  =  0.9867
Time: CPU 0.02 s, Wall: 0.19 s
 
 
Checking if the test <function prime_exclude_MR at 0xabd133c> works:
Checking how times test fails up to 100 (must be 0): 0
Checking the ratio of the number of accurate predictions to the number of all composites: 75  :  75  =  1.000
Time: CPU 0.53 s, Wall: 0.53 s

#Some preliminaries x,y = var('x,y') Small_Primes = [2,3,5,7,11,13,17,19,23,31,37,41,43] K_Small_Primes = prod(Small_Primes) # divisible by all small primes def test_correct_factorisation( f , N = 10^10, Z = 10): for z in range(Z): n = randint(floor(N/43),43*N) Fact = f(n) n_prod = prod(Fact) print 'Checking factorisation for', n print 'Acquired the following factorisation:', Fact print 'Product of the factors is correct:', n == n_prod print 'All factors are prime:', prod([ p in Primes() for p in Fact]) == true print ' ' if n != n_prod or prod([ p in Primes() for p in Fact]) != true: print 'The factorisation is wrong! Fail!' raise ArithmeticError print 'All tests passed. The algorithm appears to be correct.' def test_speed_factorisation_step( f , N, Z): for z in range(Z): n = randint(floor(N/43), 43*N) time f(n) def test_speed_factorisation( f, N = 10^10, Z = 10): print 'Testing for time needed to factorise', Z, 'random integer of length approximately ', (log(N)/log(10)).n(digits = 3) time test_speed_factorisation_step( f , N, Z) # returns [n/prod(Small_Factors) , Small_Factors] where Small_Factors are the small prime factors of n def factorise_Small_Primes(n): #print 'Removing small factors from:', n Small_Factors = [] G = gcd(K_Small_Primes, n) while G != 1 : n = n/G for p in Small_Primes : while G % p == 0 : Small_Factors = Small_Factors + [p] G = G/p G = gcd(K_Small_Primes, n) return [n, Small_Factors] # The 'combo' factorising function. I think it makes sense to get rid of the really small primes in a initial stage (for instance to make sure we need not worry that some theorems break down in pathologically small characteristic) def factorise(n): # Of course printing the in-build factor(n) is just so that we know what is going on, it's not used anywhere else. print 'Factorising:', n print 'Real factorisation:', factor(n) #just in case... if n == 0: return [] res_sp = factorise_Small_Primes(n) n = res_sp[0] Fact_sp = res_sp[1] Fact_el = factorise_elliptic(n) Fact = Fact_el + Fact_sp Fact.sort() return Fact

# Some preliminaries for the second stage max_B1 = 2^10 max_B2 = 10*max_B1 Medium_Primes = [] for p in range(max_B2): if prime_test_MR(p, prob = 0.1): Medium_Primes.append(p) # Performs the second stage, trying to multiply the Q from the first stage by all medium primes; the medium primes are taken to include the small ones, because it makes implementation a little bit easier, reduces the number of arbitrary choices and doesn't really cost much def second_stage(Q, B2 = max_B2, B1 = 1): Q = Medium_Primes[0]*Q i = 0 while i + 1 < len(Medium_Primes) and Medium_Primes[i+1] < B2 : Q = ( Medium_Primes[i+1]-Medium_Primes[i] )*Q i = i+1 return [false, Q] print len(Medium_Primes)

# The meat def get_B(p): return floor( exp(sqrt(0.5 * log(p) * log(log(p)) )) ) def get_u(p): return floor( sqrt(2*log(p)/log(log(p)) ) ) def get_max_iter(p): u = get_u(p) return 2*floor( u^u ) # the factor 2 is to be on the safe side def get_times_reuse(p): return floor( log(p)/log(get_B(p)) ) # we want to reuse points around log(p_exp)/log, so that the B-smooth points actually do work fine even when the actual powers of primes in the decomposition are high; calculation shows that we if #C(F_p) is B-smooth we need at most that many trials. I think it makes sense, and if I am wrong this does not spoil performance drastically. # The key factorising function. We assume that n is not zero, and not divisible by (really) small primes def factorise_elliptic(n): if n == 1: return [] p_exp = floor( sqrt(sqrt(n)) ) B = get_B(p_exp) max_iter = get_max_iter(p_exp) # if n is prime, we are done if prime_test_MR(n): return [n] iter = 0 Fact_inter = [false] while Fact_inter == [false] and iter < max_iter: Fact_inter = find_factor_elliptic(n, p_exp) iter = iter+1 if iter == max_iter : #lots of loops, we are out of patience p_exp *= 2 B = get_B(p_exp) max_iter = get_max_iter(p_exp) iter = 0 if B >= max_B1: print 'We know', n, 'is composite but cannot factor it! Fail!' raise ArithmeticError Fact1 = factorise_elliptic(Fact_inter[1]) Fact0 = factorise_elliptic(Fact_inter[0]) Fact = Fact0 + Fact1 #Fact.sort() return Fact def first_stage(P,B): Q = P J = 1 for j in range(1,B): # at stage j, J is lcm(1,...,j) dJ = Integer( j/gcd(J,j) ) J = dJ*J Q = dJ*Q return Q def random_curve(n): R = Integers(n) P = [randint(0,n-1),randint(0,n-1)] a = randint(0,n-1) C = false while C == false: try: C = EllipticCurve(R, [a, P[1]^2 - P[0]^3 - a*P[0]]) except ArithmeticError: P = [randint(0,n-1),randint(0,n-1)] a = randint(0,n-1) P = C(P[0],P[1]) return [P,C] def find_factor_elliptic(n, p_exp): #returns a list: [false] on fail, [n1, n2] on success B = get_B(p_exp) R = Integers(n) [P,C] = random_curve(n) times_reuse = get_times_reuse(p_exp) try: while true: Q = first_stage(P,B) if Q == C([0,1,0]): # order somewhere in the middle, aboriting return [false] else: # Q hopefully has small order; maybe it makes sense to use it as a new P; we want to do it around times_repeat times if randint(1,times_reuse) != 1 : P = Q else: second_stage(Q, B2 = 10*B) return [false] except ZeroDivisionError as err: d = Integer(err.args[0].split()[2]) g = gcd(d,n) if g > 1 and g < n : return [g, n/g] else : return [false]

# The tests for the implemented methods test_correct_factorisation( factorise , N = 10^15, Z = 5) print ' ' test_speed_factorisation( factorise , N = 10^10, Z = 100) print ' ' test_speed_factorisation( factorise , N = 10^15, Z = 100) print ' ' test_speed_factorisation( factorise , N = 10^20, Z = 100) print ' ' test_speed_factorisation( factorise , N = 10^25, Z = 100) #print ' ' #test_speed_factorisation( factorise , N = 10^30, Z = 100) #print ' ' #test_speed_factorisation( factorise , N = 10^35, Z = 100) #print ' ' #test_speed_factorisation( factorise , N = 10^40, Z = 100)

WARNING: Output truncated!  


Factorising: 14198000697931772
Real factorisation: 2^2 * 19 * 186815798656997
Checking factorisation for 14198000697931772
Acquired the following factorisation: [2, 2, 19, 186815798656997]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 41349901926119733
Real factorisation: 3^2 * 23 * 29 * 43 * 137 * 251 * 659 * 7069
Checking factorisation for 41349901926119733
Acquired the following factorisation: [3, 3, 23, 29, 43, 137, 251, 659,
7069]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 14959092803180667
Real factorisation: 3 * 139 * 4447 * 8066814533
Checking factorisation for 14959092803180667
Acquired the following factorisation: [3, 139, 4447, 8066814533]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 34824833285863877
Real factorisation: 41 * 103 * 239 * 5659 * 6097199
Checking factorisation for 34824833285863877
Acquired the following factorisation: [41, 103, 239, 5659, 6097199]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 26631404139199727
Real factorisation: 7 * 43 * 88476425711627
Checking factorisation for 26631404139199727
Acquired the following factorisation: [7, 43, 88476425711627]
Product of the factors is correct: True
All factors are prime: True
 
All tests passed. The algorithm appears to be correct.
 
Testing for time needed to factorise 100 random integer of length
approximately  10.0
Factorising: 150444613198
Real factorisation: 2 * 11 * 6838391509
Time: CPU 0.05 s, Wall: 0.05 s
Factorising: 120588274182
Real factorisation: 2 * 3 * 20098045697
Time: CPU 0.08 s, Wall: 0.08 s
Factorising: 45303514842
Real factorisation: 2 * 3 * 19 * 5233 * 75941
Time: CPU 0.94 s, Wall: 0.93 s
Factorising: 168211892001
Real factorisation: 3 * 29 * 1291 * 1497653
Time: CPU 0.20 s, Wall: 0.21 s
Factorising: 247873124522
Real factorisation: 2 * 127 * 975878443
Time: CPU 0.20 s, Wall: 0.20 s
Factorising: 11792085022
Real factorisation: 2 * 47 * 125447713
Time: CPU 0.16 s, Wall: 0.16 s
Factorising: 240703918691
Real factorisation: 197 * 2281 * 535663
Time: CPU 0.80 s, Wall: 0.80 s

...

Time: CPU 2.75 s, Wall: 2.76 s
Factorising: 210656536442021210253192963
Real factorisation: 3^5 * 7 * 1009 * 592788811 * 207052015237
Time: CPU 35.31 s, Wall: 35.70 s
Factorising: 58201813328635401144240539
Real factorisation: 167 * 3547 * 98255949328242980311
Time: CPU 0.31 s, Wall: 0.31 s
Factorising: 292969351507995956706628601
Real factorisation: 83 * 659 * 5153 * 496831921 * 2092131841
Time: CPU 12.54 s, Wall: 12.62 s
Factorising: 399220167047483432043395260
Real factorisation: 2^2 * 5 * 426990119 * 46748173936957477
Time: CPU 5.87 s, Wall: 5.89 s
Factorising: 54203371007310291194209500
Real factorisation: 2^2 * 3 * 5^3 * 371387 * 97298991810537779
Time: CPU 2.09 s, Wall: 2.10 s
Factorising: 239635518554772151038859390
Real factorisation: 2 * 5 * 11 * 293 * 367 * 2273 * 1539301 * 5790310223
Time: CPU 6.30 s, Wall: 6.33 s
Factorising: 314658643889667305191267670
Real factorisation: 2 * 5 * 53 * 181 * 1913 * 1714629344753132063
Time: CPU 0.35 s, Wall: 0.35 s
Factorising: 32459283541198153813402692
Real factorisation: 2^2 * 3^2 * 19 * 109131973 * 434841335604631
Time: CPU 2.94 s, Wall: 3.00 s
Factorising: 106392501005534743970300973
Real factorisation: 3 * 139 * 349 * 3271 * 5443 * 41061107383277
Time: CPU 0.46 s, Wall: 0.46 s
Factorising: 2225736288761785013426676
Real factorisation: 2^2 * 3 * 498368903 * 372170139322441
Time: CPU 37.37 s, Wall: 37.60 s
Factorising: 160632457915173260394533085
Real factorisation: 3^2 * 5 * 7 * 23 * 22171491775731298881233
Time: CPU 0.07 s, Wall: 0.07 s
Factorising: 283077516570026920027570871
Real factorisation: 641419 * 4790921027 * 92118009367
Time: CPU 4.47 s, Wall: 4.50 s
Factorising: 61066542388723163634314923
Real factorisation: 2789 * 17839 * 1227394890422185313
Time: CPU 0.87 s, Wall: 0.88 s
Factorising: 21553138012788775042541510
Real factorisation: 2 * 5 * 127 * 3023 * 5613951311022000631
Time: CPU 0.27 s, Wall: 0.27 s
Factorising: 199098311637733902809033347
Real factorisation: 79 * 463 * 3324817 * 1637162349509683
Time: CPU 4.22 s, Wall: 4.23 s
Factorising: 219597417719839501807457443
Real factorisation: 11 * 97 * 205808264029840207879529
Time: CPU 0.15 s, Wall: 0.16 s
Factorising: 53441149299070092816181029
Real factorisation: 3 * 13 * 43 * 281 * 113406097778973410017
Time: CPU 0.16 s, Wall: 0.16 s
Factorising: 212109683717985872172972955
Real factorisation: 5 * 59 * 135271 * 215814803 * 24629327273
Time: CPU 23.73 s, Wall: 23.88 s
Factorising: 261454993931786580104869987
Real factorisation: 11 * 285616601 * 83218677516509617
Time: CPU 43.61 s, Wall: 43.83 s
Factorising: 281564650142481547494967638
Real factorisation: 2 * 3 * 263217633529 * 178283806678337

WARNING: Output truncated!  


Factorising: 14198000697931772
Real factorisation: 2^2 * 19 * 186815798656997
Checking factorisation for 14198000697931772
Acquired the following factorisation: [2, 2, 19, 186815798656997]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 41349901926119733
Real factorisation: 3^2 * 23 * 29 * 43 * 137 * 251 * 659 * 7069
Checking factorisation for 41349901926119733
Acquired the following factorisation: [3, 3, 23, 29, 43, 137, 251, 659, 7069]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 14959092803180667
Real factorisation: 3 * 139 * 4447 * 8066814533
Checking factorisation for 14959092803180667
Acquired the following factorisation: [3, 139, 4447, 8066814533]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 34824833285863877
Real factorisation: 41 * 103 * 239 * 5659 * 6097199
Checking factorisation for 34824833285863877
Acquired the following factorisation: [41, 103, 239, 5659, 6097199]
Product of the factors is correct: True
All factors are prime: True
 
Factorising: 26631404139199727
Real factorisation: 7 * 43 * 88476425711627
Checking factorisation for 26631404139199727
Acquired the following factorisation: [7, 43, 88476425711627]
Product of the factors is correct: True
All factors are prime: True
 
All tests passed. The algorithm appears to be correct.
 
Testing for time needed to factorise 100 random integer of length approximately  10.0
Factorising: 150444613198
Real factorisation: 2 * 11 * 6838391509
Time: CPU 0.05 s, Wall: 0.05 s
Factorising: 120588274182
Real factorisation: 2 * 3 * 20098045697
Time: CPU 0.08 s, Wall: 0.08 s
Factorising: 45303514842
Real factorisation: 2 * 3 * 19 * 5233 * 75941
Time: CPU 0.94 s, Wall: 0.93 s
Factorising: 168211892001
Real factorisation: 3 * 29 * 1291 * 1497653
Time: CPU 0.20 s, Wall: 0.21 s
Factorising: 247873124522
Real factorisation: 2 * 127 * 975878443
Time: CPU 0.20 s, Wall: 0.20 s
Factorising: 11792085022
Real factorisation: 2 * 47 * 125447713
Time: CPU 0.16 s, Wall: 0.16 s
Factorising: 240703918691
Real factorisation: 197 * 2281 * 535663
Time: CPU 0.80 s, Wall: 0.80 s

...

Time: CPU 2.75 s, Wall: 2.76 s
Factorising: 210656536442021210253192963
Real factorisation: 3^5 * 7 * 1009 * 592788811 * 207052015237
Time: CPU 35.31 s, Wall: 35.70 s
Factorising: 58201813328635401144240539
Real factorisation: 167 * 3547 * 98255949328242980311
Time: CPU 0.31 s, Wall: 0.31 s
Factorising: 292969351507995956706628601
Real factorisation: 83 * 659 * 5153 * 496831921 * 2092131841
Time: CPU 12.54 s, Wall: 12.62 s
Factorising: 399220167047483432043395260
Real factorisation: 2^2 * 5 * 426990119 * 46748173936957477
Time: CPU 5.87 s, Wall: 5.89 s
Factorising: 54203371007310291194209500
Real factorisation: 2^2 * 3 * 5^3 * 371387 * 97298991810537779
Time: CPU 2.09 s, Wall: 2.10 s
Factorising: 239635518554772151038859390
Real factorisation: 2 * 5 * 11 * 293 * 367 * 2273 * 1539301 * 5790310223
Time: CPU 6.30 s, Wall: 6.33 s
Factorising: 314658643889667305191267670
Real factorisation: 2 * 5 * 53 * 181 * 1913 * 1714629344753132063
Time: CPU 0.35 s, Wall: 0.35 s
Factorising: 32459283541198153813402692
Real factorisation: 2^2 * 3^2 * 19 * 109131973 * 434841335604631
Time: CPU 2.94 s, Wall: 3.00 s
Factorising: 106392501005534743970300973
Real factorisation: 3 * 139 * 349 * 3271 * 5443 * 41061107383277
Time: CPU 0.46 s, Wall: 0.46 s
Factorising: 2225736288761785013426676
Real factorisation: 2^2 * 3 * 498368903 * 372170139322441
Time: CPU 37.37 s, Wall: 37.60 s
Factorising: 160632457915173260394533085
Real factorisation: 3^2 * 5 * 7 * 23 * 22171491775731298881233
Time: CPU 0.07 s, Wall: 0.07 s
Factorising: 283077516570026920027570871
Real factorisation: 641419 * 4790921027 * 92118009367
Time: CPU 4.47 s, Wall: 4.50 s
Factorising: 61066542388723163634314923
Real factorisation: 2789 * 17839 * 1227394890422185313
Time: CPU 0.87 s, Wall: 0.88 s
Factorising: 21553138012788775042541510
Real factorisation: 2 * 5 * 127 * 3023 * 5613951311022000631
Time: CPU 0.27 s, Wall: 0.27 s
Factorising: 199098311637733902809033347
Real factorisation: 79 * 463 * 3324817 * 1637162349509683
Time: CPU 4.22 s, Wall: 4.23 s
Factorising: 219597417719839501807457443
Real factorisation: 11 * 97 * 205808264029840207879529
Time: CPU 0.15 s, Wall: 0.16 s
Factorising: 53441149299070092816181029
Real factorisation: 3 * 13 * 43 * 281 * 113406097778973410017
Time: CPU 0.16 s, Wall: 0.16 s
Factorising: 212109683717985872172972955
Real factorisation: 5 * 59 * 135271 * 215814803 * 24629327273
Time: CPU 23.73 s, Wall: 23.88 s
Factorising: 261454993931786580104869987
Real factorisation: 11 * 285616601 * 83218677516509617
Time: CPU 43.61 s, Wall: 43.83 s
Factorising: 281564650142481547494967638
Real factorisation: 2 * 3 * 263217633529 * 178283806678337

factorisation

171 days ago by jakub.konieczny