算額（その1688）

算額（その1688）

山形県新庄市 戸沢神社 文政元年(1818)
http://www.wasan.jp/yamagata/tozawa.html

キーワード：正方形3個，長方形
#Julia, #SymPy, #算額, #和算，#数学

長方形の中に大中小 3 個の正方形を容れる。長方形の長辺，短辺の長さが与えられたとき，小正方形の一辺の長さはいかほどか。

長方形の長辺，短辺を a, b
点 A の座標を (0, y2)
点 B の座標を (x2, y1)
点 C の座標を (x3, y6)
点 D の座標を (x1, b)
点 E の座標を (x5, 0)
点 F の座標を (a, y3)
点 G の座標を (x4, y4)
点 H の座標を (a, y5)
点 I の座標を (x6, b)
とおき，以下の連立方程式を解く。
SymPy の能力上，解析解は求められないので，数値解を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms x1::positive, x2::positive, x3::positive, x4::positive, x5::positive, x6::positive, a::positive,
      y1::positive, y2::positive, y3::positive, y4::positive, y5::positive, y6::positive, b::positive
AB = (x2^2 + (y2 - y1)^2)
eq1 = AB - ((x3 - x2)^2 + (y6 - y1)^2) |> expand
eq2 = AB - ((x3 - x1)^2 + (b - y6)^2) |> expand
eq3 = AB - (x1^2 + (b - y2)^2) |> expand
BE = (x5 - x2)^2 + y1^2
eq4 = BE - ((a - x5)^2 + y3^2) |> expand
eq5 = BE - ((a - x4)^2 + (y4 - y3)^2) |> expand
eq6 = BE - ((x4 - x2)^2 + (y4 - y1)^2) |> expand
GH = (a - x4)^2 + (y5 - y4)^2
eq7 = GH - ((a - x6)^2 + (b - y5)^2) |> expand
eq8 = GH - ((x6 - x3)^2 + (b - y6)^2) |> expand
eq9 = GH - ((x4 - x3)^2 + (y6 - y4)^2) |> expand
eq10 = x2 - (b - y2)
eq11 = (a - x2) - y4
eq12 = (a - x3) - (b - y4);

# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12],
#     (x1, x2, x3, x4, x5, x6, y1, y2, y3, y4, y5, y6))

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")
println(eq12, ",  # eq12")

    2*x2*x3 - x3^2 - 2*y1*y2 + 2*y1*y6 + y2^2 - y6^2,  # eq1
    -b^2 + 2*b*y6 - x1^2 + 2*x1*x3 + x2^2 - x3^2 + y1^2 - 2*y1*y2 + y2^2 - y6^2,  # eq2
    -b^2 + 2*b*y2 - x1^2 + x2^2 + y1^2 - 2*y1*y2,  # eq3
    -a^2 + 2*a*x5 + x2^2 - 2*x2*x5 + y1^2 - y3^2,  # eq4
    -a^2 + 2*a*x4 + x2^2 - 2*x2*x5 - x4^2 + x5^2 + y1^2 - y3^2 + 2*y3*y4 - y4^2,  # eq5
    2*x2*x4 - 2*x2*x5 - x4^2 + x5^2 + 2*y1*y4 - y4^2,  # eq6
    -2*a*x4 + 2*a*x6 - b^2 + 2*b*y5 + x4^2 - x6^2 + y4^2 - 2*y4*y5,  # eq7
    a^2 - 2*a*x4 - b^2 + 2*b*y6 - x3^2 + 2*x3*x6 + x4^2 - x6^2 + y4^2 - 2*y4*y5 + y5^2 - y6^2,  # eq8
    a^2 - 2*a*x4 - x3^2 + 2*x3*x4 - 2*y4*y5 + 2*y4*y6 + y5^2 - y6^2,  # eq9
    -b + x2 + y2,  # eq10
    a - x2 - y4,  # eq11
    a - b - x3 + y4,  # eq12

function driver(a, b)
    function H(u)
        (x1, x2, x3, x4, x5, x6, y1, y2, y3, y4, y5, y6) = u
        return [
            x2^2 - (-x2 + x3)^2 + (-y1 + y2)^2 - (-y1 + y6)^2,  # eq1
            x2^2 - (b - y6)^2 - (-x1 + x3)^2 + (-y1 + y2)^2,  # eq2
            -x1^2 + x2^2 - (b - y2)^2 + (-y1 + y2)^2,  # eq3
            y1^2 - y3^2 - (a - x5)^2 + (-x2 + x5)^2,  # eq4
            y1^2 - (a - x4)^2 + (-x2 + x5)^2 - (-y3 + y4)^2,  # eq5
            y1^2 - (-x2 + x4)^2 + (-x2 + x5)^2 - (-y1 + y4)^2,  # eq6
            (a - x4)^2 - (a - x6)^2 - (b - y5)^2 + (-y4 + y5)^2,  # eq7
            (a - x4)^2 - (b - y6)^2 - (-x3 + x6)^2 + (-y4 + y5)^2,  # eq8
            (a - x4)^2 - (-x3 + x4)^2 + (-y4 + y5)^2 - (-y4 + y6)^2,  # eq9
            -b + x2 + y2,  # eq10
            a - x2 - y4,  # eq11
            a - b - x3 + y4,  # eq12
        ]
    end;

    iniv = BigFloat[1.2, 6.4, 7.6, 8.8, 9.6, 10.8, 2.4, 3.6, 3.2, 5.6, 6.8, 8.8]
    res = nls(H, ini=iniv)
    return res[1]
end;

たとえば，長方形の長辺，短辺が 11, 10 のとき，小正方形の一辺の長さ(IH)は 3.05287 である。
A = (0, 4.8)
B = (5.2, 3.2)
C = (6.8, 8.4)
D = (1.6, 10)
E = (7.8, 0)
F = (11, 2.6)
G = (8.4, 5.8)
H = (11, 7.4)
I = (9.4, 10)

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (x1, x2, x3, x4, x5, x6, y1, y2, y3, y4, y5, y6) = driver(a, b)
    len = sqrt((a - x6)^2 + (b - y5)^2)
    @printf("長方形の長辺，短辺が %g, %g のとき，小正方形の一辺の長さ(IH)は %g である。\n", a, b, len)
    @printf("A = (%g, %g)\n", 0, y2)
    @printf("B = (%g, %g)\n", x2, y1)
    @printf("C = (%g, %g)\n", x3, y6)
    @printf("D = (%g, %g)\n", x1, b)
    @printf("E = (%g, %g)\n", x5, 0)
    @printf("F = (%g, %g)\n", a, y3)
    @printf("G = (%g, %g)\n", x4, y4)
    @printf("H = (%g, %g)\n", a, y5)
    @printf("I = (%g, %g)\n", x6, b)
    plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:brown, lw=0.5)
    plot!([0, x2, x3, x1, 0], [y2, y1, y6, b, y2], color=:green, lw=0.5)
    plot!([x2, x5, a, x4, x2], [y1, 0, y3, y4, y1], color=:red, lw=0.5)
    plot!([x4, a, x6, x3, x4], [y4, y5, b, y6, y4], color=:blue, lw=0.5)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, y2, " A", :green, :left, delta=-delta)
        point(x2, y1, "B ", :green, :right, delta=-delta)
        point(x3, y6, " C", :green, :left, :vcenter)
        point(x1, b, " D", :green, :left, delta=-delta)
        point(x5, 0, "E", :green, :left, :bottom, delta=delta)
        point(a, y3, " F", :green, :left, :vcenter)
        point(x4, y4, "G", :green, :right, delta=-delta)
        point(a, y5, " H", :green, :left, :vcenter)
        point(x6, b, "I", :green, :right, delta=-delta)
        point(a, b, "(a,b)", :brown, :right, :bottom, delta=delta)
    end
end;

draw(11, 10, true)

算額（その1687）

算額（その1687）

九十九 江刺市 雨宝堂 現中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03007
https://w.atwiki.jp/sangaku/pages/205.html

キーワード：円7個，円弧1個，楕円1個
#Julia, #SymPy, #算額, #和算，#数学

円弧の中に楕円 1 個，甲円 2個，乙円 3 個，丙円 2個を容れる。甲円と乙円の直径が与えられたとき，円弧の直径を求める術を述べよ。

円弧の半径と中心座標を R, (0, R - 4r2)
甲円の半径と中心座標を r1, (x1, 0)
乙円の半径と中心座標を r2, (0, R - r2), (0, R - 3r2), (0, R - 5r2)
丙円の半径と中心座標を r3, (x3, y3)
楕円の長半径，短半径と中心座標を a, b, (0, R - 4r2); b = 2r2
楕円と円弧の接点座標を (x0, y0)
楕円と丙円の接点座標を　(x03, y03)
とおく。

丙円は円弧の直径には無関係なので，eq1〜eq6 で a, x1, x0, y0 を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive, x1::positive, r2::positive, a::positive, b::positive, x0::positive, y0::negative
b = 2r2
eq1 = x1^2 + r2^2 - (r1 + r2)^2
eq2 = 4x1^2 - 4(a^2 - b^2)*(b^2 - r1^2)/b^2  # 算法助術-公式84
eq3 = x0^2/a^2 + (y0 - R + 4r2)^2/b^2 - 1 |> factor |> numerator
eq4 = -b^2*x0/(a^2*(y0 - R + 4r2)) + x0/y0 |> factor |> numerator
eq4 = R*a^2 - 4a^2*r2 - a^2*y0 + 4r2^2*y0
eq5 = x0^2 + y0^2 - R^2;

res = solve([eq1, eq2, eq4, eq5], (a, x1, x0, y0))[2]

    (2*sqrt(2)*r2^(3/2)*sqrt(-1/(r1 - 2*r2)), sqrt(r1)*sqrt(r1 + 2*r2), sqrt((R*r1 - 2*R*r2 + 8*r2^2)*(R*r1 + 2*R*r2 - 8*r2^2))/r1, -2*r2*(-R + 4*r2)/r1)

eq13 = eq3(a => res[1], x1 => res[2], x0 => res[3], y0 => res[4]) |> simplify |> numerator
eq13 |> println

    4*r2^2*(R^2*r1^2 - 4*R^2*r1*r2 + 4*R^2*r2^2 + 16*R*r1*r2^2 - 32*R*r2^3 - 24*r1*r2^3 + 64*r2^4)

ans_R = solve(eq13, R)[2]  # 2 of 2
ans_R |> println

    2*(sqrt(6)*sqrt(r1)*r2^(3/2)*(r1 - 2*r2)^2 + 4*r2^2*(-r1^2 + 4*r1*r2 - 4*r2^2))/((r1 - 2*r2)*(r1^2 - 4*r1*r2 + 4*r2^2))

√r1, √r2 を s, t とおき，簡約化したあとにもとに戻すと以下のようになる。

@syms s, t
ans_R = ans_R(sqrt(r1) => s, sqrt(r2) => t) |> simplify |> x -> x(s => sqrt(r1), t => sqrt(r2))
ans_R |> display
ans_R |> println

    2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2)

甲円の直径が 20 寸，乙円の直径が 12 寸のとき，全円の直径は 60.3160084678767 寸である。

2*ans_R(r1 => 20/2, r2 => 12/2).evalf() |> println

    60.3160084678767

術は以下のようになっており，2*ans_R と同じである。

@syms 甲径, 乙径
A = 3甲径/8乙径
全円径 = 8(1 - √A)*乙径^2/(2乙径 - 甲径)
全円径 |> display

全円径(甲径 => 20, 乙径 => 12).evalf() |> println

    60.3160084678767

a, x0, y0 は以下の通りである。

# a
res[1] |> println
res[1](r1 => 20/2, r2 => 12/2).evalf() |> println

    2*sqrt(2)*r2^(3/2)*sqrt(-1/(r1 - 2*r2))
    29.3938769133981

# x0
res[2] |> println
res[2](r1 => 20/2, r2 => 12/2, R => ans_R).evalf() |> println

    sqrt(r1)*sqrt(r1 + 2*r2)
    14.8323969741913

# y0
res[3] |> println
res[3](r1 => 20/2, r2 => 12/2, R => ans_R).evalf() |> println

    sqrt((R*r1 - 2*R*r2 + 8*r2^2)*(R*r1 + 2*R*r2 - 8*r2^2))/r1
    29.2386551695722

算額の解はここまでである。
算額には丙円が描かれているが，丙円は算額の答え（全円の直径）には何の影響も与えない。
つまり，丙円は解答者を惑わす「お飾り」である。

以下は丙円を描くためのパラメータを求める連立方程式であるが，SymPy の能力的に解析解を求めることができないので，数値解を求める。

@syms R::positive, r1::positive, x1::positive, r2::positive, a::positive, b::positive, x0::positive, y0::negative,
      r3::positive, x3::positive, y3::positive, x03::positive, y03::positive
b = 2r2
R = 2*r2^(Sym(3)/2)*(sqrt(Sym(6))*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2)
(a, x1, x0, y0) = (2*sqrt(Sym(2))*r2^(Sym(3)/2)*sqrt(-1/(r1 - 2*r2)), sqrt(r1)*sqrt(r1 + 2*r2),
    sqrt((R*r1 - 2*R*r2 + 8*r2^2)*(R*r1 + 2*R*r2 - 8*r2^2))/r1,
    -2*r2*(-R + 4*r2)/r1);

eq6 = x3^2 + (R - r2 - y3)^2 - (r2 + r3)^2
eq7 = x3^2 + y3^2 - (R - r3)^2
eq8 = (x03 - x3)^2 + (y03 - y3)^2 - r3^2
eq9 = x03^2/a^2 + (y03 - R + 4r2)^2/b^2 - 1
eq10 = -b^2*x03/(a^2*(y03 - R + 4r2)) + (x3 - x03)/(y3 - y03);    

println(eq6, ",  # eq6")
println(eq7, ",  # eq7")
println(eq8, ",  # eq8")
println(eq9, ",  # eq9")
println(eq10, ",  # eq10")

    x3^2 - (r2 + r3)^2 + (2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) - r2 - y3)^2,  # eq6
    x3^2 + y3^2 - (2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) - r3)^2,  # eq7
    -r3^2 + (x03 - x3)^2 + (y03 - y3)^2,  # eq8
    -1 + (-2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) + 4*r2 + y03)^2/(4*r2^2) + x03^2*(-r1 + 2*r2)/(8*r2^3),  # eq9
    (-x03 + x3)/(-y03 + y3) + x03*(r1 - 2*r2)/(2*r2*(-2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) + 4*r2 + y03)),  # eq10

function H(u)
    (r3, x3, y3, x03, y03) = u
    return [
        x3^2 - (r2 + r3)^2 + (2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) - r2 - y3)^2,  # eq6
        x3^2 + y3^2 - (2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) - r3)^2,  # eq7
        -r3^2 + (x03 - x3)^2 + (y03 - y3)^2,  # eq8
        -1 + (-2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) + 4*r2 + y03)^2/(4*r2^2) + x03^2*(-r1 + 2*r2)/(8*r2^3),  # eq9
        (-x03 + x3)/(-y03 + y3) + x03*(r1 - 2*r2)/(2*r2*(-2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2) + 4*r2 + y03)),  # eq10
    ]
end;

r1 = 10;  r2 = 6
iniv = BigFloat[5, 11, 22.5, 10.2, 17.3];
res = nls(H, ini=iniv)

    ([5.1302159947337, 11.0029172476153, 22.47945720404743, 10.233847537564895, 17.407214314190472], true)

甲円の直径が 20 寸，乙円の直径が 12 寸のとき，丙円の半径は 5.1302159947337，中心座標は (11.0029172476153, 22.47945720404743) である。

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 2*r2^(3/2)*(sqrt(6)*sqrt(r1) - 4*sqrt(r2))/(r1 - 2*r2)
    (a, x1, x0, y0) = (2*sqrt(2)*r2^(3/2)*sqrt(-1/(r1 - 2*r2)), sqrt(r1)*sqrt(r1 + 2*r2), sqrt((R*r1 - 2*R*r2 + 8*r2^2)*(R*r1 + 2*R*r2 - 8*r2^2))/r1, -2*r2*(-R + 4*r2)/r1)
    b = 2r2
    @printf("r1 = %g;  r2 = %g;  R = %g;  a = %g;  x1 = %g;  x0 = %g;  y0 = %g;  b = %g\n", r1, r2, R, a, x1, x0, y0, b)
    y = R - 6r2
    x = sqrt(R^2 - y^2)
    θ = atand(y, x)
    plot()
    circle(0, 0, R, beginangle=θ, endangle=180 - θ, n=500)
    segment(-x, y, x, y, :red)
    circle(0, R - r2, r2, :blue)
    circle(0, R - 3r2, r2, :blue)
    circle(0, R - 5r2, r2, :blue)
    circle2(x1, R - 4r2, r1, :magenta)
    ellipse(0, R - 4r2, a, b)
    (r3, x3, y3, x03, y03) = [5.1302159947337, 11.0029172476153, 22.47945720404743, 10.233847537564895, 17.407214314190472]
    @printf("x3 = %g;  y3 = %g;  r3 = %g;  x03 = %g;  y03 = %g", x3, y3, r3, x03, y03)
    circle2(x3, y3, r3, :green)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, R, "R", :red, :center, :bottom, delta=delta)
        point(0, R - r2, "乙円:r2\n(0,R-r2)", :blue, :center, delta=-delta)
        point(0, R - 3r2, "乙円:r2\n(0,R-3r2)", :blue, :center, delta=-delta)
        point(0, R - 5r2, "乙円:r2\n(0,R-5r2)", :blue, :center, delta=-delta)
        point(x1, R - 4r2, "甲円:r1\n(x1,R-4r2)", :magenta, :center, delta=-delta)
        point(0, R - 2r2, "R-2r2", :black, :center, delta=-delta)
        point(0, R - 4r2, "R-4r2", :black, :center, delta=-delta)
        point(x0, y0, "(x0,y0)", :green, :left, :bottom, delta=delta)
        point(x03, y03, "(x03,y03)", :green, :right, delta=-delta, deltax=2delta)
        point(0, R - 6r2, "R-6r2", :red, :center, delta=-delta)
        dimension_line(0, R - 4r2, a, R - 4r2, "a", :green, :center, :bottom, delta=delta, deltax=3delta, length=0)
        point(x3, y3, "丙円:r3\n(x3,y3)", :green, :center, delta=-delta)
        plot!(xlims=(-R - 2delta, R + 12delta), ylims=(R - 6r2 - 5delta, R + delta))
    end
end;

draw(20/2, 12/2, true)

算額（その1686）

算額（その1686）

九十九 江刺市 雨宝堂 現中善観音堂 文政10年(1827)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03007
https://w.atwiki.jp/sangaku/pages/205.html

キーワード：円2個，正方形3個
#Julia, #SymPy, #算額, #和算，#数学

正方形の中に小さな正方形を 2 個容れ，甲円と乙円を 2 個ずつ容れる。乙円の直径が与えられたとき，甲円の直径はいかほどか。

甲円の半径を r1
乙円の半径を r2
正方形の一辺の長さを a
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r1::positive,
      r2::positive
b = a/2
eq1 = b/a- (sqrt(a^2 + b^2) - 4r1)/2r1
eq2 = a + b - sqrt(a^2 + b^2) - 2r2
res = solve([eq1, eq2], (r1, a))

# r1 甲円の半径
res[r1] |> simplify |> println

    r2*(5 + 3*sqrt(5))/10

甲円の半径 r1 は，乙円の半径 r2 の (5 + 3√5)/10 倍である。
たとえば，乙円の直径が 1 寸のとき，甲円の直径は (5 + 3√5)/10 = 1.170820393249937 寸である。

術は以下のようになっており，上記の解と一致している。

@syms 乙円径
甲円径 = (sqrt(1.8) + 1)*乙円径/2
甲円径(乙円径 => 1) |> println

    1.17082039324994

# a 正方形の一辺の長さ
res[a] |> simplify |> println

    r2*(sqrt(5) + 3)

正方形の一辺の長さ a は，乙円の半径の √5 + 3 倍である。
たとえば，乙円の直径が 1 寸のとき，正方形の一辺の長さは (√5 + 3)*(1/2) = 2.618033988749895 寸である。

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = r2*(5 + 3√5)/10
    a = r2*(√5 + 3)
    b = a/2
    @printf("乙円の直径が %g のとき，甲円の直径は %g である。なお，正方形の一辺の長さは %g である。\n", 2r2, 2r1, a)
    plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
    circle(r2, a - r2, r2)
    circle(a - r2, r2, r2)
    θ = atand(2)
    (Ax, Ay) = (b, a)
    (Bx, By) = (b + 2r1*sind(θ), a - 2r1*cosd(θ))
    (Cx, Cy) = (b, 0)
    (Dx, Dy) = (b - 2r1*sind(θ), 2r1*cosd(θ))
    (Ex, Ey) = ((2b - 2r1*sind(θ))/2, (a + 2r1*cosd(θ))/2)
    (Fx, Fy) = ((2b + 2r1*sind(θ))/2, (a - 2r1*cosd(θ))/2)
    segment(0, 0, Ax, Ay)
    segment(Cx, Cy, a, a)
    segment(Ax, Ay, Bx, By)
    segment(Ex, Ey, Fx, Fy)
    segment(Cx, Cy, Dx, Dy)
    (o1x, o1y) = ((Ax + Fx)/2, (Ay + Fy)/2)
    circle(o1x, o1y, r1, :blue)
    (o2x, o2y) = ((Ex + Cx)/2, (Ey + Cy)/2)
    circle(o2x, o2y, r1, :blue)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(Ax, Ay, "A", :black, :left, delta=-delta)
        point(Bx, By, " B", :black, :left, :vcenter)
        point(Cx, Cy, " C", :black, :left, :vcenter)
        point(Dx, Dy, "D", :black, :left, delta=-delta)
        point(Ex, Ey, "E", :black, :left, delta=-delta)
        point(Fx, Fy, " F", :black, :left, :vcenter)
        point(o1x, o1y, "甲円:r1", :blue, :center, delta=-delta)
        point(o2x, o2y, "甲円:r1", :blue, :center, delta=-delta)
        point(r2, a - r2, "乙円:r2", :red, :center, delta=-delta)
        point(a-r2, r2, "乙円:r2", :red, :center, delta=-delta)
        point(a, a, "(a,a)", :green, :right, :bottom, delta=delta)
    end
end;

draw(1/2, true)

 

算額（その1685）

算額（その1685）

七十一 岩手県一関市川崎町薄衣諏訪前 浪分神社(第4面) 明治35年(1902)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03088
https://w.atwiki.jp/sangaku/pages/277.html

キーワード：円9個，円弧2個
#Julia, #SymPy, #算額, #和算，#数学

2 個の円弧が交わる中に，大円 1 個，中円 4 個，小円 4 個を容れる。小円の直径が与えられたとき，大円の直径を求める術を述べよ。

円弧の半径と中心座標を R, (0, R - 2r2)
大円の半径と中心座標を r1, (0, 0)
中円の半径と中心座標を r2, (r1 + r2, 0)
小円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive, r2::positive,
      r3::positive, x3::positive3, y3::negative
r1 = 2r2
eq1 = (r1 + r2)^2 + (R - 2r2)^2 - (R - r2)^2
eq2 = x3^2 + (R - 2r2 - y3)^2 - (R - r3)^2
eq3 = (r1 + r2 - x3)^2 + y3^2 - (r2 + r3)^2
eq4 = x3^2 + y3^2 - (r1 + r3)^2
res = solve([eq1, eq2, eq3, eq4], (r2, R, x3, y3))[1]

    (7*r3/3, 14*r3, 5*r3, -8*r3/3)

大円の半径 r1 は，中円の半径 r2 の 2 倍である。
中円の半径 r2 は，小円の半径 r3 の 7/3 倍である。
すなわち，大円の半径 r1 は 小円の半径 r3 の 14/3 倍である。

たとえば，小円の直径が 1 のとき，大円の直径は 14/3 = 4.666666666666667 である。

function draw(r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, R, x3, y3) = (7*r3/3, 14*r3, 5*r3, -8*r3/3)
    r1 = 2r2
    @printf("小円の直径が %g のとき，大円の直径は %g である。\n", 2r3, 2r1)
    @printf("円弧の半径は %g，中心座標は (0, %g), (0, %g) である。\n", R, R - 2r2, 2r2 - R)
    y = R - 2r2
    x = sqrt(R^2 - y^2)
    θ = atand(y, x)
    plot()
    circle(0, R - 2r2, R, beginangle=180+θ, endangle=360−θ, n=500)
    circle(0, -R + 2r2, R, beginangle=θ, endangle=180−θ, n=500)
    circle(0, 0, r1, :blue)
    circle22(0, r2, r2, :magenta)
    circle2(r1 + r2, 0, r2, :magenta)
    circle4(x3, y3, r3, :green)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, 0, "大円:r1,(0,0)", :blue, :center, :bottom, delta=delta)
        point(0, r2, "中円:r2\n(0,r2)", :magenta, :center, :bottom, delta=delta)
        point(r1 + r2, 0, "中円:r2\n(r1+r2,0)", :magenta, :center, :bottom, delta=delta)
        point(x3, y3, "小円:r3,(x3,y3)", :green, :left, delta=-10delta, deltax=-3delta)
    end
end;

draw(1/2, true)

 

算額（その1684）

算額（その1684）

六十四 岩手県一関市花泉町金沢大門沢 大門神社・大門観世音菩薩 明治26年(1893)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03077
https://w.atwiki.jp/sangaku/pages/293.html

キーワード：球1個，円錐3個，3次元
#Julia, #SymPy, #算額, #和算，#数学

盤上に 3 個の円錐が底面の中心が円周上にあり底面で互いに外接している。その中心に1個の球を容れる。球は盤面に接し，またそれぞれの円錐と 1 点で外接している。円錐の底面の直径と高さが与えられたとき，球の直径はいかほどか。

算額（その1683） の類題である。

いつもどおり，山村の図は不適切なものである。球は盤面に載っていなければならない。

球の半径と中心座標を r1, (0, 0, r1)
円錐の高さと底面の半径と中心座標を h, r2, (a*cos(θ), a*sin(θ)); θ = 0°, 120°, 240°
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, r1::positive, r2::positive, h::positive
eq1 = dist2(a - r2, 0, a, h, 0, r1, r1)
eq2 = a*sind(Sym(60)) - r2
res = solve([eq1, eq2], (r1, a))[2]  # 2 of 2

    (r2^2*(-3 + 2*sqrt(3))/(3*h) + r2*sqrt(3*h^2*(7 - 4*sqrt(3)) + r2^2*(3 - 2*sqrt(3))^2)/(3*h), 2*sqrt(3)*r2/3)

# r1 球の半径
res[1] |> together

たとえば円錐の底面の半径が 1/2，高さが 3/2 のとき，球の直径は 0.214635531637327 である。

2 * res[1](r2 => 1/2, h => 3/2).evalf() |> println

    0.214635531637327

術は以下のようであり，前述の解と一致する。

下径 = 1
高 = 1.5
A = sqrt((下径/2)^2 + 高^2)
天 = 下径/2 + 高 - A
球径 = (2/sqrt(3) - 1)*((下径 - 天)*下径)/天

    0.21463553163732685

function draw(r2, h)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, a) = (r2^2*(-3 + 2*sqrt(3))/(3*h) + r2*sqrt(3*h^2*(7 - 4*sqrt(3)) + r2^2*(3 - 2*sqrt(3))^2)/(3*h), 2*sqrt(3)*r2/3)
    @printf("円錐の底面の直径が %g のとき，球の直径は %g である。\n", 2r2, 2r1)
    p1 = plot([a - r2, a + r2, a, a - r2], [0, 0, h, 0], color=:blue, lw=0.5, xlabel="x-axis", ylabel="z-axis")
    circle(0, r1, r1)
    point(0, r1, "球:r1,(0,r1)", :red, :center, delta=-0.1r1)
    point(a, 0, "a", :blue, :center, :bottom, delta=0.1r1)
    point(a + r2, 0, "a+r2", :blue, :center, :bottom, delta=0.1r1)
    p2 = plot(xlabel="x-axis", ylabel="y-axis")
    rotate(a, 0, r2, :blue, angle=120)
    circle(0, 0, r1)
    point(a, 0, "(a,0)", :blue, :center, delta=-0.1r1)
    plot(p1, p2)
end;

draw(1/2, 3/2)

算額（その1683）

算額（その1683）

四十六 岩手県一関市真滝 熊野白山滝神社 明治31年(1898)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03081
https://w.atwiki.jp/sangaku/pages/329.html

キーワード：球1個，円錐5個，3次元
#Julia, #SymPy, #算額, #和算，#数学

盤上に 5 個の円錐が底面の中心が円周上にあり底面で互いに外接している。その中心に1個の球を容れる。球は盤面に接し，またそれぞれの円錐と 1 点で外接している。球の高さと円錐の高さは等しい。円錐の底面の直径が与えられたとき，球の直径はいかほどか。

いつもどおり，山村の図は不適切なものである。なぜか，円錐とも見えないものが 6 個描かれていたりする。

球の半径と中心座標を r1, (0, 0, r1)
円錐の底面の半径と中心座標を r2, (a*cos(θ), a*sin(θ)); θ = 0°, 72°, 144°, 216°, 288°
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, r1::positive, r2::positive
eq1 = dist2(a - r2, 0, a, 2r1, 0, r1, r1)
eq2 = a*sind(Sym(36)) - r2
res = solve([eq1, eq2], (r1, a))[1]

    (sqrt(-2*sqrt(2)*r2^2*sqrt(5 - sqrt(5)) + 8*r2^2)/sqrt(5 - sqrt(5)), 2*sqrt(2)*r2/sqrt(5 - sqrt(5)))

# r1
res[1] |> simplify

円錐の底面の半径 r2 が 1/2 のとき，球の半径 r1 は 0.546151438315381 である。

res[1](r2 => 1/2).evalf() |> println

    0.546151438315381

術は以下の通りで，上の解と一致する。

円錐径 = 1
天 = sqrt(0.8) + 2
球径 = sqrt(天 - sqrt(天))*円錐径

    1.0923028766307612

function draw(r2)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, a) = (sqrt(-2*sqrt(2)*r2^2*sqrt(5 - sqrt(5)) + 8*r2^2)/sqrt(5 - sqrt(5)), 2*sqrt(2)*r2/sqrt(5 - sqrt(5)))
    @printf("円錐の底面の直径が %g のとき，球の直径は %g である。\n", 2r2, 2r1)
    p1 = plot([a - r2, a + r2, a, a - r2], [0, 0, 2r1, 0], color=:blue, lw=0.5, xlabel="x-axis", ylabel="z-axis")
    circle(0, r1, r1)
    p2 = plot(xlabel="x-axis", ylabel="y-axis")
    rotate(a, 0, r2, :blue, angle=72)
    circle(0, 0, r1)
    plot(p1, p2)
end;

draw(1/2)

 

算額（その1682）

算額（その1682）

四十五 岩手県一関市真滝 熊野白山滝神社 明治8年(1875)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03061
https://w.atwiki.jp/sangaku/pages/327.html

キーワード：円5個，斜線2本，直線上
#Julia, #SymPy, #算額, #和算，#数学

直線上に全円 1 個，大円 2 個，小円 2 個が載っている。全円と大円，全円と小円は外接しており，大円と小円は共通接線を持つ。大円と小円の直径が与えられたとき，全円の直径はいかほどか。

いつもどおり，山村の図は不適切なものである。

全円の半径と中心座標を r1, (0, r1)
大円の半径と中心座標を r2, (c, r2)
小円の半径と中心座標を r3, (b, r3)
大円と小円の共通接線と x 軸の交点座標を (a, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::negative, b::negative, c::positive, θ::positive,
      r1::positive, r2::positive, r3::positive
eq1 = c - 2sqrt(r1*r2)
eq2 = -b - 2sqrt(r1*r3)
eq3 = r3/(b - a) - r2/(c - a)
eq4 = dist2(a, 0, 0, r2, c, r2, r2);

res = solve([eq1, eq2, eq3, eq4], (r1, a, b, c))[2]  # 2 of 2

    (r2^3*(-r2^2 + 2*r2*r3 - r3^2)/(4*(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)), -r2*sqrt(r2*sqrt(-(r2^2 - 2*r2*r3 + r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) - 1)*sqrt(r2*sqrt(-(r2^2 - 2*r2*r3 + r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) + 1), r2*r3*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) + (-r2 + r3)*sqrt(r2*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) - 1)*sqrt(r2*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) + 1), r2^2*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)))

全円の半径は大円と小円の半径で以下のように表される。術で述べられている式は不適切なものである。

たとえば，大円の直径が 13，小円の直径が 3 のとき，全円の直径は 22.3675 である。

function draw(r2, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, a, b, c) = (r2^3*(-r2^2 + 2*r2*r3 - r3^2)/(4*(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)), -r2*sqrt(r2*sqrt(-(r2^2 - 2*r2*r3 + r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) - 1)*sqrt(r2*sqrt(-(r2^2 - 2*r2*r3 + r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) + 1), r2*r3*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) + (-r2 + r3)*sqrt(r2*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) - 1)*sqrt(r2*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)) + 1), r2^2*sqrt((-r2^2 + 2*r2*r3 - r3^2)/(2*r2^(5/2)*r3^(3/2) - r2^4 + 3*r2^3*r3)))
    θ = atan(r2/(c - a))
    @printf("大円の直径が %g，小円の直径が %g のとき，全円の直径は %g である。\n", 2r2, 2r3, 2r1)
    plot()
    circle(0, r1, r1)
    circle2(c, r2, r2, :magenta)
    circle2(b, r3, r3, :blue)
    abline(a, 0, tan(2θ), a, -a + 0.5r2, :green)
    abline(-a, 0, -tan(2θ), -a, a - 0.5r2, :green)
    segment(-c - r2, 0, c + r2, 0, lw=1)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, r1, "全円:r1,(0,r1)", :red, :center, :bottom, delta=delta)
        point(c, r2, "大円:r2,(c,r2)", :magenta, :center, :bottom, delta=delta)
        point(b, r3, " 小円:r3,(b,r3)", :blue, :left, :vcenter)
        point(a, 0, "a", :green, :center, delta=-2delta)
        point(b, 0, "b", :blue, :center, delta=-2delta)
        point(c, 0, "c", :magenta, :center, delta=-2delta)
        ylims!(-7delta, 2r1 + delta)
    end
end;

draw(13/2, 3/2, true)

算額（その1681）

算額（その1681）

四十五 岩手県一関市真滝 熊野白山滝神社 明治8年(1875)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図
https://w.atwiki.jp/sangaku/pages/327.html

キーワード：円1個，楕円3個，正三角形
#Julia, #SymPy, #算額, #和算，#数学

楕円 3 個が正三角形を囲み，正三角形に内接する円と楕円は外接している。楕円の短径が 1 寸のとき長径はいかほどか。

山村の図はいつもどおり，ひどいものである。

内接円の半径と中心座標を r, (0, 0)
楕円の長半径，短半径と中心座標を a, b, (0, 3r/2)
2 個の楕円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r::positive, x0::positive, y0::positive
r = 2b
eq1 = x0^2/a^2 + (y0 - 3r/2)^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - 3r/2)) - 1/√Sym(3)
eq3 = √Sym(3)*y0 - x0;

res = solve([eq1, eq2, eq3], (a, x0, y0))[1]

    (2*sqrt(6)*b, 8*sqrt(3)*b/3, 8*b/3)

楕円の長半径は短半径の 2√6 = √24 倍である。

function draw(b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = 2b
    (a, x0, y0) = b .* (2√6, 8√3/3, 8/3)
    plot([√3r, 0, -√3r, √3r], [-r, 2r, -r, -r], color=:green, lw=0.5)
    ellipse(0, 3r/2, a, b, color=:red)
    ellipse(3r/2*cosd(30), -3r/2*sind(30), a, b, color=:red, φ=-120)
    ellipse(-3r/2*cosd(30), -3r/2*sind(30), a, b, color=:red, φ=120)
    circle(0, 0, r, :blue)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(x0, y0, "(x0,y0)", :red, :right, :bottom, delta=delta/2)
        point(0, 2b, "2b", :blue, :center, delta=-delta)
        point(0, 4b, "4b", :red, :center, :bottom, delta=delta)
        point(0, 3b, " 楕円:a,b,(0,3b)", :red, :left, :vcenter)
        point(√3r, -r, " (√3r,-r)", :green, :left, :vcenter)
    end
end;

draw(1/2, true)

算額（その1680）

算額（その1680）

七十 岩手県一関市川崎町薄衣 浪分神社 文久年間
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円4個，楕円
#Julia, #SymPy, #算額, #和算，#数学

楕円と大円が交差しており，隙間に小円 3 個を容れる。楕円の長径，短径が与えられたとき，小円の直径はいかほどか。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
大円の半径と中心座標を r1, (0, r1 - b)
小円の半径と中心座標を r2, (d, 0), (0, 2r1 - b - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive, d::positive
eq1 = (2r1 - b) - (b + 2r2)
eq2 = d^2 + (r1 - b)^2 - (r1 + r2)^2
eq3 = 4d^2 - 4(a^2 - b^2)*(b^2 - r2^2)/b^2;  # 算法助術-公式84

res = solve([eq1, eq2, eq3], (r1, r2, d))[1]

    (2*a^2*b/(a^2 + 2*b^2), b*(a^2 - 2*b^2)/(a^2 + 2*b^2), 2*sqrt(2)*a*b*sqrt(a - b)*sqrt(a + b)/(a^2 + 2*b^2))

たとえば，長半径が 5/2，短半径が 2/2 のとき，小円の直径は 1.03030303030303 である。

# r2
res[2] |> println
2 * res[2](a => 5/2, b => 2/2) |> println

    b*(a^2 - 2*b^2)/(a^2 + 2*b^2)
    1.03030303030303

術は以下の通りで，同等の式である。

長径 = 5
短径 = 2
天 = 長径^2
地 = 2短径^2
実 = (天 - 地)*短径
小円径 = 実/(天 + 地)

    1.0303030303030303

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, d) = (2*a^2*b/(a^2 + 2*b^2), b*(a^2 - 2*b^2)/(a^2 + 2*b^2), 2*sqrt(2)*a*b*sqrt(a - b)*sqrt(a + b)/(a^2 + 2*b^2))
    @printf("長径 = %g，短径 = %g のとき，小円径 = %g，大円径 = %g\n", 2a, 2b, 2r2, 2r1)
    plot()
    ellipse(0, 0, a, b, color=:red)
    circle2(d, 0, r2, :blue)
    circle(0, -b + 2r1 - r2, r2, :blue)
    circle(0, r1 - b, r1, :magenta)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, r1 - b, "大円:r1,(0,r1-b)", :magenta, :center, delta=-delta)
        point(d, 0, "小円:r2,(d,0)", :blue, :center, delta=-delta)
        point(0, 2r1 - b - r2, "小円:r2,(0,2r1-b-r2)", :blue, :center, delta=-delta)
        point(a, 0, " a", :red, :left, :vcenter)
        point(0, b, "b", :red, :center, delta=-delta)
    end
end;

draw(5/2, 2/2, true)

算額（その1679）

5月1日に公開

算額（その1678）

算額（その1678）

宮城県塩竈市森山 塩竈神社 文化14年(1817) 昭和60年復元奉納
http://www.wasan.jp/miyagi/siogama4.html
キーワード：円3個，楕円
#Julia, #SymPy, #算額, #和算，#数学

楕円の中に，甲円，乙円，丙円を容れる。楕円の長径が 13 寸，短径が 5 寸，乙円の直径が 4 寸のとき，丙円の直径はいかほどか。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
甲円の半径と中心座標を r1, (d1, 0)
乙円の半径と中心座標を r2, (d2, 0)
丙円の半径と中心座標を r3, (x3, y3)
丙円と楕円の接点座標を (x0, y0)
とおき，以下の連立方程式を解く。

SymPy の能力的に，解析解を求めることができないので，数値解を求める。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a, b, r2, r1, r3, x3, y3, d1, d2, x0, y0
d1 = sqrt((a^2 - b^2)*(b^2 - r2^2))/b
d2 = sqrt((a^2 - b^2)*(b^2 - r1^2))/b
eq1 = a^4*(r1 - r2)^2 - 4b^2*(a^2*b^2 - a^2*r2*r1 - b^4)  # 算法助術-公式87
eq2 = (x3 + d1)^2 + y3^2 - (r3 + r2)^2
eq3 = (x3 - d2)^2 + y3^2 - (r3 + r1)^2
eq4 = (x0 - x3)^2 + (y0 - y3)^2 - r3^2
eq5 = x0^2/a^2 + y0^2/b^2 - 1
eq6 = -b^2*x0/(a^2*y0) + (x0 - x3)/(y0 - y3);
# solve([eq1, eq2, eq3, eq4, eq5, eq6], (r1, r3, x3, y3, x0, y0))

ans_r1 = solve(eq1, r1)[2]  # 2 of 2
ans_r1 |> println

    (2*b*sqrt((a - b)*(a + b)*(b - r2)*(b + r2)) + r2*(a^2 - 2*b^2))/a^2

ans_r1(a => 13/2, b => 5/2, r2 => 4/2).evalf() |> println

    2.47337278106509

function H(u)
    (r1, r3, x3, y3, x0, y0) = u
    return [
        a^4*(r1 - r2)^2 - 4*b^2*(a^2*b^2 - a^2*r1*r2 - b^4),  # eq1
        y3^2 - (r2 + r3)^2 + (x3 + sqrt((a^2 - b^2)*(b^2 - r2^2))/b)^2,  # eq2
        y3^2 - (r1 + r3)^2 + (x3 - sqrt((a^2 - b^2)*(b^2 - r1^2))/b)^2,  # eq3
        -r3^2 + (x0 - x3)^2 + (y0 - y3)^2,  # eq4
        -1 + y0^2/b^2 + x0^2/a^2,  # eq5
        (x0 - x3)/(y0 - y3) - b^2*x0/(a^2*y0),  # eq6
    ]
end;
a = 13/2
b = 5/2
r2 = 4/2
iniv = BigFloat[2.47337, 0.62, -1.67, 1.8, -1.7, 2.42]
res = nls(H, ini=iniv)

    ([2.473372781065089, 0.63, -1.6666666666666667, 1.783009316358785, -1.7333333333333334, 2.4094720491334933], true)

長径 = 13, 短径 = 5, 乙円の直径 = 4 のとき，丙円の直径は 1.26 である。

function draw(a, b, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    #r1 = (2*b*sqrt((a - b)*(a + b)*(b - r2)*(b + r2)) + r2*(a^2 - 2*b^2))/a^2
    (r1, r3, x3, y3, x0, y0) = [2.473372781065089, 0.63, -1.6666666666666667, 1.783009316358785, -1.7333333333333334, 2.4094720491334933]
    d2 = -sqrt((a^2 - b^2)*(b^2 - r2^2))/b
    d1 = sqrt((a^2 - b^2)*(b^2 - r1^2))/b
    @printf("長径 = %g, 短径 = %g, 乙円の直径 = %g のとき，丙円の直径は %g である。\n", 2a, 2b, 2r2, 2r3)
    plot()
    ellipse(0, 0, a, b)
    circle(d2, 0, r2)
    circle(d1, 0, r1, :blue)
    circle(x3, y3, r3, :green)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(x0, y0, "(x0,y0)", :green, :center, :bottom, delta=2delta)
        point(d1, 0, "甲円:r1,(d1,0)", :blue, :center, delta=-delta)
        point(d2, 0, "乙円:r2,(d2,0)", :red, :center, delta=-delta)
        point(x3, y3, "丙円:r3,(x3,y3)", :green, :left, delta=-delta, deltax=8delta)
        point(a, 0, " a", :black, :left, :vcenter)
        point(0, b, "b", :black, :center, :bottom, delta=delta)
    end
end;

draw(13/2, 5/2, 4/2, true)

算額（その1677）

算額（その1677）

岩手県紫波郡矢巾町岩清水 龍泉寺不動堂 明治12年(1879)
http://www.wasan.jp/iwate/ryusenji.html

キーワード：長精度計算
#Julia, #SymPy, #算額, #和算，#数学

2220446049250313080847263336181640625 の 26 乗根 = 25 の計算

Julia では，長精度数を big"数字列" のように指定するだけで，長精度計算をすることができる。

精度はデフォルトで十分なことが多いが，setprecision(BigFloat, 2000) などで指定することもできる。

big"2220446049250313080847263336181640625"^(1/big"26")

    25.00000000000000000000000000000000000000000000000000000000000000000000000000028

検算

big"25"^26

    2220446049250313080847263336181640625

なお，素因子分解も簡単に行える。

using Primes
factor(2220446049250313080847263336181640625)

    5^52

 

算額（その1676）

算額（その1676）

宮城県丸森町小斎日向 鹿島神社 明治9年(1876)

徳竹亜紀子，谷垣美保: 2022年度の算額調査，仙台高等専門学校名取キャンパス 研究紀要，第 59 号, p.9-47, 2022.
https://www.sendai-nct.ac.jp/natori-library/wp-content/uploads/2023/05/kiyo2023-2.pdf

キーワード：円5個，楕円
#Julia, #SymPy, #算額, #和算，#数学

楕円の中に大円 2 個，中円 2 個，小円 1 個を容れる。楕円の長径と短径が与えられたとき，大円の直径はいかほどか。

楕円の長半径，短半径，中心座標を a, b, (0, 0)
小円の半径と中心座標を r3, (0, 0)
中円の半径と中心座標を r2, (0, b - r2)
大円の半径と中心座標を r1, (r3 + r1, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms a, b, r1, r2, r3
eq1 = (r3 + r1)^2 + (b - r2)^2 - (r1 + r2)^2
eq2 = (r3 + r1)^2 - (a^2 - b^2)*(b^2 - r1^2)/b^2
eq3 = b - (2r2 + r3)
res = solve([eq1, eq2, eq3], (r1, r2, r3))[3];

# r1
res[1] |> println

    b*(a^4 - b^4 + b^2*sqrt(5*a^4 - 6*a^2*b^2 + b^4))/(a^2*(a^2 + 3*b^2))

# r2
res[2] |> println

    b*(3*a^2 + b^2 - sqrt(5*a^4 - 6*a^2*b^2 + b^4))/(2*(a^2 + 3*b^2))

# r3
res[3] |> println

    (-2*a^2*b + 2*b^3 + b*sqrt(5*a^4 - 6*a^2*b^2 + b^4))/(a^2 + 3*b^2)

楕円の長半径が 5，短半径が 3 のとき，大円の直径は 4.30030092052470 である。

2 * res[1](a => 5, b => 3).evalf() |> println

    4.30030092052470

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, r3) = (b*(a^4 - b^4 + b^2*sqrt(5*a^4 - 6*a^2*b^2 + b^4))/(a^2*(a^2 + 3*b^2)), b*(3*a^2 + b^2 - sqrt(5*a^4 - 6*a^2*b^2 + b^4))/(2*(a^2 + 3*b^2)), (-2*a^2*b + 2*b^3 + b*sqrt(5*a^4 - 6*a^2*b^2 + b^4))/(a^2 + 3*b^2))
    @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g;  r3 = %g\n", a, b, r1, r2, r3)
    plot()
    ellipse(0, 0, a, b)
    circle22(0, b - r2, r2, :blue)
    circle2(r3 + r1, 0, r1)
    circle(0, 0, r3, :magenta)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(r3, 0, "r3 ", :magenta, :right, :vcenter)
        point(r3 + r1, 0, "大円:r1,(r3+r1,0)", :red, :center, delta=-delta)
        point(0, b - r2, "中円:r2,(0,b-r2)", :blue, :center, delta=-delta)
        point(0, b, "b", :black, :center, :bottom, delta=delta)
        point(a, 0, " a", :black, :left, :vcenter)
    end
end;

draw(5, 3, true)

算額（その1675）

算額（その1675）

長野市若穂 清水寺観音堂
中村信弥「幻の算額」（204）
http://www.wasan.jp/maborosi/maborosi.html
キーワード：円2個，直角三角形，斜線
#Julia, #SymPy, #算額, #和算，#数学

直角三角形の中に斜線（矢）を描き，区画内に大円と小円を容れる。鈎と股の和が 42 寸，弦が 30 寸，大円の直径が 9 寸のとき，小円の直径はいかほどか。

直角三角形の 3 辺を「鈎」，「股」，「弦」とし，「中鈎」，「和」もそのまま変数とする。
大円の半径と中心座標を r1, (r1, y1)
小円の半径と中心座標を r2, (x2, r2)
斜線と斜辺（弦）の交点座標を (x0, y0)
とおき，以下の連立方程式を解く。

eq1 は「隔斜 n 円径の定理」である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms 鈎, 股, 弦, r1, r2, R, 中鈎, 和
R = (鈎 + 股 - sqrt(鈎^2 - 股^2))/2
R = (和 - 弦)/2  # 鈎 + 股 = 和
中鈎 = (和^2 - 弦^2)/(2*弦)
eq1 = (1 - 2R/中鈎) - (1 - 2r1/中鈎)*(1 - 2r2/中鈎);

ans_r2 = solve(eq1, r2)[1]
ans_r2 |> println

    (2*r1*和^2 - 2*r1*弦^2 - 和^3 + 和^2*弦 + 和*弦^2 - 弦^3)/(2*(4*r1*弦 - 和^2 + 弦^2))

ans_r2(和 => 42, 弦 => 30, r1 => 9/2) |> println

    4.00000000000000

小円の半径は 4 寸（直径は 8 寸）である。

算額の問への答えはここまでである。

以下では，図を描くためにパラメータを求める。

まず最初に，鈎，股を求める。

eq2 = 弦^2 - (鈎^2 + 股^2)
eq3 = 鈎 + 股 - 和;

res = solve([eq2, eq3], (鈎, 股))[2];

# 鈎
res[1] |> println

    和/2 - sqrt(-和^2 + 2*弦^2)/2

# 股
res[2]  |> println

    和/2 + sqrt(-和^2 + 2*弦^2)/2

次に，斜線と斜辺（弦）の交点座標 (x0, y0) と小円の中心の x 座標 x2 を求める。

@syms  x0::positive, y0::positive, y1::positive,
      x2::positive
eq4 = dist2(0, 0, x0, y0, r1, y1, r1)
eq5 = dist2(0, 0, x0, y0, x2, r2, r2)
eq6 = (鈎 - y0)/x0 - 鈎/股
eq7 = dist2(0, 鈎, 股, 0, x2, r2, r2);
res2 = solve([eq5, eq6, eq7], (x0, y0, x2))[1];

# x0
res2[1] |> println

    股*(-4*r2^3*股^2 - 4*r2^3*股*sqrt(股^2 + 鈎^2) - 4*r2^3*鈎^2 + 2*r2^2*股^2*鈎 + 2*r2^2*股*鈎*sqrt(股^2 + 鈎^2) + 4*r2^2*鈎^3 + 2*r2*股^2*鈎^2 - 股^2*鈎^3)/(鈎*(4*r2^2*股^2 + 4*r2^2*鈎^2 - 股^2*鈎^2))

# y0
res2[2] |> println

    2*r2*(2*r2^2*股^2 + 2*r2^2*股*sqrt(股^2 + 鈎^2) + 2*r2^2*鈎^2 + r2*股^2*鈎 - r2*股*鈎*sqrt(股^2 + 鈎^2) - 股^2*鈎^2)/(4*r2^2*股^2 + 4*r2^2*鈎^2 - 股^2*鈎^2)

# x2
res2[3] |> println

    (-r2*sqrt(股^2 + 鈎^2) - 股*(r2 - 鈎))/鈎

最後に 大円の中心の y 座標 y1 を求める。

eq8 = dist2(0, 鈎, 股, 0, r1, y1, r1)
ans_y1 = solve(eq8, y1)[1];

# y1
ans_y1 |> println

    (-r1*sqrt(股^2 + 鈎^2) - 鈎*(r1 - 股))/股

function draw(和, 弦, r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = (和 - 弦)/2  # 鈎 + 股 = 和
    中鈎 = (和^2 - 弦^2)/(2*弦)
    r2 = (2*r1*和^2 - 2*r1*弦^2 - 和^3 + 和^2*弦 + 和*弦^2 - 弦^3)/(2*(4*r1*弦 - 和^2 + 弦^2))
    (鈎, 股) = (和/2 - sqrt(-和^2 + 2*弦^2)/2, 和/2 + sqrt(-和^2 + 2*弦^2)/2)
    (x0, y0, x2) = (股*(-4*r2^3*股^2 - 4*r2^3*股*sqrt(股^2 + 鈎^2) - 4*r2^3*鈎^2 + 2*r2^2*股^2*鈎 + 2*r2^2*股*鈎*sqrt(股^2 + 鈎^2) + 4*r2^2*鈎^3 + 2*r2*股^2*鈎^2 - 股^2*鈎^3)/(鈎*(4*r2^2*股^2 + 4*r2^2*鈎^2 - 股^2*鈎^2)), 2*r2*(2*r2^2*股^2 + 2*r2^2*股*sqrt(股^2 + 鈎^2) + 2*r2^2*鈎^2 + r2*股^2*鈎 - r2*股*鈎*sqrt(股^2 + 鈎^2) - 股^2*鈎^2)/(4*r2^2*股^2 + 4*r2^2*鈎^2 - 股^2*鈎^2), (-r2*sqrt(股^2 + 鈎^2) - 股*(r2 - 鈎))/鈎)
    (x0, y0, y1) = (2*r1*(2*r1^2*股^2 + 2*r1^2*鈎^2 + 2*r1^2*鈎*sqrt(股^2 + 鈎^2) + r1*股*鈎^2 - r1*股*鈎*sqrt(股^2 + 鈎^2) - 股^2*鈎^2)/(4*r1^2*股^2 + 4*r1^2*鈎^2 - 股^2*鈎^2), 鈎*(-4*r1^3*股^2 - 4*r1^3*鈎^2 - 4*r1^3*鈎*sqrt(股^2 + 鈎^2) + 4*r1^2*股^3 + 2*r1^2*股*鈎^2 + 2*r1^2*股*鈎*sqrt(股^2 + 鈎^2) + 2*r1*股^2*鈎^2 - 股^3*鈎^2)/(股*(4*r1^2*股^2 + 4*r1^2*鈎^2 - 股^2*鈎^2)), (-r1*sqrt(股^2 + 鈎^2) - 鈎*(r1 - 股))/股)
    y1 = (-r1*sqrt(股^2 + 鈎^2) - 鈎*(r1 - 股))/股
    @printf("和 = %g;  弦 = %g;  大円の直径 = %g;  内接円の直径 = %g;  中鈎 = %g;  小円の直径 = %g\n",  和, 弦, 2r1, 2R, 中鈎, 2r2)
    @printf("鈎 = %g;  股 = %g;  x0 = %g;  y0 = %g;  y1 = %g;  x2 = %g\n", 鈎, 股, x0, y0, y1, x2)
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
    circle(r1, y1, r1)
    circle(x2, r2, r2, :blue)
    segment(0, 0, x0, y0, :magenta) 
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(r1, y1, "大円:r1,(r1,y1)", :red, :center, delta=-delta/2)
        point(x2, r2, "小円:r2,(x2,r2)", :blue, :center, delta=-delta/2)
        point(x0, y0, "(x0,y0)", :magenta, :left, :bottom, delta=delta/2)
        point(0, 鈎, " 鈎", :green, :left, :bottom, delta=delta/2)
        point(股, 0, " 股", :green, :left, :bottom, delta=delta/2)
    end

end;

draw(42, 30, 9/2, true)

算額（その1674）

算額（その1674）

福島県田村郡三春町平沢担橋 諏訪神社 大正15年(1926)
キーワード：円5個，外円，斜線2本
#Julia, #SymPy, #算額, #和算，#数学

外円と甲円 2 個が輪違いになっている。上の甲円の中心を通り下の甲円に接する斜線を 2 本引いてできる領域に乙円 2 個を容れる。黒積（扇形の面積から乙円の面積を差し引いた面積）を求めよ。

外円の半径と中心座標を R, (0, r1-R)
甲円の半径と中心座標を r1, (0, r1), (0, -r1)
乙円の半径と中心座標を r2, (0, r2), (0, 2r1 - r2)
斜線と外円の交点座標を (x1, y1)
斜線と甲円の交点座標を (x2, y2)
とおき，以下の連立方程式を解く。

いくつかのパラメータは自明である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy

@syms R::positive, r1::positive, r2::positive, θ::positive
R = 3r1/2
eq = 2r2 - (r1 - r2);

ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println

    3*r2

乙円の半径 r2 は甲円の半径 r1 の 1/3 である。

黒積は，甲円の面積の 1/6 から乙円の面積を差し引いたものである。

黒積 = pi*(3r2)^2/6 - pi*r2^2
黒積 |> println

    pi*r2^2/2

黒積は乙円の面積 πr2^2 の 1/2 である。
乙円の面積が 2 歩のとき，黒積は 1 歩である。

function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 3r1/2
    r2 = r1/3
    θ = 30
    @printf("R = %g;  r1 = %g;  r2 = %g;  θ = %g\n", R, r1, r2, θ)
    plot()
    circle22(0, r1, r1)
    circle(0, r1 - R, R, :blue)
    circle(0, 2r1 - r2, r2, :green)
    circlef(0, r1, r1, :black, beginangle=240, endangle=300)
    plot!([0, r1*sind(θ), -r1*sind(θ), 0], [r1,  r1 - r1*cosd(θ), r1 - r1*cosd(θ), r1], color=:black, lw=0.5, seriestype=:shape, fillcolor=colorfillcolor=:black, )
    circlef(0, r2, r2 , :white)
    circle(0, r2, r2, :green)
    x1 = R*sind(2θ)
    y1 = R*cosd(2θ) + (r1 - 2R)
    x2 = -r1*sind(θ)
    y2 = 2r1 + x2*tand(θ/2)
    segment(x1, y1, x2, y2, :magenta)
    segment(-x1, y1, -x2, y2, :magenta) 
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, 2r1, "2r1", :red, :center, :bottom, delta=delta)
        point(0, r1, "r1", :red, :center, :bottom, delta=2delta)
        point(0, 2r1 - r2, "乙円:r2\n(0,2r1-r2)", :green, :center, :bottom, delta=delta)
        point(0, r2, "乙円:r2\n(0,r2)", :green, :center, :bottom, delta=delta/2)
        point(0, -r1, "甲円,r1,(0,-r1)", :red, :center, :bottom, delta=delta)
        point(0, r1 - R, "外円,R,(0,r1-R)", :blue, :center, :bottom, delta=delta)
        point(x1, y1, " (x1, y1)", :magenta)
        point(x2, y2, "(x2,y2)  ", :magenta, :right, :vcenter)
    end

end;

draw(1/2, true)