算額（その15）

算額（その15）

岩手県東磐井郡川崎村 浪分神社 年代不詳（文久か）
http://www.wasan.jp/iwate/namiwake5.html

岩手県陸前高田市気仙町 今泉諏訪神社慶応 4 年
http://www.wasan.jp/iwate/suwa.html

下図のように，半径 1 の円の中に 2 種類の円がありそれぞれ接している。それぞれの円の半径を求めよ。

紙と鉛筆でも求めることができる。

小さい円の半径は r2 = 1/5 である。大きい円の半径は r1 = 3*r2。右の小さい円の中心を (x, ±r2) とする。

using SymPy

@syms x::positive;

r2 = 1//5
r1 = 3r2
eq1 = x^2 + r2^2 - (1 - r2)^2
eq1 |> expand |> println

   x^2 - 3/5

solve(eq1, x)[1] |> println # x について解く

   sqrt(15)/5

using Plots

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=:center)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   x = sqrt(15)/5
   r2 = 1/5
   r1 = 3r2
   println("r2 = $r2,  x = $x")
   plot()
   circle(0,  2r2, 3r2, :blue)
   circle(0, -2r2, 3r2, :blue)
   circle( x,  r2, r2, :green)
   circle( x, -r2, r2, :green)
   circle(-x,  r2, r2, :green)
   circle(-x, -r2, r2, :green)
   circle( 0,   0, r2, :green)
   if more
       hline!(-1:r2:1, linewidth=0.1)
       point(x, r2, "(x,r2)\n", :green, :center, :bottom)
       plot!([0, x], [0, r2], linewidth=0.5, linestyle=:dash)
   end
   circle(0, 0, 1)
end;

算額（その14）

算額（その14）

山形県山形市山寺 立石寺（山寺）根本中堂 大正3年
http://www.wasan.jp/yamagata/yamadera.html

下図のように，半径 1 の円の中に 3 つの正方形がある。一番小さい正方形の一辺の長さを求めよ。

最小の正方形の一辺を a，二番目に小さい正方形の一辺をその x 倍，すなわち ax とする。

using SymPy

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

sqrt2 = sqrt(Sym(2));
eq1 = (a*x + a*sqrt2)^2 +(a*x + a/sqrt2)^2 - 1;
eq2 = (a*x)^2 + (3a*x)^2 - 1;

連立方程式を解き，r, x を求める。

res = solve([eq1, eq2], (a, x))

   1-element Vector{Tuple{Sym, Sym}}:
    (4*sqrt(5)/25, 5*sqrt(2)/8)

using Plots

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=:center)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function square(x1, y1, x2, y2, color=:green)
   plot!([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1],
         linewidth=0.5, color=color)
end

function square2(a, x, color=:magenta)
   plot!([a*(x+1/sqrt(2)), a*(x+sqrt(2)), a*(x+1/sqrt(2)), a*x, a*(x+1/sqrt(2))],
         [a*x, a*(x+1/sqrt(2)), a*(x + sqrt(2)), a*(x+1/sqrt(2)), a*x],
         linewidth=0.5, color=color)
end

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   (a, x) = (4*sqrt(5)/25, 5*sqrt(2)/8)
   println("a = $a, x = $x")
   plot()
   if more
       point(0, a*x, "(0,a*x)", :blue, :right, :top)
       point(a*x, a*x, "(a*x,a*x)", :blue, :center, :top)
       point(a*x, a*x+a/√2, "(a*x,a*x+a/√2)", :blue, :right, :top)
       point(a*x+a/sqrt(2), a*x, "(a*x+a/√2)", :blue, :left, :top)
       point(0, -a*x-a*sqrt(2), "(0,-a*x-a*√2)", :blue, :center, :top)
       point(0, 3a*x, "(0,3a*x)", :blue, :center, :top)
   end
   square(-a*x - a/sqrt(2), -a*x - a*sqrt(2), a*x + a/sqrt(2), a*x, :green)
   square(-a*x, a*x, a*x, 3a*x, :brown)
   square2(a, x)
   circle(0, 0, 1)
end;

算額（その13）

算額（その13）

山形県山形市山寺 立石寺（山寺）根本中堂 大正3年
http://www.wasan.jp/yamagata/yamadera.html

下図のように，正方形の中に四分円と円がある。円はそれぞれ四分円と正方形の対角線に接する。円の半径を求めよ。

正方形の一辺を 1，円の半径と中心座標を r, (x, r - 1) とする。

using SymPy

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

四分円と円が接することから

eq1 = x^2 + (r - 1)^2 - (1 + r)^2 |> expand
eq1 |> println

   -4*r + x^2

円と正方形の対角線が接することから

x2 = x - r/sqrt(Sym(2))
y2 = r - 1 + r/sqrt(Sym(2))
eq2 = r^2 + x2^2 + y2^2 - (1 + r)^2 |> expand
println(eq2)

   sqrt(2)*r^2 + 2*r^2 - sqrt(2)*r*x - 4*r - sqrt(2)*r + x^2

連立方程式を解き，r, x を求める。

res = solve([eq1, eq2], (r, x))

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

using Plots

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=:center)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color, vertical))
end;

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   (r, x) = ((-2 - 2*sqrt(2 - sqrt(2)) + 2*sqrt(2))^2/4, -2 - 2*sqrt(2 - sqrt(2)) + 2*sqrt(2))
   println("r = $r,  x = $x")
   x2 = x - r/sqrt(2)
   y2 = r - 1 + r/sqrt(2)
   plot()
   if more
       point(x, r - 1, "(x,r-1)", :blue, :left, :top)
       point(x2, y2, "(x2,y2)", :blue, :right, :bottom)
       plot!([0, x], [0, r - 1], color=:red, linewidth=0.5, linestyle=:dash)
   end

   circle(0, 0, 1, beginangle=180, endangle=270)
   plot!([-1, 0, -1, -1, 0, 0], [-1, 0, 0, -1, -1, 0], color=:red, linewidth=0.5,
       xlims=(-1.05, 0.05), ylims=(-1.05, 0.05))
   circle(x, r - 1, r, :blue)
end;

 

算額（その12）

算額（その12）

山形県山形市山寺 立石寺（山寺）根本中堂 大正3年
http://www.wasan.jp/yamagata/yamadera.html

下図のように，正方形の中に四分円，大小の円がある。円の半径を求めよ。

これも，ちょっと面倒であるが紙と鉛筆で計算できる。
正方形の一辺の長さを 1，大円，小円の半径と中心座標を r1, (r1, y1), r2, (x2, y2) とする。

using SymPy
@syms r1::positive, r2::positive, y1::positive, x2::positive, y2::positive
eq1 = r1^2 + y1^2 - (1 - r1)^2;
eq2 = (1//2)^2 + (1 - r2)^2 - (1 + r2)^2;
eq3 = (1 - r1)^2 + y1^2 - (1 + r1)^2;

solve([eq1, eq2, eq3], (r1, r2, y1))

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (1/6, 1/16, sqrt(6)/3)

using Plots

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=:center)
   scatter!([x], [y], color=color, markerstrokewidth=0)
   annotate!(x, y, text(string, 10, position, color))
end;

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, y1) = (1//6, 1//16, sqrt(6)/3)
   println("r1 = $r1,  r2 = $r2,  y1 = $y1")
   plot()
   if more
       point(r1, y1, " (r1, y1)", :blue)
       point(1/2, 1 - r2, " (1/2, 1 - r2)", :blue)
   end

   circle(1, 0, 1, beginangle=90, endangle=180)
   plot!([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], color=:red, linewidth=0.5,
       xlims=(-0.05, 1.05), ylims=(-0.05, 1.05))
   circle(r1, y1, r1, :blue)
   circle(1/2, 1 - r2, r2, :blue)
   circle(0, 0, 1, beginangle=0, endangle=90)
end;

算額（その11）

算額（その11）

山形県山形市山寺 立石寺（山寺）根本中堂 大正3年(1914)4月
http://www.wasan.jp/yamagata/yamadera.html

東京都住吉神社 嘉永8年(1851)
http://www.wasan.jp/tokyo/sumiyosi.html

下図のように，半円の中に 3 つの円がある。円の半径を求めよ。

これは簡単。中学生でも暗算で解答できるだろう。

using SymPy
@syms x::positive, y::positive, r::positive
eq1 = x^2 + r^2 - (1-r)^2;
eq2 = x^2 +(1 - 2r)^2 - 4r^2;

solve([eq1, eq2], (x, r))

   1-element Vector{Tuple{Sym, Sym}}:
    (sqrt(3)/3, 1/3)

using Plots

function circle(ox, oy, r, color=:red)
   θ = 0:0.01:2pi
   x = r.*cos.(θ)
   y = r.*sin.(θ)
   plot!(ox .+ x, oy .+ y, color=color)
end;

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

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   x, r = (sqrt(3)/3, 1/3)
   println("$x, $r")
   plot()
   if more
       point(x, r, " (x, r)", :blue)
       point(-x, r, " (-x, r)", :blue)
       point(0, 1 - r, " (0, 1-r)", :blue)
   end

   circle(0, 1 - r, r, :blue)
   circle(x, r, r, :blue)
   circle(-x, r, r, :blue)
   for i = 0:60:180
       plot!([0, cosd(i)], [0, sind(i)], color=:black)
   end
   circle(0, 0, 1, :black)
end;

 

Julia 1.8.2 速い，Python 3.11.0 よりも

Python 3.11.0 になって，Python が速くなって嬉しいって！！

> 自分のPCでも23.2秒から12.5秒になり、2倍ぐらい早くなってるのでこれはうれしい。

フィボナッチ数列を再帰関数により求める場合の話なので，一般的なものではないことをお断りしておく。

私のマシンで比較してみました。

結論を述べておきます。なんだかんだいっても，Python に比べて Julia は 28 倍くらい速い。

Python 3.11.0

Python 3.11.0 (v3.11.0:deaf509e8f, Oct 24 2022, 14:43:23) [Clang 13.0.0 (clang-1300.0.29.30)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> def fb(n):
...     if n == 0 or n == 1:
...         return n
...     else:
...         return fb(n - 2) + fb(n - 1)
... 
>>> from time import time
>>> s = time();fb(40);print(time()-s)
102334155
11.548144340515137

ほどほど早い。（Python 3.10 で確かめようなんてつもりはさらさない）

R はどうか?

R version 4.2.2 RC (2022-10-27 r83209) -- "Innocent and Trusting"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)

> fb = function(n) {
+  if (n < 2) {
+   return(n)
+  } else {
+   return(fb(n - 2) + fb(n - 1))
+  }
+ }
> system.time(print(fb(40)))
[1] 102334155
   user  system elapsed 
 64.022   0.394  64.373 

遅い！！！

でも，gmp ライブラリを使ってずるっこすると，測定限界以下の速度になります

> library(gmp)
> system.time(print(fibnum(40)))
Big Integer ('bigz') :
[1] 102334155
   user  system elapsed 
      0       0       0 

Julia は?

   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.8.2 (2022-09-29)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> function fib(n)
           if n < 2
               return n
           else
               return fib(n - 2) + fib(n - 1)
           end
       end;

julia> @time fib(40)
  0.405616 seconds
102334155

何の文句もない。速い（あたりまえなんだけど）

ちなみに，噛ませ犬で AWK

function fb(n) {
    if (n < 2) {
        return(n)
    }
    else  {
        return(fb(n - 2) + fb(n - 1))
    }
}

BEGIN {
    print fb(40)
}

というファイルを用意して

$ time awk -f fb.awk
102334155
37.665u 0.136s 0:37.81 99.9% 0+0k 0+0io 0pf+0w

意外にも R より速い

算額（その10）

算額（その10）

山形県米沢市小野川町 小野川薬師堂 明治10年
http://www.wasan.jp/yamagata/onogawa.html

円内に半径の異なる 3 種類の円が下図のように収まっている。大円と小円は直線に接している。 外円の直径が 9 寸のとき，大円の直径はいかほどか。

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

include("julia-source.txt");

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

    (5*R/9, 20*R/49, -3*R/7, 20*R/81, 20*R/27)

大円の半径 r1 は，外円の半径 R の 5/9 倍である。
外円の直径が 9 寸のとき，大円の直径は 5 寸である。

その他のパラメータは以下のとおりである。

    R = 4.5; r1 = 2.5; r2 = 1.83673; y2 = -1.92857; r3 = 1.11111; x3 = 3.33333

function draw(R, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, y2, r3, x3) = (5*R/9, 20*R/49, -3*R/7, 20*R/81, 20*R/27)
    @printf("R = %g;  r1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  x3 = %g\n", R, r1, r2, y2, r3, x3)
    y = R - 2r1
    x = sqrt(R^2 - y^2)
    plot()
    circle(0, 0, R)
    circle(0, R - r1, r1, :green)
    circle2(r2, y2, r2, :magenta)
    circle2(x3, y + r3, r3, :blue)
    segment(-x, y, x, y, :black)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        vline!([0], color=:gray80, lw=0.5)
        point(0, R - r1, "大円:R,(0,R-r1)", :green, :center, delta=-delta/2)
        point(r2, y2, "中円:r2,(r2,y2)", :magenta, :center, delta=-delta/2)
        point(x3, y + r3, "小円:r3\n(x3,R-2r1+r3)", :blue, :center, :bottom, delta=delta/2)
    end
end;

draw(9/2, true)

算額（その9）

算額（その9）

大きな円の中に 4 種類の大きさの異なる円 A, B, C, D があり，図のように互いに接している（すべての円が相互に接しているわけではない）。

外側の円の半径を 1 としたとき，それぞれの円の半径を求めよ。

山形県鶴岡市大山 椙尾神社 文政元年8月
http://www.wasan.jp/yamagata/sugio.html

円 A, B, C, D の半径を r1, r2, r3, r4 とする。円 C, D の中心座標を (x3, y3), (x4, y4) とし，方程式を立てる。

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

cos30d = sqrt(Sym(3))/2;
sin30d = 1//2;
x3 = (r1 + r3)cos30d;
y3 = (r1 + r3)sin30d;
x4 = (1 - r4)cos30d;
y4 = -(1 - r4)sin30d;

eq1 = (x4 + r1*cos30d)^2 + (r1*sin30d-y4)^2 - r4^2;
eq2 = x4^2 + (1 - r2 -  y4)^2 - (r2 + r4)^2;
eq3 = (x3 - x4)^2 + (y3 - y4)^2 - (r4 - r3)^2;

この図形においては，r1, r2, r3, r4 は独立に決まるわけではなく，r1 を決めると残りの r2, r3, r4 が決まる。

eq1, eq2, eq3 を r2, r3, r4 について解くと，それぞれは r1 についての式になる。

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

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (-3*(r1 - 1)/(r1 + 7), -3*r1*(r1 - 1)/(7*r1 + 1), (r1 + 1)/2)

なお，算額では A と B の大きさが同じように見受けられるので，以下（r1 = r2）を解けば対応する図を得る。

solve(-3*(r1 - 1)/(r1 + 7) - r1)[1] |> println

   -5 + 2*sqrt(7)

using Plots

function circle(ox, oy, r, color=:red)
   θ = 0:0.01:2pi
   x = r.*cos.(θ)
   y = r.*sin.(θ)
   plot!(ox .+ x, oy .+ y, color=color)
end;

function point(x, y, string="", color=:green, position=:left, vertical=:center)
   # scatter!([x], [y])
   annotate!(x, y, text(string, 10, position, color))
end;

function draw(r1=2*sqrt(7) - 5; more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (-3*(r1 - 1)/(r1 + 7), -3*r1*(r1 - 1)/(7*r1 + 1), (r1 + 1)/2)
   cos30d = sqrt(Sym(3))/2;
   sin30d = 1/2;
   x3 = (r1 + r3)cos30d;
   y3 = (r1 + r3)sin30d;
   x4 = (1 - r4)*cos30d;
   y4 = -(1 - r4)*sin30d;
   println("$r1, $r2, $r3, $r4")
   plot()
   if more
       point(0, 0, "A", :green)
       point(0, 1 - r2, "B", :magenta)
       point(x3, y3, "C", :brown)
       point(0.8, 0, "D", :blue)
   end

   circle(0, 0, 1)

   circle(0, 0, r1, :green)

   circle(0, 1 - r2, r2, :magenta)
   circle((1 - r2)cos30d, -(1 - r2)sin30d, r2, :magenta)
   circle((-(1 - r2))cos30d, -(1 - r2)sin30d, r2, :magenta)

   circle(x3, y3, r3, :brown)
   circle(-x3, y3, r3, :brown)
   circle(0, -r1-r3, r3, :brown)

   circle(x4, y4, r4, :blue)
   circle(-x4, y4, r4, :blue)
   circle(0, 1-r4, r4, :blue)
end;

r1 = 0.26 のとき

   0.26, 0.30578512396694213, 0.20468085106382977, 0.63

r1 = r2 = 2*sqrt(7) - 5 のとき

   0.29150262212918143, 0.29150262212918104, 0.20377661238703057, 0.6457513110645907

 

算額（その8）

算額（その8）

円の内部に大きさの異なる 5 種類の円がある。それぞれは互いに接している。
円の半径を求めよ。

山形県新庄市神明町 総鎮守天満宮 文政元年7月（現在は戸沢神社にあり）
http://www.wasan.jp/yamagata/tozawa.html

なお，算額には上半分の図が描かれているが，ここでは下半分も描くことにする。

外側の円の半径を 2 とする。内側の一番大きい円の半径は r1 で，以下 r2, r3, r4, r5 とする。

それぞれの円の中心座標を (0, 1), (x2, r2), (x3, y3), (x4, y4), (x5, y5) とする。

互いに接することから以下の 11 本の方程式を立てる。

using SymPy
@syms r2::positive, r3::positive, r4::positive, r5::positive
@syms x2::positive, x3::positive, x4::positive, x5::positive
@syms y3::positive, y4::positive, y5::positive
eq1 = x3^2 + (y3 - 1)^2 - (1 + r3)^2;
eq2 = x5^2 + (y5 - 1)^2 - (1 + r5)^2;
eq3 = x2^2 + (1 - r2)^2 - (1 + r2)^2;
eq4 = x2^2 + r2^2 - (2 - r2)^2;
eq5 = x4^2 + y4^2 - (2 - r4)^2;
eq6 = x3^2 + y3^2 - (2 - r3)^2;
eq7 = (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2;
eq8 = (x4 - x3)^2 + (y3 - y4)^2 - (r3 + r4)^2;
eq9 = (x2 - x4)^2 + (y4 - r2)^2 - (r2 + r4)^2;
eq10 = (x2 - x5)^2 + (y5 - r2)^2 - (r2 + r5)^2;
eq11 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2;

未知数が 11 個で，方程式が 11 本なので，理論的には解けるが，SymPy の solve では時間ばかりかかって一向に解が求まらない。

eq9 と eq11 を図と見比べると，x2 = x4, x3 = x5 であることがわかる。未知数の個数は変わらないが式の数を増やす。

eq12 = x2 - x4;
eq13 = x3 - x5;

solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13], (r2, r3, r4, r5, x2, x3, x4, x5, y3, y4, y5))

解が求まった。

   1-element Vector{NTuple{11, Sym}}:
    (1/2, 1/5, 1/6, 2/15, sqrt(2), 4*sqrt(2)/5, sqrt(2), 4*sqrt(2)/5, 7/5, 7/6, 16/15)

using Plots

function circle(ox, oy, r, color=:red)
   θ = 0:0.01:2pi
   x = r.*cos.(θ)
   y = r.*sin.(θ)
   plot!(ox .+ x, oy .+ y, color=color)
end;

function point(x, y, string="", position=:left, vertical=:center)
   scatter!([x], [y])
   annotate!(x, y, text(string, 10, position))
end;

function circle2(x, y, r, color)
   circle( x,  y, r, color)
   circle( x, -y, r, color)
   circle(-x,  y, r, color)
   circle(-x, -y, r, color)
end

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4, r5, x2, x3, x4, x5, y3, y4, y5) = (1/2, 1/5, 1/6, 2/15, sqrt(2), 4*sqrt(2)/5, sqrt(2), 4*sqrt(2)/5, 7/5, 7/6, 16/15)
   plot()
   circle(0, 0, 2)
   circle(0, 1, 1, :blue)
   circle(0, -1, 1, :blue)
   circle2(x2, r2, r2, :green)
   circle2(x3, y3, r3, :magenta)
   circle2(x4, y4, r4, :brown)
   circle2(x5, y5, r5, :black)
end;

draw()
savefig("fig.png")

なお，数値的に解を求めることもできる。

この場合は，方程式は 11 本でよい。

# import Pkg; Pkg.add("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])
       v = r.zero[1]
   else
       r = nlsolve((vout,vin)->vout .= func(vin,params...), ini)
       v = r.zero
   end
   return v, r.f_converged
end;

function h(u)
   r2, r3, r4, r5, x2, x3, x4, x5, y3, y4, y5 = u
   return [
       x3^2 + (y3 - 1)^2 - (1 + r3)^2,
       x5^2 + (y5 - 1)^2 - (1 + r5)^2,
       x2^2 + (1 - r2)^2 - (1 + r2)^2,
       x2^2 + r2^2 - (2 - r2)^2,
       x4^2 + y4^2 - (2 - r4)^2,
       x3^2 + y3^2 - (2 - r3)^2,
       (x4 - x5)^2 + (y4 - y5)^2 - (r4 + r5)^2,
       (x4 - x3)^2 + (y3 - y4)^2 - (r3 + r4)^2,
       (x2 - x4)^2 + (y4 - r2)^2 - (r2 + r4)^2,
       (x2 - x5)^2 + (y5 - r2)^2 - (r2 + r5)^2,
       (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^2
   ]
end;

iniv = [0.525,  0.225,  0.175,  0.125,  1.45,  1.175,  1.45,  1.025,  1.45,  1.2,  1.1]  # 初期値

nls(h, ini=iniv)

   ([0.4999999999999999, 0.19999999995440246, 0.16666666671681116, 0.13333333457600596, 1.4142135623730967, 1.1313708498182247, 1.414213562263483, 1.1313708505474644, 1.4000000001367925, 1.166666666725507, 1.0666666672577463], true)

 

算額（その7）

算額（その7）

円の内部に，一番大きい円から一番小さい円まで4種類の大きさの円が互いに接している。4つの円の半径を求めよ。

山形県新庄市神明町 総鎮守天満宮 文政元年7月（現在は戸沢神社にあり）
http://www.wasan.jp/yamagata/tozawa.html

これは，山形県鶴尾市大山 椙尾神社 文政元年8月に奉納されたものより一段階複雑なものだ。
http://www.wasan.jp/yamagata/sugio.html
算額の問題解決編
https://blog.goo.ne.jp/r-de-r/e/903d534eb0f73787d2a98800a44cfa95

外側の円の半径を 1 として，内部の円の半径を大きい順に r1, r2, r3, r4 とする。図を描くために必要な 2 つの円の中心の x 座標も求める。

円が互いに接しているので順次方程式を立て，連立方程式を解く。

注：当初，eq4 = x2^2 + r2^2 - (1 - r2)^2 を eq4 = sqrt(x2^2 + r2^2) + r2 - 1 としていたので解が求まらなかった。

using SymPy
@syms r1::positive, r2::positive, r3::positive, r4::positive, x2::positive, x3::positive;

eq1 = x2^2 + (1- r1 - r2)^2 - (r1 + r2)^2;  # 円1と2円が接することから
eq2 = (x2 - x3)^2 + r2^2 - (r2 +r3)^2;  # 円2と円3が接することから
eq3 = x3^2 + r4^2 - (r3 + r4)^2;  # 円3と円4が接することから
eq4 = x2^2 + r2^2 - (1 - r2)^2;   # 円2と外の円が接することから
eq5 = 2r1 + 2r4 - 1;  # 円1と円4が接することから
eq6 = (1 - r1)^2 + x3^2 - (r1 + r3)^2;  # 円1と円3が接することから

solve([eq1, eq2, eq3, eq4, eq5, eq6], (r1, r2, r3, r4, x2, x3))  # 方程式を解く

解は以下の通り。

   1-element Vector{NTuple{6, Sym}}:
    (3/7, 2/7, 1/5, 1/14, sqrt(21)/7, 2*sqrt(21)/35)

using Plots

function circle(ox, oy, r, color=:red)
   θ = 0:0.01:2pi
   x = r.*cos.(θ)
   y = r.*sin.(θ)
   plot!(ox .+ x, oy .+ y, color=color)
end;

function point(x, y, string="", position=:left, vertical=:center)
   scatter!([x], [y])
   annotate!(x, y, text(string, 10, position))
end;

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3, r4, x2, x3) = (3/7, 2/7, 1/5, 1/14, sqrt(21)/7, 2*sqrt(21)/35)
   plot()
   circle(0, 0, 1)  # 外側の円
   circle(0, 1 - r1, r1, :blue)  # 一番大きい円
   circle(0, r1 - 1, r1, :blue)
   circle( x2,  r2, r2, :green)  # 二番目に大きい円
   circle( x2, -r2, r2, :green)
   circle(-x2,  r2, r2, :green)
   circle(-x2, -r2, r2, :green)
   circle( x3, 0, r3, :brown)  # 三番目に大きい円
   circle(-x3, 0, r3, :brown)
   circle(0,  r4, r4, :magenta)  # 一番小さい円
   circle(0, -r4, r4, :magenta)
   if more
       point(0, 1-r1, " (0, 1-r1)")
       point(x2, r2, " (x2, r2)")
       point(x3, 0, " (x3, 0)", :left, :top)
       point(0, r4, "   (0, r4)", :left)
   end
end;

draw(true)
savefig("fig1.png")

draw()
savefig("fig2.png")

 

算額（その6）

算額（その6）

香川県高松市　石清尾八幡宮 文化三年
http://www.wasan.jp/kagawa/iwaseo.html

千葉県富津市 吾妻神社 明治十年
http://www.wasan.jp/chiba/azuma.html

• 外側の円の半径を 1 とする。
• 内部の三角形は正三角形である。
• 円 A の半径を 0.17 とする。
• それぞれの円と正三角形は互いに接している。
• 円 B の中心 (x1, y1)，半径 r1 を求める。

using SymPy
@syms x, y, r::positive, r1::positive, r2::positive, x1::positive, y1::positive, x2::positive, y2::positive;

eq1 = (1 - r2 - y1)^2 + x1^2 - (r1 + r2)^2;  # 円 A，円 B が接する
eq2 = y2 - y1 - (x2 - x1)/√Sym(3);  # 円 B の中心と接点 C の直線の方程式
eq3 = y2 - (1 - 2r2) + x2√Sym(3);   # 接点 C を通る正三角形の辺の方程式
res = solve([eq2, eq3], (x2, y2));  # 接点 x2, y2 の座標
eq4 = (res[y2] - y1)^2 + (res[x2] - x1)^2 - r1^2;  #  BC が半径 r1
eq5 = sqrt(x1^2 + y1^2) + r1 - 1;  # 外側の円と円 B が接する
a = 0.17;
EQ1 = eq1(r2 => a);  # 方程式は繊細で，r2 = 0.17 でないと解けない？
EQ4 = eq4(r2 => a);
res2 = solve([EQ1, EQ4, eq5], (x1, y1, r1))  # x1, y1, r1 を求める

   1-element Vector{Tuple{Sym, Sym, Sym}}:
    (0.396484958572158, 0.575520654261460, 0.301126373472639)

using Plots

function circle(ox, oy, r; color=:red)
   θ = 0:0.01:2pi
   x = r.*cos.(θ)
   y = r.*sin.(θ)
   plot!(ox .+ x, oy .+ y, color=color)
end;

function triangle(ox, oy, r; color=:blue)
   θ = [90, 210, 330, 90]
   x = r .* cosd.(θ)
   y = r .* sind.(θ)
   plot!(ox .+ x, oy .+ y, color=color)
end;

function point(x, y, string="", position=:left)
   scatter!([x], [y])
   annotate!(x, y, text(string, 10, position))
end;

function draw(more=false)
   pyplot(size=(500, 500), aspectratio=1, label="", fontfamily="IPAMincho")
   r, r2 = (1, 0.17)
   x1, y1, r1 = (0.396484958572158, 0.575520654261460, 0.301126373472639)
   plot()
   circle(0, 0, r)
   triangle(0, 0, r　-　2r2)
   circle(0, r　-　r2, r2, color=:black)
   circle((r　-　r2)*cosd(330), (r　-　r2)*sind(330), r2, color=:black)
   circle((r　-　r2)*cosd(210), (r　-　r2)*sind(210), r2, color=:black)
   circle(x1, y1, r1)
   circle(-x1, y1, r1)
   circle(x1+2*r1*sind(30), y1-2*r1*cosd(30), r1)
   circle(-(x1+2*r1*sind(30)), y1-2*r1*cosd(30), r1)
   circle(-r1, -r1 - (r - 2r2)*sind(30), r1)
   circle(r1, -r1 - (r - 2r2)*sind(30), r1)
   if more
       point(0, 0, " O")
       point(0, 1 - r2, " A")
       point(x1, y1, "B ", :right)
       point(x1, y1, " (x1, y1)")
       x2 = x1/4 + sqrt(3)*(1/4 - r2/2 - y1/4)
       y2 = (3*y1 + 1 - sqrt(3)*x1)/4 - r2/2
       point(x2, y2, "C ", :right)
       point(x2, y2, " (x2, y2)")
       plot!([x1, x2], [y1, y2])
       annotate!((x1 + x2)/2, (y1 + y2)/1.8, text("r1", 10, :right))
   end
end;