2019年5月19日日曜日

190519

Ruby


F_2(q) and F_3(q)

B. Mazur, Perturbations, deformations and variations (and "near-misses") in geometry, physics, and number theory
に載っているF_2 およびF_3 を計算してみた。 

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 A(a, b, c, s, t, u, n)
  mul(mul(power(a, s, n), power(b, t, n), n), power(c, u, n), n)
end

def A126858(a, b, c, n)
  a1 = A(a, b, c, 3, 0, 0, n)
  a2 = A(a, b, c, 1, 1, 0, n)
  a3 = A(a, b, c, 0, 0, 1, n)
  (0..n).map{|i| (5 * a1[i] - 3 * a2[i] - 2 * a3[i]) / 51840}
end

def F_3(a, b, c, n)
  a1 = A(a, b, c, 4, 1, 0, n)
  a2 = A(a, b, c, 6, 0, 0, n)
  a3 = A(a, b, c, 2, 2, 0, n)
  a4 = A(a, b, c, 0, 3, 0, n)
  a5 = A(a, b, c, 3, 0, 1, n)
  a6 = A(a, b, c, 1, 1, 1, n)
  a7 = A(a, b, c, 0, 0, 2, n)
  (0..n).map{|i| (15 * a1[i] - 6 * a2[i] - 12 * a3[i] + 7 * a4[i] + 4 * a5[i] - 12 * a6[i] + 4 * a7[i]) / 35831808r}
end

n = 100
a = E_2k(1, n)
b = E_2k(2, n)
c = E_2k(3, n)
p A126858(a, b, c, n)
p F_3(a, b, c, n)

出力結果
[0, 0, 1, 8, 30, 80, 180, 336, 620, 960, 1590, 2200, 3416, 4368, 6440, 7920, 11160, 13056, 18333, 20520, 27860, 31360, 41052, 44528, 59760, 62400, 80990, 87120, 109872, 113680, 147960, 148800, 188976, 196416, 240210, 243040, 311910, 303696, 376580, 385840, 467400, 459200, 578592, 556248, 685960, 689040, 814200, 795616, 1000928, 940800, 1142575, 1138320, 1346436, 1289808, 1597320, 1502160, 1815520, 1781440, 2071470, 1984760, 2474640, 2269200, 2709152, 2665488, 3108960, 2941120, 3587760, 3307656, 3968412, 3841920, 4433520, 4174800, 5145660, 4667328, 5518550, 5396600, 6208440, 5785472, 7030296, 6408480, 7665680, 7318080, 8336202, 7813288, 9580480, 8641440, 10093820, 9750960, 11196240, 10338240, 12560670, 11313120, 13383056, 12745216, 14427120, 13505200, 16390080, 14601216, 17071257, 16382520, 18802050]
[(0/1), (0/1), (1/12), (20/3), (102/1), (2288/3), (3773/1), (14232/1), (133616/3), (119904/1), (584517/2), (1927900/3), (4013432/3), (2569296/1), (14394518/3), (8365192/1), (14426496/1), (23381600/1), (151885575/4), (58125708/1), (269849564/3), (395149888/3), (195967551/1), (828880856/3), (398774464/1), (544543680/1), (4586626939/6), (1018905048/1), (1396485648/1), (5453010736/3), (2448671142/1), (3121442400/1), (4128375808/1), (5187906080/1), (13495484403/2), (25020191728/3), (10739989086/1), (13073496336/1), (49818496739/3), (60139107448/3), (25150090080/1), (89915259200/3), (37377803512/1), (44031127764/1), (163015794472/3), (63689848272/1), (77792286390/1), (270800856112/3), (329864101376/3), (126331201920/1), (1829900154535/12), (174829519992/1), (209303207172/1), (237830107920/1), (284352406578/1), (320567905752/1), (1141395757696/3), (1285519158688/3), (1009346752689/2), (1692885599372/3), (665140467376/1), (738055938000/1), (2591163670808/3), (959464286328/1), (1115821891584/1), (3692188459456/3), (1432795280676/1), (1570557014940/1), (1818030330092/1), (1995407460928/1), (2295038099676/1), (2504195399000/1), (2888008439184/1), (3131191652832/1), (21533911590655/6), (11708697709580/3), (4449903597384/1), (14425877365952/3), (5497041209758/1), (5909603218032/1), (20169906407168/3), (7246380240288/1), (16391092220607/2), (26370350052964/3), (29945483611456/3), (10645119258336/1), (36058977478421/3), (12870138877544/1), (14466994073664/1), (15402512851936/1), (34752009094569/2), (18413605256592/1), (62040662272400/3), (65985397699840/3), (24568961270988/1), (78053912323960/3), (29202477969408/1), (30753913227744/1), (137364461667811/4), (36345626985348/1), (40420330609290/1)]

0 件のコメント:

コメントを投稿

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