Rubyで三角数の約数を探せ！〜Rubyでオイラープロジェクトを解こう！Problem12

Problem 12 - Project Eulerより

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...
Let us list the factors of the first seven triangle numbers:
1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
We can see that 28 is the first triangle number to have over five divisors.
What is the value of the first triangle number to have over five hundred divisors?
三角数の数列は自然数を足し合わせていくことで作られる。故に7番目の三角数は1 + 2 + 3 + 4 + 5 + 6 + 7 = 28となる。最初の10項は1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...となる。
最初の7つの三角数の因数を並べてみよう。
1: 1
3: 1,3
6: 1,2,3,6
10: 1,2,5,10
15: 1,3,5,15
21: 1,3,7,21
28: 1,2,4,7,14,28
そうすると28は約数が5つを超える最初の三角数であることがわかる。
では約数が500を超える最初の三角数の値は何か。

ちょっと悩んだけど以下の戦略で

factorメソッドで順番に三角数の因数*1を求める
Array#combinationメソッドで因数同士の全て組み合わせを求める
その組み合わせから約数*2を求める

 require "mathn"
 def triangle_number(divs)
   i = 1; tri = 1
   loop do
     divisors = []
     factors = factor(tri)
     1.upto(factors.length-1) do |n|
       divisors += factors.combination(n).to_a
     end
     result = divisors.map { |c| c.inject(:*) }.uniq
     return result if result.length > divs
     i += 1; tri += i
   end
 end

 def factor(n)
   return [1, n] if Prime.prime?(n)
   2.upto(n-1) do |i|
     return factor(n/i) << i if n.modulo(i).zero?
   end
 end

 triangle_number(500).max # => 76576500

ちょっと遅いです…
結果が出たということで…

(追記:2009/1/22)
Problem 21 - Project Eulerを解いていたら
ここでバカみたいに複雑なことやっているのに気付いた

方針：

三角数は (1..n).inject(:+)で求まる
約数は三角数を割ってあまりが出ないもので求まる

 def tr_with_divs(limit)
   i = 1
   loop do
     tr = (1..i).inject(:+)
     return tr if div_num(tr) > limit
     i += 1
   end
 end

 def div_num(tr)
   num = 0
   (1..tr).each do |i|
     num += 1 if tr % i == 0
   end
   num
 end

 tr_with_divs 500 # =>

ところが…
答えが出てこない！
他の方法が必要だ…

でWikipediaを調べてみると…

約数の個数 - Wikipediaより

自然数 n の全ての正の約数の個数を d(n) で表す。
n の素因数分解が
$n = p1^{a1}p2^{a2}...pm^{am}$
と表せるとき、d(n) は以下の式で求められる。
$d(n) = (a1 + 1)(a2 +1)...(am + 1)$

じゃあ今度はこの方向で…

 require "mathn"
 def tr_with_divs(limit)
   i = 1
   loop do
     tr = (1..i).inject(:+)
     return tr if div_num(tr) > limit
     i += 1
   end
 end

 def div_num(tr)
   cnt = Hash.new(0)
   f = factor(tr)
   f.each { |e| cnt[e] += 1 }
   cnt.values.map { |e| e + 1 }.inject(:*)
 end

 def factor(n)
   return [n] if Prime.prime?(n)
   (2...n).each do |i|
     return factor(n/i) << i if n % i == 0
   end
 end

 t = Time.now
 tr_with_divs 500 # => 76576500
 Time.now - t  3 # => 330.58157

やっぱり時間が掛かる…

*1:ここでは三角数を積の形で表したときの要素とします

*2:ここでは三角数を割りきれる数とします