ユークリッド互除法って速いのかしら

nが与えられたとき、

1/x + 1/y = 1/n

を満たす、x、yのユニークな組み合わせ数を調べる。初めて組み合わせ数が1000を超えるnはいくつか？

RubyのComparableみたいなものを使ってみたかったこともあって、分数クラスを作ってブルートフォースでやってみた。gcdとかlcmはあえて再発明というか自分で書いてみた。

module MyMath
  def gcd(a, b)
    if (b == 0) then a
    else
      gcd(b, (a % b))
    end
  end

  def lcm(a, b)
    a * b / gcd(a, b)
  end
end

class Fraction
  include Comparable
  include MyMath
  attr_reader :num, :den

  def initialize(num, den)
    @num = num # numerator
    @den = den # denominator
  end

# reduction(2/4) -> 1/2

  def reduction
    g = gcd(self.num, self.den)
    g == 1 ? self : Fraction.new(self.num / g, self.den / g)
  end

# convert_fraction(2/3, 1/4) -> 8/12, 3/12

  def convert_fraction(f1, f2)
    l = lcm(f1.den, f2.den)
    a, b = [f1, f2].map do |f|
      num = f.num * (l / f.den)
      Fraction.new(num, l)
    end
    return a, b
  end

  def +(other)
    a, b = convert_fraction(self, other)
    sum = Fraction.new(a.num + b.num, a.den).reduction
  end

  def -(other)
    a, b = convert_fraction(self, other)
    sum = Fraction.new(a.num - b.num, a.den).reduction
  end

  def <=>(other)
    (self.num.to_f / self.den) <=> (other.num.to_f / other.den)
  end

  def ==(other)
    (self.num == other.num) && (self.den == other.den)
  end

  def inspect
    self.num.to_s + "/" + self.den.to_s
  end

  def to_s
    self.inspect
  end
end

def solutions(n)
  count = 0
  b = Fraction.new(1, n)

  (n+1).upto(n*2) do |i|
    b1 = Fraction.new(1, i)
    b2 = (b - b1).reduction
    if b2.num == 1 then
      print "#{b1.reduction} + #{b2.reduction} = #{b}\n"
      count += 1
    end
  end
  count
end

# 1.upto(100000) do |i|
#   c = solutions(i)
#   print "#{i}, #{c}\n"
#   if c >= 1000 then break end
# end

puts solutions(ARGV.shift.to_i)

Comparableは分数同士の大小を比較するのに書いたけど、結局考えてみたら、そんな演算子は不要なのだった。ともあれ、これで、

% ruby gcd.rb 4
1/5 + 1/20 = 1/4
1/6 + 1/12 = 1/4
1/8 + 1/8 = 1/4
3
% ruby gcd.rb 12
1/13 + 1/156 = 1/12
1/14 + 1/84 = 1/12
1/15 + 1/60 = 1/12
1/16 + 1/48 = 1/12
1/18 + 1/36 = 1/12
1/20 + 1/30 = 1/12
1/21 + 1/28 = 1/12
1/24 + 1/24 = 1/12
8

という感じでちゃんとできている。しかし、ともかく遅くて1000個なんて無理。クラスとかやめて、Cで書き直したらだいぶ速くなったけど、それでも無理。どうも大きな数字のときのgcdの計算には速いアルゴリズムがあるとかわかったけど、そんな程度の話じゃない。しかし、考えてみたらgcdとか関係ないし、もっと単純に書き換えた。

n = 1

while true
  n += 1
  count = 0
  (n+1).upto(n*2) do |x|
    if (x - n) == 1 ||
      (x * n) % (x - n) == 0 then
      count += 1
    end
  end
  print "#{n}: #{count}\n"
  if count > 1000 then break end
end

これでも20分ほどかかったけど、ちゃんと答えが出た。