Miller-Rabin 素数判定と Pollard のρ法

最近 Project Euler の問題を解いているのだけど、使い回せるコードを探すのが大変になってきたのでここにまとめておくことにする。
前に他のプログラム中で使ったものだが、とりあえず Miller-Rabin 素数判定と Pollard のρ法(またはモンテカルロ法)による因数分解と、このふたつを使った素因数分解のコードを多少修正したもの。言語は Ruby。

# Miller-Rabin 法による素数判定
class Miller_Rabin
  class << self
    def mod_pow(b, e, n)
      x = 1
      loop do
        x = x * b % n if e.odd?
        return x if (e >>= 1) == 0
        b = b ** 2 % n
      end
    end

    def init_sd
      @s, @d = 0, @n >> 1
      while @d.even? do
        @s += 1
        @d >>= 1
      end
    end

    def witness?(x)
      x = mod_pow(x, @d, @n)
      return true if x == 1 || x == @n - 1
      @s.times.any? do
        x = x ** 2 % @n
        x == @n - 1
      end
    end

    def probably_prime?(n, k = 20)
      @n = n
      return true if @n == 2 || @n == 3
      return false if @n < 2 || @n.even?
      init_sd
      k.times.all? {witness?(rand(@n - 3) + 2)}
    end
  end
end


# Pollard のρ法による因数分解(Brent のアルゴリズムを使用)
class Rho
  class << self
    M = 64

    def f(x)
      (x ** 2 + @c) % @n
    end

    def init_c
      begin
        @c = rand(@n - 1) + 1
      end while @c == @n - 2
    end

    def brent
      k = 1
      x = y = rand(@n)
      loop do
        k.times {x = f(x)}
        m = [k, M].min
        (k / m).times do
          prod = m.times.inject(1) do |q|
            x = f(x)
            q * (x - y) % @n
          end
          g = @n.gcd(prod)
          return g if g > 1
        end
        y = x
        k <<= 1
      end
    end

    # 非自明な因数をひとつ返す
    def factor(n)
      @n = n
      loop do
        init_c
        q = brent
        return q unless q == @n
      end
    end
  end
end


# 素因数分解
class PF
  class << self
    def factorize(n)
      return if n == 1
      return @d[n] += 1 if Miller_Rabin.probably_prime?(n)
      q = Rho.factor(n)
      factorize(q)
      factorize(n / q)
    end

    # Miller-Rabin 法とρ法により素因数分解を求める
    def prime_division(n)
      @d = Hash.new(0)
      factorize(n)
      @d.sort
    end
  end
end


# 以下使用例

# フェルマー数 F_6
F6 = 2 ** (2 ** 6) + 1

# 素数判定
p Miller_Rabin.probably_prime?(F6)

# 証人数を10として素数判定
p Miller_Rabin.probably_prime?(F6, 10)

# 非自明な因数をひとつ見つける
p Rho.factor(F6)

# 素因数分解
p PF.prime_division(F6)