算額（その291）

算額（その291）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(281)
長野県中野市片塩 大徳寺 天保9年(1838)

上下の甲円は外接しそれぞれ左右の甲円にも外接している。乙円 2 個は互いに外接し，下の甲円と容円および左右の甲円に外接している。容円は上の甲円に内接している。
甲円の直径が 125 寸のとき，容円の直径を求めよ。

半径，中心座標を以下のように決める。
甲円 r1, (0, r1), (x1, 0)
乙円 r2, (r2, y2)
容円 r3, (0, 2r1 - r3)

r1 = 125/2 であるが，r1 も未知数として以下の連立方程式を解き(x1, r2, y2, r3) を求める（それぞれの式に r1 が含まれる）。

include("julia-source.txt");

using SymPy

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

r1 = 125/2
eq1 = x1^2 + r1^2 - 4r1^2
eq2 = r2^2 + (2r1 -r3 -y2)^2 - (r2 + r3)^2
eq3 = r2^2 + (y2 + r1)^2 - (r1 + r2)^2
eq4 = (x1 - r2)^2 + y2^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3, eq4], (x1, r2, y2, r3))

   1-element Vector{NTuple{4, Sym}}:
    (sqrt(3)*r1, r1*(-92*sqrt(-11013 + 12723*sqrt(3)) - 86*sqrt(-3671 + 4241*sqrt(3)) + 4499 + 8998*sqrt(3))/13497, r1*(-13497 - 4499*sqrt(3) + 86*sqrt(-11013 + 12723*sqrt(3)) + 276*sqrt(-3671 + 4241*sqrt(3)))/13497, r1*(-sqrt(-58736/167281 + 67856*sqrt(3)/167281) - 62*sqrt(3)/409 + 627/409))

   x1 = 108.25318;  r2 = 24.13163;  y2 = 20.70279;  r3 = 42.34994
   容円の直径(2r3) = 84.69989

容円の半径は以下の式に甲円の半径を掛けたものになる。

@syms d
apart(res[1][4](r1 => 1), d) |> simplify

「算額への招待 追補1」 http://www.wasan.jp/syotai/syotai.html では以下の解が示されている。

A = ((1 + 2√Sym(3)) - 2sqrt(1+√Sym(3)))/3
@syms d
eq = ((12+2A+A^2) - 2(A + 2)sqrt(2A + 1))/(8 + 4A + A^2)

簡約化すると以下のようになるが，SymPy で求めた式より若干複雑である。

apart(eq, d) |> expand |> simplify   

いずれにせよ，甲円の直径が 125 寸のとき，容円の直径は 84.69988610132283 寸 である。

2*(125/2 * (-4*sqrt(-3671 + 4241*sqrt(3)) - 62*sqrt(3) + 627)/409)

   84.69988610132283

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 125 / 2
   (x1, r2, y2, r3) = (
       sqrt(3)*r1,
       r1*(-92*sqrt(-11013 + 12723*sqrt(3)) - 86*sqrt(-3671 + 4241*sqrt(3)) + 4499 + 8998*sqrt(3))/13497,
       r1*(-13497 - 4499*sqrt(3) + 86*sqrt(-11013 + 12723*sqrt(3)) + 276*sqrt(-3671 + 4241*sqrt(3)))/13497,
       r1*(-sqrt(-58736/167281 + 67856*sqrt(3)/167281) - 62*sqrt(3)/409 + 627/409)
   )
   @printf("x1 = %.5f;  r2 = %.5f;  y2 = %.5f;  r3 = %.5f\n", x1, r2, y2, r3)
   @printf("容円の直径(2r3) = %.5f\n", 2r3)
   plot()
   circle(0, r1, r1)
   circle(0, -r1, r1)
   circle(x1, 0, r1)
   circle(-x1, 0, r1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, 2r1 - r3, r3, :orange)
   if more
       point(0, r1, " r1", :red)
       point(0, -r1, " -r1", :red)
       point(x1, 0, " x1", :red)
       point(r2, y2, "(r2,y2)", :blue, :center, :bottom)
       point(0, 2r1 - r3, " 2r1-r3", :orange)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

Julia の挙動不審

Julia で 2 * *(4 + 6)  を入力すると，エラーなしで  2 * (4 + 6) を計算する。

julia> 2 * *(4 + 6)

20

算額（その290）

算額（その290）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(249)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

長方形内に右下の頂点で長辺に接する円弧の一部があり，他に大円，中円，小円がある。それぞれは図のように互いに接している。
中円の直径が 500 寸のとき，大円の直径を求めよ。

長方形の長辺の長さを a とする。
その一部が円弧である円の半径，中心座標を r0, (0, r0) とする。
長方形の短辺の長さは sqrt(r0^2 - a^2) である。
大円の半径，中心座標を r1, (a - r1, r1)
中円の半径，中心座標を r2, (x2, x2)
小円の半径，中心座標を r3, (x3, y3) および (a - r3, r3) とする。
以下の連立方程式を nlsolve() で解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r0::positive, r1::positive, r2::positive,
     x2::positive, r3::positive, x3::positive, y3::positive;

r2 = 500 // 2
eq1 = (a - r1)^2 + (r0 - r1)^2 - (r0 + r1)^2
eq2 = x3^2 + (r0 - y3)^2 - (r0 + r3)^2
eq3 = x2^2 + (r0 - r2)^2 - (r0 + r2)^2
eq4 = 2(r1 - r3)^2 - (r1 + r3)^2
eq5 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq6 = (a - r1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq7 = (x3 - x2)^2 + (y3 - r2)^2 - (r2 + r3)^2;

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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, r0, r1, x2, r3, x3, y3) = u
   return [
       (a - r1)^2 + (r0 - r1)^2 - (r0 + r1)^2,  # eq1
       x3^2 - (r0 + r3)^2 + (r0 - y3)^2,  # eq2
       x2^2 + (r0 - 250)^2 - (r0 + 250)^2,  # eq3
       2*(r1 - r3)^2 - (r1 + r3)^2,  # eq4
       (r1 - 250)^2 - (r1 + 250)^2 + (a - r1 - x2)^2,  # eq5
       (-r1 + y3)^2 - (r1 + r3)^2 + (a - r1 - x3)^2,  # eq6
       -(r3 + 250)^2 + (-x2 + x3)^2 + (y3 - 250)^2,  # eq7
      ]
end;

iniv = [big"3100.0", 3900, 450, 2000, 75, 2100, 550]
res = nls(H, ini=iniv);

names = [" a", "r0", "r1", "x2", "r3", "x3", "y3"]
for (i, name) in enumerate(names)
   @printf("%s = %.6f\n", name, res[1][i])
end

    a = 3086.670209
   r0 = 3867.992287
   r1 = 449.500799
   x2 = 1966.721202
   r3 = 77.122145
   x3 = 2118.355531
   y3 = 539.855012

中円の半径は 449.500799 で，直径はほぼ 899 寸である。

using Plots

function diamond(x0, y0, a)
   plot!([1, 0, -1, 0, 1] ./ √2 .* a .+ x0, [0, 1, 0, -1, 0] ./ √2 .* a .+ y0)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r0, r1, x2, r3, x3, y3) = res[1]
   r2 = 500 / 2
   b = r0 - sqrt(r0^2 - a^2)
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(0, r0, r0, beginangle=270, endangle=360)
   circle(a - r1, r1, r1, :blue)
   circle(x2, r2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(a - r3, r3, r3, :magenta)
   if more
       point(a, b, "(a,b)")
       point(0, r0, " r0", :red)
       point(a - r1, r1, "大円:r1\n(a-r1,r1)", :blue, :center, :bottom)
       point(x2, r2, "中円:r2 \n(x2,r2) ", :green, :right, :top)
       point(x3, y3, "小円:r3  \n(x3,y3)  ", :magenta, :right, :bottom)
       point(a - r3, r3, "小円:r3  \n(a-r3,r3)  ", :magenta, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false, xlims=(0, 1.1Float64(a)), ylims=(-20, 1.1Float64(b)))
   end
end;

算額（その289）

算額（その289）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(250)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

長方形内に斜線を描き，大円，小円が2個ずつ入っている。大円は長方形と斜線に接している。小円は長方形，大円，斜線に接している。
長方形の短辺の長さが 360 寸，小円の直径が 40 寸のとき，長方形の長辺の長さを求めよ。

長方形の長辺の長さを a とする。短辺の長さは大円の直径に等しい。

大円の半径と中心座標を r1, (r1, r1)
小円の半径と中心座標を r2, (x, 2r1 - r2) とする。
斜線と長方形の長辺の交点の座標を (b, 2r1), (a - b, 0) とする。
以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r1::positive, r2::positive, x::positive;

r1 = 360 // 2
r2 = 40 // 2
eq1 = (x - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = distance(b, 2r1, a - b, 0, r1, r1) - r1^2
eq3 = distance(b, 2r1, a - b, 0, x, 2r1 - r2) - r2^2

res = solve([eq1, eq2, eq3], (a, b, x))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (735, 315, 300)
    (960, 240, 300)

最初の組が適解である。すなわち，長方形の長辺は 735 寸である。

using Plots

function diamond(x0, y0, a)
   plot!([1, 0, -1, 0, 1] ./ √2 .* a .+ x0, [0, 1, 0, -1, 0] ./ √2 .* a .+ y0)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, x) = res[1]
   println("a = $a")
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(a - r1, r1, r1)
   circle(x, 2r1 - r2, r2, :blue)
   circle(a - x, r2, r2, :blue)
   segment(a - b, 0, b, 2r1, :green)
   if more
       point(r1, r1, "(r1,r1)", :red)
       point(x, 2r1 - r2, "(x,2r1-r2)  ", :blue, :right)
       point(b, 2r1, "(b,2r1)", :black, :left, :bottom)
       point(a, 0, " a", :black, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その288）

算額（その288）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(250)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

円内に 2 本の平行な弦と 3 個の正方形が入っている。正方形の一辺の長さが 1 寸のとき，円の径を求めよ。

円の半径を R とする。正方形の一辺の長さを a とすれば，対角線の長さの半分は a/√2 = √2a/2 で，これを b とする。
⊿OAB において，OB = R, OA = 3b - R, AB = 2b

以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, a::positive, b::positive

b = a/√Sym(2)
eq = (3b-R)^2 + (2b)^2 - R^2
solve(eq, R)[1] |> println  # 半径 R

   13*sqrt(2)*a/12

直径は a*13√2/6 で，a = 1 寸のとき，円の直径は 3.0641293851417064 寸。

13√2/6

   3.0641293851417064

using Plots

function diamond(x0, y0, a)
   plot!([1, 0, -1, 0, 1] ./ √2 .* a .+ x0, [0, 1, 0, -1, 0] ./ √2 .* a .+ y0)
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   b = a/√2
   R = a*13*sqrt(2)/12
   println("R = $R  2R = $(2R)")
   plot()
   circle(0, 0, R)
   diamond(0, R - b, a)
   diamond(b, R - 3b, a)
   diamond(-b, R - 3b, a)
   y = R - 2b
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :blue)
   y = R - 4b
   x = sqrt(R^2 - y^2)
   segment(-x, y, x, y, :blue)
   if more
       point(0, R - b, "  R-b")
       point(0, R - 3b, "  A")
       point(0, R - 3b, "R-3b  ", :green, :right)
       point(2b, 0, "2b", :green, :center)
       point(R, 0, " R", :green)
       point(0, 0, " O")
       point(2b, R - 3b, " B")
       plot!([0, 2b, 0, 0], [0, R - 3b, R - 3b, 0], color=:gray, linestyle=:dot)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その287）

算額（その287）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(251)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

鉤股弦で，
股×弦 - 鉤×弦 = 425 平方寸
弦 - 鉤 = 18 寸
のとき，鉤，股，弦を求めよ。

以下の連立方程式を解く。
(鉤^2 + 股^2) - 弦^2 を使わないのは，虚数解を除くため。

include("julia-source.txt");

using SymPy
@syms 鉤::positive, 股::positive, 弦::positive;

eq1 = (股*弦 - 鉤*弦) - 425
eq2 = (弦 - 鉤) - 18
eq3 = (鉤^2 + 股^2) - 弦^2
eq3 = sqrt(鉤^2 + 股^2) - 弦
res = solve([eq1, eq2, eq3], (鉤, 股, 弦));

解は 2 組ある。

(鉤, 股, 弦) = res[1];
鉤.evalf(), 股.evalf(), 弦.evalf()

   (7.00000000000000, 24.0000000000000, 25.0000000000000)

(鉤, 股, 弦) = res[2];
鉤.evalf(), 股.evalf(), 弦.evalf()

   (25.3934598546653, 35.1875625010877, 43.3934598546653)

---

算額への招待 第4章 で術および現代解法が紹介されているが，方程式の係数に誤りがあるため，a=7, 25.3934598546653 の片方しか求めていない。

@syms a, b, c, A, B
A = 425
B = 18
eq1 = a^2 + b^2 - c^2
eq2 = b*c - a*c - A
eq3 = c - a - B;

solve(eq3, c)[1] |> println  # c

   a + 18

eq12 = eq2(c => a + 18)
eq12 |> println

   -a*(a + 18) + b*(a + 18) - 425

solve(eq12, b)[1] |> println  # b

   (a^2 + 18*a + 425)/(a + 18)

eq11 = eq1(c => a + 18, b => (a^2 + 18*a + 425)/(a + 18))
eq111 = eq11 |> simplify
eq111|> println

   (a^4 - 446*a^2 - 8028*a + 75649)/(a^2 + 36*a + 324)

factor(numerator(eq111)) |> println

   (a - 7)*(a^3 + 7*a^2 - 397*a - 10807)

pyplot(size=(200, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(eq111, xlims=(0, 30))
hline!([0])

res = solve(eq111)
res[1].evalf() |> println
res[2].evalf() |> println
res[3].evalf() |> println
res[4].evalf() |> println

   7.00000000000000
   -16.1967299273327 - 12.7768525871672*I
   -16.1967299273327 + 12.7768525871672*I
   25.3934598546653

件の現代解法では定数項が + 10807 となっており，もう一つの解が見逃されている。

A^2 - B^4 + (2A*B - 4B^3)a + (2A - 4B^2)a^2 + a^4 |> factor |> println

   (a - 7)*(a^3 + 7*a^2 - 397*a - 10807)

SymPy の solve() は，ここに展開したような式変形と求解を間違いなくやってくれている。

 

算額（その286）

算額（その286）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(252)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

等脚台形内に，甲，乙，丙の 5 個の円がある。それぞれ台形に内接したり，互いに外接しあっている。
甲円の直径が 1 寸のとき，乙円の直径を求めよ。

台形の上底と下底の長さを 2b, 2a，高さを c とする。
甲円の半径と中心座標 r1, (0, r1)
乙円の半径と中心座標 r2, (0, r2 + 2r1)
上の丙円の半径と中心座標 r3, (0, r3 + 2r1 + 2r2)
右の丙円の半径と中心座標 r3, (x3, r3)
c = 2(r1 + r2 + r3) である。

以下の連立方程式を解く。なお，円の中心から直線までの平方距離は分数式になるので分子 = 0 を解けば数式は若干短くなる。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive,
     r3::positive, x3::positive;

r1 = 1//2
c = 2(r1 + r2 + r3)
eq1 = x3^2 + (r3 - r1)^2 - (r3 + r1)^2
eq2 = x3^2 + (r2 + 2r1 - r3)^2 - (r2 + r3)^2
eq3 = distance(a, 0, b, c, 0, c - r3) - r3^2
eq3 = -2*a*b*r3 - 2*b^2*r2 - b^2 + 2*r2*r3^2 + 2*r3^3 + r3^2
eq4 = distance(a, 0, b, c, 0, r2 + 2r1) - r2^2
eq4 = 4*a^2*r2*r3 + 4*a^2*r3^2 + 4*a*b*r2^2 + 4*a*b*r2*r3 + 2*a*b*r2 + 4*a*b*r3 + 2*b^2*r2 + b^2 - 4*r2^4 - 8*r2^3*r3 - 4*r2^3 - 4*r2^2*r3^2 - 4*r2^2*r3 - r2^2
eq5 = distance(a, 0, b, c, x3, r3) - r3^2
eq5 = 2*a^2*r2 + a^2 + 2*a*b*r3 - 4*a*r2*x3 - 2*a*r3*x3 - 2*a*x3 - 2*b*r3*x3 - 2*r2*r3^2 + 2*r2*x3^2 - 2*r3^3 - r3^2 + 2*r3*x3^2 + x3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5])#, (a, r1, r2, r3, x3, r6, x6, y6))

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")

   x3^2 + (r3 - 1/2)^2 - (r3 + 1/2)^2,  # eq1
   x3^2 - (r2 + r3)^2 + (r2 - r3 + 1)^2,  # eq2
   -2*a*b*r3 - 2*b^2*r2 - b^2 + 2*r2*r3^2 + 2*r3^3 + r3^2,  # eq3
   4*a^2*r2*r3 + 4*a^2*r3^2 + 4*a*b*r2^2 + 4*a*b*r2*r3 + 2*a*b*r2 + 4*a*b*r3 + 2*b^2*r2 + b^2 - 4*r2^4 - 8*r2^3*r3 - 4*r2^3 - 4*r2^2*r3^2 - 4*r2^2*r3 - r2^2,  # eq4
   2*a^2*r2 + a^2 + 2*a*b*r3 - 4*a*r2*x3 - 2*a*r3*x3 - 2*a*x3 - 2*b*r3*x3 - 2*r2*r3^2 + 2*r2*x3^2 - 2*r3^3 - r3^2 + 2*r3*x3^2 + x3^2,  # eq5

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)
   (r2, r3, x3, a, b) = u
   return [
       x3^2 + (r3 - 1/2)^2 - (r3 + 1/2)^2,  # eq1
       x3^2 - (r2 + r3)^2 + (r2 - r3 + 1)^2,  # eq2
       -2*a*b*r3 - 2*b^2*r2 - b^2 + 2*r2*r3^2 + 2*r3^3 + r3^2,  # eq3
       4*a^2*r2*r3 + 4*a^2*r3^2 + 4*a*b*r2^2 + 4*a*b*r2*r3 + 2*a*b*r2 + 4*a*b*r3 + 2*b^2*r2 + b^2 - 4*r2^4 - 8*r2^3*r3 - 4*r2^3 - 4*r2^2*r3^2 - 4*r2^2*r3 - r2^2,  # eq4
       2*a^2*r2 + a^2 + 2*a*b*r3 - 4*a*r2*x3 - 2*a*r3*x3 - 2*a*x3 - 2*b*r3*x3 - 2*r2*r3^2 + 2*r2*x3^2 - 2*r3^3 - r3^2 + 2*r3*x3^2 + x3^2,  # eq5
      ]
end;

iniv = [1.3, 0.7, 1.2, 2.1, 0.5]
res = nls(H, ini=iniv);
println(res);

   ([1.3090169943749475, 0.6909830056250524, 1.1755705045849463, 2.1266270208801, 0.5020285397155682], true)

乙円の半径 = 1.3090169943749475 なので，直径はほぼ 2.618 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (r2, r3, x3, a, b) = res[1]
   c = 2(r1 + r2 + r3)
   plot([a, b, -b, -a, a], [0,c, c, 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1, :blue)
   circle(0, r2 + 2r1, r2, :orange)
   circle(0, c - r3, r3)
   circle(x3, r3, r3)
   circle(-x3, r3, r3)
   if more
       str = @sprintf("r2 ≒ %.5f", 2r2)
       point(0, r1, " r1\n 甲:r1", :blue)
       point(0, r2 + 2r1, " r2+2r1\n 乙:r2\n $str", :orange)
       point(0, c - r3, " r3+2r2+2r1\n 丙:r3", :red)
       point(x3, r3, " 丙:r3\n (x3,r3)", :red)
       point(b, c, "(b,2(r1+r2+r3))", :black, :left, :bottom)
       point(a, 0, "a ", :black, :right, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その285）

算額（その285）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(253)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

外円の中に弦を引き，等円を 6 個入れる。そのうち 3 個は弦の下，2 個は弦の下に置き，それぞれが弦に接し，互いに外接する。また弦の下の 3 個の円の左右の等円は外円に内接する。残り 1 個は弦の上の2個に外接し，外円に内接する。

等円の径が 1 寸のとき，外円の径を求めよ。

外円，等円の半径を R, r，弦の位置を y とする（y < 0, 外円の中心から弦までの距離は -y）。
等円の中心座標は図に示すようになる。
r を変数のままにして，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive, y::negative;

eq1 = r^2 + ((R - r) - (r + y))^2 - 4r^2
eq2 = 4r^2 + (y - r)^2 - (R - r)^2

res = solve([eq1, eq2], (R, y))

   1-element Vector{Tuple{Sym, Sym}}:
    (3*r*(4 - 5*sqrt(3)/2)/5 + 18*r/5, r*(4 - 5*sqrt(3)/2))

解の式を簡約化する。

res[1][1] |> simplify |> println  # 外円の半径
res[1][2] |> simplify |> println  # 弦の位置（y 座標）

   3*r*(4 - sqrt(3))/2
   r*(8 - 5*sqrt(3))/2

等円の直径が 1 寸のとき，外円の直径は 3*(4 - sqrt(3))/2 ≒ 3.401923788646684 ≒ 3寸4分

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (R, y) = (3*r*(4 - sqrt(3))/2, r*(8 - 5*sqrt(3))/2)
   x = sqrt(R^2 - y^2)
   plot()
   circle(0, 0, R, :black)
   circle(0, R - r, r)
   circle(r, r + y, r)
   circle(-r, r + y, r)
   circle(0, y - r, r)
   circle(2r, y - r, r)
   circle(-2r, y - r, r)
   segment(-x, y, x, y, :blue)
   if more
       annotate!(0.65, 1.05, text(@sprintf("R = %.3f\ny = %.3f\nx = %.3f\nr = %.3f", R, y, x, r), 10, :left))
       point(x, y, "(x,y)", :blue, :right, :bottom)
       point(2r, y - r, "(2r,y-r)", :red)
       point(r, r + y, "(r,r+y)", :red)
       point(0, R - r, "(0,R-r)", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その284）

算額（その284）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(247)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

菱形の長，平（横の対角線と縦の対角線），面（一辺の長さ）の和に，面と長の再乗冪（3乗）を加えると　901430720，菱形の面積は 42504 になる。このとき，長，平，面を求めよ。

「面と長の再乗冪」は「面の再乗冪と長の再乗冪」ではない。
連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms 長::positive, 平::positive, 面::positive;

eq1 = (長 + 平 + 面) + (面 + 長^3) - 901430720
eq2 = 長 * 平 / 2 - 42504
eq3 = (長^2 + 平^2)/4 - 面^2;

res = solve([eq1, eq2, eq3], (長, 平, 面))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (966, 88, 485)

長 = 966，平 = 88，面 = 485 である。

算額（その283）

算額（その283）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(246)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

大円内に，正三角形，元円，亨円，利円をそれぞれ 3 個，貞円を 6 個が入っている。それぞれは，大円，正三角形，隣り合う円に接している。

大円の半径，中心座標 R, (0, 0)
元円の半径，中心座標 r1, (x1, r1)
亨円の半径，中心座標 r2, (x2, r2)
利円の半径，中心座標 r3, (0, R - r3)
貞円の半径，中心座標 r4, (x4, R - 2r3 + r4
ここで，
r2 = 16 // 2
r3 = (1//2 - √3/4)R
x1 = (R - r1)√3/2
x2 = (R - 2r1 - r2)√3/2
である。

図のように記号を定め，連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, x1::positive, r1::positive, x2::positive, r2::positive,
   r3::positive, x4::positive, r4::positive;

r2 = 16 // 2
r3 = (1//2 - sqrt(Sym(3))/4)R
x1 = (R - r1)sqrt(Sym(3))/2
x2 = (R - 2r1 - r2)sqrt(Sym(3))/2
eq0 = 2r2 + r2 + 2r1 - R
eq1 = (x1 - x2)^2 - 4r1*r2
eq2 = x4^2 - 4r3*r4
eq3 = x1^2 + r1^2 - (R - r1)^2
eq4 = x4^2 + (R - 2r3 + r4)^2 - (R - r4)^2;

変数は R, r1, r4, x4 の 4 個であるが，条件式は5個以上考えられる。この内，適切なものを選択しないと solve() で解が得られないことがある。eq0, eq1, eq2, eq4 の 4 連立方程式を解くことで 2 通りの解が得られるが，最初の解が適切である。また，得られた式が複雑なものもあるが，simplify(), apart() などで単純な数式になる。また evalf() で数値として求めてもよい。

res = solve([eq0, eq1, eq2, eq4], (R, r1, r4, x4))

   2-element Vector{NTuple{4, Sym}}:
    (-64*sqrt(3)*sqrt(4*sqrt(3) + 7)/3 + 152/3 + 128*sqrt(4*sqrt(3) + 7)/3, -32*sqrt(3)*sqrt(4*sqrt(3) + 7)/3 + 40/3 + 64*sqrt(4*sqrt(3) + 7)/3, (41361027804795656 + 23879800537058368*sqrt(3) + 26321303246238494*sqrt(4*sqrt(3) + 7) + 15196611514637565*sqrt(3)*sqrt(4*sqrt(3) + 7))/(6*(1385331749802026 + 799821658665135*sqrt(3))*sqrt(4*sqrt(3) + 7)), sqrt(-1700*sqrt(3)/9 + 1216/(9*sqrt(4*sqrt(3) + 7)) + 3400/9))
    (-128*sqrt(4*sqrt(3) + 7)/3 + 152/3 + 64*sqrt(3)*sqrt(4*sqrt(3) + 7)/3, -64*sqrt(4*sqrt(3) + 7)/3 + 40/3 + 32*sqrt(3)*sqrt(4*sqrt(3) + 7)/3, (-41361027804795656 - 23879800537058368*sqrt(3) + 26321303246238494*sqrt(4*sqrt(3) + 7) + 15196611514637565*sqrt(3)*sqrt(4*sqrt(3) + 7))/(6*(1385331749802026 + 799821658665135*sqrt(3))*sqrt(4*sqrt(3) + 7)), sqrt(-1700*sqrt(3)/9 - 1216/(9*sqrt(4*sqrt(3) + 7)) + 3400/9))

貞円の半径は 4.5 である。元の単位で，丁円の直径は 9 寸である。

names = ("R", "r1", "r4", "x4")
i = 1
for (j, name) in enumerate(names)
   println("$name = $(res[i][j].evalf())")
end

   R = 72.0000000000000
   r1 = 24.0000000000000
   r4 = 4.50000000000000
   x4 = 9.31748562369075

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, r1, r4, x4) = res[1]
   r2 = 16 // 2
   r3 = (1//2 - sqrt(Sym(3))/4)R
   x1 = (R - r1)sqrt(Sym(3))/2
   x2 = (R - 2r1 - r2)sqrt(Sym(3))/2
   θ = 0:60:360
   x = [R*cosd(deg) for deg in θ]
   y = [R*sind(deg) for deg in θ]
   @printf("R = %.6f;  r1 = %.6f;  r4 = %.6f;  v4 = %.6f\n", R, r1, r4, x4)
   @printf("2r4 = %.6f\n", 2r4)
   plot([x[1], x[4], x[5], x[2], x[3], x[6], x[1]], [y[1], y[4], y[5], y[2], y[3], y[6], y[1]],
       color=:black, lw=0.5)
   circle(0, 0, R, :black)
   rotate(x1, r1, r1)
   rotate(x2, r2, r2, :blue)
   rotate(0, R - r3, r3, :green)
   rotate(x4, R - 2r3 + r4, r4, :magenta)
   rotate(-x4, R - 2r3 + r4, r4, :magenta)
   if more
       point(0, R, " R", :black, :left, :bottom)
       point(R, 0, "R ", :black, :right, :bottom)
       point(x1, r1, " 元:r1\n (x1,r1)", :red, :left, :top)
       point(x2, r2, " 亨:r2\n (x2,r2)", :blue, :left, :bottom)
       point(0, R - r3, "利:r3 \nR-r3 ", :green, :right, :top)
       point(x4, R - 2r3 + r4, " 貞:r4\n (x4,R-2r3+r4)", :green, :left, :top)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その282）

算額（その282）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(246)
長野県諏訪市中洲 諏訪大社上社 明治12年(1897)

外円の中に菱形，その中に正方形がある。それらに等円が図のように接している。等円の径が 8453 寸のとき，外円の径を求めよ。

正方形の一辺の長さを 2a，等円の半径を r，外円の半径，中心座標を R, (0, 0) とする。
上の等円の半径，中心座標は r, (0, R - r)
右の等円の半径，中心座標は r, (a + r, 0) である。

図のように記号を定め，連立方程式を解く。
等円の半径 r を未知数のまま解く。

include("julia-source.txt");

using SymPy
@syms R, r, a

r = 8453 // 2

eq1 = (R-2r) / R - a / (R - a)
eq2 = a/sqrt(a^2 + (R - a)^2) - r / (R - a - r)
res = solve([eq1, eq2], (R, a))

   1-element Vector{Tuple{Sym, Sym}}:
    (r + sqrt(3)*r + 2*sqrt(r^2), sqrt(3)*r)

外円の半径は　r + √3r + 2*sqrt(r^2) = (3 + √3) * r である。
直径が 8453 寸のとき，外円の直径は 40000.02547637972 寸である。約 40000 寸。

R = 8453 / 2; (3 + √3) * r * 2

   40000.02547637971

   R = 20000.012738;  a = 7320.512738

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 8453/2
   (R, a) = (r + sqrt(3)*r + 2*sqrt(r^2), sqrt(3)*r)
   @printf("R = %.6f;  a = %.6f\n", R, a)
   # @printf("2r = %.6f\n", 2r)
   plot([R, 0, -R, 0, R], [0, R - 2r, 0, - R + 2r, 0], color=:blue, lw=0.5)
   plot!([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:green, lw=0.5)
   circle(0, 0, R, :black)
   circle(0, R - r, r)
   circle(0, -R + r, r)
   circle(a + r, 0, r, :red)
   circle(-a - r, 0, r, :red)
   if more
       point(0, R - r, " R-r")
       point(0, a, " a", :green, :left, :bottom)
       point(a, 0, "a ", :green, :right, :bottom)
       point(a + r, 0, " a+r", :green, :left, :bottom)
       point(0, R, " R", :green, :left, :bottom)
       point(0, R - 2r, "   R-2r")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その279）

算額（その279）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(239)
長野県松本市筑摩 筑摩神社 明治6年(1873)

甲円 4 個，乙円 2 個が互いに外接し，3 個の円は直線に接している。
乙円の径が 1 寸のとき，甲円の径を求めよ。

甲円の半径を r1 とする。右の甲円の中心は (√3r1, 2r1) である。
乙円の半径を r2 とする。右の乙円の中心を (x2, r2) である。
以下の連立方程式を r2 を変数のままで解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, y1::positive, r2::positive, x2::positive;

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

   2-element Vector{Tuple{Sym, Sym}}:
    ((2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)^2/(4*r2), 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)
    ((-2*sqrt(6)*r2/3 + 2*sqrt(3)*r2/3)^2/(4*r2), -2*sqrt(6)*r2/3 + 2*sqrt(3)*r2/3)

   r1 = 0.971405;  x1 = 1.682522;  y1 = 1.942809;  x2 = 1.393847
   2r1 = 1.942809

最初の解が適解である。

甲円の半径は (2*sqrt(Sym(3))*r2/3 + 2*sqrt(Sym(6))*r2/3)^2/(4*r2) を簡約化して r2*(2*sqrt(Sym(2)) + 3)/3 である。
すなわち，「甲円の半径 = (√8/3 + 1) × 乙円の半径」である。
よって，乙円の直径が 1 寸の場合，甲円の直径は √8/3 + 1 = 1.9428090415820636 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1//2
   (r1, x2) = ((2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)^2/(4*r2), 2*sqrt(3)*r2/3 + 2*sqrt(6)*r2/3)
   y1 = 2r1
   x1 = sqrt(Sym(3))r1
   @printf("r1 = %.6f;  x1 = %.6f;  y1 = %.6f;  x2 = %.6f\n", r1, x1, y1, x2)
   @printf("2r1 = %.6f\n", 2r1)
   plot()
   circle(0, r1, r1, :blue)
   circle(0, 3r1, r1, :blue)
   circle(x1, y1, r1, :blue)
   circle(-x1, y1, r1, :blue)
   circle(x2, r2, r2)
   circle(-x2, r2, r2)
   if more
       point(0, r1, " r1")
       point(0, 2r1, " 2r1")
       point(0, 3r1, " 3r1")
       point(0, 4r1, " 4r1")
       point(x1, y1, " 甲円:r1\n (x1,y1)", :blue, :center, :top)
       point(x2, r2, " 乙円:r2\n (x2,r1)", :red, :center, :top)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その281）

算額（その281）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(240)
長野県松本市筑摩 筑摩神社 明治6年(1873)

等脚台形に対角線を引き，時計回りに 90° 回転させ上半分だけを考える。対角線で区切られた2つの区域に対角線及び辺に接する等円を入れる。元の台形の高さ a を 12 寸，等円の径 r を 4 寸とするとき，元の台形の底辺の半分 b はいかほどか。

図のように記号を定め，連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r::positive, a::positive, b::positive, c::positive, d::positive, y1::positive, dummy;
(r, a) = (2, 12)
eq1 = distance(0, 0, a, c, r, y1) - r^2
eq2 = distance(0, 0, a, c, a - r, r) - r^2
eq3 = distance(0, b, d, 0, r, y1) - r^2
eq4 = distance(0, b, d, 0, a - r, r) - r^2;

res = solve([eq1, eq2, eq3, eq4], (b, c, d, y1))

   3-element Vector{NTuple{4, Sym}}:
    (7, 5, 28/3, 3)
    (13/4 - sqrt(65)/4, 5, 26 - 2*sqrt(65), 3)
    (sqrt(65)/4 + 13/4, 5, 2*sqrt(65) + 26, 3)

最初の組が適解である。したがって，元の台形の底辺の半分は 7 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, c, d, y1) = res[1]
   @printf("b = %.6f;  c = %.6f;  d = %.6f;  y1 = %.6f\n", b, c, d, y1)
   # @printf("2r = %.6f\n", 2r)
   plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:black, lw=0.5)
   segment(0, 0, a, c)
   segment(0, b, d, 0)
   circle(r, y1, r, :blue)
   circle(a - r, r, r, :blue)
   if more
       point(r, y1, " (r,y1)")
       point(a - r, r, " (a-r,r)")
       point(a, 0, " a", :green, :left, :bottom)
       point(0, b, " b", :green, :left, :bottom)
       point(a, c, " (a,c)")
       point(d, 0, "  d", :green, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その280）

算額（その280）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(240)
長野県松本市筑摩 筑摩神社 明治6年(1873)

鉤股弦（直角三角形）に 2 本の斜線を入れ，区切られた領域に 2 個の等円を入れる。鉤が 5 寸のとき，等円の直径を求めよ。

図のように記号を定め，連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r::positive, x1::positive, y1::positive,
   x::positive, y::positive;

eq1 = distance(0, 5, a, 0, r, y1) - r^2
eq2 = distance(0, 5, a, 0, x1, r) - r^2
eq3 = distance(0, 0, x, y, r, y1) - r^2
eq4 = distance(0, 0, x, y, x1, r) - r^2
eq5 = distance(0, 5, a + b, 0, x1, r) - r^2
eq6 = (5 - y)*(a + b) - 5x
eq7 = x*y + x*(5-y)/2 + (a + b - x)*y/2 - 5(a + b)/2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

nlsove() で解く。

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")

   -r^2 + (y1 - 5*(a^2 - a*r + 5*y1)/(a^2 + 25))^2 + (-a*(a*r - 5*y1 + 25)/(a^2 + 25) + r)^2,  # eq1
   -r^2 + (r - 5*(a^2 - a*x1 + 5*r)/(a^2 + 25))^2 + (-a*(a*x1 - 5*r + 25)/(a^2 + 25) + x1)^2,  # eq2
   -r^2 + (r - x*(r*x + y*y1)/(x^2 + y^2))^2 + (-y*(r*x + y*y1)/(x^2 + y^2) + y1)^2,  # eq3
   -r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq4
   -r^2 + (r - 5*(a^2 + 2*a*b - a*x1 + b^2 - b*x1 + 5*r)/(a^2 + 2*a*b + b^2 + 25))^2 + (x1 - (a + b)*(a*x1 + b*x1 - 5*r + 25)/(a^2 + 2*a*b + b^2 + 25))^2,  # eq5
   -5*x + (5 - y)*(a + b),  # eq6
   -5*a/2 - 5*b/2 + x*y + x*(5 - y)/2 + y*(a + b - x)/2,  # eq7

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, r, x1, y1, x, y) = u
   return [
       -r^2 + (y1 - 5*(a^2 - a*r + 5*y1)/(a^2 + 25))^2 + (-a*(a*r - 5*y1 + 25)/(a^2 + 25) + r)^2,  # eq1
       -r^2 + (r - 5*(a^2 - a*x1 + 5*r)/(a^2 + 25))^2 + (-a*(a*x1 - 5*r + 25)/(a^2 + 25) + x1)^2,  # eq2
       -r^2 + (r - x*(r*x + y*y1)/(x^2 + y^2))^2 + (-y*(r*x + y*y1)/(x^2 + y^2) + y1)^2,  # eq3
       -r^2 + (r - y*(r*y + x*x1)/(x^2 + y^2))^2 + (-x*(r*y + x*x1)/(x^2 + y^2) + x1)^2,  # eq4
       -r^2 + (r - 5*(a^2 + 2*a*b - a*x1 + b^2 - b*x1 + 5*r)/(a^2 + 2*a*b + b^2 + 25))^2 + (x1 - (a + b)*(a*x1 + b*x1 - 5*r + 25)/(a^2 + 2*a*b + b^2 + 25))^2,  # eq5
       -5*x + (5 - y)*(a + b),  # eq6
       -5*a/2 - 5*b/2 + x*y + x*(5 - y)/2 + y*(a + b - x)/2,  # eq7
      ]
end;

iniv = [big"2.8", 2.9, 0.9, 3.3, 1.6, 3.4, 2.1]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[2.800221924622118970432095313760788507835397409564486177867680366139659180899345, 2.91487027133284102199601750983229645815568299593371469024204485373220221457499, 0.8974814920915777221723012652934296052881492161222889500985443442946459740211824, 3.326236908516905783403296227571694673430896734617166839237861420836889759024065, 1.560761971632353997076535469661088458153271795581832971123527071275699538470573, 3.431977890341429946507935106849530222260823283947106463592744755845496558390283, 1.997443109692506364316550909933324776725666683744229847556829344314537003376505], true)

   a = 2.800222;  b = 2.914870;  r = 0.897481;  x1 = 3.326237;  y1 = 1.560762;  x = 3.431978;  y = 1.997443
   2r = 1.794963

等円の直径は 1.794963 寸である。

算額の術では 1.830127 寸としているが，図に描く限り上の方が正しそう。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r, x1, y1, x, y) = res[1]
   @printf("a = %.6f;  b = %.6f;  r = %.6f;  x1 = %.6f;  y1 = %.6f;  x = %.6f;  y = %.6f\n", a, b, r, x1, y1, x, y)
   @printf("2r = %.6f\n", 2r)
   plot([0, a + b, 0, 0, x], [0, 0, 5, 0, y], color=:black, lw=0.5)
   segment(0, 5, a, 0)
   circle(r, y1, r, :blue)
   circle(x1, r, r, :blue)
   if more
       point(x, y, " (x,y)", :green, :left, :bottom)
       point(x1, r, " (x1,r)")
       point(r, y1, " (r,y1)")
       point(a, 0, "a ", :green, :right, :bottom)
       point(a + b, 0, " a+b", :green, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その278）

算額（その278）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額(239)
長野県松本市筑摩 筑摩神社 明治6年(1873)

等脚台形に2本の斜線を引き，分割された区画に2個の等円が入っている。上底，下底がそれぞれ2寸4分，1尺2寸のとき，等円の径をもとめよ。

小数点を避けるため，長さを「分」単位で表す。すなわち，上底，下底をそれぞれ 24（分），120（分）とし，y軸を対称軸とする。
等円の半径を r，台形の高さを y1, 上底と下底の長さを 2a, 2b，斜線と脚の交点座標を(x2, y2) とし，図のように記号を定め，方程式を解く。

include("julia-source.txt");

using SymPy
@syms y1::positive, x2::positive, y2::positive, r::positive;
@syms y1, x2, y2, r, a, b

(a, b) = (24, 120) .// 2

eq1 = (y1 - y2)*(b - a) - y1*(x2 - 12)
eq2 = distance(a, y1, b, 0, 0, y1 - r) - r^2
eq3 = distance(-b, 0, x2, y2, 0, y1 - r) - r^2
eq4 = distance(-b, 0, x2, y2, 0, r) - r^2;

eq2, eq3, eq4 は円の中心から斜線までの平方距離であるが，分数式になっておりそのままでは solve() で解けないので，分子=0 として解く。

@syms d
apart(eq2, d) |> numerator |>  println
apart(eq3, d) |> numerator |>  println
apart(eq4, d) |> numerator |>  println

   -r^2*y1^2 + 1152*r*y1 + 144*y1^2
   -r^2*y2^2 - 2*r*x2^2*y1 - 240*r*x2*y1 + 120*r*x2*y2 - 7200*r*y1 + 7200*r*y2 + x2^2*y1^2 + 120*x2*y1^2 - 120*x2*y1*y2 + 3600*y1^2 - 7200*y1*y2 + 3600*y2^2
   -r^2*y2^2 - 120*r*x2*y2 - 7200*r*y2 + 3600*y2^2

eq1 = (y1 - y2)*(b - a) - y1*(x2 - 12)
eq2 = -r^2*y1^2 + 1152*r*y1 + 144*y1^2
eq3 = -r^2*y2^2 - 2*r*x2^2*y1 - 240*r*x2*y1 + 120*r*x2*y2 - 7200*r*y1 + 7200*r*y2 + x2^2*y1^2 + 120*x2*y1^2 - 120*x2*y1*y2 + 3600*y1^2 - 7200*y1*y2 + 3600*y2^2
eq4 = -r^2*y2^2 - 120*r*x2*y2 - 7200*r*y2 + 3600*y2^2
solve([eq1, eq2, eq3, eq4], (y1, x2, y2, r))

   9-element Vector{NTuple{4, Sym}}:
    (-90, 180/7, -450/7, -20)
    (-20, -60, -50, -60)
    (0, x2, 0, r)
    (20, -60, 50, 60)
    (90, 180/7, 450/7, 20)
    (-24*sqrt(5), 0, -30*sqrt(5), -12*sqrt(5))
    (-24*sqrt(5), 60, 0, -12*sqrt(5))
    (24*sqrt(5), 0, 30*sqrt(5), 12*sqrt(5))
    (24*sqrt(5), 60, 0, 12*sqrt(5))

5 番目の組が適解である。
すなわち等円の半径が 20 であり，元の単位での直径は 20/10*2 = 4 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b) = (60, 12)
   (y1, x2, y2, r) = (90, 180/7, 450/7, 20)
   @printf("y1 = %.6f;  x2 = %.6f;  y2 = %.6f;  r = %.6f\n", y1, x2, y2, r)
   plot([a, b, -b, -a, a], [0, y1, y1, 0, 0], color=:black, lw=0.5)
   segment(-a, 0, x2, y2)
   segment(a, 0, -x2, y2)
   circle(0, y1 - r, r, :blue)
   circle(0, r, r, :blue)
   if more
       point(0, y1 - r, " y1-r")
       point(0, r, " r")
       point(x2, y2, " (x2,y2)")
       point(0, y1, " y1", :green, :left, :top)
       point(a, 0, " a", :green, :left, :bottom)
       point(b, y1, " (b, y1)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;