算額（その366）

算額（その366）

山形県山形市小白川町 天満社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

正方形内に二等辺三角形，大円，中円，小円が入っている。大円の直径を与えられたときに，中円の直径を求めよ。

正方形の左下隅の座標を (0, 0)
大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (r1 - r3, r3)
二等辺三角形の底辺の長さを 2b とする。

算額（その361）の簡単なものである。
https://blog.goo.ne.jp/r-de-r/e/ff37569366bf01f983d6cec055ecc257

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, b::positive;

eq1 = 2(r1 - r3)^2 - (r1 + r3)^2
eq2 = distance(0, 2r1, b, 0, 0, r2) - r2^2
eq3 = distance(0, 2r1, b, 0, r1 - r3, r3) - r3^2

res = solve([eq1, eq2, eq3], (r2, r3, b))

   3-element Vector{Tuple{Sym, Sym, Sym}}:
    (r1/2, r1*(3 - 2*sqrt(2)), sqrt(2)*r1/2)
    (-11*r1/4 - 19*r1*sqrt(203753 - 142008*sqrt(2))/1156 - 3*sqrt(2)*r1*sqrt(203753 - 142008*sqrt(2))/578 + 3*sqrt(2)*r1/2, r1*(3 - 2*sqrt(2)), r1*(-2 + 3*sqrt(2))/2)
    (-11*r1/4 + 3*sqrt(2)*r1*sqrt(203753 - 142008*sqrt(2))/578 + 19*r1*sqrt(203753 - 142008*sqrt(2))/1156 + 3*sqrt(2)*r1/2, r1*(3 - 2*sqrt(2)), r1*(-2 + 3*sqrt(2))/2)

1番目のものが適解である。
r2 = r1/2 で，中円の直径は大円の直径の 1/2 である。

 r2 = 0.5;  r3 = 0.171573;  b = 0.707107;  r2/r1 = 0.5

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r2, r3, b) = (r1/2, r1*(3 - 2*sqrt(2)), sqrt(2)*r1/2)
   @printf("r2 = %g;  r3 = %g;  b = %g;  r2/r1 = %g\n", r2, r3, b, r2/r1)
   plot([-r1, r1, r1, -r1, -r1], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
   plot!([-b, b, 0, -b], [0, 0, 2r1, 0], color=:green, lw=0.5)
   circle(0, r1, r1)
   circle(0, r2, r2, :blue)
   circle(r1 - r3, r3, r3, :magenta)
   circle(r3 - r1, r3, r3, :magenta)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(r1, 0, " r1", :black, :left, :bottom)
       point(b, 0, "b ", :black, :right, :bottom)
       point(0, 2r1, " 2r1", :black, :left, :bottom)
       point(0, r1, " r1", :red)
       point(0, r2, " r2", :blue)
       point(r1 - r3, r3, "(r1-r3,r3)", :magenta, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その365）

算額（その365）

山形県山形市小白川町 天満社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

正方形内に大小の正三角形と甲円，乙円がある。甲円の直径が与えられたときに，乙円の直径を求めよ。

正方形の辺の長さを a とする。

甲円の半径と中心座標を r1, (x1, r1)
容円の半径と中心座標を r2, (a - r2, y2)
とする。

以下の連立方程式を解く。a = 1 としても乙円と甲円の比を求める場合には一般性を失わない。

include("julia-source.txt");

using SymPy

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

a = 1
eq1 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x1, r1) - r1^2
eq2 = distance(0, (sqrt(Sym(3)) - 1)a, (sqrt(Sym(3)) - 1)a, 0, x1, r1) - r1^2
eq3 = distance(a, 0, a/2, sqrt(Sym(3))a/2, a - r2, y2) - r2^2
eq4 = distance(a, a, (sqrt(Sym(3)) - 1)a, 0, a - r2, y2) - r2^2

res = solve([eq1, eq2, eq3, eq4], (r1, x1, r2, y2))

   4-element Vector{NTuple{4, Sym}}:
    ((-11 - 19*sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + (19*sqrt(3)/3 + 11)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3 + sqrt(3), -1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16) + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5), 1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))
    ((-11 - 19*sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + (19*sqrt(3)/3 + 11)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3 + sqrt(3), -1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5) + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16), 1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))
    ((-11 - 19*sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + sqrt(3) + (19*sqrt(3)/3 + 11)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3, -1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16) + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5), 1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))
    ((-11 - 19*sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + sqrt(3) + (19*sqrt(3)/3 + 11)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3, -1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5) + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16), 1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))

1 番目のものが適解である。

甲円の半径

res[1][1] |> println

   (-11 - 19*sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + (19*sqrt(3)/3 + 11)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3 + sqrt(3)

乙円の半径

res[1][3] |> println

   (-22 - 10*sqrt(3))*(1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16) + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5)

甲円と乙円の式は複雑で長いが，二重根号を外して簡約化すると短い式になる。

r1 = res[1][1] |> sympy.sqrtdenest |> simplify
r2 = res[1][3] |> sympy.sqrtdenest |> simplify
r1 |> println
r2 |> println

   -sqrt(2) - 1/2 + sqrt(3)/2 + sqrt(6)/2
   -sqrt(2)/2 - 1/4 + sqrt(3)/4 + sqrt(6)/4

乙円の半径 / 甲円の半径 を求める。

r2/r1 |> simplify |> println

   1/2

乙円の直径は甲円の直径の 1/2 である。

   a = 1;  r1 = 0.176557;  r2 = 0.0882784;  r2/r1 = 0.5

x1, y2 も簡単な式になる。

res[1][2] |> sympy.sqrtdenest |> simplify |> println
res[1][4] |> sympy.sqrtdenest |> simplify |> println

   -sqrt(6) - sqrt(3)/2 + 3/2 + 3*sqrt(2)/2
   -sqrt(2)/4 + 1/4 + sqrt(3)/4

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (r1, x1, r2, y2) = res[1]
   @printf("a = %g;  r1 = %g;  r2 = %g;  r2/r1 = %g\n", a, r1, r2, r2/r1)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   plot!([a, 0, (√3 - 1)a, a], [a, (√3 - 1)a, 0, a], color=:blue, lw=0.5)
   plot!([0, a, a/2, 0], [0, 0, √3a/2], color=:red, lw=0.5)
   circle(x1, r1, r1, :magenta)
   circle(a - r2, y2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       r3m1a = (√3 - 1)a
       point(0, r3m1a, "   (√3-1)a", :blue, :left, :vcenter)
       point(r3m1a, 0, " (√3-1)a", :blue, :left, :bottom)
       point(a, 0, " a", :blue, :left, :bottom)
       point(0, a, " a", :blue, :left, :bottom)
       point(a/2, √3a/2, "(a/2,√3a/2)", :red, :right, :bottom)
       point(x1, r1, "甲円:r1,(x1,r1)", :magenta, :center, delta=-delta/2)
       point(a - r2, y2, "乙円:r2\n(a-r2,y2)", :green, :center, delta=2.3delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その364）

算額（その364）

山形県寒河江市大字柴橋字台下 伊豆神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

円内に三角形（甲斜，乙斜，丙斜）と容円がある。甲斜，乙斜，丙斜の長さがそれぞれ 31寸，26寸，8寸，容円の直径が 9寸のとき，矢の長さを求めよ。

甲斜，乙斜，丙斜 の長さをそれぞれ 甲，乙，丙とする。
乙斜と丙斜の交点座標を (x, y) とする。
外円の半径と中心座標を r0, (0, 0)
容円の半径と中心座標を r1, (x1, y1)
小円の半径と中心座標を r3, (x3, r3)
とする。
甲斜と y 軸の交点座標は sqrt(r0^2 + (甲/2)^2) である。
以下の連立方程式を解く。

三角形は確定しているので，甲斜に対しての相対的頂点座標 (x, y) を求めておく。

include("julia-source.txt");

using SymPy

@syms 甲::positive, 乙::positive, 丙::positive,
     r0::positive, r1::positive, x1::negative, y1::positive,
     x::negative, y::positive;
@syms x, y

(甲, 乙, 丙) = (31, 26, 8)
eq1 = (甲//2 - x)^2 + y^2 - 乙^2
eq2 = (x + 甲//2)^2 + y^2 - 丙^2
solve([eq1, eq2], (x, y))[1]

   (-306/31, -91*sqrt(15)/62)

他の座標を求める。

r1 = 9//2
y0 = sqrt(r0^2 - (甲//2)^2)
x = -306//31
y = y0 - 91*sqrt(15)/62
eq1 = x1^2 + y1^2 - (r0 - r1)^2
eq2 = distance(x, y, 甲//2, y0, x1, y1) - r1^2
eq3 = distance(x, y, -甲//2, y0, x1, y1) - r1^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=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)
   (r0, x1, y1) = u
   return [
     x1^2 + y1^2 - (r0 - 9/2)^2,  # eq1
     (-1385976137815419915050000000000000*x1/6496359999999999653635817057725169 + 6185822500000000000000000000000000*y1/6496359999999999653635817057725169 - 3092911250000000000000000000000000*sqrt(4*r0^2 - 961)/6496359999999999653635817057725169 + 21482630136139008683275000000000000/6496359999999999653635817057725169)^2 + (310537499999999653635817057725169*x1/6496359999999999653635817057725169 - 1385976137815419915050000000000000*y1/6496359999999999653635817057725169 + 692988068907709957525000000000000*sqrt(4*r0^2 - 961)/6496359999999999653635817057725169 - 9626662499999989262710328789480239/12992719999999999307271634115450338)^2 - 81/4,  # eq2
     (307505195230503210650000000000000*x1/615039999999999653635817057725169 + 304502500000000000000000000000000*y1/615039999999999653635817057725169 - 152251250000000000000000000000000*sqrt(4*r0^2 - 961)/615039999999999653635817057725169 + 4766330526072799765075000000000000/615039999999999653635817057725169)^2 + (310537499999999653635817057725169*x1/615039999999999653635817057725169 + 307505195230503210650000000000000*y1/615039999999999653635817057725169 - 153752597615251605325000000000000*sqrt(4*r0^2 - 961)/615039999999999653635817057725169 + 9626662499999989262710328789480239/1230079999999999307271634115450338)^2 - 81/4,  # eq3
   ]
end;

iniv = [big"22.0", -8.5, 15]
res = nls(H, ini=iniv);
println(res);

 (BigFloat[24.05739177340256001701446698368955710351888843775921149169035805392494843246005, -8.425291930678338067678088530312220772526057834341420944064372854532599673572294, 17.6495333893335722390792552249254553187347022853655274345347471331828351144709], true)

   r0 = 24.0574;  x1 = -8.42529;  y1 = 17.6495
   x = -9.87097;  y = 12.714
   矢 = 5.658802182377457
   甲斜長 = 31
   乙斜長 = 26.000000000000004
   丙斜長 = 8.0

矢の長さは r0 - y0 = 5.658802182377457 であるが，術では 5 寸と言い切っており，疑問が残る。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, x1, y1) = Float64.(res[1])
   @printf("r0 = %g;  x1 = %g;  y1 = %g\n", r0, x1, y1)
   (甲, 乙, 丙) = (31, 26, 8)
   r1 = 9//2
   y0 = sqrt(r0^2 - (甲//2)^2)
   x = -306//31
   y = y0 - 91*sqrt(15)/62
   @printf("x = %g;  y = %g\n", x, y)
   矢 = r0 - y0
   println("矢 = $矢")
   println("甲斜長 = $甲")
   println("乙斜長 = $(sqrt((甲/2 - x)^2 + (y0 - y)^2))")
   println("丙斜長 = $(sqrt((x + 甲/2)^2 + (y0 - y)^2))")
   plot()
   circle(0, 0, r0)
   circle(x1, y1, r1, :blue)
   plot!([x, 甲/2, -甲/2, x], [y, y0, y0, y], color=:black, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, y0, " y0=sqrt(r0^2-甲^2/4)", :red, :left, :bottom, delta=delta/2)
       point(0, r0, " r0", :red, :left, :bottom, delta=delta/2)
       point(甲/2, y0, " (甲/2,y0)")
       point(x, y, "(x,y)", delta=-delta/2)
       point(x1, y1, "(x1,y1)", :blue, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;