算額（その618）
三重県四日市市 神明神社 寛政2年(1790)
「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216
今，「角」から「軫」まで二十八宿の名前を冠する 28 種の香がある。図のように環状に置き，初日に「角」を薫し，二日目に「氐」，三日目に「尾」，四日目に「女」のように n 日目に n*(n + 1)/2 番目の香を薫す。なお 28 番目の「軫」を過ぎれば，29 番目は「角」になる。7 番目の「箕」が焚かれるのは何日目か。
n 日目に焚かれる香の順番を an とすれば，a1 = 1, a2 = 2, a3 = 6, ..., an = n*(n + 1)/2 である。
「箕」の番号は 7, 35, 63, ..., 7 + 28k である。
両者が一致すると，その日に「箕」が焚かれるので，n*(n + 1)//2 = 7 + 28k を満たす n を求めればよい。
include("julia-source.txt");
using SymPy
@syms n, an, k, 箕
an = n*(n + 1)//2
箕 = 7 + 28k
res = solve(an - 箕, n)
res |> println
Sym[-sqrt(224*k + 57)/2 - 1/2, sqrt(224*k + 57)/2 - 1/2]
n は正の整数なので，少なくとも 2 番目のものが適解である。
さらに，224*k + 57 は平方数でなくてはならない。
k を順次増やして計算すればよい。k = 3， n = sqrt(224*k + 57)/2 - 1/2 = 13 は手計算でも見つかるが，そのあとはちょっと面倒かもしれないので，プログラムを書いてみる。
sqrt(224*3 + 57)/2 - 1/2
13.0
for k = 1:100
k2 = 224*k + 57
if isapprox(isqrt(k2), sqrt(k2))
println("k = $k, k2 = $k2, $(Int(sqrt(224*k + 57)/2 - 1/2))日目")
end
end
k = 3, k2 = 729, 13日目
k = 8, k2 = 1849, 21日目
k = 21, k2 = 4761, 34日目
k = 32, k2 = 7225, 42日目
k = 86, k2 = 19321, 69日目
当然といえば当然だが，簡単な数列ではない。
オンライン整数列大辞典 https://oeis.org/welcome にもない。
ブルートフォースで解を求めるのも簡単でよい。
二十八宿 = [
"角", "亢", "氐", "房", "心", "尾", "箕", # 東方青龍
"斗", "牛", "女", "虚", "危", "室", "壁", # 北方玄武
"奎", "婁", "胃", "昴", "畢", "觜", "参", # 西方白虎
"井", "鬼", "柳", "星", "張", "翼", "軫" # 南方朱雀
]
for n = 1:100
an = Int(n*(n + 1)//2)
suffix = mod(an - 1, 28) + 1
if suffix == 7
println(n, "日目，通し番号 ", an, " 番目の香，通し番号を 28 で割ったあまり ", suffix, "。香の名前は ", 二十八宿[suffix])
end
end
13日目，通し番号 91 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
21日目，通し番号 231 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
34日目，通し番号 595 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
42日目，通し番号 903 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
69日目，通し番号 2415 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
77日目，通し番号 3003 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
90日目，通し番号 4095 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
98日目，通し番号 4851 番目の香，通し番号を 28 で割ったあまり 7。香の名前は 箕
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r) = (1, 0.1)
plot(showaxis=false)
rotate(0, R, 1.1r, :blue, angle=360//28)
for i = 1:28
θ = (i + 6)*360/28
x = (R - 0.01r)*cosd(θ)
y = (R - 0.01r)*sind(θ)
annotate!(x, y, text(二十八宿[i], :center, :vcenter, 18))
end
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
#vline!([0], color=:black, lw=0.5)
#hline!([0], color=:black, lw=0.5)
end
end;
