2017年12月31日日曜日

171231

Ruby


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

実はτ(1), τ(2), ... と順に求めなくても
τ(p) (pは素数) を求める方法はあります。

http://www.ipgp.fr/~rozier/math/raman.html
に載っている式をもとに実装すると次のようになります。

require 'prime'

def A259825(n)
  return -1 if n == 0
  return 0 if n % 4 == 1 || n % 4 == 2
  s = 0
  (1..Math.sqrt(n / 3).to_i).each{|a|
    (0..a).each{|b|
      c4 = n + b * b
      if c4 % (4 * a) == 0
        c = c4 / (4 * a)
        if a == c
          if a == b
            s += 4
          elsif b == 0
            s += 6
          else
            s += 12
          end
        elsif a < c
          if a == b
            s += 12
          elsif b == 0
            s += 12
          else
            s += 24
          end
        end
      end
    }
  }
  s
end

def A(n)
  m = 4 * n
  (1..(2 * Math.sqrt(n)).to_i).inject(m ** 5 * A259825(m) / 2){|s, i| s + (m - i * i) ** 5 * A259825(m - i * i)} / 12
end

# nは素数
def B(n)
  A(n) - (462 * n ** 6 + 330 * n ** 4 - 165 * n ** 3 + 55 * n ** 2 - 11 * n + 1)
end

p B(104729)

出力結果
2643202687128887204371152330

このコードは以前のコードに比べ遅いわけではないのですが、
Hurwitz class number の計算アルゴリズムがまずいので、
思ったほど高速になっておらず改善が必要です。

ちなみにPARI で書けば、τ(n) はすぐ求まります。

2017年12月24日日曜日

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

2017年12月9日土曜日

171209(2)

Ruby


Number of Gaussian primes a+b*i such that a^2 + b^2 <= 10^n

A091098 のデータを用いて計算してみた。

require 'prime'

def A(n)
  ary = []
  cnt = 0
  (0..n).each{|i|
    cnt += 1 if i % 4 == 3 && i.prime?
    ary << cnt
  }
  ary
end

ary1 = [0] +
[1,11,80,609,4783,39175,332180,2880504,25423491,
 227523275,2059020280,18803924340,173032709183,
 1602470783672,14922284735484,139619168787795,
 1311778575685086,12369977142579584,
 117028833597800689,1110409801150582707]
ary2 =
[1,3,10,31,100,316,1000,3162,10000,31622,100000,
 316227,1000000,3162277,10000000,31622776,
 100000000,316227766,1000000000,3162277660,
 10000000000,31622776601,100000000000]

n = 8
ary = A(10 ** n)
p (1..2 * n).map{|i| 8 * ary1[i] + 4 * ary[ary2[i]] + 4}

出力結果
[16, 100, 668, 4928, 38404, 313752, 2658344, 23046512, 203394764, 1820205436, 16472216912, 150431552012, 1384262129028, 12819767598972, 119378281788240, 1116953361826164]

171209

Ruby


Number of Gaussian primes a+b*i such that sqrt(a^2 + b^2) <= 10^n

A091098 とA091099 のデータを用いて計算してみた。

ary1 = [0] +
[1,11,80,609,4783,39175,332180,2880504,25423491,
 227523275,2059020280,18803924340,173032709183,
 1602470783672,14922284735484,139619168787795,
 1311778575685086,12369977142579584,
 117028833597800689,1110409801150582707]
ary3 = [0] +
[2,13,87,619,4808,39322,332398,2880950,25424042,
 227529235,2059034532,18803987677,173032827655,
 1602470967129,14922285687184,139619172246129,
 1311778581969146,12369977145161275,
 117028833678543917,1110409801410336132]
p (1..10).map{|i| 8 * ary1[2 * i] + 4 * ary3[i] + 4}

出力結果
[100, 4928, 313752, 23046512, 1820205436, 150431552012, 12819767598972, 1116953361826164, 98959817242332844, 8883278410114778600]