数学ガール読了/分割数を求めるプログラム

数学ガールを読み終わったので、記念に分割数Pnを求めるプログラムをrubyで書いてみました。

ソース

# m以上の数を使った場合のnの分割数
def partition_number(n, m=1)
  pn = 1
  (n-m).step((n+1)/2, -1) do |i|
    pn += partition_number(i, n-i)
  end
  pn
end

if $0 == __FILE__
  print(partition_number((ARGV.shift || 15).to_i))
end

このプログラムの処理は、

P(n,m) = {\sum_{k=\frac{n}{2}}^{n-m}{P(k,n-k)} + 1}
を求めています。
ただし、1 ≦ m ≦ n/2 で、n/2は正確には ceil(n/2) です。*1
P(n,m)は 正の整数nをm以上の整数を使って分割する場合の数(ただし、該当する組み合わせがない場合は1)になっているはずです。

P(20)までは、このページの表と答えが一致していることを確認しましたが、それ以上の場合でも正しいかは確認できていません。多分、正しいと思うのですが (^_^;)

計算量とか

ループの再帰呼び出しになっているので、計算量はほぼ階乗のオーダーで、数が増えると計算時間も爆発的に増えます。
以下にこのプログラムで計算した分割数P(n)の値と計算時間T(n)を乗っけます。

n P(n) T(n)[sec]
9 30 0.001
15 176 0.007
20 627 0.002
30 5604 0.018
50 204226 0.690
60 966467 3.282
70 4087968 13.804
80 15796476 53.457
90 56634173 191.500

ところで、

上記のプログラムはとてもシンプルですが、最初はもっと泥臭いやり方でこんな感じでした。

require 'set'

# 整数の分割を表すクラス
class Partition
  attr_reader :n, :part
  
  def self.divider(n, div)
    pa = Array.new(n+1, 0)
    pa[div] = n / div unless div == 0
    
    Partition.new(pa)
  end
  
  def initialize(part)
    @n = part.size - 1
    @part = part
    
    # assertion
    #unless self.check?
    #  raise ArgumentError.new("Illegal aruguments : [#{@part.join(', ')}] != #{@n}")
    #end
  end
  
  # Setによる重複解除去のため、eql?とhashを定義する
  def eql?(other)
    @n == other.n && @part == other.part
  end
  
  def hash
    @n ^ part.hash
  end
  
  def check?
    sum = 0
    @part.each_with_index {|x, i|
      sum += x * i
    }
    sum == @n
  end
  
  def +(other)
    Partition.new(Array.new(@n + other.n + 1) { |i|
      (@part[i] || 0) + (other.part[i] || 0)
    })
  end
end

# Setに集合の積をとるメソッドを追加
class Set
  def cross_join other
    newset = Set.new
    self.each { |a|
      other.each { |b|
        newset << (yield(a, b))
      }
    }
    newset
  end
end

class Fixnum 
  @@partitions_cache = {}
  def partitions
    return Set.new([Partition.divider(self, self)]) if self <= 1
    
    # return cache
    return @@partitions_cache[self] if @@partitions_cache[self]
    
    pset = Set.new

    pset << Partition.divider(self, 1)
    pset << Partition.divider(self, self)
    
    (self-1).step((self+1)/2, -1) do |i|
      pset += i.partitions.cross_join((self-i).paritons) { |a, b| a + b }
    end
    
    # set to cache
    @@partitions_cache[self] = pset
  end
  
  def partition_number
    self.partitions.size
  end
  
  def self.clear_partitions_cache
    @@partitions_cache = {}
    nil
  end
end

このプログラムでは分割式クラスを定義して重複も含めて片っ端から列挙して、Setを使って重複を取り除いてます。最初のプログラムではP(30)が一瞬で求められるのに対して、こちらは1分以上かかってしまいます。
これは、『当たり前のところから始めたかった』からです。*2
数学ガールでは、数学の問題を考えるのに『当たり前のところから始める』ことが大切だと書かれていますが、これはプログラミングにも共通でいえることです。
最初から完璧なものを作ろうとすると行き詰ってしまいます。そこで、泥臭かったり、非効率だったりするけど確実に正しい処理を作ってから、正しい状態を保ちながらより良い形にリファクタリングしてと上手く行くことが多いと思います。



このエントリーももっとリファクタリングした方がいい気がしますが、時間が遅くなってしまったのでもう寝ます。(-_-)zzZZ

*1:ceilは数学的にはどう書くの?

*2:「ミルカさんの解説をちゃんと理解していなかった」からとも言えますが。