require "primes/utils/version"

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

module Primes
  module Utils
    private
    factor_platforms = %w/linux bsd unix/
    os_has_factor = factor_platforms.any? {|p| RUBY_PLATFORM.match(p)}

    public
    if os_has_factor   # for platforms with cli 'factor' command

      def prime?
        `factor #{self.abs}`.split(' ').size == 2
      end

      def factors(p=0) # p is unused dummy variable for method consistency
        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

      puts "Using cli 'factor' for prime?/factors/prime_division"

    else  # use pure ruby versions for platforms without cli command 'factor'

      def prime?
        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 factors(p=13)      # P13 is default prime generator here
        seeds  = [2, 3, 5, 7, 11, 13, 17, 19]
        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.reduce(:*)
        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 [] if n == 0
        return [[n,1]] if primes.include? n
        primes.each {|p| while n%p ==0; factors << p; n /= p end }

        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 if n > 1
        h=Hash.new(0); factors.each {|f| h[f] +=1}; h.to_a.sort
      end

      puts "Using pure ruby versions"
    end
  
    # Replace slow ruby library method prime_division with faster version
    alias  prime_division  factors

    def primenth(p=11) # return value of nth prime
      # Return value of nth prime
      # Uses sozP11 Sieve of Zakiya (SoZ) as default prime sieve

      seeds  = [2, 3, 5, 7, 11, 13, 17, 19]
      primes = []

      n = self.abs # the desired nth prime
      return seeds[n-1] if n < 14
      numb = 22*n  # approx value need to check primes up to

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

      num = numb-1 | 1              # if N even number then subtract 1
      k=num/mod; modk = mod*k; r=1  # kth residue group; base num value
      while num >= modk+residues[r]; r += 1 end # compute extra prms size
      maxprms = k*rescnt + r-1      # max size of prime candidates array
      prms=Array.new(maxprms,true)  # set all prime candidates to True

      # array of residues offsets to compute nonprimes positions in prms
      pos =[]; rescnt.times {|i| pos[residues[i]] = i-1};

      # sieve (SoZ) to eliminate nonprimes from prms
      sqrtN = Math.sqrt(num).to_i
      modk,r,k=0,0,0
      prms.each do |prime|
        r +=1; if r > rescnt; r=1; modk +=mod; k +=1 end
        next unless prime
        res_r = residues[r]
        prime = modk + res_r
        break if prime > sqrtN
        prmstep = prime * rescnt
        residues[1..-1].each do |ri|
          # compute (prime * (modk + ri)) position index in prms
          kk,rr  = (res_r * ri).divmod mod  # residues product res[group|track]
          nonprm = (k*(prime + ri) + kk)*rescnt + pos[rr]
          while nonprm < maxprms; prms[nonprm]=nil; nonprm +=prmstep end
        end
      end
      # the prms array now has all the primes positions for primes r1..N
      # find the nth prime and output it
      count = primes.size
      modk,r=0,0
      prms.each do |prime|
        r +=1; if r > rescnt; r=1; modk +=mod end
        count +=1 if prime
        return modk+residues[r] if count == n
      end
    end

    alias  nthprime  primenth       # to make life easier

    def primes(start_num=0)
      # Find primes between a number range: end_num - start_num
      # Use as: end_num.primes(start_num) or end_num.primes
      # If start_num omitted, will find all primes <= end_num
      # If start_num > self, values are switched to continue
      # Output is an array of the primes within the given range
      # To find count of these primes do: end_num.primes(start_num).size
      # This method uses the P13 Strictly Prime (SP) Prime Generator

      num = self.abs;  start_num = start_num.abs
      num, start_num = start_num, num  if start_num > num

      primes = [2,3,5,7,11,13]      # P13 excluded primes lists

      return primes.select {|p| p >= start_num && p <= num} if num <= 13

      mod = 30030                   # P13 modulus value 2*3*5*7*11*13

      residues=[1]; 17.step(mod,2) {|i| residues << i if mod.gcd(i) == 1}
      rescnt = residues.size        # number of residues
      residues << mod+1             # to make algorithm easier

      num = num-1 | 1               # if N even number then subtract 1
      k=num/mod; modk = mod*k; r=1  # kth residue group; base num value
      while num >= modk+residues[r]; r += 1 end # compute extra prms size
      maxprms = k*rescnt + r-1      # max size of prime candidates array
      prms=Array.new(maxprms,true)  # set all prime candidates to True

      # hash of residues offsets to compute nonprimes positions in prms
      pos =[]; rescnt.times {|i| pos[residues[i]] = i-1};

      # sieve (SoZ) to eliminate nonprimes from prms
      sqrtN= Math.sqrt(num).to_i
      modk,r,k=0,0,0
      prms.each do |prime|
        r +=1; if r > rescnt; r=1; modk +=mod; k +=1 end
        next unless prime
        res_r = residues[r]
        prime = modk + res_r
        break if prime > sqrtN
        prmstep = prime * rescnt
        residues[1..-1].each do |ri|
          # compute (prime * (modk + ri)) position index in prms
          kk,rr  = (res_r * ri).divmod mod  # residues product res[group|track]
          nonprm = (k*(prime + ri) + kk)*rescnt + pos[rr] # 1st prime multiple
          while nonprm < maxprms; prms[nonprm]=nil; nonprm +=prmstep end
        end
      end
      # the prms array now has all the primes positions for primes r1..N
      primes = start_num <= 13 ? primes.drop_while {|p| p < start_num} : []
      k = start_num / mod
      modk = mod*k
      m = rescnt*k
      r = 0
      while m < maxprms
        r +=1; if r > rescnt; r=1; modk +=mod end
        if prms[m]
          prime = modk + residues[r]
	  primes << prime if prime >= start_num
        end
        m +=1
      end
      primes
    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.even?
      k.times do
        a = 2 + rand(n-4)
        x = a.to_bn.mod_exp(d,n)       #x = (a**d) mod n
        next if x == 1 or x == n-1
        (s-1).times do
          x = x.mod_exp(2,n)           #x = (x**2) mod n
          return false if x == 1
          break if x == n-1
        end
        return false if x != n-1
      end
      true  # with high probability
    end
  end
end

class Integer; include Primes::Utils end