算額（その382）
東京都渋谷区 金王八幡宮（金王八幡神社） 嘉永3年(1850)
https://www.naruto-u.ac.jp/journal/info-edu/j05006.pdf
二十八宿のうち15宿の名前をつけた円がある。角と亢の円周の和が16寸，心，尾，箕の円周の和が30寸，虚，危，室，壁，奎 の周の和が63寸であるとき，角の円周を求めよ。
答は 7 寸 7 分 6 厘 3 毛 3 糸 2 忽 1 微と残りは微少数（7.763321...）である。
include("julia-source.txt");
using Plots
function draw()
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
f(n) = -69//10960*n^2 + 1079//2192*n + 997//137
names = ["角", "亢", "氐", "房", "心", "尾", "箕", "斗", "牛", "女", "虚", "危", "室", "壁", "奎"]
plot()
for (i, name) in enumerate(names)
point(i, 1, name, :blue, :center, mark=false)
circle(i, 0, f(i)/30)
end
plot!(showaxis=false)
end;
金王八幡宮といえば，冲方丁の「天地明察」にも出てくる。その中に，この問題を参考にしたと思われる算額問題がある。これについては，末尾に記す。
15宿の名前を順番に書いておく。角が一番小さく，奎が一番大きい。
角，亢，氐，房，心，尾，箕，斗，牛，女，虚，危，室，壁，奎
天文（学）に，しかも昔の天文（学）に興味がなければ，読み方も分かりづらいので，この算額の題意を，現代数学風に書き下す。
数列 a[n], n = 1, 2, ..., 15 がある。
a[i] < a[i+1], i = 1, 2, ..., 14 である。
a[1] + a[2] = 16
a[5] + a[6] + a[7] = 30
a[11] + a[12] + a[13] + a[14] + a[15] = 63
がわかっているとき，a[1] を求めよ。
等差数列ではなさそうだが，次に簡単な形式の数列といえば a[n] = p×n^2 + q×n + r であろうか。
第1の条件 a[1] + a[2] = 16 は
(p×1^2 + q×1 + r) + (p×2^2 + q×2 + r)
= 5p + 3q + 2r
のようになる。
2, 3 番目の条件における p, q, r の係数を求める。
using SymPy
function coef(s, e)
coef2 = coef1 = coef0 = 0
for n in s:e
coef2 += n^2
coef1 += n
coef0 += 1
end
return (coef2, coef1, coef0)
end
coef(1, 2) |> println
coef(5, 7) |> println
coef(11, 15) |> println
(5, 3, 2)
(110, 18, 3)
(855, 65, 5)
以下の3元連立方程式を解けばよい。
@syms p, q, r
eq1 = 5p + 3q + 2r - 16
eq2 = 110p + 18q + 3r - 30
eq3 = 855p + 65q + 5r - 63
res = solve([eq1, eq2, eq3], (p, q, r))
Dict{Any, Any} with 3 entries:
p => -69/10960
q => 1079/2192
r => 997/137
分数形式を小数点付きで表す。
res[p].evalf(), res[q].evalf(), res[r].evalf()
(-0.00629562043795620, 0.492244525547445, 7.27737226277372)
a[1] は 1^2p + 1q + r なので，7.76332116788321 が答えである。
(res[p] + res[q] + res[r]).evalf() |> println
7.76332116788321
行列演算で行うならば，以下のようにする。
A = Sym[5 3 2; 110 18 3; 855 65 5]
b = [16, 30, 63];
pqr = A \ [16.0, 30, 63]
pqr |> println
Sym[-0.00629562043795620, 0.492244525547445, 7.27737226277372]
pqr' * [1;1;1] |> println
7.76332116788321
演算子 \ は Julia 特有のものなので，一般言語的には以下のようにする。
inv(A) * b |> println
Sym[-69/10960, 1079/2192, 997/137]
---------------------------------
さて，冲方丁の天地明察に出てくる問題である。
冲方 丁. 天地明察（特別合本版） (Japanese Edition) (p.237). Kindle 版.
『今、図の如く、大小の十四宿の星の名を持つ円が並んでいる。角星と亢星の周の長さを足すと九寸である。また房星と心星と尾星の周の長さを足すと十八寸である。さらに女星、虚星、危星、室星、壁星の五つの星の周の長さを足すと四十五寸である。角星の周の長さは何寸であるか問う』
変更点は以下の通り。
- 15宿から14宿に縮小された
- 条件数は変わらず3個であるが，係数行列と定数ベクトルは異なる
- 2 番目の条件式で，金王八幡宮のものは「心，尾，箕の円周の和」であるが天地明察のものでは「房星と心星と尾星の周の長さを足す」となっている
- p*n^2+q*n+r だが，（結果的に）p = 0 である。つまり，簡単な階差数列である
coef(1, 2) |> println
coef(4, 6) |> println
coef(10, 14) |> println
(5, 3, 2)
(77, 15, 3)
(730, 60, 5)
以下の3元連立方程式を解けばよい。
@syms p, q, r
eq1 = 5p + 3q + 2r - 9
eq2 = 77p + 15q + 3r - 18
eq3 = 730p + 60q + 5r - 45
$730 p + 60 q + 5 r - 45$
res = solve([eq1, eq2, eq3])
Dict{Any, Any} with 3 entries:
p => 0
q => 3/7
r => 27/7
分数形式を小数点付きで表す。なお，p = 0 なので，一般項は
a[n] = 27//7 + 3/7*n
つまり，階差数列である。これを村瀬は「……招差術か」と言っている。
res[p].evalf(), res[q].evalf(), res[r].evalf()
(0, 0.428571428571429, 3.85714285714286)
a[1] は 1^2p + 1q + r なので，7.76332116788321 が答えである。
res[r] + res[q]
30/7
1^2*res[p] + 1*res[q] + 1*res[r] |> println
30/7
因縁の "30/7" が出てきた。
一般項 a(n) を求める関数を定義しておこう。
func_a(n) = (3n + 27)//7;
for i in 1:14
println("func_a($i) = $(func_a(i))")
end
func_a(1) = 30//7
func_a(2) = 33//7
func_a(3) = 36//7
func_a(4) = 39//7
func_a(5) = 6//1
func_a(6) = 45//7
func_a(7) = 48//7
func_a(8) = 51//7
func_a(9) = 54//7
func_a(10) = 57//7
func_a(11) = 60//7
func_a(12) = 9//1
func_a(13) = 66//7
func_a(14) = 69//7
func_a(1) + func_a(2)
9//1
func_a(4) + func_a(5) + func_a(6)
18//1
func_a(10) + func_a(11) + func_a(12) + func_a(13) + func_a(14)
45//1
角星の周の長さは 30/7
func_a(1) |> println
30//7
行列演算で行うならば，以下のようにする。
A = Sym[5 3 2; 77 15 3; 730 60 5]
b = [9, 18, 45];
A, b
(Sym[5 3 2; 77 15 3; 730 60 5], [9, 18, 45])
pqr = A \ [9, 18, 45]
pqr |> println
Sym[0, 3/7, 27/7]
pqr' * [1;1;1] |> println
30/7
演算子 \ は Julia 特有のものなので，一般言語的には以下のようにする。
inv(A) * b |> println
Sym[0, 3/7, 27/7]
---
なお，この数列の問題は，初版本（単行本）では解答不能ではないものの，なぜこのような答えになる問題を出したのかという状態であった。
236ページ
『今、図の如く、大小の十五宿の星の名を持つ円が並んでいる。角星と亢星の周の長さを足すと十寸である。また心星と尾星と箕星の周の長さを足すと二十七寸五分である。さらに虚星、危星、室星、壁星、奎星の五つの星の周の長さを足すと四十寸である。角星の周の長さは何寸であるか問う』
角星 + 亢星 = 10寸
心星 + 尾星 + 箕星 = 27寸5分
虚星 + 危星 + 室星 + 壁星 + 奎星 = 40寸
using SymPy
A = Sym[5 3 2; 110 18 3; 855 65 5]
b = [10, 27.5, 40];
pqr = A\b
pqr |> println
Sym[-0.0942062043795621, 1.64119525547445, 2.77372262773723]
g(n) = -0.0942062043795621*n^2 + 1.64119525547445*n + 2.77372262773723
using SymPy
solve(diff(g))
8.71065375302662
g(7), g(8), g(9), g(10)
(9.645985401459837, 9.874087591240855, 9.91377737226275, 9.765054744525518)
9 番目が最大でそれ以降では減少する。そういう意味では，この問題は「無術」なのであろう。
角星の周の長さは 4.32071167883212 ということになる。
pqr' * [1;1;1]
4.32071167883212
しかし，この数値も特段の意味を持つものではないようだ。
ということで，改訂版（電子版）では著者の意図したように問題（の数値）が変更された。
