180917

Ruby


The sphere packing problem in dimension 8

Maryna Viazovska氏による講演が以下にアップロードされている。

https://www.youtube.com/watch?v=GSMoIb-D4p4
https://www.youtube.com/watch?v=0a2vdHlbJKc

前者の動画に出てくる
φ = 518400q + 31104000q^2 + 870912000q^3 + ... ,
ψ = q^(-1) + 144 - 5120q^(1/2) + 70524q - 626688q^(3/2) + ...
（ただし、q = exp(2π i t))
の係数を出力してみた。

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 bernoulli(n)
  ary = []
  a = []
  (0..n).each{|i|
    a << 1r / (i + 1)
    i.downto(1){|j| a[j - 1] = j * (a[j - 1] - a[j])}
    ary << a[0] # Bn = a[0]
  }
  ary
end

def E_2k(k, n)
  a = -4 * k / bernoulli(2 * k)[-1]
  b = a.denominator
  c = a.numerator
  [b] + (1..n).map{|i| c * sigma(2 * k - 1, i)}
end

# 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 I(ary, n)
  a = [1]
  i = 0
  while i < n
    a << -(0..i).inject(0){|s, j| s + ary[1 + i - j] * a[j]}
    i += 1
  end
  a
end

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

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

def A(n)
  a = E_2k(1, n)
  b = E_2k(2, n)
  c = E_2k(3, n)
  d = mul(a, b, n)
  e = (1..n).map{|i| d[i] - c[i]}
  r = A000594(n)
  mul(mul(e, e, n - 1), I(r, n - 1), n - 1)
end

def A00_4(n)
  ary = [1] + [0] * n
  (1..n).each{|i|
    j = i * i
    break if j > n
    ary[j] = 2
  }
  power(ary, 4, n)
end

def A01_4(n)
  ary = [1] + [0] * n
  (1..n).each{|i|
    j = i * i
    break if j > n
    ary[j] = (-1) ** (i % 2) * 2
  }
  power(ary, 4, n)
end

def A10_4(n)
  ary = [1] + [0] * n
  (1..n).each{|i|
    j = i * (i + 1)
    break if j > n
    ary[j] = 1
  }
  power(ary, 4, n)
end

# -2～n-2 
def A319294(n)
  a = A00_4(n)
  b = A01_4(n)
  c = A10_4(n)
  d = (0..n).map{|i| a[i] + b[i]}
  e = I(mul(c, c, n), n)
  f = mul(d, e, n).map{|i| i / 2}
  g = [1] + (1..n).map{|i| b[i] - 16 * c[i - 1]}
  h = mul(g, I(mul(a, a, n), n), n).map{|i| i * 128}
  ([1, 0] + (2..n).map{|i| f[i] + h[i - 2]})[0..n]
end

n = 50
p a = A(n)
p A319294(n + 2)

# A281373
p a.map{|i| i / 518400}

出力結果
[518400, 31104000, 870912000, 15697152000, 210303475200, 2265242112000, 20595617280000, 163319763456000, 1156167168134400, 7434092726323200, 43996666201804800, 242175877622784000, 1250286900875712000, 6095989790740684800, 28230535439298662400, 124774907252150476800, 528508060756732224000, 2152911394014850944000, 8460229761257773363200, 32157626555717287833600, 118510579413170806579200, 424336487248949755392000, 1478950905909909993984000, 5025889410756351414681600, 16677888376529592125280000, 54116367018561184100198400, 171914123063311861284864000, 535276309312266116216832000, 1635212754461027241135244800, 4905795736572272712782438400, 14466381382727713044935884800, 41963590052233149723052032000, 119830807528853675615488512000, 337089191768542107424905984000, 934712156683604751871271731200, 2556376529361173779444614604800, 6899624930489021508910336704000, 18386710788968496778380304896000, 48402789412667653557030560256000, 125927725910609020521612819456000, 323923060378851737609682025228800, 824149276407205470840652886016000, 2074805192196777843429565441536000, 5170241957113295217494592168960000, 12757119890783853676424057997350400, 31177470401176483400583413586124800, 75492971663046873095674098385920000, 181164766163746719677990358085632000, 430983257182115546274947305704403200, 1016666528235478923452749998848640000]
[1, 0, 144, -5120, 70524, -626688, 4265600, -24164352, 119375370, -529539072, 2151757440, -8125793280, 28827864296, -96885780480, 310514729472, -954123868160, 2823202073655, -8074060259328, 22387521828480, -60344692402176, 158484892943628, -406368240128000, 1019049374174976, -2503142549913600, 6030967901500134, -14270261339381760, 33196945984012800, -76000499056847872, 171386626286516040, -381008035914817536, 835625106694349824, -1809276617175275520, 3869780788718208819, -8180990261004369920, 17103874636380714624, -35380669867248697344, 72446742810305722680, -146905267940338876416, 295116441246442028160, -587553474790036682752, 1159710012298896697386, -2270074082571419443200, 4408097952551024955136, -8493898810748820572160, 16245174113200620330624, -30847020761299021213696, 58167128640245120017920, -108947149043693266262016, 202731394027206007021980, -374870460273855469240320, 688938058997189908185600, -1258627893976024198476800, 2286168896136927086572356]
[1, 60, 1680, 30280, 405678, 4369680, 39729200, 315045840, 2230260741, 14340456648, 84870112272, 467160257760, 2411818867430, 11759239565472, 54457051387536, 240692336520352, 1019498573990610, 4152992658207660, 16319887656747248, 62032458633713904, 228608370781579488, 818550322625288880, 2852914556153375760, 9695002721366418624, 32171852578182083575, 104391140082101049576, 331624465785709608960, 1032554609012859020480, 3154345591167105017622, 9463340541227377918176, 27905828284582779793472, 80948283279770736348480, 231155107115844281665680, 650249212516477830680760, 1803071289898928919504768, 4931281885341770407879272, 13309461671468019886015310, 35468192108349723723727440, 93369578342337294670197840, 242916137944847647611135840, 624851582520933135821145882, 1589794128871924133566074240, 4002324830626500469578637040, 9973460565419165157204074400, 24608641764629347369645173606, 60141725310911426312853807072, 145626874350013258286408368800, 349469070531918826539333252480, 831372023885253754388401438473, 1961162284404859034438175152100]