算額（その341）

算額（その341）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf

7月 高木神社
長方形の中に，春円 1 個，夏円 2 個，秋円 1 個，冬円 1 個を入れる。5 個の円の直径の和が 8 寸のとき，長方形の縦の長さ（春円の直径）を求めよ。

春円の半径，中心座標を r1, (a - r1, r1)
夏円の半径，中心座標を r2, (x2, r2), (r2, r2)
秋円の半径，中心座標を r3, (x3, r3)
冬円の半径，中心座標を r4, (x3, 2r1 - r4)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, r2::positive, x2::positive, r3::positive, x3::positive, r4::positive

eq1 = r1 + 2r2 + r3 + r4 - 8//2
eq2 = (a - r1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand
eq3 = (a - r1 - x3)^2 + (r4 - r1)^2 - (r1 + r4)^2 |> expand
eq4 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand
eq5 = (x2 - x3)^2 + (2r1 - r4 - r2)^2 - (r2 + r4)^2 |> expand
eq6 = (x3 - r2)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand
eq7 = r3 + r4 - r1;

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (a, r1, r2, x2, r3, x3, r4))

   2-element Vector{NTuple{7, Sym}}:
    (2 - sqrt(3), 3/2, 1/2, 1/2, 3/8, 1/2 + sqrt(3)/2, 9/8)
    (2 + 2*sqrt(3), 3/2, 1/2, 1/2 + sqrt(3), 3/8, 1/2 + sqrt(3)/2, 9/8)

   a = 5.4641;  r1 = 1.5;  r2 = 0.5;  x2 = 2.23205;  r3 = 0.375;  x3 = 1.36603;  r4 = 1.125
   春円の直径 = 3;  夏円の直径 = 1;  秋円の直径 = 0.75;  冬の直円径 = 2.25
   長方形の縦の長さ = 春円の直径 = 3

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1, r2, x2, r3, x3, r4) = res[2]
   @printf("a = %g;  r1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g\n", a, r1, r2, x2, r3, x3, r4)
   @printf("春円の直径 = %g;  夏円の直径 = %g;  秋円の直径 = %g;  冬の直円径 = %g\n", 2r1, 2r2, 2r3, 2r4)
   @printf("長方形の縦の長さ = 春円の直径 = %g\n", 2r1)
   plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   circle(a - r1, r1, r1)
   circle(x2, r2, r2, :green)
   circle(r2, r2, r2, :green)
   circle(x3, r3, r3, :blue)
   circle(x3, 2r1 - r4, r4, :orange)
   if more
       point(a - r1, r1, " 春円:r1,(a-r1,r1)", :red, :center)
       point(x2, r2, " 夏円:r2\n(x2,r2)", :green, :center)
       point(r2, r2, " 夏円:r2\n(r2,r2)", :green, :center)
       point(x3, r3, " 秋円:r3\n(x3,r3)", :blue, :center)
       point(x3, 2r1 - r4, " 冬円:r4,(x3,2r1-r4)", :orange, :center)
       point(a, 0, "a ", :black, :right, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

 

算額（その340）

算額（その340）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf

5月 愛宕神社

求めるものは違うが，算額（その142）「三重県四日市市 神明神社 寛政2年」と基本的に同じ。

乾円，坤円の直径がそれぞれ 4 寸，1 寸のとき，地円の直径を求めよ。

なお，乾円の直径と坤円の直径がある値より大きくなれば，鉤が股より大きくなる。算額のような，乾円径 = 4, 坤円径 = 1 の場合だと，算額の図とは異なるイメージのものになる。

上図のように記号を定め，式を立て，解く。solve() ではなく，nlsolve() による。

連立方程式の段階では r4，r5 は記号として扱う。

include("julia-source.txt");

using SymPy

@syms r1, r2, r3, x2, x3, x4, x5, x, y
@syms r4, r5

eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2 |> expand  # 天地
eq2 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2 |> expand  # 地人
eq3 = (x4 - r1)^2 + (r1 - r4)^2 - (r1 + r4)^2 |> expand  # 天乾
eq4 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2 |> expand  # 地乾
eq5 = (x5 - x2)^2 + (r2 - r5)^2 - (r2 + r5)^2 |> expand  # 地坤
eq6 = (x3 - x5)^2 + (r3 - r5)^2 - (r3 + r5)^2 |> expand  # 人坤
eq7 = r1*(x + y + sqrt(x^2 + y^2)) - x*y |> expand  # 直角三角形の面積
eq8 = r1*(x - x2) - r2*(x - r1) |> expand  # 半径の関係（三角形の相似）
eq9 = r1*(x - x3) - r3*(x - r1) |> expand; # 半径の関係（三角形の相似）

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

   r1^2 - 4*r1*r2 - 2*r1*x2 + x2^2, # eq1
   -4*r2*r3 + x2^2 - 2*x2*x3 + x3^2, # eq2
   r1^2 - 4*r1*r4 - 2*r1*x4 + x4^2, # eq3
   -4*r2*r4 + x2^2 - 2*x2*x4 + x4^2, # eq4
   -4*r2*r5 + x2^2 - 2*x2*x5 + x5^2, # eq5
   -4*r3*r5 + x3^2 - 2*x3*x5 + x5^2, # eq6
   r1*x + r1*y + r1*sqrt(x^2 + y^2) - x*y, # eq7
   r1*r2 + r1*x - r1*x2 - r2*x, # eq8
   r1*r3 + r1*x - r1*x3 - r3*x, # eq9

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)
   (r1, r2, r3, x2, x3, x4, x5, x, y) = u
   return [
       r1^2 - 4*r1*r2 - 2*r1*x2 + x2^2, # eq1
       -4*r2*r3 + x2^2 - 2*x2*x3 + x3^2, # eq2
       r1^2 - 4*r1*r4 - 2*r1*x4 + x4^2, # eq3
       -4*r2*r4 + x2^2 - 2*x2*x4 + x4^2, # eq4
       -4*r2*r5 + x2^2 - 2*x2*x5 + x5^2, # eq5
       -4*r3*r5 + x3^2 - 2*x3*x5 + x5^2, # eq6
       r1*x + r1*y + r1*sqrt(x^2 + y^2) - x*y, # eq7
       r1*r2 + r1*x - r1*x2 - r2*x, # eq8
       r1*r3 + r1*x - r1*x3 - r3*x, # eq9
   ]
end;

(r4, r5) = (4, 1) .// 2  # r4, r5 に数値を代入
iniv = [big"40.0", 25, 18, 101, 145, 72, 125, 230, 95]
res = nls(H, ini=iniv)
println(res);
   (BigFloat[17.99999999999999999999959359538046894328544006940578449319282835040346165571312, 4.500000000000000000000021850969763478122224875168800267359829481564552397400594, 1.125000000000000000000015511954584626415427036519923717249791558076525292630489, 35.99999999999999999999950081078087049912383423570819427141543308984033365075256, 40.49999999999999999999954859392954415808055400655154524243081491576853407899673, 29.99999999999999999999948510964983470350422331631829474897253215082231657323331, 38.99999999999999999999950866369699465815370090785716009867301337616840188238613, 41.99999999999999999999957349741006947409069773160858398799318001605314871943797, 143.9999999999999993599636251260970577496755779165141698135744945781538654670987], true)

   鉤 = 144;  股 = 42
   天円径 = 36;  地円径 = 9;  人円径 = 2.25
   乾円径 = 4;  坤円径 = 1

地円の半径 = 4.5。元の単位では，直径 9 寸。

なお，術では「乾円径と坤円径の積の平方根を2倍したものに乾円径と坤円径を加える」とのこと。
2sqrt(4 * 1) + 4 + 1 = 9

2sqrt(4 * 1) + 4 + 1

   9.0

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x2, x3, x4, x5, x, y) = res[1]
   @printf("鉤 = %g;  股 = %g\n", y, x)
   @printf("天円径 = %g;  地円径 = %g;  人円径 = %g\n", 2r1, 2r2, 2r3)
   @printf("乾円径 = %g;  坤円径 = %g\n", 2r4, 2r5)
   plot([0, x, 0, 0], [0, 0, y, 0], color=:black, lw=0.5)
   circle(r1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :magenta)
   circle(x4, r4, r4, :orange)
   circle(x5, r5, r5, :green)
   if more == true
       point(r1, r1, "天:(r1,r1)", :red, :center, :bottom)
       point(x2, r2, "地:(x2,r2)", :blue, :center, :bottom)
       point(x3, r3, "人:(x3,r3)", :magenta, :center, :bottom)
       point(x4, r4, "乾:(x4,r4)", :orange, :center, :bottom)
       point(x5, r5, "坤:(x5,r5)", :green, :center, :bottom)
       # point(0, y, " y", :black, :left, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-2, 45), ylims=(-2, 20))
   else
       plot!(showaxis=false)
   end
end;

算額（その339）

算額（その339）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf

3月 田村大元神社
外円の中に 2 個の等円があり，その交わった部分に共通弦を設ける。
外円の直径が 9 寸，弦の長さが 3 寸のとき，等円の直径を求めよ。

外円と等円の半径を r，r1 とし，以下の方程式を解く。

include("julia-source.txt");

using SymPy

@syms r::positive, r1::positive, 弦::positive;

eq = (r - r1)^2 + (弦//2)^2 - r1^2
r1 = solve(eq, r1)[1]
2r1 |> println

   r + 弦^2/(4*r)

等円の直径は，「外円の半径に弦の二乗を外円の直径の2倍で割ったものを足す」

外円の直径が 9寸，弦の長さが 3寸のとき，等円径は 5寸 である。

r = 9//2
弦 = 3
等円径 = r + 弦^2/(4*r)
@printf("等円径 = %g\n", 等円径)

   等円径 = 5

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r1) = (9//2, 5//2)
   plot()
   circle(0, 0, r, :blue)
   circle(0, r - r1, r1, :red)
   circle(0, r1 - r, r1, :red)
   segment(-3//2, 0, 3//2, 0, :green, lw=2)
   if more
       point(3/2, 0, "  3/2", :magenta, :left, :bottom)
       point(r, 0, "r0=9/2 ", :blue, :right, :bottom)
       point(0, r - r1, " r-r1", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その338）

算額（その338）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
2月 諏訪神社

一辺の長さが 9寸の正方形の中に図のように直角三角形を入れる。中央に小さな正方形ができるが，この面積が直角三角形1個の面積と等しくなるようにするとき(すなわち外側の正方形を 5 等分する），直角三角形の鉤，股（それぞれ直角三角形の直角を挟む短い方の辺と長い方の辺）の長さを求めよ。

外側の正方形の一辺の長さを a とし，正方形の中心を原点とする。第1象限内にある頂点 A の座標を (x, y) とする。鉤は AB, 股は AC である。以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive, y::positive, AB::positive, AC::positive, a::positive

eq1 = (a//2 - x)^2 + (a//2 - y)^2 - AB^2
eq2 = (a//2 + x)^2 + (a//2 - y)^2 - AC^2
eq3 = AB*AC//2 - a^2//5
eq4 = AB^2 + AC^2 - a^2

solve([eq1, eq2, eq3, eq4], (x, y, AB, AC))

   2-element Vector{NTuple{4, Sym}}:
    (3*a/10, a/10, sqrt(5)*a/5, 2*sqrt(5)*a/5)
    (3*a/10, 9*a/10, sqrt(5)*a/5, 2*sqrt(5)*a/5)

最初の方が適解。
x = 3*a/10，y = a/10，鉤 = a/√5，股 = 2a/√5 である。

a = 9寸 のとき，鉤 = 4.024922359499621寸，股 = 8.049844718999243寸である。

a = 9
println("鉤 = $(a/√5)")
println("股 = $(2a/√5)")

   鉤 = 4.024922359499621
   股 = 8.049844718999243

直角三角形の面積は a/√5 * 2a/√5 / 2 = 2a^2/5 = 16.2, 正方形の面積の 1/5 は 81/5 = 16.2

a/√5 * 2a/√5 / 2, 2a^2/5 / 2, a^2/5

   (16.2, 16.2, 16.2)

図の全体を描くには，⊿ABC を90°回転して4個描けばよい。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 9
   (x, y) = (3*a/10, a/10)
   plot()
   # rect(-a/2, -a/2, a/2, a/2, :black)
   A = [0.0 -1.0; 1.0 0]
   xy = [[x, a/2, -a/2, x], [y, a/2, a/2, y]]
   for i = 1:4
       plot!(xy[1], xy[2], color=:blue, lw=1)
       xy = A * xy
   end
   if more
       point(x, y, " A:(x,y)\n(3*a/10,a/10)")
       point(a/2, a/2, "B:(a/2,a/2) ", :green, :right, :bottom)
       point(-a/2, a/2, " C:(-a/2,a/2)", :green, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その337）

算額（その337）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
1月 龍穏院

長方形に三角形を内接させると，図のように甲，乙，丙の直角三角形ができる。長方形の面積が 16歩，甲と乙の面積がそれぞれ 4歩，2歩のとき，三角形の面積を求めよ。

中学レベルだと長方形の面積（直）から甲，乙，丙の3個の直角三角形の面積を引いた残りが求める三角形の面積とするであろう。

図のように，横 = a，縦 = b の長方形があり，c, d が横，縦を区分するとき，以下のようになる。

include("julia-source.txt");

using SymPy

@syms a, b, c, d, 甲, 乙, 丙, 直, 三角

直 = a*b
甲 = a*d/2
乙 = b*c/2
丙 = (a - c)*(b - d)/2;
三角 = 直 - 甲 - 乙 - 丙 |> expand;

これを SymPy で簡約化すると，a*b/2 - c*d/2 になる。見慣れない式になるが，図のように等積変形を行うと，このような式になる。

三角 = 直 - 甲 - 乙 - 丙
⊿ODF = □OAEB - ⊿OAD　- ⊿OBF - ⊿DEF
     = □OAEB - ⊿OAD　- ⊿OBX - ⊿DEX
     = ⊿BEX
     = a*(b - x)/2
     = a*(b - c*d/a)/2
     = (a*b - a*c*d/a)/2
     = a*b/2 - c*d/2
     = a*b/2 - (2乙/b)*(2甲/a)/2
     = 直/2 - 2甲*乙/直

三角 |> println

   a*b/2 - c*d/2

更に，甲，乙の面積の計算式から c, d を代入すると，以下のようになる。

直/2 - (2甲/a)*(2乙/b)/2 |> expand |> println

   a*b/2 - c*d/2

整理すると算額の術に書かれている式になる。
直/2 - 2甲*乙/直
直，甲，乙だけを用いるトリッキーな式であるが，実際には適用しにくい（甲，乙，丙のどれを使ってもよいのでないのは明らかだが）。最初に述べた中学レベルの解法のほうが間違いもない。

直/2 - 2甲*乙/直 |> expand |> println

   a*b/2 - c*d/2

s(直, 甲, 乙) = 直/2 - 2甲*乙/直;  # 関数定義

s(16, 4, 2)

   7.0

using Plots

function polygon(xys, color=:blue; lw=0.5, alpha=0.2, fill=false)
   (x, y) = ([], [])
   for xy in xys
       append!(x, xy[1])
       append!(y, xy[2])
   end
   append!(x, x[1])
   append!(y, y[1])
   if fill
       plot!(x, y, linecolor=color, lw=lw, seriestype=:shape, fillcolor=color, alpha=alpha) 
   else
       plot!(x, y, color=color, lw=lw) 
   end
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, d) = (7, 5, 2, 3)
   #(a, b, c, d) = (8, 2, 2, 1)
   x = d/a * c
   plot()
   rect(0, 0, a, b, :black, fill=false)
   polygon([(0, 0), (a, 0), (a, d)], :cadetblue4, alpha=0.8, fill=true)  # 甲
   polygon([(0, 0), (c, b), (0, b)], :chocolate1, alpha=0.8, fill=true)  # 乙
   polygon([(0, 0), (c, x), (0, b)], :chocolate1, alpha=0.8, fill=true)  # 乙
   polygon([(a, d), (a, b), (c, b)], :coral3, alpha=0.8, fill=true)  # 丙
   polygon([(a, d), (a, b), (c, x)], :coral3, alpha=0.8, fill=true)  # 丙
   polygon([(0, 0), (c, b), (a, d)], :blue4, alpha=0.3, fill=false)  # 三角
   segment(c, 0, c, b, :red, linestyle=:dash)
   if more
       annotate!([
               (0, 0, ("O ", 10, :right, :top)),
               (a, 0, (" A:a", 10, :left, :top)),
               (0, b, ("B:b ", 10, :right, :bottom)),
               (a, b, (" E:(a,b)", 10, :bottom)),
               (a, d, (" D:d", 10, :left)),
               (c, 0, (" C:c", 10, :left, :top)),
               (c, b, ("F", 10, :bottom)),
               (c, x, (" X:(c,x), x=c*d/a", 10, :white, :left, :top)),
               (4a/5, d/2, ("甲")),
               (0.5, b/2, ("乙")),
               (a - 1, (b+d)/2, ("丙"))
               ])
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       plot!(xlims=(-0.8, a + 0.5), ylims=(-0.5, b + 0.5))
   else
       plot!(showaxis=false)
   end
end;

 

算額（その336）

算額（その336）

愛知県名古屋市熱田区 熱田神宮 文化3年(1806)
http://www.wasan.jp/aichi/atuta1.html

二等辺三角形の中に直交する斜線甲と斜線乙があり，面積の等しい円が 3 個内接している。斜線乙が 2 寸のとき，円の直径はいくらか。

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

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, d::positive,
     h::positive, w::positive, r::positive,
     x1::negative, y1::positive, x2::positive, y2::positive,
     X1::negative, X2::positive, Y2::positive, Y3::positive,
     dummy;

a = sqrt((x1 + w)^2 + y1^2)
b = sqrt((x2 - x1)^2 + (y2 - y1)^2)
c = sqrt((w - x2)^2 + y2^2)
d = sqrt(h^2 + w^2)
eq1 = a + 2 - 2w - 2r
eq2 = b + 2 - c - 2r
eq3 = 2a - (2w + a + 2)r
eq4 = 2b - (b + c + 2)r
eq5 = distance(0, h, w, 0, 0, Y3) - r^2
eq6 = distance(0, h, w, 0, X2, Y2) - r^2
eq7 = distance(-w, 0, x2, y2, X2, Y2) - r^2
eq8 = distance(-w, 0, x2, y2, 0, Y3) - r^2
eq9 = distance(-w, 0, x2, y2, X1, r) - r^2
eq10 = distance(x1, y1, w, 0, X1, r) - r^2
eq11 = distance(x1, y1, w, 0, X2, Y2) - r^2;

println(eq1, ",  # eq1")
println(eq2, ",  # eq2")
println(eq3, ",  # eq3")
println(eq4, ",  # eq4")
println(eq5, ",  # eq5")
println(eq6, ",  # eq6")
println(eq7, ",  # eq7")
println(eq8, ",  # eq8")
println(eq9, ",  # eq9")
println(eq10, ",  # eq10")
println(eq11, ",  # eq11")

   -2*r - 2*w + sqrt(y1^2 + (w + x1)^2) + 2,  # eq1
   -2*r - sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2,  # eq2
   -r*(2*w + sqrt(y1^2 + (w + x1)^2) + 2) + 2*sqrt(y1^2 + (w + x1)^2),  # eq3
   -r*(sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2) + 2*sqrt((-x1 + x2)^2 + (-y1 + y2)^2),  # eq4
   h^2*w^2*(-Y3 + h)^2/(h^2 + w^2)^2 - r^2 + (Y3 - h*(Y3*h + w^2)/(h^2 + w^2))^2,  # eq5
   -r^2 + (X2 - w*(X2*w - Y2*h + h^2)/(h^2 + w^2))^2 + (Y2 - h*(-X2*w + Y2*h + w^2)/(h^2 + w^2))^2,  # eq6
   -r^2 + (X2 - (X2*(w^2 + 2*w*x2 + x2^2 + y2^2) - y2*(X2*y2 - Y2*w - Y2*x2 + w*y2))/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (Y2 - y2*(X2*w + X2*x2 + Y2*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq7
   -r^2 + y2^2*(Y3*w + Y3*x2 - w*y2)^2/(w^2 + 2*w*x2 + x2^2 + y2^2)^2 + (Y3 - y2*(Y3*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq8
   -r^2 + (X1 - (X1*(w^2 + 2*w*x2 + x2^2 + y2^2) - y2*(X1*y2 - r*w - r*x2 + w*y2))/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (r - y2*(X1*w + X1*x2 + r*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq9
   -r^2 + (X1 - (X1*w^2 - 2*X1*w*x1 + X1*x1^2 - r*w*y1 + r*x1*y1 + w*y1^2)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (r - y1*(-X1*w + X1*x1 + r*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq10
   -r^2 + (X2 - (X2*(w^2 - 2*w*x1 + x1^2 + y1^2) - y1*(X2*y1 + Y2*w - Y2*x1 - w*y1))/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (Y2 - y1*(-X2*w + X2*x1 + Y2*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq11

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)
   (h, w, r, x1, y1, x2, y2, X1, X2, Y2, Y3) = u
   return [
       -2*r - 2*w + sqrt(y1^2 + (w + x1)^2) + 2,  # eq1
       -2*r - sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2,  # eq2
       -r*(2*w + sqrt(y1^2 + (w + x1)^2) + 2) + 2*sqrt(y1^2 + (w + x1)^2),  # eq3
       -r*(sqrt(y2^2 + (w - x2)^2) + sqrt((-x1 + x2)^2 + (-y1 + y2)^2) + 2) + 2*sqrt((-x1 + x2)^2 + (-y1 + y2)^2),  # eq4
       h^2*w^2*(-Y3 + h)^2/(h^2 + w^2)^2 - r^2 + (Y3 - h*(Y3*h + w^2)/(h^2 + w^2))^2,  # eq5
       -r^2 + (X2 - w*(X2*w - Y2*h + h^2)/(h^2 + w^2))^2 + (Y2 - h*(-X2*w + Y2*h + w^2)/(h^2 + w^2))^2,  # eq6
       -r^2 + (X2 - (X2*w^2 + 2*X2*w*x2 + X2*x2^2 + Y2*w*y2 + Y2*x2*y2 - w*y2^2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (Y2 - y2*(X2*w + X2*x2 + Y2*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq7
       -r^2 + y2^2*(Y3*w + Y3*x2 - w*y2)^2/(w^2 + 2*w*x2 + x2^2 + y2^2)^2 + (Y3 - y2*(Y3*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq8
       -r^2 + (X1 - (X1*(w^2 + 2*w*x2 + x2^2 + y2^2) - y2*(X1*y2 - r*w - r*x2 + w*y2))/(w^2 + 2*w*x2 + x2^2 + y2^2))^2 + (r - y2*(X1*w + X1*x2 + r*y2 + w^2 + w*x2)/(w^2 + 2*w*x2 + x2^2 + y2^2))^2,  # eq9
       -r^2 + (X1 - (X1*(w^2 - 2*w*x1 + x1^2 + y1^2) - y1*(X1*y1 + r*w - r*x1 - w*y1))/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (r - y1*(-X1*w + X1*x1 + r*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq10
       -r^2 + (X2 - (X2*w^2 - 2*X2*w*x1 + X2*x1^2 - Y2*w*y1 + Y2*x1*y1 + w*y1^2)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2 + (Y2 - y1*(-X2*w + X2*x1 + Y2*y1 + w^2 - w*x1)/(w^2 - 2*w*x1 + x1^2 + y1^2))^2,  # eq11
   ]
end;

iniv = [big"4.3", 1.3, 0.4, -0.3, 0.9, 0.5, 2.2, -0.2, 0.4, 1.2, 2.4]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[4.285714285714285714285714285714285714285714285714285714285574246330323963257686, 1.249999999999999999999999999999999999999999999999999999999958894684739139553026, 0.4999999999999999999999999999999999999999999999999999999999816989607612493273477, -0.3499999273461119877023840155306407964015055211277256589619837238743919315567991, 1.199999945509583990776788011647980597301129140845794244221458224297416065202069, 0.5499999999999999999999999999999999999999999999999999999999544520988613601022169, 2.399999999999999999999999999999999999999999999999999999999895684919240762427152, -0.2499999999999999999999999999999999999999999999999999999999983023984933316200398, 0.3499999999999999999999999999999999999999999999999999999999877259401613410750471, 1.299999999999999999999999999999999999999999999999999999999966810926003067216378, 2.499999999999999999999999999999999999999999999999999999999924978510268800929172], true)

names = ["h", "w", "r", "x1", "y1", "x2", "y2", "X1", "X2", "Y2", "Y3"]
for (i,  name) in enumerate(names)
   @printf("%2s = %.6f\n", name, res[1][i])
end

    h = 4.285714
    w = 1.250000
    r = 0.500000
   x1 = -0.350000
   y1 = 1.200000
   x2 = 0.550000
   y2 = 2.400000
   X1 = -0.250000
   X2 = 0.350000
   Y2 = 1.300000
   Y3 = 2.500000

円の直径は 2r = 1 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (h, w, r, x1, y1, x2, y2, X1, X2, Y2, Y3) = res[1]
   plot([w, 0, -w, w], [0, h, 0, 0], color=:black, lw=0.5)
   circle(X1, r, r)
   circle(X2, Y2, r)
   circle(0, Y3, r)
   segment(-w, 0, x2, y2)
   segment(x1, y1, w, 0)
   if more
       point(0, Y3, " Y3")
       point(X2, Y2, "(X2,Y2)\n", :green, :center, :bottom)
       point(X1, r, "(X1,r)\n", :green, :center, :bottom)
       point(-w, 0, "-w")
       point(w, 0, "w")
       point(0, h, " h")
       point(x1, y1, "(x1,y1)", :green, :right, :bottom)
       point(x2, y2, " (x2,y2)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その335）

算額（その335）

石川県能登町 白山神社
【算額に挑戦】石川県の算額 －その１－
http://blog.livedoor.jp/enjoy_math/archives/51203926.html

等脚台形に甲円，乙円が入っている。それぞれは台形の辺に3点で内接し，また界斜にも内接している。
上頭は 12寸，甲円，乙円の直径はそれぞれ 24寸，15寸である。界斜の長さはいかほどか。

元の図を反時計回りに 90° 回転したほうが若干わかりやすいか。
甲円，乙円の半径を r2, r1 とする。等脚台形の脚を延長し，交点を原点 O とする。∠AOC を θとする。AC = 上頭/2 = 6 である。
AO = a, OC = b として，a, b を求める。

include("julia-source.txt");

using SymPy

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

(r1, r2, c) = (15//2, 24//2, 6)
sinθ = c/b
eq1 = (a + r1)sinθ - r1
eq2 = a^2 + c^2 - b^2

solve([eq1, eq2], (a, b))[1]

   (80/3, 82/3)

(a, b) = (80//3, 82//3)
(sinθ, tanθ, cosθ) = (c/b, c/a, a/b)

   (9//41, 9//40, 40//41)

界斜と脚の交点 (x1, y1), (x2, y2) を求める。

@syms x1::positive, y1::positive, x2::positive, y2::negative
y1 = x1*tanθ
y2 = -x2*tanθ
eq3 = distance(x1, y1, x2, y2, a + r1, 0) - r1^2
eq4 = distance(x1, y1, x2, y2, 2b, 0) - r2^2
solve([eq3, eq4])

   2-element Vector{Dict{Any, Any}}:
    Dict(x1 => 5200/123 - 40*sqrt(10)/41, x2 => 40*sqrt(10)/41 + 5200/123)
    Dict(x1 => 40*sqrt(10)/41 + 5200/123, x2 => 5200/123 - 40*sqrt(10)/41)

2 番目のものが適解である（1 番目のものと対称）。

(x1, x2) = (5200//123 + 40*sqrt(Sym(10))/41, 5200//123 - 40*sqrt(Sym(10))/41)
(y1, y2) = (x1, -x2) .* tanθ
(x1, y1, x2, y2) |> println

   (40*sqrt(10)/41 + 5200/123, 9*sqrt(10)/41 + 390/41, 5200/123 - 40*sqrt(10)/41, -390/41 + 9*sqrt(10)/41)

界斜の長さは，2 点間の距離を求めればよい。

界斜の長さ = sqrt((x1 - x2)^2 + (y1 - y2)^2)
@printf("界斜の長さ = %.6g\n", 界斜の長さ)

   界斜の長さ = 20

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot()
   segment(0, 0, (2b + r2), (2b + r2)*tanθ, :green)
   segment(0, 0, (2b + r2), -(2b + r2)*tanθ, :green)
   segment((2b + r2), (2b + r2)*tanθ, (2b + r2), -(2b + r2)*tanθ)
   segment(a, c, a, -c)
   circle(a + r1, 0, r1)
   segment(a + r1, 0, (a + r1)*cosθ^2, (a + r1)*cosθ*sinθ)
   circle(2b, 0, r2)
   segment(2b, 0, 2b*cosθ^2, 2b*cosθ*sinθ)
   segment(x1, y1, x2, y2, :blue)
   if more
       point(a, 0, "a ", :green, :right)
       point(a, c, "(a,c) ", :green, :right, :bottom)
       point(a + r1, 0, " a+r1")
       point(2b, 0, "2b ", :green, :right)
       point(2b + r2, 0, "2b+r2 ", :green, :right)
       point(2b + r2, (2b + r2)*tanθ, "(2b+r2,(2b+r2)*tanθ", :green, :right, :bottom)
       point(x1, y1, "(x1,y1)", :blue, :right, :bottom)
       point(x2, y2, "(x2,y2)", :blue, :right, :top)
       point(0, 0, " O", :red)
       point(a, 0, " A", :red, :left)
       point(a, c, " C", :red, :left, :bottom)
       annotate!(0.5, 15, text("a = 80/3, b = (0,0)-(a,c) = 82/3, c = 6\n" *
               "∠AOC をθとおく\n2b+r2 = 200/3, (2b + r2)*tanθ = 15\n",
               :black, :left, :top, 9))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その334）

算額（その334）

遠藤寛子:「算法少女」
https://www.chikumashobo.co.jp/product/9784480090133/

久留米藩の殿様の有馬頼徸が，千葉あきと中根宇多に出した問題。
『円のうちに，大円二個，小円二個が接した形があるが，それらの大円小円は，またおたがいに接している。いま，いちばん外側の円の直径を七寸，内に接している大きい方の円の直径を三寸としたら，小円の直径はいかほどか』

図は示されていないが以下のようなものを考えた。

include("julia-source.txt");

using SymPy

@syms r::positive, r1::positive, y1::positive, r2::positive, y2::negative;

(r, r1) = (7//2, 3//2)

eq1 = (r1 - r2)^2 + (y1 - y2)^2 - (r1 + r2)^2
eq2 = r1^2 + y1^2 - (r - r1)^2
eq3 = r2^2 + y2^2 - (r - r2)^2

res1 = solve([eq1, eq2, eq3], (r2, y1, y2))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (21*sqrt(14)/169 + 315/338, (7 + 6*sqrt(14))*sqrt(79 - 12*sqrt(14))/130, -7*sqrt(79/676 - 3*sqrt(14)/169))

小円の半径は 1.39689233800150 で，元の単位での直径は 2寸7分9厘3毛7糸...である。
あきの出した答えは「2寸8分」。きれいな数値だ。

res1[1][1].evalf() |> println

   1.39689233800150

1.39689233800150*2

   2.793784676003

y1, y2 は二重根号を外すと簡約化できる。

res1[1][2] |> sympy.sqrtdenest |> simplify |> println

   sqrt(7)/2

res1[1][3] |> sympy.sqrtdenest |> simplify |> println

   -21*sqrt(2)/13 + 7*sqrt(7)/26

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r1) = (7//2, 3//2)
   (r2, y1, y2) = ((42*sqrt(14) + 315)/338, sqrt(7)/2, (7*sqrt(7) - 42*sqrt(2))/26)
   @printf("小円の直径 = %.6f\n", 2r2)
   plot()
   circle(0, 0, r)
   circle(r1, y1, r1, :blue)
   circle(-r1, y1, r1, :blue)
   circle(r2, y2, r2, :orange)
   circle(-r2, y2, r2, :orange)
   if more
       point(r1, y1, "大円:r1\n(r1,y1)", :blue)
       point(r2, y2, "小円:r2\n(r2,y2)", :orange)
       point(r, 0, "r ", :red, :right)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

ネットを散策して，もう一つの図の可能性がわかった。映画版では以下のような図が示されていたとのこと。

街角の数学
http://streetwasan.web.fc2.com/math17.12.13.html

こちらは，大円，小円の中心が x 軸，y 軸上にあるので非常に簡単で eq は 10r2 = 14 になり爆速で答えが出る（お話にならない）。しかもきれいな数値で，あきが答えた通り，直径は 2.8 = 2寸8 分である。
しかし，有馬の殿様がこんな簡単な問題を出したとは思えない。やはり，最初の図のほうが解きがいがあるのでは?

using SymPy
@syms r, r1, r2
(r, r1) = (7//2, 3//2)
eq = (r - r1)^2 + (r - r2)^2 - (r1 + r2)^2 |> expand
eq |> println
solve(eq)[1] |> println

   14 - 10*r2
   7/5

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r1) = (7//2, 3//2)
   r2 = 7/5
   @printf("小円の直径 = %.6f\n", 2r2)
   plot()
   circle(0, 0, r)
   circle(0, r - r1, r1, :blue)
   circle(0, r1 - r, r1, :blue)
   circle(r - r2, 0, r2, :orange)
   circle(r2 - r, 0, r2, :orange)
   if more
       point(0, r - r1, " r-r1", :blue)
       point(r - r2, 0, " r-r2", :orange)
       point(r, 0, "r ", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その333）

算額（その333）

遠藤寛子:「算法少女」
https://www.chikumashobo.co.jp/product/9784480090133/
国立国会図書館:江戸の数学，コラム 算額
https://www.ndl.go.jp/math/s1/c5.html

半円の中に直角三角形と等円 2 個が入っている。半円と等円の直径の比を求めよ。

半円の半径と中心座標を r, (0, 0)
等円の半径と中心座標を r1, (x1, y1), (x2, r1) とおく。
作図に必要なパラメータも含めて一挙に solve() で解くのは無理なので，まずは r，r1 を求める。a は 2 個の等円の共通接線となる (-r, 0) と（x，y）を結ぶ弦の長さである。

include("julia-source.txt");

using SymPy

@syms r::positive, r1::positive, a::positive;

eq1 = a^2 + (2r - 4r1) - (2r)^2
eq2 = (2r - 4r1)a - (a + (2r - 4r1) + 2r)r1
eq3 = a + (2r - 4r1) - 2r - 2r1

res1 = solve([eq1, eq2, eq3], (r, r1, a))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (13/10, 2/5, 12/5)

等円の半径は半円の半径の 4/13 倍である。

これらを既知として，図を描くためのパラメータを求める。

@syms r::positive, r1::positive, a::positive,
     x1::negative, y1::positive, x2::positive,
     x::positive, y::positive;

(r, r1, a) = (13//10, 2//5, 12//5)
eq4 = x1^2 + y1^2 - (r - r1)^2
eq5 = 2r*y - (2r - 4r1)a
eq6 = x^2 + y^2  - r^2
eq7 = distance(-r, 0, x, y, x1, y1) - r1^2
eq8 = distance(-r, 0, x, y, x2, r1) - r1^2;
res2 = solve([eq4, eq5, eq6, eq7, eq8], (x1, y1, x2, x, y))

   1-element Vector{NTuple{5, Sym}}:
    (-9/26, 54/65, 7/10, 119/130, 12/13)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, r1, a) = (13/10, 2/5, 12/5)
   (x1, y1, x2, x, y) = (-9/26, 54/65, 7/10, 119/130, 12/13)
   @printf("半円の半径 = %g;  円の半径 = %g\n",  r, r1)
   plot([-r, r, x, -r], [0, 0, y, 0], color=:black, lw=0.5)
   circle(0, 0, r, beginangle=0, endangle=180)
   circle(x1, y1, r1, :blue)
   circle(x2, r1, r1, :blue)
   if more
       point(r, 0, " r", :red, :left, :bottom)
       point(x1, y1, "等円:r1,(x1,y1)", :blue, :center)
       point(x2, r1, "等円:r1,(x2,r1)", :blue, :center)
       point(x, y, " (x,y)", :black, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その332）

算額（その332）

石川県能登町 白山神社
【算額に挑戦】石川県の算額 ―その2―
http://blog.livedoor.jp/enjoy_math/archives/51206109.html

直線上に甲円と乙円が載っており，互いに接している。
大円と小円も同じ直線上にあり，互いに接している。
また，甲円と大円，乙円と小円も互いに接している。
甲円，乙円，大円の直径がそれぞれ 64寸，54寸，11寸のとき，小円の直径を求めよ。

甲円の直径，中心座標を r1, (0, r1)
乙円の直径，中心座標を r2, (x2, r2)
大円の直径，中心座標を r3, (x3, r3)
小円の直径，中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
なお，同一直線に載っている2つの円 r1, r2 の中心間の水平距離は 2√(r1*r2) なので，eq1 は見かけは長ったらしいが実際には x2^2 = 3456 という式になっている。ほかも同様である。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x2::positive,
     r3::positive, x3::positive, r4::positive, x4::positive;
#@syms r1, r2, x2, r3, x3, r4, x4;
(r1, r2, r3) = (64, 54, 11) .// 2
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq4 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4], (x2, x3, r4, x4))

   1-element Vector{NTuple{4, Sym}}:
    (24*(18 - sqrt(66))/sqrt(65 - 6*sqrt(66)), 8*(-11 + 3*sqrt(66))/sqrt(65 - 6*sqrt(66)), 211232/1849 - 24960*sqrt(66)/1849, 48*sqrt(4290/1849 - 396*sqrt(66)/1849))

それぞれの式は，必要に応じて分母の有理化や二重根号を外すなどで簡約化できる場合もある。

@syms dummy
apart(res[1][1], dummy) |> sympy.sqrtdenest |> simplify |> println
apart(res[1][2], dummy) |> sympy.sqrtdenest |> simplify |> println
apart(res[1][3], dummy) |> println
apart(res[1][4], dummy) |> sympy.sqrtdenest |> simplify |> println

   24*sqrt(6)
   8*sqrt(11)
   211232/1849 - 24960*sqrt(66)/1849
   -528*sqrt(6)/43 + 864*sqrt(11)/43

小円の直径は 2(211232 - 24960*sqrt(66))/1849 = 9.146567247470445 である。

2(211232 - 24960*sqrt(66))/1849

   9.146567247470445

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (64, 54, 11) .// 2
   x2 = 24*sqrt(6)
   x3 = 8*sqrt(11)
   r4 = 211232/1849 - 24960*sqrt(66)/1849
   x4 = -528*sqrt(6)/43 + 864*sqrt(11)/43
   @printf("小円の直径 = %.6f\n", 2r4)
   plot()
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   circle(x4, r4, r4, :orange)
   if more
       point(0, r1, " 甲円\n r1", :red, :left, :vcenter)
       point(x2, r2, " 乙円:r2\n (x2,r2)", :blue, :left, :vcenter)
       point(x3, r3, " 大円:r3   \n(x3,r3)   ", :green, :right, :bottom)
       point(x4, r4, "    小円:r4,(x4,r4)", :orange, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その331）

算額（その331）

静岡県磐田市蒲田 医王寺 安永8年(1779)
http://www.wasan.jp/sizuoka/ioji2.html
数学者による和算額の内容解析
https://www.iouji.net/wp-content/uploads/sugaku.pdf

長方形 ABCD があり，線分 BE の二乗と 平(AB) の和が 237 寸，AF は 13寸6分，FC は 3寸4分である。長方形の二辺の長さ（平と長）はいかほどか。

AB，BC，CE の長さを a, c, e とおき，直線の交点 F の座標を (x, y) とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, c::positive, e::positive,
     x::positive, y::positive;
eq1 = (c^2 + e^2) + a - 237
eq2 = (a - y)^2 + x^2 - 13.6^2
eq3 = (c - x)^2 + y^2 - 3.4^2
eq4 = y/e - 13.6/17
eq5 = a^2 + c^2 - 17^2;
res = solve([eq1, eq2, eq3, eq4, eq5], (a, c, e, x, y))

   1-element Vector{NTuple{5, Sym}}:
    (8.00000000000000, 15.0000000000000, 2.00000000000000, 12.0000000000000, 1.60000000000000)

長方形の縦，横の長さは 8寸，15寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, c, e, x, y) = (8, 15, 2, 12, 1.6)
   @printf("平 = %d;  長 = %d;  e = %d;  x = %d;  y = %.1f\n", a, c, e, x, y)
   plot([0, c, c, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   segment(0, a, c, 0, :blue)
   segment(0, 0, c, e, :red)
   if more
       point(x, y, "(x,y)\n", :green, :center, :bottom)
       point(0, a, " a", :black, :left, :bottom)
       point(c, 0, " c", :black, :left, :bottom)
       point(c, e, " e", :black, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その330）

算額（その330）

静岡県磐田市蒲田 医王寺 安永8年(1779)
http://www.wasan.jp/sizuoka/ioji2.html
【算額に挑戦】静岡県の算額 ―その1―
http://blog.livedoor.jp/enjoy_math/archives/cat_50038397.html
数学者による和算額の内容解析
https://www.iouji.net/wp-content/uploads/sugaku.pdf

鈎股弦（直角三角形）内に界弦（斜線;図の青線）を引き，二分割された各領域それぞれに直径が同じ円を入れる。円は鈎，股，弦（全弦），界弦に接している。弦は界弦より 7 寸長く，鈎は 8 寸である。このとき，界弦，股，弦，円の直径を求めよ。

界弦と股の交点の x 座標を a，右側の円の中心の x 座標を x とおく。以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鉤::positive, 股::positive, 全弦::positive,
     界弦::positive, r::positive, a::positive,
     x::positive, dummy;

鉤 = 8
eq1 = 全弦 - 界弦 - 7
eq2 = 鉤^2 + a^2 - 界弦^2 
eq3 = 鉤 + a - 界弦 - 2r
eq4 = 鉤^2 + 股^2 - 全弦^2
eq6 = numerator(apart(distance(0, 鉤, a, 0, x, r) - r^2, dummy))
eq7 = (界弦 + (股 - a) + 全弦)r - (股 - a)*鉤;

解を求める変数の指定順によっては（一定の時間内に）解が求まらないこともあるので注意。

二番目の組のものが適解である。

res = solve([eq1, eq2, eq3, eq4, eq6, eq7], (股, 全弦, 界弦, r, a, x))

   2-element Vector{NTuple{6, Sym}}:
    (15, 17, 10, 2, 6, 2)
    (15, 17, 10, 2, 6, 7)

股 = 15;  全弦 = 17;  界弦 = 10;  円の直径 = 2;  a = 6;  x = 7

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (股, 全弦, 界弦, r, a, x) = (15, 17, 10, 2, 6, 7)
   @printf("股 = %d;  全弦 = %d;  界弦 = %d;  円の直径 = %d;  a = %d;  x = %d\n",
       股, 全弦, 界弦, r, a, x)
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   circle(r, r, r)
   circle(x, r, r)
   segment(0, 鉤, a, 0, :blue)
   if more
       point(a, 0, " a", :blue)
       point(股, 0, " 股", :black)
       point(0, 鉤, " 鉤", :black, :left, :bottom)
       point(r, r, " (r,r)", :red)
       point(x, r, " (x,r)", :red)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
       plot!(ylims=(-0.6, 8.5))
   else
       plot!(showaxis=false)
   end
end;

 

算額（その329）

算額（その329）

四二 浦和市西堀（現さいたま市桜区西堀） 氷川神社 嘉永5年(1852)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県浦和市西堀 氷川神社 嘉永5年（文化9年の算額の修復）
http://www.wasan.jp/saitama/uhikawa.html

正三角形の中に大円，中円，小円が入っている。中円の直径が520寸であるとき，小円の直径を求めよ。

正三角形の一辺の長さを 2a とする。座標原点を正三角形の「底辺/2」に置く。
大円の半径，中心座標は a/√3，(0, a/√3)
中円の半径，中心座標を r2, (x2, r2)
小円の半径，中心座標を r3, (x3, r3)
とおき，連立方程式を a を未知数のまま r2, x2, r3, x3 について解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r2::positive, x2::positive, r3::positive, x3::positive;

r1 = a/sqrt(Sym(3))
eq1 = r2/(a - x2) - 1/sqrt(Sym(3))
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2

res = solve([eq1, eq2, eq3, eq4], (a, x2, r3, x3))

   2-element Vector{NTuple{4, Sym}}:
    (3*sqrt(3)*r2, 2*sqrt(3)*r2, 3*r2*(2 - sqrt(3))/2, 3*r2*(-1 + sqrt(3)))
    (3*sqrt(3)*r2, 2*sqrt(3)*r2, 3*r2*(sqrt(3) + 2)/2, 3*r2*(1 + sqrt(3)))

最初のものが適解である。

r2 = 520
3*r2*(2 - sqrt(3))/2

   209.0003700962758

中円の直径が 520寸のとき，小円の直径は 209寸あまりである。
算額の答えは 109寸あまりということであるが，記載誤りである。

   大円径 = 1560; 中円径 = 520; 小円径 = 209

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 520/2
   (a, x2, r3, x3) = (3*sqrt(3)*r2, 2*sqrt(3)*r2, 3*r2*(2 - sqrt(3))/2, 3*r2*(-1 + sqrt(3)))
   r1 = a/√3
   @printf("大円径 = %g;  中円径 = %g;  小円径 = %g\n", 2r1, 2r2, 2r3)
   plot([a, 0, -a, a], [0, a*sqrt(3), 0, 0], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   if more
       point(0, r1, " 大円:r1,(0,r1)", :red, :vcenter, :left)
       point(x2, r2, "中円:r2,(x2,r2) ", :blue, :right, :vcenter)
       point(x3, r3, " 小円:r3,(x3,r3)", :green, :left, :vcenter)
       point(a, 0, "a", :black, :left, :bottom)
       point(0, √3a, " (0,√3a)", :black)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その328）

算額（その328）

埼玉県吉見町 安楽寺 文政5年(1822)
http://www.wasan.jp/saitama/anrakuji.html

菱形の中に直径の同じ 4 個の円を入れる。菱長（長い方の対角線）の二乗と菱横（短い方の対角線）の二乗の和が 400歩，外積（菱形の面積から 4 個の円の面積の和を除いた面積; 図で灰色の部分の面積）が 45.44 歩である。
菱長，菱横，円の直径，菱面（菱形の一辺の長さ）を求めよ。

菱長，菱横，円の直径，菱面 を 2a, 2b, 2r とおく。

方程式を解く前に，ここで用いられている円積率（π/4 の近似値，つまり円の面積 ≒ 円積率 * 直径の二乗）は当時広く採用されていた 0.75（すなわち π ≒ 3 とする）ではないので，採用した円積率がいくつなのか調べていく途中で，「吉田光由『塵劫記』(1627)が3.16を使っている」という記述を見つけた（https://www.ndl.go.jp/math/s1/c4.html）。つまり，3.16/4 = 0.79 である。
これによれば「外積が 45.44 歩」という数値の意味がわかる。

a > b なので 2 組目のものが適解である。
すなわち，菱長 = 16寸，菱横 = 12寸, 円径の直径 = 4寸，菱面 = 10寸 である。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, r::positive, 菱面::positive;

円積率 = 0.79  # = 3.16/4
eq1 = (2a)^2 + (2b)^2 - 400
eq2 = a + b - sqrt(a^2 + b^2) - 2r
eq3 = 2a*b - 4(円積率*(2r)^2) - 45.44
eq4 = a^2 + b^2 - 菱面^2

res = solve([eq1, eq2, eq3, eq4], (a, b, r, 菱面))

   2-element Vector{NTuple{4, Sym}}:
    (6.00000000000000, 8.00000000000000, 2.00000000000000, 10.0000000000000)
    (8.00000000000000, 6.00000000000000, 2.00000000000000, 10.0000000000000)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, r, 菱面) = (8, 6, 2, 10)
   plot([a, 0, -a, 0, a], [0, b, 0, -b, 0], linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=:gray80)
   circlef(r, r, r, :red2)
   circlef(-r, r, r, :red2)
   circlef(r, -r, r, :red2)
   circlef(-r, -r, r, :red2)
   circle4(r, r, r, :black)
   if more
       point(a, 0, "a", :green, :left, :bottom)
       point(0, b, "   b", :green, :left, :center)
       point(r, r, "      (r,r)", :black, :left, :center)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その327）

算額（その327）

埼玉県吉見町 安楽寺 文政5年(1822)
http://www.wasan.jp/saitama/anrakuji.html

鉤股弦に中鉤を引き，二分割された領域に大円と小円が鉤，股，弦，中鉤に内接している。小円の直径は大円の直径より 4 寸短く，弦は股より10 寸長い。
鉤，股，弦，大円の直径，小円の直径を求めよ。

一度に（描画に必要な変数も含めて）全ての変数を求めるのは solve() にとっては重荷なので，要求されている変数だけを求める連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, 鉤::positive, 股::positive, 弦::positive, 長弦::positive, 中鉤::positive;

eq1 = 2r1 -2r2 - 4
eq2 = 中鉤 * 弦 / 2 - 鉤 * 股 / 2
eq3 = 長弦 + 中鉤 - 股 - 2r1
eq4 = (弦 - 長弦) + 中鉤 - 鉤 - 2r2
eq5 = 鉤^2 + 股^2 - 弦^2
eq6 = 弦 - 股 - 10
eq7 = 中鉤/長弦 - 鉤/股

res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (鉤, 股, 弦, r1, r2, 中鉤, 長弦))

   1-element Vector{NTuple{7, Sym}}:
    (30, 40, 50, 8, 6, 24, 32)

鉤 = 30寸，股 = 40寸，弦 = 50寸，大円直径 = 16寸，小円直径 = 12寸である。

図を描くために，大円と小円の中心座標 (x1, r1), (股 - r2, y2) の 2変数 x1, y2 を求める連立方程式を解く（鉤，股，r1, r2 は前段階の結果を引き継ぐ）。
結果は 2 通り得られるが，最初のものが適解である。

@syms r1::positive, r2::positive, 鉤::positive, 股::positive, 弦::positive, 長弦::positive, 中鉤::positive,
   x1::positive, y2::positive;

(鉤, 股, 弦, r1, r2, 中鉤, 長弦) = (30, 40, 50, 8, 6, 24, 32)
eq8 = distance(0, 0, 股, 鉤, x1, r1) - r1^2
eq9 = distance(0, 0, 股, 鉤, 股 - r2, y2) - r2^2
res2 = solve([eq8, eq9], (x1, y2))

   2-element Vector{Tuple{Sym, Sym}}:
    (24, 18)
    (24, 33)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, r1, r2, 中鉤, 長弦) = (30, 40, 50, 8, 6, 24, 32)
   (x1, y2) = (24, 18)
   (x, y) = 長弦 .* (4/5, 3/5)  # 弦と中鉤の交点座標（弦と中鉤は直交する）
   plot([0, 股, 股, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   circle(x1, r1, r1)
   circle(股 - r2, y2, r2, :blue)
   segment(x, y, 股, 0, :green)
   if more
       point(x1, r1, "\n大円:r1\n(x1, r1)", :red, :center)
       point(股 - r2, y2, "\n小円:r2\n(股-r2,y2)", :blue, :center)
       point(股, 0, " A", :black, :left, :bottom)
       point(股, 鉤, " B", :black)
       point(0, 0, "   O", :black, :left, :bottom)
       point(x, y, "C ", :black, :right, :bottom)
       annotate!(10, 20, text("AB: 鉤\nAO: 股\nBO: 弦\nAC: 中鉤\nCO: 長弦", :left, 10))
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;