Project

General

Profile

Feature #16468

Updated by steveb3210 (Stephen Blackstone) 6 months ago

The miller-rabin algorithm is a non-deterministic primality test, however it is known that below 2**64, you can always get a deterministic answer by only checking a=[2,3,5,7,11,13,17,19,23, 29, 31, 37] 

 Given that Prime.prime? would never respond in a reasonable amount of time for larger numbers, we can gain much more utility and performance by    switching.. 

 ``` 
                                                user       system        total          real 
 miller_rabin: random set                     0.150000     0.000000     0.150000 (    0.152212) 
 Prime.prime?: random set                     0.270000     0.000000     0.270000 (    0.281257) 

                                                user       system        total          real 
 miller_rabin: 16 digits                      0.010000     0.000000     0.010000 (    0.000300) 
 Prime.prime? 16 digits                       2.200000     0.020000     2.220000 (    2.368247) 

                                                user       system        total          real 
 miller_rabin: 2-10000                        0.030000     0.000000     0.030000 (    0.035752) 
 Prime.prime? 2-10000                         0.020000     0.000000     0.020000 (    0.022948) 
 ``` 


 ``` 
 require 'benchmark' 
 require 'prime' 

 def modpow(base, power, mod) 
   result = 1 
   while power > 0 
     result = (result * base) % mod if power & 1 == 1 
     base = (base * base) % mod 
     power >>= 1; 
   end 
   result 
 end 

 def miller_rabin(n) 
   return false if n < 2 # 0, 1 
   return true    if n < 4 # 2, 3 
   return false if n % 2 == 0 # 4, 6, 8, 10 
   d = (n-1) 
   s = 0 
   while (d % 2 == 0) 
     s +=1 
     d /= 2 
   end 
  
  
   # https://arxiv.org/pdf/1509.00864.pdf 
   # 
   # if n < 18,446,744,073,709,551,616 = 2**64, it is enough to test a = 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, and 37. 
   # 
   [2,3,5,7,11,13,17,19,23, 29, 31, 37].each do |a| 
     x = modpow(a,d,n) 
     next if x == 1 or x == (n - 1) or n == a 
     skip = false 
     1.upto(s-1) do |k| 
       x = x*x % n 
       if x == 1 
         return false 
       elsif x == (n - 1) 
         skip = true 
         break 
       end 
     end 
     next if skip 
     return false 
   end 
   return true 
 end 



 test_set = [] 

 while test_set.length < 50_000 
   test_set << rand(1_000_000) 
 end 


 c = 1 

 Benchmark.bm(40) do |bm|   
   bm.report("miller_rabin: random set") do  
     test_set.each do |x| 
       miller_rabin(x) 
     end     
   end 

   bm.report("Prime.prime?: random set") do 
     test_set.each do |x| 
       Prime.prime?(x) 
     end 
   end 
 end 

 puts 


 Benchmark.bm(40) do |bm|   
   bm.report("miller_rabin: 16 digits") do  
     miller_rabin(1000000000100011) 
   end 

   bm.report("Prime.prime? 16 digits") do 
     Prime.prime?(1000000000100011) 
   end 
 end 

 puts 


 Benchmark.bm(40) do |bm|   
   bm.report("miller_rabin: 2-10000") do  
     (2..10000).each do |x| 
       miller_rabin(x) 
     end 
    
   end 

   bm.report("Prime.prime? 2-10000") do 
     (2..10000).each do |x| 
       Prime.prime?(x) 
     end 
   end 
 end 

 ```

Back