2017年7月13日木曜日

170713

Ruby


j-invariant

次の求め方は多項式同士のかけ算をするより速いと思う。

require 'prime'

def power0(a, n)
  return 1 if n == 0
  k = power0(a, n >> 1)
  k *= k
  return k if n & 1 == 0
  return k * a
end

# x > 0
def sigma(x, i)
  sum = 1
  pq = i.prime_division
  pq.each{|a, n| sum *= (power0(a, (n + 1) * x) - 1) / (power0(a, x) - 1)}
  sum
end

def A000521(n)
  s3 = [0] + (1..n + 1).map{|i| sigma(3, i)}
  s5 = [0] + (1..n + 1).map{|i| sigma(5, i)}
  ary = [1]
  (0..n).each{|i| ary << (1..i + 1).inject(0){|s, j| s += (504 * s5[j] - 240 * (i - j) * s3[j]) * ary[-j]} / (i + 1)}
  ary
end

n = 50
p A000521(n)

出力結果
[1, 744, 196884, 21493760, 864299970, 20245856256, 333202640600, 4252023300096, 44656994071935, 401490886656000, 3176440229784420, 22567393309593600, 146211911499519294, 874313719685775360, 4872010111798142520, 25497827389410525184, 126142916465781843075, 593121772421445058560, 2662842413150775245160, 11459912788444786513920, 47438786801234168813250, 189449976248893390028800, 731811377318137519245696, 2740630712513624654929920, 9971041659937182693533820, 35307453186561427099877376, 121883284330422510433351500, 410789960190307909157638144, 1353563541518646878675077500, 4365689224858876634610401280, 13798375834642999925542288376, 42780782244213262567058227200, 130233693825770295128044873221, 389608006170995911894300098560, 1146329398900810637779611090240, 3319627709139267167263679606784, 9468166135702260431646263438600, 26614365825753796268872151875584, 73773169969725069760801792854360, 201768789947228738648580043776000, 544763881751616630123165410477688, 1452689254439362169794355429376000, 3827767751739363485065598331130120, 9970416600217443268739409968824320, 25683334706395406994774011866319670, 65452367731499268312170283695144960, 165078821568186174782496283155142200, 412189630805216773489544457234333696, 1019253515891576791938652011091437835, 2496774105950716692603315123199672320, 6060574415413720999542378222812650932, 14581598453215019997540391326153984000]

3 件のコメント:

  1. http://oeis.org/A000521
    での
    a(n) = (1/n)*(Sum_{r in Z} A027652(n - r^2) + Sum_{r>0, r odd} ((-1)^n * A027652(4*n - r^2) - A027652(16*n - r^2))) for n > 0. - Seiichi Manyama, Jun 11 2017
    から求めるプログラムを見せてください。

    返信削除
    返信
    1. 元々次の論文を見て、j-invariant の計算をしようと思ったのですが、普通に計算しても十分速かったので作っていないです。(ちなみに、GAI さんが書いている公式を使う方が速いみたいです。)
      https://projecteuclid.org/download/pdf_1/euclid.em/1064858788

      削除
  2. http://oeis.org/A027652
    での出力プログラム(-1から10000までの数表提示)を教えてください。

    返信削除