算額（その202）

算額（その202）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（105）
長野県伊那市東春近 春近神社 文政4年(1821)

直線上に大円，中円，小円がある。それぞれの円は直線，垂直線，斜線に接している。3 個の円の径の和が 54 寸，小円を含む鉤股弦の和が 30 寸のとき，小円の径を求めよ。

大円，中円，小円の半径をそれぞれ r1, r2, r3 とおく。また，鉤股の長さをそれぞれ x, y とおく。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms r1::positive, r2::positive, r3::positive, x::positive, y::positive;

eq1 = distance(0, y, x, 0, r1, r1) - r1^2
eq2 = distance(0, y, x, 0, -r2, r2) - r2^2
eq3 = distance(0, y, x, 0, r3, r3) - r3^2
eq4 = 2(r1 + r2 + r3) - 54
eq5 = x + y + sqrt(x^2 + y^2) - 30;
# solve([eq1, eq2, eq3, eq4, eq5]);

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

   -r1^2 + (r1 - x*(r1*x - r1*y + y^2)/(x^2 + y^2))^2 + (r1 - y*(-r1*x + r1*y + x^2)/(x^2 + y^2))^2,
   -r2^2 + (-r2 - x*(-r2*x - r2*y + y^2)/(x^2 + y^2))^2 + (r2 - y*(r2*x + r2*y + x^2)/(x^2 + y^2))^2,
   -r3^2 + (r3 - x*(r3*x - r3*y + y^2)/(x^2 + y^2))^2 + (r3 - y*(-r3*x + r3*y + x^2)/(x^2 + y^2))^2,
   2*r1 + 2*r2 + 2*r3 - 54,
   x + y + sqrt(x^2 + y^2) - 30,

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r1, r2, r3, x, y) = u
   return [
       -r1^2 + (r1 - x*(r1*x - r1*y + y^2)/(x^2 + y^2))^2 + (r1 - y*(-r1*x + r1*y + x^2)/(x^2 + y^2))^2,
       -r2^2 + (-r2 - x*(-r2*x - r2*y + y^2)/(x^2 + y^2))^2 + (r2 - y*(r2*x + r2*y + x^2)/(x^2 + y^2))^2,
       -r3^2 + (r3 - x*(r3*x - r3*y + y^2)/(x^2 + y^2))^2 + (r3 - y*(-r3*x + r3*y + x^2)/(x^2 + y^2))^2,
       2*r1 + 2*r2 + 2*r3 - 54,
       x + y + sqrt(x^2 + y^2) - 30,
   ]
end;

iniv = [big"20.0", 8, 3, 12, 8]
res = nls(H, ini=iniv)

   (BigFloat[14.99999999999999999984666619116852447418727617700909974481543102428426836661943, 9.999999999999999998427131201068448021956451753434173348206627211441946563294198, 2.000000000000000001726202607763027503856272069556726906977939919601503399470133, 5.000000000000000001419534990100076452230824423574926396608810570956077428130958, 11.99999999999999999932144013561540591403353792037404330533967686003722496194269], true)

(r1, r2, r3, x, y) = res[1]

   14.99999999999999999984666619116852447418727617700909974481543102428426836661943
   9.999999999999999998427131201068448021956451753434173348206627211441946563294198
   2.000000000000000001726202607763027503856272069556726906977939919601503399470133
   5.000000000000000001419534990100076452230824423574926396608810570956077428130958
   11.99999999999999999932144013561540591403353792037404330533967686003722496194269

r1 = 15.000;  r2 = 10.000;  r3 = 2.000;  鉤 = x = 5.000;  股 = y = 12.000;  弦 = 13.000

小円の半径は 2 である。元の単位でいえば，直径 4 寸。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function abline(x0, y0, slope, minx, maxx)
   f(x0) = slope * x0 + y0
   # println("slope = $slope $x0, $(f(x0)), $y0, $(f(y0))")
   segment(minx, f(minx), maxx, f(maxx))
end

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, x, y) = res[1]
   @printf("r1 = %.3f;  r2 = %.3f;  r3 = %.3f;  鉤 = x = %.3f;  股 = y = %.3f;  弦 = %.3f\n",
       r1, r2, r3, x, y, sqrt(x^2 + y^2))
   plot(ylims=(-2, 31))
   circle(r1, r1, r1)
   circle(-r2, r2, r2, :blue)
   circle(r3, r3, r3, :green)
   abline(0, y, -y/x, -0.4r2, 2r1)
   hline!([0], color=:black, lw=0.5)
   vline!([0], color=:black, lw=0.5)
   if more == true
       point(r1, r1, "大円:(r1, r1)", :red, :center)
       point(-r2, r2, "中円:(-r2, r2)", :blue, :center)
       point(r3, r3, "  小円:(r3, r3)", :green, :left, :bottom)
       point(x, 0, "  x", :green, :left)
       point(0, y, "y ", :blue, :right)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その201）

算額（その201）

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

直線と弧に挟まれて元，利，亨，貞の4円がある。元円の径は 9 寸，亨円の径は 8 寸である。貞円の径はいかほどか。

弧は半径が R，中心座標 (0, R) の円の一部とする。元，亨の円の中心の x 座標を x1, x3 とする。利円の半径を r2, 中心座標を (x2, y2) とする。貞円の半径を r4, 中心座標を (x4, r4) とする。
以下のように方程式を立て，解く。ただし，solve() では無理なようなので，nlsolver() を用いる。

using SymPy

@syms R::positive, x1::positive, x2::positive, y2::positive, r2::positive, x3::positive, x4::positive, r4::positive;

eq1 = x4^2 + (R - r4)^2 - (r4 + R)^2
eq2 = x2^2 + (R - y2)^2 - (r2 + R)^2
eq3 = x1^2 + (R - 9)^2 - (9 + R)^2
eq4 = (x4 - x2)^2 + (r4 - y2)^2 - (r4 + r2)^2
eq5 = (x2 - x1)^2 + (y2 - 9)^2 - (r2 + 9)^2
eq6 = (x3 - x1)^2 + (9 -8)^2 - (9 +  8)^2
eq7 = (x4 - x3)^2 + (r4 - 8)^2 - (r4 + 8)^2
eq8 = (x3 - x2)^2 + (y2 - 8)^2 - (r2 + 8)^2;
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])

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

   x4^2 + (R - r4)^2 - (R + r4)^2,
   x2^2 - (R + r2)^2 + (R - y2)^2,
   x1^2 + (R - 9)^2 - (R + 9)^2,
   -(r2 + r4)^2 + (r4 - y2)^2 + (-x2 + x4)^2,
   -(r2 + 9)^2 + (-x1 + x2)^2 + (y2 - 9)^2,
   (-x1 + x3)^2 - 288,
   (r4 - 8)^2 - (r4 + 8)^2 + (-x3 + x4)^2,
   -(r2 + 8)^2 + (-x2 + x3)^2 + (y2 - 8)^2,

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (R, x1, x2, y2, r2, x3, x4, r4) = u
   return [
x4^2 + (R - r4)^2 - (R + r4)^2,
x2^2 - (R + r2)^2 + (R - y2)^2,
x1^2 + (R - 9)^2 - (R + 9)^2,
-(r2 + r4)^2 + (r4 - y2)^2 + (-x2 + x4)^2,
-(r2 + 9)^2 + (-x1 + x2)^2 + (y2 - 9)^2,
(-x1 + x3)^2 - 288,
(r4 - 8)^2 - (r4 + 8)^2 + (-x3 + x4)^2,
-(r2 + 8)^2 + (-x2 + x3)^2 + (y2 - 8)^2,
   ]
end;

iniv = [big"100.0", 60, 70, 20, 6, 78, 105, 28]
res = nls(H, ini=iniv)

   (BigFloat[72.00000000000000000825964442488647097716026341883355173366655519563804792349667, 50.91168824543142175984557839685392111956016953456497828078367372174802581654967, 61.09402589451770611111762160903542988124497556281032888237971124359306058389206, 21.59999999999999999909987194204532742690715360782740479871527281754573245830236, 7.199999999999999999201109949817982282417360698313820741808391250977189279801073, 67.88225099390856234546584308737029806240472447241742343849815521828381897192468, 101.8233764908628435172626133333358387527768508723559120900733596797539734806522, 36.00000000000000000036157749114823774743490558123508432409135127356873277697285], true)

   R = 72.000;  x1 = 50.912;  x2 = 61.094;  y2 = 21.600;  r2 = 7.200;  x3 = 67.882;  x4 = 101.823;  r4 = 36.000

貞円の半径は 36。元の単位では 丁円の径は 36 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x1, x2, y2, r2, x3, x4, r4) = res[1]
   #(R, x1, x2, y2, r2, x3, x4, r4) = [100, 60, 70, 20, 6, 78, 105, 28]
   @printf("R = %.3f;  x1 = %.3f;  x2 = %.3f;  y2 = %.3f;  r2 = %.3f;  x3 = %.3f;  x4 = %.3f;  r4 = %.3f\n",
   R, x1, x2, y2, r2, x3, x4, r4)
   plot(xlims=(-2, 150), ylims=(-2, 75))
   circle(0, R, R)
   circle(x1, 9, 9, :blue)
   circle(x2, y2, r2, :green)
   circle(x3, 8, 8, :orange)
   circle(x4, r4, r4, :purple)
   if more == true
       point(0, R, " R")
       point(x1, 9, "元:(x1,9)", :blue, :right)
       point(x2, y2, "利:(x2,y2,r2)", :blue, :right)
       point(x3, 8, "亨:(x3,8)")
       point(x4, r4, "貞:(x4,r4)")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   end
end;

算額（その200）

算額（その200）

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

問6 鉤股（直角三角形）に斜線を入れ，区分された領域に大円，小円を入れる。鉤股の和（直角を挟む二辺の長さの和）は 42 寸，弦（斜辺）の長さは 30 寸，小円の径は 8 寸である。大円の径はいかほどか。

鉤股弦それぞれの長さを b, c, a とする。大円，小円の中心座標と半径を (r, y, r)，(x, 4, 4），斜線と斜辺の交点座標を (x2, y2) とおき，以下の連立方程式を解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms a::positive, b::positive, c::positive, x::positive, y::positive, x2::positive, y2::positive, r::positive;

a = 30
eq1 = distance(0, b, c, 0, r, y) - r^2
eq2 = distance(0, 0, x2, y2, r, y) - r^2
eq3 = distance(0, b, c, 0, x, 4) - 4^2
eq4 = distance(0, 0, x2, y2, x, 4) - 4^2
eq5 = b^2 + c^2 - 30^2
eq6 = b + c - 42
eq7 = b/(b - y2) - c/x2;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq8])

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

   -r^2 + (-b*(b*y + c^2 - c*r)/(b^2 + c^2) + y)^2 + (-c*(b^2 - b*y + c*r)/(b^2 + c^2) + r)^2,
   -r^2 + (r - x2*(r*x2 + y*y2)/(x2^2 + y2^2))^2 + (y - y2*(r*x2 + y*y2)/(x2^2 + y2^2))^2,
   (-b*(4*b + c^2 - c*x)/(b^2 + c^2) + 4)^2 + (-c*(b^2 - 4*b + c*x)/(b^2 + c^2) + x)^2 - 16,
   (x - x2*(x*x2 + 4*y2)/(x2^2 + y2^2))^2 + (-y2*(x*x2 + 4*y2)/(x2^2 + y2^2) + 4)^2 - 16,
   b^2 + c^2 - 900,
   b + c - 42,
   b/(b - y2) - c/x2,

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)
   (b, c, x, y, x2, y2, r) = u
   return [
       -r^2 + (-b*(b*y + c^2 - c*r)/(b^2 + c^2) + y)^2 + (-c*(b^2 - b*y + c*r)/(b^2 + c^2) + r)^2,
       -r^2 + (r - x2*(r*x2 + y*y2)/(x2^2 + y2^2))^2 + (y - y2*(r*x2 + y*y2)/(x2^2 + y2^2))^2,
       (-b*(4*b + c^2 - c*x)/(b^2 + c^2) + 4)^2 + (-c*(b^2 - 4*b + c*x)/(b^2 + c^2) + x)^2 - 16,
       (x - x2*(x*x2 + 4*y2)/(x2^2 + y2^2))^2 + (-y2*(x*x2 + 4*y2)/(x2^2 + y2^2) + 4)^2 - 16,
       b^2 + c^2 - 900,
       b + c - 42,
       b/(b - y2) - c/x2, 
   ]
end;

iniv = [big"17.0", 25, 17, 8, 23, 8, 5]
res = nls(H, ini=iniv)
   (BigFloat[17.99999999999999999999999999997392261968940866058339734161588411500118676158968, 24.0000000000000000000000000000260773803105913394166026583841158849915477118184, 11.99999999999999999999996243820109743337714823617693357054397757466207088386991, 9.000000000000000000000022567842421108529450387827899267586371106337607442867201, 11.99999999999999999999999448634115543134366334798061420939477084496247619303071, 9.000000000000000000000005638682909565557685988715612295652539939416822847399611, 4.499999999999999999999988758603463582049975552632475204639736619751622970598919], true)

using Printf
(b, c, x, y, x2, y2, r) = res[1]
@printf("b = %.3f;  c = %.3f;  x = %.3f;  y = %.3f;  x2 = %.3f;  y2 = %.3f;  r = %.3f\n",
       b, c, x, y, x2, y2, r)

   b = 18.000;  c = 24.000;  x = 12.000;  y = 9.000;  x2 = 12.000;  y2 = 9.000;  r = 4.500

大円の半径は 4.5 である。元の単位でいえば，大円の径は 9 寸である。

例によって，算額の図とはかなり異なるものになる。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (b, c, x, y, x2, y2, r) = res[1]
   @printf("b = %.3f;  c = %.3f;  x = %.3f;  y = %.3f;  x2 = %.3f;  y2 = %.3f;  r = %.3f\n",
       b, c, x, y, x2, y2, r)
   plot([0, c, 0, 0], [0, 0, b, 0], color=:black, lw=0.5)
   circle(x, 4, 4)
   circle(r, y, r, :blue)
   segment(0, 0, x2, y2)
   if more == true
       point(c, 0, " c", :green, :left, :bottom)
       point(0, b, " b", :green, :left, :bottom)
       point(x2, y2, " (x2,y2)", :green, :left, :bottom)
       point(x, 4, " 小円:(x,4)", :red, :center)
       point(r, y, " 大円:(r,y)", :blue, :center)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

draw(false)

   b = 18.000;  c = 24.000;  x = 12.000;  y = 9.000;  x2 = 12.000;  y2 = 9.000;  r = 4.500

 

 

算額（その199）

算額（その199）

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

問5 正方形内に4本の斜線を引き，区画された部分に半径の等しい円を5個入れる。正方形の一辺の長さが与えられたとき，円の径を求めよ。

正方形の一辺の長さを x，円の半径を r とする。
⊿ACD は直角三角形（∠ADC = π/3）で，AB = a，AD = x，BC = 2r である。
以下の 2 つの方程式を r, a に対して解く。

using SymPy

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

eq1 = a^2 + (a + 2r)^2 - x^2
eq2 = 2r*(x + 2a + 2r) + 4r^2 - x^2
res = solve([eq1, eq2], (r, a))

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

println("r = ", res[1][1])
println("a = ", res[1][2])

   r = x*(-1/4 + sqrt(3)/4)
   a = x/2

円の直径は「(√3 - 1)/2 × 正方形の一辺の長さ」である。

draw(10, true)

   AD = x = 10.00000;  BC = 2r = 3.66025;  AB = CD = a = 5.00000

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(x, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = x*(-1/4 + sqrt(3)/4)
   a = x/2
   @printf("AD = x = %.5f;  BC = 2r = %.5f;  AB = CD = a = %.5f\n", x, 2r, a)
   plot([0, x, x, 0, 0], [0, 0, x, x, 0])
   circle(x/2, x/2, r)
   cosine = (a + 2r)/x
   sine   = a/x
   segment(0, 0, x - a*sine, a*cosine, :green)
   segment(x, 0, x - a*cosine, x - a*sine, :green)
   segment(x, x, a*sine, x - a*cosine, :green)
   segment(0, x, a*cosine, a*sine, :green)
   circle(x/2 + 2r*cosine, x/2 + 2r*sine, r)
   circle(x/2 - 2r*sine, x/2 + 2r*cosine, r)
   circle(x/2 - 2r*cosine, x/2 - 2r*sine, r)
   circle(x/2 + 2r*sine, x/2 - 2r*cosine, r)
   if more == true
       point(x, x, "A ", :blue, :right, :bottom)
       point(x - a*cosine, x - a*sine, " B")
       point(a*sine, x - a*cosine, " C")
       point(0, x, " D", :blue, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その198）

算額（その198）

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

問4 全円の中に大円，中円，小円が入っている。大円，中円，小円の径が与えられたとき，全円の径を求めよ。

図のように，全円，大円，中円，小円の半径を R, r1, r2, r3 とし，小円の中心座標を (x3, y3) とする。連立方程式を立て，R, x3, y3 について解く。

using SymPy

@syms r1::positive, r2::positive, r3::positive, R::positive,
     x3::positive, y3::negative;

eq1 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq2 = x3^2 + y3^2 - (R - r3)^2
eq3 = x3^2 + (R - 2r1 - r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (R, x3, y3));

println("R = ", res[1][1])
println("x3 = ", res[1][2])
println("y3 = ", res[1][3])

   R = r1^2*(r1 + r2 + r3)/(r1^2 + r1*r2 - r2*r3)
   x3 = 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(r1 + r2)
   y3 = -(r1^4 + 2*r1^3*r2 + r1^2*r2^2 - 3*r1^2*r2*r3 - 3*r1*r2^2*r3 - r1*r2*r3^2 + r2^2*r3^2)/((r1 + r2)*(r1^2 + r1*r2 - r2*r3))

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(r1, r2, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = r1^2*(r1 + r2 + r3)/(r1^2 + r1*r2 - r2*r3)
   x3 = 2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(r1 + r2)
   y3 = -(r1^4 + 2*r1^3*r2 + r1^2*r2^2 - 3*r1^2*r2*r3 - 3*r1*r2^2*r3 - r1*r2*r3^2 + r2^2*r3^2)/((r1 + r2)*(r1^2 + r1*r2 - r2*r3))
   @printf("R = %.5f,  x3 = %.5f,  y3 = %.5f\n", R, x3, y3)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   circle(0, R - 2r1 - r2, r2, :orange)
   if more == true
       point(0, R - r1, " R-r1", :blue)
       point(0, R - 2r1 - r2, " R-2r1-r2", :orange)
       point(x3, y3, "(x3,y3)", :green)
       point(0.4r1, R - 0.6r1, "大円", :blue, mark=false)
       point(0.4r2, R - 2r1 - 0.6r2, "中円", :orange, mark=false)
       point(x3 + 0.4r3, y3 + 0.4r3, "小円", :green, mark=false)
       point(0.75R, 0.75R, "全円", :red, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

draw(6, 2, 4, false)

   R = 10.80000,  x3 = 6.00000,  y3 = -3.20000

算額（その197）

算額（その197）

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

問2 大円の中に 3 個の小円が互いに外接して入っている。小円の径が 5 寸のとき，大円の径はいかほどか。

大円半径を r とおき方程式を解く。

using SymPy

@syms r::positive;

eq1 = (√3(r - 5)/2)^2 + (3(r - 5)/2)^2 - (2*5)^2 
res = solve(eq1)[1]
res.evalf() |> println

   10.7735026918963

小円の径が 5 寸のとき，大円の径は 10寸7分7厘3毛5糸である。
算額では 17.77352 寸としているが，末尾桁は 0 または 17.773502 寸の写し間違いであろう。17.773502 寸としても，正しい末尾は 3 である。

なお，eq1 を簡約化すると 3r^2 - 30r - 25 = 0 である。
術では「4 を 3 で割り，平方根を求めて 1 を加え，小円の径を掛ける」とあるが，これは二次方程式の正の解を計算することを表している。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:blue, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = res[1]
   println("r = $r")
   plot()
   circle(0, 0, r)
   circle(0, r - 5, 5, :blue)
   circle(√3(r - 5)/2, (5 - r)/2, 5, :blue)
   circle(-√3(r - 5)/2, (5 - r)/2, 5, :blue)
   if more == true
       point(0, r - 5, " r-5")
       point(0, r, " r", :red)
       point(√3(r - 5)/2, (5 - r)/2, "√3(r-5)/2,(5-r)/2", :blue, :center, :top)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その196）

算額（その196）

中村信弥「改訂増補 長野県の算額」（91,165）
http://www.wasan.jp/zoho/zoho.html
長野県下高井郡木島平村穂高 一川谷大元神社 文化8年(1811)
長野県飯山市下木島 鳥出神社 天保14年(1843)

問1 2つの大円が交わっている部分に中円，それ以外の部分に小円が入っている。大円，中円の径が与えられたとき，小円の径を求めよ。

大円，中円，小円の半径を r1, r2, r3 とおく。
小円の中心の y 座標を y3 とおく。
以下の連立方程式を r3, y3 について解く。

using SymPy

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

eq1 = r3^2 + (r1 - r2 + y3)^2 - (r1 + r3)^2 |> expand
eq2 = r3^2 + (r2 - r1 + y3)^2 - (r1 - r3)^2 |> expand
solve([eq1, eq2], (r3, y3))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(r2)*(r1 - r2)*sqrt(2*r1 - r2)/r1, sqrt(r2)*sqrt(2*r1 - r2))

式を変形して，小円の径は「（1 - 中円の径/大円の径）× 平方根(2×大円の径×中円の径 - 中円の径の二乗）」である。当然であるが，文中の「径」は「直径」と読み替えても「半径」と読み替えてもよい。

中円の径が大円の径に比して小さいと，「小円が中円より大きい」ということになる（その限界値を求めるのも一興である）。

r1 = 10;  r2 = 4.56310987307924;  r3 = 4.5631098730792345;  y3 = 8.392867552141613

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, y3) = (sqrt(r2)*(r1 - r2)*sqrt(2*r1 - r2)/r1, sqrt(r2)*sqrt(2*r1 - r2))
   println("r1 = $r1;  r2 = $r2;  r3 = $(r3);  y3 = $y3")
   plot()
   circle(0, 0, r2)
   circle(0, r1 - r2, r1, :green)
   circle(0, r2 - r1, r1, :green)
   circle(r3, y3, r3, :blue)
   circle(-r3, y3, r3, :blue)
   circle(r3, -y3, r3, :blue)
   circle(-r3, -y3, r3, :blue)
   if more == true
       point(0, r1 - r2, "r1-r2")
       point(0, r2 - r1, "r2-r1")
       point(r3, y3, "(r3,y3)", :blue)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その195）

算額（その195）

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

長方形内に2本の斜線を引き，4区分されたそれぞれの領域に大円，甲円，乙円，丙円を入れる。甲円，丙円の径がそれぞれ 126寸，66寸のとき，乙円の径はいかほどか。

長方形の長辺，短辺の長さをそれぞれ x, y とおく。その他，
大円の中心座標，半径を (x - r0, r0, r0)
甲円の中心座標，半径を (r1, r1, r1) ただし r1 = 63
乙円の中心座標，半径を (x2, y - r2, r2)
丙円の中心座標，半径を (x3, r3, r3) ただし r3 = 33
として，各円の中心から2本の斜線までの距離とそれぞれの円の半径の関係についての連立方程式を nlsolve() で解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms a::positive, b::positive, c::positive, x::positive, y::positive,
     x2::positive, r0::positive, r2::positive, x3::positive;

r1 = 126 // 2
r3 = 66 // 2
r0 = y / 2
eq1 = distance(a, 0, c, y, x- r0, r0) - r0^2
eq2 = distance(0, y, b, 0, x- r0, r0) - r0^2
eq3 = distance(a, 0, c, y, x2, y - r2) - r2^2
eq4 = distance(0, y, b, 0, x2, y - r2) - r2^2
eq5 = distance(a, 0, c, y, r1, r1) - r1^2
eq6 = distance(0, y, b, 0, r1, r1) - r1^2
eq7 = distance(a, 0, c, y, x3, r3) - r3^2
eq8 = distance(0, y, b, 0, x3, r3) - r3^2;

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

   -y^2/4 + (y/2 - y*(a^2 - 2*a*c + c^2 + y^2 + (a - c)*(a + c - 2*x + y))/(2*(a^2 - 2*a*c + c^2 + y^2)))^2 + (x - y/2 - (y^2*(a + c - 2*x + y)/2 + (2*x - y)*(a^2 - 2*a*c + c^2 + y^2)/2)/(a^2 - 2*a*c + c^2 + y^2))^2,
   -y^2/4 + (y/2 - y*(2*b^2 - 2*b*x + b*y + y^2)/(2*(b^2 + y^2)))^2 + (-b*(2*b*x - b*y + y^2)/(2*(b^2 + y^2)) + x - y/2)^2,
   -r2^2 + (x2 - (x2*(a^2 - 2*a*c + c^2 + y^2) + y*(a*r2 - c*r2 + c*y - x2*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-r2 + y - y*(a^2 - a*c - a*x2 + c*x2 - r2*y + y^2)/(a^2 - 2*a*c + c^2 + y^2))^2,
   -r2^2 + (-b*(b*x2 + r2*y)/(b^2 + y^2) + x2)^2 + (-r2 + y - y*(b^2 - b*x2 - r2*y + y^2)/(b^2 + y^2))^2,
   (63 - (63*a^2 - 126*a*c + a*y^2 - 63*a*y + 63*c^2 + 63*c*y)/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - 63*a + 63*c + 63*y)/(a^2 - 2*a*c + c^2 + y^2) + 63)^2 - 3969,
   (-b*(63*b + y^2 - 63*y)/(b^2 + y^2) + 63)^2 + (-y*(b^2 - 63*b + 63*y)/(b^2 + y^2) + 63)^2 - 3969,
   (x3 - (x3*(a^2 - 2*a*c + c^2 + y^2) + y*(a*y - 33*a + 33*c - x3*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - a*x3 + c*x3 + 33*y)/(a^2 - 2*a*c + c^2 + y^2) + 33)^2 - 1089,
   (-b*(b*x3 + y^2 - 33*y)/(b^2 + y^2) + x3)^2 + (-y*(b^2 - b*x3 + 33*y)/(b^2 + y^2) + 33)^2 - 1089,

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)
   (a, b, c, x, y, x2, r2, x3) = u
   return [
       -y^2/4 + (y/2 - y*(2*a^2 - 2*a*c - 2*a*x + a*y + 2*c*x - c*y + y^2)/(2*(a^2 - 2*a*c + c^2 + y^2)))^2 + (x - y/2 - (a^2*x - a^2*y/2 - 2*a*c*x + a*c*y + a*y^2/2 + c^2*x - c^2*y/2 + c*y^2/2)/(a^2 - 2*a*c + c^2 + y^2))^2,
       -y^2/4 + (y/2 - y*(2*b^2 - 2*b*x + b*y + y^2)/(2*(b^2 + y^2)))^2 + (-b*(2*b*x - b*y + y^2)/(2*(b^2 + y^2)) + x - y/2)^2,
       -r2^2 + (x2 - (x2*(a^2 - 2*a*c + c^2 + y^2) + y*(a*r2 - c*r2 + c*y - x2*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-r2 + y - y*(a^2 - a*c - a*x2 + c*x2 - r2*y + y^2)/(a^2 - 2*a*c + c^2 + y^2))^2,
       -r2^2 + (-b*(b*x2 + r2*y)/(b^2 + y^2) + x2)^2 + (-r2 + y - y*(b^2 - b*x2 - r2*y + y^2)/(b^2 + y^2))^2,
       (63 - (63*a^2 - 126*a*c + a*y^2 - 63*a*y + 63*c^2 + 63*c*y)/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - 63*a + 63*c + 63*y)/(a^2 - 2*a*c + c^2 + y^2) + 63)^2 - 3969,
       (-b*(63*b + y^2 - 63*y)/(b^2 + y^2) + 63)^2 + (-y*(b^2 - 63*b + 63*y)/(b^2 + y^2) + 63)^2 - 3969,
       (x3 - (x3*(a^2 - 2*a*c + c^2 + y^2) + y*(a*y - 33*a + 33*c - x3*y))/(a^2 - 2*a*c + c^2 + y^2))^2 + (-y*(a^2 - a*c - a*x3 + c*x3 + 33*y)/(a^2 - 2*a*c + c^2 + y^2) + 33)^2 - 1089,
       (-b*(b*x3 + y^2 - 33*y)/(b^2 + y^2) + x3)^2 + (-y*(b^2 - b*x3 + 33*y)/(b^2 + y^2) + 33)^2 - 1089,     
   ]
end;

iniv = [big"85.0", 315, 310, 420, 170, 175, 45, 185]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[84.0000000000000000000001735702499203810234362830123112659631373822954941885675, 314.9999999999999999999947583595807515206167302592172476055552893595369151008652, 307.9999999999999999999947727498833028947963179438839311145144054791595619767401, 419.9999999999999999999951927515441018951317502337587509125424356893062721470855, 168.0000000000000000000006020767375679450188933432128189504750790778747497805281, 175.999999999999999999997350763673708796543835009537017244584622091582997453893, 44.00000000000000000000017089021421936614679777263092841262631968462778651775099, 182.9999999999999999999975156955872708477881618570615211838481194365997674781264], true)

   a = 84.00000, b = 315.00000, c = 308.00000
   x = 420.00000, y = 168.00000, x2 = 176.00000, r2 = 44.00000, x3 = 183.00000
   乙円の直径 = 88.00000

乙円の半径は 44 なので，直径は 88 寸である。
図に描いてみると，算額の図は正確なものではない（与えられた条件を満たすものではない）ことがわかる。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c, x, y, x2, r2, x3) = res[1]
   r0 = y/2
   @printf("a = %.5f, b = %.5f, c = %.5f\nx = %.5f, y = %.5f, x2 = %.5f, r2 = %.5f, x3 = %.5f\n",
       a, b, c, x, y, x2, r2, x3)
   @printf("乙円の直径 = %.5f\n", 2r2)
   plot([0, x, x, 0, 0], [0, 0, y, y, 0], color=:black, lw=0.5)
   circle(x - r0, r0, r0)
   circle(r1, r1, r1, :blue)
   circle(x2, y - r2, r2, :green)
   circle(x3, r3, r3, :orange)
   segment(a, 0, c, y)
   segment(0, y, b, 0)
   if more == true
       point(x - r0, r0, "大:(x-r0,r0,r0)", :red, :center)
       point(r1, r1, "甲:(r1,r1,r1)", :blue, :center)
       point(x2, y - r2, "乙:(x2,y-r2,r2)", :green, :center)
       point(x3, r3, "丙:(x3,r3,r3)", :orange, :center)
       point(0, y, " y", :black, :left, :bottom)
       point(x, 0, " x", :black, :left, :bottom)
       point(a, 0, "a ", :black, :right, :bottom)
       point(b, 0, " b", :black, :left, :bottom)
       point(c, y, " c", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

算額（その194）

算額（その194）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（89）
長野県中野市田上 田上観音堂 文化6年(1809)

問3. 大円の中に，甲，乙，丙，丁，戊，己の 5 個の円が外円に内接すると同時に隣同士の円は外接している。さらに，子，丑，寅，卯，辰の 5 個の円が甲円と乙，丙，丁，戊，己円と外接している。大円の径が 76 寸のとき，辰円の径はいかほどか。

大円の中心座標，半径を (0, R0, R0) ただし R0 = 2
甲円の中心座標，半径を (0, R1, R1) および (0, 3R1, R1) ただし R1 = 1
乙円の中心座標，半径を (X2, Y2, R2)
丙円の中心座標，半径を (X3, Y3, R3)
丁円の中心座標，半径を (X4, Y4, R4)
戊円の中心座標，半径を (X5, Y5, R5)
己円の中心座標，半径を (X6, Y6, R6)
子円の中心座標，半径を (x2, y2, r2)
丑円の中心座標，半径を (x3, y3, r3)
寅円の中心座標，半径を (x4, y4, r4)
卯円の中心座標，半径を (x5, y5, r5)
辰円の中心座標，半径を (x6, y6, r6)
とおき方程式を解く。なお，解として得られるのは大円の半径を 2 としたときの各円の半径なので，元の大円の大きさ，単位（寸）での直径はそれぞれを 38 倍する。大円の半径は 2 の倍数のほうが良いが，大きさは任意に指定できない。SymPy のせいであろうが，16 以上にすると平方根の計算においてエラーになる。

using SymPy

@syms R0::positive, R1::positive,
     X2::positive, Y2::positive, R2::positive,
     x2::positive, y2::positive, r2::positive;

R0 = 2  # 大円の半径
R1 = R0 // 2  # 甲円の半径
X2 = R0 - R2
Y2 = 2R1
eq1 = X2^2 + R1^2 - (R1 + R2)^2 
R2 = solve(eq1, R2)[1]  # 乙円の半径

$\frac{2}{3}$

R0 - R0//3, Y2, R2  # 乙円の中心座標と半径

   (4//3, 2//1, 2/3)

eq2 = x2^2 + (y2 - R1)^2 - (R1 + r2)^2
eq3 = x2^2 + (3R1 - y2)^2 - (R1 + r2)^2
eq4 = (R0 - R2 - x2)^2 + (y2 - 2R1)^2 - (r2 + R2)^2;
(r2, x2, y2) = solve([eq2, eq3, eq4], (r2, x2, y2))[1]  # 子円の中心座標と半径

   (2/15, 8/15, 2)

乙円から始めて，丙円以降の中心座標と半径を求める漸化式を求める。

@syms Xn, Yn, Rn, Xnp1, Ynp1, Rnp1,
     xn, yn, rn, xnp1, ynp1, rnp1;
eq5 = (Xn - Xnp1)^2 + (Yn - Ynp1)^2 - (Rn + Rnp1)^2 |> expand
eq6 = Xnp1^2 + (R1 - Ynp1)^2 - (R1 + Rnp1)^2 |> expand
eq7 = Xnp1^2 + (2R1 - Ynp1)^2 - (R0 - Rnp1)^2 |> expand;
res = solve([eq5, eq6, eq7], (Xnp1, Ynp1, Rnp1));

println(res[2])

   出力省略

2番めの組が適切解である。これに基づいて，現在の円の中心座標と半径を与えて，次の円の中心座標と半径を求める関数を定義する。

NEXTCIRCLE(Xn, Yn, Rn) = return ((-4*Rn^2*Xn - 4*Rn*Xn + sqrt(2)*Rn*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2) + 4*Xn^3 + 4*Xn*Yn^2 - 12*Xn*Yn + 16*Xn + 3*sqrt(2)*Yn*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2) - 4*sqrt(2)*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2))/(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16), 3*(-Rn^3 - 3*Rn^2*Yn + 4*Rn^2 + Rn*Xn^2 + Rn*Yn^2 + 3*Xn^2*Yn + 4*Xn^2 - 2*sqrt(2)*Xn*sqrt(-Rn^4 - 2*Rn^3 + 2*Rn^2*Xn^2 + 2*Rn^2*Yn^2 - 6*Rn^2*Yn + 8*Rn^2 + 2*Rn*Xn^2 + 2*Rn*Yn^2 - Xn^4 - 2*Xn^2*Yn^2 + 6*Xn^2*Yn - Yn^4 + 6*Yn^3 - 8*Yn^2) + 3*Yn^3 - 4*Yn^2)/(2*(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16)), -sqrt(2)*Xn*sqrt(-(Rn^2 - 2*Rn - Xn^2 - Yn^2 + 2*Yn)*(Rn^2 + 4*Rn - Xn^2 - Yn^2 + 4*Yn))/(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16) - (Rn^3 + 3*Rn^2*Yn - 4*Rn^2 - Rn*Xn^2 - Rn*Yn^2 - 3*Xn^2*Yn - 4*Xn^2 - 3*Yn^3 + 4*Yn^2)/(2*(Rn^2 + 6*Rn*Yn - 8*Rn + 8*Xn^2 + 9*Yn^2 - 24*Yn + 16)));

下側の甲円と乙円，丙円から始めて，丑円以降の中心座標と半径を求める漸化式を求める。

eq8 = xnp1^2 + (R1 -ynp1)^2 - (R1 + rnp1)^2 |> expand
eq9 = (Xn - xnp1)^2 + (Yn - ynp1)^2 - (Rn + rnp1)^2 |> expand
eq10 = (Xnp1 - xnp1)^2 + (ynp1 - Ynp1)^2 - (Rnp1 + rnp1)^2 |> expand
res2 = solve([eq8, eq9, eq10], (xnp1, ynp1, rnp1));

println(res2[1])

    出力省略

1番めの組が適切解である。これに基づいて，丑円以降の円の中心座標と半径を求める関数を定義する。

function nextcircle(Xn, Yn, Rn, Xnp1, Ynp1, Rnp1)
   term = sqrt(-Rn^4*Rnp1^2 + 2*Rn^4*Rnp1 + Rn^4*Xnp1^2 + Rn^4*Ynp1^2 - 2*Rn^4*Ynp1 + 2*Rn^3*Rnp1^3 - 2*Rn^3*Rnp1^2 - 2*Rn^3*Rnp1*Xnp1^2 - 2*Rn^3*Rnp1*Ynp1^2 + 4*Rn^3*Rnp1*Ynp1 - 4*Rn^3*Rnp1 - 2*Rn^3*Xnp1^2 - 2*Rn^3*Ynp1^2 + 4*Rn^3*Ynp1 - Rn^2*Rnp1^4 - 2*Rn^2*Rnp1^3 + 2*Rn^2*Rnp1^2*Xn^2 - 2*Rn^2*Rnp1^2*Xn*Xnp1 + 2*Rn^2*Rnp1^2*Xnp1^2 + 2*Rn^2*Rnp1^2*Yn^2 - 2*Rn^2*Rnp1^2*Yn*Ynp1 - 2*Rn^2*Rnp1^2*Yn + 2*Rn^2*Rnp1^2*Ynp1^2 - 2*Rn^2*Rnp1^2*Ynp1 + 8*Rn^2*Rnp1^2 - 4*Rn^2*Rnp1*Xn^2 + 4*Rn^2*Rnp1*Xn*Xnp1 + 2*Rn^2*Rnp1*Xnp1^2 - 4*Rn^2*Rnp1*Yn^2 + 4*Rn^2*Rnp1*Yn*Ynp1 + 4*Rn^2*Rnp1*Yn + 2*Rn^2*Rnp1*Ynp1^2 - 8*Rn^2*Rnp1*Ynp1 - 2*Rn^2*Xn^2*Xnp1^2 - 2*Rn^2*Xn^2*Ynp1^2 + 4*Rn^2*Xn^2*Ynp1 + 2*Rn^2*Xn*Xnp1^3 + 2*Rn^2*Xn*Xnp1*Ynp1^2 - 4*Rn^2*Xn*Xnp1*Ynp1 - Rn^2*Xnp1^4 - 2*Rn^2*Xnp1^2*Yn^2 + 2*Rn^2*Xnp1^2*Yn*Ynp1 + 2*Rn^2*Xnp1^2*Yn - 2*Rn^2*Xnp1^2*Ynp1^2 + 2*Rn^2*Xnp1^2*Ynp1 - 2*Rn^2*Yn^2*Ynp1^2 + 4*Rn^2*Yn^2*Ynp1 + 2*Rn^2*Yn*Ynp1^3 - 2*Rn^2*Yn*Ynp1^2 - 4*Rn^2*Yn*Ynp1 - Rn^2*Ynp1^4 + 2*Rn^2*Ynp1^3 + 2*Rn*Rnp1^4 - 2*Rn*Rnp1^3*Xn^2 - 2*Rn*Rnp1^3*Yn^2 + 4*Rn*Rnp1^3*Yn - 4*Rn*Rnp1^3 + 2*Rn*Rnp1^2*Xn^2 + 4*Rn*Rnp1^2*Xn*Xnp1 - 4*Rn*Rnp1^2*Xnp1^2 + 2*Rn*Rnp1^2*Yn^2 + 4*Rn*Rnp1^2*Yn*Ynp1 - 8*Rn*Rnp1^2*Yn - 4*Rn*Rnp1^2*Ynp1^2 + 4*Rn*Rnp1^2*Ynp1 + 2*Rn*Rnp1*Xn^2*Xnp1^2 + 2*Rn*Rnp1*Xn^2*Ynp1^2 - 4*Rn*Rnp1*Xn^2*Ynp1 + 4*Rn*Rnp1*Xn^2 - 8*Rn*Rnp1*Xn*Xnp1 + 2*Rn*Rnp1*Xnp1^2*Yn^2 - 4*Rn*Rnp1*Xnp1^2*Yn + 4*Rn*Rnp1*Xnp1^2 + 2*Rn*Rnp1*Yn^2*Ynp1^2 - 4*Rn*Rnp1*Yn^2*Ynp1 + 4*Rn*Rnp1*Yn^2 - 4*Rn*Rnp1*Yn*Ynp1^2 + 4*Rn*Rnp1*Ynp1^2 + 2*Rn*Xn^2*Xnp1^2 + 2*Rn*Xn^2*Ynp1^2 - 4*Rn*Xn^2*Ynp1 - 4*Rn*Xn*Xnp1^3 - 4*Rn*Xn*Xnp1*Ynp1^2 + 8*Rn*Xn*Xnp1*Ynp1 + 2*Rn*Xnp1^4 + 2*Rn*Xnp1^2*Yn^2 - 4*Rn*Xnp1^2*Yn*Ynp1 + 4*Rn*Xnp1^2*Ynp1^2 - 4*Rn*Xnp1^2*Ynp1 + 2*Rn*Yn^2*Ynp1^2 - 4*Rn*Yn^2*Ynp1 - 4*Rn*Yn*Ynp1^3 + 8*Rn*Yn*Ynp1^2 + 2*Rn*Ynp1^4 - 4*Rn*Ynp1^3 + Rnp1^4*Xn^2 + Rnp1^4*Yn^2 - 2*Rnp1^4*Yn - 2*Rnp1^3*Xn^2 - 2*Rnp1^3*Yn^2 + 4*Rnp1^3*Yn - Rnp1^2*Xn^4 + 2*Rnp1^2*Xn^3*Xnp1 - 2*Rnp1^2*Xn^2*Xnp1^2 - 2*Rnp1^2*Xn^2*Yn^2 + 2*Rnp1^2*Xn^2*Yn*Ynp1 + 2*Rnp1^2*Xn^2*Yn - 2*Rnp1^2*Xn^2*Ynp1^2 + 2*Rnp1^2*Xn^2*Ynp1 + 2*Rnp1^2*Xn*Xnp1*Yn^2 - 4*Rnp1^2*Xn*Xnp1*Yn - 2*Rnp1^2*Xnp1^2*Yn^2 + 4*Rnp1^2*Xnp1^2*Yn - Rnp1^2*Yn^4 + 2*Rnp1^2*Yn^3*Ynp1 + 2*Rnp1^2*Yn^3 - 2*Rnp1^2*Yn^2*Ynp1^2 - 2*Rnp1^2*Yn^2*Ynp1 + 4*Rnp1^2*Yn*Ynp1^2 - 4*Rnp1^2*Yn*Ynp1 + 2*Rnp1*Xn^4 - 4*Rnp1*Xn^3*Xnp1 + 2*Rnp1*Xn^2*Xnp1^2 + 4*Rnp1*Xn^2*Yn^2 - 4*Rnp1*Xn^2*Yn*Ynp1 - 4*Rnp1*Xn^2*Yn + 2*Rnp1*Xn^2*Ynp1^2 - 4*Rnp1*Xn*Xnp1*Yn^2 + 8*Rnp1*Xn*Xnp1*Yn + 2*Rnp1*Xnp1^2*Yn^2 - 4*Rnp1*Xnp1^2*Yn + 2*Rnp1*Yn^4 - 4*Rnp1*Yn^3*Ynp1 - 4*Rnp1*Yn^3 + 2*Rnp1*Yn^2*Ynp1^2 + 8*Rnp1*Yn^2*Ynp1 - 4*Rnp1*Yn*Ynp1^2 + Xn^4*Xnp1^2 + Xn^4*Ynp1^2 - 2*Xn^4*Ynp1 - 2*Xn^3*Xnp1^3 - 2*Xn^3*Xnp1*Ynp1^2 + 4*Xn^3*Xnp1*Ynp1 + Xn^2*Xnp1^4 + 2*Xn^2*Xnp1^2*Yn^2 - 2*Xn^2*Xnp1^2*Yn*Ynp1 - 2*Xn^2*Xnp1^2*Yn + 2*Xn^2*Xnp1^2*Ynp1^2 - 2*Xn^2*Xnp1^2*Ynp1 + 2*Xn^2*Yn^2*Ynp1^2 - 4*Xn^2*Yn^2*Ynp1 - 2*Xn^2*Yn*Ynp1^3 + 2*Xn^2*Yn*Ynp1^2 + 4*Xn^2*Yn*Ynp1 + Xn^2*Ynp1^4 - 2*Xn^2*Ynp1^3 - 2*Xn*Xnp1^3*Yn^2 + 4*Xn*Xnp1^3*Yn - 2*Xn*Xnp1*Yn^2*Ynp1^2 + 4*Xn*Xnp1*Yn^2*Ynp1 + 4*Xn*Xnp1*Yn*Ynp1^2 - 8*Xn*Xnp1*Yn*Ynp1 + Xnp1^4*Yn^2 - 2*Xnp1^4*Yn + Xnp1^2*Yn^4 - 2*Xnp1^2*Yn^3*Ynp1 - 2*Xnp1^2*Yn^3 + 2*Xnp1^2*Yn^2*Ynp1^2 + 2*Xnp1^2*Yn^2*Ynp1 - 4*Xnp1^2*Yn*Ynp1^2 + 4*Xnp1^2*Yn*Ynp1 + Yn^4*Ynp1^2 - 2*Yn^4*Ynp1 - 2*Yn^3*Ynp1^3 + 2*Yn^3*Ynp1^2 + 4*Yn^3*Ynp1 + Yn^2*Ynp1^4 + 2*Yn^2*Ynp1^3 - 8*Yn^2*Ynp1^2 - 2*Yn*Ynp1^4 + 4*Yn*Ynp1^3)
   return ((Rn^3*Rnp1*Xnp1 - Rn^3*Xnp1 - Rn^2*Rnp1^2*Xn - Rn^2*Rnp1^2*Xnp1 + 2*Rn^2*Rnp1*Xn - Rn^2*Rnp1*Xnp1 + Rn^2*Xn*Ynp1^2 - 2*Rn^2*Xn*Ynp1 + Rn^2*Xnp1^3 - Rn^2*Xnp1*Yn*Ynp1 + Rn^2*Xnp1*Yn + Rn^2*Xnp1*Ynp1^2 - Rn^2*Xnp1*Ynp1 + 2*Rn^2*Xnp1 + Rn*Rnp1^3*Xn - Rn*Rnp1^2*Xn + 2*Rn*Rnp1^2*Xnp1 - Rn*Rnp1*Xn^2*Xnp1 - Rn*Rnp1*Xn*Xnp1^2 - Rn*Rnp1*Xn*Ynp1^2 + 2*Rn*Rnp1*Xn*Ynp1 - 2*Rn*Rnp1*Xn - Rn*Rnp1*Xnp1*Yn^2 + 2*Rn*Rnp1*Xnp1*Yn - 2*Rn*Rnp1*Xnp1 + Rn*Xn^2*Xnp1 + Rn*Xn*Xnp1^2 - Rn*Xn*Ynp1^2 + 2*Rn*Xn*Ynp1 - 2*Rn*Xnp1^3 + Rn*Xnp1*Yn^2 + 2*Rn*Xnp1*Yn*Ynp1 - 4*Rn*Xnp1*Yn - 2*Rn*Xnp1*Ynp1^2 + 2*Rn*Xnp1*Ynp1 - Rn*Ynp1*term + Rn*term - Rnp1^3*Xn + Rnp1^2*Xn^3 + Rnp1^2*Xn*Yn^2 - Rnp1^2*Xn*Yn*Ynp1 - Rnp1^2*Xn*Yn + Rnp1^2*Xn*Ynp1 + 2*Rnp1^2*Xn + Rnp1^2*Xnp1*Yn^2 - 2*Rnp1^2*Xnp1*Yn - 2*Rnp1*Xn^3 + Rnp1*Xn^2*Xnp1 + Rnp1*Xn*Xnp1^2 - 2*Rnp1*Xn*Yn^2 + 2*Rnp1*Xn*Yn*Ynp1 + 2*Rnp1*Xn*Yn + Rnp1*Xn*Ynp1^2 - 4*Rnp1*Xn*Ynp1 - Rnp1*Xnp1*Yn^2 + 2*Rnp1*Xnp1*Yn + Rnp1*Yn*term - Rnp1*term - Xn^3*Ynp1^2 + 2*Xn^3*Ynp1 + Xn^2*Xnp1*Yn*Ynp1 - Xn^2*Xnp1*Yn - Xn^2*Xnp1*Ynp1 + Xn*Xnp1^2*Yn*Ynp1 - Xn*Xnp1^2*Yn - Xn*Xnp1^2*Ynp1 - Xn*Yn^2*Ynp1^2 + 2*Xn*Yn^2*Ynp1 + Xn*Yn*Ynp1^3 - Xn*Yn*Ynp1^2 - 2*Xn*Yn*Ynp1 - Xn*Ynp1^3 + 2*Xn*Ynp1^2 - Xnp1^3*Yn^2 + 2*Xnp1^3*Yn + Xnp1*Yn^3*Ynp1 - Xnp1*Yn^3 - Xnp1*Yn^2*Ynp1^2 - Xnp1*Yn^2*Ynp1 + 2*Xnp1*Yn^2 + 2*Xnp1*Yn*Ynp1^2 - 2*Xnp1*Yn*Ynp1 - Yn*term + Ynp1*term)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)), (Rn^3*Rnp1*Ynp1 - Rn^3*Rnp1 - Rn^3*Ynp1 + Rn^3 - Rn^2*Rnp1^2*Yn - Rn^2*Rnp1^2*Ynp1 + 2*Rn^2*Rnp1^2 + 2*Rn^2*Rnp1*Yn - Rn^2*Rnp1*Ynp1 - Rn^2*Rnp1 - Rn^2*Xn*Xnp1*Ynp1 + Rn^2*Xn*Xnp1 + Rn^2*Xnp1^2*Yn + Rn^2*Xnp1^2*Ynp1 - Rn^2*Yn + Rn^2*Ynp1^3 - Rn^2*Ynp1^2 + Rn^2*Ynp1 + Rn*Rnp1^3*Yn - Rn*Rnp1^3 - Rn*Rnp1^2*Yn + 2*Rn*Rnp1^2*Ynp1 - Rn*Rnp1^2 - Rn*Rnp1*Xn^2*Ynp1 + Rn*Rnp1*Xn^2 - 4*Rn*Rnp1*Xn*Xnp1 - Rn*Rnp1*Xnp1^2*Yn + Rn*Rnp1*Xnp1^2 - Rn*Rnp1*Yn^2*Ynp1 + Rn*Rnp1*Yn^2 - Rn*Rnp1*Yn*Ynp1^2 + Rn*Rnp1*Ynp1^2 + Rn*Xn^2*Ynp1 - Rn*Xn^2 + 2*Rn*Xn*Xnp1*Ynp1 + 2*Rn*Xn*Xnp1 - Rn*Xnp1^2*Yn - 2*Rn*Xnp1^2*Ynp1 - Rn*Xnp1^2 + Rn*Xnp1*term + Rn*Yn^2*Ynp1 - Rn*Yn^2 + Rn*Yn*Ynp1^2 - 2*Rn*Ynp1^3 + Rn*Ynp1^2 - Rnp1^3*Yn + Rnp1^3 + Rnp1^2*Xn^2*Yn + Rnp1^2*Xn^2*Ynp1 - Rnp1^2*Xn*Xnp1*Yn + Rnp1^2*Xn*Xnp1 + Rnp1^2*Yn^3 - Rnp1^2*Yn^2 + Rnp1^2*Yn - Rnp1^2*Ynp1 - 2*Rnp1*Xn^2*Yn - Rnp1*Xn^2*Ynp1 - Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1*Yn + 2*Rnp1*Xn*Xnp1 - Rnp1*Xn*term + Rnp1*Xnp1^2*Yn - Rnp1*Xnp1^2 - 2*Rnp1*Yn^3 + Rnp1*Yn^2*Ynp1 + Rnp1*Yn^2 + Rnp1*Yn*Ynp1^2 - Rnp1*Ynp1^2 + Xn^3*Xnp1*Ynp1 - Xn^3*Xnp1 - Xn^2*Xnp1^2*Yn - Xn^2*Xnp1^2*Ynp1 + 2*Xn^2*Xnp1^2 + Xn^2*Yn - Xn^2*Ynp1^3 + Xn^2*Ynp1^2 + Xn^2*Ynp1 + Xn*Xnp1^3*Yn - Xn*Xnp1^3 + Xn*Xnp1*Yn^2*Ynp1 - Xn*Xnp1*Yn^2 + Xn*Xnp1*Yn*Ynp1^2 - 2*Xn*Xnp1*Yn - Xn*Xnp1*Ynp1^2 - 2*Xn*Xnp1*Ynp1 + Xn*term - Xnp1^2*Yn^3 + Xnp1^2*Yn^2 + Xnp1^2*Yn + Xnp1^2*Ynp1 - Xnp1*term + Yn^3 - Yn^2*Ynp1 - Yn*Ynp1^2 + Ynp1^3)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)), sqrt(-(Rn^2 - 2*Rn - Xn^2 - Yn^2 + 2*Yn)*(Rnp1^2 - 2*Rnp1 - Xnp1^2 - Ynp1^2 + 2*Ynp1)*(Rn^2 - 2*Rn*Rnp1 + Rnp1^2 - Xn^2 + 2*Xn*Xnp1 - Xnp1^2 - Yn^2 + 2*Yn*Ynp1 - Ynp1^2))*(Xn*Ynp1 - Xn - Xnp1*Yn + Xnp1)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)) - (Rn^3*Xnp1^2 + Rn^3*Ynp1^2 - 2*Rn^3*Ynp1 + Rn^3 - Rn^2*Rnp1*Xn*Xnp1 - Rn^2*Rnp1*Yn*Ynp1 + Rn^2*Rnp1*Yn + Rn^2*Rnp1*Ynp1 - Rn^2*Rnp1 + Rn^2*Xn*Xnp1 - Rn^2*Xnp1^2 + Rn^2*Yn*Ynp1 - Rn^2*Yn - Rn^2*Ynp1^2 + Rn^2*Ynp1 - Rn*Rnp1^2*Xn*Xnp1 - Rn*Rnp1^2*Yn*Ynp1 + Rn*Rnp1^2*Yn + Rn*Rnp1^2*Ynp1 - Rn*Rnp1^2 - Rn*Xn^2*Xnp1^2 - Rn*Xn^2*Ynp1^2 + 2*Rn*Xn^2*Ynp1 - Rn*Xn^2 + Rn*Xn*Xnp1^3 + Rn*Xn*Xnp1*Ynp1^2 - 2*Rn*Xn*Xnp1*Ynp1 + 2*Rn*Xn*Xnp1 - Rn*Xnp1^2*Yn^2 + Rn*Xnp1^2*Yn*Ynp1 + Rn*Xnp1^2*Yn - Rn*Xnp1^2*Ynp1 - Rn*Xnp1^2 - Rn*Yn^2*Ynp1^2 + 2*Rn*Yn^2*Ynp1 - Rn*Yn^2 + Rn*Yn*Ynp1^3 - Rn*Yn*Ynp1^2 - Rn*Ynp1^3 + Rn*Ynp1^2 + Rnp1^3*Xn^2 + Rnp1^3*Yn^2 - 2*Rnp1^3*Yn + Rnp1^3 - Rnp1^2*Xn^2 + Rnp1^2*Xn*Xnp1 - Rnp1^2*Yn^2 + Rnp1^2*Yn*Ynp1 + Rnp1^2*Yn - Rnp1^2*Ynp1 + Rnp1*Xn^3*Xnp1 - Rnp1*Xn^2*Xnp1^2 + Rnp1*Xn^2*Yn*Ynp1 - Rnp1*Xn^2*Yn - Rnp1*Xn^2*Ynp1^2 + Rnp1*Xn^2*Ynp1 - Rnp1*Xn^2 + Rnp1*Xn*Xnp1*Yn^2 - 2*Rnp1*Xn*Xnp1*Yn + 2*Rnp1*Xn*Xnp1 - Rnp1*Xnp1^2*Yn^2 + 2*Rnp1*Xnp1^2*Yn - Rnp1*Xnp1^2 + Rnp1*Yn^3*Ynp1 - Rnp1*Yn^3 - Rnp1*Yn^2*Ynp1^2 - Rnp1*Yn^2*Ynp1 + Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1^2 - Rnp1*Ynp1^2 - Xn^3*Xnp1 + 2*Xn^2*Xnp1^2 - Xn^2*Yn*Ynp1 + Xn^2*Yn + Xn^2*Ynp1 - Xn*Xnp1^3 - Xn*Xnp1*Yn^2 + 4*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - Xn*Xnp1*Ynp1^2 - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn*Ynp1 + Xnp1^2*Yn + Xnp1^2*Ynp1 - Yn^3*Ynp1 + Yn^3 + 2*Yn^2*Ynp1^2 - Yn^2*Ynp1 - Yn*Ynp1^3 - Yn*Ynp1^2 + Ynp1^3)/(2*(Rn^2*Xnp1^2 + Rn^2*Ynp1^2 - 2*Rn^2*Ynp1 + Rn^2 - 2*Rn*Rnp1*Xn*Xnp1 - 2*Rn*Rnp1*Yn*Ynp1 + 2*Rn*Rnp1*Yn + 2*Rn*Rnp1*Ynp1 - 2*Rn*Rnp1 + 2*Rn*Xn*Xnp1 - 2*Rn*Xnp1^2 + 2*Rn*Yn*Ynp1 - 2*Rn*Yn - 2*Rn*Ynp1^2 + 2*Rn*Ynp1 + Rnp1^2*Xn^2 + Rnp1^2*Yn^2 - 2*Rnp1^2*Yn + Rnp1^2 - 2*Rnp1*Xn^2 + 2*Rnp1*Xn*Xnp1 - 2*Rnp1*Yn^2 + 2*Rnp1*Yn*Ynp1 + 2*Rnp1*Yn - 2*Rnp1*Ynp1 - Xn^2*Ynp1^2 + 2*Xn^2*Ynp1 + 2*Xn*Xnp1*Yn*Ynp1 - 2*Xn*Xnp1*Yn - 2*Xn*Xnp1*Ynp1 - Xnp1^2*Yn^2 + 2*Xnp1^2*Yn + Yn^2 - 2*Yn*Ynp1 + Ynp1^2)));
end

   乙 (1.33333, 2.00000, 0.66667)
   丙 (1.33333, 1.00000, 0.33333)
   丁 (1.09091, 0.54545, 0.18182)
   戊 (0.88889, 0.33333, 0.11111)
   己 (0.74074, 0.22222, 0.07407)
   子 (0.53333, 2.00000, 0.13333) 元の単位で 5.06667 寸
   丑 (1.04348, 1.30435, 0.08696) 元の単位で 3.30435 寸
   寅 (1.02564, 0.76923, 0.05128) 元の単位で 1.94872 寸
   卯 (0.88889, 0.47619, 0.03175) 元の単位で 1.20635 寸
   辰 (0.75789, 0.31579, 0.02105) 元の単位で 0.80000 寸

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R0 = 2  # 大円の半径
   R1 = R0 // 2  # 甲円の半径
   X2, Y2, R2 = (4//3, 2//1, 2//3)
   XnA = zeros(6)
   YnA = zeros(6)
   RnA = zeros(6)
   plot()
   circle(0, R0, R0, :black)
   circle(0, R1, R1, :blue)
   circle(0, 3R1, R1, :green)
   circle(X2, Y2, R2, :green)
   circle(-X2, Y2, R2, :green)
   (XnA[1], YnA[1], RnA[1]) = (Xn, Yn, Rn) = (X2, Y2, R2)
   large = Char["乙丙丁戊己　"...]
   for i = 1:5
       @printf("%s (%.5f, %.5f, %.5f)\n", large[i], Xn, Yn, Rn)
       circle(Xn, Yn, Rn, :green)
       circle(-Xn, Yn, Rn, :green)
       more && point(Xn, Yn, "  " * large[i], mark=false)
       i == 5 && break
       (XnA[i+1], YnA[i+1], RnA[i+1]) = (Xn, Yn, Rn) = NEXTCIRCLE(Xn, Yn, Rn)
   end
   (xn, yn, rn) = x2.evalf(), y2.evalf(), r2.evalf()
   circle(xn, yn, rn, :blue, fill=true)
   circle(-xn, yn, rn, :blue, fill=true)
   more && point(xn, yn, "子　 ", :blue, :right, mark=false)
   small = Char["子丑寅卯辰"...]
   @printf("%s (%.5f, %.5f, %.5f) 元の単位で %.5f 寸\n", small[1], xn, yn, rn, 2rn * 76 / 4)
   for i = 1:4
       (xn, yn, rn) = nextcircle(XnA[i], YnA[i], RnA[i], XnA[i+1], YnA[i+1], RnA[i+1])
       @printf("%s (%.5f, %.5f, %.5f) 元の単位で %.5f 寸\n", small[i+1], xn, yn, rn, 2rn * 76 / 4)
       circle(xn, yn, rn, :blue, fill=true)
       circle(-xn, yn, rn, :blue, fill=true)
       more && i < 5 && point(xn, yn, small[i+1] * "　", :blue, :right, :bottom, mark=false)
   end
   if more == true
       point(0, 0, "O ", :green, :right, :bottom)
       point(0, R1, "R1 ", :green, :right, :bottom)
       point(0, R0, "R0=2R1", :green, :center, :bottom)
       point(0, R0+R1, "3R1 ", :green, :right, :bottom)
       point(0, 2R0, "2R0 ", :green, :right, :bottom)
       point(0, 3R1, " 甲", :green, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

 

算額（その193）

算額（その193）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（88）
長野県中野市田上 田上観音堂 文化6年(1809)

問2. 大円の中に甲，乙，丙，丁を入れる。それぞれは互いに隣同士の円は外接し，大円に内接する正三角形の底辺及び大円に接している。大円の径が 196 寸，丙円の径が 17寸6分5厘のとき，丁円の径はいかほどか。

大，甲，乙，丙，丁の円の半径を R, r1, r2, r3, r4，それぞれの円の中心座標を (0, 0), (0, -R/2-r1), (x2, -R/2-r2), (x3, -R/2-r3), (x4, -R/2-r4) と置く。ただし，R = 19600, r3 = 1765 とする。

using SymPy

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

R = 19600
r3 = 1764

eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq3 = (x4 - x3)^2 + (r3 - r4)^2 - (r3 + r4)^2
eq4 = x2^2 + (R//2 + r2)^2 - (R - r2)^2
eq5 = x3^2 + (R//2 + r3)^2 - (R - r3)^2
eq6 = x4^2 + (R//2 + r4)^2 - (R - r4)^2

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

   4-element Vector{NTuple{6, Sym}}:
    (4900, 3675, 675, 4900*sqrt(3), 7840*sqrt(3), 9100*sqrt(3))
    (4900, 3675, 3675, 4900*sqrt(3), 7840*sqrt(3), 4900*sqrt(3))
    (828100/9, 675, 675, 9100*sqrt(3), 7840*sqrt(3), 9100*sqrt(3))
    (828100/9, 675, 3675, 9100*sqrt(3), 7840*sqrt(3), 4900*sqrt(3))

r1 = 4900.00000, r2 = 3675.00000, r4 = 675.00000, x2 = 8487.04896, x3 = 13579.27833, x4 = 15761.66235

r2 > r4 であるから，1 番目の解が適切である。
丙円の半径は r4 = 675，もとの単位でいえば丙円の直径は 6寸7分5厘 である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r4, x2, x3, x4) = res[1]
   @printf("r1 = %.5f, r2 = %.5f, r4 = %.5f, x2 = %.5f, x3 = %.5f, x4 = %.5f\n", r1, r2, r4, x2, x3, x4)
   plot([√3R/2, 0, -√3R/2, √3R/2], [-R/2, R, -R/2, -R/2], color=:blue, lw=0.5)
   circle(0, 0, R, :black)
   circle(0, -R/2 - r1, r1, :blue)
   circle(x2, -R/2 - r2, r2, :green)
   circle(-x2, -R/2 - r2, r2, :green)
   circle(x3, -R/2 - r3, r3, :brown)
   circle(-x3, -R/2 - r3, r3, :brown)
   circle(x4, -R/2 - r4, r4, :red)
   circle(-x4, -R/2 - r4, r4, :red)
   if more == true
       point(0, -R, " -R")
       point(0, -R/2, " -R/2", :green, :left, :bottom) 
       point(0, -R/2-r1, " 甲")
       point(x2, -R/2-r2, " 乙")
       point(x3, -R/2-r3, " 丙")
       point(x4, -R/2-r4, " 丁")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

算額（その192）

算額（その192）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（88）
長野県中野市田上 田上観音堂 文化6年(1809)

問1. 長方形内に 2 本の斜線を入れ，それぞれに接する甲円と乙円を入れる。長方形の短辺の長さが 10 寸，甲円の直径が 8 寸のとき，長方形の長辺の長さはいかほどか。

長方形の長辺の長さを x，乙円の半径を r，斜線が x 軸，y軸と交わる座標を x1, y1 とする。各円の中心から斜線までの距離が半径に等しいことを条件として，方程式を解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

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

eq1 = distance(0, y1, x, 10, r, r) - r^2 |> simplify
eq2 = distance(0, 10, x1, 0, r, r) - r^2 |> simplify
eq3 = distance(0, y1, x, 10, x - 4, 4) - 4^2 |> simplify
eq4 = distance(0, 10, x1, 0, x - 4, 4) - 4^2 |> simplify

res = solve([eq1, eq2, eq3, eq4], (x, x1, y1, r))

   1-element Vector{NTuple{4, Sym}}:
    (8, 24, 20/3, 4)

今までに経験のないことであるが，solve で得られる解が不適切解である。得られる乙円の半径が 4 という解は甲円と乙円が重なるという解である。

そこで，nlsolve() で解いてみる。

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

   x*(2*r^2*y1 - 20*r^2 - 2*r*x*y1 - 2*r*y1^2 + 20*r*y1 + x*y1^2)/(x^2 + y1^2 - 20*y1 + 100),
   20*x1*(r^2 - r*x1 - 10*r + 5*x1)/(x1^2 + 100),
   4*x*(5*x + 12*y1 - 120)/(x^2 + y1^2 - 20*y1 + 100),
   20*(5*x^2 - 6*x*x1 - 40*x + x1^2 + 24*x1)/(x1^2 + 100),

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, x1, y1, r) = u
   return [
       x*(2*r^2*y1 - 20*r^2 - 2*r*x*y1 - 2*r*y1^2 + 20*r*y1 + x*y1^2)/(x^2 + y1^2 - 20*y1 + 100),
       20*x1*(r^2 - r*x1 - 10*r + 5*x1)/(x1^2 + 100),
       4*x*(5*x + 12*y1 - 120)/(x^2 + y1^2 - 20*y1 + 100),
       20*(5*x^2 - 6*x*x1 - 40*x + x1^2 + 24*x1)/(x1^2 + 100),
   ]
end;

iniv = [big"13.8", 7.8, 4.3, 2.6]
res = nls(H, ini=iniv)
println(res)
   (BigFloat[13.75900072948533162286601728523747073950320887218339386802707109168781010348274, 7.807200583588265253825888517990309019730833318263487969214336156987235887470571, 4.267083029381111827960970371439448927535572362184485605389630903733447101473092, 2.560249817628667040685677767677303452771782038991463293198420137207969259035087], true)

x = 13.75900, x1 = 7.80720, y1 = 4.26708, r = 2.56025
であり，x = 13.75900 は元の単位でいえば 13寸7分5厘9毛 である。算額の答えでは 13寸7分5厘6毛あまりあり となっている。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, x1, y1, r) = res[1]
   @printf("x = %.5f, x1 = %.5f, y1 = %.5f, r = %.5f\n", x, x1, y1, r)
   plot([0, x, x, 0, 0], [0, 0, 10, 10, 0], color=:blue, lw=0.5)
   circle(r, r, r)
   circle(x - 4, 4, 4, :blue)
   segment(0, 10, x1, 0, :green)
   segment(0, y1, x, 10, :green)
   if more == true
       point(x - 4, 4, " 甲: (x-4,4)")
       point(r, r, " 乙: (r,r)")
       point(x1, 0, "x1")
       point(0, y1, "y1 ", :green, :right)
       point(x, 0, " x")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5, xlims=(-1, 14.5),ylims=(-0.5, 10.5))
   else
      plot!(showaxis=false)
   end
end;

算額（その191）

算額（その191）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（88）
長野県中野市田上 田上観音堂 文化6年(1809)

大円内に，天斜，地斜，左斜，右斜で構成される四辺形及びその対角線（乾斜，坤斜）を描く。天斜，乾斜，坤斜で囲まれる領域に小円を描く。
小円と大円の直径の（値の）積が 780 寸（本当は寸^2=歩？），また，天斜，左斜，右斜の長さがそれぞれ 39 寸，60 寸，25 寸のとき，地斜の長さはいかほどか。

四辺形の頂点を (x1, y1), (x2, y2), (x3, y3), (x4, y4)，小円の中心座標と半径を (x5, y5), r5 とおく。そのままでは未知数が 12 個，条件が 11 個で解けないが，図は回転しても一般性を失わないので，x1 = 0 になるように座標軸を定めれば方程式が解ける。また， x1 = 0 なら， y1 は大円の円周上になるので，未知数も方程式も 1 個ずつ減る。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, r5, R;

x1 = 0
y1 = -R
eq1 = R * r5 - 780//4
eq2 = x2^2 + y2^2 - R^2
eq3 = x3^2 + y3^2 - R^2
eq4 = x4^2 + y4^2 - R^2
eq5 = (x4 - x1)^2 + (y4 - y1)^2 - 39^2
eq6 = (x4 - x3)^2 + (y4 - y3)^2 - 60^2
eq7 = (x1 - x2)^2 + (y1 - y2)^2 - 25^2
eq8 = distance(x4, y4, x1, y1, x5, y5) - r5^2 |> simplify
eq9 = distance(x4, y4, x2, y2, x5, y5) - r5^2 |> simplify
eq10 = distance(x3, y3, x1, y1, x5, y5) - r5^2 |> simplify;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10])

nlsolve() で解く。

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

   R*r5 - 195,
   -R^2 + x2^2 + y2^2,
   -R^2 + x3^2 + y3^2,
   -R^2 + x4^2 + y4^2,
   x4^2 + (R + y4)^2 - 1521,
   (-x3 + x4)^2 + (-y3 + y4)^2 - 3600,
   x2^2 + (-R - y2)^2 - 625,
   (-r5^2*(R^2 + 2*R*y4 + x4^2 + y4^2)^2 + x4^2*(R*x4 - R*x5 + x4*y5 - x5*y4)^2 + (x4*(R^2 + R*y4 + R*y5 + x4*x5 + y4*y5) - x5*(R^2 + 2*R*y4 + x4^2 + y4^2))^2)/(R^2 + 2*R*y4 + x4^2 + y4^2)^2,
   (-r5^2*(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2 + (x2 - x4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2 + (y2 - y4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2)/(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2,
   (-R^2*r5^2 + R^2*x3^2 - 2*R^2*x3*x5 + R^2*x5^2 - 2*R*r5^2*y3 + 2*R*x3^2*y5 - 2*R*x3*x5*y3 - 2*R*x3*x5*y5 + 2*R*x5^2*y3 - r5^2*x3^2 - r5^2*y3^2 + x3^2*y5^2 - 2*x3*x5*y3*y5 + x5^2*y3^2)/(R^2 + 2*R*y3 + x3^2 + y3^2),

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x2, y2, x3, y3, x4, y4, x5, y5, r5, R) = u
   return [
R*r5 - 195,
-R^2 + x2^2 + y2^2,
-R^2 + x3^2 + y3^2,
-R^2 + x4^2 + y4^2,
x4^2 + (R + y4)^2 - 1521,
(-x3 + x4)^2 + (-y3 + y4)^2 - 3600,
x2^2 + (-R - y2)^2 - 625,
(-r5^2*(R^2 + 2*R*y4 + x4^2 + y4^2)^2 + x4^2*(R*x4 - R*x5 + x4*y5 - x5*y4)^2 + (x4*(R^2 + R*y4 + R*y5 + x4*x5 + y4*y5) - x5*(R^2 + 2*R*y4 + x4^2 + y4^2))^2)/(R^2 + 2*R*y4 + x4^2 + y4^2)^2,
(-r5^2*(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2 + (x2 - x4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2 + (y2 - y4)^2*(x2*y4 - x2*y5 - x4*y2 + x4*y5 + x5*y2 - x5*y4)^2)/(x2^2 - 2*x2*x4 + x4^2 + y2^2 - 2*y2*y4 + y4^2)^2,
(-R^2*r5^2 + R^2*x3^2 - 2*R^2*x3*x5 + R^2*x5^2 - 2*R*r5^2*y3 + 2*R*x3^2*y5 - 2*R*x3*x5*y3 - 2*R*x3*x5*y5 + 2*R*x5^2*y3 - r5^2*x3^2 - r5^2*y3^2 + x3^2*y5^2 - 2*x3*x5*y3*y5 + x5^2*y3^2)/(R^2 + 2*R*y3 + x3^2 + y3^2),
   ]
end;

iniv = [big"-23.0", -22, -15, 28, 31, -9, 3, -22, 6, 32]
res = nls(H, ini=iniv)
println(res)
   (BigFloat[-23.07692307692307692307692307691865904585941003474930311365314461268297834238257, -22.88461538461538461538461538450049187950077181984566989711738965555646289279357, -15.50769230769230769230769230792675686699338875941266776383234390713069575421587, 28.5615384615384615384615384616459508480950144867108280676139955837025986961933, 31.19999999999999999999999999998122624448363949495361597070768953084622475373686, -9.099999999999999999999999999863248885904729542964274601211974716937721151767532, 3.599999999999999999999999999995741859144646799042250528533190648410558984504582, -22.29999999999999999999999999988215259415310464508394927693531668395582734424483, 6.000000000000000000000000000009441163418536999038475318876384580066842683370904, 32.49999999999999999999999999990015502464806610960030256664669973254443033639196], true)

x2 = -23.07692, y2 = -22.88462, x3 = -15.50769, y3 = 28.56154
x4 = 31.20000, y4 = -9.10000, x5 = 3.60000, y5 = -22.30000, r5 = 6.00000, R = 32.50000
となり，求める地斜の長さは三平方の定理より 52 と計算される。
地斜 = sqrt((x3 - x2)^2 + (y3 - y2)^2) = 52.00000

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360)
  θ = beginangle:0.1:endangle
  x = r.*cosd.(θ)
  y = r.*sind.(θ)
  plot!(ox .+ x, oy .+ y, color=color, linewidth=0.5)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x2, y2, x3, y3, x4, y4, x5, y5, r5, R) = res[1]
   x1 = 0
   y1 = -R
   @printf("x1 = %.5f, y1 = %.5f, x2 = %.5f, y2 = %.5f, x3 = %.5f, y3 = %.5f\nx4 = %.5f, y4 = %.5f, x5 = %.5f, y5 = %.5f, r5 = %.5f, R = %.5f\n", x1, y1, x2, y2, x3, y3, x4, y4, x5, y5, r5, R)
   @printf("地斜 = sqrt((x3 - x2)^2 + (y3 - y2)^2) = %.5f\n", sqrt((x3 - x2)^2 + (y3 - y2)^2))
   plot([x1, x2, x3, x4, x1], [y1, y2, y3, y4, y1], color=:blue, lw=0.5)
   circle(0, 0, R)
   circle(x5, y5, r5)
   segment(x1, y1, x3, y3, :green)
   segment(x2, y2, x4, y4, :green)
   if more == true
       point(x1, y1, "(x1,y1)")
       point(x2, y2, "(x2,y2) ", :green, :right)
       point(x3, y3, "  (x3,y3)")
       point(x4, y4, " (x4,y4)")
       point(x5, y5, "(x5,y5;r5)", :green, :center)
       point(x5, y5, "小円", :red, :center, :bottom)
       point(R, 0, " R: 大円", :red)
       point((x3 + x4)/2, (y3 + y4)/2, "   左斜", :blue, mark=false)
       point((x2 + x3)/2, (y2 + y3)/2, " 地斜", :blue, mark=false)
       point((x1 + x2)/2, (y1 + y2)/2, "右斜", :blue, :right, mark=false)
       point((x1 + x4)/2, (y1 + y4)/2, "天斜", :blue, mark=false)
       point((x1 + x3)/2, (y1 + y3)/2, " 乾斜", :green, mark=false)
       point((x2 + x4)/2, (y2 + y4)/2, "    坤斜", :green, mark=false)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

   x1 = 0.00000, y1 = -32.50000, x2 = -23.07692, y2 = -22.88462, x3 = -15.50769, y3 = 28.56154
   x4 = 31.20000, y4 = -9.10000, x5 = 3.60000, y5 = -22.30000, r5 = 6.00000, R = 32.50000
   地斜 = sqrt((x3 - x2)^2 + (y3 - y2)^2) = 52.00000

算額（その190）

算額（その190）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（87）
長野県中野市田上 田上観音堂 文化6年(1809)
県内の算額1（64）
長野県下高井郡木島平村往郷 水穂神社 寛政12年(1800)
県内の算額1（46）
長野県長野市若穂 清水寺観音堂 寛政6年(1794)

第6問 鉤股の中に，全，甲，乙，丙，丁の正方形を入れる。丙，丁の一辺が 192 寸, 108 寸 のとき，全の一辺はいかほどか。

全，甲，乙，丙，丁 の正方形の一辺を a, b, c, d, e とする。c = 192, e = 108 で，構成される複数の直角三角形の鉤股の比を t =　AC / AB = (a + d + e - c)/(a + b + c - e) として，連立方程式を解く。

using SymPy

@syms a::positive, b::positive, d::positive;
c = 192
e = 108
t = (a + d + e - c)/(a + b + c - e)
eq1 = t - (b - c)/c
eq2 = t - e/(d - e)
eq3 = t - (a - c)/(b + c)

res = solve([eq1, eq2, eq3], (a, b, d))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (588, 336, 252)

全，甲，乙の一辺の長さはそれぞれ 588 寸，336 寸，252 寸である。

ちなみに，鉤股弦は 3:4:5, t = 3/4

using Plots
using Printf

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function rect(x1, y1, x2, y2, color=:pink)
   plot!([x1, x2, x2, x1, x1], [y1, y1, y2,  y2, y1], color=color, linewidth=0.5, seriestype=:shape, fillcolor=color)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, d) = res[1]
   (c, e) = (192, 108)
   t = (a + d + e - c)/(a + b + c - e)
   println("t = $t")
   plot()
   rect(0, 0, a, a)
   rect(a, 0, a + b, b)
   rect(a + b, 0, a + b + c, c)
   rect(0, a, d, a + d)
   rect(0, a + d, e, a + d + e)
   plot!([0, a + b + c + c/t, 0, 0], [0, 0, a + d + e + e*t, 0], color=:black, lw=0.5)
   if more == true
       point(a, 0, " a", :black, :left, :top)
       point(a + b, 0, " a+b", :black, :left, :top)
       point(a + b + c, 0, " a+b+c", :black, :left, :top)
       point(0, a, "a ", :black, :right, :bottom)
       point(0, a + d, "a+d ", :black, :right, :bottom)
       point(0, a + d + e, "a+d+e ", :black, :right, :bottom)
       plot!([e, a + b + c, e, e], [c, c, a + d + e, c], color=:red, lw=1)
       point(e, c, " A", :black, :left, :bottom)
       point(a + b + c, c, " B", :black, :left, :bottom)
       point(e, a + d + e, " C", :black, :left, :bottom)
       hline!([0], color=:black, lw=0.5, xlims=(-175, 1400))
       vline!([0], color=:black, lw=0.5, ylims=(-70, 1050))
   else
      plot!(showaxis=false)
   end
end;

   t = 3/4

 

算額（その189）

算額（その189）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（86）
長野県中野市田上 田上観音堂 文化6年(1809)

第5問 内外2個の菱形により決まる黒の部分の面積が 144 寸（本当は平方寸）
外側，内側の菱形の辺の長さが 17 寸，10 寸のとき共通対角線の長さはいかほどか。

内側の菱形の OB，OA の長さを a，b とする。BC の長さを c とする。

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

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive
eq1 = a^2 + b^2 - 10^2  # ⊿ AOB についてピタゴラスの定理
eq2 = (a + c)^2 + b^2 - 17^2  # ⊿ AOC についてピタゴラスの定理
eq3 = 2b*(a + c) - 2a*b - 144  # 外側の菱形の面積から内側のひし形の面積を引いたものが黒積

solve([eq1, eq2, eq3], (a, b, c))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (6, 8, 9)

共通対角線の長さは 2b = 16

using Plots
using Printf

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, b, c) = (6, 8, 9)
   plot([15, 0, -a - c, 0, a + c], [0, b, 0, -b, 0], linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=:black)
   plot!([a, 0, -a, 0, a], [0, b, 0, -b, 0], linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=:white)
   if more == true
       point(0, 0, " O", :green, :left, :bottom)
       point(a, 0, " B: a", :cyan, :left, :bottom)
       point(a + c, 0, " C: a+c", :red, :left, :bottom)
       point(0, b, " A: b", :green, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       segment(a + c, 0, a, 0, :gray, linewidth=1)
       segment(-a - c, 0, -a, 0, :gray, linewidth=1)
   else
      plot!(showaxis=false)
   end
end;

算額（その188）

算額（その188）

中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（86）
長野県中野市田上 田上観音堂 文化6年(1809)

第4問 鉤股（直角三角形）内に大中小の3円を入れる。大円と中円の径を掛けると 72，大円と小円の径を掛けると 48 であるとき，それぞれの径を求めよ。
注：算額の図は「算額（その187）https://blog.goo.ne.jp/r-de-r/e/2053d2ea5ef9100ee13ecf80457132b9」と違い垂直線，水平線に接するという条件がない場合の解である。

大円，中円，小円の半径を r1,r2,r3 とする。中円，小円の中心座標を (x2, r2), (r3, y3) とする。
中円，小円の中心から弦への距離，鉤股弦の面積の関係式を解く。

using SymPy

function distance(x1, y1, x2, y2, x0, y0)
   p1, p2 = sympy.Point(x1, y1), sympy.Point(x2, y2)
   l = sympy.Line(p1, p2)
   l.distance(sympy.Point(x0, y0))^2  # 二乗距離を返す！！
end;

@syms X::positive, Y::positive, x2::pisitive, y3::positive,
     r1::positive, r2::positive, r3::positive, l::positive;

eq1 = 2r1 * 2r2 - 72
eq2 = 2r1 * 2r3 - 48
eq3 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq5 = distance(0, Y, X, 0, r3, y3) - r3^2
eq6 = distance(0, Y, X, 0, x2, r2) - r2^2
eq7 = X^2 + Y^2 - l^2
eq8 = (X + Y + l)*r1 - X * Y;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])

solve() では解けないので，nlsolve() を使用する。

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

   4*r1*r2 - 72,
   4*r1*r3 - 48,
   (-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2,
   (-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2,
   -r3^2 + (-X*(X*r3 + Y^2 - Y*y3)/(X^2 + Y^2) + r3)^2 + (-Y*(X^2 - X*r3 + Y*y3)/(X^2 + Y^2) + y3)^2,
   -r2^2 + (-X*(X*x2 + Y^2 - Y*r2)/(X^2 + Y^2) + x2)^2 + (-Y*(X^2 - X*x2 + Y*r2)/(X^2 + Y^2) + r2)^2,
   X^2 + Y^2 - l^2,
   -X*Y + r1*(X + Y + l),

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r1, r2, r3, X, Y, l, x2, y3) = u
   return [
4*r1*r2 - 72,
4*r1*r3 - 48,
(-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2,
(-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2,
-r3^2 + (-X*(X*r3 + Y^2 - Y*y3)/(X^2 + Y^2) + r3)^2 + (-Y*(X^2 - X*r3 + Y*y3)/(X^2 + Y^2) + y3)^2,
-r2^2 + (-X*(X*x2 + Y^2 - Y*r2)/(X^2 + Y^2) + x2)^2 + (-Y*(X^2 - X*x2 + Y*r2)/(X^2 + Y^2) + r2)^2,
X^2 + Y^2 - l^2,
-X*Y + r1*(X + Y + l),        
      ]
end;

iniv = [big"6.1", 4.2, 3.3, 24, 17, 31, 15, 13]
res = nls(H, ini=iniv)
println(res)

   (BigFloat[5.748772323345360879553386071530534002107771135435234241146154053030068382725282, 3.131103301291523318846022596310963240657188641774061152337086444494206214742888, 2.08740220086101554589734764786891092285620654321669230693296468944296425038439, 24.38365332206977707890963182337988451889403111451717974855438930320079233053448, 16.62684846651738071918895370581186617641229049594184966902009096464686858816747, 29.51295714189643603899180810391979664284976102613029457759886594916383321385045, 14.23405369758393117236351841678872247352580238769692473134187249862639602572867, 12.67697555362087005366317143755402346987899215067675675336938195307101120791379], true)

using Printf
@printf("大円径 = %6.3f; 中円径 = %6.3f; 小円径 = %6.3f\n", 2res[1][1], 2res[1][2], 2res[1][3])

   大円径 = 11.498; 中円径 =  6.262; 小円径 =  4.175

大円，中円，小円の径は元の単位で，11.498 寸，6.262 寸，4.175 寸である。

using Plots
using Printf

function circle(ox, oy, r, color=:red; beginangle=0, endangle=360, fill=false)
   θ = beginangle:0.1:endangle
   x = r.*cosd.(θ)
   y = r.*sind.(θ)
   if fill
       plot!(ox .+ x, oy .+ y, linecolor=:black, linewidth=0.5, seriestype=:shape, fillcolor=color)
   else
       plot!(ox .+ x, oy .+ y, color=color, linewidth=0.25)
   end
end;

function point(x, y, string="", color=:green, position=:left, vertical=:top; mark=true)
   mark && scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function segment(x1, y1, x2, y2, color=:black; linestyle=:solid, linewidth=0.5)
   plot!([x1, x2], [y1, y2], color=color, linestyle=linestyle, linewidth=linewidth)
end;

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, X, Y, l, x2, y3) = res[1]
   @printf("X = %6.3f; Y = %6.3f; l = %6.3f; x2 = %6.3f; y3 = %6.3f\n",
       X, Y, l, x2, y3)
   @printf("r1 = %6.3f; r2 = %6.3f; r3 = %6.3f\n", r1, r2, r3)
   plot([0, X, 0, 0], [0, 0, Y, 0], color=:black, lw=0.5)
   circle(r1, r1, r1, :red)
   circle(x2, r2, r2, :blue)
   circle(r3, y3, r3, :green)
   if more == true
       point(r1, r1, "大円(r1,r1,r1)")
       point(x2, r2, "中円(x2,r2,r2)")
       point(r3, y3, "小円(r3,y3,r3)")
       point(0, 0, " O", :green, :left, :bottom)
       point(X, 0, " X", :green, :left, :bottom)
       point(0, Y, " Y", :green, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
      plot!(showaxis=false)
   end
end;

   X = 24.384; Y = 16.627; l = 29.513; x2 = 14.234; y3 = 12.677
   r1 =  5.749; r2 =  3.131; r3 =  2.087