|
|
@@ -0,0 +1,277 @@ |
|
|
#!/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) |