仕事の部屋

2020年3月1日日曜日

200301(2)

Ruby

フーリエ変換による畳み込みを用いた多項式の乗算(2)

フーリエ変換による畳み込みを用いないで、普通に計算してみた。

# 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

def A(k, n)
  ary = []
  x = [1]
  a = [1] * k
  (n + 1).times{|i|
    ary << x
    x = mul(x, a, k * (i + 1))
  }
  ary.flatten
end

def B(k, n)
  ary = []
  x = [1]
  a = (1..k).to_a
  (n + 1).times{|i|
    ary << x
    x = mul(x, a, k * (i + 1))
  }
  ary.flatten
end

def A329708(n)
  ary = []
  (0..n).each{|i|
    a = (1..i + 1).to_a
    ary << mul(a, a, 2 * i)
  }
  ary.flatten
end

(1..5).each{|i| p A(i, 10)}
(1..5).each{|i| p B(i, 10)}
p A329708(10)

出力結果
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 2, 1, 1, 3, 3, 1, 1, 4, 6, 4, 1, 1, 5, 10, 10, 5, 1, 1, 6, 15, 20, 15, 6, 1, 1, 7, 21, 35, 35, 21, 7, 1, 1, 8, 28, 56, 70, 56, 28, 8, 1, 1, 9, 36, 84, 126, 126, 84, 36, 9, 1, 1, 10, 45, 120, 210, 252, 210, 120, 45, 10, 1]
[1, 1, 1, 1, 1, 2, 3, 2, 1, 1, 3, 6, 7, 6, 3, 1, 1, 4, 10, 16, 19, 16, 10, 4, 1, 1, 5, 15, 30, 45, 51, 45, 30, 15, 5, 1, 1, 6, 21, 50, 90, 126, 141, 126, 90, 50, 21, 6, 1, 1, 7, 28, 77, 161, 266, 357, 393, 357, 266, 161, 77, 28, 7, 1, 1, 8, 36, 112, 266, 504, 784, 1016, 1107, 1016, 784, 504, 266, 112, 36, 8, 1, 1, 9, 45, 156, 414, 882, 1554, 2304, 2907, 3139, 2907, 2304, 1554, 882, 414, 156, 45, 9, 1, 1, 10, 55, 210, 615, 1452, 2850, 4740, 6765, 8350, 8953, 8350, 6765, 4740, 2850, 1452, 615, 210, 55, 10, 1]
[1, 1, 1, 1, 1, 1, 2, 3, 4, 3, 2, 1, 1, 3, 6, 10, 12, 12, 10, 6, 3, 1, 1, 4, 10, 20, 31, 40, 44, 40, 31, 20, 10, 4, 1, 1, 5, 15, 35, 65, 101, 135, 155, 155, 135, 101, 65, 35, 15, 5, 1, 1, 6, 21, 56, 120, 216, 336, 456, 546, 580, 546, 456, 336, 216, 120, 56, 21, 6, 1, 1, 7, 28, 84, 203, 413, 728, 1128, 1554, 1918, 2128, 2128, 1918, 1554, 1128, 728, 413, 203, 84, 28, 7, 1, 1, 8, 36, 120, 322, 728, 1428, 2472, 3823, 5328, 6728, 7728, 8092, 7728, 6728, 5328, 3823, 2472, 1428, 728, 322, 120, 36, 8, 1, 1, 9, 45, 165, 486, 1206, 2598, 4950, 8451, 13051, 18351, 23607, 27876, 30276, 30276, 27876, 23607, 18351, 13051, 8451, 4950, 2598, 1206, 486, 165, 45, 9, 1, 1, 10, 55, 220, 705, 1902, 4455, 9240, 17205, 29050, 44803, 63460, 82885, 100110, 112035, 116304, 112035, 100110, 82885, 63460, 44803, 29050, 17205, 9240, 4455, 1902, 705, 220, 55, 10, 1]
[1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 4, 3, 2, 1, 1, 3, 6, 10, 15, 18, 19, 18, 15, 10, 6, 3, 1, 1, 4, 10, 20, 35, 52, 68, 80, 85, 80, 68, 52, 35, 20, 10, 4, 1, 1, 5, 15, 35, 70, 121, 185, 255, 320, 365, 381, 365, 320, 255, 185, 121, 70, 35, 15, 5, 1, 1, 6, 21, 56, 126, 246, 426, 666, 951, 1246, 1506, 1686, 1751, 1686, 1506, 1246, 951, 666, 426, 246, 126, 56, 21, 6, 1, 1, 7, 28, 84, 210, 455, 875, 1520, 2415, 3535, 4795, 6055, 7140, 7875, 8135, 7875, 7140, 6055, 4795, 3535, 2415, 1520, 875, 455, 210, 84, 28, 7, 1, 1, 8, 36, 120, 330, 784, 1652, 3144, 5475, 8800, 13140, 18320, 23940, 29400, 34000, 37080, 38165, 37080, 34000, 29400, 23940, 18320, 13140, 8800, 5475, 3144, 1652, 784, 330, 120, 36, 8, 1, 1, 9, 45, 165, 495, 1278, 2922, 6030, 11385, 19855, 32211, 48879, 69675, 93600, 118800, 142740, 162585, 175725, 180325, 175725, 162585, 142740, 118800, 93600, 69675, 48879, 32211, 19855, 11385, 6030, 2922, 1278, 495, 165, 45, 9, 1, 1, 10, 55, 220, 715, 1992, 4905, 10890, 22110, 41470, 72403, 118360, 182005, 264220, 363165, 473694, 587400, 693450, 780175, 837100, 856945, 837100, 780175, 693450, 587400, 473694, 363165, 264220, 182005, 118360, 72403, 41470, 22110, 10890, 4905, 1992, 715, 220, 55, 10, 1]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[1, 1, 2, 1, 4, 4, 1, 6, 12, 8, 1, 8, 24, 32, 16, 1, 10, 40, 80, 80, 32, 1, 12, 60, 160, 240, 192, 64, 1, 14, 84, 280, 560, 672, 448, 128, 1, 16, 112, 448, 1120, 1792, 1792, 1024, 256, 1, 18, 144, 672, 2016, 4032, 5376, 4608, 2304, 512, 1, 20, 180, 960, 3360, 8064, 13440, 15360, 11520, 5120, 1024]
[1, 1, 2, 3, 1, 4, 10, 12, 9, 1, 6, 21, 44, 63, 54, 27, 1, 8, 36, 104, 214, 312, 324, 216, 81, 1, 10, 55, 200, 530, 1052, 1590, 1800, 1485, 810, 243, 1, 12, 78, 340, 1095, 2712, 5284, 8136, 9855, 9180, 6318, 2916, 729, 1, 14, 105, 532, 2009, 5922, 13993, 26840, 41979, 53298, 54243, 43092, 25515, 10206, 2187, 1, 16, 136, 784, 3388, 11536, 31864, 72592, 137638, 217776, 286776, 311472, 274428, 190512, 99144, 34992, 6561, 1, 18, 171, 1104, 5364, 20664, 65100, 170928, 378414, 710828, 1135242, 1538352, 1757700, 1673784, 1303452, 804816, 373977, 118098, 19683, 1, 20, 210, 1500, 8085, 34704, 122520, 363120, 915570, 1980440, 3692140, 5941320, 8240130, 9804240, 9924120, 8433072, 5893965, 3280500, 1377810, 393660, 59049]
[1, 1, 2, 3, 4, 1, 4, 10, 20, 25, 24, 16, 1, 6, 21, 56, 111, 174, 219, 204, 144, 64, 1, 8, 36, 120, 310, 648, 1124, 1608, 1905, 1840, 1376, 768, 256, 1, 10, 55, 220, 690, 1772, 3830, 7040, 11085, 14970, 17203, 16660, 13280, 8320, 3840, 1024, 1, 12, 78, 364, 1335, 4032, 10324, 22776, 43743, 73580, 108558, 140316, 158089, 153672, 126960, 86784, 46848, 18432, 4096, 1, 14, 105, 560, 2345, 8106, 23849, 60860, 136395, 270690, 478051, 753144, 1058715, 1325030, 1469835, 1434076, 1215984, 880320, 528640, 250880, 86016, 16384, 1, 16, 136, 816, 3836, 14896, 49336, 142256, 362086, 821456, 1672056, 3066896, 5081916, 7614096, 10308616, 12583696, 13793761, 13493856, 11673536, 8813056, 5694976, 3055616, 1294336, 393216, 65536, 1, 18, 171, 1140, 5940, 25560, 93900, 300960, 854190, 2169740, 4970250, 10323720, 19517700, 33666840, 53050140, 76370880, 100343385, 120066930, 130377315, 127816740, 112317120, 87578880, 59742720, 34928640, 16957440, 6488064, 1769472, 262144, 1, 20, 210, 1540, 8805, 41544, 167400, 589200, 1840050, 5156600, 13076140, 30190200, 63754850, 123554400, 220231800, 361542480, 546902925, 762066900, 977024850, 1150145700, 1239350265, 1217172600, 1083118800, 866419200, 616358400, 384159744, 205332480, 91095040, 31784960, 7864320, 1048576]
[1, 1, 2, 3, 4, 5, 1, 4, 10, 20, 35, 44, 46, 40, 25, 1, 6, 21, 56, 126, 234, 369, 504, 594, 574, 465, 300, 125, 1, 8, 36, 120, 330, 768, 1544, 2728, 4275, 5920, 7256, 7848, 7386, 5880, 3900, 2000, 625, 1, 10, 55, 220, 715, 1972, 4730, 10040, 19085, 32670, 50553, 70860, 89905, 102820, 105490, 96224, 76775, 52250, 29375, 12500, 3125, 1, 12, 78, 364, 1365, 4332, 11974, 29376, 64818, 129740, 236958, 396516, 609389, 860772, 1117050, 1329584, 1446498, 1430532, 1276546, 1016220, 709125, 422500, 206250, 75000, 15625, 1, 14, 105, 560, 2380, 8526, 26579, 73600, 183645, 417060, 868266, 1665804, 2956345, 4865630, 7437615, 10566136, 13946849, 17084340, 19380690, 20294820, 19525821, 17148254, 13626235, 9672600, 6020000, 3193750, 1378125, 437500, 78125, 1, 16, 136, 816, 3876, 15456, 53536, 164656, 456586, 1154096, 2680616, 5756096, 11479216, 21334096, 37042456, 60192656, 91636211, 130755056, 174842536, 218927296, 256329136, 280028816, 284582936, 267947216, 232466026, 184497760, 132647200, 85218000, 47962500, 22950000, 8875000, 2500000, 390625, 1, 18, 171, 1140, 5985, 26280, 100020, 337680, 1027710, 2852660, 7284870, 17229240, 37932570, 78053760, 150575760, 272977200, 465881355, 749445750, 1137244185, 1628385660, 2199912615, 2802614400, 3363549840, 3797152560, 4023870210, 3991747284, 3693744342, 3173605864, 2516661270, 1827606600, 1202770500, 707040000, 363628125, 158531250, 55546875, 14062500, 1953125, 1, 20, 210, 1540, 8855, 42444, 175950, 646200, 2138175, 6452600, 17924140, 46156200, 110794850, 249009400, 525822300, 1046166480, 1965440925, 3492711900, 5878567350, 9379622700, 14195606265, 20383802100, 27768280050, 35873674200, 43918845525, 50898216744, 55755208980, 57617580040, 56031446210, 51115460520, 43571112676, 34530075200, 25279752375, 16956022500, 10305506250, 5590462500, 2648984375, 1060937500, 339843750, 78125000, 9765625]
[1, 1, 4, 4, 1, 4, 10, 12, 9, 1, 4, 10, 20, 25, 24, 16, 1, 4, 10, 20, 35, 44, 46, 40, 25, 1, 4, 10, 20, 35, 56, 70, 76, 73, 60, 36, 1, 4, 10, 20, 35, 56, 84, 104, 115, 116, 106, 84, 49, 1, 4, 10, 20, 35, 56, 84, 120, 147, 164, 170, 164, 145, 112, 64, 1, 4, 10, 20, 35, 56, 84, 120, 165, 200, 224, 236, 235, 220, 190, 144, 81, 1, 4, 10, 20, 35, 56, 84, 120, 165, 220, 264, 296, 315, 320, 310, 284, 241, 180, 100, 1, 4, 10, 20, 35, 56, 84, 120, 165, 220, 286, 340, 381, 408, 420, 416, 395, 356, 298, 220, 121]

200301

Ruby

フーリエ変換による畳み込みを用いた多項式の乗算(1)

Numo::Pocketfftをgemコマンドでインストールしておく。

require 'numo/pocketfft'
include Numo

def A(k, n)
  ary = []
  x = Int32[1]
  a = Int32[*[1] * k]
  (n + 1).times{
    ary << x.to_a
    # Numo::PocketfftはDFloat
    x = Int32.cast(Numo::Pocketfft.fftconvolve(x, a).map{|i| i.round})
  }
  ary.flatten
end

def B(k, n)
  ary = []
  x = Int64[1]
  a = Int64[1..k]
  (n + 1).times{
    ary << x.to_a
    # Numo::PocketfftはDFloat
    x = Int64.cast(Numo::Pocketfft.fftconvolve(x, a).map{|i| i.round})
  }
  ary.flatten
end

def A329708(n)
  ary = []
  (0..n).each{|i|
    a = Int32[1..i + 1]
    # Numo::PocketfftはDFloat
    ary << Int32.cast(Numo::Pocketfft.fftconvolve(a, a).map{|i| i.round}).to_a
  }
  ary.flatten
end

(1..5).each{|i| p A(i, 10)}
(1..5).each{|i| p B(i, 10)}
p A329708(10)

2020年2月29日土曜日

200229(2)

Crystal

A332800(2)

Conjecture についての計算は以下の通りです。

def a332800_conjecture(n)
  return 1 if n == 0
  a = (1..n).map{|i| [i]}
  len = 1
  while len < n
    b = Array(Array(Int32)).new
    a.each{|c|
      (1..n).each{|num|
        if !c.includes?(num)
          i = c.clone + [num]
          if len > 1
            if (i[-3] % i[-2]) < (i[-2] % i[-1])
              b << i
            end
          else
            b << i
          end
        end
      }
    }
    a = b
    len += 1
  end
  a.size
end

(0..19).each{|i| p a332800_conjecture(i)}

出力結果
1
1
2
3
4
6
7
9
11
13
14
19
20
22
25
29
30
36
37
42

200229

Crystal

A332800(1)

Crystal で計算してみた。途中で計算が止まる。

def a332800(n)
  return 1 if n == 0
  a = (1..n).map{|i| [i]}
  len = 1
  while len < n
    b = Array(Array(Int32)).new
    a.each{|c|
      (1..n).each{|num|
        if !c.includes?(num)
          i = c.clone + [num]
          if len > 1
            if (i[-3] % i[-2]) <= (i[-2] % i[-1])
              b << i
            end
          else
            b << i
          end
        end
      }
    }
    a = b
    len += 1
  end
  a.size
end

(0..20).each{|i| p a332800(i)}

出力結果
1
1
2
4
9
21
44
109
241
530
1176
3180
6456
14835
34672
81877
179434
479275
Program was killed

2020年2月27日木曜日

200227

Crystal

A318790

Crystal で計算してみた。

def check(d, a, i)
  return true if i == 0
  j = 1
  d_max = [i, d - 1].min
  while (a[i] - a[i - j]).abs >= d && j < d_max
    j += 1
  end
  (a[i] - a[i - j]).abs >= d
end

def solve(d, len, a = Array(Int32).new)
  b = Array(Array(Int32)).new
  if a.size == len
    b << a
  else
    (1..len).each{|m|
      s = a.size
      if s == 0 || (s > 0 && !a.includes?(m))
        if check(d, a + [m], s)
          b += solve(d, len, a + [m])
        end
      end
    }
  end
  b
end

def a318790(n)
  (1..n).map{|i| solve(i, i * i + 1).size // 2}
end

p a318790(4)

出力結果
[1, 7, 20, 37]

2020年2月26日水曜日

200226

Crystal

深さが不定の配列の作成

Ruby と違い型の設定が要る。

alias Ary = Int32 | String | Array(Ary)

a = Array(Ary).new
5.times{|i|
  a << i
  a << i.to_s
  p a
  b = Array(Ary).new
  b << a
  a = b
}

出力結果
[0, "0"]
[[0, "0"], 1, "1"]
[[[0, "0"], 1, "1"], 2, "2"]
[[[[0, "0"], 1, "1"], 2, "2"], 3, "3"]
[[[[[0, "0"], 1, "1"], 2, "2"], 3, "3"], 4, "4"]

2020年2月24日月曜日

200224

Crystal

A320843(2)

Crystal で計算してみた。

class Count
  @@cnt = 0
  def initialize
    @@cnt += 1
  end
  def self.instances
    @@cnt
  end
  def self.to_zero
    @@cnt = 0
  end
end

def search(a, num, n)
  if num == n + 1
    Count.new
  else
    (1..n).each{|i|
      if a[i] == 0
        if num % i == 0 || i % num == 0
          a[i] = num
          search(a, num + 1, n)
          a[i] = 0
        end
      end
    }
  end
end

def a320843(n)
  a = [0] * (n + 1)
  Count.to_zero
  search(a, 1, n)
  Count.instances
end

p (0..20).map{|i| a320843(i)}

出力結果
[1, 1, 2, 3, 8, 10, 36, 41, 132, 250, 700, 750, 4010, 4237, 10680, 24679, 87328, 90478, 435812, 449586, 1939684]

2020年2月15日土曜日

200215

Ruby

素数の個数(8)

150709分をRuby2.7で

ruby --jit ○○.rb

とするだけで、10分くらいで計算が終わった。

2020年2月9日日曜日

200209

Crystal

素数の和(2)

以前のRuby のコードをCrystal のコードにかえて出力してみた。

require "big"

def sum_of_primes(n)
  m = Math.sqrt(n).to_i
  keys = (1..m).map{|i| n // i}
  keys += (1..keys[-1] - 1).to_a.reverse
  h = Hash(Int32, BigInt).new
  keys.each{|i|
    j = i.to_big_i
    h[i] = j * (j + 1) // 2 - 1
  }
  (2..m).each{|i|
    if h[i] > h[i - 1] # このときiは素数
      hp = h[i - 1]
      i2 = i * i
        keys.each{|j|
          break if j < i2
          h[j] -= i * (h[j // i] - hp)
        }
    end
  }
  h[n]
end

time = Time.local
p (0..9).map{|i| sum_of_primes(10 ** i)}
puts Time.local - time

出力結果
[0, 17, 1060, 76127, 5736396, 454396537, 37550402023, 3203324994356, 279209790387276, 24739512092254535]
00:00:02.812922000

2020年2月8日土曜日

200208(2)

Crystal

A006880とA046731(2)

10^n までの素数の数と和を出力してみた。

require "big"

def a006880(n)
  max = 10 ** n - 1
  ary = Array.new(max + 1, true)

  i = 2
  while i * i <= max
    if ary[i]
      (2 * i).step(to: max, by: i){|j| ary[j] = false}
    end
    i += 1
  end

  s = 0.to_big_i
  (2..max).each{|i| s += 1 if ary[i]}
  p s
end

def a046731(n)
  max = 10 ** n - 1
  ary = Array.new(max + 1, true)

  i = 2
  while i * i <= max
    if ary[i]
      (2 * i).step(to: max, by: i){|j| ary[j] = false}
    end
    i += 1
  end

  s = 0.to_big_i
  (2..max).each{|i| s += i if ary[i]}
  p s
end

n = 8
p (0..n).map{|i| a006880(i)}
p (0..n).map{|i| a046731(i)}

出力結果
0
4
25
168
1229
9592
78498
664579
5761455
[0, 4, 25, 168, 1229, 9592, 78498, 664579, 5761455]
0
17
1060
76127
5736396
454396537
37550402023
3203324994356
279209790387276
[0, 17, 1060, 76127, 5736396, 454396537, 37550402023, 3203324994356, 279209790387276]

登録: 投稿 (Atom)