primalitytest6.rb

primalitytest6.rb

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

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

class Integer

   def primz?
      residues = [1,7,11,13,17,19,23,29,31,37,41,43,47,49,53,59,61]
      res1 = [1,13,17,29,37,41,49,53]
      res2 = [7,19,31,43]
      res3 = [11,23,47,59]

      mod=60; rescnt=16

      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 ( ! residues.include?(n%mod) || n < 2 )

      m= n%12; res=[]
      res = res1 if (m == 1 or m == 5)
      res = res2 if m == 7
      res = res3 if m == 11
      return false if res.size == 0

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

   def primz7?
      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 ( ! residues.include?(n%mod) || n < 2 )

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

   def primz7a?
      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 < 2

      sqrtN = Math.sqrt(n).to_i
      p=11          # first prime candidate 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    # next prime candidates from next kth residue group
      end
      return true
   end

   def primz7b?
      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  # use prime candidates from next kth residue group
      end
      return true
   end

   def primepa?(p=17)
     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 < 10009

     # 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 !residues.include?(n%mod) || n < 2

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

   def primep?(p=17)
     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
     return true  if primes.include? n
     return false if !residues.include?(n%mod) || n < 2

     sqrtN = Math.sqrt(n).to_i
     modk,r=0,1; p=residues[1] # firts prime candidate
     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 factorp(p=17)
     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] # firts prime candidate
     while p <= sqrtN
       if n%p == 0
         factors << p; modk,r=0,0; n /= p; sqrtN = Math.sqrt(n).to_i
       end
       r +=1; if r > rescnt; r=1; modk +=mod end
       p = modk+residues[r]   # next prime candidate
     end
     factors << n
     factors.sort # return n if prime, or its prime factors
   end

   # Miller-Rabin prime test in Ruby
   # From: http://en.wikipedia.org/wiki/Miller-Rabin_primality_test

   def primemr?
     n = self.abs()
     return true  if n == 2
     return false if n == 1 || n & 1 == 0
     return false if n > 3 && n % 6 != 1 && n % 6 != 5     # added

     # cf. http://betterexplained.com/articles/another-look-at-prime-numbers/ and
     # http://everything2.com/index.pl?node_id=1176369

     d = n-1
     d >>= 1 while d & 1 == 0
     20.times do                               # 20 = k from above
       a = rand(n-2) + 1;   t = d
       y = ModMath.pow(a,t,n)                  # implemented below
       while t != n-1 && y != 1 && y != n-1
         y = (y * y) % n;   t <<= 1
       end
       return false if y != n-1 && t & 1 == 0
     end
     return true
   end

  private
  module ModMath
     def ModMath.pow(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
  end
end

=begin
module ModMath
   def ModMath.pow(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
end
=end

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

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("primz7?      ") do prime.primz7?     end
     t.report("primz7a?     ") do prime.primz7a?    end
     t.report("primz7b?     ") do prime.primz7b?    end
     t.report("primep? 13   ") do prime.primep? 13  end
     t.report("primep? 17   ") do prime.primep? 17  end
     t.report("primepa? 13  ") do prime.primepa? 13 end
     t.report("primepa? 17  ") do prime.primepa? 17 end
     t.report("factorp 13   ") do prime.factorp 13  end
     t.report("factorp 17   ") do prime.factorp 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)