Skip to content

Instantly share code, notes, and snippets.

@mjc
Forked from jzakiya/primalitytest6.rb
Last active August 29, 2015 14:27

Revisions

  1. @jzakiya jzakiya created this gist Jul 9, 2015.
    277 changes: 277 additions & 0 deletions primalitytest6.rb
    Original file line number Diff line number Diff line change
    @@ -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)