算額（その619）

算額（その619）

和算図形問題あれこれ
令和 4 年 4 月の問題 - No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html

正方形の一辺の長さを a とする。図のように，正方形内に 3 本の線を引く。

甲，乙，丙の三角形の面積が 111 平方寸，68 平方寸, 8 平方寸のとき，正方形の一辺の長さはいかほどか。

周りにある直角三角形の面積は計算できる。中央の三角形の面積を x とおき，以下の連立方程式を解けばよい。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, x::positive;

eq1 = a*b/2 - 68
eq2 = a*c/2 - 111
eq3 = (a - b)*(a - c)/2 - 8
eq4 = 111 + 8 + 68 + x - a^2
res = solve([eq1, eq2, eq3, eq4], (a, b, c, x))

   1-element Vector{NTuple{4, Sym}}:
    (sqrt(187 - sqrt(4777))*(sqrt(31191) + 11*sqrt(1887))/444, 2*sqrt(352869 - 1887*sqrt(4777))/111, sqrt(1221/4 - 111*sqrt(4777)/68), sqrt(4777))

a = sqrt(187 - sqrt(4777))*(sqrt(31191) + 11*sqrt(1887))/444 となり，SymPy ではこれ以上簡約化できない。
しかし，裏をかいて a^2 を求めて簡約化すると以下のようになる。この平方根をとれば a = sqrt(sqrt(4777) + 187) である。

res[1][1]^2 |> simplify |> println

   sqrt(4777) + 187

a = sqrt(sqrt(4777) + 187) = 16.003619739999756 である。

sqrt(sqrt(4777) + 187) 

   16.003619739999756

その他のパラメータは以下の通り。

   a = 16.0036;  b = 8.49808;  c = 13.8719;  x = 69.1158

function draw(more=false)
   pyplot(size=(300, 300), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, x) = res[1]
   @printf("正方形の一辺の長さ a = %g;  b = %g;  c = %g;  x = %g\n", a, b, c, x)
   plot([0, a, a, 0, 0,  a, c, 0], [0, 0, a, a, 0,  b, a, 0], color=:black, lw=0.5)
   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)
       point(a, 0, " a", :red, :left, :bottom, delta=delta/2)
       point(a, b, " (a,b)", :red, :left, :vcenter)
       point(c, a, "(c,a)", :red, :center, :bottom, delta=delta/2)
       point(0, a, "  a", :red, :center, :bottom, delta=delta/2)
       point(5, 12, "甲: 111", :green, :center, :vcenter, mark=false)
       point(12.5, 2.5, "乙: 68", :green, :center, :vcenter, mark=false)
       point(16, 13, "丙: 8 ", :green, :right, :vcenter, mark=false)
       point(11, 9, "x", :green, :center, :vcenter, mark=false)
   end
end;

算額（その618）

算額（その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;

算額（その617）

算額（その617）

三重県四日市市 神明神社 寛政2年(1790)
「三重県に現存する算額の研究」福島完(2007/2/13)
https://mie-u.repo.nii.ac.jp/?action=repository_uri&item_id=7216

問題文2
菱形の中に，大円 1 個，小円 2 個が入っている。菱形の面積から3個の円の面積を除いた面積(A)，菱形の対角線の長さの和(B)と，大円と小円の直径の差(C)が与えられているとき，小円の直径を求めよ。

菱形の対角線 2a, 2b (a > b)，大円と小円の半径をそれぞれ r1, r2 とする。

以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, y1::positive, r2::positive, y2::negative, a::positive, b::positive,
   A::positive, B::positive, C::positive;
円周率 = 3.16  # 円積率 3.16/4 = 0.79
eq1 = r2^2 + (y1 - y2)^2 - (r1 + r2)^2
eq2 = distance(0, b, a, 0, 0, y1) - r1^2
eq3 = distance(0, -b, a, 0, r2, y2) - r2^2
eq4 = 2a*b - 円周率*r1^2 - 2円周率*r2^2 - A
eq5 = 2a + 2b - B
eq6 = 2r1 - 2r2 - C;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, b, r1, y1, r2, y2) = u
   return [
       r2^2 - (r1 + r2)^2 + (y1 - y2)^2,  # eq1
       a^2*b^2*(b - y1)^2/(a^2 + b^2)^2 - r1^2 + (-b*(a^2 + b*y1)/(a^2 + b^2) + y1)^2,  # eq2
       -r2^2 + (-a*(a*r2 + b^2 + b*y2)/(a^2 + b^2) + r2)^2 + (-b*(-a^2 + a*r2 + b*y2)/(a^2 + b^2) + y2)^2,  # eq3
       2*a*b - 3.16*r1^2 - 6.32*r2^2 - 53.7,  # eq4
       2*a + 2*b - 32,  # eq5
       2*r1 - 2*r2 - 5.2,  # eq6
   ]
end;

(A, B, C) = (53.7, 32, 5.2)
iniv = BigFloat[8.62, 7.38, 4.23, 1.81, 1.63, -3.83]
res = nls(H, ini=iniv)
  (BigFloat[8.618869516619060933261217496129591164127618589189709440245477710397030837896864, 7.381130483380939066738782503870408835872381410810290559754522289595245106897612, 4.234238957828194708369843972428936199044440731565763308078691436517181518285074, 1.806377116536225577423953366326376812684128330422196332266830152186767547254349, 1.634238957828194619552002002416412965153907284300138308078691438199894855227677, -3.829959998582366768106243609407652459638817222378998651702439496268955544637784], true)

菱形の面積から3個の円の面積を除いた面積(A) = 53.7
菱形の対角線の長さの和(B) = 32
大円と小円の直径の差(C) = 5.2
のとき，
小円の直径 = 3.26848
その他のパラメータは以下の通りである。
a = 8.61887; b = 7.38113;  r1 = 4.23424;  y1 = 1.80638;  r2 = 1.63424;  y2 = -3.82996

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B, C) = (53.7, 32, 5.2)
   (a, b, r1, y1, r2, y2) = res[1]
   @printf("菱形の面積から3個の円の面積を除いた面積(A) = %g\n", A)
   @printf("菱形の対角線の長さの和(B) = %g\n", B)
   @printf("大円と小円の直径の差(C) = %g\n", C)
   @printf("小円の直径 = %g;  a = %g; b = %g;  r1 = %g;  y1 = %g;  r2 = %g;  y2 = %g\n", 2r2, a, b, r1, y1, r2, y2)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], color=:black, lw=0.5)
   circle(0, y1, r1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   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)
       point(0, y1, " 大円:r1,(0,y1)", :red, :left, :vcenter)
       point(r2, y2, "小円:r2\n(r2,y2)", :blue, :center, delta=-delta/2)
       point(0, b, " b", :black, :left, :bottom)
       point(a, 0, " a", :black, :left, :bottom)
   end
end;