仕事の部屋: 9月 2015

2015年9月29日火曜日

150929

Ruby

オイラー関数のベキ(3)

φ(x)^n の1≦n≦10のときについて、
オンライン整数列大辞典の
A010815（http://oeis.org/A010815/list）、
A002107（http://oeis.org/A002107/list）、
A010816（http://oeis.org/A010816/list）、
A000727（http://oeis.org/A000727/list）、
A000728（http://oeis.org/A000728/list）、

A000729（http://oeis.org/A000729/list）、

A000730（http://oeis.org/A000730/list）、

A000731（http://oeis.org/A000731/list）、

A010817（http://oeis.org/A010817/list）、

A010818（http://oeis.org/A010818/list）

と比較し、答え合わせしてみる。

# 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 phi_k(k, 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)
  }
  # k乗
  power(ary, k, n)
end

def A010815(n)
  phi_k(1, n)
end
ary = A010815(92)

# OEIS A010815のデータ
ary0 =
[1,-1,-1,0,0,1,0,1,0,0,0,0,-1,0,0,-1,0,0,0,0,0,0,
 1,0,0,0,1,0,0,0,0,0,0,0,0,-1,0,0,0,0,-1,0,0,0,0,0,
 0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
 -1,0,0,0,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1]
# 一致の確認
p ary == ary0

def A002107(n)
  phi_k(2, n)
end
ary = A002107(77)

# OEIS A002107のデータ
ary0 =
[1,-2,-1,2,1,2,-2,0,-2,-2,1,0,0,2,3,-2,2,0,0,-2,
 -2,0,0,-2,-1,0,2,2,-2,2,1,2,0,2,-2,-2,2,0,-2,0,-4,
 0,0,0,1,-2,0,0,2,0,2,2,1,-2,0,2,2,0,0,-2,0,-2,0,
 -2,2,0,-4,0,0,-2,-1,2,0,2,0,0,0,-2]
# 一致の確認
p ary == ary0

def A010816(n)
  phi_k(3, n)
end
ary = A010816(98)

# OEIS A010816のデータ
ary0 =
[1,-3,0,5,0,0,-7,0,0,0,9,0,0,0,0,-11,0,0,0,0,0,13,
 0,0,0,0,0,0,-15,0,0,0,0,0,0,0,17,0,0,0,0,0,0,0,0,
 -19,0,0,0,0,0,0,0,0,0,21,0,0,0,0,0,0,0,0,0,0,-23,
 0,0,0,0,0,0,0,0,0,0,0,25,0,0,0,0,0,0,0,0,0,0,0,0,
 -27,0,0,0,0,0,0,0]
# 一致の確認
p ary == ary0

def A000727(n)
  phi_k(4, n)
end
ary = A000727(81)

# OEIS A000727のデータ
ary0 =
[1,-4,2,8,-5,-4,-10,8,9,0,14,-16,-10,-4,0,-8,14,
 20,2,0,-11,20,-32,-16,0,-4,14,8,-9,20,26,0,2,-28,
 0,-16,16,-28,-22,0,14,16,0,40,0,-28,26,32,-17,0,
 -32,-16,-22,0,-10,32,-34,-8,14,0,45,-4,38,8,0,0,
 -34,-8,38,0,-22,-56,2,-28,0,0,-10,20,64,-40,-20,
 44]
# 一致の確認
p ary == ary0

def A000728(n)
  phi_k(5, n)
end
ary = A000728(56)

# OEIS A000728のデータ
ary0 =
[1,-5,5,10,-15,-6,-5,25,15,-20,9,-45,-5,25,20,10,
 15,20,-50,-35,-30,55,-50,15,80,1,50,-35,-45,-15,5,
 -50,-25,-55,85,51,50,10,-40,65,10,-10,-115,50,
 -115,-100,85,80,-30,5,20,45,70,65,45,-55,-100]
# 一致の確認
p ary == ary0

def A000729(n)
  phi_k(6, n)
end
ary = A000729(68)

# OEIS A000729のデータ
ary0 =
[1,-6,9,10,-30,0,11,42,0,-70,18,-54,49,90,0,-22,
 -60,0,-110,0,81,180,-78,0,130,-198,0,-182,-30,90,
 121,84,0,0,210,0,-252,-102,-270,170,0,0,-69,330,0,
 -38,420,0,-190,-390,0,-108,0,0,0,-300,99,442,210,
 0,418,-294,0,0,-510,378,-540,138,0]
# 一致の確認
p ary == ary0

def A000730(n)
  phi_k(7, n)
end
ary = A000730(45)

# OEIS A000730のデータ
ary0 =
[1,-7,14,7,-49,21,35,41,-49,-133,98,-21,126,112,
 -176,-105,-126,140,-35,147,259,98,-420,-224,238,
 -455,273,-14,322,406,-35,-7,-637,-196,245,-181,
 -574,462,147,924,217,-329,-140,-7,-371,-777]
# 一致の確認
p ary == ary0

def A000731(n)
  phi_k(8, n)
end
ary = A000731(60)

# OEIS A000731のデータ
ary0 =
[1,-8,20,0,-70,64,56,0,-125,-160,308,0,110,0,-520,
 0,57,560,0,0,182,-512,-880,0,1190,-448,884,0,0,0,
 -1400,0,-1330,1000,1820,0,-646,1280,0,0,-1331,
 -2464,380,0,1120,0,2576,0,0,-880,1748,0,-3850,0,
 -3400,0,2703,4160,-2500,0,3458]
# 一致の確認
p ary == ary0

def A010817(n)
  phi_k(9, n)
end
ary = A010817(40)

# OEIS A010817のデータ
ary0 =
[1,-9,27,-12,-90,135,54,-99,-189,-85,657,-162,
 -135,-171,-810,702,495,837,-673,-900,243,-1053,
 -297,1566,2700,-1764,81,-1188,-1377,270,-2043,
 3321,-756,3726,3015,-4563,-3348,504,-351,-1350,
 -468]
# 一致の確認
p ary == ary0

def A010818(n)
  phi_k(10, n)
end
ary = A010818(45)

# OEIS A010818のデータ
ary0 =
[1,-10,35,-30,-105,238,0,-260,-165,140,1054,-770,
 -595,0,-715,2162,455,0,-2380,-1820,2401,-680,1495,
 3080,1615,-6958,-1925,0,0,5100,-1442,8330,-5355,
 1330,0,-16790,0,8190,8265,0,1918,0,8415,-10230,
 -7140,-9362]
# 一致の確認
p ary == ary0

出力結果

true
true
true
true
true
true
true
true
true
true

2015年9月27日日曜日

150927(2)

Ruby

オイラー関数のベキ(2)

φ(x)^n の最初の方を出力してみた。
また、「ラマヌジャンの遺した関数」の3.8 において、
φ(x)^n の最初の方の，たとえば500個の，係数の中の0の数を求めることは，
単純なコンピュータ・プログラムの問題である
と書いてあったので、この個数も出力してみた。

# -*- coding: cp932 -*-

# 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 phi_k(k, 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)
  }
  # k乗
  power(ary, k, n)
end

(1..35).each{|k| p [k, phi_k(k, 50)]}
puts "500次までに現れる0の個数"
(1..35).each{|k| p [k, phi_k(k, 500).count(0)]}

出力結果

[1, [1, -1, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
[2, [1, -2, -1, 2, 1, 2, -2, 0, -2, -2, 1, 0, 0, 2, 3, -2, 2, 0, 0, -2, -2, 0, 0, -2, -1, 0, 2, 2, -2, 2, 1, 2, 0, 2, -2, -2, 2, 0, -2, 0, -4, 0, 0, 0, 1, -2, 0, 0, 2, 0, 2]]
[3, [1, -3, 0, 5, 0, 0, -7, 0, 0, 0, 9, 0, 0, 0, 0, -11, 0, 0, 0, 0, 0, 13, 0, 0, 0, 0, 0, 0, -15, 0, 0, 0, 0, 0, 0, 0, 17, 0, 0, 0, 0, 0, 0, 0, 0, -19, 0, 0, 0, 0, 0]]
[4, [1, -4, 2, 8, -5, -4, -10, 8, 9, 0, 14, -16, -10, -4, 0, -8, 14, 20, 2, 0, -11, 20, -32, -16, 0, -4, 14, 8, -9, 20, 26, 0, 2, -28, 0, -16, 16, -28, -22, 0, 14, 16, 0, 40, 0, -28, 26, 32, -17, 0, -32]]
[5, [1, -5, 5, 10, -15, -6, -5, 25, 15, -20, 9, -45, -5, 25, 20, 10, 15, 20, -50, -35, -30, 55, -50, 15, 80, 1, 50, -35, -45, -15, 5, -50, -25, -55, 85, 51, 50, 10, -40, 65, 10, -10, -115, 50, -115, -100, 85, 80, -30, 5, 20]]
[6, [1, -6, 9, 10, -30, 0, 11, 42, 0, -70, 18, -54, 49, 90, 0, -22, -60, 0, -110, 0, 81, 180, -78, 0, 130, -198, 0, -182, -30, 90, 121, 84, 0, 0, 210, 0, -252, -102, -270, 170, 0, 0, -69, 330, 0, -38, 420, 0, -190, -390, 0]]
[7, [1, -7, 14, 7, -49, 21, 35, 41, -49, -133, 98, -21, 126, 112, -176, -105, -126, 140, -35, 147, 259, 98, -420, -224, 238, -455, 273, -14, 322, 406, -35, -7, -637, -196, 245, -181, -574, 462, 147, 924, 217, -329, -140, -7, -371, -777, 588, -560, -196, -489, 1246]]
[8, [1, -8, 20, 0, -70, 64, 56, 0, -125, -160, 308, 0, 110, 0, -520, 0, 57, 560, 0, 0, 182, -512, -880, 0, 1190, -448, 884, 0, 0, 0, -1400, 0, -1330, 1000, 1820, 0, -646, 1280, 0, 0, -1331, -2464, 380, 0, 1120, 0, 2576, 0, 0, -880, 1748]]
[9, [1, -9, 27, -12, -90, 135, 54, -99, -189, -85, 657, -162, -135, -171, -810, 702, 495, 837, -673, -900, 243, -1053, -297, 1566, 2700, -1764, 81, -1188, -1377, 270, -2043, 3321, -756, 3726, 3015, -4563, -3348, 504, -351, -1350, -468, -891, 7074, 1611, 2700, -2423, -1512, -3267, -5265, -1800, 3510]]
[10, [1, -10, 35, -30, -105, 238, 0, -260, -165, 140, 1054, -770, -595, 0, -715, 2162, 455, 0, -2380, -1820, 2401, -680, 1495, 3080, 1615, -6958, -1925, 0, 0, 5100, -1442, 8330, -5355, 1330, 0, -16790, 0, 8190, 8265, 0, 1918, 0, 8415, -10230, -7140, -9362, -7315, 10010, 0, 14260, 14641]]
[11, [1, -11, 44, -55, -110, 374, -143, -462, 55, 495, 1287, -2069, -902, 1210, -275, 3795, -1507, -2431, -3575, -385, 8690, -1661, 1143, 1265, -4290, -12716, 2299, 11440, 3905, 8635, -10472, 6105, -20548, -1540, 8690, -24904, 29634, 25003, 8470, -23320, -18183, -4741, 2420, -19195, 2200, 18271, 5643, 52382, -14520, 990, -12287]]
[12, [1, -12, 54, -88, -99, 540, -418, -648, 594, 836, 1056, -4104, -209, 4104, -594, 4256, -6480, -4752, -298, 5016, 17226, -12100, -5346, -1296, -9063, -7128, 19494, 29160, -10032, -7668, -34738, 8712, -22572, 21812, 49248, -46872, 67562, 2508, -47520, -76912, -25191, 67716, 32076, 7128, 29754, 36784, -51072, 45144, -122398, -53460, 11286]]
[13, [1, -13, 65, -130, -65, 728, -871, -715, 1560, 845, 78, -6513, 2730, 8605, -4355, 2483, -13299, -2275, 11440, 10010, 19734, -41834, -11375, 12870, -2730, 14911, 33201, 25155, -70070, -36595, -28925, 64389, 13650, 52780, 72215, -173693, 87867, -62920, -111865, -25870, 89908, 284505, -37895, -88660, -59995, -57759, -172081, 117390, -121550, 61490, 301041]]
[14, [1, -14, 77, -182, 0, 924, -1547, -506, 3003, 0, -1729, -8372, 9177, 13090, -15625, 0, -17017, 10556, 30107, 0, 7084, -89206, 11571, 69160, 0, 27132, 0, -19096, -153502, 0, 93093, 165242, 0, -38962, 0, -420838, 257439, 0, -76153, 218750, 168245, 397320, -638066, -399126, 0, 67158, 92092, 530348, 0, 0, 424879]]
[15, [1, -15, 90, -245, 105, 1107, -2485, 195, 4860, -2420, -3990, -8190, 19695, 13755, -38475, 3990, -9750, 34020, 43015, -46605, -13860, -127385, 106485, 165240, -79275, -16380, -92340, -35840, -151995, 188550, 315783, 90090, -271215, -307485, 20475, -505440, 915385, 209340, -284130, 337645, -294225, 269325, -1707970, -70305, 1297620, 574210, 492765, 251370, -847245, -1102725, 438129]]
[16, [1, -16, 104, -320, 260, 1248, -3712, 1664, 6890, -7280, -5568, -4160, 33176, 4640, -74240, 29824, 14035, 54288, 27040, -142720, 1508, -110240, 289536, 222720, -380770, -83200, -123904, 142912, 7640, 408000, 386048, -530816, -755943, -294320, 716560, -194688, 1843200, -250432, -1688960, 324480, -988858, 1025440, -2283008, 1963520, 3857360, -1118240, -965120, -2204800, -2004730, -976720, 3450304]]
[17, [1, -17, 119, -408, 476, 1309, -5236, 4233, 8602, -15470, -4250, 5236, 45815, -21182, -117776, 101065, 46767, 36685, -36771, -267036, 143514, -18241, 486285, 81753, -1007250, 104006, 165767, 579292, 78829, 187510, 60214, -1706885, -616539, 956879, 2210985, -450109, 1609730, -2449615, -4158149, 2521015, -86275, 3856314, -3279895, 3835982, 4922435, -9624346, -2911454, -2670955, 2759372, 4788186, 8763296]]
[18, [1, -18, 135, -510, 765, 1242, -7038, 8280, 9180, -27710, 3519, 20196, 50370, -68850, -153765, 244782, 52785, -71010, -130525, -343620, 517293, 54978, 498780, -390150, -1835865, 1161270, 896751, 793730, -633420, -906660, -75582, -2589984, 1523745, 3589380, 2472615, -3740850, -767039, -4649670, -4222800, 11166210, 1718937, 4728294, -10403235, 2349450, 5331285, -23826622, 6799950, 7601040, 14132100, 7375230, -989604]]
[19, [1, -19, 152, -627, 1140, 988, -9063, 14212, 7410, -44270, 22781, 38114, 36176, -137256, -154850, 480605, -46493, -316065, -153406, -254525, 1156948, -184927, 88483, -1051042, -2381650, 3838874, 1417039, -542146, -2649911, -2171510, 2131306, -2282489, 5694509, 5022973, -2849050, -10735551, -1364789, -88825, 1790369, 24936550, -6189326, -6533150, -25397471, 6233254, 19679725, -38263549, 39975430, 19318117, 5441714, -16667750, -38728042]]
[20, [1, -20, 170, -760, 1615, 476, -11210, 22440, 1615, -64600, 60002, 51680, -9520, -213180, -83980, 803528, -379525, -692360, 119700, 80920, 1899830, -1235360, -755990, -1200040, -1981435, 8388956, -361760, -5068440, -4585935, -788120, 9421910, -2949160, 8315255, 768740, -16070560, -13715512, 10200340, 16428540, 5608780, 28549800, -40127031, -28504580, -23992440, 43014080, 60146875, -77059100, 80840440, -5335960, -56174545, -45601520, -56122066]]
[21, [1, -21, 189, -910, 2205, -378, -13321, 33345, -10395, -86870, 122703, 46683, -98287, -264915, 96390, 1163064, -1113588, -1066527, 1042055, 536025, 2287467, -3603805, -1391733, 478170, -562555, 13742379, -7889805, -12745348, -1009470, 6926850, 21064883, -13691664, 4004343, -8355270, -30343950, 5444712, 39441969, 30902949, -22027005, 4899895, -91000161, -20752011, 44496424, 116926740, 78003135, -226108554, 116695530, -60695838, -129044790, 45585540, 12284811]]
[22, [1, -22, 209, -1078, 2926, -1672, -15169, 47234, -31350, -107426, 218680, -266, -234707, -237006, 405878, 1444806, -2415413, -1091398, 3018169, 523050, 1618309, -7344304, -134905, 5365866, 5852, 17297588, -24278276, -18767364, 17865419, 19729952, 27154743, -49645442, -3483403, -4925446, -31064495, 61081922, 61961867, -55594, -108218187, -18341400, -88377696, 77121660, 185417067, 120639398, -58813391, -545883338, 294675997, 2967910, -128962680, 344990030, 286748]]
[23, [1, -23, 230, -1265, 3795, -3519, -16445, 64285, -64515, -120175, 354706, -123763, -407560, -48530, 817190, 1464341, -4376693, -135355, 6303955, -1282710, -682088, -11372603, 5678585, 13479425, -5451115, 16579596, -48805655, -11515065, 61570080, 21234520, 7731427, -119019296, 17214120, 46163645, -22347260, 134763417, 9991982, -115146395, -208431980, 72814665, 62601078, 221269499, 248467505, -164778900, -396525635, -820862617, 1040502519, 303756860, -310026775, 615008615, -577640492]]
[24, [1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920, 534612, -370944, -577738, 401856, 1217160, 987136, -6905934, 2727432, 10661420, -7109760, -4219488, -12830688, 18643272, 21288960, -25499225, 13865712, -73279080, 24647168, 128406630, -29211840, -52843168, -196706304, 134722224, 165742416, -80873520, 167282496, -182213314, -255874080, -145589976, 408038400, 308120442, 101267712, -17125708, -786948864, -548895690, -447438528, 2687348496, 248758272, -1696965207, 611981400, -1740295368]]
[25, [1, -25, 275, -1700, 6050, -9405, -15550, 107525, -182875, -81675, 756655, -801550, -662975, 1220175, 1361350, -209440, -9601900, 8608900, 14889050, -19948500, -6262465, -7057550, 38788925, 19716425, -69119875, 23579969, -82427400, 98068850, 191984400, -192983175, -128640655, -199535875, 424794500, 284736850, -398127725, 141193525, -454458725, -190957250, 306041450, 867327675, 248309215, -781592900, -498119600, -1140582625, 331248600, 951687880, 4441237625, -2000633400, -5357206250, 1875749425, -1810005816]]
[26, [1, -26, 299, -1950, 7475, -13754, -12220, 132756, -276575, 0, 1010100, -1486030, -519961, 2486300, 829725, -2215486, -11643060, 18523050, 16317925, -42861650, 0, 11010090, 59644221, -5743400, -138219900, 79631474, -64447293, 190305726, 197368275, -523231800, -99201921, 0, 873014519, 160509700, -1222259675, 305612684, -511748120, 355081636, 1089106200, 759784300, -698434100, -2494548030, 0, 0, 2468263525, 2158297050, 3767134280, -8192279446, -9594939900, 9772308650, 1157638599]]
[27, [1, -27, 324, -2223, 9126, -19278, -5967, 159030, -399087, 151593, 1270971, -2500875, 74970, 4203522, -1004157, -4796037, -11750778, 32885190, 10452375, -77533092, 27104868, 43070625, 63798840, -69960267, -215939061, 236414349, -37046646, 237487433, 85921371, -1008703449, 286178139, 474257484, 1224628470, -608265126, -2606289075, 1473758712, 10716732, 1192052160, 1464312006, -1064732643, -2341040562, -3747291822, 3750796530, 2904349500, 3309301413, -30203550, -1508652000, -15903743310, -7320117141, 31078702161, 1804201776]]
[28, [1, -28, 350, -2520, 11025, -26180, 4158, 184600, -554400, 401100, 1496964, -3920280, 1444625, 6224400, -4972350, -7121296, -8308965, 50796900, -8971200, -121968000, 94011435, 80598288, 20282500, -175228200, -254651775, 554394204, -88470242, 133725200, -133867445, -1515976700, 1369368000, 1033213272, 879369575, -2314141200, -3980458300, 5107622936, 781936407, 910872900, 222868800, -5087226200, -2187994564, -1973664000, 11996468950, 3892051800, -3360359625, -6947241084, -7654070088, -15374268000, 9869832000, 62945277200, -22144800636]]
[29, [1, -29, 377, -2842, 13195, -34684, 19285, 206973, -745706, 782275, 1621564, -5803161, 4026360, 8149841, -12056025, -7428263, 254504, 69194580, -49156653, -167517050, 224634319, 94868280, -112333182, -288914501, -172722550, 1061590530, -420678727, -212254364, -271663242, -1795846200, 3413927270, 967351840, -978399506, -4524424502, -3825362450, 12725565255, -657413528, -3315667498, -2494031813, -9144119775, 5035400442, 3578388387, 20759047670, -5893341662, -22777727375, -8797952807, -1616838769, 3681888314, 36274904939, 80118476175, -103728485744]]
[30, [1, -30, 405, -3190, 15660, -45036, 40745, 222750, -974835, 1334580, 1547469, -8174520, 8380245, 9200250, -23243355, -2643380, 14704740, 82050570, -116275500, -195804810, 442809990, 25147930, -371898000, -313802910, 125394405, 1688931000, -1364323095, -737497840, 158838945, -1653918750, 6309965146, -1076120370, -4802530500, -5257026620, -47508120, 24290828532, -10318690855, -13910780610, -1364146515, -7675874280, 24023587500, 5245892100, 18578165165, -33913100250, -47491260255, 23066963660, 27308816280, 29296187730, 36536324105, 44794138050, -245483273862]]
[31, [1, -31, 434, -3565, 18445, -57505, 70091, 227447, -1241550, 2102730, 1139498, -11000164, 15185009, 8060465, -39266925, 11975548, 33735905, 79961555, -212042635, -176681400, 762467041, -231771190, -762218948, -59474275, 687626655, 2193123086, -3317871844, -996975686, 1948141060, -1380351105, 9298132746, -6788705842, -9608214462, -500684100, 8195740630, 35945695112, -37464288578, -27703334946, 18233392140, 3452352820, 50110166442, -18648267155, -5426442244, -68394536935, -47776054550, 125566716471, 51003097355, 3739338358, -34904617090, -37022572120, -360301225363]]
[32, [1, -32, 464, -3968, 21576, -72384, 109120, 215296, -1542684, 3135712, 217248, -14153856, 25215616, 2704192, -60182656, 43083520, 52111434, 50631680, -328746320, -68928128, 1172526144, -825260672, -1202344640, 768450816, 1395312728, 2106191584, -6556788640, 100803072, 5930017280, -2446579520, 11062835584, -17191741952, -11762712973, 14959546528, 16880156208, 38804082816, -90268764128, -27924752384, 79537896320, 11907174400, 59476916292, -96945887328, -36151163360, -58860723712, 11806125312, 301046180224, -41929489920, -152538292480, -142811793922, -40571481632, -282758291200]]
[33, [1, -33, 495, -4400, 25080, -89991, 159896, 179025, -1871100, 4485140, -1452132, -17376120, 39295135, -9791100, -84751920, 99240515, 58170882, -19058490, -443044250, 173982600, 1617603372, -1940174742, -1456327455, 2535117750, 1782559240, 783282456, -10878123663, 4453210080, 12274486950, -8274669480, 10738181916, -30949593906, -4649643405, 44427298525, 12632491905, 22022095887, -166221342056, 21965133795, 201364886250, -36031327805, 13370409459, -232741171476, 2257921270, 73994101500, 103475478645, 461792586948, -436927039032, -414115099695, -36601699750, 242404463730, 31750490970]]
[34, [1, -34, 527, -4862, 28985, -110670, 224774, 109616, -2214454, 6202790, -4160563, -20224152, 58212420, -33421286, -109592285, 190464702, 30766260, -138537388, -509479358, 590424450, 1978075528, -3762986788, -1052300952, 5565993766, 814913275, -2348596580, -15108004380, 14570701530, 19246124150, -24572149690, 10398722367, -42736429098, 19198424003, 83788921020, -29646853725, -16183052760, -240028781751, 173007827374, 367539275460, -257643143340, -99099230160, -344243246886, 252100115707, 356338645124, 24177692825, 446603464544, -1249233123069, -441921135810, 765164623614, 771325632000, 152306691108]]
[35, [1, -35, 560, -5355, 33320, -134792, 306425, -3960, -2553740, 8337140, -8264431, -22012620, 82596080, -73399200, -128071540, 327108117, -64156470, -305606350, -455573580, 1186922730, 2056498976, -6416740355, 788588570, 9925006980, -3311307160, -7223330464, -16617037990, 32704289240, 21589314305, -58378128620, 19258974280, -43471191470, 61914761755, 115717397495, -142117772440, -48251561329, -258673133355, 462387945810, 485243715915, -824196278360, -158905982564, -216715690380, 829306675115, 595583270100, -665090670370, 246500217957, -2246704399660, 504524131615, 2608529471870, 601239063870, -912722944252]]
500次までに現れる0の個数
[1, 464]
[2, 243]
[3, 469]
[4, 158]
[5, 0]
[6, 212]
[7, 0]
[8, 250]
[9, 0]
[10, 151]
[11, 0]
[12, 0]
[13, 0]
[14, 172]
[15, 2]
[16, 0]
[17, 0]
[18, 0]
[19, 0]
[20, 0]
[21, 0]
[22, 0]
[23, 0]
[24, 0]
[25, 0]
[26, 80]
[27, 0]
[28, 0]
[29, 0]
[30, 0]
[31, 0]
[32, 0]
[33, 0]
[34, 0]
[35, 0]

150927

Ruby

k角数定理？(1)

「組合せ論プロムナード」の第5講において以下の内容が載っていた。

「ヤコビの3重積公式」の変数を置き換えると、

「ガウスの4角数定理」および「オイラーの5角数定理」が導ける。

さて、同様の変換を行うことで、

以下のようなk角数定理？のようなものが導ける。

Π(1 - q^((k - 2)n))(1 - q^((k - 2)n - 1))(1 - q^((k - 2)n - k + 3))

= Σ(-1)^m q^(m((k - 2) * m + k - 4) / 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 lhs(k, n)
  ary = [1]
  # 無限積のうち必要なところだけ取り出す
  # 1 - q^((k - 2)n)
  (k - 2).step(n, k - 2){|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary, n)
  }
  # 1 - q^((k - 2)n - 1)
  (k - 3).step(n, k - 2){|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary, n)
  }
  # 1 - q^((k - 2)n - k + 3)
  1.step(n, k - 2){|i|
    b_ary = Array.new(i + 1, 0)
    b_ary[0], b_ary[-1] = 1, -1
    ary = mul(ary, b_ary, n)
  }
  s = ary.size
  # 形を整えておく
  ary += [0] * (n + 1 - s) if s < n + 1
  ary
end

# k角数
def k_number(k, n)
  n * ((k - 2) * n + k - 4) / 2
end

def rhs(k, n)
  ary = Array.new(n + 1, 0)
  i = 0
  j = k_number(k, i)
  while j <= n
    ary[j] = (-1) ** i
    i += 1
    j = k_number(k, i)
  end
  i = -1
  j = k_number(k, i)
  while j <= n
    ary[j] = (-1) ** i
    i -= 1
    j = k_number(k, i)
  end
  ary
end

n = 150
(5..10).each{|k|
  # lhs = rhs となることを確認する（一致したらT、しないならF）
  p (0..n).map{|i| lhs(k, i) == rhs(k, i) ? 'T' : 'F'}
  p [k, lhs(k, n)]
}

出力結果

["T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"]
[5, [1, -1, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]]
["T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"]
[6, [1, -1, 0, -1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
["T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"]
[7, [1, -1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]]
["T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"]
[8, [1, -1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
["T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"]
[9, [1, -1, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
["T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T", "T"]
[10, [1, -1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

2015年9月24日木曜日

150924

Ruby

オイラー関数のベキ(1)

一昨日ダイソンが見つけた公式を確認したが、クラインが見つけた公式も確認しておく。

# m次まで計算
def K(m, n)
  # i<=j<=kとし、n = [-i, k].maxとする
  ary = Array.new(m + 1, 0)
  (-n).upto(0){|i|
    n.downto(0){|k|
      j = 1 - i - k
      if i <= j && j <= k
        m0 = i * j + i * k + k * j
        if -m0 <= m
          # i = j = kとはならないので、以下の条件はちょうど2つ同じものがある場合を表している
          if i == j || i == k || j == k
            # i,j,kの大小関係を置かないときは3通り
            s = 0.5 * (2 + 9 * (3 * i * j * k - m0))
          else
            # i,j,kの大小関係を置かないときは6通り
            s = 2 + 9 * (3 * i * j * k - m0)
          end
          p [-m0, [i, j, k], s]
          ary[-m0] += s
        end
      end
    }
  }
  ary.map(&:to_i)
end
ary = K(200, 200)
p ary

# 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
  mul(ary, power(ary, n - 1, m), m)
end

def phi_8(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)
  }
  # 8乗
  power(ary, 8, n)
end
ary0 = phi_8(200)

# 一致の確認
p ary == ary0

出力結果

[200, [-16, 8, 9], -29302]

[192, [-15, 4, 12], -17710]

[185, [-15, 5, 11], -20608]

[180, [-15, 6, 10], -22678]

[177, [-15, 7, 9], -23920]

[176, [-15, 8, 8], -12167.0]

[196, [-14, 1, 14], -3526]

[184, [-14, 2, 13], -8170]

[174, [-14, 3, 12], -12040]

[166, [-14, 4, 11], -15136]

[160, [-14, 5, 10], -17458]

[156, [-14, 6, 9], -19006]

[154, [-14, 7, 8], -19780]

[197, [-13, -1, 15], 7040]

[182, [-13, 0, 14], 1640]

[169, [-13, 1, 13], -3040]

[158, [-13, 2, 12], -7000]

[149, [-13, 3, 11], -10240]

[142, [-13, 4, 10], -12760]

[137, [-13, 5, 9], -14560]

[134, [-13, 6, 8], -15640]

[133, [-13, 7, 7], -8000.0]

[186, [-12, -2, 15], 11396]

[170, [-12, -1, 14], 6068]

[156, [-12, 0, 13], 1406]

[144, [-12, 1, 12], -2590]

[134, [-12, 2, 11], -5920]

[126, [-12, 3, 10], -8584]

[120, [-12, 4, 9], -10582]

[116, [-12, 5, 8], -11914]

[114, [-12, 6, 7], -12580]

[196, [-11, -4, 16], 20774]

[177, [-11, -3, 15], 14960]

[160, [-11, -2, 14], 9758]

[145, [-11, -1, 13], 5168]

[132, [-11, 0, 12], 1190]

[121, [-11, 1, 11], -2176]

[112, [-11, 2, 10], -4930]

[105, [-11, 3, 9], -7072]

[100, [-11, 4, 8], -8602]

[97, [-11, 5, 7], -9520]

[96, [-11, 6, 6], -4913.0]

[190, [-10, -5, 16], 23312]

[170, [-10, -4, 15], 17732]

[152, [-10, -3, 14], 12710]

[136, [-10, -2, 13], 8246]

[122, [-10, -1, 12], 4340]

[110, [-10, 0, 11], 992]

[100, [-10, 1, 10], -1798]

[92, [-10, 2, 9], -4030]

[86, [-10, 3, 8], -5704]

[82, [-10, 4, 7], -6820]

[80, [-10, 5, 6], -7378]

[186, [-9, -6, 16], 25004]

[165, [-9, -5, 15], 19712]

[146, [-9, -4, 14], 14924]

[129, [-9, -3, 13], 10640]

[114, [-9, -2, 12], 6860]

[101, [-9, -1, 11], 3584]

[90, [-9, 0, 10], 812]

[81, [-9, 1, 9], -1456]

[74, [-9, 2, 8], -3220]

[69, [-9, 3, 7], -4480]

[66, [-9, 4, 6], -5236]

[65, [-9, 5, 5], -2744.0]

[184, [-8, -7, 16], 25850]

[162, [-8, -6, 15], 20900]

[142, [-8, -5, 14], 16400]

[124, [-8, -4, 13], 12350]

[108, [-8, -3, 12], 8750]

[94, [-8, -2, 11], 5600]

[82, [-8, -1, 10], 2900]

[72, [-8, 0, 9], 650]

[64, [-8, 1, 8], -1150]

[58, [-8, 2, 7], -2500]

[54, [-8, 3, 6], -3400]

[52, [-8, 4, 5], -3850]

[161, [-7, -7, 15], 10648.0]

[140, [-7, -6, 14], 17138]

[121, [-7, -5, 13], 13376]

[104, [-7, -4, 12], 10010]

[89, [-7, -3, 11], 7040]

[76, [-7, -2, 10], 4466]

[65, [-7, -1, 9], 2288]

[56, [-7, 0, 8], 506]

[49, [-7, 1, 7], -880]

[44, [-7, 2, 6], -1870]

[41, [-7, 3, 5], -2464]

[40, [-7, 4, 4], -1331.0]

[120, [-6, -6, 13], 6859.0]

[102, [-6, -5, 12], 10640]

[86, [-6, -4, 11], 7904]

[72, [-6, -3, 10], 5510]

[60, [-6, -2, 9], 3458]

[50, [-6, -1, 8], 1748]

[42, [-6, 0, 7], 380]

[36, [-6, 1, 6], -646]

[32, [-6, 2, 5], -1330]

[30, [-6, 3, 4], -1672]

[85, [-5, -5, 11], 4096.0]

[70, [-5, -4, 10], 6032]

[57, [-5, -3, 9], 4160]

[46, [-5, -2, 8], 2576]

[37, [-5, -1, 7], 1280]

[30, [-5, 0, 6], 272]

[25, [-5, 1, 5], -448]

[22, [-5, 2, 4], -880]

[21, [-5, 3, 3], -512.0]

[56, [-4, -4, 9], 2197.0]

[44, [-4, -3, 8], 2990]

[34, [-4, -2, 7], 1820]

[26, [-4, -1, 6], 884]

[20, [-4, 0, 5], 182]

[16, [-4, 1, 4], -286]

[14, [-4, 2, 3], -520]

[33, [-3, -3, 7], 1000.0]

[24, [-3, -2, 6], 1190]

[17, [-3, -1, 5], 560]

[12, [-3, 0, 4], 110]

[9, [-3, 1, 3], -160]

[8, [-3, 2, 2], -125.0]

[16, [-2, -2, 5], 343.0]

[10, [-2, -1, 4], 308]

[6, [-2, 0, 3], 56]

[4, [-2, 1, 2], -70]

[5, [-1, -1, 3], 64.0]

[2, [-1, 0, 2], 20]

[1, [-1, 1, 1], -8.0]

[0, [0, 0, 1], 1.0]

[1, -8, 20, 0, -70, 64, 56, 0, -125, -160, 308, 0, 110, 0, -520, 0, 57, 560, 0, 0, 182, -512, -880, 0, 1190, -448, 884, 0, 0, 0, -1400, 0, -1330, 1000, 1820, 0, -646, 1280, 0, 0, -1331, -2464, 380, 0, 1120, 0, 2576, 0, 0, -880, 1748, 0, -3850, 0, -3400, 0, 2703, 4160, -2500, 0, 3458, 0, 0, 0, -1150, -456, -5236, 0, 0, -4480, 6032, 0, 6160, 0, -3220, 0, 4466, 0, 0, 0, -7378, -1456, -3920, 0, 0, 4096, 2200, 0, 0, 7040, 812, 0, -4030, 0, 5600, 0, -4913, -9520, 0, 0, -10400, 3584, 10640, 0, 10010, -7072, 0, 0, 8750, 0, 992, 0, -4930, 0, -5720, 0, -11914, 0, 0, 0, -3723, 11200, 4340, 0, 12350, 0, -8584, 0, 0, 10640, 0, 0, 1190, -8000, -21560, 0, 8246, -14560, 0, 0, 17138, 0, 3640, 0, -2590, 5168, 14924, 0, 0, -10240, 0, 0, 12710, 0, -19780, 0, -17600, 0, -7000, 0, -7700, 10648, 20900, 0, 0, 19712, -15136, 0, 0, -3040, 23800, 0, 0, 0, -12040, 0, -12167, -8960, 0, 0, -22678, 0, 1640, 0, 17680, -20608, 36400, 0, 0, 0, 23312, 0, -17710, 0, 0, 0, 17248, 7040, 0, 0, -29302]

true

2015年9月22日火曜日

150922(2)

Ruby

フリーマン・ダイソンによるτ関数に関する公式

「ラマヌジャンの遺した関数」を読んでいると、τ関数に関する公式が出てきたので
自分で確かめたくなった。
オンライン整数列大辞典の
A000594（http://oeis.org/A000594/list）
と比較し、答え合わせしてみる。

def ary(i, m)
  i.step(m, 5).to_a + (i - 5).step(-m, -5).to_a
end

def f(a, b, c, d, e)
  (a - b) * (a - c) * (a - d) * (a - e) * (b - c) * (b - d) * (b - e) * (c - d) * (c - e) * (d - e)
end

def tau_function(n)
  sum = 0
  # -mからmまで調べれば十分
  m = Math.sqrt(10 * n).to_i + 1
  ary(1, m).each{|a|
    ary(2, m).each{|b|
      ary(3, m).each{|c|
        ary(4, m).each{|d|
          ary(5, m).each{|e|
            if a + b + c + d + e == 0
              if a * a + b * b + c * c + d * d + e * e == 10 * n
                f = f(a, b, c, d, e)
                p [n, [a, b, c, d, e], f]
                sum += f
              end
            end
          }
        }
      }
    }
  }
  # 288 = 1!2!3!4!
  sum / 288
end

def A000594(n)
  (1..n).map{|i| tau_function(i)}
end
ary = A000594(28)

# OEIS A000594のデータ
ary0 =
[1,-24,252,-1472,4830,-6048,-16744,84480,-113643,
 -115920,534612,-370944,-577738,401856,1217160,
 987136,-6905934,2727432,10661420,-7109760,
 -4219488,-12830688,18643272,21288960,-25499225,
 13865712,-73279080,24647168]
# 一致の確認
p ary == ary0

出力結果
[1, [1, 2, -2, -1, 0], 288]
[2, [1, -3, 3, -1, 0], -6912]
[3, [1, -3, -2, 4, 0], 36288]
[3, [-4, 2, 3, -1, 0], 36288]
[4, [1, 2, 3, -1, -5], -64512]
[4, [1, -3, -2, -1, 5], -64512]
[4, [-4, 2, -2, 4, 0], -294912]
[5, [1, 2, 3, -6, 0], 36288]
[5, [1, 2, -2, 4, -5], 489888]
[5, [6, -3, -2, -1, 0], 36288]
[5, [-4, 2, -2, -1, 5], 489888]
[5, [-4, -3, 3, 4, 0], 338688]
[6, [1, -3, 3, 4, -5], -870912]
[6, [-4, -3, 3, -1, 5], -870912]
[7, [1, 2, -2, -6, 5], -2483712]
[7, [1, 2, -7, 4, 0], -266112]
[7, [6, 2, -2, -1, -5], -2483712]
[7, [-4, 2, 3, 4, -5], 338688]
[7, [-4, 7, -2, -1, 0], -266112]
[7, [-4, -3, -2, 4, 5], 338688]
[8, [1, 2, -7, -1, 5], 2239488]
[8, [1, 7, -2, -1, -5], 2239488]
[8, [1, -3, 3, -6, 5], 6386688]
[8, [6, 2, -2, -6, 0], 7077888]
[8, [6, -3, 3, -1, -5], 6386688]
[9, [1, 7, -2, -6, 0], -4953312]
[9, [1, -8, 3, 4, 0], 684288]
[9, [6, 2, -7, -1, 0], -4953312]
[9, [6, -3, 3, -6, 0], -17006112]
[9, [6, -3, -2, 4, -5], -3592512]
[9, [-4, 2, 3, -6, 5], -3592512]
[9, [-4, -3, 8, -1, 0], 684288]
[10, [1, 7, -7, -1, 0], 3161088]
[10, [1, -3, 8, -1, -5], -6918912]
[10, [1, -3, -7, 4, 5], -11354112]
[10, [1, -8, 3, -1, 5], -6918912]
[10, [-4, 7, 3, -1, -5], -11354112]
[11, [1, -3, 8, -6, 0], 13039488]
[11, [1, -8, -2, 4, 5], 12737088]
[11, [6, 2, 3, -6, -5], 6386688]
[11, [6, -3, -2, -6, 5], 6386688]
[11, [6, -3, -7, 4, 0], 36324288]
[11, [6, -8, 3, -1, 0], 13039488]
[11, [-4, 2, 8, -1, -5], 12737088]
[11, [-4, 2, -7, 4, 5], 9237888]
[11, [-4, 7, 3, -6, 0], 36324288]
[11, [-4, 7, -2, 4, -5], 9237888]
[11, [-4, -3, -2, 9, 0], -741312]
[11, [-9, 2, 3, 4, 0], -741312]
[12, [1, 7, 3, -6, -5], -22643712]
[12, [1, -3, -2, 9, -5], 6386688]
[12, [6, -3, -7, -1, 5], -22643712]
[12, [6, -8, -2, 4, 0], -37158912]
[12, [-4, 2, 8, -6, 0], -37158912]
[12, [-9, 2, 3, -1, 5], 6386688]
[13, [1, 2, 3, 4, -10], 288288]
[13, [1, 2, 8, -6, -5], 17978688]
[13, [1, -8, 8, -1, 0], -8128512]
[13, [6, 2, -7, 4, -5], -28540512]
[13, [6, -8, -2, -1, 5], 17978688]
[13, [-4, 2, -2, 9, -5], -14126112]
[13, [-4, 7, -2, -6, 5], -28540512]
[13, [-4, 7, -7, 4, 0], -95622912]
[13, [-4, -3, 8, 4, -5], -6918912]
[13, [-4, -3, -2, -1, 10], 288288]
[13, [-4, -8, 3, 4, 5], -6918912]
[13, [-9, 2, -2, 4, 5], -14126112]
[14, [1, 7, -7, 4, -5], 86220288]
[14, [1, -3, -7, 9, 0], -37158912]
[14, [-4, 7, -7, -1, 5], 86220288]
[14, [-4, -3, 3, 9, -5], 8805888]
[14, [-9, 7, 3, -1, 0], -37158912]
[14, [-9, -3, 3, 4, 5], 8805888]
[15, [1, -3, -2, -6, 10], -22643712]
[15, [1, -8, -2, 9, 0], 34898688]
[15, [6, 2, 3, -1, -10], -22643712]
[15, [6, 2, -7, -6, 5], 17791488]
[15, [6, 7, -2, -6, -5], 17791488]
[15, [6, -8, 3, 4, -5], 26345088]
[15, [-4, 2, -7, 9, 0], 118879488]
[15, [-4, -3, 8, -6, 5], 26345088]
[15, [-9, 2, 8, -1, 0], 34898688]
[15, [-9, 7, -2, 4, 0], 118879488]
[16, [1, 2, 3, -11, 5], -1677312]
[16, [1, 2, -7, 9, -5], -75866112]
[16, [1, 7, 3, -1, -10], 67212288]
[16, [1, 7, -7, -6, 5], -64576512]
[16, [1, -3, -7, -1, 10], 67212288]
[16, [6, 2, -2, 4, -10], 66060288]
[16, [6, 7, -7, -1, -5], -64576512]
[16, [11, -3, -2, -1, -5], -1677312]
[16, [-4, 2, -2, -6, 10], 66060288]
[16, [-4, -8, 8, 4, 0], 301989888]
[16, [-9, 7, -2, -1, 5], -75866112]
[17, [1, 2, 8, -1, -10], -48498912]
[17, [1, 7, -2, 4, -10], -183218112]
[17, [1, -8, 8, 4, -5], -305690112]
[17, [1, -8, -2, -1, 10], -48498912]
[17, [6, 2, 3, -11, 0], 14702688]
[17, [6, 7, -7, -6, 0], 50083488]
[17, [6, -3, 3, 4, -10], -46230912]
[17, [6, -3, 8, -6, -5], -28540512]
[17, [6, -8, 3, -6, 5], -28540512]
[17, [11, -3, -2, -6, 0], 14702688]
[17, [-4, 2, -7, -1, 10], -183218112]
[17, [-4, -3, 3, -6, 10], -46230912]
[17, [-4, -8, 3, 9, 0], -352864512]
[17, [-4, -8, 8, -1, 5], -305690112]
[17, [-9, 2, -2, 9, 0], -138311712]
[17, [-9, -3, 8, 4, 0], -352864512]
[18, [1, 7, 3, -11, 0], -33530112]
[18, [1, -8, 3, 9, -5], 325721088]
[18, [11, -3, -7, -1, 0], -33530112]
[18, [-4, -3, -7, 9, 5], -103514112]
[18, [-9, 7, 3, 4, -5], -103514112]
[18, [-9, -3, 3, 9, 0], 408146688]
[18, [-9, -3, 8, -1, 5], 325721088]
[19, [1, 2, 8, -11, 0], 21909888]
[19, [1, 2, -2, 9, -10], 148313088]
[19, [1, 2, -7, -6, 10], 78962688]
[19, [1, 2, -12, 4, 5], 3564288]
[19, [1, -3, 8, 4, -10], 502020288]
[19, [1, -8, 8, -6, 5], 339026688]
[19, [6, 2, -2, -11, 5], -85542912]
[19, [6, 7, -2, -1, -10], 78962688]
[19, [6, -8, 8, -1, -5], 339026688]
[19, [6, -8, -7, 4, 5], 7495488]
[19, [11, 2, -2, -6, -5], -85542912]
[19, [11, -8, -2, -1, 0], 21909888]
[19, [-4, 7, 3, 4, -10], 137225088]
[19, [-4, 7, 8, -6, -5], 7495488]
[19, [-4, 12, -2, -1, -5], 3564288]
[19, [-4, -3, -7, 4, 10], 137225088]
[19, [-4, -8, 3, -1, 10], 502020288]
[19, [-4, -8, -2, 9, 5], 382269888]
[19, [-9, 2, 8, 4, -5], 382269888]
[19, [-9, 2, -2, -1, 10], 148313088]
[20, [1, 7, -2, -11, 5], 282175488]
[20, [1, -3, 3, 9, -10], -525837312]
[20, [6, 2, -12, 4, 0], -37158912]
[20, [6, -3, 3, -11, 5], 78962688]
[20, [6, -3, -7, 9, -5], 166053888]
[20, [6, -8, 8, -6, 0], -346816512]
[20, [11, 2, -7, -1, -5], 282175488]
[20, [11, -3, 3, -6, -5], 78962688]
[20, [-4, 2, 8, 4, -10], -501645312]
[20, [-4, 12, -2, -6, 0], -37158912]
[20, [-4, -8, -2, 4, 10], -501645312]
[20, [-9, 2, 3, 9, -5], -312947712]
[20, [-9, 7, 3, -6, 5], 166053888]
[20, [-9, -3, 3, -1, 10], -525837312]
[20, [-9, -3, -2, 9, 5], -312947712]
[21, [1, 7, -12, 4, 0], 71705088]
[21, [1, 12, -2, -6, -5], 71251488]
[21, [1, -8, 3, -6, 10], -452656512]
[21, [6, 2, -12, -1, 5], 71251488]
[21, [6, 7, -2, -11, 0], -183218112]
[21, [6, -3, 8, -1, -10], -452656512]
[21, [6, -8, -2, 9, -5], -522510912]
[21, [11, 2, -7, -6, 0], -183218112]
[21, [-4, 2, 3, 9, -10], 407822688]
[21, [-4, 12, -7, -1, 0], 71705088]
[21, [-9, 2, 8, -6, 5], -522510912]
[21, [-9, -3, -2, 4, 10], 407822688]
[22, [1, 7, -12, -1, 5], -212838912]
[22, [1, 12, -7, -1, -5], -212838912]
[22, [1, -3, 8, -11, 5], -862912512]
[22, [1, -8, -7, 9, 5], -391053312]
[22, [1, -13, 3, 4, 5], -3290112]
[22, [11, -3, -7, 4, -5], -312947712]
[22, [11, -8, 3, -1, -5], -862912512]
[22, [-4, 7, 3, -11, 5], -312947712]
[22, [-4, 7, -7, 9, -5], -64576512]
[22, [-4, -3, 13, -1, -5], -3290112]
[22, [-9, 7, 8, -1, -5], -391053312]
[22, [-9, 7, -7, 4, 5], -64576512]
[23, [1, 2, -2, -11, 10], -305690112]
[23, [1, 2, -12, 9, 0], -46230912]
[23, [1, 12, -7, -6, 0], 106178688]
[23, [1, -8, -7, 4, 10], 471132288]
[23, [6, 7, 3, -6, -10], -238298112]
[23, [6, 7, -12, -1, 0], 106178688]
[23, [6, -3, 8, -11, 0], 810425088]
[23, [6, -3, -2, 9, -10], 485388288]
[23, [6, -3, -7, -6, 10], -238298112]
[23, [6, -3, -12, 4, 5], -44416512]
[23, [6, -8, -7, 9, 0], 449100288]
[23, [6, -13, 3, 4, 0], 29023488]
[23, [11, 2, -2, -1, -10], -305690112]
[23, [11, -8, 3, -6, 0], 810425088]
[23, [11, -8, -2, 4, -5], 968018688]
[23, [-4, 2, 8, -11, 5], 968018688]
[23, [-4, 7, 8, -1, -10], 471132288]
[23, [-4, 12, 3, -6, -5], -44416512]
[23, [-4, -3, 13, -6, 0], 29023488]
[23, [-9, 2, 3, -6, 10], 485388288]
[23, [-9, 7, 8, -6, 0], 449100288]
[23, [-9, 12, -2, -1, 0], -46230912]
[24, [1, -3, 3, -11, 10], 1109541888]
[24, [1, -3, 13, -6, -5], -66189312]
[24, [6, 2, 8, -6, -10], 891813888]
[24, [6, 7, 3, -11, -5], 325721088]
[24, [6, -8, -2, -6, 10], 891813888]
[24, [6, -13, 3, -1, 5], -66189312]
[24, [11, -3, 3, -1, -10], 1109541888]
[24, [11, -3, -7, -6, 5], 325721088]
[24, [-9, 2, -7, 9, 5], 804722688]
[24, [-9, 7, -2, 9, -5], 804722688]
[25, [1, 2, -12, -1, 10], 209297088]
[25, [1, 7, 8, -6, -10], -720431712]
[25, [1, 12, -2, -1, -10], 209297088]
[25, [1, -13, 8, 4, 0], -174650112]
[25, [6, 2, 8, -11, -5], -1210521312]
[25, [6, 7, -7, 4, -10], 137225088]
[25, [6, -8, -7, -1, 10], -720431712]
[25, [6, -13, -2, 4, 5], 42977088]
[25, [11, 2, -2, -11, 0], 583041888]
[25, [11, -3, -2, 4, -10], -880955712]
[25, [11, -8, -2, -6, 5], -1210521312]
[25, [11, -8, -7, 4, 0], -778643712]
[25, [-4, 2, 3, -11, 10], -880955712]
[25, [-4, 2, 13, -6, -5], 42977088]
[25, [-4, 7, 8, -11, 0], -778643712]
[25, [-4, 7, -2, 9, -10], -877960512]
[25, [-4, 7, -7, -6, 10], 137225088]
[25, [-4, 7, -12, 4, 5], 196466688]
[25, [-4, 12, -7, 4, -5], 196466688]
[25, [-4, -3, -2, 14, -5], 1116288]
[25, [-4, -8, 8, 9, -5], 92671488]
[25, [-4, -8, 13, -1, 0], -174650112]
[25, [-9, 2, -7, 4, 10], -877960512]
[25, [-9, -8, 8, 4, 5], 92671488]
[25, [-14, 2, 3, 4, 5], 1116288]
[26, [1, 7, 8, -11, -5], 968018688]
[26, [1, -3, -12, 9, 5], 2053029888]
[26, [1, -8, 13, -1, -5], 576108288]
[26, [1, -13, 3, 9, 0], 166053888]
[26, [1, -13, 8, -1, 5], 576108288]
[26, [11, -3, 3, -11, 0], -1803174912]
[26, [11, -8, -7, -1, 5], 968018688]
[26, [-9, 7, -7, 9, 0], -1024192512]
[26, [-9, 12, 3, -1, -5], 2053029888]
[26, [-9, -3, 8, 9, -5], -352864512]
[26, [-9, -3, 13, -1, 0], 166053888]
[26, [-9, -8, 3, 9, 5], -352864512]
[27, [1, 12, -2, -11, 0], -302968512]
[27, [1, -3, -12, 4, 10], -2428538112]
[27, [1, -8, 13, -6, 0], -376451712]
[27, [6, 2, -7, 9, -10], -1720922112]
[27, [6, 7, -12, 4, -5], -273030912]
[27, [6, -3, -2, -11, 10], -1154829312]
[27, [6, -3, -12, 9, 0], -2142770112]
[27, [6, -13, 8, -1, 0], -376451712]
[27, [11, 2, 3, -6, -10], -1154829312]
[27, [11, 2, -12, -1, 0], -302968512]
[27, [-4, 2, -12, 9, 5], -2357776512]
[27, [-4, 12, 3, -1, -10], -2428538112]
[27, [-4, 12, -7, -6, 5], -273030912]
[27, [-4, -3, 8, 9, -10], 295783488]
[27, [-4, -3, -7, 14, 0], -90683712]
[27, [-9, 7, -2, -6, 10], -1720922112]
[27, [-9, 12, 3, -6, 0], -2142770112]
[27, [-9, 12, -2, 4, -5], -2357776512]
[27, [-9, -8, 3, 4, 10], 295783488]
[27, [-14, 7, 3, 4, 0], -90683712]
[28, [1, 7, -7, 9, -10], 1833689088]
[28, [1, -3, 13, -1, -10], -685504512]
[28, [1, -3, -7, 14, -5], 270885888]
[28, [1, -13, 3, -1, 10], -685504512]
[28, [1, -13, -2, 9, 5], -1803174912]
[28, [6, 7, -7, -11, 5], -85542912]
[28, [6, -8, 8, 4, -10], -346816512]
[28, [11, 2, 3, -11, -5], 1549836288]
[28, [11, 7, -7, -6, -5], -85542912]
[28, [11, -3, -2, -11, 5], 1549836288]
[28, [-4, 2, -12, 4, 10], 2543321088]
[28, [-4, 12, -2, 4, -10], 2543321088]
[28, [-4, -8, 8, -6, 10], -346816512]
[28, [-4, -8, -2, 14, 0], 272498688]
[28, [-9, 2, 13, -1, -5], -1803174912]
[28, [-9, 7, -7, -1, 10], 1833689088]
[28, [-14, 2, 8, 4, 0], 272498688]
[28, [-14, 7, 3, -1, 5], 270885888]
true

150922

Ruby

Ulam spiral

素数を塗りつぶす場合、p角数を塗りつぶす場合両方に対応。

# -*- coding: cp932 -*-

require 'prime'

# mはp角数か？
def p_number?(p, m)
  i = 0
  a = ((p - 2) * i - p + 4) * i / 2
  while a < m
    i += 1
    a = ((p - 2) * i - p + 4) * i / 2
  end
  a == m
end

# (x, y)に埋まる数字を算出
def cell(n, x, y, start = 1)
  y, x = y - n / 2, x - (n - 1) / 2
  l = 2 * [x.abs, y.abs].max
  if y >= x
    d = l * 3 + x + y
  else 
    d = l - x - y
  end
  (l - 1) * (l - 1) + d + start - 1
end

# condition = 0 のとき素数を、それ以外のときcondition角数を塗りつぶす
def show_spiral(condition, n, symbol, start = 1)
  puts "\nN : #{condition} : #{n}"
  n.times{|y|
    n.times{|x|
      m = cell(n, x, y, start)
      if condition == 0
        print m.prime? ? symbol[0] : symbol[1]
      else
        print p_number?(condition, m) ? symbol[0] : symbol[1]
      end
    }
    puts
  }
end

show_spiral(0, 40, "■□")
show_spiral(3, 40, "■□")
show_spiral(4, 40, "■□")
show_spiral(5, 40, "■□")
show_spiral(6, 40, "■□")

出力結果（実際より縦長になっていますが、コピペしてメモ帳に貼り付けるときれいな図になります。）

N : 0 : 40
□□□■□□□□□□□□□□□□□■□□□■□□□□□□□■□□□■□□□□□□
□□□□□□■□□□□□■□□□■□■□□□■□□□□□□□□□□□□□■□□□
□■□□□□□■□■□□□□□■□□□■□■□□□□□□□□□□□□□□□□□■
■□□□□□■□■□□□□□□□□□□□□□□□□□□□□□■□□□□□■□□□
□□□□□□□■□■□□□□□■□□□■□□□□□□□□□□□■□□□□□■□□
□□□□□□□□□□□□□□□□□□■□□□■□■□□□■□□□□□□□□□□□
□■□□□□□□□□□□□□□□□■□□□■□□□□□□□■□□□■□■□□□□
■□□□□□■□□□■□□□□□□□□□□□■□■□□□□□■□□□□□□□□□
□■□□□■□■□□□□□■□■□□□□□■□□□□□■□□□□□□□□□□□■
■□■□□□□□□□□□□□■□□□□□□□□□□□■□□□■□□□□□■□■□
□□□■□□□□□■□□□■□□□□□□□■□□□□□■□□□□□□□□□□□□
□□□□■□□□■□□□□□□□□□■□□□■□■□□□■□■□■□□□□□□□
□■□■□□□□□□□■□□□□□■□□□□□□□□□■□■□□□■□□□□□■
□□□□□□■□□□□□■□□□■□■□□□□□□□□□□□□□□□■□□□□□
□□□□□□□□□□□□□□□□□□□■□■□□□□□■□□□■□□□■□■□□
■□□□■□□□■□□□■□■□□□■□□□□□□□■□□□■□■□□□□□□□
□□□□□□□□□□□□□□□□□□□■□■□□□■□□□□□□□□□□□□□□
□□■□□□□□□□■□□□■□■□□□□□■□■□■□□□□□■□■□■□□□
□□□■□■□■□■□■□■□■□■□□□■□□□□□□□■□□□□□□□□□■
□□□□□□□□□□□□□□□□□□■□■□■□□□□□□□□□□□■□□□□□
□□□□□□□□□□□□□■□□□■□□■■□■□■□■□□□■□■□■□□□□
□□□□□□■□□□□□□□■□■□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□■□□□■□□□□□□□□□□□□□□□□□□□□
□□□□■□□□■□■□□□■□■□□□■□□□■□■□□□■□□□■□■□□□
□■□□□□□■□□□■□□□■□□□□□■□□□□□■□■□□□■□□□□□□
□□□□□□□□□□□□□□□□■□□□□□□□□□□□■□□□□□□□□□□□
□■□□□□□□□□□■□■□□□□□■□□□■□□□■□□□□□□□■□■□□
■□■□□□□□■□□□■□□□□□□□□□□□■□□□□□□□■□□□□□■□
□□□■□□□□□■□□□□□■□□□■□■□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□■□■□□□■□□□□□■□□□■□□□□□■□□□
□□□■□■□■□■□□□□□□□□□■□■□□□□□■□□□□□■□■□□□■
□□□□□□■□□□■□□□□□□□□□□□■□■□□□□□□□□□□□□□□□
□■□□□■□■□□□□□■□□□□□■□□□■□■□□□□□□□□□□□■□□
□□■□■□□□□□□□■□□□□□□□□□■□□□□□□□■□□□□□■□□□
□□□□□□□□□□□□□■□■□□□■□■□□□□□□□□□■□□□□□■□□
□□□□□□□□□□■□□□■□□□□□■□□□□□■□□□□□□□□□■□■□
□□□□□□□■□■□□□□□■□□□□□□□□□□□□□□□□□■□□□□□□
■□□□■□□□□□□□■□□□□□□□□□□□■□□□■□□□□□■□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□■□□□■
■□□□■□■□□□■□□□□□■□□□□□□□□□□□■□□□□□□□□□□□

N : 3 : 40
□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□■□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□■□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□■□□□□□□□□□□
□□□■□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□
□□□□□□□□□■□□□□□□□□□□□□□□■□□□□□□□□□■□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□■□■□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□■□□□□■□□□□□□□□□□□
□□□□□□□□□□□□□□■□□□□□■□□□□□□□□□□□□□□□□□□□
□□□□□□□□■□□□□□□□□□■■□□■□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□■
□□■□□□□□□□□□□■□□□■□□□□□□□□□■□□□□□■□□□□□□
□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□■□□□□□□□□□□□□□
□□□□□□□■□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□■□□□□□□□
□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□■□
□■□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□

N : 4 : 40
■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□

N : 5 : 40
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□
□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□■□□□□□□□□□□
□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□■□□
□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□■□□■□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□■□□□□□□□□□□□□□
□□□□□□□■□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□■□□□□■□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□
□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□
□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□
□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□

N : 6 : 40
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□■□□□□□□□□□□
□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□■□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□■□□□□□□□□□■■□□■□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■
□□□□□□□□□□□□□■□□□□□□□□□□□□□■□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□■□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□
□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□■□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□■□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□■□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□
□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□□

2015年9月21日月曜日

150921(3)

Ruby

階段状に現れるフィボナッチ数列

（n が増えると形が崩れてくるが）きれいに並んでいる。

def f(m, n)
  l = m - n
  ary = Array.new(l, 1)
  p ary
  n.times{|i|
    new_ary = ary.clone
    new_ary.shift
    (1..new_ary.size - 1).each{|i|
      new_ary[i] += new_ary[i - 1]
    }
    new_ary.push(new_ary[-1])
    ary = new_ary
    p [0] * (i + 1) + ary
  }
end

l = 4
n = 6
m = n + l
f(m, n)

出力結果
[1, 1, 1, 1]
[0, 1, 2, 3, 3]
[0, 0, 2, 5, 8, 8]
[0, 0, 0, 5, 13, 21, 21]
[0, 0, 0, 0, 13, 34, 55, 55]
[0, 0, 0, 0, 0, 34, 89, 144, 144]
[0, 0, 0, 0, 0, 0, 89, 233, 377, 377]

150921(2)

Ruby

縦読み、横読みの一般化

以下のサイトで掲題について質問してみた。
（http://ja.stackoverflow.com/questions/16776/%E5%A4%A7%E3%81%8D%E3%81%AA%E8%A1%8C%E5%88%97%E3%81%8B%E3%82%89%E6%8C%87%E5%AE%9A%E3%81%95%E3%82%8C%E3%81%9F%E6%96%87%E5%AD%97%E5%88%97%E3%82%92%E6%8E%A2%E3%81%97%E5%87%BA%E3%81%99%E3%81%AB%E3%81%AF）

orangecat さんの回答のコードを少しアレンジしてみた。

# -*- coding: cp932 -*-

W = 4
ary1 = [
'す', 'す', 'ま', '番兵',
'き', 'す', 'す', '番兵',
'す', 'き', 'だ', '番兵'
]
ary2 = [
'す', 'す', 'ま', '番兵',
'き', 'す', 'す', '番兵',
'よ', 'め', 'だ', '番兵'
]

def dpsearch(ary, str)
  # 番兵行を上に追加
  ary = Array.new(W, '番兵') + ary
  ss = Array.new(ary.size, [])
  last = str.size - 1
  W.upto(ary.size - 1){|pos|
    pre = (ss[pos - 1] + ss[pos - W]).uniq + [-1]
    # ary[pos] == str[s]のとき、str[0..s]まで一致している
    ss[pos] = pre.map{|s| s + 1}.select{|s| ary[pos] == str[s]}
    if ss[pos].include?(last)
      # posは番兵行を含んでいるが、元々はpos - W
      return pos - W
    end
  }
  nil
end

[ary1, ary2, ary2 * 1000 + ary1].each{|ary|
  pos = dpsearch(ary, 'すすき')
  if pos
    puts "Found (last char #{ary[pos]} at #{pos % W}, #{pos / W})"
  else
    puts "Not found"
  end
}

出力結果
Found (last char き at 1, 2)
Not found
Found (last char き at 1, 3002)

150921

Ruby

Half-Catalan number

オンライン整数列大辞典の
A000992（http://oeis.org/A000992/list）
と比較し、答え合わせしてみる。

def A000992(n)
  ary = [0, 1]
  i = 2
  while i <= n
    s = 0
    (1..i / 2).each{|j| s += ary[j] * ary[i - j]}
    ary[i] = s
    i += 1
  end
  ary[1..-1]
end
ary = A000992(32)

# OEIS A000992のデータ
ary0 =
[1,1,1,2,3,6,11,24,47,103,214,481,1030,2337,5131,
 11813,26329,60958,137821,321690,734428,1721998,
 3966556,9352353,21683445,51296030,119663812,
 284198136,666132304,1586230523,3734594241,
 8919845275]
# 一致の確認
p ary == ary0

2015年9月20日日曜日

150920

Ruby

分割が絡んだ係数について

Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables

に載っているTable 24.2のM3を求めるコードを書いてみた。

（ただし、順番が異なっている箇所があります。）

# 階乗
def f(ary)
  ary.map{|i| (1..i).inject(:*)}.inject(:*)
end

# nの分割を列挙
def partion_ary(n)
  ary = [[n]]
  # nの分割になっていないものは以下の操作を繰り返す
  f_ary = []
  (n - 1).downto(1){|i| f_ary.push([i])}
  while f_ary.size > 0
    b_ary = []
    f_ary.each{|c|
      i = [n - c.inject(:+), c[-1]].min
      (1..i).each{|j|
        d = c.clone
        d.push(j)
        # nの分割になったら、aryに加える
        d.inject(:+) == n ? ary.push(d) : b_ary.push(d)
      }
    }
    f_ary = b_ary
  end
  ary
end

def multinomials_and_partitions(n)
  m = (1..n).inject(:*)
  ary = []
  partion_ary(n).each{|a|
    ary << m / (f(a) * f(a.group_by(&:to_i).map{|i| i.last.size}))
  }
  ary
end

(1..15).each{|i| p multinomials_and_partitions(i)}

出力結果
[1]
[1, 1]
[1, 3, 1]
[1, 4, 3, 6, 1]
[1, 5, 10, 10, 15, 10, 1]
[1, 6, 15, 10, 15, 60, 15, 20, 45, 15, 1]
[1, 7, 21, 35, 21, 105, 105, 70, 35, 210, 105, 35, 105, 21, 1]
[1, 8, 28, 56, 35, 28, 168, 210, 280, 280, 56, 420, 840, 280, 105, 70, 560, 420, 56, 210, 28, 1]
[1, 9, 36, 84, 126, 36, 252, 378, 504, 1260, 315, 280, 84, 756, 1890, 1260, 1260, 2520, 126, 1260, 3780, 840, 945, 126, 1260, 1260, 84, 378, 36, 1]
[1, 10, 45, 120, 210, 126, 45, 360, 630, 840, 2520, 1260, 2100, 1575, 120, 1260, 3780, 2520, 3150, 12600, 1575, 6300, 2800, 210, 2520, 9450, 4200, 12600, 12600, 945, 252, 3150, 12600, 2100, 4725, 210, 2520, 3150, 120, 630, 45, 1]
[1, 11, 55, 165, 330, 462, 55, 495, 990, 1320, 4620, 2310, 4620, 6930, 1386, 5775, 165, 1980, 6930, 4620, 6930, 27720, 6930, 34650, 23100, 17325, 15400, 330, 4620, 20790, 9240, 34650, 69300, 5775, 17325, 69300, 15400, 462, 6930, 34650, 11550, 69300, 46200, 10395, 462, 6930, 34650, 4620, 17325, 330, 4620, 6930, 165, 990, 55, 1]
[1, 12, 66, 220, 495, 792, 462, 66, 660, 1485, 1980, 7920, 3960, 9240, 13860, 5544, 27720, 8316, 5775, 220, 2970, 11880, 7920, 13860, 55440, 13860, 83160, 55440, 83160, 8316, 138600, 51975, 69300, 15400, 495, 7920, 41580, 18480, 83160, 166320, 27720, 51975, 415800, 138600, 103950, 138600, 184800, 792, 13860, 83160, 27720, 207900, 277200, 17325, 207900, 415800, 61600, 10395, 924, 16632, 103950, 27720, 277200, 138600, 62370, 792, 13860, 83160, 9240, 51975, 495, 7920, 13860, 220, 1485, 66, 1]
[1, 13, 78, 286, 715, 1287, 1716, 78, 858, 2145, 2860, 12870, 6435, 17160, 25740, 10296, 60060, 36036, 6006, 45045, 36036, 286, 4290, 19305, 12870, 25740, 102960, 25740, 180180, 120120, 180180, 36036, 360360, 270270, 360360, 108108, 200200, 450450, 75075, 715, 12870, 77220, 34320, 180180, 360360, 60060, 135135, 1081080, 360360, 540540, 36036, 900900, 1801800, 675675, 450450, 600600, 200200, 1287, 25740, 180180, 60060, 540540, 720720, 90090, 675675, 2702700, 600600, 450450, 270270, 1801800, 1201200, 1716, 36036, 270270, 72072, 900900, 900900, 45045, 1351350, 1801800, 200200, 135135, 1716, 36036, 270270, 60060, 900900, 360360, 270270, 1287, 25740, 180180, 17160, 135135, 715, 12870, 25740, 286, 2145, 78, 1]
[1, 14, 91, 364, 1001, 2002, 3003, 1716, 91, 1092, 3003, 4004, 20020, 10010, 30030, 45045, 18018, 120120, 72072, 24024, 105105, 168168, 42042, 126126, 364, 6006, 30030, 20020, 45045, 180180, 45045, 360360, 240240, 360360, 72072, 840840, 630630, 840840, 504504, 42042, 560560, 2522520, 630630, 378378, 504504, 1051050, 525525, 1001, 20020, 135135, 60060, 360360, 720720, 120120, 315315, 2522520, 840840, 1261260, 168168, 2522520, 5045040, 3783780, 2522520, 756756, 6306300, 2802800, 1576575, 6306300, 525525, 1401400, 2002, 45045, 360360, 120120, 1261260, 1681680, 210210, 1891890, 7567560, 1681680, 2522520, 126126, 945945, 12612600, 12612600, 4729725, 2102100, 3153150, 8408400, 1401400, 3003, 72072, 630630, 168168, 2522520, 2522520, 252252, 4729725, 12612600, 2102100, 1576575, 3783780, 12612600, 5605600, 135135, 3432, 84084, 756756, 168168, 3153150, 2522520, 105105, 6306300, 6306300, 560560, 945945, 3003, 72072, 630630, 120120, 2522520, 840840, 945945, 2002, 45045, 360360, 30030, 315315, 1001, 20020, 45045, 364, 3003, 91, 1]
[1, 15, 105, 455, 1365, 3003, 5005, 6435, 105, 1365, 4095, 5460, 30030, 15015, 50050, 75075, 30030, 225225, 135135, 45045, 225225, 360360, 180180, 25740, 630630, 210210, 126126, 455, 8190, 45045, 30030, 75075, 300300, 75075, 675675, 450450, 675675, 135135, 1801800, 1351350, 1801800, 1081080, 180180, 1401400, 6306300, 1576575, 1891890, 2522520, 630630, 6306300, 4729725, 3783780, 1891890, 2627625, 1365, 30030, 225225, 100100, 675675, 1351350, 225225, 675675, 5405400, 1801800, 2702700, 360360, 6306300, 12612600, 9459450, 6306300, 3783780, 210210, 18918900, 8408400, 9459450, 37837800, 4729725, 5675670, 3783780, 21021000, 23648625, 15765750, 7882875, 1401400, 3003, 75075, 675675, 225225, 2702700, 3603600, 450450, 4729725, 18918900, 4204200, 6306300, 630630, 2837835, 37837800, 37837800, 28378350, 12612600, 3783780, 23648625, 94594500, 21021000, 23648625, 47297250, 2627625, 21021000, 21021000, 5005, 135135, 1351350, 360360, 6306300, 6306300, 630630, 14189175, 37837800, 6306300, 9459450, 378378, 14189175, 94594500, 63063000, 23648625, 7882875, 4729725, 47297250, 63063000, 7007000, 6435, 180180, 1891890, 420420, 9459450, 7567560, 630630, 23648625, 47297250, 6306300, 4729725, 28378350, 63063000, 21021000, 2027025, 6435, 180180, 1891890, 360360, 9459450, 6306300, 225225, 23648625, 18918900, 1401400, 4729725, 5005, 135135, 1351350, 225225, 6306300, 1801800, 2837835, 3003, 75075, 675675, 50050, 675675, 1365, 30030, 75075, 455, 4095, 105, 1]

ちなみに、面白い性質を持っている。
例えば、n = 4 のとき、
1×(x + 3) + 4×(x + 3)(x + 2) + 3×(x + 3)(x + 2) + 6×(x + 3)(x + 2)(x + 1) + 1×(x + 3)(x + 2)(x + 1)x
= (x + 3)^4
となっている。

2015年9月19日土曜日

150919

Ruby

二重根号が外れてきれいになる式

√(( 5＋2√3)( 4＋√13)) - √(( 4＋ √3)( 5＋√13)) = 1

√(( 7＋3√3)( 5＋√22)) - √(( 5＋ √3)( 7＋√22)) = 2

√((23＋8√8)( 7＋√17)) - √(( 7＋2√8)(23＋√17)) = 4

√((13＋9√2)( 5＋√ 7)) - √(( 5＋3√2)(13＋√ 7)) = 2

√((16＋9√3)(11＋√13)) - √((11＋6√3)(16＋√13)) = 1

を以下のコードで見つけ、

http://www.artofproblemsolving.com/community/c4h1142129_curious_equation

で証明した方法で全て（手計算で）確認した。

def square_number?(n)
  (Math.sqrt(n).to_i) ** 2 == n
end

M = 10
N = 30
(1..M).each{|f2|
  (f2 + 1..M).each{|d|
    (2..M).each{|x|
      if !square_number?(x)
        (x + 1..N).each{|y|
          if x.gcd(y) == 1 && !square_number?(y)
            (1..N).each{|e2|
              (e2 + 1..N).each{|c|
                m0 = Math.sqrt(c + d * Math.sqrt(x)) * Math.sqrt(e2 + Math.sqrt(y))
                m1 = Math.sqrt(e2 + f2 * Math.sqrt(x)) * Math.sqrt(c + Math.sqrt(y))
                z = m0 - m1
                p [[c, d, e2, f2, x, y], z.to_i] if (z - z.to_i).abs < 0.0000001 && z > 0
              }
            }
          end
        }
      end
    }
  }
}

出力結果
[[5, 2, 4, 1, 3, 13], 1]
[[7, 3, 5, 1, 3, 22], 2]
[[23, 8, 7, 2, 8, 17], 4]
[[13, 9, 5, 3, 2, 7], 2]
[[16, 9, 11, 6, 3, 13], 1]

2015年9月15日火曜日

150915

Ruby

ラマヌジャンが見つけた等式

http://www.nn.iij4u.or.jp/~hsat/misc/math/ramanujan.html
上記ページに載っている等式の両辺をそれぞれ計算してみた。

p (2 ** (1.0 / 3) - 1) ** (1.0 / 3)
p (1.0 / 9) ** (1.0 / 3) - (2.0 / 9) ** (1.0 / 3) + (4.0 / 9) ** (1.0 / 3)

p (4 ** (1.0 / 5) + 1) ** (1.0 / 2)
p (16 ** (1.0 / 5) + 8 ** (1.0 / 5) + 2 ** (1.0 / 5) - 1) / (5 ** (1.0 / 2))

p ((3 + 2 * 5 ** (1.0 / 4)) / (3 - 2 * 5 ** (1.0 / 4))) ** (1.0 / 4)
p (5 ** (1.0 / 4) + 1) / (5 ** (1.0 / 4) - 1)

p (5 ** (1.0 / 3) - 4 ** (1.0 / 3)) ** (1.0 / 2)
p (2 ** (1.0 / 3) + 20 ** (1.0 / 3) - 25 ** (1.0 / 3)) / 3

出力結果

0.6381858208606442
0.6381858208606441
1.5229930764034663
1.522993076403466
5.037559141801549
5.037559141801561
0.3501069760923045
0.3501069760923044

2015年9月14日月曜日

150914(2)

Ruby

xx + 27yy 型の素数

tsujimotter さんの本日のブログ

(http://tsujimotter.hatenablog.com/entry/xx-27yy)

について、自分で実験してみることにした。

require 'prime'

# 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

# 1 - q^(6n)
def f1(n)
  f1 = Array.new(6 * n + 1, 0)
  f1[0] = 1
  f1[-1] = -1
  f1
end

# 1 - q^(18n)
def f2(n)
  f2 = Array.new(18 * n + 1, 0)
  f2[0] = 1
  f2[-1] = -1
  f2
end

ary = [0, 1]
n = 1
while 6 * n <= 300
  ary = mul(ary, f1(n), 300)
  n += 1
end
n = 1
while 18 * n <= 300
  ary = mul(ary, f2(n), 300)
  n += 1
end
p ary

Prime.each(300){|i| p [i, ary[i]]}

出力結果

[0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

[2, 0]

[3, 0]

[5, 0]

[7, -1]

[11, 0]

[13, -1]

[17, 0]

[19, -1]

[23, 0]

[29, 0]

[31, 2]

[37, -1]

[41, 0]

[43, 2]

[47, 0]

[53, 0]

[59, 0]

[61, -1]

[67, -1]

[71, 0]

[73, -1]

[79, -1]

[83, 0]

[89, 0]

[97, -1]

[101, 0]

[103, -1]

[107, 0]

[109, 2]

[113, 0]

[127, 2]

[131, 0]

[137, 0]

[139, -1]

[149, 0]

[151, -1]

[157, 2]

[163, -1]

[167, 0]

[173, 0]

[179, 0]

[181, -1]

[191, 0]

[193, -1]

[197, 0]

[199, -1]

[211, -1]

[223, 2]

[227, 0]

[229, 2]

[233, 0]

[239, 0]

[241, -1]

[251, 0]

[257, 0]

[263, 0]

[269, 0]

[271, -1]

[277, 2]

[281, 0]

[283, 2]

[293, 0]

赤字の素数がxx + 27yy 型の素数となる。

150914

C

Number of knight's tours on a m×n chessboard(3)

さらに、3×21 のときも求まる。

m = 3 のとき

kto3s =

[0, 0, 0, 16, 0, 0, 104, 792, 1120, 6096, 21344, 114496,

257728, 1292544, 3677568, 17273760, 46801984,

211731376, 611507360, 2645699504, 7725948608]

m = 6 のとき

kto6s =

[0, 0, 0, 1488, 37568, 6637920, 779938932]

m = 7 のとき

kto7s =

[0, 0, 104, 12756, 1245736, 779938932, 165575218320]

2015年9月13日日曜日

150913

C

Number of knight's tours on a m×n chessboard(2)

さらに、3×19, 3×20 のときも求まる。

7×7 のときも含めて現在分かっている値をまとめておく。

m = 3 のとき

kto3s =

[0, 0, 0, 16, 0, 0, 104, 792, 1120, 6096, 21344,

114496, 257728, 1292544, 3677568, 17273760,

46801984, 211731376, 611507360, 2645699504]

m = 6 のとき

kto6s =

[0, 0, 0, 1488, 37568, 6637920, 779938932]

m = 7 のとき

kto7s =

[0, 0, 104, 12756, 1245736, 779938932, 165575218320]

2015年9月12日土曜日

150912

C

Number of knight's tours on a m×n chessboard(1)

yoh2 さんのコード（VS2013用）を用いる。

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <assert.h>


struct board {
 // 盤面の大きさ
 int num_rows, num_cols;

 //counts[行インデックス][列インデックス]:
 //そのマスから/マスへ移動可能なマスがいくつあるか。
 //解の探索で、利き筋のマスが訪問済になると減少する。
 //また、訪問済のマスは -1 とする。
 int **counts;

 // カウントが1のマスの数。これが3以上になると
 // そのすべてを訪問することはできないため解なしとなる。
 int num_ones;
};


// マスの位置を表す構造体。
struct board_index {
 int row, col;
};


// ナイトの利き筋。
const struct board_index g_knight_movement[] = {
 { 1, 2 }, { -1, 2 }, { 1, -2 }, { -1, -2 },
 { 2, 1 }, { -2, 1 }, { 2, -1 }, { -2, -1 },
};
// g_knight_movement の要素数。
#define NUM_KNIGHT_MOVEMENTS  ((sizeof(g_knight_movement) / sizeof(g_knight_movement[0])))


// 指定されたインデックスが盤面の範囲に入っているか。
static bool in_board(const struct board * brd, struct board_index index)
{
 return ((0 <= index.row) && (index.row < brd->num_rows))
  && ((0 <= index.col) && (index.col < brd->num_cols));
}


// 指定したマスの利き筋の数の初期値を求める。
int initial_count_at(const struct board * brd, int row, int col)
{
 int count = 0;
 for (size_t i = 0; i < NUM_KNIGHT_MOVEMENTS; i++) {
  struct board_index next = { // (row, col) から一回動いた位置
   .row = g_knight_movement[i].row + row,
   .col = g_knight_movement[i].col + col,
  };
  count += in_board(brd, next);
 }
 return count;
}


// init_board_result() の戻り値。
enum init_board_result {
 IB_SUCCESS,         // 初期化成功
 IB_ERROR,           // エラー発生
 IB_NO_SOLUTION,     // 初期化時に解がないことが確定して初期化を中断した
};


// 盤面の情報を初期化する。
enum init_board_result init_board(struct board * brd, int num_rows, int num_cols)
{
 enum init_board_result err_retval;
 brd->num_rows = num_rows;
 brd->num_cols = num_cols;
 brd->num_ones = 0;
 brd->counts = (int **)malloc(sizeof(brd->counts[0]) * num_rows);
 if (brd->counts == NULL) {
  err_retval = IB_ERROR;
  goto err_alloc_baseptr;
 }
 int row;
 for (row = 0; row < num_rows; row++) {
  brd->counts[row] = (int *)malloc(sizeof(brd->counts[0][0]) * num_cols);
  if (brd->counts[row] == NULL) {
   err_retval = IB_ERROR;
   goto err_alloc_row;
  }
  for (int col = 0; col < num_cols; col++) {
   brd->counts[row][col] = initial_count_at(brd, row, col);
   switch (brd->counts[row][col]) {
   case 1:
    brd->num_ones++;
    break;

   case 0:
    err_retval = IB_NO_SOLUTION;
    goto no_solution;

   default:
    break;
   }
  }
 }
 return IB_SUCCESS;

 // error or no solution
no_solution:
 free(brd->counts[row]);
err_alloc_row:
 while (row-- > 0) {
  free(brd->counts[row]);
 }
 free(brd->counts);
err_alloc_baseptr:
 return err_retval;
}


// 盤面を破棄する。
void destray_board(struct board * brd)
{
 for (int row = 0; row < brd->num_rows; row++) {
  free(brd->counts[row]);
 }
 free(brd->counts);
}


// 後述の update_next_counts_and_enum_available_indices() とペアで使う。
// 上記関数で列挙されたマスのカウントを元に戻す。
void restore_next_counts(
struct board * brd,
 const struct board_index nexts[],
 size_t num_nexts)
{
 for (size_t i = 0; i < num_nexts; i++) {
  // カウント戻し + 1のマス数の調整
  switch (brd->counts[nexts[i].row][nexts[i].col]++) {
  case 1:
   brd->num_ones--;
   break;

  case 0:
   brd->num_ones++;
   break;
  }
 }
}


// 現在位置から到達できるマスのカウントを減らしつつ、
// 指定された位置から次に動ける候補一覧を列挙する。
// nexts[] に候補が格納され、候補数が返される。解がないと判断されたら
// カウントの変更はキャンセルされ、0が返されら。
size_t update_next_counts_and_enum_available_indices(
struct board * brd,
struct board_index index,
 bool is_last_step,
struct board_index nexts[])
{
 size_t num_nexts = 0;
 for (size_t i = 0; i < NUM_KNIGHT_MOVEMENTS; i++) {
  struct board_index next = {
   .row = g_knight_movement[i].row + index.row,
   .col = g_knight_movement[i].col + index.col,
  };
  if (!in_board(brd, next)) {
   // 盤面の範囲外は候補たり得ないので何もしない
   continue;
  }
  int count = brd->counts[next.row][next.col];
  if (count < 0) {
   // 訪問済のマスは触らない。
   continue;
  }
  // 以降、すべて正の数となる。(0となるパターンは存在しない)
  assert(count != 0);
  // カウントを減らす。
  brd->counts[next.row][next.col] = count - 1;
  if (count == 1) {
   // 1のマス (今回0になった) が見付かった時点で、それが唯一の候補。
   // ただし、それが最後の空きマスでなければ、さらに
   // 次の一歩を移動できないので解なしとなる。
   if (is_last_step) {
    num_nexts = 1;
    nexts[0] = next;
    brd->num_ones--;
    // 他のマスはすべて訪問済なので、これ以上見る必要はない。
    break;
   }
   else {
    // 解なし。変更を元に戻して0を返す。
    brd->counts[next.row][next.col] = count;
    restore_next_counts(brd, nexts, num_nexts);
    return 0;
   }
  }
  else {
   nexts[num_nexts++] = next;
   if (count == 2) {
    // 減算の結果、1になるので、記録してある1の数を増やす。
    brd->num_ones++;
   }
  }
 }
 assert(num_nexts > 0);
 return num_nexts;
}


long long solve_knight_tour_impl(
struct board * brd,
 int depth,
struct board_index index)
{
 if (depth == brd->num_rows * brd->num_cols - 1) {
  // 解がひとつ見付かった
  return 1;
 }

 struct board_index nexts[NUM_KNIGHT_MOVEMENTS]; // 移動先候補は最大でナイトの利き筋の数。
 size_t num_nexts = update_next_counts_and_enum_available_indices(
  brd, index,
  depth == brd->num_rows * brd->num_cols - 2,
  nexts);
 if (num_nexts == 0) {
  // この探索では解がない。
  return 0;
 }
 // 現在位置を無効化する。
 // 後で戻せるように値の退避もしておく。
 int saved_count = brd->counts[index.row][index.col];
 if (saved_count == 1) {
  brd->num_ones--;
 }
 brd->counts[index.row][index.col] = -1;

 long long found = 0;
 // ここで、現在の1のマスの数を調べる。
 if (brd->num_ones >= 3) {
  // 3以上なら解なしとなるので何もしない。
 }
 else {
  if (brd->num_ones >= 2) {
   // カウントが1の移動先候補それぞれについて解を探索。
   for (size_t i = 0; i < num_nexts; i++) {
    if (brd->counts[nexts[i].row][nexts[i].col] == 1) {
     found += solve_knight_tour_impl(brd, depth + 1, nexts[i]);
    }
   }
  }
  else {
   // 移動先候補それぞれについて解を探索。
   for (size_t i = 0; i < num_nexts; i++) {
    found += solve_knight_tour_impl(brd, depth + 1, nexts[i]);
   }
  }
 }
 // [A] の変更元に戻す。
 brd->counts[index.row][index.col] = saved_count;
 if (saved_count == 1) {
  brd->num_ones++;
 }
 restore_next_counts(brd, nexts, num_nexts);

 return found;
}


// num_rows × num_cols の盤面のナイトツアーを満たす解の数を返す。
// エラー発生時は -1 。
// なお、回転・線対称・逆周りはすべて別個の解として数えている。
long long solve_knight_tour(size_t num_rows, size_t num_cols)
{
 struct board brd;
 switch (init_board(&brd, num_rows, num_cols)) {
 case IB_SUCCESS:
  break;

 case IB_ERROR:
  return -1;

 case IB_NO_SOLUTION:
  return 0;
 }
 long long found = 0;

 // 各座標を視点とした解の個数を合計する。
 // なお、縦横2等分すれば、それぞれ線対称の位置になるので、
 // その 1/4 の部分だけ計算し、4倍することで全体の結果が求められる。
 // ただし、奇数の場合は、中央の行・列の追加計算も必要。
 for (int row = 0; row < num_rows / 2; row++) {
  for (int col = 0; col < num_cols / 2; col++) {
   found += solve_knight_tour_impl(&brd, 0, (struct board_index){ .row = row, .col = col });
  }
 }
 found *= 4;
 // 奇数の場合の追加パターン
 if (num_rows & 1) {
  long long extra_found = 0;
  // 行数が奇数 → 真ん中の行を追加計算。これも半分だけ計算して2倍する。
  for (int col = 0; col < num_cols / 2; col++) {
   extra_found += solve_knight_tour_impl(&brd, 0, (struct board_index){ .row = num_rows / 2, .col = col });
  }
  found += extra_found * 2;
 }
 if (num_cols & 1) {
  long long extra_found = 0;
  // 列数が奇数 → 真ん中の列を追加計算。これも半分だけ計算して2倍する。
  for (int row = 0; row < num_rows / 2; row++) {
   extra_found += solve_knight_tour_impl(&brd, 0, (struct board_index){ .row = row, .col = num_cols / 2 });
  }
  found += extra_found * 2;
 }
 if ((num_rows & 1) && (num_cols & 1)) {
  // どちらも奇数 → 中央のマスを追加計算。
  found += solve_knight_tour_impl(&brd, 0, (struct board_index){ .row = num_rows / 2, .col = num_cols / 2 });
 }
 destray_board(&brd);
 return found;
}


// 引数に盤面のサイズを指定する。
// 不正な引数が渡された場合のエラー処理はしていない(手抜き)。
int main(int argc, char *argv[])
{
 if (argc < 3) {
  printf("Usage: %s #-rows #-cols\n", argv[0]);
  return 1;
 }
 int num_rows = strtol(argv[1], NULL, 0);
 int num_cols = strtol(argv[2], NULL, 0);

 // 探索本体
 long long num_solutions = solve_knight_tour(num_rows, num_cols);
 if (num_solutions == -1) {
  perror("");
  return 1;
 }
 printf("# of solutions = %lld\n", num_solutions);
 return 0;
}

m = 3 のとき

C:\Users\Seiichi>20150912.exe 3 1

# of solutions = 0

C:\Users\Seiichi>20150912.exe 3 2

# of solutions = 0

C:\Users\Seiichi>20150912.exe 3 3

# of solutions = 0

C:\Users\Seiichi>20150912.exe 3 4

# of solutions = 16

C:\Users\Seiichi>20150912.exe 3 5

# of solutions = 0

C:\Users\Seiichi>20150912.exe 3 6

# of solutions = 0

C:\Users\Seiichi>20150912.exe 3 7

# of solutions = 104

C:\Users\Seiichi>20150912.exe 3 8

# of solutions = 792

C:\Users\Seiichi>20150912.exe 3 9

# of solutions = 1120

C:\Users\Seiichi>20150912.exe 3 10

# of solutions = 6096

C:\Users\Seiichi>20150912.exe 3 11

# of solutions = 21344

C:\Users\Seiichi>20150912.exe 3 12
# of solutions = 114496

C:\Users\Seiichi>20150912.exe 3 13

# of solutions = 257728

C:\Users\Seiichi>20150912.exe 3 14

# of solutions = 1292544

C:\Users\Seiichi>20150912.exe 3 15

# of solutions = 3677568

C:\Users\Seiichi>20150912.exe 3 16

# of solutions = 17273760

C:\Users\Seiichi>20150912.exe 3 17

# of solutions = 46801984

C:\Users\Seiichi>20150912.exe 3 18

# of solutions = 211731376

m = 6 のとき

C:\Users\Seiichi>20150912.exe 6 1

# of solutions = 0

C:\Users\Seiichi>20150912.exe 6 2

# of solutions = 0

C:\Users\Seiichi>20150912.exe 6 3

# of solutions = 0

C:\Users\Seiichi>20150912.exe 6 4

# of solutions = 1488

C:\Users\Seiichi>20150912.exe 6 5

# of solutions = 37568

C:\Users\Seiichi>20150912.exe 6 6

# of solutions = 6637920

C:\Users\Seiichi>20150912.exe 6 7

# of solutions = 779938932

m = 7 のとき

C:\Users\Seiichi>20150912.exe 7 1

# of solutions = 0

C:\Users\Seiichi>20150912.exe 7 2

# of solutions = 0

C:\Users\Seiichi>20150912.exe 7 3

# of solutions = 104

C:\Users\Seiichi>20150912.exe 7 4

# of solutions = 12756

C:\Users\Seiichi>20150912.exe 7 5

# of solutions = 1245736

C:\Users\Seiichi>20150912.exe 7 6

# of solutions = 779938932

2015年9月10日木曜日

150910(2)

コード用(1)

これから、以下を使用する。

コード用　コード用　コード用　コード用　コード用　コード用　コード用　コード用　コード用

150910

Ruby

Number of knight's tours on a 3×k chessboard(2)

Cに比べて遅いため、h = 14 までしか計算していない。

@knight_move =
[[2, -1], [2, 1], [-2, -1], [-2, 1],
[1, -2], [1, 2], [-1, -2], [-1, 2]]

# startが[x, y]
def search(x, y, w, h, used, depth)
  return 0 if x < 0 || w <= x || y < 0 || h <= y || used[x + y * w]
  return 1 if depth == w * h
  cnt = 0
  used[x + y * w] = true
  @knight_move.each{|m|
    cnt += search(x + m[0], y + m[1], w, h, used, depth + 1)
  }
  used[x + y * w] = false
  return cnt
end

def directed_open_tours3(h)
  return 0 if h < 3
  total = 0
  (0..h / 2 - 1).each{|y|
    used = Array.new(3 * h, false)
    total += search(0, y, 3, h, used, 1) * 4
    p [y, 0, total]
    used = Array.new(3 * h, false)
    total += search(1, y, 3, h, used, 1) * 2
    p [y, 1, total]
  }
  if h % 2 == 1
    y = h / 2
    used = Array.new(3 * h, false)
    total += search(0, y, 3, h, used, 1) * 2
    used = Array.new(3 * h, false)
    total += search(1, y, 3, h, used, 1) * 1
  end
  total
end

(1..14).each{|h|
  p directed_open_tours3(h)
}

出力結果

[0, 0, 0]

[0, 1, 0]

[0, 0, 8]

[0, 1, 16]

[1, 0, 16]

[1, 1, 16]

[0, 0, 0]

[0, 1, 0]

[1, 0, 0]

[1, 1, 0]

[0, 0, 0]

[0, 1, 0]

[1, 0, 0]

[1, 1, 0]

[2, 0, 0]

[2, 1, 0]

[0, 0, 32]

[0, 1, 32]

[1, 0, 32]

[1, 1, 88]

[2, 0, 104]

[2, 1, 104]

104

[0, 0, 328]

[0, 1, 464]

[1, 0, 544]

[1, 1, 640]

[2, 0, 656]

[2, 1, 656]

[3, 0, 752]

[3, 1, 792]

792

[0, 0, 688]

[0, 1, 688]

[1, 0, 688]

[1, 1, 880]

[2, 0, 912]

[2, 1, 912]

[3, 0, 912]

[3, 1, 992]

1120

[0, 0, 1792]

[0, 1, 2816]

[1, 0, 3440]

[1, 1, 4272]

[2, 0, 4464]

[2, 1, 4528]

[3, 0, 4784]

[3, 1, 4944]

[4, 0, 5584]

[4, 1, 6096]

6096

[0, 0, 10416]

[0, 1, 10416]

[1, 0, 10416]

[1, 1, 14416]

[2, 0, 15584]

[2, 1, 15584]

[3, 0, 15584]

[3, 1, 16624]

[4, 0, 20224]

[4, 1, 20224]

21344

[0, 0, 38632]

[0, 1, 51320]

[1, 0, 58840]

[1, 1, 70072]

[2, 0, 73880]

[2, 1, 74584]

[3, 0, 80360]

[3, 1, 84544]

[4, 0, 97168]

[4, 1, 102592]

[5, 0, 110144]

[5, 1, 114496]

114496

[0, 0, 121968]

[0, 1, 121968]

[1, 0, 121968]

[1, 1, 158800]

[2, 0, 167664]

[2, 1, 167664]

[3, 0, 167664]

[3, 1, 180416]

[4, 0, 223168]

[4, 1, 223168]

[5, 0, 223168]

[5, 1, 243456]

257728

[0, 0, 412512]

[0, 1, 528176]

[1, 0, 597840]

[1, 1, 727952]

[2, 0, 762928]

[2, 1, 769072]

[3, 0, 813936]

[3, 1, 853040]

[4, 0, 988656]

[4, 1, 1040560]

[5, 0, 1116400]

[5, 1, 1177232]

[6, 0, 1257680]

[6, 1, 1292544]

1292544

2015年9月8日火曜日

150908

C

Number of knight's tours on a 3×k chessboard(1)

一日半かけて計算した。

#include <stdio.h>

int search(int x, int y, int w, int h, long long used, int depth){
  int cnt = 0;
  if (x < 0 || w <= x || y < 0 || h <= y || (used & (1LL << (x + y * w))) > 0) return 0;
  if (depth == w * h) return 1;
  used += 1LL << (x + y * w);
  cnt += search(x + 2, y - 1, w, h, used, depth + 1);
  cnt += search(x + 2, y + 1, w, h, used, depth + 1);
  cnt += search(x - 2, y - 1, w, h, used, depth + 1);
  cnt += search(x - 2, y + 1, w, h, used, depth + 1);
  cnt += search(x + 1, y - 2, w, h, used, depth + 1);
  cnt += search(x + 1, y + 2, w, h, used, depth + 1);
  cnt += search(x - 1, y - 2, w, h, used, depth + 1);
  cnt += search(x - 1, y + 2, w, h, used, depth + 1);
  used -= 1LL << (x + y * w);
  return cnt;
}

int directed_open_tours3(int h){
  if (h < 3) return 0;
  int y;
  long long used;
  int total = 0;
  for (y = 0; y < h / 2; y++){
    used = 0LL;
    total += search(0, y, 3, h, used, 1) * 4;
    used = 0LL;
    total += search(1, y, 3, h, used, 1) * 2;
  }
  if (h % 2 == 1){
    y = h / 2;
    used = 0LL;
    total += search(0, y, 3, h, used, 1) * 2;
    used = 0LL;
    total += search(1, y, 3, h, used, 1) * 1;
  }
  return total;
}

int main(void){
  int h;
  for (h = 1; h < 17; h++){
    printf("%d\n", directed_open_tours3(h));
  }
  return 0;
}

出力結果

104

792

1120

6096

21344

114496

257728

1292544

3677568

17273760

登録: 投稿 (Atom)