171224

Ruby


ラマヌジャンのタウ函数のいろいろな計算方法(1)

この記事は
インテジャーズ Advent Calendar 2017 (https://adventar.org/calendars/2342)

の12/24 分として書いております。
前回はコロちゃんぬさんの「積分を使って自然数のべき乗の和を求める私」でした。

ラマヌジャンのタウ函数τ(n) は、

q*Procuct_{n>0} (1-q^n)^24 = Sum_{n>0} τ(n)*q^n

で求められます。

基本的な性質として、

① gcd(m, n) = 1 のとき、τ(m*n) ＝ τ(m)*τ(n).
② p が素数で r > 0 のとき、τ(p^(r + 1)) = τ(p)*τ(p^r) − p^11*τ(p^(r − 1))．

があります。

また、インテジャーズ にはτ(n)(http://integers.hatenablog.com/archive/category/%CF%84%28n%29)

に関するいろいろな事柄が載っています。

しかし、実際τ(n) をどのように求めればよいでしょうか？

性質①②を使えるときはいいのですが、τ(p)  (pは素数) を求めるには、

どうすればいいでしょうか？

方法1

定義に従う。

(Product_{n>0} (1-q^n))^24 のp - 1 次以下の項を計算する。

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  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[0..m]
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  k = power(ary, n >> 1, m)
  k = mul(k, k, m)
  return k if n & 1 == 0
  return mul(k, ary, m)
end

def A(n)
  ary = [1]
  (1..n).each{|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary, n)
  }
  power(ary, 24, n)
end

n = 19
p A(n - 1)

出力結果
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920, 534612, -370944, -577738, 401856, 1217160, 987136, -6905934, 2727432, 10661420]

方法2
ヤコビの公式を用いる。（「数学の微笑み」山下純一著 p.289 ）

(Product_{n>0} (1-q^n))^3 の8乗のp - 1 次以下の項を計算する。

# m次以下を取り出す
def mul(f_ary, b_ary, m)
  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[0..m]
end

# m次以下を取り出す
def power(ary, n, m)
  return [1] if n == 0
  k = power(ary, n >> 1, m)
  k = mul(k, k, m)
  return k if n & 1 == 0
  return mul(k, ary, m)
end

def A(n)
  ary = Array.new(n + 1, 0)
  # ヤコビの公式の必要なところだけ取り出す
  i = 0
  j, k = 2 * i + 1, i * (i + 1) / 2
  while k <= n
    i & 1 == 1? ary[k] = -j : ary[k] = j
    i += 1
    j, k = 2 * i + 1, i * (i + 1) / 2
  end
  power(ary, 8, n)
end

n = 19
p A(n - 1)

出力結果
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920, 534612, -370944, -577738, 401856, 1217160, 987136, -6905934, 2727432, 10661420]

方法3
{-24, -24, -24, ... } のEuler transform を利用する。

171119 分において、
f(n) = -m, g(n) = 1
とすると、次のようになる。

Product_{n>0} (1-q^n)^m = a(0) + a(1)*q + a(2)*q^2 + ... .
このとき、
a(n) = -(m/n) * Sum_{k=1..n} b(k)*a(n-k).
(ただし、b(n) = Sum_{d|n} d)

この方法はm乗しているのに、計算量が少なくてすむので、良い方法と言えます。

以下は、b(n) について公式を使わず単純に約数を足していったコードです。

def s(n)
  s = 0
  (1..n).each{|i| s += i if n % i == 0}
  s
end

def A(m, n)
  ary = [1]
  a = [0] + (1..n).map{|i| s(i)}
  (1..n).each{|i| ary << (1..i).inject(0){|s, j| s - m * a[j] * ary[-j]} / i}
  ary
end

n = 19
p A(24, n - 1)

出力結果
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920, 534612, -370944, -577738, 401856, 1217160, 987136, -6905934, 2727432, 10661420]

この方法が速いとしばらく思っていたのですが、
今年（2017）9月に「ラマヌジャン《ゼータ関数論文集》」p.112
に次の方法が載っていて見つけて、衝撃を受けました。
(もちろん方法3 も載っていました。さすがですね、Ramanujan は！)

方法4
次の等式を利用する。

1*(n - 1)*τ(n - 0) - 3*(n - 10)*τ(n - 1) + 5(n - 28)*τ(n - 3) - 7*(n-55)*τ(n - 6) + ... = 0
ここで、
数列0, 1, 3, 6, ... の第r項はr*(r - 1)/2 であり、
数列1, 10, 28, 55, ... の第r項は1+9*r*(r - 1)/2.

def A(n)
  ary = [0, 1]
  (2..n).each{|i|
    s, t, u = 0, 1, 0
    (1..n).each{|j|
      t += 9 * j
      u += j
      break if i <= u
      s += (-1) ** (j % 2 + 1) * (2 * j + 1) * (i - t) * ary[-u]
    }
    ary << s / (i - 1)
  }
  ary[1..-1]
end

n = 458329
p A(n)[-1]

出力結果
-11695495424911987900947041440697

方法5
方法4 に対し、性質①②を使う。

require 'prime'

def A(n)
  ary = [0, 1]
  (2..n).each{|i|
    a = i.prime_division
    # i = p^q
    if a.size == 1
      s, t, u = 0, 1, 0
      (1..n).each{|j|
        t += 9 * j
        u += j
        break if i <= u
        s += (-1) ** (j % 2 + 1) * (2 * j + 1) * (i - t) * ary[-u]
      }
      ary << s / (i - 1)
    else
      s = 1
      a.each{|j|
        k, l = j[0], j[1]
        if l == 1
          s *= ary[k]
        else
          s *= ary[k ** (l - 1)] * ary[k] - k ** 11 * ary[k ** (l - 2)]
        end
      }
      ary << s
    end
  }
  ary[1..-1]
end

n = 458329
p A(n)[-1]

出力結果
-11695495424911987900947041440697