算額（その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
