算額（その267）

算額（その267）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（64）
長野県下高井郡木島平村往郷 水穂神社 寛政12年(1800)

甲，乙，丙の 3 個の正方形がある。甲，乙の面積の和を A，乙と丙の面積の和を B とする。このとき丙の一辺の長さを求めよ。ただし正方形の一辺の差は等しいとする。

甲，乙，丙の 3 個の正方形の一辺の長さを，甲，乙，丙とする。

include("julia-source.txt");

using SymPy
@syms 甲::positive, 乙::positive, 丙::positive, A::positive, B::positive;

eq1 = 甲^2 + 乙^2 - A
eq2 = 乙^2 + 丙^2 - B
eq3 = (甲 - 乙) - (乙 - 丙)
res = solve([eq1, eq2, eq3], (甲, 乙, 丙));

以下の 2 組の解が得られる。

res[1][1] |> println
res[1][2] |> println
res[1][3] |> println

    sqrt(2)*(-A + B - 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B - sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B))
    -(B/2 + sqrt(-A^2 + 6*A*B - B^2)/4)*sqrt(-2*A + 14*B - 2*sqrt(-A^2 + 6*A*B - B^2))/(A - 5*B)
    sqrt(-A/8 + 7*B/8 - sqrt(-A^2 + 6*A*B - B^2)/8)

res[2][1] |> println
res[2][1] |> println
res[2][3] |> println

    sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B))
    sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B))
    sqrt(-A/8 + 7*B/8 + sqrt(-A^2 + 6*A*B - B^2)/8)

それぞれを使って，実際に甲，乙，丙を与えて A, B を求め，数式により正しく甲，乙，丙が返ってくるか確かめる。

function func1(甲, 乙, 丙)  # 1 組目の解
    A = 甲^2 + 乙^2
    B = 乙^2 + 丙^2
    return (sqrt(2)*(-A + B - 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B - sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B)),
-(B/2 + sqrt(-A^2 + 6*A*B - B^2)/4)*sqrt(-2*A + 14*B - 2*sqrt(-A^2 + 6*A*B - B^2))/(A - 5*B),
sqrt(-A/8 + 7*B/8 - sqrt(-A^2 + 6*A*B - B^2)/8))
end;

function func2(甲, 乙, 丙)  # 2 組目の解
    A = 甲^2 + 乙^2
    B = 乙^2 + 丙^2
    return (sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B)),
sqrt(2)*(-A + B + 2*sqrt(-A^2 + 6*A*B - B^2))*sqrt(-A + 7*B + sqrt(-A^2 + 6*A*B - B^2))/(4*(A - 5*B)),
sqrt(-A/8 + 7*B/8 + sqrt(-A^2 + 6*A*B - B^2)/8))
end;

func1(8, 6, 4) |> println
func1(9.34, 8.084, 6.828) |> println

    (8.000000000000002, 6.0, 4.0)
    (9.340000000000003, 8.084000000000001, 6.828000000000001)

func2(8, 6, 4) |> println
func2(9.34, 8.084, 6.828) |> println

    (-9.899494936611665, -9.899494936611665, 7.0710678118654755)
    (-12.32062855539441, -12.32062855539441, 10.544376321053798)

1 組目の解 func1 のほうが適切解である。

算額（その266）

算額（その266）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（171）
長野県下高井郡木島平村上木島 天然寺 弘化2年(1845)

図のように，正三角形の内側に 3 個の正方形がある。それぞれの正方形の 2 つの頂点は正三角形の辺の上にあり，他の 2 個の頂点は他の正方形の頂点と一致している。正三角形の一辺が 8 寸のとき，正方形の一辺の長さを求めよ。

正三角形の一辺の長さを a，正方形の一辺の長さを x とすると，⊿ABC は一辺の長さが x の正三角形，∠ OAD は 30°，OB = a/2 ゆえ 2(x + √3/2x) = a。x について解くと x = (2 - √3)a。a に 8 を代入すると x = 2.14359353944898。

include("julia-source.txt");

using SymPy
@syms x::positive, a::positive;

eq1 = 2(x + sqrt(Sym(3))/2 * x) - a
res = solve(eq1, x)[1]
res |> println

   -sqrt(3)*a + 2*a

res(a => 8).evalf()

   2.14359353944898

正方形の一辺の長さは 2.143593 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 8
   x = (2 - √3)a
   @printf("正三角形の一辺の長さ = %.6f;  正方形の一辺の長さ = %.6f\n", a, x)
   plot([a/2, 0, -a/2, a/2], [0, a√3/2, 0, 0], color=:blue, lw=0.5)
   xy = x .* [ # 図形の頂点座標，一筆書き
           0                        sind(30)                 # 1
           cosd(60)                 sind(30)+sind(60)        # 2
           cosd(30) + cosd(60)      sind(60)                 # 3
           cosd(30)                 0                        # 4
           0                        sind(30)                 # 5
           -cosd(30)                0                        # 6
           -cosd(30) - cosd(60)     sind(60)                 # 7
           -cosd(60)                sind(30) + sind(60)      # 8
           -cosd(60)                sind(30) + sind(60) + 1  # 9
           cosd(60)                 sind(30) + sind(60) + 1  # 10
           cosd(60)                 sind(30) + sind(60)      # 11
           -cosd(60)                sind(30) + sind(60)      # 12
           0                        sind(30)                 # 13
       ]
   plot!(xy[:, 1], xy[:, 2], color=:red, lw=0.5)
   if more
       point(0, 0, " O", :red, :left, :bottom)
       point(x*cosd(30), 0, " A", :red, :left, :bottom)
       point(x*cosd(30) + x, 0, " B", :red, :left, :bottom)
       point(x*(cosd(30)+cosd(60)), x*sind(60), " C", :red, :left, :bottom) 
       point(0, x*sind(30), "  D", :red, :left, :bottom)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その265）

算額（その265）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（91）
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)

直角三角形（鉤股弦）内に円と正方形がある。三角形の面積は 294 平方寸　，鉤の3乗と股の3上の和に円の直径（径）と正方形の一辺の長さ（方面）の差を掛けると 62426 となる（数字の掛け算で次元は無視）。このとき，円の直径を求めよ。

include("julia-source.txt");

using SymPy
@syms 鉤::positive, 股::positive, 弦::positive,
   径::positive, 方面::positive;

eq1 = (鉤 + 股 - 弦) - 径
eq2 = (鉤 * 股 / 2) - 294;
eq3 = (鉤^3 + 股^3)*(径 - 方面) - 62426
eq4 = (鉤^2 + 股^2) - 弦^2
eq5 = (鉤 - 方面)/方面 - 方面/(股 - 方面)
res = solve([eq1, eq2, eq3, eq4, eq5], (鉤, 股, 弦, 径, 方面))

   2-element Vector{NTuple{5, Sym}}:
    (21, 28, 35, 14, 12)
    (28, 21, 35, 14, 12)

鉤 < 股 なので (21, 28, 35, 14, 12) が適切である。

鉤 = 21;  股 = 28;  弦 = 35;  径 = 14;  方面 = 12

円の直径は 14 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 径, 方面) = (21, 28, 35, 14, 12)
   r = 径/2
   @printf("鉤 = %.0f;  股 = %.0f;  弦 = %.0f;  径 = %.0f;  方面 = %.0f\n", 鉤, 股, 弦, 径, 方面)
   plot([0, 股, 0, 0], [0, 0, 鉤, 0], color=:black, lw=0.5)
   plot!([0, 方面, 方面, 0, 0], [0, 0, 方面, 方面, 0], color=:black, lw=0.5)
   circle(径/2, 径/2, 径/2)
   if more
       point(r, r, "半径=$(Int(r))")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その264）

算額（その264）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（135）
長野県長野市元善町 善光寺 天保3年(1832)

藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html

正方形内に，甲，乙，丙，丁の 4 円が入っている。甲円か乙円の径が与えられているとき，正方形の一辺の長さを求めよ。

正方形の一辺の長さを x とする。
甲円の半径，中心座標を r1, (x - r1, x - r1)
乙円の半径，中心座標を r2, (r2, r2)
丙円の半径，中心座標を r3, (x - r3, r3)
丁円の半径，中心座標を r4, (x4, r4) とする。

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

甲円の径を与えるか乙円の径を与えるかは本質的なものではない。どちらにどのような値を与えても，得られる解（図）は相似である。つまり，各円の径及び正方形の一辺の長さは，どれか一つの数値の定数倍である。

以下でまず 甲円の半径を 1 とした場合の解を求める。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
   r4::positive, x4::positive, x::positive;

r1 = 1
eq1 = 2(x - r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (x - r1 - r3)^2 - (r1 + r3)^2
eq3 = (x - r1 - x4)^2 + (x - r1 - r4)^2 - (r1 + r4)^2
eq4 = (r2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq5 = (x - r3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2

res = solve([eq1, eq2, eq3, eq4, eq5], (r2, r3, r4, x4, x));

7 組の解が得られるが，そのうちの 4 番目のものが適切である。

res[4]

   ((-32362318*sqrt(420*sqrt(2) + 594) - 45767229*sqrt(210*sqrt(2) + 297) + 1115429952 + 788728083*sqrt(2))/(-560916*sqrt(6)*sqrt(70*sqrt(2) + 99) - 793255*sqrt(3)*sqrt(70*sqrt(2) + 99) + 19333056 + 13670535*sqrt(2)), (-64396059*sqrt(210*sqrt(2) + 297) - 45534890*sqrt(420*sqrt(2) + 594) + 1569448152 + 1109767431*sqrt(2))/(-560916*sqrt(6)*sqrt(70*sqrt(2) + 99) - 793255*sqrt(3)*sqrt(70*sqrt(2) + 99) + 19333056 + 13670535*sqrt(2)), (-68*sqrt(420*sqrt(2) + 594) - 96*sqrt(210*sqrt(2) + 297) + 1656*sqrt(2) + 2342)/(-sqrt(3)*sqrt(70*sqrt(2) + 99) + 9*sqrt(2) + 13), (-118630*sqrt(420*sqrt(2) + 594) - 167768*sqrt(210*sqrt(2) + 297) + 4088811 + 2891226*sqrt(2))/(-447*sqrt(6)*sqrt(70*sqrt(2) + 99) - 632*sqrt(3)*sqrt(70*sqrt(2) + 99) + 15405 + 10893*sqrt(2)), -2*sqrt(210*sqrt(2) + 297) + 18*sqrt(2) + 26)

そのまま数値解として表示すると以下のようになる。

((-32362318*sqrt(420*sqrt(2) + 594) - 45767229*sqrt(210*sqrt(2) + 297) + 1115429952 + 788728083*sqrt(2))/(-560916*sqrt(6)*sqrt(70*sqrt(2) + 99) - 793255*sqrt(3)*sqrt(70*sqrt(2) + 99) + 19333056 + 13670535*sqrt(2)), (-64396059*sqrt(210*sqrt(2) + 297) - 45534890*sqrt(420*sqrt(2) + 594) + 1569448152 + 1109767431*sqrt(2))/(-560916*sqrt(6)*sqrt(70*sqrt(2) + 99) - 793255*sqrt(3)*sqrt(70*sqrt(2) + 99) + 19333056 + 13670535*sqrt(2)), (-68*sqrt(420*sqrt(2) + 594) - 96*sqrt(210*sqrt(2) + 297) + 1656*sqrt(2) + 2342)/(-sqrt(3)*sqrt(70*sqrt(2) + 99) + 9*sqrt(2) + 13), (-118630*sqrt(420*sqrt(2) + 594) - 167768*sqrt(210*sqrt(2) + 297) + 4088811 + 2891226*sqrt(2))/(-447*sqrt(6)*sqrt(70*sqrt(2) + 99) - 632*sqrt(3)*sqrt(70*sqrt(2) + 99) + 15405 + 10893*sqrt(2)), -2*sqrt(210*sqrt(2) + 297) + 18*sqrt(2) + 26)

   (0.5887908616460656, 0.4184622124074972, 0.36337064107190264, 1.513883688904384, 2.712235388919648)

数式のまま取り扱い，簡約化すると以下のようになる。

r2 = (-32362318*sqrt(420*sqrt(Sym(2)) + 594) - 45767229*sqrt(210*sqrt(Sym(2)) + 297) + 1115429952 + 788728083*sqrt(Sym(2)))/(-560916*sqrt(Sym(6))*sqrt(70*sqrt(Sym(2)) + 99) - 793255*sqrt(Sym(3))*sqrt(70*sqrt(Sym(2)) + 99) + 19333056 + 13670535*sqrt(Sym(2)))
r3 = (-64396059*sqrt(210*sqrt(Sym(2)) + 297) - 45534890*sqrt(420*sqrt(Sym(2)) + 594) + 1569448152 + 1109767431*sqrt(Sym(2)))/(-560916*sqrt(Sym(6))*sqrt(70*sqrt(Sym(2)) + 99) - 793255*sqrt(Sym(3))*sqrt(70*sqrt(Sym(2)) + 99) + 19333056 + 13670535*sqrt(Sym(2)))
r4= (-68*sqrt(420*sqrt(Sym(2)) + 594) - 96*sqrt(210*sqrt(Sym(2)) + 297) + 1656*sqrt(Sym(2)) + 2342)/(-sqrt(Sym(3))*sqrt(70*sqrt(Sym(2)) + 99) + 9*sqrt(Sym(2)) + 13)
x4 = (-118630*sqrt(420*sqrt(Sym(2)) + 594) - 167768*sqrt(210*sqrt(Sym(2)) + 297) + 4088811 + 2891226*sqrt(Sym(2)))/(-447*sqrt(Sym(6))*sqrt(70*sqrt(Sym(2)) + 99) - 632*sqrt(Sym(3))*sqrt(70*sqrt(Sym(2)) + 99) + 15405 + 10893*sqrt(Sym(2)))
x  = -2*sqrt(210*sqrt(Sym(2)) + 297) + 18*sqrt(Sym(2)) + 26;

apart(r2, x) |> println
apart(r3, x) |> println
apart(r4, x) |> println
apart(x4, x) |> println
apart(x , x) |> sympy.sqrtdenest |> println

   -6*sqrt(6) - 8*sqrt(3) + 10*sqrt(2) + 15
   -12*sqrt(3) - 8*sqrt(6) + 14*sqrt(2) + 21
   -27*sqrt(3) - 19*sqrt(6) + 33*sqrt(2) + 47
   -27*sqrt(6) - 38*sqrt(3) + 47*sqrt(2) + 67
   -10*sqrt(6) - 14*sqrt(3) + 18*sqrt(2) + 26

次に乙円の半径を 1 にしたときの解を示す。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
   r4::positive, x4::positive, x::positive;

r2 = 1
eq1 = 2(x - r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (x - r1 - r3)^2 - (r1 + r3)^2
eq3 = (x - r1 - x4)^2 + (x - r1 - r4)^2 - (r1 + r4)^2
eq4 = (r2 - x4)^2 + (r2 - r4)^2 - (r2 + r4)^2
eq5 = (x - r3 - x4)^2 + (r3 - r4)^2 - (r3 + r4)^2

res2 = solve([eq1, eq2, eq3, eq4, eq5], (r1, r3, r4, x4, x));

同じく 7 組の解が得られるが，2 番目のものが適切である。

r2 = 1
res2[2]

   ((-105731*sqrt(2) - 25282*sqrt(18 - 12*sqrt(2)) + 35765*sqrt(9 - 6*sqrt(2)) + 149534)/(-1831*sqrt(2) - 436*sqrt(6)*sqrt(3 - 2*sqrt(2)) + 623*sqrt(3)*sqrt(3 - 2*sqrt(2)) + 2594), (-252671*sqrt(2) - 60424*sqrt(18 - 12*sqrt(2)) + 85457*sqrt(9 - 6*sqrt(2)) + 357334)/(-1831*sqrt(2) - 436*sqrt(6)*sqrt(3 - 2*sqrt(2)) + 623*sqrt(3)*sqrt(3 - 2*sqrt(2)) + 2594), (-58789*sqrt(2) - 14058*sqrt(18 - 12*sqrt(2)) + 19885*sqrt(9 - 6*sqrt(2)) + 83143)/(-1831*sqrt(2) - 436*sqrt(6)*sqrt(3 - 2*sqrt(2)) + 623*sqrt(3)*sqrt(3 - 2*sqrt(2)) + 2594), (-610*sqrt(2) - 144*sqrt(18 - 12*sqrt(2)) + 210*sqrt(9 - 6*sqrt(2)) + 867)/(-49*sqrt(2) - 11*sqrt(6)*sqrt(3 - 2*sqrt(2)) + 18*sqrt(3)*sqrt(3 - 2*sqrt(2)) + 71), -2*sqrt(2) + 2*sqrt(3)*sqrt(3 - 2*sqrt(2)) + 6)

r1 = (-105731*sqrt(Sym(2)) - 25282*sqrt(18 - 12*sqrt(Sym(2))) + 35765*sqrt(9 - 6*sqrt(Sym(2))) + 149534)/(-1831*sqrt(Sym(2)) - 436*sqrt(Sym(6))*sqrt(3 - 2*sqrt(Sym(2))) + 623*sqrt(Sym(3))*sqrt(3 - 2*sqrt(Sym(2))) + 2594)
r3 = (-252671*sqrt(Sym(2)) - 60424*sqrt(18 - 12*sqrt(Sym(2))) + 85457*sqrt(9 - 6*sqrt(Sym(2))) + 357334)/(-1831*sqrt(Sym(2)) - 436*sqrt(Sym(6))*sqrt(3 - 2*sqrt(Sym(2))) + 623*sqrt(Sym(3))*sqrt(3 - 2*sqrt(Sym(2))) + 2594)
r4 = (-58789*sqrt(Sym(2)) - 14058*sqrt(18 - 12*sqrt(Sym(2))) + 19885*sqrt(9 - 6*sqrt(Sym(2))) + 83143)/(-1831*sqrt(Sym(2)) - 436*sqrt(Sym(6))*sqrt(3 - 2*sqrt(Sym(2))) + 623*sqrt(Sym(3))*sqrt(3 - 2*sqrt(Sym(2))) + 2594)
x4 = (-610*sqrt(Sym(2)) - 144*sqrt(18 - 12*sqrt(Sym(2))) + 210*sqrt(9 - 6*sqrt(Sym(2))) + 867)/(-49*sqrt(Sym(2)) - 11*sqrt(Sym(6))*sqrt(3 - 2*sqrt(Sym(2))) + 18*sqrt(Sym(3))*sqrt(3 - 2*sqrt(Sym(2))) + 71)
x  = -2*sqrt(Sym(2)) + 2*sqrt(Sym(3))*sqrt(3 - 2*sqrt(Sym(2))) + 6;

apart(r1, x) |> println
apart(r3, x) |> println
apart(r4, x) |> println
apart(x4, x) |> println
apart(x , x) |> sympy.sqrtdenest |> println

   -10*sqrt(2) - 8*sqrt(3) + 6*sqrt(6) + 15
   -20*sqrt(3) - 24*sqrt(2) + 14*sqrt(6) + 35
   -5*sqrt(3) - 5*sqrt(2) + 3*sqrt(6) + 9
   -2*sqrt(3) - sqrt(2) + sqrt(6) + 5
   -2*sqrt(3) - 2*sqrt(2) + 2*sqrt(6) + 6

数式としてみると r1 = 1 としたときと違うように見えるが比を取ってみると同じであることがわかる。

図を描くプログラムも，相似であることを利用すれば単純になるが，以下のプログラムでは r1 = 1 とする場合と r2 = 1 とする場合を分けて書いた。なお，r1, r2 が 1 以外の場合には単に定数を掛けてやればよい。

using Plots

function draw(more=false; which=:r1)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   if which == :r1
       r1 = 1
       r2 = -6*sqrt(6) - 8*sqrt(3) + 10*sqrt(2) + 15
       r3 = -12*sqrt(3) - 8*sqrt(6) + 14*sqrt(2) + 21
       r4 = -27*sqrt(3) - 19*sqrt(6) + 33*sqrt(2) + 47
       x4 = -27*sqrt(6) - 38*sqrt(3) + 47*sqrt(2) + 67
       x  = -10*sqrt(6) - 14*sqrt(3) + 18*sqrt(2) + 26
   else
       r2 = 1
       r1 = -10*sqrt(2) - 8*sqrt(3) + 6*sqrt(6) + 15
       r3 = -20*sqrt(3) - 24*sqrt(2) + 14*sqrt(6) + 35
       r4 = -5*sqrt(3) - 5*sqrt(2) + 3*sqrt(6) + 9
       x4 = -2*sqrt(3) - sqrt(2) + sqrt(6) + 5
       x  = -2*sqrt(3) - 2*sqrt(2) + 2*sqrt(6) + 6
   end
   @printf("r1 = %.5f;  r2 = %.5f;  r3 = %.5f;  r4 = %.5f;  x4 = %.5f;  x = %.5f\n", r1, r2, r3, r4, x4, x)
   plot([0, x, x, 0, 0], [0, 0, x, x, 0], color=:black, lw=0.5)
   circle(x - r1, x - r1, r1)
   circle(r2, r2, r2)
   circle(x - r3, r3, r3)
   circle(x4, r4, r4)
   if more
       point(x - r1, x - r1, " 甲円:r1\n (x-r1,x-r1")
       point(r2, r2, " 乙円:r2\n (r2,r2)")
       point(x - r3, r3, " 丙円:r3\n (x-r3,r3)")
       point(x4, r4, " 丁円:r4\n (x4,r4)")
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

甲円の半径が 1 の場合

draw(true)

   r1 = 1.00000;  r2 = 0.58879;  r3 = 0.41846;  r4 = 0.36337;  x4 = 1.51388;  x = 2.71224

乙円の半径が 1 の場合

draw(true; which=:r2)

   r1 = 1.69840;  r2 = 1.00000;  r3 = 0.71071;  r4 = 0.61715;  x4 = 2.57117;  x = 4.60645

 

算額（その263）

算額（その263）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（94）
長野県長野市三輪 美和神社 文化10年(1813)

図のような，横 4丈2尺，縦  2丈4尺の長方形がある。長方形の内側に図のように糸を張ったところ，糸の長さが 15丈になった。このとき，等弦の本数（折れ線の数）を求めよ。

弦の本数を n とする。各弦は直角三角形の斜辺である。直角を挟む辺の長さは片方は 24尺であり，もう一方は 42/n 尺であるので，弦の長さは sqrt((42/n)^2 + 24^2) である。これが n 本あり，その合計が 15丈 = 150尺なので，以下の方程式を解けばよい。

include("julia-source.txt");

using SymPy

@syms n::positive;

eq = n * sqrt((42/n)^2 + 24^2) - 150
solve(eq)[1] |> println

   6

6 本である。

using Printf
using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   plot([0, 42, 42, 0, 0], [0, 0, 24, 24, 0], linecolor=:black, linewidth=0.5)
   plot!([0, 5, 10, 15, 20, 25], [24, 0, 24, 0, 24, 0], linestyle=[:solid, :solid, :solid, :solid, :solid, :solid], linecolor=:blue, linewidth=1)
   plot!([25, 27.5], [0, 12], linecolor=:blue, linestyle=:dash)
   if more == true
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5, xlims=(-2, 60))
   else
       plot!(showaxis=false)
   end
end;

算額（その262）

算額（その262）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（71）
長野県長野市 久保寺観音堂 享和3年(1803)

鉤股弦（直角三角形）の中に円と正方形が入っている。
方の一辺と円の直径の積が 672 平方寸，鉤股弦の周と中鉤の積が 5644.8 平方寸であるとき，方面（正方形の一辺）を求めよ。
注：この 2 つの数値を，小数でなくするために 10 倍するのは，正しくない。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms  鉤::positive, 股::positive, 弦::positive,
      中鉤::positive, 長弦::positive, 短弦::positive,
      方面::positive, 径::positive;

eq1 = 鉤^2 + 股^2 - 弦^2
eq2 = 短弦^2 + 中鉤^2 - 鉤^2
eq3 = 長弦^2 + 中鉤^2 - 股^2
eq4 = 弦 - 長弦 - 短弦
eq5 = 方面 * 径 - 672.0
eq6 = (鉤 + 股 + 弦) * 中鉤 - 5644.8
eq7 = 鉤 + 股 - 弦 - 径
eq8 = 方面 * (方面 + (股 - 方面)/2 + (鉤 - 方面)/2) - 鉤 * 股 / 2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])  #, (鉤, 股, 弦, 中鉤, 長弦, 短弦, 方面, 径))

   2-element Vector{Dict{Any, Any}}:
    Dict(股 => 56.0000000000000, 弦 => 70.0000000000000, 径 => 28.0000000000000, 短弦 => 25.2000000000000, 長弦 => 44.8000000000000, 方面 => 24.0000000000000, 鉤 => 42.0000000000000, 中鉤 => 33.6000000000000)
    Dict(股 => 42.0000000000000, 弦 => 70.0000000000000, 径 => 28.0000000000000, 短弦 => 44.8000000000000, 長弦 => 25.2000000000000, 方面 => 24.0000000000000, 鉤 => 56.0000000000000, 中鉤 => 33.6000000000000)

2 つの解が求まる。最初の解は「鉤 < 股」で，もう一方は裏返しに反転して回転する「鉤 > 股」のもので，本質的に同じものである。

なお，「solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8], (鉤, 股, 弦, 中鉤, 長弦, 短弦, 方面, 径))」とすると解が求まらない（有限の時間内にできなさそう）。

for i = 1:2
   println("***** i = $i *****")
   println((res[i][鉤], res[i][股], res[i][弦], res[i][中鉤], res[i][長弦], res[i][短弦], res[i][方面], res[i][径]))
end

   ***** i = 1 *****
   (42.0000000000000, 56.0000000000000, 70.0000000000000, 33.6000000000000, 44.8000000000000, 25.2000000000000, 24.0000000000000, 28.0000000000000)
   ***** i = 2 *****
   (56.0000000000000, 42.0000000000000, 70.0000000000000, 33.6000000000000, 25.2000000000000, 44.8000000000000, 24.0000000000000, 28.0000000000000)

   鉤 = 42.00;  股 = 56.00;  弦 = 70.00;  中鉤 = 33.60;  長弦 = 44.80;  短弦 = 25.20; 方面 = 24.00;  径 = 28.00

答え：方面は 24 寸である。

using Printf
using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 中鉤, 長弦, 短弦, 方面, 径) = (42.0000000000000, 56.0000000000000, 70.0000000000000, 33.6000000000000, 44.8000000000000, 25.2000000000000, 24.0000000000000, 28.0000000000000)
   # (鉤, 股, 弦, 中鉤, 長弦, 短弦, 方面, 径) = (56.0000000000000, 42.0000000000000, 70.0000000000000, 33.6000000000000, 25.2000000000000, 44.8000000000000, 24.0000000000000, 28.0000000000000)
   @printf("鉤 = %.2f;  股 = %.2f;  弦 = %.2f;  中鉤 = %.2f;  長弦 = %.2f;  短弦 = %.2f; 方面 = %.2f;  径 = %.2f\n", 鉤, 股, 弦, 中鉤, 長弦, 短弦, 方面, 径)
   plot([股, 股, 0, 股, 長弦*股/弦], [ 0, 鉤, 0, 0, 長弦*鉤/弦], linecolor=:black, linewidth=0.5)
   circle(股 - 径/2, 径/2, 径/2)
   plot!([股 - 方面, 股 - 方面, 股], [0, 方面, 方面], color=:black, lw=0.5)
   if more == true
       point(股 - 径/2, 径/2, "\n円:径/2\n(股-径/2,径/2)", :red, :center, :top)
       point(股 - 方面, 方面, "(股-方面,方面)", :red, :right, :bottom)
       point(股, 鉤/2, " 鉤", mark=false)
       point(股/2, 0, "股  ", :green, :right, :bottom, mark=false)
       point((股 + 長弦*股/弦)/2, 長弦*鉤/弦/2, "  中鉤", mark=false)
       point(長弦*股/弦/2, 長弦*鉤/弦/2, "長弦   \n", :green, :right, mark=false)
       point((股 + 長弦*股/弦)/2, (鉤 + 長弦*鉤/弦)/2, "短弦\n", :green, :right, :bottom, mark=false)
       point(股/2, 鉤/2, "弦 = 長弦 + 短弦   ", :green, :right, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その261）

算額（その261）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（71）
長野県長野市 久保寺観音堂 享和3年(1803)

十七 群馬県高崎市八幡町 文化7年(1810)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

キーワード：円2個，直角三角形，中鈎
#Julia, #SymPy, #算額, #和算

鈎股弦（直角三角形）の中鈎を隔てて，大円，小円が入っている。大円，小円の直径が 68.8 寸，51.6 寸であるとき，鈎，股，弦を求めよ。

鈎，股，弦，中鈎，短弦，長弦 をそのまま「鈎」，「股」，「弦」，「中鈎」，「短弦」，「長弦」
大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (股 - r2, y2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

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

eq1 = 鈎^2 + 股^2 - 弦^2

# eq2 = 短弦^2 + 中鈎^2 - 鈎^2
# eq3 = 長弦^2 + 中鈎^2 - 股^2
eq2 = 短弦/鈎 - 鈎/弦
eq3 = 中鈎/股 - 鈎/弦

eq4 = 弦 - 長弦 - 短弦

# eq5 = dist2(長弦*股/弦, 長弦*鈎/弦, 股, 0, x1, r1, r1)
# eq6 = dist2(長弦*股/弦, 長弦*鈎/弦, 股, 0, 股 - r2, y2, r2)
eq5 = 長弦 + 中鈎 - 股 - 2r1
eq6 = 短弦 + 中鈎 - 鈎 - 2r2

eq7 = dist2(0, 0, 股, 鈎, x1, r1, r1)
eq8 = dist2(0, 0, 股, 鈎, 股 - r2, y2, r2);

SymPy の性能上，一度に解けないので，まず eq1 〜 eq6 を解いて 鈎, 股, 弦, 中鈎, 長弦, 短弦 を求める。

res = solve([eq1, eq2, eq3, eq4, eq5, eq6], (鈎, 股, 弦, 中鈎, 長弦, 短弦))[2]  # 2 of 2

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

2 組の解が得られるが，2 番目のものが適解である。弦を表す式は以下のように簡約化できる。

ans_鈎 = res[1] |> factor
ans_鈎 |> println

   (r1^2 + r2^2)*(r1^2 + r1*r2 + r1*sqrt(r1^2 + r2^2) + 2*r2^2 + 2*r2*sqrt(r1^2 + r2^2))/(r1*(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2)))

ans_股 = res[2] |> factor
ans_股 |> println

   (r1^2 + r2^2)*(r1^2 + r1*r2 + r1*sqrt(r1^2 + r2^2) + 2*r2^2 + 2*r2*sqrt(r1^2 + r2^2))/(r2*(r1^2 + r2^2 + r2*sqrt(r1^2 + r2^2)))

ans_弦 = res[3] |> factor
ans_弦 |> println

   (r1^2 + r2^2)*(r1 + r2 + sqrt(r1^2 + r2^2))/(r1*r2)

大円，小円の直径が 68.8 寸，51.6 寸であるとき，鈎 = 129 寸，股 = 172 寸，弦 = 215 寸である。

ans_鈎(r1 => 68.8/2, r2 => 51.6/2).evalf() |> println
ans_股(r1 => 68.8/2, r2 => 51.6/2).evalf() |> println
ans_弦(r1 => 68.8/2, r2 => 51.6/2).evalf() |> println

   129.000000000000
   172.000000000000
   215.000000000000

鈎，股の解が得られたので，eq7, eq8 の鈎，股に代入し，連立方程式を解いて x1, y2 を求める。

eq17 = eq7(鈎 => res[1], 股 => res[2]) |> simplify |> numerator;
eq18 = eq8(鈎 => res[1], 股 => res[2]) |> simplify |> numerator;
res2 = solve([eq17, eq18], (x1, y2))[3]  # 3 of 4

   (r1^2/r2 + r1*sqrt(r1^2 + r2^2)/r2, (r1^4 + 2*r1^3*r2 + r1^3*sqrt(r1^2 + r2^2) + 4*r1^2*r2^2 + 3*r1^2*r2*sqrt(r1^2 + r2^2) + 2*r1*r2^3 + 2*r1*r2^2*sqrt(r1^2 + r2^2) + 2*r2^4 + 2*r2^3*sqrt(r1^2 + r2^2) - r2*sqrt(r1^6 + 9*r1^4*r2^2 + 4*r1^4*r2*sqrt(r1^2 + r2^2) + 16*r1^2*r2^4 + 12*r1^2*r2^3*sqrt(r1^2 + r2^2) + 8*r2^6 + 8*r2^5*sqrt(r1^2 + r2^2)))/(r1^3 + 2*r1*r2^2 + 2*r1*r2*sqrt(r1^2 + r2^2)))

4 組の解が得られるが 3 番目のものが適解である。

x
res2[1](r1 => 68.8/2, r2 => 51.6/2) |>  println

   103.200000000000

y
res2[2](r1 => 68.8/2, r2 => 51.6/2) |>  println

   77.4000000000000

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = sqrt(r1^2 + r2^2)
   (鈎, 股, 弦, 中鈎, 長弦, 短弦) = ((r1^4 + r1^3*r2 + r1^3*t + 3*r1^2*r2^2 + 2*r1^2*r2*t + r1*r2^3 + r1*r2^2*t + 2*r2^4 + 2*r2^3*t)/(r1*(r1^2 + r2^2 + r2*t)), (r1^4 + r1^3*r2 + r1^3*t + 3*r1^2*r2^2 + 2*r1^2*r2*t + r1*r2^3 + r1*r2^2*t + 2*r2^4 + 2*r2^3*t)/(r2*(r1^2 + r2^2 + r2*t)), (r1^2*(r1 + t) + r1*r2*(r1 + r2) + r2^2*(r2 + t))/(r1*r2), (r1^3 + 2*r1^2*r2 + r1^2*t + r1*r2^2 + r1*r2*t + 2*r2^3 + 2*r2^2*t)/(r1^2 + r2^2 + r2*t), r1*(r1 + r2 + t)/r2, r2*(r1 + r2 + t)/r1)
   (x1, y2) = (r1^2/r2 + r1*t/r2, (r1^4 + 2*r1^3*r2 + r1^3*t + 4*r1^2*r2^2 + 3*r1^2*r2*t + 2*r1*r2^3 + 2*r1*r2^2*t + 2*r2^4 + 2*r2^3*t - r2*sqrt(r1^6 + 9*r1^4*r2^2 + 4*r1^4*r2*t + 16*r1^2*r2^4 + 12*r1^2*r2^3*t + 8*r2^6 + 8*r2^5*t))/(r1^3 + 2*r1*r2^2 + 2*r1*r2*t))
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  中鈎 = %g;  長弦 = %g;  短弦 = %g; x1 = %g;  y2 = %g\n", 鈎, 股, 弦, 中鈎, 長弦, 短弦, x1, y2)
   # println("2r1 + 4r2 = $(2r1 + 4r2)")
   plot([股, 股, 0, 股, 長弦*股/弦], [ 0, 鈎, 0, 0, 長弦*鈎/弦], linecolor=:black, linewidth=0.5)
   circle(x1, r1, r1)
   circle(股 - r2, y2, r2, :blue)
   if more == true
       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(x1, r1, "大円:r1\n(x1,r1)", :red, :center, delta=-delta)
       point(股 - r2, y2, "小円:r2\n(股-r2,y2)", :blue, :center, delta=-delta)
       point(股, 鈎/2, " 鈎", mark=false)
       point(股/2, 0, "股  ", :green, :right, :bottom, mark=false)
       point((股 + 長弦*股/弦)/2, 長弦*鈎/弦/2, "  中鈎", mark=false)
       point(長弦*股/弦/2, 長弦*鈎/弦/2, "長弦   \n", :green, :right, mark=false)
       point((股 + 長弦*股/弦)/2, (鈎 + 長弦*鈎/弦)/2, "短弦\n", :green, :right, :bottom, mark=false)
       point(股/2, 鈎/2, "弦 = 長弦 + 短弦   ", :green, :right, mark=false)
       xlims!(-3delta, 股 + 5delta)
   end
end;

draw(68.8/2, 51.6/2, true)

 

算額（その260）

算額（その260）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（70）
長野県長野市 久保寺観音堂 享和3年(1803)

圭（二等辺三角形）内に，大円 1 個，小円 2 個が入っている。
中鉤（圭の高さ）が 36 寸，大円と小円 2 個の径の和が 44寸のとき，大円の径を求めよ。

大円と小円の半径をそれぞれ r1, r2，圭の底辺の長さを 2x とする。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms l::positive, x::positive, r1::positive, r2::positive;

l = sqrt(36^2 + x^2)
eq1 = 2r1 + 4r2 - 44  # 大円の径と小円の径の 2 倍（2 個ぶん）の和が 44 寸
eq2 = (l + x)r1 - 36x  # 面積との関係
eq3 = (36 + l + x)r2 - 36x  # 面積との関係

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

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (10, 6, 15)

大円の径 = 20 寸，小円の径 = 12 寸，底辺の長さ = 30 寸である。

using Printf
using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x) = res[1]
   @printf("r1 = %.5f;  r2 = %.5f;  x = %.5f\n", r1, r2, x)
   println("2r1 + 4r2 = $(2r1 + 4r2)")
   plot([x, 0, -x, x], [0, 36, 0, 0], linecolor=:black, linewidth=0.5)
   circle(0, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(-r2, r2, r2, :blue)
   if more == true
       point(0, 36, " 36")
       point(x, 0, " x", :green, :left, :bottom)
       point(0, r1, " 大円:r1")
       point(r2, r2, " 小円:r1\n (r2,r2)")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

A = 中鉤, B = 径の和 として A, B を変数のまま eq1, eq2, eq3 を解く。4 組の解が得られるが，4 番目のものが適切である。

using SymPy

@syms l::positive, x::positive, r1::positive, r2::positive;
@syms A, B

l = sqrt(A^2 + x^2)
eq1 = 2r1 + 4r2 - B  # 大円の径と小円の径の 2 倍（2 個ぶん）の和が 44 寸
eq2 = (l + x)r1 - A*x
eq3 = r1*(x - r2) - r2*x

res = solve([eq1, eq2, eq3], (r1, r2, x))[4]

   ((-3*A*(4*A - B) + 3*sqrt(A*(4*A - B))*(2*A - B) + (3*A - B)*(B + sqrt((B^2*(3*A - B)^2 - 2*B*(3*A - B)*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B)) + 9*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B))^2)/(3*A - B)^2)))/(4*(3*A - B)), (3*A*(4*A - B) + 3*sqrt(A*(4*A - B))*(-2*A + B) + (3*A - B)*(B - sqrt((B^2*(3*A - B)^2 - 2*B*(3*A - B)*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B)) + 9*(A*(4*A - B) - sqrt(A*(4*A - B))*(2*A - B))^2)/(3*A - B)^2)))/(8*(3*A - B)), (A*(4*A - B) + sqrt(A*(4*A - B))*(-2*A + B))/(2*(3*A - B)))

同じ項が何箇所にも現れるのでそれをまとめると文字数としては少なくなるが，複雑な式であることに変わりはない。

A = 36; B = 44
t1 = A*(4*A - B)
t2 = (2*A - B)
t3 = (3*A - B)
t4 = sqrt(t1)*t2
t5 = (t1 - t4);

((-3t1 + 3t4 + t3*(B + sqrt((B^2*t3^2 - 2B*t3*t5 + 9t5^2)/t3^2)))/(4t3),
(3t1 + 3sqrt(t1)*(-2A + B) + t3*(B - sqrt((B^2*t3^2 - 2B*t3*t5 + 9t5^2)/t3^2)))/(8t3),
(t1 + sqrt(t1)*(-2A + B))/(2t3))

   (10.0, 6.0, 15.0)

 

算額（その259）

算額（その259）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（70）
長野県長野市 久保寺観音堂 享和3年(1803)

鉤股弦（直角三角形）において，
   (1) 鉤^3 + 股^3 = 223.894125
   (2) 中鉤^3 + 短弦^3 = 48.361131
のとき，鉤，股，弦を求めよ。

「中鉤」，「短弦」は図の如し。
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms 鉤::positive, 股::positive, 弦::positive, 中鉤::positive, 短弦::positive, 長弦::positive;

eq1 = 鉤^3 + 股^3 - 223.894125
eq2 = 中鉤^3 + 短弦^3 - 48.361131
eq3 = 鉤^2 + 股^2 - 弦^2
eq4 = 短弦^2 + 中鉤^2 - 鉤^2
eq5 = 長弦^2 + 中鉤^2 - 股^2
eq6 = 長弦 + 短弦 - 弦

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

   1-element Vector{NTuple{6, Sym}}:
    (4.05000000000000, 5.40000000000000, 6.75000000000000, 3.24000000000000, 2.43000000000000, 4.32000000000000)

   鉤 = 4.05000;  股 = 5.40000;  弦 = 6.75000;  中鉤 = 3.24000;  短弦 = 2.43000;  長弦 = 4.32000

鉤 = 4.05, 股 = 5.40, 弦 = 6.75寸 である。

using Printf
using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, 弦, 中鉤, 短弦, 長弦) = res[1]
   @printf("鉤 = %.5f;  股 = %.5f;  弦 = %.5f;  中鉤 = %.5f;  短弦 = %.5f;  長弦 = %.5f\n", 鉤, 股, 弦, 中鉤, 短弦, 長弦)
   plot([股, 股, 0, 股, 長弦*股/弦], [ 0, 鉤, 0, 0, 長弦*鉤/弦], linecolor=:black, linewidth=0.5)
   if more == true
       point(股, 鉤/2, " 鉤", mark=false)
       point(股/2, 0, " 股", :green, :center, :bottom, mark=false)
       point((股 + 長弦*股/弦)/2, 長弦*鉤/弦/2, "  中鉤", mark=false)
       point(長弦*股/弦/2, 長弦*鉤/弦/2, " 長弦", mark=false)
       point((股 + 長弦*股/弦)/2, (鉤 + 長弦*鉤/弦)/2, " 短弦", mark=false)
       point(股/2, 鉤/2, "弦 = 長弦 + 短弦   ", :green, :right, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その258）

算額（その258）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（70）
長野県長野市 久保寺観音堂 享和3年(1803)

菱形の中に大円 1 個と小円 2 個が入っている。菱形の辺の長さの和と菱形の横の対角線の長さの積が 160 平方寸，菱形の辺の長さの和と大円の直径の積が 96 平方寸のとき，小円の直径を求めよ。

菱形の中心を原点とし，長い方と短い方の対角線の長さを 2x, 2y とおく。
大円の半径，中心座標を r1, (0, 0)，
小円の半径，中心座標を r2, (r2, 0) とおく。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, x::positive, y::positive, sinθ::positive;

xy = sqrt(x^2 + y^2)
sinθ = y / xy
eq1 = 2x * 4xy - 160
eq2 = 2r1 * 4xy - 96
eq3 = x*sinθ - r1
eq4 = (x - r2)sinθ - r2;

# res = solve([eq1, eq2, eq3, eq4], (r1, r2, x, y))

なぜか solve() では解けないので nlsolve() を用いる。

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, x, y) = u
   return [
       8*x*sqrt(x^2 + y^2) - 160,  # eq1
       8*r1*sqrt(x^2 + y^2) - 96,  # eq2
       -r1 + x*y/sqrt(x^2 + y^2),  # eq3
       -r2 + y*(-r2 + x)/sqrt(x^2 + y^2),  # eq4
   ]
end;

iniv = [1.0, 1, 1, 1]
res = nls(H, ini=iniv);
println(res);

   ([2.4, 1.5, 4.0, 2.9999999999999996], true)

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

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, x, y) = res[1]
   @printf("r1 = %.5f;  r2 = %.5f;  x = %.5f;  y = %.5f\n", r1, r2, x, y)
   println("小円の直径 = $(2r2)")
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], linecolor=:black, linewidth=0.5)
   circlef(0, 0, r1, :red)
   circlef(r2, 0, r2, :yellow)
   circlef(-r2, 0, r2, :yellow)
   circle(0, 0, r1, :black)
   circle(r2, 0, r2, :black)
   circle(-r2, 0, r2, :black)
   if more == true
       point(0, y, " y", :black, :left, :bottom)
       point(x, 0, " y", :black, :left, :bottom)
       point(0.2, 1.8, "大円:r1", :black, mark=false)
       point(r2, 0, "小円:r2\n", :black, :center, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その257）

算額（その257）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額2（112）
長野県長野市信更町田野口 清水神社 文政11年(1828)

正方形内に斜線で区切られた部分に円と正方形がある。円の径と内側の正方形の一辺の長さが等しい。外側の正方形の一辺の長さが 100 寸のとき，円の径を求めよ。

外側の正方形の対角座標を (0, 0), (x, x) とする。
円の半径と中心座標を r, (x - r, x - r)
内側の正方形の対角座標を (0, 0), (r, r) とする。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms x::positive, r::positive;

x = 100
eq = 2((x - r) - 2r)^2 - r^2
res = solve(eq, r)[1]
2res |> println

   2*x*(6 - sqrt(2))/17

最初の解が適切である。
外側の正方形の一辺の長さに (12 - 2√2)/17 を掛けて円の直径を得る。

円の直径は 53寸9分5厘あまりである。

100 * (12 - 2√2)/17

   53.950428677963586

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = 100
   r = x * (6 - √2)/17
   println("円の直径 = $(2r)")
   plot([0, x, x, 0, 0], [0, 0, x, x, 0], linecolor=:black, linewidth=0.5)
   rect(0, 0, 2r, 2r, :aliceblue)
   circlef(x - r, x - r, r, :pink)
   segment(4r - x, x, x, 4r - x)
   if more == true
       point(x, x, "(x,x)", :black, :right, :top)
       point(2r, 2r, " (2r,2r)", :black, :right, :top)
       point(x - r, x - r, "(x-r,x-r)", :red, :center)
       point(x, 4r-x, "(x,4r-x) ", :black, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その256）

算額（その256）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額4（186)
長野県長野市北尾張部 連開庵 嘉永6年(1853)

正方形内に四分円，半円，天円，地円，人円が入っており，互いに内接・外接している。
地円，人円の径がそれぞれ 1寸4分4厘，1 寸のとき天円の径を求めよ。

四分円の半径を r0 とする。半円の半径は r0/2 である。
天円の半径と中心座標を r1, (x1, r0 - r1)
地円の半径と中心座標を r2, (x2, r2) とする。r2 = 1.44 である。
人円の半径と中心座標を r3, (r3, y3) とする。r3 = 1 である。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

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

r2 = 144//200
r3 = 100//200
eq1 = (r0 - x1)^2 + r1^2 - (r0 - r1)^2
eq2 = (r0 - x1)^2 + (r0 - r1 - r0//2)^2 - (r1 + r0//2)^2
eq3 = (r0 - r3)^2 + (r0 - y3)^2 - (r0 + r3)^2
eq4 = (r0 - x2)^2 + (r0 - r2)^2 - (r0 + r2)^2
eq5 = (x2 - r3)^2 + (y3 - r2)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, x1, x2, y3))

   4-element Vector{NTuple{5, Sym}}:
    (3*sqrt(2)/5 + 181/100 + 11*sqrt(120*sqrt(2) + 241)/100, 3*sqrt(2)/20 + 181/400 + 11*sqrt(120*sqrt(2) + 241)/400, 3*sqrt(2)/5 + 181/100 + 11*sqrt(120*sqrt(2) + 241)/100 + sqrt(660*sqrt(2)*sqrt(120*sqrt(2) + 241) + 18120*sqrt(2) + 34561 + 1991*sqrt(120*sqrt(2) + 241))/100, -sqrt(120*sqrt(2) + 241)/100 + 49/100 + 3*sqrt(2)/5, sqrt(3*sqrt(2)/250 + 241/10000) + 71/100 + 3*sqrt(2)/5)
    (3*sqrt(2)/5 + 181/100 + 11*sqrt(120*sqrt(2) + 241)/100, 3*sqrt(2)/20 + 181/400 + 11*sqrt(120*sqrt(2) + 241)/400, -sqrt(660*sqrt(2)*sqrt(120*sqrt(2) + 241) + 18120*sqrt(2) + 34561 + 1991*sqrt(120*sqrt(2) + 241))/100 + 3*sqrt(2)/5 + 181/100 + 11*sqrt(120*sqrt(2) + 241)/100, -sqrt(120*sqrt(2) + 241)/100 + 49/100 + 3*sqrt(2)/5, sqrt(3*sqrt(2)/250 + 241/10000) + 71/100 + 3*sqrt(2)/5)
    (-11*sqrt(120*sqrt(2) + 241)/100 + 3*sqrt(2)/5 + 181/100, -11*sqrt(120*sqrt(2) + 241)/400 + 3*sqrt(2)/20 + 181/400, -11*sqrt(120*sqrt(2) + 241)/100 - sqrt(-1991*sqrt(120*sqrt(2) + 241) - 660*sqrt(2)*sqrt(120*sqrt(2) + 241) + 18120*sqrt(2) + 34561)/100 + 3*sqrt(2)/5 + 181/100, sqrt(120*sqrt(2) + 241)/100 + 49/100 + 3*sqrt(2)/5, -sqrt(3*sqrt(2)/250 + 241/10000) + 71/100 + 3*sqrt(2)/5)
    (-11*sqrt(120*sqrt(2) + 241)/100 + 3*sqrt(2)/5 + 181/100, -11*sqrt(120*sqrt(2) + 241)/400 + 3*sqrt(2)/20 + 181/400, -11*sqrt(120*sqrt(2) + 241)/100 + sqrt(-1991*sqrt(120*sqrt(2) + 241) - 660*sqrt(2)*sqrt(120*sqrt(2) + 241) + 18120*sqrt(2) + 34561)/100 + 3*sqrt(2)/5 + 181/100, sqrt(120*sqrt(2) + 241)/100 + 49/100 + 3*sqrt(2)/5, -sqrt(3*sqrt(2)/250 + 241/10000) + 71/100 + 3*sqrt(2)/5)

二番目の解が適切である。

天円の直径は 2寸4分4厘あまりである。

2res[2][2].evalf()

   2.44388710953848

   r0 = 4.88777;  r1 = 1.22194;  x1 = 1.43160;  x2 = 1.13587;  y3 = 1.76119
   天円の直径 = 2.44388710953848

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, x1, x2, y3) = res[2]
   @printf("r0 = %.5f;  r1 = %.5f;  x1 = %.5f;  x2 = %.5f;  y3 = %.5f\n", r0, r1, x1, x2, y3)
   println("天円の直径 = $(2r1.evalf())")
   plot([0, r0, r0, 0, 0], [0, 0, r0, r0, 0], linecolor=:black, linewidth=0.5)
   circle(r0, r0, r0, beginangle=180, endangle=270)
   circle(r0, r0//2, r0//2, :blue, beginangle=90, endangle=270)
   circle(x1, r0 - r1, r1, :orange)
   circle(x2, r2, r2, :magenta)
   circle(r3, y3, r3, :green)
   if more == true
       point(r0, r0, "(r0,r0)", :black, :right, :bottom)
       point(r0, r0/2, "(r0,r0/2)", :blue, :right)
       point(x1, r0 - r1, " 天円:r1\n (x1,r0-r1)", :orange)
       point(x2, r2, " 地円:r2\n (x2,r2)", :magenta)
       point(r3, y3, "\n 人円:r3\n (r3,y3)", :green, :center, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その255）

算額（その255）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県外(285)
東京都あきる野市 二宮神社 寛政6年(1794)

勾股弦（直角三角形）に 5 個の等円（直径が同じ円）が弦（直角に対する変）及び鉤（直角を挟む短い方の辺），股（同じく長い方の辺）に接している。内接する全円の径が 44 寸のとき，等円の径を求めよ。

全円の半径を r1 = 22 寸 とする。
等円の半径を r2 とする。

以下の方程式をとき，まず勾，股，弦の長さを求める。

include("julia-source.txt");

using SymPy

# 全円半径 r1 = 22 寸、
# 等円半径 r2
@syms 鉤, 股, 弦, r1::positive, r2::positive,  x::positive, y::positive
弦 = 110
eq1 = 鉤 + 股 - 弦 - 44  # 鉤股弦と内接円の径の関係
eq2 = 鉤^2 + 股^2 - 弦^2  # ピタゴラスの定理
(鉤, 股) = solve([eq1, eq2], (鉤, 股))[1]

   (66, 88)

eq3 = distance(0, 鉤, 股, 0, r2, y) - r2^2
eq4 = distance(0, 鉤, 股, 0, x, r2) - r2^2
eq5= (x - r2)^2 + (y - r2)^2 -(8r2)^2;

res = solve([eq3, eq4, eq5], (r2, x, y))

   4-element Vector{Tuple{Sym, Sym, Sym}}:
    (110/13, 814/13, 638/13)
    (660/53, 4884/53, 3828/53)
    (-1540/191 + 880*sqrt(15)/191, 21428/191 - 2640*sqrt(15)/191, 440*sqrt(15)/191 + 11836/191)
    (-2310/491 + 1980*sqrt(15)/491, 660*sqrt(15)/491 + 42438/491, 37026/491 - 3960*sqrt(15)/491)

names = ["等円半径", "内側の股の x 座標", "内側の鉤の y 座標"]
for i = 1:4
   println(i)
   for j = 1:3
       println("$(names[j]) = $(res[i][j]) = $(res[i][j].evalf())")
   end
end

   1
   等円半径 = 110/13 = 8.46153846153846
   内側の股の x 座標 = 814/13 = 62.6153846153846
   内側の鉤の y 座標 = 638/13 = 49.0769230769231
   2
   等円半径 = 660/53 = 12.4528301886792
   内側の股の x 座標 = 4884/53 = 92.1509433962264
   内側の鉤の y 座標 = 3828/53 = 72.2264150943396
   3
   等円半径 = -1540/191 + 880*sqrt(15)/191 = 9.78128452702894
   内側の股の x 座標 = 21428/191 - 2640*sqrt(15)/191 = 58.6561464189132
   内側の鉤の y 座標 = 440*sqrt(15)/191 + 11836/191 = 70.8906422635145
   4
   等円半径 = -2310/491 + 1980*sqrt(15)/491 = 10.9134562637285
   内側の股の x 座標 = 660*sqrt(15)/491 + 42438/491 = 91.6378187545762
   内側の鉤の y 座標 = 37026/491 - 3960*sqrt(15)/491 = 44.1730874725430

一番目の解が適切である。

   (r2, x, y) = (110/13, 814/13, 638/13)

   (8.461538461538462, 62.61538461538461, 49.07692307692308)

   等円の半径 = 8.46154;  等円の直径 = 16.92308;  x = 62.61538;  y = 49.07692

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鉤, 股, r2, x, y) = (66, 88, 110/13, 814/13, 638/13)
   @printf("等円の半径 = %.5f;  等円の直径 = %.5f;  x = %.5f;  y = %.5f\n", r2, 2r2, x, y)
   plot([0, 88, 0, 0], [0, 0, 66, 0], color=:black, lw=0.5)
   circle(22, 22, 22, :red)
   for i = 0:4
       circle(r2*(1 + 2i*4/5) , y - 2i*r2*3/5, r2, :blue)
   end
   if more == true
       plot!([r2, x, r2, r2], [r2, r2, y, r2], color=:black, lw=0.5)
       circle(r2, r2, r2, :blue)
       point(股, 0, "88", :black, :left, :bottom)
       point(0, 鉤, "  66", :black)
       point(22, 22, " 全円\n(r1,r1)", :red) 
       point(r2, y, " 等円 (r2,y)\n", :blue, :center, :bottom)
       point(x, r2, "\n (x,r2)", :blue, :center, :top)
       point(r2, r2, " (r2,r2)", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その249）

算額（その249）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額3（230)
長野県諏訪市下諏訪 諏訪大社下社 慶応4年(1868)

長方形内に甲，乙，丙，丁の円が 1 個ずつ入っている。丙円の径が 5寸0分6厘2毛5糸，丁円の径が 1寸4分4厘のとき，甲円の径はいかほどか。

注：以下の説明図は算額の図に近くなるように，r3 = 2, r4 = 1 としたものである。

長方形の長辺，短辺の長さをそれぞれ x，y とする。左下隅角を原点とする。
甲円の半径，中心座標を r1, (x − r1, r1) とする。y = 2r1 である。
乙円の半径，中心座標を r2, (r2, r2) とする。
丙円の半径，中心座標を r3, (x3, y - r3) とする。r3 = 50625 / 2
丁円の半径，中心座標を r4, (x4, r4) とする。r4 = 14400 / 2

以下の連立方程式を nlsolve() で解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, x3::positive,
   r4::positive, x4::positive, x::positive;

この段階では r3, r4 も未知数とする

y = 2r1
eq1 = (x3 - r2)^2 + (y - r3 - r2)^2 - (r2 + r3)^2      # 乙円と丙円が外接
eq2 = (x - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2      # 甲円と乙円が外接
eq3 = (x4 - r2)^2 + (r2 - r4)^2 - (r2 + r4)^2          # 乙円と丁円が外接
eq4 = (x - r1 - x3)^2 + (y - r3 - r1)^2 - (r1 + r3)^2  # 甲円と丙円が外接
eq5 = (x - r1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2;     # 甲円と丁円が外接

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

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

   (-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2,  # eq1
   (r1 - r2)^2 - (r1 + r2)^2 + (-r1 - r2 + x)^2,  # eq2
   (-r2 + x4)^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq3
   (r1 - r3)^2 - (r1 + r3)^2 + (-r1 + x - x3)^2,  # eq4
   (r1 - r4)^2 - (r1 + r4)^2 + (-r1 + x - x4)^2,  # eq5

算額の通り，r3 = 50625 / 2, r4 = 14400 / 2 とすると，算額の図とは似ても似つかない図が得られる。しかし，求まる r1, r2 は算額の通りである。

算額の多くは，実際に図を描くためのすべてのパラメータを計算しないので，実際とは異なる図を掲示してしまう。
算額では r1 = r4*(sqrt(2sqrt(r3/r4) + 1/4) + 1/2)^2
算額の図は r3/r4≒2 の場合のもの（前述の説明図）であるが，問では 5.0625/1.4400≒3.5 なので違って当然なわけである。

r3 = 50625 / 2
r4 = 14400 / 2
f(r3, r4) = r4*(sqrt(2sqrt(r3/r4) + 1/4) + 1/2)^2
f(50625 / 2, 14400 / 2), f(50625 / 2, 50625 / 4)

   (45000.0, 64331.362230694234)

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)
   (x, r1, r2, x3, x4) = u
   return [
       (-r2 + x3)^2 - (r2 + r3)^2 + (2*r1 - r2 - r3)^2,  # eq1
       (r1 - r2)^2 - (r1 + r2)^2 + (-r1 - r2 + x)^2,  # eq2
       (-r2 + x4)^2 + (r2 - r4)^2 - (r2 + r4)^2,  # eq3
       (r1 - r3)^2 - (r1 + r3)^2 + (-r1 + x - x3)^2,  # eq4
       (r1 - r4)^2 - (r1 + r4)^2 + (-r1 + x - x4)^2,  # eq5
   ]
end;

(r3, r4) = (50625 // 2, 14400 // 2)
(r3, r4) = (2, 1)
iniv = [big"200000.0", 89000, 50000, 50000, 60000]
res = nls(H, ini=iniv);
println(res);
   (BigFloat[125000.0000000000000000000000000000000078311250578101097708795336203908490939767, 45000.00000000000000000000000000000000240863585814089291883850249994466664153376, 20000.00000000000000000000000000000000145859703738035599114200719600180134501038, 12500.00000000000000000000000000000000319845103367027695329127846269509477692267, 44000.00000000000000000000000000000000423640766899408968775137416801857989402613], true)

names = ["x", "r1", "r2", "x3", "x4"]
for i in 1:length(names)
   println("$(names[i]) = $(Float64(res[1][i]))")
end

   x = 125000.0
   r1 = 45000.0
   r2 = 20000.0
   x3 = 12500.0
   x4 = 44000.0

甲円の半径は 45000 なので，元の単位での直径は 9 寸。

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   # res, r3, r4 は　グローバル変数
   (x, r1, r2, x3, x4) = res[1]
   @printf("r3 = %.5f;  r4 = %.5f\n", r3, r4)
   @printf("x = %.5f;  r1 = %.5f;  r2 = %.5f;  x3 = %.5f;  x4 = %.5f\n", x, r1, r2, x3, x4)
   y = 2r1
   plot([0, x, x, 0, 0], [0, 0, y, y, 0], color=:black, lw=0.5)
   circle(x - r1, r1, r1)
   circle(r2, r2, r2, :blue)
   circle(x3, y - r3, r3, :green)
   circle(x4, r4, r4, :magenta)
   if more == true
       point(x − r1, r1, " 甲円:r1\n (x−r1,r1)\n", :red, :center, :bottom)
       point(r2, r2, "\n乙円:r2\n (r2,r2)\n", :blue, :center, :bottom)
       point(x3, y - r3, "丙円:r3\n(x3,y-r3)\n", :green, :center, :bottom)
       point(x4, r4, "   丁円:r4\n      (x4,r4)", :magenta, :left, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       point(0, y, " y", :black, :left, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;