Skip to content

Instantly share code, notes, and snippets.

@mjc
Forked from jzakiya/primalitytest6.rb
Last active Aug 29, 2015
Embed
What would you like to do?
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment