算額（その432）

算額（その432）

神奈川県横須賀市 浦賀叶大明神 文化11年(1814)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円に等円 2 個が内接し，その上下に甲円と乙円が大円に内接し，等円に外接している。大円の直径が 1 尺，甲円の直径が 5 寸のとき乙円の直径はいかほどか。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (0, r0 - r1)
等円の半径と中心座標を r2, (r2, y2)
乙円の半径と中心座標を r3, (0, r3 - r0)
とし，以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive,
     y2::negative, r3::positive;
eq1 = r2^2 + (r0 - r1 - y2)^2 - (r1 + r2)^2
eq2 = r2^2 + (y2 - r3 + r0)^2 - (r2 + r3)^2
eq3 = r2^2 + y2^2 - (r0 - r2)^2
res = solve([eq1, eq2, eq3], (r2, r3, y2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*r0*r1*(r0 - r1)/(r0 + r1)^2, r0*(r0 - r1)/(r0 + 3*r1), (r0^2 - 3*r0*r1)/(r0 + r1))

術および中村の現代解法は乙円の半径についてかなり複雑な式を提示している。
SymPy での結果は，乙円の半径は，r0*(r0 - r1)/(r0 + 3*r1) すなわち，「大円の半径と大円の半径と甲円の半径の差をかけたものを大円の半径と甲円の半径の3倍の和で割る」となり，藤井貞夫が「長野県和算研究会報No.5」で指摘したものになった。

大円の直径が 1 尺，甲円の直径が 5 寸のとき乙円の直径は 2 寸である。

2res[1][2](r0 => 10/2, r1 => 5/2).evalf() |> float |> println

   2.0

   r2 = 2.22222;  r3 = 1;  y2 = -1.66667;  乙円の直径 = 2

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r0 = 10//2
   r1 = 5//2
   (r2, r3, y2) = (4*r0*r1*(r0 - r1)/(r0 + r1)^2, r0*(r0 - r1)/(r0 + 3*r1), (r0^2 - 3*r0*r1)/(r0 + r1))
   @printf("r2 = %g;  r3 = %g;  y2 = %g;  乙円の直径 = %g\n", r2, r3, y2, 2r3)
   plot()
   circle(0, 0, r0, :blue)
   circle(0, r0 - r1, r1, :red)
   circle(r2, y2, r2, :green)
   circle(-r2, y2, r2, :green)
   circle(0, r3 - r0, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r0 - r1, " 甲円:r1,(0,r0-r1)", :red, :left, :vcenter)
       point(r2, y2, "(等円:r2,(r2,y2)", :green, :center, delta=-delta)
       point(0, r3 - r0, " 乙円:r3\n (0,r3-r0)", :black, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その431）

算額（その431）

東京都浅草 池田邸内 伽山

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円 2 個，中円 5 個，小円 4 個が図のように互いに内接・外接している。大円の直径が 21 寸のとき，小円の直径はいかほどか。

大円の半径と中心座標を r1, (0, r2), (0, -r2)
中円の半径と中心座標を r2, (2r2, 0)
小円の半径と中心座標を r3, (x3, x4)
とし，以下の連立方程式の解を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     x3::positive;
r2 = r1//2
eq1 = (2r2 - x3)^2 + x3^2 - (r2 + r3)^2
eq2 = x3^2 + (x3 + r2)^2 - (r1 + r3)^2
res = solve([eq1, eq2], (r3, x3))

   1-element Vector{Tuple{Sym, Sym}}:
    (3*r1/14, 4*r1/7)

大円の直径が 21 寸のとき，小円の直径は 21 * 3 / 14 = 4.5 寸である。

21 * 3 / 14

   4.5

r3 = 2.25;  x3 = 6;  小円の直径 = 4.5

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 21//2
   r2 = r1//2
   (r3, x3) = (3*r1/14, 4*r1/7)
   @printf("r3 = %g;  x3 = %g;  小円の直径 = %g\n", r3, x3, 2r3)
   plot()
   circle(0, r2, r1, :red)
   circle(0, -r2, r1, :red)
   circle(0, 0, r2, :green)
   circle42(0, 2r2, r2, :green)
   circle4(x3, x3, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r2, " r2\n 大円:r0", :red, :left, :top, delta=-delta)
       point(0, 2r2, " 2r2\n 中円:r2", :green, :left, :vcenter)
       point(x3, x3, " 小円:r3,(x3,x3)", :black, :left, :vcenter)
       point(2r2, 0, "2r2", :green, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

#算額

算額（その430）

算額（その430）

東京都浅草 池田邸内 伽山

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

大円の中に甲円 2 個と乙円 1 個，等円 3 個が入っている。等円は 1 本の元と長さの等しい 3 本の面（隣り合う 3 本の弦）に接している。面の長さが 1 寸のとき，乙円の直径を求めよ。

大円の半径と中心座標を r0, (0, 0)
甲円の半径と中心座標を r1, (r1, y1) ただし y1 = a + r1
乙円の半径と中心座標を r2, (0, r0 - r2)
等円の半径と中心座標を r3, (0, a - r3), (2r3, a - r3)
弦と y 軸の交点座標を a
真ん中の面と y 軸の交点座標を b
面の長さの半分を c とし，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r0::positive, r1::positive, r2::positive,
     r3::positive, a::negative;
c = 1//2
b = a - 2r3
d = sqrt(r0^2 - a^2)
y1 = a + r1
eq1 = r1^2 + (r0 - r2 - y1)^2 - (r1 + r2)^2
eq2 = r1^2 + y1^2 - (r0 - r1)^2
eq3 = r0^2 - b^2 - c^2
eq4 = distance(c, b, d, a, 2r3, b+ r3) - r3^2
eq5 = (d - c)^2 + (a - b)^2 - (2c)^2;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, r2, r3, a))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r0, r1, r2, r3, a) = u
   return [
       r1^2 - (r1 + r2)^2 + (-a + r0 - r1 - r2)^2,  # eq1
       r1^2 + (a + r1)^2 - (r0 - r1)^2,  # eq2
       r0^2 - (a - 2*r3)^2 - 1/4,  # eq3
       -r3^2 + (2*r3 - 2*r3*(-4*a^2 + 4*r0^2 + 4*r3*sqrt(-a^2 + r0^2) + 2*r3 - 4*sqrt(-a^2 + r0^2) + 1)/(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))^2 + (a - r3 - (2*r3*(8*r3^2 + 8*r3*sqrt(-a^2 + r0^2) - 4*r3 - 2*sqrt(-a^2 + r0^2) + 1) + (a - 2*r3)*(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))/(-4*a^2 + 4*r0^2 + 16*r3^2 - 4*sqrt(-a^2 + r0^2) + 1))^2,  # eq4
       4*r3^2 + (sqrt(-a^2 + r0^2) - 1/2)^2 - 1,  # eq5
   ]
end;

iniv = [big"65.0", 31, 27, 10, -32] ./ 40
res = nls(H, ini=iniv);
println(res);
   (BigFloat[2.176250899482821511100052865997767880198073191589329947230101745924818974139391, 0.9777711616568012201157129796477534133807314665018133584786185122552240330104753, 1.070019583175477769385089300350675253824657209000148549887247280185448490007896, 0.2236067977499789696409173668731276235440618359611525724270897245410529337440141, -1.670820393249936908922752100619382870632185507883457717281269173623140154153112], true)

r0 = 2.17625;  r1 = 0.977771; , r2 = 1.07002; , r3 = 0.223607; , a = -1.67082
乙円の直径 = 2.14004;  甲円の直径 = 1.95554;  等円の直径 = 0.447214
面 = 1;  弦 = 2d = 2.78885;  矢 = r0 + a = 0.505431

面の長さが 1 寸のとき，乙円の直径は 2.14004 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 169.3/2
   (r0, r1, r2, r3, a) = res[1]
   c = 1//2
   b = a - 2r3
   d = sqrt(r0^2 - a^2)
   y1 = a + r1
   @printf("r0 = %g;  r1 = %g; , r2 = %g; , r3 = %g; , a = %g\n", r0, r1, r2, r3, a)
   @printf("乙円の直径 = %g;  甲円の直径 = %g;  等円の直径 = %g\n", 2r2, 2r1, 2r3)
   @printf("面 = %g;  弦 = 2d = %g;  矢 = r0 + a = %g\n",
       sqrt((sqrt(r0^2 - a^2) - c)^2 + (a - b)^2), 2d, r0 + a)
   plot()
   circle(0, 0, r0, :red)
   circle(r1, y1, r1, :green)
   circle(-r1, y1, r1, :green)
   circle(0, r0 - r2, r2, :magenta)
   circle(0, a - r3, r3, :orange)
   circle(2r3, a - r3, r3, :orange)
   circle(-2r3, a - r3, r3, :orange)
   segment(-c, b, c, b, :blue)
   segment(-d, a, d, a, :blue)
   segment(c, b, d, a, :blue)
   segment(-c, b, -d, a, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, y1, "甲円:r1,(r1,y1)", :green, :center, :top, delta=-delta)
       point(0, r0 - r2, " 乙円:r2,(0,r0-r2)", :magenta, :center, :top, delta=-delta)
       point(0, a - r3, "a-r3 ", :black, :right, :vcenter)
       point(2r3, a - r3, "(2r3,a-r3)", :black, :center, :bottom, delta=delta)
       point(c, b, "(c,b)")
       point(d, a, "(d,a)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

術では，乙円径は 1.4221...寸としているが，この値は間違っている。
術の式通りに計算すると下記のように 2.14004...寸となる。

@syms 角, 元, 氏, 房, 面

面 = 1
角 = sqrt(80)
元 = sqrt(角 + 10)
氏 = 元 - sqrt(角 + 5.8 - sqrt(12.8))
房 = sqrt(2氏 * 元) - 氏
乙円径 = (元 * 房 * 面) / (2元 -(2sqrt((元 -2房)元) + 房))
@printf("面 = %g;  角 = %g;  元 = %g;  氏 = %g;  房 = %g\n", 面, 角, 元, 氏, 房)  
@printf("乙円径 = %g\n", 乙円径)  

   面 = 1;  角 = 8.94427;  元 = 4.3525;  氏 = 1.01086;  房 = 1.95554
   乙円径 = 2.14004

中村は，術と中村の導いた乙円径を表す式が同じであることを示しているが，(7)の分母の式において（またしても）符号を誤っているのでそれに基づいたのでは正しい答えにならない。"+ 房" は "- 房" である。

@printf("大円径 = %g\n", 元*面)

@printf("h(矢) = %g\n", 氏*面/2)

@printf("(7) の分子 = 甲円径 = %g\n", 房*面)
@printf("(7) の分母 = %g\n", (2元 - 2sqrt(元*(元 - 2房)) - 房)面)

x = (R*d/((2元 - 2sqrt(元*(元 - 2房)) - 房)面))
@printf("乙円径 = %g\n", x.evalf())

   大円径 = 4.3525
   h(矢) = 0.505431
   (7) の分子 = 甲円径 = 1.95554
   (7) の分母 = 3.97726
   乙円径 = 2.14004

以下は，中村による現代解法である。

@syms R,  # 大円径
     r,  # 等円径
     d,  # 甲円径
     x,  # 乙円径
     h,  # 矢
     a   # 面

a = 1

r = a/sqrt(Sym(5))
@printf("等円径 = %g\n", r.evalf())

R = sqrt(10 + 4sqrt(Sym(5)))a
@printf("大円径 = %g\n", R.evalf())

h = (sqrt(10 + 4sqrt(Sym(5))) - sqrt(4sqrt(Sym(5)) + 29//5 - 8sqrt(Sym(5))/5))a/2
@printf("h(矢) = %g\n", h.evalf())

eq3 = (d/2)^2 - h*(R - d) + h^2;
d = solve(eq3, d)[1]
@printf("甲円径 = %g\n", d.evalf())

x = (R*(2R - d) + 2R*sqrt(R*(R - 2d)))/(d + 4R)
@printf("乙円径 = %g\n", x.evalf())

   等円径 = 0.447214
   大円径 = 4.3525
   h(矢) = 0.505431
   甲円径 = 1.95554
   乙円径 = 2.14004

 

 

#算額

算額（その429）

算額（その429）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

2 個の大円が交わっており，正方形 3 個，小円 2 個が入っている。大円の直径が 169.3 寸のとき，小円の直径を求めよ。

正方形の一辺の長さを 2a
大円の半径と中心座標を r1, (x1, 0), (-x1, 0)
小円の半径と中心座標を r2, (0, a + r2)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     a::positive;

eq1 = x1^2 + (a + r2)^2 - (r1 - r2)^2
eq2 = a^2 + (x1 + a)^2 - r1^2
eq3 = (r1 - 2x1 + 2a)^2 + a^2 - r1^2
res = solve([eq1, eq2, eq3], (r2, x1, a))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (2*r1*(4 - sqrt(6))/15, r1*(2*sqrt(6) + 7)/25, r1*(-4/25 + 6*sqrt(6)/25))

小円の半径は大円の半径の 2(4 - √6)/15 倍である。
小円の直径が 169.3 のとき，大円の直径は 35.00018487290774 である。

res[1][1](r1 => 169.3/2).evalf() |> println
res[1][2](r1 => 169.3/2).evalf() |> println
res[1][3](r1 => 169.3/2).evalf() |> println

   17.5000924364539
   40.2899445381277
   36.2198336143830

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 169.3/2
   (r2, x1, a) = (2*r1*(4 - sqrt(6))/15, r1*(2*sqrt(6) + 7)/25, r1*(-4/25 + 6*sqrt(6)/25))
   @printf("r2 = %g;  x1 = %g;  a = %g\n", r2, x1, a)
   @printf("小円の直径 = %.7g\n", 2r2)
   plot()
   circle(x1, 0, r1, :red)
   circle(-x1, 0, r1, :red)
   rect(-a, -a, a, a, :blue)
   circle(0, a + r2, r2, :green)
   rect(r1 - x1, -a, r1 - x1 + 2a, a, :blue)
   rect(-r1 + x1, -a, -r1 + x1 - 2a, a, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1 - x1, 0, " r1-x1", :red, :left, delta=-delta)
       point(x1 - r1, 0, "x1-r1", :red, :right, delta=-delta)
       point(x1, 0, "x1 ", :red, :right, :bottom, delta=delta/2)
       point(-a, a, "(-a,a)", :blue)
       point(r1 - x1 + 2a, a, "(r1-x1+2a,a)", :blue, :right, :top, delta=-delta)
       point(0, a + r2, " a+r2", :green)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

#算額

算額（その428）

算額（その428）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

外円の中に長方形が内接している。一辺の長さが 41.86 寸の正方形が 6個，大円が 2 個，小円が 4 個入っている。小円の直径を求めよ。

正方形の一辺の長さを a
大円の半径と中心座標を r1, (r0 - r1)
小円の半径と中心座標を r2, (r0 - 2r1 + r2, y)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, r0::positive, r1::positive,
     r2::positive, y::positive;
@syms a, r0, r1, r2, y
eq1 = (sqrt(Sym(2))a/2 + a)^2 + (r0 - sqrt(Sym(2))a + a)^2 - r0^2
eq2 = (r0 - 2r1)^2 + (r0 - sqrt(Sym(2))a)^2 - r0^2
eq3 = (r0 - 2r1 + r2)^2 + y^2 - (r0 - r2)^2
eq4 = (r1 - r2)^2 + y^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (r0, r1, r2, y));

2 組の解が得られるが，2 番目のものが適解である。

names = ["r0", "r1", "r2", "y"]
for (i, name) in enumerate(names)
   println("\n", name)
   println(simplify(res[2][i]))
   println(res[1][i](a => 4186).evalf())
end

   r0
   a*(5 + 7*sqrt(2))/4
   15592.3214511641

   r1
   (-23922439916777731675282139353831061069516954103552970349841905*a^2 + 16915719487681281681493317896073504500369711065425523483418778*sqrt(2)*a^2 - 6986054008647705185868918429789392*sqrt(5)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4) + 4939886163250256057660601952272630*sqrt(10)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4))/(8*a*(4882633868649679478319412976906635182914530485393454527631129 - 3452543518573294933348514588454588356006400466465731641999825*sqrt(2)))
   1681.32809749726

   r2
   a*(205 + 151*sqrt(2))/1168
   1500.02961796760

   y
   -sqrt(-6848627247131292427665721780*sqrt(2)*a^2 + 9685421536530988262785792422*a^2 - 4*sqrt(10)*sqrt(1144573653449326586061314076037948929840593626875394274 - 809335791921480250432937095033206217367486276799679753*sqrt(2))*sqrt(a^4))/(8*sqrt(8442762866313367764895228377 - 5969934874720155329318365400*sqrt(2)))
   3176.18761647798

小円の半径は a*(205 + 151*√2)/1168 なので，直径は a*(205 + 151*√2)/584 である。

a = 41.86 を代入すると，小円の直径は 30 寸あまりあり。

41.86*(205 + 151*√2)/584

   30.00059235935206

res2 = Float64.([res[2][i](a => 4186).evalf()/100 for i in 1:4])

   4-element Vector{Float64}:
    155.92321451164108
     16.81328097497264
     15.00029617967603
     31.761876164779796

r0 = 155.923;  r1 = 16.8133;  r2 = 15.0003;  y = 31.7619
小円の直径 = 30.0006

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r2, y) = res2
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  y = %g\n", r0, r1, r2, y)
   @printf("小円の直径 = %g\n", 2r2)
   a = 41.86
   (x0 , y0) = (r0 - 2r1, r0 - √2a)
   plot([x0, x0, -x0, -x0, x0], [-y0, y0, y0, -y0, -y0], color=:orange, lw=0.5)
   plot!([0, √2a/2, 0, -√2a/2, 0], [r0 - √2a, r0 - √2a + √2a/2, r0, r0 - √2a + √2a/2, r0 - √2a], color=:red, lw=0.5)
   rect(√2a/2, r0 - √2a, √2a/2 + a, r0 - √2a + a, :red)
   rect(-√2a/2, r0 - √2a, -√2a/2 - a, r0 - √2a + a, :red)
   plot!([0, √2a/2, 0, -√2a/2, 0], -[r0 - √2a, r0 - √2a + √2a/2, r0, r0 - √2a + √2a/2, r0 - √2a], color=:red, lw=0.5)
   rect(√2a/2, -r0 + √2a, √2a/2 + a, -r0 + √2a - a, :red)
   rect(-√2a/2, -r0 + √2a, -√2a/2 - a, -r0 + √2a - a, :red)
   circle(0, 0, r0, :black)
   circle(r0 - r1, 0, r1, :blue)
   circle(-r0 + r1, 0, r1, :blue)
   circle4(r0 - 2r1 + r2, y, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta)
       point(r0 - r1, 0, "r0-r1", :blue, :center, delta=-delta/2)
       point(r0 - 2r1 + r2, y, "(r0-2r1+r2,y) ", :green, :right, :vcenter)
       point(0, r0 - √2a, "r0-√2a ", :orange, :right, :top, delta=-delta/2)
       point(a/√2, r0-√2a, "(a/√2,r0-√2a)", :red, :center, delta=-3delta)
       point(a/√2 + a, r0-√2a, "(a/√2+a,r0-√2a)", :red, :center, delta=-delta/2)
       point(a/√2 + a, r0 - √2a + a, "(a/√2+a,r0-√2a+a)", :red, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その427）

算額（その427）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

大円は鈎股弦（直角三角形）の，鈎と弦に接しており，中円，小円は弦と大円に接しかつそれぞれも互いに接している。鈎，股がそれぞれ　3000 寸，4000 寸，小円の直径が 456 寸のとき，大円の直径を求めよ。

鈎，股，弦の長さをそれぞれ鈎，股，弦
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, x2)
小円の半径と中心座標を r3, (x31, y31), (x32, y32)
として，連立方程式を解く。

solve() の能力的に一度では解けないので，まず最初に大円と中円の半径を求める。このためには，図を反時計回りに回転させ，弦（下図の紫の線）が水平になるように回転させて方程式を解くとよい。eq13 はr1 と鈎股弦の面積の関係式である（上図を参照）。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, r2::positive,
     r3::positive, x3::positive;

eq11 = x3^2 + (r1 - 2r2 + r3)^2 - (r1 - r3)^2
eq12 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2
eq13 = r1*(鈎 + 股) + (r1 - 2r2)*弦 - 鈎*股
res1 = solve([eq11, eq12, eq13], (r1, r2, x3))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    ((-2*r3*弦^2 - 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 - sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), -sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) + 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 - 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 - sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) + 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 + 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 + sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), -sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) - 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))
    ((-2*r3*弦^2 + 弦*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2) + 股^2*鈎 + 股*鈎^2)/(-弦^2 + 股^2 + 2*股*鈎 + 鈎^2), (-r3*弦 + 股*鈎/2 + sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/2)/(-弦 + 股 + 鈎), sqrt(4*r3^2*弦/(弦 - 股 - 鈎) - 2*r3*股*鈎/(弦 - 股 - 鈎) - 2*r3*sqrt(4*r3^2*弦^2 - 4*r3*股^2*鈎 - 4*r3*股*鈎^2 + 股^2*鈎^2)/(弦 - 股 - 鈎)))

4 組の解が得られるが，そのうちの最初の 2 組が適解である。1 番目と 2 番目は，小円の位置（x3）が対称なもので実質的に同じである。

println("r1 = ", res1[2][1](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))
println("r2 = ", res1[2][2](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))
println("x3 = ", res1[2][3](鈎 => 3000, 股 => 4000, 弦 => 5000, r3 => 456//2))

   r1 = 1250
   r2 = 300
   x3 = 120*sqrt(19)

算額の解としては大円の半径が 1250，直径は 2500 でここまででよい。

次いで，図を描くために，r1, r2 が既知として，小円の中心座標を求める。

using SymPy
@syms r1::positive, r2::positive, x2::positive, y2::positive,
     r3::positive, x31::positive, y31::positive,
     x32::positive, y32::positive;

(r1, r2, r3) = (1250, 300, 456//2)
x2 = r1 + (r1 - r2)*3//5
y2 = r1 + (r1 - r2)*4//5
eq1 = (x2 - x31)^2 + (y2 - y31)^2 - (r2 + r3)^2
eq2 = (x2 - x32)^2 + (y2 - y32)^2 - (r2 + r3)^2
eq5 = (x31 - r1)^2 + (y31 - r1)^2 - (r1 - r3)^2
eq6 = (x32 - r1)^2 + (y32 - r1)^2 - (r1 - r3)^2
res2 = solve([eq1, eq2, eq5, eq6], (x31, y31, x32, y32))

   4-element Vector{NTuple{4, Sym}}:
    (8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5)
    (8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19))
    (96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5)
    (96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19))

4 組の解が得られるが，2 番目と 3 番目のものが適解である（小円の位置関係が逆になるだけ）。

(8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5, 96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19)) |> println
(96*sqrt(19) + 8884/5, 9762/5 - 72*sqrt(19), 8884/5 - 96*sqrt(19), 72*sqrt(19) + 9762/5) |> println

   (1358.3457014200953, 2266.2407239349286, 2195.2542985799046, 1638.5592760650716)
   (2195.2542985799046, 1638.5592760650716, 1358.3457014200953, 2266.2407239349286)

2 つの解をまとめると，以下のようになる。

   鈎 = 3000;  股 = 4000;  r3 = 228
   r1 = 1250;  r2 = 300;  x31 = 2195.25;  y31 = 1638.56;  x32 = 1358.35;  y32 = 2266.24
   大円の直径 = 2500

using Plots

function draw0(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x3) = (1250, 300, 456//2, 120*sqrt(19))
   plot()
   circle(0, 0, r1)
   circle(0, r1 - r2, r2, :blue)
   circle(x3, r1 - 2r2 + r3, r3, :green)
   circle(-x3, r1 - 2r2 + r3, r3, :green)
   x = sqrt(r1^2 - (r1 - 2r2)^2)
   segment(-x, r1 - 2r2, x, r1 - 2r2, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1-r2, " r1-r2", :blue)
       point(0, r1-2r2, " r1-2r2", :magenta, delta=-delta)
       point(x3, r1 - 2r2 + r3, "(x3,r1-2r2+r3)", :green, :center, :top, delta=-delta)
       point(r1, 0, "r1 ", :red, :right, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, r3) = (3000, 4000, 456//2)
   (r1, r2) = (1250, 300)
   (x31, y31, x32, y32) = res2[3]
   x2 = r1 + (r1 - r2)*3/5
   y2 = r1 + (r1 - r2)*4/5
   @printf("鈎 = %g;  股 = %g;  r3 = %g\n", 鈎, 股, r3)
   @printf("r1 = %g;  r2 = %g;  x31 = %g;  y31 = %g;  x32 = %g;  y32 = %g\n", r1, r2, x31, y31, x32, y32)
   @printf("大円の直径 = %.8g\n", 2r1)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, y2, r2, :blue)
   circle(x31, y31, r3, :green)
   circle(x32, y32, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, r1, " 大円:r1,(r1,r1)", :red, :left, :vcenter)
       point(r1 + (r1 - r2)*3/5, r1 + (r1 - r2)*4/5, " 中円:r2,(x2,y2)", :blue, :left, :vcenter)
       point(x31, y31, " 小円:r3,(x31,y31)", :green, :left, :vcenter)
       point(x32, y32, " 小円:r3,(x32,y32)", :green, :left, :bottom, delta=delta)
       point(股, 0, "股", :black, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その426）

算額（その426）

長野県東筑摩郡筑北村青柳 碩水寺豊川稲荷社 慶応4年(1868)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

等脚台形内に中円が内接している。大円，小円の直径は等脚台形の下底と上底に等しい。それぞれの円の中心は等脚台形の右下角，右上角にあり，互いに接している。中円，小円の直径がそれぞれ 6 寸，4 寸のとき，大円の直径を求めよ。

大円の半径，中心座標を r1, (r1, 0)
中円の半径，中心座標を r2, (0, r2)
小円の半径，中心座標を r3, (r3, 2r2)
とおき，以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive;

eq1 = (r1 - r3)^2 + (2r2)^2 - (r1 + r3)^2
res = solve([eq1], (r1))

   Dict{Any, Any} with 1 entry:
     r1 => r2^2/r3

大円の半径は，中円の半径の2乗を小円の半径で割って得られる。

中円，小円の半径，が 3 寸，2 寸のとき，大円の半径は 9/2 寸であり，直径は 9 寸である。

res[r1](r2 => 3, r3 => 2) |> println

   9/2

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (6, 4) ./ 2
   r1 = r2^2/r3
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   @printf("大円の直径 = %g\n", 2r1)
   plot([r1, r3, -r3, -r1, r1], [0, 2r2, 2r2, 0, 0], color=:black, lw=0.5)
   circle(r1, 0, r1)
   circle(0, r2, r2, :blue)
   circle(r3, 2r2, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1, 0, "大円:r1,(r1,0)", :red, :center, :top, delta=-delta)
       point(0, r2, "中円:r2,(0,r2) ", :blue, :right, :vcenter)
       point(r3, 2r2, "小円:r3,(r3,2r2)", :green, :center, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その425）

算額（その425）

長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

長方形内に，甲，乙，丙，丁，戊，己の 6 個の円がある。それぞれ長方形に内接したり，互いに外接しあっている。
己円の直径を丁円と戊円の直径で表わせ。

長方形の長辺の長さを a とする（短辺の長さは 甲円の直径に等しい）。
甲円の半径と中心座標 r1, (a - r1, r1)
乙円の半径と中心座標 r2, (r2, r2)
丙円の半径と中心座標 r3, (x3, 2r1 - r3)
丁円の半径と中心座標 r4, (a - r4, r4)
戊円の半径と中心座標 r5, (r5, r5)
己円の半径と中心座標 r6, (x6, y6)

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive,
     x3::positive, r4::positive, r5::positive, r6::positive,
     x6::positive, y6::positive;

eq1 = (a - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (2r1 - r3 - r1)^2 - (r1 + r3)^2
eq3 = 2(r1 - r4)^2 - (r1 + r4)^2
eq4 = (a - r1 - x6)^2 + (r1 - y6)^2 - (r1 + r6)^2
eq5 = (x3 - r2)^2 + (2r1 - r3 - r2)^2 - (r2 + r3)^2
eq6 = 2(r2 - r5)^2 - (r2 + r5)^2
eq7 = (r2 - x6)^2 + (r2 - y6)^2 - (r2 + r6)^2
eq8 = (x3 - x6)^2 + (2r1 - r3 - y6)^2 - (r3 + r6)^2
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (a, r1, r2, r3, x3, r6, x6, y6))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, r1, r2, r3, x3, r6, x6, y6) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (a - r1 - r2)^2,  # eq1
       (r1 - r3)^2 - (r1 + r3)^2 + (a - r1 - r3)^2,  # eq2
       2*(r1 - r4)^2 - (r1 + r4)^2,  # eq3
       -(r1 + r6)^2 + (r1 - y6)^2 + (a - r1 - x6)^2,  # eq4
       (-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2,  # eq5
       2*(r2 - r5)^2 - (r2 + r5)^2,  # eq6
       -(r2 + r6)^2 + (r2 - x6)^2 + (r2 - y6)^2,  # eq7
       -(r3 + r6)^2 + (x3 - x6)^2 + (2*r1 - r3 - y6)^2,  # eq8
   ]
end;

(r4, r5) = (2, 1)
iniv = [big"180.0", 59, 35, 21, 55, 5, 65, 62] ./2
res = nls(H, ini=iniv);
println(res);

  (BigFloat[33.97056274847714058562026469051637694283606250452337687812015685588878974154529, 11.65685424949238019520675489683879231427868750150779229270671895196292991384848, 5.828427124746190097603377448419396157139343750753896146353359475981464956924239, 5.828427124746190097603377448419396157139343750753896146353359475981464956924308, 5.82842712474619009760948182491957344489198207867602126024818280529793350621324, 1.093836321356054313599895375610329707159527028223527974012807033844387485155129, 9.563017928136325881606859521228462607119160473284264318693911918118542430249862, 11.65685424949238019520479919353753625575412910164858126097863477985370829598065], true)

r4 = 2;  r5 = 1
a = 33.9706;  r1 = 11.6569;  r2 = 5.82843;  r3 = 5.82843;  x3 = 5.82843;  r6 = 1.09384;  x6 = 9.56302;  y6 = 11.6569

己円の直径は 2.18767 であるが，術では以下のように 7.31370849898476 であり，答えが違う。中村が術を読み違えているのか，式に起こすときに誤記があるのか。

(丁, 戊) = (2, 1) .* 2
天 = 丁 * 戊
地 = 丁 + 戊
法 = 2(sqrt((4地 * 戊 + 丁^2)天) + 2戊^2) + 地*丁
法*丁/(3 + √8)天

   7.31370849898476

中村は下記の x1 式を簡約化して x2 を提示し術と同じとしている。
しかし，x1 は正しいが x2 は式の記述ミスがある。正しくは x3（分母の式の一箇所の - が + でなければならない）。

@syms d, e
x1 = (3 + 2sqrt(Sym(2)))*d^2*e*((d^2 + d*e + 4e^2) - 2(d + 2e)*sqrt(d*e))/(d^4 - 2d^3*e - 7d^2*e^2 - 8d*e^3 + 16e^4)
x2 = (3 + 2sqrt(Sym(2)))*d^2*e / (d^2 + d*e + 4e^2 - 2(d + 2e)*sqrt(d*e))
x3 = (3 + 2sqrt(Sym(2)))*d^2*e / (d^2 + d*e + 4e^2 + 2(d + 2e)*sqrt(d*e))
x1(d => 4, e => 2).evalf() |> println
x2(d => 4, e => 2).evalf() |> println
x3(d => 4, e => 2).evalf() |> println

   2.18767264271211
   -35.4929704984046
   2.18767264271211

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r4, r5) = (2, 1)
   (a, r1, r2, r3, x3, r6, x6, y6) = res[1]
   @printf("r4 = %g;  r5 = %g\n", r4, r5)
   @printf("a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  r6 = %g;  x6 = %g;  y6 = %g\n", a, r1, r2, r3, x3, r6, x6, y6)
   @printf("己円の直径 = %g\n", 2r6)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(x3, 2r1 - r3, r3, :green)
   circle(a - r4, r4, r4, :magenta)
   circle(r5, r5, r5, :orange)
   circle(x6, y6, r6, :purple)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a - r1, r1, " 甲円:r1,(a-r1,r1)", :red, :left, :vcenter)
       point(r2, r2, " 乙円:r2,(r2,r2)", :blue, :center, :top, delta=-delta)
       point(x3, 2r1 - r3, " 丙円:r3,(x3,2r1-r3)", :green, :center, :top, delta=-delta)
       point(a - r4, r4, "丁円:r4,(a-r4,r4) ", :magenta, :right, :vcenter)
       point(r5, r5, "    戊円:r5,(r5,r5) ", :orange, :left, :vcenter)
       point(x6, y6, "  己円:r6,(x5,y5) ", :purple, :left, :vcenter)
       point(a, 2r1, "(a,2r1) ", :black, :right, :top, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その424）

算額（その424）

長野県諏訪市中洲 諏訪大社上社 明治12年(1879)

中村信弥(1999)：算額への招待
http://www.wasan.jp/syotai/syotai.html

直線上に大円 1 個，小円 3 個が載っている。大円と小円 2 個は直線に接し，大円と小円 2 個は接している。小円 3 個は互いに接している。このとき，大円の直径を 3 個の小円で囲まれた部分の面積（黒積）で表わせ。

大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (x2, r2), (x2 - r2, (1 + √3)r2)
黒積は √3r2^2 - πr2^2/2 = (√3 - π/2)r2^2 である。

以下の連立方程式を立て，数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::negative,
     黒積::positive;

まず，小円の半径 r2 を黒積で表す。

eq3 = (sqrt(Sym(3)) - PI/2)r2^2 - 黒積
r2 = solve(eq3, r2)[1]
println("r2: ", r2)

   r2: sqrt(2)*sqrt(黒積)/sqrt(-pi + 2*sqrt(3))

以下の二元連立方程式で大円の半径 r1 と小円の中心の x 座標 x2 を求める。

eq1 = (x2 - r2)^2 + (r1 - (1 + sqrt(Sym(3)))r2)^2 - (r1 + r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, x2))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-104256*sqrt(2)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) - 60192*sqrt(6)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) - 8688*sqrt(2)*pi^3*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) - 5016*sqrt(6)*pi^3*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 627*sqrt(2)*pi^4*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 362*sqrt(6)*pi^4*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 90288*sqrt(2)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 52128*sqrt(6)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 45144*sqrt(2)*pi^2*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 26064*sqrt(6)*pi^2*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + (-4*sqrt(黒積)*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3))/(-6 + sqrt(3)*pi) - (2*sqrt(6)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 12*sqrt(2)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)))/(36 - 3*pi^2))*(-69840*pi^2 - 40320*sqrt(3)*pi^2 - 2910*pi^4 - 1680*sqrt(3)*pi^4 - 83808 - 48384*sqrt(3) + 168*pi^5 + 97*sqrt(3)*pi^5 + 120960*pi + 69840*sqrt(3)*pi + 20160*pi^3 + 11640*sqrt(3)*pi^3))/(3*(-pi + 2*sqrt(3))*(-28*sqrt(3)*pi - 48*pi + 4*sqrt(3)*pi^2 + 7*pi^2 + 48*sqrt(3) + 84)^2), -4*sqrt(黒積)*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3))/(-6 + sqrt(3)*pi) - (2*sqrt(6)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 12*sqrt(2)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)))/(36 - 3*pi^2))
    ((-104256*sqrt(2)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) - 60192*sqrt(6)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) - 8688*sqrt(2)*pi^3*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) - 5016*sqrt(6)*pi^3*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 627*sqrt(2)*pi^4*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 362*sqrt(6)*pi^4*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 90288*sqrt(2)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 52128*sqrt(6)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 45144*sqrt(2)*pi^2*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 26064*sqrt(6)*pi^2*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + (4*sqrt(黒積)*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3))/(-6 + sqrt(3)*pi) - (2*sqrt(6)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 12*sqrt(2)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)))/(36 - 3*pi^2))*(-69840*pi^2 - 40320*sqrt(3)*pi^2 - 2910*pi^4 - 1680*sqrt(3)*pi^4 - 83808 - 48384*sqrt(3) + 168*pi^5 + 97*sqrt(3)*pi^5 + 120960*pi + 69840*sqrt(3)*pi + 20160*pi^3 + 11640*sqrt(3)*pi^3))/(3*(-pi + 2*sqrt(3))*(-28*sqrt(3)*pi - 48*pi + 4*sqrt(3)*pi^2 + 7*pi^2 + 48*sqrt(3) + 84)^2), 4*sqrt(黒積)*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3))/(-6 + sqrt(3)*pi) - (2*sqrt(6)*pi*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)) + 12*sqrt(2)*sqrt(黒積)*sqrt(-pi + 2*sqrt(3)))/(36 - 3*pi^2))

2番めのものが適解である。

res[2][1] |> simplify |> println  # 大円の半径 r1

   -sqrt(黒積)*(3*(12 - pi^2)*sqrt(-pi + 2*sqrt(3))*(-sqrt(3)*pi + 6)*(-104256*sqrt(2)*pi - 60192*sqrt(6)*pi - 8688*sqrt(2)*pi^3 - 5016*sqrt(6)*pi^3 + 627*sqrt(2)*pi^4 + 362*sqrt(6)*pi^4 + 90288*sqrt(2) + 52128*sqrt(6) + 45144*sqrt(2)*pi^2 + 26064*sqrt(6)*pi^2) + 2*(sqrt(-pi + 2*sqrt(3))*(-sqrt(3)*pi + 6)*(sqrt(6)*pi + 6*sqrt(2)) + 6*(12 - pi^2)*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3)))*(-11640*sqrt(3)*pi^3 - 20160*pi^3 - 69840*sqrt(3)*pi - 120960*pi - 97*sqrt(3)*pi^5 - 168*pi^5 + 48384*sqrt(3) + 83808 + 1680*sqrt(3)*pi^4 + 2910*pi^4 + 40320*sqrt(3)*pi^2 + 69840*pi^2))/(9*(12 - pi^2)*(pi - 2*sqrt(3))*(-sqrt(3)*pi + 6)*(-28*sqrt(3)*pi - 48*pi + 4*sqrt(3)*pi^2 + 7*pi^2 + 48*sqrt(3) + 84)^2)

res[2][2] |> factor |> println  # 小円の中心の x 座標 x2

   2*sqrt(黒積)*(sqrt(2)*sqrt(-pi + 2*sqrt(3)) + 2*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3)))/(-6 + sqrt(3)*pi)

黒積を与えて r1 を求める

res[2][1](黒積 => 1.2345).evalf()*2

   25.6917381798322

SymPy では r1 を求める式はこれ以上簡約化できないが，術では sqrt(黒積*(448sqrt(3) + 776)/(18sqrt(3) - 9pi)) という式が示されている。

黒積 = 1.2345
sqrt(黒積*(448sqrt(3) + 776)/(18sqrt(3) - 9pi))

   25.691738179835568

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   黒積 = 1.2345
   r2 = sqrt(2)*sqrt(黒積)/sqrt(-pi + 2*sqrt(3))
   r1 = -sqrt(黒積)*(3*(12 - pi^2)*sqrt(-pi + 2*sqrt(3))*(-sqrt(3)*pi + 6)*(-104256*sqrt(2)*pi - 60192*sqrt(6)*pi - 8688*sqrt(2)*pi^3 - 5016*sqrt(6)*pi^3 + 627*sqrt(2)*pi^4 + 362*sqrt(6)*pi^4 + 90288*sqrt(2) + 52128*sqrt(6) + 45144*sqrt(2)*pi^2 + 26064*sqrt(6)*pi^2) + 2*(sqrt(-pi + 2*sqrt(3))*(-sqrt(3)*pi + 6)*(sqrt(6)*pi + 6*sqrt(2)) + 6*(12 - pi^2)*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3)))*(-11640*sqrt(3)*pi^3 - 20160*pi^3 - 69840*sqrt(3)*pi - 120960*pi - 97*sqrt(3)*pi^5 - 168*pi^5 + 48384*sqrt(3) + 83808 + 1680*sqrt(3)*pi^4 + 2910*pi^4 + 40320*sqrt(3)*pi^2 + 69840*pi^2))/(9*(12 - pi^2)*(pi - 2*sqrt(3))*(-sqrt(3)*pi + 6)*(-28*sqrt(3)*pi - 48*pi + 4*sqrt(3)*pi^2 + 7*pi^2 + 48*sqrt(3) + 84)^2)
   x2 = 2*sqrt(黒積)*(sqrt(2)*sqrt(-pi + 2*sqrt(3)) + 2*sqrt(-2*pi - sqrt(3)*pi + 6 + 4*sqrt(3)))/(-6 + sqrt(3)*pi)
   @printf("黒積 = %g;  r1 = %g;  r2 = %g;  x2 = %g\n", 黒積, r1, r2, x2)
   @printf("大円の直径 = %g\n", 2r1)
   plot()
   circlef(x2 - r2, (1 + 1/√3)r2, r2/√3, :black)
   circlef(x2, r2, r2, :paleturquoise)
   circlef(x2 - 2r2, r2, r2, :paleturquoise)
   circlef(x2 -  r2, (1 + √3)r2, r2, :paleturquoise)
   circlef(0, r1, r1, :pink)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1, " r1    大円", :black, :left, :vcenter)
       point(x2, r2, " 小円:r2,(x2,r2)", :black, :left, :vcenter)
       point(x2 - r2, (1 + √3)r2, " 小円:r2,(x2-r2,(1+√3)r2))", :black, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その423）

算額（その423）

長野県諏訪市 諏訪大明神 天保12年(1841)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

直角三角形（鈎股弦）で，その面積が 210 歩，周と弦の和が 99 寸のとき，各辺の長さを求めよ。

鈎，股，弦の長さを 鈎, 股, 弦とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     a::positive, b::positive;

(a, b) = (210, 99)
eq1 = 股*鈎/2 - a
eq2 = 鈎 + 股 + 2弦 - b
eq3 = 鈎^2 + 股^2 - 弦^2

res = solve([eq1, eq2, eq3], (鈎, 股, 弦))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (-b/6 - sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3)
    (-b/6 - sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3)
    (-b/6 + sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3)
    (-b/6 + sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3)

鈎 < 股 なので，3 番目の組が適解である。

(a, b) = (210, 99)
(-b/6 - sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3) |> println
(-b/6 - sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 + 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b + 2*sqrt(12*a + b^2))^2)/6 - sqrt(12*a + b^2)/3, 2*b/3 + sqrt(12*a + b^2)/3) |> println
(-b/6 + sqrt(12*a + b^2)/3 - sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 + sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3) |> println
(-b/6 + sqrt(12*a + b^2)/3 + sqrt(-24*a + 5*b^2 - 4*b*sqrt(12*a + b^2))/6, -b/6 - sqrt(-72*a + (b - 2*sqrt(12*a + b^2))^2)/6 + sqrt(12*a + b^2)/3, 2*b/3 - sqrt(12*a + b^2)/3) |> println

   (-102.91912585224469, -4.080874147755303, 103.0)
   (-4.080874147755303, -102.91912585224469, 103.0)
   (20.0, 21.0, 29.0)
   (21.0, 20.0, 29.0)

(a, b) = (210, 99) のとき，(鈎, 股, 弦) = (20, 21, 29) である。

 

算額（その422）

算額（その422）

長野県 諏訪社 文化2年(1805)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

鈎股弦に大円と菱形が内接し，菱形に小円が内接している。小円の直径が　83.7 寸のとき，大円の直径を求めよ。

鈎，股，弦の長さを 鈎, 股, 弦
大円の半径と中心座標を r1, (股 - r1, r1)
小円の半径と中心座標を r2, (x2, r2)
菱形の一辺と股の交点座標を (a, 0) とする（図参照）
菱形の高さは 2r2 である（小円は菱形の上辺と下辺に接している）。
以下の連立方程式を立て，数値解を求める。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     r1::positive, r2::positive, x2::positive,
     a::positive, b::positive;

eq1 = 鈎^2 + 股^2 - 弦^2
eq2 = 鈎 + 股 - 弦 - 2r1
eq3 = distance(0, 0, 股, 鈎, x2, r2) - r2^2
eq4 = distance(a, 0, 股, 2r2, x2, r2) - r2^2
eq5 = (股 - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq6 = 2r2/(股 - a) - 鈎/股;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鈎, 股, 弦, r1, x2, a))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (鈎, 股, 弦, r1, x2, a) = u
   return [
       -弦^2 + 股^2 + 鈎^2,  # eq1
       -2*r1 - 弦 + 股 + 鈎,  # eq2
       -r2^2 + (r2 - 鈎*(r2*鈎 + x2*股)/(股^2 + 鈎^2))^2 + (x2 - 股*(r2*鈎 + x2*股)/(股^2 + 鈎^2))^2,  # eq3
       -r2^2 + (r2 - 2*r2*(a^2 - a*x2 - a*股 + 2*r2^2 + x2*股)/(a^2 - 2*a*股 + 4*r2^2 + 股^2))^2 + (x2 - (2*r2^2*(a - 2*x2 + 股) + x2*(a^2 - 2*a*股 + 4*r2^2 + 股^2))/(a^2 - 2*a*股 + 4*r2^2 + 股^2))^2,  # eq4
       (r1 - r2)^2 - (r1 + r2)^2 + (-r1 - x2 + 股)^2,  # eq5
       2*r2/(-a + 股) - 鈎/股,  # eq6
   ]
end;

r2 = 83.7/2
iniv = [big"180.0", 330.0, 378.0, 72.0, 150.0, 192.0]

res = nls(H, ini=iniv);
println(res);

   (BigFloat[177.7149360644993599352411051684090983213846737911418418853773103721415723113209, 347.4194693522861942938268302046185973814888150527085239226608605207901268261221, 390.2343990288826671178217820524254737950070308416933506657933666693103715275904, 67.4500031939514435556230766603011109539332290010785075711224021118106638049274, 173.709734676143097146913415102309298777291728644346174703749593832642708454533, 183.7922007121731238055940109044496826135959195331646665299783932361428838118555], true)

r2 = 41.85
r1 = 67.45;  x2 = 173.71;  鈎 = 177.715;  股 = 347.419;  弦 = 390.234;  a = 183.792

小円の直径が 83.7 寸のとき，大円の直径は 134.9 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 83.7/2
   (鈎, 股, 弦, r1, x2, a) = res[1]
   @printf("r2 = %g\n", r2)
   @printf("r1 = %g;  x2 = %g;  鈎 = %g;  股 = %g;  弦 = %g;  a = %g\n", r1, x2, 鈎, 股, 弦, a)
   @printf("大円の直径 = %g\n", 2r1)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(股 - r1, r1, r1)
   circle(x2, r2, r2, :blue)
   segment(a, 0, 股, 2r2)
   segment(股 - a, 2r2, 股, 2r2)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(a, 0, " a", :black, :right, :bottom, delta=delta)
       point(股 - a, 2r2, "(股-a,2r2)", :black, :right, :bottom, delta=delta)
       point(股 - r1, r1, "大円:r1,(股-r1,r1)", :red, :center, :bottom, delta=2delta)
       point(x2, r2, "小円:r2,(x2,r2)", :blue, :center, :bottom, delta=2delta)
       point(股, 2r2, " (股,2r2)", :black, :left, :vcenter)
       point(股, 鈎, " (股,鈎)", :black, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その421）

算額（その421）

長野県伊那市高遠町藤沢 諏訪大明神社 寶歴元年(1751)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

問題　直角三角形があり，その面積を S，(鈎 + 股 + 弦)/(股 - 鈎) = A として，各辺の長さを求めよ。

算額および中村も，S, A をもとにして鈎の 4 次方程式を解けばよいとしているが，ピタゴラスの定理を第 3 の条件式とすれば，SymPy で連立方程式を解くことで解が（一度に）求まる。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive, S::positive, A::positive;

eq1 = 股*鈎/2 - S
eq2 = (鈎 + 股 + 弦)/(股 - 鈎) - A
eq3 = 鈎^2 + 股^2 - 弦^2

res = solve([eq1, eq2, eq3], (鈎, 股, 弦))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (sqrt(S)*(-A*sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - A*sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), sqrt(S)*(A*sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - A*sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), 2*sqrt(S)*sqrt((A^3 - A - 2*sqrt(2*A^2 + 1))/(A^2 - 4))/sqrt(A))
    (sqrt(S)*(A*sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + A*sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), sqrt(S)*(-A*sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + A*sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((A^5 - 5*A^3 - 2*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((2*A^7 - 9*A^5 - 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 + 14*A^2*sqrt(2*A^2 + 1) + 4*A + 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), -2*sqrt(S)*sqrt(-(-A^3 + A + 2*sqrt(2*A^2 + 1))/(A^2 - 4))/sqrt(A))
    (sqrt(S)*(-A*sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + A*sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), sqrt(S)*(A*sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + A*sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), 2*sqrt(S)*sqrt((A^3 - A + 2*sqrt(2*A^2 + 1))/(A^2 - 4))/sqrt(A))
    (sqrt(S)*(A*sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - A*sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), sqrt(S)*(-A*sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - A*sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) + sqrt((A^5 - 5*A^3 + 2*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)) - sqrt((2*A^7 - 9*A^5 + 4*A^4*sqrt(2*A^2 + 1) + 3*A^3 - 14*A^2*sqrt(2*A^2 + 1) + 4*A - 8*sqrt(2*A^2 + 1))/(A^4 - 8*A^2 + 16)))/(sqrt(A)*(A^2 + 1)), -2*sqrt(S)*sqrt((A^3 - A + 2*sqrt(2*A^2 + 1))/(A^2 - 4))/sqrt(A))

3番めの組が適解である。

res[3][1] |> simplify |> display

res[3][2] |> simplify |> display

res[3][3] |> simplify |> display

解が適正であることを，例で見てみよう。
鈎, 股, 弦がそれぞれ，3.9, 5.2, 6.5 のとき，S = 10.14, A = 12 である。

(鈎, 股, 弦) = (3.9, 5.2, 6.5)
println("S = $(股*鈎/2),  A = $((鈎 + 股 + 弦)/(股 - 鈎))")

   S = 10.14,  A = 11.999999999999996

res[3][1], res[3][2], res[3][3] に S = 10.14, A = 12 を代入すると，鈎，股，弦が求まる。

println([res[3][i](S => 10.14, A => 12).evalf() for i in 1:3])

   Sym[3.90000000000000, 5.20000000000000, 6.50000000000000]

 

算額（その420）

算額（その420）

茨城県笠間市福原 吾國山上頂神（吾国山上頂神前） 文化7年(1810)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

長径 8 寸，短径 6 寸の菱形内に，甲円 2 個，乙円 2 個が入っている。甲円の直径が 4 寸のとき，乙円の直径を求めよ。

菱形の中心を原点とする。
甲円の半径と中心座標を r1, (x1, 0)
乙円の半径と中心座標を r2, (0, y2)
とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x1::positive, y2::positive;

eq1 = x1^2 + y2^2 - (r1 + r2)^2
eq2 = distance(0, 3, 4, 0, x1, 0) - r1^2
eq3 = distance(0, 3, 4, 0, 0, y2) - r2^2;
res = solve([eq1, eq2, eq3], (r2, x1, y2))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (4*(-16*sqrt(15)*r1^(5/2) - 900*sqrt(15)*sqrt(r1) + 16*r1^3 + 60*r1^2 + 900*r1 + 3375)/(9*(4*r1^2 + 225)), -(5*r1 - 12)/3, 20*sqrt(15)*sqrt(r1)/9 - 20*r1/9 - 16/3)
    (4*(4*r1 - 15)/9, -(5*r1 - 12)/3, 20*r1/9 - 16/3)

最初の組が適解である。

r1 = 4/2
(4*(-16*sqrt(15)*r1^(5/2) - 900*sqrt(15)*sqrt(r1) + 16*r1^3 + 60*r1^2 + 900*r1 + 3375)/(9*(4*r1^2 + 225)), -(5*r1 - 12)/3, 20*sqrt(15)*sqrt(r1)/9 - 20*r1/9 - 16/3)

   (0.48493231101926837, 0.6666666666666666, 2.3938346112259135)

r2 = 0.484932;  x1 = 0.666667;  y2 = 2.39383

乙円の直径は 0.969865 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 4//2
   (r2, x1, y2) = (4*(-16*sqrt(15)*r1^(5/2) - 900*sqrt(15)*sqrt(r1) + 16*r1^3 + 60*r1^2 + 900*r1 + 3375)/(9*(4*r1^2 + 225)), -(5*r1 - 12)/3, 20*sqrt(15)*sqrt(r1)/9 - 20*r1/9 - 16/3)
   @printf("r2 = %g;  x1 = %g;  y2 = %g\n", r2, x1, y2)
   @printf("乙円の直径 = %g\n", 2r2)
   plot([4, 0, -4, 0, 4], [0, 3, 0, -3, 0], color=:black, lw=0.5)
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle(0, y2, r2, :blue)
   circle(0, -y2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, 0, " x1", :red)
       point(0, y2, " y2", :blue)
       point(4, 0, " 4", :black, :left, :bottom)
       point(0, 3, " 3", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その419）

算額（その419）

東京都江東区木場（深川） 大丸屋内稲荷絵馬堂 文化7年(1810)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

弓形の中に正三角形 1 個，大円 2 個，中円 4 個，小円 4 個が入っている。大円の直径が 1 尺のとき，小円の直径を求めよ。

弓形を構成する外円の中心を原点とする。
弓形を構成する弦と外円の中心の距離を a とする（図参照）
外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (x1, a + r1)
中円の半径と中心座標を r2, (x21, a + r2), (x22, y22)
小円の半径と中心座標を r3, (x31, a + r3), (x32, y32)
として，以下の連立方程式の数値解を求める。

未知数は r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32 の11個であるが，条件式（方程式）は 12 個になった。
以下の eq4 〜 eq9 のうちどれか一つを使わないようにする（eq9 を除いた）。

include("julia-source.txt");

using SymPy

@syms r0::positive, a::positive, r1::positive, x1::positive,
     r2::positive, x21::positive, x22::positive, y22::positive,
     r3::positive, x31::positive, x32::positive, y32::positive;

eq1 = x1^2 + (a + r1)^2 - (r0 - r1)^2
eq2 = x21^2 + (a + r2)^2 - (r0 - r2)^2
eq3 = x22^2 + y22^2 - (r0 - r2)^2
eq4 = (x21 - x1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq5 = (x1 - x22)^2 + (y22 - a - r1)^2 - (r1 + r2)^2
eq6 = (x31 - x1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq7 = (x1 - x32)^2 + (y32 - a - r1)^2 - (r1 + r3)^2
eq8 = (x21 - x31)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq9 = (x22 - x32)^2 + (y22 - y32)^2 - (r2 + r3)^2
eq10 = distance(0, r0, (r0 - a)/sqrt(Sym(3)), a, x1, a + r1) - r1^2
eq11 = distance(0, r0, (r0 - a)/sqrt(Sym(3)), a, x22, y22) - r2^2
eq12 = distance(0, r0, (r0 - a)/sqrt(Sym(3)), a, x32, y32) - r3^2;

 

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32) = u
   return [
       x1^2 + (a + r1)^2 - (r0 - r1)^2,  # eq1
       x21^2 + (a + r2)^2 - (r0 - r2)^2,  # eq2
       x22^2 + y22^2 - (r0 - r2)^2,  # eq3
       (r1 - r2)^2 - (r1 + r2)^2 + (-x1 + x21)^2,  # eq4
       -(r1 + r2)^2 + (x1 - x22)^2 + (-a - r1 + y22)^2,  # eq5
       (r1 - r3)^2 - (r1 + r3)^2 + (-x1 + x31)^2,  # eq6
       -(r1 + r3)^2 + (x1 - x32)^2 + (-a - r1 + y32)^2,  # eq7
       (r2 - r3)^2 - (r2 + r3)^2 + (x21 - x31)^2,  # eq8
       #-(r2 + r3)^2 + (x22 - x32)^2 + (y22 - y32)^2,  # eq9
       -r1^2 + (a/4 - r0/4 + r1/4 + sqrt(3)*x1/4)^2 + (sqrt(3)*a/4 - sqrt(3)*r0/4 + sqrt(3)*r1/4 + 3*x1/4)^2,  # eq10
       -r2^2 + (-r0/4 + sqrt(3)*x22/4 + y22/4)^2 + (-sqrt(3)*r0/4 + 3*x22/4 + sqrt(3)*y22/4)^2,  # eq11
       -r3^2 + (-r0/4 + sqrt(3)*x32/4 + y32/4)^2 + (-sqrt(3)*r0/4 + 3*x32/4 + sqrt(3)*y32/4)^2,  # eq12
   ]
end;
r1 = 10//2
iniv = [big"100.0", 60, 32, 8, 53, 15, 100, 3, 47, 12, 95] ./ 3.2
iniv = [big"31.25", 18.75, 10.0, 2.5, 16.5625, 4.6875, 31.25, 0.9375, 14.6875, 3.75, 29.6875]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[25.49038105676657998080848135617869247549515697170824483702794809282177979870005, 12.74519052838329013851885544688169255263494103177866816560880000647679820897323, 10.24519052838328991096467576691497724140415363411739473165687289639408890927506, 2.470438416976302806278212382647506529008215409707534746382983461692078901157198, 17.27432762617317838252967445536201528319013903197721182970015560305202172736248, 4.539957388152644446167925002781799424476871763565147607613821577005623652425758, 22.56782103024110685639605573917644547722212244239117758536157640732734540923325, 0.851900255445560589676464780553194462776632582898376305125453038963129343031564, 14.37290237744824272406559296990917376944931384289685950439663092059990005842643, 4.58897484763492825949410623312361446943111492490923485552995147479476165882802, 19.24584397689835787931127610394921161788065249923856835686360465925916856690406], true)

r0 = 25.4904;  a = 12.7452;  x1 = 10.2452;  r2 = 2.47044;  x21 = 17.2743;  x22 = 4.53996;  y22 = 22.5678;  r3 = 0.8519;  x31 = 14.3729;  x32 = 4.58897;  y32 = 19.2458

小円の直径は 1.7038 寸である。

using Plots

function draw(more=false)
   pyplot(size=(1000, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 10//2
   (r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32) = res[1]
   @printf("r0 = %g;  a = %g;  x1 = %g;  r2 = %g;  x21 = %g;  x22 = %g;  y22 = %g;  r3 = %g;  x31 = %g;  x32 = %g;  y32 = %g\n", r0, a, x1, r2, x21, x22, y22, r3, x31, x32, y32)
   @printf("小円の直径 = %g\n", 2r3)
   plot()
   circle(0, 0, r0, beginangle=0, endangle=180)
   circle(x1, a + r1, r1, :blue)
   circle(-x1, a + r1, r1, :blue)
   circle(x21, a + r2, r2, :green)
   circle(-x21, a + r2, r2, :green)
   circle(x22, y22, r2, :green)
   circle(-x22, y22, r2, :green)
   circle(x31, a + r3, r3, :magenta)
   circle(-x31, a + r3, r3, :magenta)
   circle(x32, y32, r3, :magenta)
   circle(-x32, y32, r3, :magenta)
   segment(0, r0, (r0 - a)/√3, a)
   segment(0, r0, -(r0 - a)/√3, a)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, a + r1, " 大円:r1,(x1,a+r1)", :blue, :left, :vcenter)
       point(x21, a + r2, "中円:r2,(x21,a+r2) ", :green, :right, :vcenter)
       point(x22, y22, " 中円:r2,(x22,y22)", :green, :left, :vcenter)
       point(x31, a + r3, "小円:r3,(x31,a+r3)  ", :magenta, :right, :vcenter)
       point(x32, y32, " 小円:r3,(x32,y32)", :magenta, :left, :vcenter)
       point(0, a, " a", :black, :left, :top, delta=-delta)
       point(√(r0^2 - a^2), a, "(√(r0^2-a^2),a)", :black, :right, :top, delta=-delta)
       point(0, r0, " r0", :red, :left, :bottom)
       point(-(r0 - a)/√3, a, "(-(r0-a)/√3, a)", :black, :center, :top, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その418）

算額（その418）

福島県田村市船引町北移東鳥堂 東鳥堂（東鳥神社）観音堂 文久2年(1862)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

団扇の中に弦の中点を通る2線分がある。3 個の等円は外円に接し，弦と2線分にも接している。弦と矢の和が 16 寸，弦と矢の積が 48 平方寸のとき，等円の直径を求めよ。

等円の半径と中心座標を r, (x, a + r) とし，以下の方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r::positive, x::positive,
     a::positive, b::positive, 矢::positive, 弦::positive;

矢 = r0 - a
弦 = 2sqrt(r0^2 - a^2)
eq1 = 弦 + 矢 - 16
eq2 = 弦*矢 - 48
eq3 = distance(0, a, b, sqrt(r0^2 - b^2), x, a + r) - r^2
eq4 = distance(0, a, b, sqrt(r0^2 - b^2), 0, r0 - r) - r^2
eq5 = x^2 + (a + r)^2 - (r0 - r)^2;
# solve([eq1, eq2, eq3, eq4, eq5], (r, r0, x, a, b))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r, r0, x, a, b) = u
   return [
       -a + r0 + 2*sqrt(-a^2 + r0^2) - 16,  # eq1
       2*(-a + r0)*sqrt(-a^2 + r0^2) - 48,  # eq2
       -r^2 + (-b*(-a*r + b*x + r*sqrt(-b^2 + r0^2))/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2) + x)^2 + (a + r - (a*(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2) - (a - sqrt(-b^2 + r0^2))*(-a*r + b*x + r*sqrt(-b^2 + r0^2)))/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2))^2,  # eq3
       b^2*(a^2 + a*r - a*r0 - a*sqrt(-b^2 + r0^2) - r*sqrt(-b^2 + r0^2) + r0*sqrt(-b^2 + r0^2))^2/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2)^2 - r^2 + (-r + r0 - (-a^2*r + a^2*r0 + a*b^2 + 2*a*r*sqrt(-b^2 + r0^2) - 2*a*r0*sqrt(-b^2 + r0^2) + b^2*r - b^2*r0 - r*r0^2 + r0^3)/(a^2 - 2*a*sqrt(-b^2 + r0^2) + r0^2))^2,  # eq4
       x^2 + (a + r)^2 - (-r + r0)^2,  # eq5
   ]
end;

iniv = [big"1.6", 5, 3, 2.5, 2.5]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[1.499999999999999999999999999999999999999999999999999999999999999999999999999948, 6.499999999999999999999999999999999999999999999999999999999999999999999999999862, 3.000000000000000000000000000000000000000000000000000000000000000000000000000138, 2.499999999999999999999999999999999999999999999999999999999999999999999999999862, 2.594733192202055198398672253319262240463466167190260192229005823351113326367367], true)

r = 1.5;  r0 = 6.5;  x = 3;  a = 2.5;  b = 2.59473
矢 = 4;  弦 = 12
等円の直径 = 3

弦と矢の和が 16 寸，弦と矢の積が 48 平方寸のとき，弦は 12 寸，矢は 4 寸，等円の直径は 3 寸である。

算額に描かれている図は，「弦は 12 寸，矢は 4 寸」に対応していない。実際に図を描いてみると，とても「団扇」とは言えないものである。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r0, x, a, b) = res[1]
   矢 = r0 - a
   弦 = 2sqrt(r0^2 - a^2)
   @printf("r = %g;  r0 = %g;  x = %g;  a = %g;  b = %g\n", r, r0, x, a, b)
   @printf("矢 = %g;  弦 = %g\n", 矢, 弦)
   @printf("等円の直径 = %g\n", 2r)
   plot()
   circle(0, 0, r0)
   circle(0, r0 - r, r, :blue)
   circle(x, a + r, r, :blue)
   circle(-x, a + r, r, :blue)
   segment(-sqrt(r0^2 - a^2), a, sqrt(r0^2 - a^2), a)
   segment(0, a, b, sqrt(r0^2 - b^2))
   segment(0, a, -b, sqrt(r0^2 - b^2))
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x, a + r, "r,(x,a+r)", :blue, :center, delta=-delta)
       point(0, a, " a", :black)
       point(b, √(r0^2-b^2), "(b,√(r0^2-b^2))", :red, :left, :bottom, delta=delta)
       point(√(r0^2-a^2), a, "(√(r0^2-a^2),a)", :red, :right, :top, delta=-delta)
       point(0, r0 - r, " r0-r", :blue, :left, :vcenter)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;