2015年11月29日日曜日

151129(4)

Ruby


break と exit の違い

どういうわけか「たのしいRuby」の索引にはexitがない!

# -*- coding: cp932 -*-

n = 4
m0 = 9
puts "breakのとき"
i, j = 0, 0
while i < n
  j = 0 if j == n || j % n == m0 % n
  while j < n
    m = i * n + j
    p m
    break if m == m0
    j += 1
  end
  i += 1
end

puts "exitのとき"
i, j = 0, 0
while i < n
  j = 0 if j == n || j % n == m0 % n
  while j < n
    m = i * n + j
    p m
    exit if m == m0
    j += 1
  end
  i += 1
end

出力結果
breakのとき
0
1
2
3
4
5
6
7
8
9
12
13
14
15
exitのとき
0
1
2
3
4
5
6
7
8
9

151129(3)

Ruby


puts と ハッシュ

「puts と 配列」と「puts と ハッシュ」とでは違いがある。

irb(main):001:0> puts [2, 3, 5][0]
2
=> nil
irb(main):002:0> puts {2 => 0, 3 => 1, 5 => 2}[2]
SyntaxError: (irb):2: syntax error, unexpected =>, expecting '}'
puts {2 => 0, 3 => 1, 5 => 2}[2]
          ^
        from C:/Ruby22-x64/bin/irb:11:in `<main>'
irb(main):003:0> puts({2 => 0, 3 => 1, 5 => 2}[2])
0
=> nil

151129(2)

Ruby


正方形の形をした領域内のSelf-avoiding walk(3)

やっともう一つ値が求まった。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, n, used)
  return 0 if n + 1 <= x.abs || n + 1 <= y.abs || used[x + n + (y + n) * (2 * n + 1)]
  return 1 if x.abs == n || y.abs == n
  cnt = 0
  used[x + n + (y + n) * (2 * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], n, used)
  }
  used[x + n + (y + n) * (2 * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
(0..4).each{|n|
  used = Array.new((2 * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, n, used)
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
1
4
92
115516
9731587524

151129

Ruby


二等辺三角形の形をした領域内のSelf-avoiding walk(4)

左右上下に進む場合について
たどり着くy 座標で分類してみた。
(実行時間は3時間ほどかかる。)

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, m, n, used)
  ary0 = Array.new(2 * n + 1, 0)
  return ary0 if x < 0 || m * n + 1 <= x || (m * y).abs > x || used[x + (y + n) * (m * n + 1)]
  if x == m * n
    # y = -n が0
    ary0[y + n] = 1
    return ary0
  end
  c_ary = Array.new(2 * n + 1, 0)
  used[x + (y + n) * (m * n + 1)] = true
  move.each{|mo|
    ary1 = search(move, x + mo[0], y + mo[1], m, n, used)
    c_ary = (0..2 * n).map{|i| c_ary[i] + ary1[i]}
  }
  used[x + (y + n) * (m * n + 1)] = false
  return c_ary
end

move = move3
puts "move = #{move3}のとき"
m = 0
puts "m = 0 のとき"
(0..7).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  a_ary = search(move, 0, 0, m, n, used)
  p [a_ary.inject(:+), a_ary]
}
m = 1
puts "m = 1 のとき"
(0..7).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  a_ary = search(move, 0, 0, m, n, used)
  p [a_ary.inject(:+), a_ary]
}
m = 2
puts "m = 2 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  a_ary = search(move, 0, 0, m, n, used)
  p [a_ary.inject(:+), a_ary]
}
m = 3
puts "m = 3 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  a_ary = search(move, 0, 0, m, n, used)
  p [a_ary.inject(:+), a_ary]
}
(4..5).each{|m|
  puts "m = #{m}のとき"
  (0..3).each{|n|
    used = Array.new((m * n + 1) * (2 * n + 1), false)
    a_ary = search(move, 0, 0, m, n, used)
    p [a_ary.inject(:+), a_ary]
  }
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
m = 0 のとき
[1, [1]]
[1, [0, 1, 0]]
[1, [0, 0, 1, 0, 0]]
[1, [0, 0, 0, 1, 0, 0, 0]]
[1, [0, 0, 0, 0, 1, 0, 0, 0, 0]]
[1, [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]]
[1, [0, 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]]
m = 1 のとき
[1, [1]]
[1, [0, 1, 0]]
[3, [0, 1, 1, 1, 0]]
[15, [0, 3, 3, 3, 3, 3, 0]]
[153, [0, 25, 25, 19, 15, 19, 25, 25, 0]]
[4169, [0, 589, 589, 475, 311, 241, 311, 475, 589, 589, 0]]
[346211, [0, 42903, 42903, 35397, 25433, 18503, 15933, 18503, 25433, 35397, 42903, 42903, 0]]
[86336255, [0, 9578225, 9578225, 7879989, 5801937, 4565323, 3908305, 3712247, 3908305, 4565323, 5801937, 7879989, 9578225, 9578225, 0]]
m = 2 のとき
[1, [1]]
[1, [0, 1, 0]]
[9, [0, 3, 3, 3, 0]]
[345, [0, 83, 65, 49, 65, 83, 0]]
[95265, [0, 18039, 14933, 10329, 8663, 10329, 14933, 18039, 0]]
m = 3 のとき
[1, [1]]
[1, [0, 1, 0]]
[29, [0, 10, 9, 10, 0]]
[10199, [0, 2473, 1881, 1491, 1881, 2473, 0]]
[83918867, [0, 15568099, 12387799, 9603105, 8800861, 9603105, 12387799, 15568099, 0]]
m = 4のとき
[1, [1]]
[1, [0, 1, 0]]
[95, [0, 33, 29, 33, 0]]
[317897, [0, 75725, 58651, 49145, 58651, 75725, 0]]
m = 5のとき
[1, [1]]
[1, [0, 1, 0]]
[313, [0, 109, 95, 109, 0]]
[9947973, [0, 2346776, 1836159, 1582103, 1836159, 2346776, 0]]

2015年11月28日土曜日

151128(4)

Ruby


正方形の形をした領域内のSelf-avoiding walk(2)

経路の長さで分けてみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, n, used, depth)
  ary0 = Array.new((2 * n - 1) * (2 * n - 1) + 1, 0)
  return ary0 if n + 1 <= x.abs || n + 1 <= y.abs || used[x + n + (y + n) * (2 * n + 1)]
  if x.abs == n || y.abs == n
    ary0[depth] = 1
    return ary0
  end
  c_ary = Array.new((2 * n - 1) * (2 * n - 1) + 1, 0)
  used[x + n + (y + n) * (2 * n + 1)] = true
  move.each{|mo|
    ary1 = search(move, x + mo[0], y + mo[1], n, used, depth + 1)
    c_ary = (0..(2 * n - 1) * (2 * n - 1)).map{|i| c_ary[i] + ary1[i]}
  }
  used[x + n + (y + n) * (2 * n + 1)] = false
  return c_ary
end

N = 3
move = move3
puts "move = #{move3}のとき"
(0..N).each{|n|
  used = Array.new((2 * n + 1) * (2 * n + 1), false)
  a_ary = search(move, 0, 0, n, used, 0)
  p [a_ary.inject(:+), a_ary]
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
[1, [1, 0]]
[4, [0, 4]]
[92, [0, 0, 4, 16, 8, 16, 8, 16, 8, 16]]
[115516, [0, 0, 0, 4, 24, 72, 80, 176, 200, 472, 568, 1368, 1584, 3664, 4000, 8696, 8392, 16528, 12704, 21504, 10616, 14352, 4528, 4576, 888, 520]]

151128(3)

Ruby


(辺が斜めの)正方形の形をした領域内のSelf-avoiding walk(2)

経路の長さで分けてみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, n, used, depth)
  ary0 = Array.new((n - 1) * (n - 1) + n * n + 1, 0)
  return ary0 if n + 1 <= (x + y).abs || n + 1 <= (y - x).abs || used[x + n + (y + n) * (2 * n + 1)]
  if (x + y).abs == n || (y - x).abs == n
    ary0[depth] = 1
    return ary0
  end
  c_ary = Array.new((n - 1) * (n - 1) + n * n + 1, 0)
  used[x + n + (y + n) * (2 * n + 1)] = true
  move.each{|mo|
    ary1 = search(move, x + mo[0], y + mo[1], n, used, depth + 1)
    c_ary = (0..(n - 1) * (n - 1) + n * n).map{|i| c_ary[i] + ary1[i]}
  }
  used[x + n + (y + n) * (2 * n + 1)] = false
  return c_ary
end

N = 5
move = move3
puts "move = #{move3}のとき"
(0..N).each{|n|
  used = Array.new((2 * n + 1) * (2 * n + 1), false)
  a_ary = search(move, 0, 0, n, used, 0)
  p [a_ary.inject(:+), a_ary]
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
[1, [1, 0]]
[4, [0, 4]]
[12, [0, 0, 12, 0, 0, 0]]
[148, [0, 0, 0, 28, 0, 40, 0, 40, 0, 40, 0, 0, 0, 0]]
[15132, [0, 0, 0, 0, 60, 0, 184, 0, 376, 0, 976, 0, 2072, 0, 3904, 0, 4680, 0, 2880, 0, 0, 0, 0, 0, 0, 0]]
[12269676, [0, 0, 0, 0, 0, 124, 0, 576, 0, 1776, 0, 5592, 0, 19368, 0, 62080, 0, 185208, 0, 496744, 0, 1130232, 0, 2072216, 0, 2868480, 0, 2792536, 0, 1791656, 0, 704600, 0, 138488, 0, 0, 0, 0, 0, 0, 0, 0]]

151128(2)

Ruby


(辺が斜めの)正方形の形をした領域内のSelf-avoiding walk(1)

点の集合Dを
格子点で-n ≦ x + y ≦ n, -n ≦ y - x ≦ nを満たしているもの
とする。
今、Dの点(0, 0)から出発し、
自身が通った点を通ることなく、左右上下に進みながら
Dの境界上の格子点のいずれか一つへたどることとする。
このたどり方が何通りあるか出力してみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, n, used)
  return 0 if n + 1 <= (x + y).abs || n + 1 <= (y - x).abs || used[x + n + (y + n) * (2 * n + 1)]
  return 1 if (x + y).abs == n || (y - x).abs == n
  cnt = 0
  used[x + n + (y + n) * (2 * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], n, used)
  }
  used[x + n + (y + n) * (2 * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
(0..5).each{|n|
  used = Array.new((2 * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, n, used)
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
1
4
12
148
15132
12269676

151128

Ruby


正方形の形をした領域内のSelf-avoiding walk(1)

点の集合Dを
格子点で-n ≦ x ≦ n, -n ≦ y ≦ nを満たしているもの
とする。
今、Dの点(0, 0)から出発し、
自身が通った点を通ることなく、左右上下に進みながら
Dの境界上の格子点のいずれか一つへたどることとする。
このたどり方が何通りあるか出力してみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, n, used)
  return 0 if n + 1 <= x.abs || n + 1 <= y.abs || used[x + n + (y + n) * (2 * n + 1)]
  return 1 if x.abs == n || y.abs == n
  cnt = 0
  used[x + n + (y + n) * (2 * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], n, used)
  }
  used[x + n + (y + n) * (2 * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
(0..3).each{|n|
  used = Array.new((2 * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, n, used)
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
1
4
92
115516

2015年11月27日金曜日

151127

Ruby


二等辺三角形の形をした領域内のSelf-avoiding walk(3)

左右上下に進む場合について
もう少しだけ調べてみた。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || - m * y > x || m * y > x || used[x + (y + n) * (m * n + 1)]
  return 1 if x == m * n
  cnt = 0
  used[x + (y + n) * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + (y + n) * (m * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
n = 2
puts "n = 2 のとき"
(0..16).each{|m|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
n = 3
puts "n = 3 のとき"
(0..5).each{|m|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
n = 4
puts "n = 4 のとき"
(0..3).each{|m|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
n = 5
puts "n = 5 のとき"
(0..2).each{|m|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
n = 2 のとき
1
3
9
29
95
313
1033
3411
11265
37205
122879
405841
1340401
4427043
14621529
48291629
159496415
n = 3 のとき
1
15
345
10199
317897
9947973
n = 4 のとき
1
153
95265
83918867
n = 5 のとき
1
4169
236871681

2015年11月25日水曜日

151125(2)

Ruby


二等辺三角形の形をした領域内のSelf-avoiding walk(2)

点の集合Dを
格子点で0 ≦ x ≦ mn, x / m ≧ y ≧ - x / mを満たしているもの
とする。
今、Dの点(0, 0)から出発し、
自身が通った点を通ることなく、左右上下に進みながら
(mn, n - 1), (mn, n -2), … , (mn, -(n - 1))
のいずれか一つへたどることとする。
このたどり方が何通りあるか出力してみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || - m * y > x || m * y > x || used[x + (y + n) * (m * n + 1)]
  return 1 if x == m * n
  cnt = 0
  used[x + (y + n) * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + (y + n) * (m * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
m = 0
puts "m = 0 のとき"
(0..7).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 1
puts "m = 1 のとき"
(0..7).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 2
puts "m = 2 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 3
puts "m = 3 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
(4..5).each{|m|
  puts "m = #{m}のとき"
  (0..3).each{|n|
    used = Array.new((m * n + 1) * (2 * n + 1), false)
    p search(move, 0, 0, m, n, used)
  }
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
m = 0 のとき
1
1
1
1
1
1
1
1
m = 1 のとき
1
1
3
15
153
4169
346211
86336255
m = 2 のとき
1
1
9
345
95265
m = 3 のとき
1
1
29
10199
83918867
m = 4のとき
1
1
95
317897
m = 5のとき
1
1
313
9947973

151125

Ruby


二等辺三角形の形をした領域内のSelf-avoiding walk(1)

点の集合Dを
格子点で0 ≦ x ≦ mn, x / m ≧ y ≧ - x / mを満たしているもの
とする。
今、Dの点(0, 0)から出発し、
自身が通った点を通ることなく、右か、上下に進みながら
(mn, n - 1), (mn, n -2), … , (mn, -(n - 1))
のいずれか一つへたどることとする。
このたどり方は
n > 0のとき、
(1 × 3 × 5× … × (2n - 1))^m
通りあるがこれを確かめてみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || - m * y > x || m * y > x || used[x + (y + n) * (m * n + 1)]
  return 1 if x == m * n
  cnt = 0
  used[x + (y + n) * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + (y + n) * (m * n + 1)] = false
  return cnt
end

move = move2
puts "move = #{move2}のとき"
m = 0
puts "m = 0 のとき"
(0..9).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 1
puts "m = 1 のとき"
(0..9).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 2
puts "m = 2 のとき"
(0..6).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 3
puts "m = 3 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (2 * n + 1), false)
  p search(move, 0, 0, m, n, used)
}

出力結果
move = [[1, 0], [0, 1], [0, -1]]のとき
m = 0 のとき
1
1
1
1
1
1
1
1
1
1
m = 1 のとき
1
1
3
15
105
945
10395
135135
2027025
34459425
m = 2 のとき
1
1
9
225
11025
893025
108056025
m = 3 のとき
1
1
27
3375
1157625

2015年11月24日火曜日

151124

Ruby


Dyck path とSelf-avoiding walk の融合(7)

F(3m, 3) について、もう一つ値が求まった。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

# startがx, y、goalがmn, n
def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || y < 0 || m * y > x || used[x + y * (m * n + 1)]
  return 1 if x == m * n && y == n
  cnt = 0
  used[x + y * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + y * (m * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
(0..11).each{|m|
  n = 3
  used = Array.new((m * n + 1) * (n + 1), false)
  p search(move, 0, 0, m, n, used)
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
1
7
44
284
1872
12384
81856
540736
3571712
23592704
155842560
1029427200

2015年11月23日月曜日

151123(6)

C


Dyck path とSelf-avoiding walk の融合(6)

m > 0, (mn + 2) (n + 1) / 2 < 64 のとき、
F(mn, n) を出力してみる。

#include <stdio.h>

int triangular_number(int x, int y, int m, int n){
 int i;
 i = x - m * y + (m * (2 * n - y + 1) + 2) * y / 2;
 return i;
}

// startがx, y、goalの座標がmn, n
long long search(int x, int y, int m, int n, long long used){
 long long cnt = 0LL;
 int z;
 z = triangular_number(x, y, m, n);
 if (x < 0 || m * n + 1 <= x || y < 0 || m * y > x || (used & (1LL << z)) > 0) return 0;
 if (x == m * n && y == n) return 1;
 used += 1LL << z;
 cnt += search(x + 1, y, m, n, used);
 cnt += search(x - 1, y, m, n, used);
 cnt += search(x, y + 1, m, n, used);
 cnt += search(x, y - 1, m, n, used);
 used -= 1LL << z;
 return cnt;
}

int main(void){
 int m;
 int n;
 long long used;
 m = 1;
 for (n = 0; n < 10; n++){
  used = 0LL;
  printf("%lld\n", search(0, 0, m, n, used));
 }
 m = 2;
 for (n = 0; n < 7; n++){
  used = 0LL;
  printf("%lld\n", search(0, 0, m, n, used));
 }
 m = 3;
 for (n = 0; n < 6; n++){
  used = 0LL;
  printf("%lld\n", search(0, 0, m, n, used));
 } 
 m = 4;
 for (n = 0; n < 5; n++){
  used = 0LL;
  printf("%lld\n", search(0, 0, m, n, used));
 }
 m = 5;
 for (n = 0; n < 5; n++){
  used = 0LL;
  printf("%lld\n", search(0, 0, m, n, used));
 }
 return 0;
}

出力結果
1
1
2
7
40
369
5680
150707
6993712
567670347
1
1
4
44
1282
105022
25769912
1
1
8
284
44360
33970104
1
1
16
1872
1592536
1
1
32
12384
57870848

151123(5)

C


Dyck path とSelf-avoiding walk の融合(5)

F(3m, 3) を出力してみる。

#include <stdio.h>

int triangular_number(int x, int y, int m, int n){
 int i;
 i = x - m * y + (m * (2 * n - y + 1) + 2) * y / 2;
 return i;
}

// startがx, y、goalの座標がmn, n
long long search(int x, int y, int m, int n, long long used){
 long long cnt = 0LL;
 int z;
 z = triangular_number(x, y, m, n);
 if (x < 0 || m * n + 1 <= x || y < 0 || m * y > x || (used & (1LL << z)) > 0) return 0;
 if (x == m * n && y == n) return 1;
 used += 1LL << z;
 cnt += search(x + 1, y, m, n, used);
 cnt += search(x - 1, y, m, n, used);
 cnt += search(x, y + 1, m, n, used);
 cnt += search(x, y - 1, m, n, used);
 used -= 1LL << z;
 return cnt;
}

int main(void){
 int m;
 int n = 3;
 long long used;
 for (m = 0; m < 11; m++){
  used = 0LL;
  printf("%lld\n", search(0, 0, m, n, used));
 }
 return 0;
}

出力結果
1
7
44
284
1872
12384
81856
540736
3571712
23592704
155842560

151123(4)

Ruby


Dyck path とSelf-avoiding walk の融合(4)

原点(0, 0) から点(i, j) への最短経路で対角線を超えないものが
i x j 格子点上のDyck path であったが、
x 軸方向もy 軸方向も逆方向に進めるが、自身が通った点は通れないとしてみる。
この条件を満たす経路の個数をF(i, j) と書く。
F(mn, n) を出力してみる。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

# startがx, y、goalがmn, n
def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || y < 0 || m * y > x || used[x + y * (m * n + 1)]
  return 1 if x == m * n && y == n
  cnt = 0
  used[x + y * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + y * (m * n + 1)] = false
  return cnt
end

move = move3
puts "move = #{move3}のとき"
m = 0
puts "m = 0 のとき"
(0..8).each{|n|
  used = Array.new((m * n + 1) * (n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 1
puts "m = 1 のとき"
(0..8).each{|n|
  used = Array.new((m * n + 1) * (n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 2
puts "m = 2 のとき"
(0..6).each{|n|
  used = Array.new((m * n + 1) * (n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 3
puts "m = 3 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (n + 1), false)
  p search(move, 0, 0, m, n, used)
}
m = 4
puts "m = 4 のとき"
(0..4).each{|n|
  used = Array.new((m * n + 1) * (n + 1), false)
  p search(move, 0, 0, m, n, used)
}
(5..8).each{|m|
  puts "m = #{m} のとき"
  (0..3).each{|n|
    used = Array.new((m * n + 1) * (n + 1), false)
    p search(move, 0, 0, m, n, used)
  }
}

出力結果
move = [[1, 0], [-1, 0], [0, 1], [0, -1]]のとき
m = 0 のとき
1
1
1
1
1
1
1
1
1
m = 1 のとき
1
1
2
7
40
369
5680
150707
6993712
m = 2 のとき
1
1
4
44
1282
105022
25769912
m = 3 のとき
1
1
8
284
44360
m = 4 のとき
1
1
16
1872
1592536
m = 5 のとき
1
1
32
12384
m = 6 のとき
1
1
64
81856
m = 7 のとき
1
1
128
540736
m = 8 のとき
1
1
256
3571712

151123(3)

Ruby


Dyck path とSelf-avoiding walk の融合(3)

原点(0, 0) から点(i, j) への最短経路で対角線を超えないものが
i x j 格子点上のDyck path であったが、
y 軸方向だけ逆方向も進めるが、自身が通った点は通れないとしてみる。
この条件を満たす経路の個数をE(i, j) と書く。
このとき、
E(mn, n) = (n!)^m
ちなみに、
D(n, n) = E(n, n) = n!
となっているが、これは対称性より明らか。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

# startがx, y、goalがmn, n
def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || y < 0 || m * y > x || used[x + y * (m * n + 1)]
  return 1 if x == m * n && y == n
  cnt = 0
  used[x + y * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + y * (m * n + 1)] = false
  return cnt
end

M = 3
N = 7
[move2].each{|move|
  puts "move = #{move}のとき"
  (0..M).each{|m|
    puts "m = #{m}のとき"
    (0..N).each{|n|
      # 実行時間の都合上
      break if m > 2 && n > 5
      used = Array.new((m * n + 1) * (n + 1), false)
      p search(move, 0, 0, m, n, used)
    }
  }
}

出力結果
move = [[1, 0], [0, 1], [0, -1]]のとき
m = 0のとき
1
1
1
1
1
1
1
1
m = 1のとき
1
1
2
6
24
120
720
5040
m = 2のとき
1
1
4
36
576
14400
518400
25401600
m = 3のとき
1
1
8
216
13824
1728000

151123(2)

Ruby


Dyck path とSelf-avoiding walk の融合(2)

原点(0, 0) から点(i, j) への最短経路で対角線を超えないものが
i x j 格子点上のDyck path であったが、
x 軸方向だけ逆方向も進めるが、自身が通った点は通れないとしてみる。
この条件を満たす経路の個数をD(i, j) と書く。
このとき、
D(mn, n) = 1 × (m + 1) × (2m + 1) × … × (m(n - 1) + 1)
(ただし、D(0, 0) = 1)
特に、m = 1 のとき、
D(mn, n) = D(n, n) = n!

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

# startがx, y、goalがmn, n
def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || y < 0 || m * y > x || used[x + y * (m * n + 1)]
  return 1 if x == m * n && y == n
  cnt = 0
  used[x + y * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + y * (m * n + 1)] = false
  return cnt
end

M = 3
N = 8
[move1].each{|move|
  puts "move = #{move}のとき"
  (0..M).each{|m|
    puts "m = #{m}のとき"
    (0..N).each{|n|
      used = Array.new((m * n + 1) * (n + 1), false)
      p search(move, 0, 0, m, n, used)
    }
  }
}

出力結果
move = [[1, 0], [-1, 0], [0, 1]]のとき
m = 0のとき
1
1
1
1
1
1
1
1
1
m = 1のとき
1
1
2
6
24
120
720
5040
40320
m = 2のとき
1
1
3
15
105
945
10395
135135
2027025
m = 3のとき
1
1
4
28
280
3640
58240
1106560
24344320

151123

Ruby


Dyck path とSelf-avoiding walk の融合(1)

原点(0, 0) から点(i, j) への最短経路で対角線を超えないものを
i x j 格子点上のDyck path といい、その個数をC(i, j) と書く。
一般に、
C(mn, n) = binomial((m + 1)n, n) / (mn + 1)
が成り立つ。

# -*- coding: cp932 -*-

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

# startがx, y、goalがmn, n
def search(move, x, y, m, n, used)
  return 0 if x < 0 || m * n + 1 <= x || y < 0 || m * y > x || used[x + y * (m * n + 1)]
  return 1 if x == m * n && y == n
  cnt = 0
  used[x + y * (m * n + 1)] = true
  move.each{|mo|
    cnt += search(move, x + mo[0], y + mo[1], m, n, used)
  }
  used[x + y * (m * n + 1)] = false
  return cnt
end

M = 3
N = 10
[move0].each{|move|
  puts "move = #{move}のとき"
  (0..M).each{|m|
    puts "m = #{m}のとき"
    (0..N).each{|n|
      used = Array.new((m * n + 1) * (n + 1), false)
      p search(move, 0, 0, m, n, used)
    }
  }
}

出力結果
move = [[1, 0], [0, 1]]のとき
m = 0のとき
1
1
1
1
1
1
1
1
1
1
1
m = 1のとき
1
1
2
5
14
42
132
429
1430
4862
16796
m = 2のとき
1
1
3
12
55
273
1428
7752
43263
246675
1430715
m = 3のとき
1
1
4
22
140
969
7084
53820
420732
3362260
27343888

2015年11月22日日曜日

151122(3)

Ruby


直角二等辺三角形の形をした領域内のSelf-avoiding walk(4)

move1のときの規則性は偶然(?)
オンライン整数列大辞典の
A208058(http://oeis.org/A208058/list)
と一致している。
このことを確認するとともに、その和が
A054091(http://oeis.org/A054091/list)
と一致することも確認しておく。

=begin
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
=end

# move2のとき規則性は簡単
def search_move2(n)
  ary = [1, n - 1, (n - 2) * (n - 2)]
  if n > 3
    n0 = (1..n - 2).inject(:*)
    ary += (1..n - 4).map{|i| n0 * (i + 1) / (1..i).inject(:*)}.reverse + [n0]
  end
  ary[0..n - 1]
end

# move1のときも規則性は簡単
def search_move1(n)
  search_move2(n).reverse
end

# 0~nまで
def f(n)
  ary = []
  i = 1
  while ary.size < n + 1
    ary += search_move1(i)
    i += 1
  end
  ary[0..n]
end

# 0~nまで
def g(n)
  (1..n + 1).map{|i| search_move1(i).inject(:+)}
end

ary = f(54)
# OEIS A208058のデータ
ary0 =
[1,1,1,1,2,1,2,4,3,1,6,12,9,4,1,24,48,36,16,5,1,
 120,240,180,80,25,6,1,720,1440,1080,480,150,36,7,
 1,5040,10080,7560,3360,1050,252,49,8,1,40320,
 80640,60480,26880,8400,2016,392,64,9,1]
# 一致の確認
p ary == ary0

ary = g(21)
# OEIS A054091のデータ
ary0 =
[1,2,4,10,32,130,652,3914,27400,219202,1972820,
 19728202,217010224,2604122690,33853594972,
 473950329610,7109254944152,113748079106434,
 1933717344809380,34806912206568842,
 661331331924808000,13226626638496160002]
# 一致の確認
p ary == ary0

出力結果
true
true

151122(2)

Ruby


直角二等辺三角形の形をした領域内のSelf-avoiding walk(3)

進むパターンを増やし、さらに分布も調べてみた。

move0 = [[1, 0], [0, 1]]
move1 = [[1, 0], [-1, 0], [0, 1]]
move2 = [[1, 0], [0, 1], [0, -1]]
move3 = [[1, 0], [-1, 0], [0, 1], [0, -1]]

# startがx, y、goalの座標の和がn - 1
def search(move, x, y, n, used)
  ary0 = Array.new(n, 0)
  return ary0 if x < 0 || n <= x || y < 0 || n <= y || used[x + y * n]
  if x + y == n - 1
    ary0[x] = 1
    return ary0
  end
  c_ary = Array.new(n, 0)
  used[x + y * n] = true
  move.each{|m|
    ary1 = search(move, x + m[0], y + m[1], n, used)
    c_ary = (0..n - 1).map{|i| c_ary[i] + ary1[i]}
  }
  used[x + y * n] = false
  return c_ary
end

# move2のとき規則性は簡単
def search_move2(n)
  ary = [1, n - 1, (n - 2) * (n - 2)]
  if n > 3
    n0 = (1..n - 2).inject(:*)
    ary += (1..n - 4).map{|i| n0 * (i + 1) / (1..i).inject(:*)}.reverse + [n0]
  end
  ary[0..n - 1]
end

# move1のときも規則性は簡単
def search_move1(n)
  search_move2(n).reverse
end

N = 10
[move0, move1, move2, move3].each{|move|
  (1..N).each{|n|
    used = Array.new(n * n, false)
    a_ary = search(move, 0, 0, n, used)
    p [a_ary.inject(:+), a_ary]
  }
}

(1..N).each{|n|
  p search_move1(n)
}

(1..N).each{|n|
  p search_move2(n)
}

出力結果
[1, [1]]
[2, [1, 1]]
[4, [1, 2, 1]]
[8, [1, 3, 3, 1]]
[16, [1, 4, 6, 4, 1]]
[32, [1, 5, 10, 10, 5, 1]]
[64, [1, 6, 15, 20, 15, 6, 1]]
[128, [1, 7, 21, 35, 35, 21, 7, 1]]
[256, [1, 8, 28, 56, 70, 56, 28, 8, 1]]
[512, [1, 9, 36, 84, 126, 126, 84, 36, 9, 1]]
[1, [1]]
[2, [1, 1]]
[4, [1, 2, 1]]
[10, [2, 4, 3, 1]]
[32, [6, 12, 9, 4, 1]]
[130, [24, 48, 36, 16, 5, 1]]
[652, [120, 240, 180, 80, 25, 6, 1]]
[3914, [720, 1440, 1080, 480, 150, 36, 7, 1]]
[27400, [5040, 10080, 7560, 3360, 1050, 252, 49, 8, 1]]
[219202, [40320, 80640, 60480, 26880, 8400, 2016, 392, 64, 9, 1]]
[1, [1]]
[2, [1, 1]]
[4, [1, 2, 1]]
[10, [1, 3, 4, 2]]
[32, [1, 4, 9, 12, 6]]
[130, [1, 5, 16, 36, 48, 24]]
[652, [1, 6, 25, 80, 180, 240, 120]]
[3914, [1, 7, 36, 150, 480, 1080, 1440, 720]]
[27400, [1, 8, 49, 252, 1050, 3360, 7560, 10080, 5040]]
[219202, [1, 9, 64, 392, 2016, 8400, 26880, 60480, 80640, 40320]]
[1, [1]]
[2, [1, 1]]
[4, [1, 2, 1]]
[12, [2, 4, 4, 2]]
[48, [6, 12, 12, 12, 6]]
[288, [30, 60, 54, 54, 60, 30]]
[2768, [254, 508, 438, 368, 438, 508, 254]]
[45408, [3750, 7500, 6484, 4970, 4970, 6484, 7500, 3750]]
[1289000, [96834, 193668, 168712, 128582, 113408, 128582, 168712, 193668, 96834]]
[63234984, [4365698, 8731396, 7629762, 5878060, 5012576, 5012576, 5878060, 7629762, 8731396, 4365698]]
[1]
[1, 1]
[1, 2, 1]
[2, 4, 3, 1]
[6, 12, 9, 4, 1]
[24, 48, 36, 16, 5, 1]
[120, 240, 180, 80, 25, 6, 1]
[720, 1440, 1080, 480, 150, 36, 7, 1]
[5040, 10080, 7560, 3360, 1050, 252, 49, 8, 1]
[40320, 80640, 60480, 26880, 8400, 2016, 392, 64, 9, 1]
[1]
[1, 1]
[1, 2, 1]
[1, 3, 4, 2]
[1, 4, 9, 12, 6]
[1, 5, 16, 36, 48, 24]
[1, 6, 25, 80, 180, 240, 120]
[1, 7, 36, 150, 480, 1080, 1440, 720]
[1, 8, 49, 252, 1050, 3360, 7560, 10080, 5040]
[1, 9, 64, 392, 2016, 8400, 26880, 60480, 80640, 40320]

151122

C


直角二等辺三角形の形をした領域内のSelf-avoiding walk(2)

以下の方法は速いが、nが11までしか求まらない。

#include <stdio.h>

int triangular_number(int x, int y, int n){
 int i;
 i = x + (2 * n - y + 1) * y / 2;
 return i;
}

// startがx, y、goalの座標の和がn - 1
long long search(int x, int y, int n, long long used){
 long long cnt = 0LL;
 int z;
 z = triangular_number(x, y, n);
 if (x < 0 || n <= x || y < 0 || n <= y || (used & (1LL << z)) > 0) return 0;
 if (x + y == n - 1) return 1;
 used += 1LL << z;
 cnt += search(x + 1, y, n, used);
 cnt += search(x - 1, y, n, used);
 cnt += search(x, y + 1, n, used);
 cnt += search(x, y - 1, n, used);
 used -= 1LL << z;
 return cnt;
}

int main(void){
 int n;
 long long used;
 for (n = 1; n < 12; n++){
  used = 0LL;
  printf("%lld\n", search(0, 0, n, used));
 }
 return 0;
}

出力結果
1
2
4
12
48
288
2768
45408
1289000
63234984
5356226560

2015年11月21日土曜日

151121(3)

Ruby


直角二等辺三角形の形をした領域内のSelf-avoiding walk(1)

点の集合Dを
格子点でx ≧ 0, y ≧ 0, 0 ≦ x + y ≦ n - 1を満たしているもの
とする。
今、Dの点(0, 0)から出発し、自身が通った点を通ることなく、
(0, n - 1), (1, n -2), … , (n - 1, 0)
のいずれか一つへたどることとする。
このたどり方が何通りあるか出力してみる。

@move0 = [[1, 0], [0, 1]]
@move = [[1, 0], [-1, 0], [0, 1], [0, -1]]

def binomial_theorem(x, y, n, used)
  return 0 if x < 0 || n <= x || y < 0 || n <= y || used[x + y * n]
  return 1 if x + y == n - 1
  cnt = 0
  used[x + y * n] = true
  @move0.each{|m|
    cnt += binomial_theorem(x + m[0], y + m[1], n, used)
  }
  used[x + y * n] = false
  return cnt
end

# startがx, y、goalの座標の和がn - 1
def search(x, y, n, used)
  return 0 if x < 0 || n <= x || y < 0 || n <= y || used[x + y * n]
  return 1 if x + y == n - 1
  cnt = 0
  used[x + y * n] = true
  @move.each{|m|
    cnt += search(x + m[0], y + m[1], n, used)
  }
  used[x + y * n] = false
  return cnt
end

# 確認用
(1..10).each{|n|
  used = Array.new(n * n, false)
  p binomial_theorem(0, 0, n, used)
}

(1..10).each{|n|
  used = Array.new(n * n, false)
  p search(0, 0, n, used)
}

出力結果
1
2
4
8
16
32
64
128
256
512
1
2
4
12
48
288
2768
45408
1289000
63234984

151121(2)

Ruby


Self-avoiding walk(6)

3 x n について、漸化式を利用してみた。
オンライン整数列大辞典の
A006192(http://oeis.org/A006192/list)
と比較し、答え合わせしてみる。

def A006192(n)
  ary = [1, 4, 12, 38]
  (5..n).each{|i| ary[i - 1] = 4 * ary[i - 2] - 3 * ary[i - 3] + 2 * ary[i - 4] + ary[i - 5]}
  ary[0..n - 1]
end
ary = A006192(24)

# OEIS A006192のデータ
ary0 =
[1,4,12,38,125,414,1369,4522,14934,49322,162899,
 538020,1776961,5868904,19383672,64019918,
 211443425,698350194,2306494009,7617832222,
 25159990674,83097804242,274453403399,906458014440]
# 一致の確認
p ary == ary0

151121

C


Self-avoiding walk(5)

m x n を調べた。
(実行時間は1時間半ほどかかる。)
なお、3 ≦ m ≦ 6 のときについて、
オンライン整数列大辞典の
A006192(http://oeis.org/A006192/list)、
A007786(http://oeis.org/A007786/list)、
A007787(http://oeis.org/A007787/list)、
A145403(http://oeis.org/A145403/list)
に詳しく載っている。
また、次のページも詳しい。
http://www.iwriteiam.nl/Crook_path.html

#include <stdio.h>

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

int main(void){
 int i;
 int j;
 long long used;
 for (i = 2; i < 10; i++){
  for (j = 1; j < 7; j++){
   used = 0LL;
   printf("%lld\n", search(0, 0, i, j, used, 1));
  }
 }
 return 0;
}

出力結果
1
2
4
8
16
32
1
4
12
38
125
414
1
8
38
184
976
5382
1
16
125
976
8512
79384
1
32
414
5382
79384
1262816
1
64
1369
29739
752061
20562673
1
128
4522
163496
7110272
336067810
1
256
14934
896476
67005561
5493330332

2015年11月18日水曜日

151118(2)

C


Self-avoiding walk(4)

4 x n を調べた。
(実行時間は6時間弱ほどかかる。)
なお、オンライン整数列大辞典の
A007786(http://oeis.org/A007786/list)
に詳しく載っている。

#include <stdio.h>

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

int main(void){
 int i;
 long long used;
 for (i = 1; i < 16; i++){
  used = 0LL;
  printf("%lld\n", search(0, 0, 4, i, used, 1));
 }
 return 0;
}

出力結果
1
8
38
184
976
5382
29739
163496
896476
4913258
26932712
147657866
809563548
4438573234
24335048679

151118

C


Self-avoiding walk(3)

3 x n を調べた。
(実行時間は1時間弱ほどかかる。)
なお、オンライン整数列大辞典の
A006192(http://oeis.org/A006192/list)
に詳しく載っている。

#include <stdio.h>

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

int main(void){
 int i;
 long long used;
 for (i = 1; i < 21; i++){
  used = 0LL;
  printf("%lld\n", search(0, 0, 3, i, used, 1));
 }
 return 0;
}

出力結果
1
4
12
38
125
414
1369
4522
14934
49322
162899
538020
1776961
5868904
19383672
64019918
211443425
698350194
2306494009
7617832222

2015年11月17日火曜日

151117(2)

C


Self-avoiding walk(2)

2 x n を調べた。
結果が2^(n - 1) (http://mathworld.wolfram.com/Self-AvoidingWalk.html)
となることを確かめよう。

#include <stdio.h>

// startがx, y、goalがw - 1, h - 1
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 (x == w - 1 && y == h - 1) return 1;
 used += 1LL << (x + y * w);
 cnt += search(x + 1, y, w, h, used, depth + 1);
 cnt += search(x - 1, y, w, h, used, depth + 1);
 cnt += search(x, y + 1, w, h, used, depth + 1);
 cnt += search(x, y - 1, w, h, used, depth + 1);
 used -= 1LL << (x + y * w);
 return cnt;
}

int main(void){
 int i;
 long long used;
 for (i = 1; i < 21; i++){
  used = 0LL;
  printf("%d\n", search(0, 0, 2, i, used, 1));
 }
 return 0;
}

出力結果
1
2
4
8
16
32
64
128
256
512
1024
2048
4096
8192
16384
32768
65536
131072
262144
524288

151117

C


Self-avoiding walk(1)

『フカシギの数え方』 おねえさんといっしょ! みんなで数えてみよう!
の動画で有名な問題を解いてみた。
なお、オンライン整数列大辞典の
A007764(http://oeis.org/A007764/list)
に詳しく載っている。

#include <stdio.h>

// startがx, y、goalがw - 1, h - 1
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 (x == w - 1 && y == h - 1) return 1;
 used += 1LL << (x + y * w);
 cnt += search(x + 1, y, w, h, used, depth + 1);
 cnt += search(x - 1, y, w, h, used, depth + 1);
 cnt += search(x, y + 1, w, h, used, depth + 1);
 cnt += search(x, y - 1, w, h, used, depth + 1);
 used -= 1LL << (x + y * w);
 return cnt;
}

int main(void){
 int i;
 long long used;
 for (i = 1; i < 8; i++){
  used = 0LL;
  printf("%d\n", search(0, 0, i, i, used, 1));
 }
 return 0;
}

出力結果
1
2
12
184
8512
1262816
575780564

2015年11月16日月曜日

151116

Ruby


Number of weakly unimodal partitions of n(2)

f(10 ** i) を調べてみた。

def f(n)
  return 1 if n == 0
  s = 0
  ps = Array.new(n + 1){0}
  ps[0] = 1
  (1..n).each{|num|
    old_ary = ps.clone
    # 最大値がnum
    (num..n).each{|i|
      j = ps[i - num]
      # 左のあがる階段がj通り、右のさがる階段がold_ary[n - i]通り
      s += j * old_ary[n - i]
      ps[i] += j
    }
  }
  s
end

(0..5).each{|i| p f(10 ** i)}

出力結果
1
209
881103407093
621820350227357514973418077551453276312298427
19002522919705812331177485379942859319055411162011514367363395916246742323537767948739063616351819551047189156629598790992584402183987248909495666154869
48678030871829155901022590522684849620391551473549000799273927405133763444866764659809512433747622053270638883565798902546892005839120274437671143669750850751160633528202728499786015906676582056664426896989396840534790145116651069254954897615303378028235812765793586406800773111878610378487167931235857873906002076196204137971760798472390952909509209509726649183557835082279874649003714814090330292729697886962484793304921164447358209863494970638027349114178340713474610030285902396235276379

2015年11月15日日曜日

151115(2)

Ruby


Number of weakly unimodal partitions of n(1)

例えば、n = 4 のとき、
[ 1 1 1 1 ],
[ 2 2 ],
[ 1 1 2 ],
[ 1 2 1],
[ 2 1 1 ],
[ 1 3 ],
[ 3 1 ],
[ 4 ]
の8通りある。

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

def f(n)
  return 1 if n == 0
  s = 0
  ps = Array.new(n + 1){0}
  ps[0] = 1
  (1..n).each{|num|
    old_ary = ps.clone
    # 最大値がnum
    (num..n).each{|i|
      j = ps[i - num]
      # 左のあがる階段がj通り、右のさがる階段がold_ary[n - i]通り
      s += j * old_ary[n - i]
      ps[i] += j
    }
  }
  s
end

def A001523(n)
  (0..n).map{|i| f(i)}
end
ary = A001523(38)

# OEIS A001523のデータ
ary0 =
[1,1,2,4,8,15,27,47,79,130,209,330,512,784,1183,
 1765,2604,3804,5504,7898,11240,15880,22277,31048,
 43003,59220,81098,110484,149769,202070,271404,
 362974,483439,641368,847681,1116325,1464999,
 1916184,2498258]
# 一致の確認
p ary == ary0

151115

Ruby


Dixon's identity

Dixon's identity を確認するコードを書いてみた。

def ncr(n, r)
  r = [r, n - r].min
  return 1 if r == 0
  return n if r == 1
  numerator = (n - r + 1..n).to_a
  denominator = (1..r).to_a
  (2..r).each{|p|
    pivot = denominator[p - 1]
    if pivot > 1
      offset = (n - r) % p
      (p - 1).step(r - 1, p){|k|
        numerator[k - offset] /= pivot
        denominator[k] /= pivot
      }
    end
  }
  result = 1
  (0..r - 1).each{|k|
    result *= numerator[k] if numerator[k] > 1
  }
  return result
end

def f(n)
  m = ncr(2 * n, n)
  s = m * m * m
  (0..n - 1).each{|i|
    (n - i) % 2 == 0 ? s0 = 2 : s0 = -2
    m = ncr(2 * n, i)
    s += s0 * m * m * m
  }
  s
end
ary = (0..100).map{|i| f(i)}
p ary

def g(n)
  ncr(3 * n, n) * ncr(2 * n, n)
end
ary0 = (0..100).map{|i| g(i)}

# 一致の確認
p ary == ary0

出力結果
[1, 6, 90, 1680, 34650, 756756, 17153136, 399072960, 9465511770, 227873431500, 5550996791340, 136526995463040, 3384731762521200, 84478098072866400, 2120572665910728000, 53494979785374631680, 1355345464406015082330, 34469858696831179429500, 879619727485803060256500, 22514366432046593564460000, 577831214478475823831865900, 14866378592908813372327325400, 383331414957648741501332688000, 9904298260191196595161087296000, 256376887255990870197659395110000, 6647750135792940867877229051444256, 172644824532516079698894427850821536, 4490186382903298862950669893074864640, 116939573426172903295569869740819625280, 3049327996713402817207903975524583094400, 79607789567531236214574346454361782651136, 2080571532547465690702652742505463614012416, 54432139997018169779222721652071650604875610, 1425432294246982705017331107505766145041177820, 37362109442632851178222433008152866950473778500, 980137754325264739153760087469388189649326812960, 25733153725826742295097099333513427886580242390260, 676129685550027699309454100092387144689767216053800, 17777902978506545489876283429576007693615707520173600, 467765731229343820187277358876812510123932028242832000, 12315686996104586105755778762527877375925475388598463020, 324457176864656573634444032288702901528992030374462542120, 8552867757486695393424970238902882608332272909530900345000, 225584778148789162220956705576428220142319245674420199040000, 5953061113440928872866052731353160003652422655736636347600000, 157178452094435902773657469301001507533471817141538271091744000, 4152006974608028932625878026488081410723506942789897542979264000, 109730270338441253547622797277670526373036819972871986672308480000, 2901279777984880331429724637396325141316257378084815835114447910000, 76743139816667950790962853694727450947519994203318056513168924100000, 2030807663084593981010775419611355697953653094605883738674081337103840, 53761358112015687488073030485490006896692325060431653483330466954103040, 1423761410651421368720750633530303030278157321589316089475330894446493600, 37719286813050009048579857830565358101545002334650667407111133678504507200, 999638712247086659559071129336793852567694712078807811118500826706045376000, 26501496253675160909090667202741659438188047976817213626598113156534880983040, 702813596120264168420100264216585729858529081007646807547722792977161961682240, 18644353293199621014505522244637588198370766848891683084807087573623142344534400, 494751524904085543615600166389508984331795961672670263262354070607072137577852800, 13132732232443982829670053253200410845233538310415033704522124712465804513701120000, 348695928658441817432456030627892908625692498038703169911236614657154552846287571200, 9261011514538879615581805880900871966633347705233238234250972770789560827623546227200, 246028745022037041192303334172111614264463441118214955380257372360309560800446550656000, 6537707150411741663005756625967313613818198591044569139341139971805776961073695340032000, 173769574576495682336776543913345076449313300929448789599626442561155990662132131755196250, 4619844079024722843588346052452165457336228601231068577649145511745421966180983155392703900, 122852217341117492036248772573886372340679379304913236859634164255863796803154353716269340900, 3267666462581739485214547437140127613605951857653517594350705995132576818511167260798075140000, 86933778740026044115943597885403161498827202773533829459739435577253316405624302354856749042500, 2313304017763894040968051757215233465718634387603234794004282523268687808323575481062130127325000, 61569766689937960254467232238567732176563555570135401815278063174410038614352534681265053609228000, 1639043128490957055010609938005720172022638941338998290419173768250967971015999354113676876441968000, 43641420336729707234397085212756472867338899126740342211913301549321201589378547617343908035355872500, 1162226138365172938879933484560700528501970704873097370244006742217407822698709836884633653698936895000, 30957323692831286498361413085833414734714654312050763773519142333789464028274984807887630520183881575000, 824736124322297825688612963447227057508515419837631867741241705523974030474875755265921557196920134631680, 21975733717096185796473626515288941743924095257190195678078868755258439600807864956709993625683724072130080, 585661272497616515044830110588768035917936632077689799138681405299103734411140283658246937560176292905595200, 15610703030636665846776793125239864373560149202954437762446036511266445496710354622007000225223318387231190400, 416168686146022827044820195509919665834471659091100681768269757209503477286119019709906216547453246798193472000, 11096487752511595522482872981716523639998320427696869834535650154495770372651104099596352535774485430175706753580, 295916487746676394175695527224358066635565635109344692496473794380659496878718652627644477636689318362463451889480, 7892573306437975093912735849556272407855791344513664061915171071419909900391047466127567016562144907067845845011800, 210538889733995326906395070634527081992574112904425463395271760408287705410808827225334698141906262923797923536961600, 5617055837928659338170191616631196596528072102360840913287968947967709827605550131118091967859106291865782186200825000, 149880931040957026874529904637212715163773136100448038243478807142270298522458628038440481400094901810883970179780242000, 3999851087769055453811143112178812448372753719708224985984974721269975992586370222784644193729736525125072268551149952000, 106757816983895439658994204191090459878486419137991228916459848486056354871488771622004732463638706936949015681424668160000, 2849788534514062176579024264201167056451569678637302538630261745782618487386933870371587588509615989229368107390592013520000, 76081906500907553147477119195886983641940845845402142134091838326051045515429226810690522325732671789210260531484470786400000, 2031443260542010119077007584099957194338814614328418974344795855037656286406786599717148568706192509121387445272510543678944000, 54247704004735895386090224503988423394192501938670360863508124163700900923904258122817457540081220468188278480591944401278464000, 1448805941667125372403791971777715395588045257943813217532610050473131929802523339870378813437542519711006657672236416813444800000, 38698144413464674445829484765423084201405831818668498310550180210452305718478360149170444296764273074181678349187734135669833600000, 1033760752303553842937255724765783702455843429433964728522211921330371462402629258670234458104146561388619664728641100654539104000000, 27618421889465228630876382419127612543085112945307461438488825443046409476861202917842908764569462858257996574040003060988349096960000, 737947690602566698072784902905870017174971432544507503188163987635460736549388430306764855863733935654795857679150889079441950057510000, 19719650106090269411886524406166744438123399647642206560378368733253681899310525351475171400767204106885463120020077377111711074817700000, 527009907687419540754903211816118904773454471132968541590919944400953355449003948449857277498137936662440446489932655167483200557195900000, 14085860158942600030216758021360117160948284181001167141634260650471762166454553039069953312255097835740526411846011462878675565550764000000, 376523493564631064367712071965768747782444205128669798396168767743500485766630075466163294008566118208045715304490994009624725072511252178400]
true

2015年11月14日土曜日

151114(2)

C


Number of directed Hamiltonian paths in mxn grid graph(2)

先程のコードを使って、
m = 3, 4, 5, 6 のときも調べる。
以下、結果のみ記す。

m = 3 のとき
2
16
40
124
264
672
1376
3156
6380
13804
27756
58076
116476
239036
478780
971644
1944828

m = 4 のとき
2
28
124
552
2012
7220
24020
77968
244376
750244
2256892
6685588

m = 5 のとき
2
44
264
2012
8648
53992
219444
1205608
4869340

m = 6 のとき
2
64
672
7220
53992
458696
3240068

151114

C


Number of directed Hamiltonian paths in mxn grid graph(1)

m = 2 のときを調べる。
なお、オンライン整数列大辞典の
A137882 のCOMMENTSには
n > 1 のとき、2 * (n^2 - n + 2)
となることが載っている。

#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 + 1, y, w, h, used, depth + 1);
 cnt += search(x - 1, y, w, h, used, depth + 1);
 cnt += search(x, y + 1, w, h, used, depth + 1);
 cnt += search(x, y - 1, w, h, used, depth + 1);
 used -= 1LL << (x + y * w);
 return cnt;
}

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

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

出力結果
2
8
16
28
44
64
88
116
148
184
224
268
316
368
424
484
548
616
688
764
844
928
1016
1108
1204
1304
1408
1516
1628
1744

2015年11月13日金曜日

151113

Ruby


Number of smooth weakly unimodal compositions of n into positive parts such that the first and last part are 1(1)

例えば、n = 9 のとき、
[ 1 1 1 1 1 1 1 1 1 ],
[ 1 1 1 1 1 1 2 1 ],
[ 1 1 1 1 1 2 1 1 ],
[ 1 1 1 1 2 1 1 1 ],
[ 1 1 1 1 2 2 1 ],
[ 1 1 1 2 1 1 1 1 ],
[ 1 1 1 2 2 1 1 ],
1 2 3 2 1 ],
[ 1 1 2 1 1 1 1 1 ],
[ 1 1 2 2 1 1 1 ],
[ 1 1 2 2 2 1 ],
[ 1 2 1 1 1 1 1 1 ],
[ 1 2 2 1 1 1 1 ],
[ 1 2 2 2 1 1 ]
の14通りある。

オンライン整数列大辞典の
A001522(http://oeis.org/A001522/list)
と比較し、答え合わせしてみる。
(なお、以下のコードにおいて、
「和がn でk のぼるのぼりかた」 = 「n をk 個の相異なる数に分割する場合の数」
を利用したが、
これは「縦横逆に見る」だけのことである。)

# nをk個の相異なる数(最小値がl)に分割
def g(n, l, k)
  return 1 if n == l && k == 1
  # 先頭がiのものにlを前から追加
  (l + 1..n - l).inject(0){|s, i| s + g(n - l, i, k - 1)}
end

# nをk個の相異なる数に分割
def g0(n, k)
  (1..n).inject(0){|s, l| s + g(n, l, k)}
end

def f(n)
  # 1, 1, … , 1の真っ平らなもの
  s = 1
  # 山型
  (1..n).each{|i|
    (1..i).each{|k|
      # 和がiでkのぼるものと和がn - iでk + 1のぼるものを組み合わせたもの
      s += g0(i, k) * g0(n - i, k + 1) 
    }
  }
  s
end

def A001522(n)
  (0..n).map{|i| f(i)}
end
ary = A001522(56)

# OEIS A001522のデータ
ary0 =
[1,1,1,1,2,3,5,7,10,14,19,26,35,47,62,82,107,139,
 179,230,293,372,470,591,740,924,1148,1422,1756,
 2161,2651,3244,3957,4815,5844,7075,8545,10299,
 12383,14859,17794,21267,25368,30207,35902,42600,
 50462,59678,70465,83079,97800,114967,134956,
 158205,185209,216546,252859]
# 一致の確認
p ary == ary0

2015年11月9日月曜日

151109(4)

Ruby


Number of partitions of n into fourth powers

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

def A046042(n)
  fourth_powers = (1..(n ** (1.0 / 4)).to_i).map{|i| i * i * i * i}
  ps = Array.new(n + 1){0}
  ps[0] = 1
  fourth_powers.each{|num|
    (num..n).each{|i|
      ps[i] += ps[i - num]
    }
  }
  ps[1..-1]
end
ary = A046042(102)

# OEIS A046042のデータ
ary0 =
[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,
 2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,4,4,
 4,4,4,4,4,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,5,5,5,5,
 5,5,5,5,5,6,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,8,9,9,9,
 9,9,9]
# 一致の確認
p ary == ary0

151109(3)

Ruby


Number of partitions of n into cubes

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

def A003108(n)
  cubes = (1..(n ** (1.0 / 3)).to_i).map{|i| i * i * i}
  ps = Array.new(n + 1){0}
  ps[0] = 1
  cubes.each{|num|
    (num..n).each{|i|
      ps[i] += ps[i - num]
    }
  }
  ps
end
ary = A003108(86)

# OEIS A003108のデータ
ary0 =
[1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,
 4,4,4,5,5,5,5,5,6,6,6,7,7,7,7,7,8,8,8,9,9,9,9,9,
 10,10,10,11,11,11,12,12,13,13,13,14,14,14,15,15,
 17,17,17,18,18,18,19,19,21,21,21,22,22,22,23,23,
 25,26,26,27,27,27,28]
# 一致の確認
p ary == ary0

151109(2)

Ruby


Number of partitions of n into squares

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

def A001156(n)
  squares = (1..Math.sqrt(n).to_i).map{|i| i * i}
  ps = Array.new(n + 1){0}
  ps[0] = 1
  squares.each{|num|
    (num..n).each{|i|
      ps[i] += ps[i - num]
    }
  }
  ps
end
ary = A001156(67)

# OEIS A001156のデータ
ary0 =
[1,1,1,1,2,2,2,2,3,4,4,4,5,6,6,6,8,9,10,10,12,13,
 14,14,16,19,20,21,23,26,27,28,31,34,37,38,43,46,
 49,50,55,60,63,66,71,78,81,84,90,98,104,107,116,
 124,132,135,144,154,163,169,178,192,201,209,220,
 235,247,256]
# 一致の確認
p ary == ary0

151109

Ruby


Number of palindromic and unimodal compositions of n

線対称な山型に分割してみる。
例えば、n = 6 のとき、
1 + 2 + 2 + 1,
1 + 4 + 1,
などがある。
この山を横に切ると、
4 + 2,
3 + 1 + 1 + 1,
…,
となる。
この操作からわかるように、
p(n | 和因子が線対称な山型に並ぶ)
= p(n | 和因子は偶数) + p(n | 和因子は奇数)
が成り立つ。

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

def A035363(n)
  even = []
  2.step(n, 2){|i| even << i}
  ps0 = Array.new(n + 1){0}
  ps0[0] = 1
  even.each{|num|
    (num..n).each{|i|
      ps0[i] += ps0[i - num]
    }
  }
  ps0
end

def A000009(n)
  odd = []
  1.step(n, 2){|i| odd << i}
  ps1 = Array.new(n + 1){0}
  ps1[0] = 1
  odd.each{|num|
    (num..n).each{|i|
      ps1[i] += ps1[i - num]
    }
  }
  ps1
end

def A096441(n)
  ps0, ps1 = A035363(n), A000009(n)
  (1..n).map{|i| ps0[i] + ps1[i]}
end

ary = A035363(69)
# OEIS A035363のデータ
ary0 =
[1,0,1,0,2,0,3,0,5,0,7,0,11,0,15,0,22,0,30,0,42,0,
 56,0,77,0,101,0,135,0,176,0,231,0,297,0,385,0,490,
 0,627,0,792,0,1002,0,1255,0,1575,0,1958,0,2436,0,
 3010,0,3718,0,4565,0,5604,0,6842,0,8349,0,10143,0,
 12310,0]
# 一致の確認
p ary == ary0

ary = A000009(55)
# OEIS A000009のデータ
ary0 =
[1,1,1,2,2,3,4,5,6,8,10,12,15,18,22,27,32,38,46,
 54,64,76,89,104,122,142,165,192,222,256,296,340,
 390,448,512,585,668,760,864,982,1113,1260,1426,
 1610,1816,2048,2304,2590,2910,3264,3658,4097,4582,
 5120,5718,6378]
# 一致の確認
p ary == ary0

ary = A096441(55)
# OEIS A096441のデータ
ary0 =
[1,2,2,4,3,7,5,11,8,17,12,26,18,37,27,54,38,76,54,
 106,76,145,104,199,142,266,192,357,256,472,340,
 621,448,809,585,1053,760,1354,982,1740,1260,2218,
 1610,2818,2048,3559,2590,4485,3264,5616,4097,7018,
 5120,8728,6378]
# 一致の確認
p ary == ary0

出力結果
true
true
true

2015年11月8日日曜日

151108(5)

Ruby


Polite number(5)

f(10 ** i) を調べてみた。

require 'prime'

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

def f(n)
  return 0 if n == 1
  factor = Prime.prime_division(n)
  o_factor = factor.clone
  n1 = 1
  if o_factor[0][0] == 2
    n1 = power(2, o_factor[0][1])
    o_factor.shift
  end
  a0, a1 = 1, 1
  o_factor.each{|i|
    p, q = i[0], i[1]
    a0 *= q + 1
    a1 *= power(p, q + 1) / (p - 1)
  }
  ((2 * n1 + 1) * a1 - a0) / 2 - n
end

(0..25).each{|i| p f(10 ** i)}

出力結果
0
4
38
324
2884
26942
259746
2548792
25244072
251220570
2506103254
25030517060
250152586860
2500762937398
25003814693162
250019073478128
2500095367415248
25000476837125426
250002384185725470
2500011920928823996
25000059604644513236
250000298023223352654
2500001490116118336178
25000007450580594826664
250000037252902980424824
2500000186264514914707082

n = 10 ** i のとき
f(n) = (2^i + 1 / 2) (5^(i + 1) - 1) / (5 - 1) - (i + 1) / 2 - n
より、
f(n) ≒ 2^i 5^(i + 1) / 4 - n = n / 4
となっている。

151108(4)

Ruby


Polite number(4)

末項の和(nは含まない)を計算してみた。

require 'prime'

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

def f(n)
  return 0 if n == 1
  factor = Prime.prime_division(n)
  o_factor = factor.clone
  n1 = 1
  if o_factor[0][0] == 2
    n1 = power(2, o_factor[0][1])
    o_factor.shift
  end
  a0, a1 = 1, 1
  o_factor.each{|i|
    p, q = i[0], i[1]
    a0 *= q + 1
    a1 *= power(p, q + 1) / (p - 1)
  }
  ((2 * n1 + 1) * a1 - a0) / 2 - n
end

def A(n)
  (1..n).map{|i| f(i)}
end
p A(100)

出力結果
[0, 0, 2, 0, 3, 3, 4, 0, 9, 4, 6, 5, 7, 5, 19, 0, 9, 13, 10, 6, 25, 7, 12, 9, 20, 8, 31, 7, 15, 28, 16, 0, 37, 10, 35, 21, 19, 11, 43, 10, 21, 36, 22, 9, 69, 13, 24, 17, 35, 26, 55, 10, 27, 44, 51, 11, 61, 16, 30, 46, 31, 17, 90, 0, 59, 52, 34, 12, 73, 48, 36, 37, 37, 20, 108, 13, 65, 60, 40, 18, 98, 22, 42, 58, 75, 23, 91, 13, 45, 102, 75, 15, 97, 25, 83, 33, 49, 43, 132, 38]

151108(3)

Polite number(3)

n を連続する正の整数に分割しよう。

n = 90 (= 2^1 3^2 5^1)のとき、
n の奇数の約数は1, 3, 5, 9, 15, 45 の6個である。

90
= 29 + 30 + 31
= 16 + 17 + 18 + 19 +20
= 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13 + 14
= 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12 + 13
= 21 + 22 + 23 + 24

さて初項の和(90も含む)は
|90 / 1 - 1 / 2| + 1 / 2
+ |90 / 3 - 3 / 2| + 1 / 2
+ |90 / 5 - 5 / 2| + 1 / 2
+ |90 / 9 - 9 / 2| + 1 / 2
+ |90 / 15 - 15 / 2| + 1 / 2
+ |90 / 45 - 45 / 2| + 1 / 2
となる。

また、末項の和(90も含む)は
90 / 1 + 1 / 2 - 1 / 2
+ 90 / 3 + 3 / 2 - 1 / 2
+ 90 / 5 + 5 / 2 - 1 / 2
+ 90 / 9 + 9 / 2 - 1 / 2
+ 90 / 15 + 15 / 2 - 1 / 2
+ 90 / 45 + 45 / 2 - 1 / 2
= 90(1 / 1 + 1 / 3 + … + 1 / 45)
 + (1 + 3 + … + 45) / 2
 - (2 + 1)(1 + 1) / 2
= 2^1 3^2 5^1 (1 + 1 / 3 + 1 / 3^2)(1 + 1 / 5^1)
 + (1 + 3 + 3^2)(1 + 5^1) / 2
 - (2 + 1)(1 + 1) / 2
= 2^1 (1 + 3 + 3^2)(1 + 5^1)
 + (1 + 3 + 3^2)(1 + 5^1) / 2
 - (2 + 1)(1 + 1) / 2
= (2^1 + 1 / 2)(1 + 3 + 3^2)(1 + 5^1)
 - (2 + 1)(1 + 1) / 2
= (2^1 + 1 / 2) (3^(+ 1) - 1) / (3 - 1) (5^(1 + 1) - 1) / (5 - 1)
 - (2 + 1)(1 + 1) / 2

一般に
n = 2^q1 3^q2 5^q3 …
と表されるとき、
末項の和(nも含む)は
(2^q1 + 1 / 2) (3^(q2 + 1) - 1) / (3 - 1) (5^(q3 + 1) - 1) / (5 - 1) …
- (q2 + 1)(q3 + 1) … / 2
となる。

151108(2)

Ruby


Polite number(2)

politeness について、
オンライン整数列大辞典の
A069283(http://oeis.org/A069283/list)
と比較し、答え合わせしてみる。

require 'prime'

def politeness(n)
  return 0 if n < 2
  factor = Prime.prime_division(n)
  o_factor = factor.clone
  o_factor.shift if o_factor[0][0] == 2
  # 奇数の約数の個数 - 1
  o_factor.inject(1){|a, n| a * (n[1] + 1)} - 1
end

def A069283(n)
  (0..n).map{|i| politeness(i)}
end
ary = A069283(100)

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

151108

Ruby


約数の出力

Polite number を考えていて、約数の出力について気になった。
haruyaさんのコード(http://d.hatena.ne.jp/haruya12/20150918/1442570781)
がシンプルで勉強になった。

require 'prime'

def d(m, k)
  return m unless @factor[k]
  (0..@factor[k][1]).map{|e| d(m * @factor[k][0] ** e, k + 1)}.flatten
end

n = 630
# 順番を逆にする方が後の結果が綺麗になる
@factor = Prime.prime_division(n).reverse
p @factor
p d(1, 0)
p d(1, n % 2 == 0 ? 1 : 0)
p d(1, n % 3 == 0 ? 2 : 0)
p d(1, n % 5 == 0 ? 3 : 0)
p d(1, n % 7 == 0 ? 4 : 0)

出力結果
[[7, 1], [5, 1], [3, 2], [2, 1]]
[1, 2, 3, 6, 9, 18, 5, 10, 15, 30, 45, 90, 7, 14, 21, 42, 63, 126, 35, 70, 105, 210, 315, 630]
[1, 2, 3, 6, 9, 18, 5, 10, 15, 30, 45, 90]
[1, 2, 3, 6, 9, 18]
[1, 2]
1

2015年11月7日土曜日

151107(3)

Ruby


Polite number(1)

2の冪を除けばよく、
f(n) = n + log(n + log(n, 2), 2).to_i
(https://en.wikipedia.org/wiki/Polite_number)
と表される。
オンライン整数列大辞典の
A138591(http://oeis.org/A138591/list)
と比較し、答え合わせしてみる。

include Math

def f(n)
  n + log(n + log(n, 2), 2).to_i
end

def A138591(n)
  (1..n).map{|i| f(i)}
end
ary = A138591(71)

# OEIS A138591のデータ
ary0 =
[1,3,5,6,7,9,10,11,12,13,14,15,17,18,19,20,21,22,
 23,24,25,26,27,28,29,30,31,33,34,35,36,37,38,39,
 40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,
 56,57,58,59,60,61,62,63,65,66,67,68,69,70,71,72,
 73,74,75,76,77]
# 一致の確認
p ary == ary0

151107(2)

Ruby


(-1)^k (n / k) の和

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

def s(n)
  (1..n).inject(0){|s, k| s - (n / k) * (-1) ** k}
end

def A059851(n)
  (0..n).map{|i| s(i)}
end
ary = A059851(74)

# OEIS A059851のデータ
ary0 =
[0,1,1,3,2,4,4,6,4,7,7,9,7,9,9,13,10,12,12,14,12,
 16,16,18,14,17,17,21,19,21,21,23,19,23,23,27,24,
 26,26,30,26,28,28,30,28,34,34,36,30,33,33,37,35,
 37,37,41,37,41,41,43,39,41,41,47,42,46,46,48,46,
 50,50,52,46,48,48]
# 一致の確認
p ary == ary0

151107

Ruby


n / k の和

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

def s(n)
  (1..n).inject(0){|s, k| s + n / k}
end

def A006218(n)
  (0..n).map{|i| s(i)}
end
ary = A006218(59)

# OEIS A006218のデータ
ary0 =
[0,1,3,5,8,10,14,16,20,23,27,29,35,37,41,45,50,52,
 58,60,66,70,74,76,84,87,91,95,101,103,111,113,119,
 123,127,131,140,142,146,150,158,160,168,170,176,
 182,186,188,198,201,207,211,217,219,227,231,239,
 243,247,249]
# 一致の確認
p ary == ary0

2015年11月4日水曜日

151104

Ruby


Kolakoski sequence(2)

短くしてみた。

def A000002(n)
  ary = [1, 2, 2]
  i = 2
  while ary.size < n
    ary += [1 + i % 2] * ary[i]
    i += 1
  end
  ary[0..n - 1]
end
ary = A000002(108)

# OEIS A000002のデータ
ary0 =
[1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,2,1,2,1,1,2,
 1,2,2,1,1,2,1,1,2,1,2,2,1,2,2,1,1,2,1,2,2,1,2,1,1,
 2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,1,2,
 2,1,2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,
 2,1,2,1,1,2,1,2,2]
# 一致の確認
p ary == ary0

2015年11月3日火曜日

151103(2)

Ruby


Kolakoski sequence(1)

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

def A000002(n)
  ary = [1, 2]
  b_ary = [2]
  # 次の文字
  str = 1
  ary += b_ary
  while ary.size < n
    f_ary, b_ary = b_ary, []
    (0..f_ary.size - 1).each{|i|
      b_ary += [str] * f_ary[i]
      # 次の文字
      str = 1 + str % 2
    }
    ary += b_ary
  end
  ary[0..n - 1]
end
ary = A000002(108)

# OEIS A000002のデータ
ary0 =
[1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,2,1,2,1,1,2,
 1,2,2,1,1,2,1,1,2,1,2,2,1,2,2,1,1,2,1,2,2,1,2,1,1,
 2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,1,2,
 2,1,2,1,1,2,2,1,2,2,1,1,2,1,2,2,1,2,2,1,1,2,1,1,2,
 2,1,2,1,1,2,1,2,2]
# 一致の確認
p ary == ary0

151103

Ruby


Gauss circle problem(4)

f(10 ** i) を調べてみた。
iが9でも計算が合ったが、速度がかなり遅くなった。
(実行時間は7時間弱ほどかかる。)

require 'bigdecimal/math'
include BigMath

def f(r)
  q = r * r
  m = (BigDecimal.new('2').sqrt(12) * r / 2).to_i
  # 正方形の内部またはx軸,y軸上の格子点
  s = (2 * m + 1) * (2 * m + 1) + (r - m) * 4
  # 上記以外
  (m + 1..r).each{|x| s += 8 * (BigDecimal.new((q - x * x).to_s).sqrt(12)).to_i}
  s
end

(0..9).each{|i| p f(10 ** i)}

出力結果
5
317
31417
3141549
314159053
31415925457
3141592649625
314159265350589
31415926535867961
3141592653589764829

2015年11月2日月曜日

151102

Ruby


Gauss circle problem(3)

f(10 ** i) を調べてみた。
計算は速いが、iが9で合わなくなった。
Math.sqrtの計算の精度の問題だろう。

def f(r)
  q = r * r
  m = (Math.sqrt(2) * r / 2).to_i
  # 正方形の内部またはx軸,y軸上の格子点
  s = (2 * m + 1) * (2 * m + 1) + (r - m) * 4
  # 上記以外
  (m + 1..r).each{|x| s += 8 * Math.sqrt(q - x * x).to_i}
  s
end

# f(10 ** 9)の計算結果は間違い
(0..9).each{|i| p f(10 ** i)}

出力結果
5
317
31417
3141549
314159053
31415925457
3141592649625
314159265350589
31415926535867961
3141592653589765277

2015年11月1日日曜日

151101

Ruby


Gauss circle problem(2)

Gauss circle problem(1)のコードはnを大きくすると遅くなるので、
別の方法で求めるようにしてみた。
オンライン整数列大辞典の
A000328(http://oeis.org/A000328/list)
と比較し、答え合わせしてみる。

def f(r)
  q = r * r
  m = (Math.sqrt(2) * r / 2).to_i
  # 正方形の内部またはx軸,y軸上の格子点
  s = (2 * m + 1) * (2 * m + 1) + (r - m) * 4
  # 上記以外
  (m + 1..r).each{|x| s += 8 * Math.sqrt(q - x * x).to_i}
  s
end

def A000328(n)
  (0..n).map{|i| f(i)}
end
ary = A000328(46)

# OEIS A000328のデータ
ary0 =
[1,5,13,29,49,81,113,149,197,253,317,377,441,529,
 613,709,797,901,1009,1129,1257,1373,1517,1653,
 1793,1961,2121,2289,2453,2629,2821,3001,3209,3409,
 3625,3853,4053,4293,4513,4777,5025,5261,5525,5789,
 6077,6361,6625]
# 一致の確認
p ary == ary0