Improved Primality Testing and Factorization in Ruby

primeszp.rb

#!/usr/local/bin/ruby -w

require 'rational' if RUBY_VERSION =~ /^(1.8)/   # for 'gcd' method

=begin
Author:  Jabari Zakiya
Version: June 26, 2013
Version: August 1, 2013  -- shorter/faster/adjustable primemr?
Version: August 28, 2013 -- two new methods use cli command "factor"
Version: September 13, 2013 -- refactored/simplified "primality"
Version: October 18, 2013 -- refactored/simplified "primality"
Description:
For an Integer n, methods primzp*? returns 'true' if n
is prime, or 'false' if not.  The method factorzp
returns an array of the prime factors of n or n if it
is prime.  These methods are significantly faster and
simpler the the standard Ruby methods 'prime' and
'prime_division' found in the library file prime.rb.

Discussion of developing this code can be found below:
http://www.scribd.com/doc/150217723/Improved-Primality-Testing-and-Factorization-in-Ruby?post_id=791539872_10151726037699873#_=_

http://www.4shared.com/dir/7467736/97bd7b71/sharing.html
=end

class Integer

   def primzp5?
      residues = [1,7,11,13,17,19,23,29,31]
      mod=30; rescnt=8

      n = self.abs
      return true  if [2,3,5].include? n
      return false if not residues.include?(n%mod) || n == 1

      sqrtN = Math.sqrt(n).to_i
      modk,r=0,1;  p=7        # first test prime pj
      while p <= sqrtN
        return false if n%p == 0
        r +=1; if r > rescnt; r=1; modk +=mod end
        p = modk+residues[r]  # next prime candidate
      end
      return true
   end

   def primzp5a?
      residues = [1,7,11,13,17,19,23,29,31]
      mod=30; rescnt=8

      n = self.abs
      return true  if [2,3,5].include? n
      return false if not residues.include?(n%mod) || n == 1

      sqrtN = Math.sqrt(n).to_i
      p=7         # first test prime pj
      while p <= sqrtN
        return false if
          n%(p)   == 0 or n%(p+4) ==0 or n%(p+6) == 0 or n%(p+10)==0 or
          n%(p+12)== 0 or n%(p+16)==0 or n%(p+22)== 0 or n%(p+24)==0
        p += mod  # first prime candidate for next kth residues group
      end
      return true
   end

   def primzp7?
      residues = [1,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,
         89,97,101,103,107,109,113,121,127,131,137,139,143,149,151,157,163,
         167,169,173,179,181,187,191,193,197,199,209,211]
      mod=210; rescnt=48

      n = self.abs
      return true  if [2, 3, 5, 7].include? n
      return false if not residues.include?(n%mod) || n == 1

      sqrtN = Math.sqrt(n).to_i
      modk,r=0,1;  p=11        # first test prime pj
      while p <= sqrtN
        return false if n%p == 0
        r +=1; if r > rescnt; r=1; modk +=mod end
        p = modk+residues[r]   # next prime candidate
      end
      return true
   end
  
   def primzp7a?
      residues = [1,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,
         89,97,101,103,107,109,113,121,127,131,137,139,143,149,151,157,163,
         167,169,173,179,181,187,191,193,197,199,209,211]
      mod=210; rescnt=48

      n = self.abs
      return true  if [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
                       47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103,
                       107, 109, 113, 127, 131, 137, 139, 149, 151, 157,163,
                       167, 173, 179, 181, 191, 193, 197, 199, 211].include? n
      return false if not residues.include?(n%mod) || n == 1

      sqrtN = Math.sqrt(n).to_i
      p=11          # first test prime pj
      while p <= sqrtN
        return false if
          n%(p)    == 0 or n%(p+2)  ==0 or n%(p+6)  == 0 or n%(p+8)  ==0 or
          n%(p+12) == 0 or n%(p+18) ==0 or n%(p+20) == 0 or n%(p+26) ==0 or
          n%(p+30) == 0 or n%(p+32) ==0 or n%(p+36) == 0 or n%(p+42) ==0 or
          n%(p+48) == 0 or n%(p+50) ==0 or n%(p+56) == 0 or n%(p+60) ==0 or
          n%(p+62) == 0 or n%(p+68) ==0 or n%(p+72) == 0 or n%(p+78) ==0 or
          n%(p+86) == 0 or n%(p+90) ==0 or n%(p+92) == 0 or n%(p+96) ==0 or
          n%(p+98) == 0 or n%(p+102)==0 or n%(p+110)== 0 or n%(p+116)==0 or
          n%(p+120)== 0 or n%(p+126)==0 or n%(p+128)== 0 or n%(p+132)==0 or
          n%(p+138)== 0 or n%(p+140)==0 or n%(p+146)== 0 or n%(p+152)==0 or
          n%(p+156)== 0 or n%(p+158)==0 or n%(p+162)== 0 or n%(p+168)==0 or
          n%(p+170)== 0 or n%(p+176)==0 or n%(p+180)== 0 or n%(p+182)==0 or
          n%(p+186)== 0 or n%(p+188)==0 or n%(p+198)== 0 or n%(p+200)==0
        p += mod    # first prime candidate for next kth residues group
      end
      return true
   end

   def primzp7b?
      residues = [1,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,
         89,97,101,103,107,109,113,121,127,131,137,139,143,149,151,157,163,
         167,169,173,179,181,187,191,193,197,199,209,211]
      mod=210; rescnt=48

      n = self.abs
      return true  if [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43,
                       47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103,
                       107, 109, 113, 127, 131, 137, 139, 149, 151, 157,163,
                       167, 173, 179, 181, 191, 193, 197, 199, 211].include? n
      return false if not residues.include?(n%mod) || n == 1

      sqrtN = Math.sqrt(n).to_i
      modk=0
      while (11+modk) <= sqrtN
        return false if
          n%(11+modk) == 0 or n%(13+modk) ==0 or n%(17+modk) == 0 or n%(19+modk) ==0 or
          n%(23+modk) == 0 or n%(29+modk) ==0 or n%(31+modk) == 0 or n%(37+modk) ==0 or
          n%(41+modk) == 0 or n%(43+modk) ==0 or n%(47+modk) == 0 or n%(53+modk) ==0 or
          n%(59+modk) == 0 or n%(61+modk) ==0 or n%(67+modk) == 0 or n%(71+modk) ==0 or
          n%(73+modk) == 0 or n%(79+modk) ==0 or n%(83+modk) == 0 or n%(89+modk) ==0 or
          n%(97+modk) == 0 or n%(101+modk)==0 or n%(103+modk)== 0 or n%(107+modk)==0 or
          n%(109+modk)== 0 or n%(113+modk)==0 or n%(121+modk)== 0 or n%(127+modk)==0 or
          n%(131+modk)== 0 or n%(137+modk)==0 or n%(139+modk)== 0 or n%(143+modk)==0 or
          n%(149+modk)== 0 or n%(151+modk)==0 or n%(157+modk)== 0 or n%(163+modk)==0 or
          n%(167+modk)== 0 or n%(169+modk)==0 or n%(173+modk)== 0 or n%(179+modk)==0 or
          n%(181+modk)== 0 or n%(187+modk)==0 or n%(191+modk)== 0 or n%(193+modk)==0 or
          n%(197+modk)== 0 or n%(199+modk)==0 or n%(209+modk)== 0 or n%(211+modk)==0
        modk += mod  # modulus for next kth residues group prime candidates
      end
      return true
   end

   def primzpa?(p=13)      # P13 is default prime generator here
     seeds  = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
     return 'PRIME OPTION NOT A SEEDS PRIME' if !seeds.include? p

     n = self.abs

     # find primes <= Pn, compute modPn then Prime Gen residues for Pn
     primes = seeds[0..seeds.index(p)]; mod = primes.inject {|a,b| a*b }
     mod = 30 if p > 5 and n < mod+2  # for Pp > P5 and n within Pp residues
     residues=[1]; 3.step(mod,2) {|i| residues << i if mod.gcd(i) == 1}
     residues << mod+1; rescnt = residues.size-1

     return true  if primes.include? n
     return false if not residues.include?(n%mod) || n == 1

     sqrtN = Math.sqrt(n).to_i
     modk = 0;  p=residues[1]             # first prime candidate pj
     res = residues[1..-1].map {|r| r-p } # residues distance from first prime
     while p <= sqrtN
       return false if res.map {|r| n%(r+p)}.include? 0
       p += mod   # first prime candidate for next residues group
     end
     return true
   end

   def primzp?(p=13)       # P13 is default prime generator here
     seeds  = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
     return 'PRIME OPTION NOT A SEEDS PRIME' if !seeds.include? p
     
     n = self.abs
     p = 5 if n < 510513   # use P5 prime generator for small numbers

     # find primes <= Pn, compute modPn then Prime Gen residues for Pn
     primes = seeds[0..seeds.index(p)]; mod = primes.inject {|a,b| a*b }
     residues=[1]; 3.step(mod,2) {|i| residues << i if mod.gcd(i) == 1}
     residues << mod+1; rescnt = residues.size-1

     return true  if primes.include? n
     return false if not residues.include?(n%mod) || n == 1

     sqrtN = Math.sqrt(n).to_i
     modk,r=0,1; p=residues[1] # first test prime pj
     while p <= sqrtN
       return false if n%p == 0
       r +=1; if r > rescnt; r=1; modk +=mod end
       p = modk+residues[r]    # next prime candidate
     end
     return true
   end

   def factorzp(p=13)      # P13 is default prime generator here
     seeds  = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41]
     return 'PRIME OPTION NOT A SEEDS PRIME' if !seeds.include? p

     # find primes <= Pn, compute modPn then Prime Gen residues for Pn
     primes = seeds[0..seeds.index(p)]; mod = primes.inject {|a,b| a*b }
     residues=[1]; 3.step(mod,2) {|i| residues << i if mod.gcd(i) == 1}
     residues << mod+1; rescnt = residues.size-1

     n = self.abs
     factors = []

     return factors << n if primes.include? n
     primes.each {|p| while n%p ==0; factors << p; n /= p end }
     return factors if n == 1  # for when n is product of only seed primes

     sqrtN= Math.sqrt(n).to_i
     modk,r=0,1; p=residues[1] # first test prime pj
     while p <= sqrtN
       if n%p == 0
         factors << p; r -=1; n /= p; sqrtN = Math.sqrt(n).to_i
       end
       r +=1; if r > rescnt; r=1; modk +=mod end
       p = modk+residues[r]    # next (or current) prime candidate
     end
     factors << n
     factors.sort # return n if prime, or its prime factors
   end

   # This produces the same output format as lib method prime_division
   def prime_division_new(p=13) # P13 is default prime generator here
     h=Hash.new(0); factorzp(p).each {|f| h[f] +=1}; h.to_a.sort
   end

   # These two method use the [Un/L]inux cli command "factor"
   # They will perform consistently fast for all *nix based Ruby versions
   def factors
     factors = `factor #{self.abs}`.split(' ')[1..-1].map {|i| i.to_i}
     h = Hash.new(0); factors.each {|f| h[f] +=1}; h.to_a.sort
   end
     
   def primality?
     # return true if number is prime or false otherwise
     `factor #{self.abs}`.split(' ').size == 2
   end

   # Miller-Rabin prime test in Ruby
   # From: http://en.wikipedia.org/wiki/Miller-Rabin_primality_test
   # Ruby Rosetta Code: http://rosettacode.org/wiki/Miller-Rabin_primality_test
   # I modified the Rosetta Code, as shown below
   
   require 'openssl'
   def primemr?(k=20) # increase k for more reliability
     n = self.abs     
     return true  if n == 2 or n == 3
     return false if n % 6 != 1 && n % 6 != 5 or n == 1
     
     d = n - 1
     s = 0
     (d >>= 1; s += 1) while d & 1 == 0  # while d even
     k.times do
       a = 2 + rand(n-4)
       x = OpenSSL::BN::new(a.to_s).mod_exp(d,n)   #x = (a**d) % n
       next if x == 1 or x == n-1
       (s-1).times do
         x = x.mod_exp(2,n)                        #x = (x**2) % n
         return false if x == 1
         break if x == n-1
       end
       return false if x != n-1
     end
     true  # probably
   end
end

def tm; s=Time.now; yield; Time.now-s end  # tm { 10001.primzp?}

require 'benchmark'
require 'prime'
def primetests(prime)
  Benchmark.bmbm(14) do |t|
     t.report("prime tests for P = #{prime}") do    end
     t.report("Miller-Rabin ") do prime.primemr?    end
     t.report("primality?   ") do prime.primality?  end
     t.report("primzp7?     ") do prime.primzp7?    end
     t.report("primzp7a?    ") do prime.primzp7a?   end
     t.report("primzp7b?    ") do prime.primzp7b?   end
     t.report("primzp? 13   ") do prime.primzp? 13  end
     t.report("primzp? 17   ") do prime.primzp? 17  end
     t.report("primzpa? 13  ") do prime.primzpa? 13 end
     t.report("primzpa? 17  ") do prime.primzpa? 17 end
     t.report("factors      ") do prime.factors     end
     t.report("factorzp 13  ") do prime.factorzp 13 end
     t.report("factorzp 17  ") do prime.factorzp 17 end
     t.report("prime? [ruby lib] ") do prime.prime? end
     t.report("prime_division [ruby lib]") do prime.prime_division end
  end
end

prime = 20_000_000_000_000_003  # 17 digits
primetests(prime)

Thanks for the implementation!

However, most of those methods return true for 1:

2.2.4 :001 > load './primzp.rb'
 => true
2.2.4 :002 > 1.primzp5?
 => true
2.2.4 :003 > 1.primzp5a?
 => true
2.2.4 :004 > 1.primzp7?
 => true
2.2.4 :005 > 1.primzp7a?
 => true
2.2.4 :006 > 1.primzp7b?
 => true
2.2.4 :007 > 1.primzpa?
 => true
2.2.4 :008 > 1.primzp?
 => true
2.2.4 :009 > 1.primemr?
 => false