2016年1月3日日曜日

160103(2)

Ruby


A190945(100)(2)

二つコードを公開します。

一つ目は、A190945(1), A190945(2), … と求める方法です。
二つ目は、OEIS に載っていたTani Akinari さんのMaxima のコードを
Ruby に置きなおしたものです。
一つ目と異なり、A190945(n) を直接求める方法です。
(n が)1から順に求める必要がないので、検証用として使用しました。

def C(n, r)
  r = [r, n - r].min
  return 1 if r == 0
  return n if r == 1
  numerator = (n - r + 1..n).to_a
  denominator = (1..r).to_a
  (2..r).each{|p|
    pivot = denominator[p - 1]
    if pivot > 1
      offset = (n - r) % p
      (p - 1).step(r - 1, p){|k|
        numerator[k - offset] /= pivot
        denominator[k] /= pivot
      }
    end
  }
  result = 1
  (0..r - 1).each{|k|
    result *= numerator[k] if numerator[k] > 1
  }
  return result
end

=begin
隣接をm個除去しておいてd + m個隣接を増やす入れ方を考える。
つまり、
隣接箇所(k)からm箇所選んでに数字を入れ
非隣接箇所(s - k)からmm(= i - d - 2m)箇所数字を入れ
入れた数字(m + mm)にd + m個数字を重複的に追加する方法の数である。

調べる範囲を整理しておく。
上記より、
① 0 ≦ m ≦ k
② 0 ≦ mm ≦ s - k
③ m + mm ≠ 0
④ d + m ≧ 0
となるが、②については、①および④より、
s - k
≧ (i - 1) * i / 2 + 1 - (i - 1) * (i - 2) / 2 = i
≧ i - (d + m) - m = mm
となるので、0 ≦ mm でよい。
=end
n = 100
# iまでの数字を使った並びで隣接がj箇所ある並び方がa[j]
a = [1]
(2..n).each{|i|
  b = []
  s = (i - 1) * i / 2 + 1
  (0..s - 1).each{|j|
    x = 0
    (0..(i - 1) * (i - 2) / 2).each{|k|
      d = j - k
      (0..k).each{|m|
        if d + m >= 0
          mm = i - d - 2 * m
          # a[k] * C(k, m) * C(s - k, mm) * H(m + mm, d + m)
          x += a[k] * C(k, m) * C(s - k, mm) * C(i - 1, d + m) if 0 <= mm && m + mm != 0
        end
      }
    }
    b[j] = x
  }
  a = b
}

p a[0]

def C(n, r)
  r = [r, n - r].min
  return 1 if r == 0
  return n if r == 1
  numerator = (n - r + 1..n).to_a
  denominator = (1..r).to_a
  (2..r).each{|p|
    pivot = denominator[p - 1]
    if pivot > 1
      offset = (n - r) % p
      (p - 1).step(r - 1, p){|k|
        numerator[k - offset] /= pivot
        denominator[k] /= pivot
      }
    end
  }
  result = 1
  (0..r - 1).each{|k|
    result *= numerator[k] if numerator[k] > 1
  }
  return result
end

def mul(f_ary, b_ary)
  s1, s2 = f_ary.size, b_ary.size
  ary = Array.new(s1 + s2 - 1, 0)
  (0..s1 - 1).each{|i|
    (0..s2 - 1).each{|j|
      ary[i + j] += f_ary[i] * b_ary[j]
    }
  }
  ary
end

def f(n)
  return 1 if n < 2
  (1..n).inject(:*)
end

def A(n)
  ary = [1]
  (1..n).each{|i|
    ary = mul(ary, [0] + (1..i).map{|j| (-1) ** (i - j) * C(i - 1, i - j) / f(j).to_r})
  }
  (n..n * (n + 1) / 2).inject(0){|s, i| s + (1..i).inject(:*) * ary[i]}.to_i
end

p A(100)

0 件のコメント:

コメントを投稿

注: コメントを投稿できるのは、このブログのメンバーだけです。