和算の心（その009）

和算の心（その009）

円の直径，弦，矢の関係式を「径矢弦の定理」と呼ぶことがある。

直角三角形 OBD において，ピタゴラスの定理（三平方の定理）を適用する。
OB² + BD² = OD²
(径/2 - 矢)² + (弦/2)² = (径/2)²

方程式:(径/2 - 矢)^2 + (弦/2)^2 - (径/2)^2 は，径，矢，弦のうちの 2 つが与えられると，残りの 1 つが計算できる。

include("julia-source.txt");
@syms 径::posiitive, 矢::positive, 弦::positiiive
eq = (径/2 - 矢)^2 + (弦/2)^2 - (径/2)^2

径，矢がわかっているときに，弦を求める。
2 通りの解が得られるが，正のものが適解である。
弦 = 2*sqrt(矢)*sqrt(径 - 矢)

ans_弦 = solve(eq, 弦)
ans_弦 |> println

   Sym{PyCall.PyObject}[-2*sqrt(矢)*sqrt(径 - 矢), 2*sqrt(矢)*sqrt(径 - 矢)]

径，弦が分かっているときに，矢を求める。
2 通りの解が得られる。小さい方が「矢」，大きい方は「大矢」と呼ばれることもある。「矢 + 大矢 = 径」である。
矢 = 径/2 - sqrt(-弦^2 + 径^2)/2
大矢 = 径/2 + sqrt(-弦^2 + 径^2)/2

ans_矢 = solve(eq, 矢)
ans_矢 |> println

   Sym{PyCall.PyObject}[径/2 - sqrt(-弦^2 + 径^2)/2, 径/2 + sqrt(-弦^2 + 径^2)/2]

矢，弦が分かっているときに，径を求める。
径 = 弦^2/(4*矢) + 矢

ans_径 = solve(eq, 径)
ans_径 |> println

   Sym{PyCall.PyObject}[弦^2/(4*矢) + 矢]

function draw(more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   弦 = 0.8
   径 = 1
   矢 = 径/2 - sqrt(-弦^2 + 径^2)/2
   矢2 = 径/2 + sqrt(-弦^2 + 径^2)/2
   @printf("弦 = %g, 径 = %g のとき，矢 = %g, 大矢 = %g\n", 弦, 径, 矢, 矢2)
   径 = 1
   矢 = 0.2
   弦 = 2*sqrt(矢*(径 - 矢))
   @printf("径 = %g, 矢 = %g のとき，弦 = %g\n", 径, 矢, 弦)
   矢 = 0.2
   弦 = 0.8
   径 = 弦^2/(4*矢) + 矢
   @printf("矢 = %g, 弦 = %g のとき，径 = %g\n", 矢, 弦, 径)
   plot()
   circle(0, 0, 径/2)
   segment(-弦/2, 径/2 - 矢, 弦/2, 径/2 - 矢)
   segment(-径/2, 0, 径/2, 0, :blue)
   segment(0, 径/2, 0, 径/2 - 矢, :green)
   segment(0, 0, 0, 径/2 - 矢, :orange)
   segment(0, 0, 弦/2, 径/2 - 矢, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       #hline!([0], color=:gray80, lw=0.5)
       #vline!([0], color=:gray80, lw=0.5)
       point(径/4, 径/4 - 矢/2, "OD:径/2", :magenta, :center, delta=-delta, mark=false)
       point(-径/4, 0, "EF:径（直径）", :blue, :center, delta=-delta, mark=false)
       point(-弦/4, 径/2 - 矢, " CD=弦", :black, :center, delta=-delta, mark=false)
       point(弦/4, 径/2 - 矢, " BD=弦/2", :black, :center, delta=-delta, mark=false)
       point(0, 径/2 - 矢/2, " AB:矢", :green, :left, delta=-delta/2, mark=false)
       point(0, 径/2 - 3矢/2, "OB:径/2-矢 ", :orange, :right, delta=-delta/2, mark=false)
       point(0, 径/2, "A", :black, :center, :bottom, delta=delta/2)
       point(0, 径/2 - 矢, "B ", :black, :right, :bottom, delta=delta/2)
       point(-弦/2, 径/2 - 矢, "C", :black, :right, :bottom, delta=delta/2)
       point(弦/2, 径/2 - 矢, "D", :black, :left, :bottom, delta=delta/2)
       point(-径/2, 0, "E ", :black, :right, :bottom, delta=delta/2)
       point(径/2, 0, " F", :black, :left, :bottom, delta=delta/2)
       point(0, 0, "O", :black, :center, delta=-delta)
       xlims!(-径/2 - 5delta, 径/2 + 5delta) 
   end
end;

draw(true)

 

さぬきうどん 山よし 佐文店

まんのう町佐文 さぬきうどん 山よし 佐文店

多度津町西港町に本店があるのかな？
懸垂幕に注目

算額（その1394）

算額（その1394）

十五 武州金鑚村 金鑚寺 文化11年(1814)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：四辺形，面積最大
#Julia, #SymPy, #算額, #和算

甲斜，乙斜，丙斜，丁斜の 4 辺からなる不等辺四角形がある。甲斜，乙斜，丙斜が与えられたとき，面積が最大となるときの丁斜を求める術を述べよ。

甲斜，乙斜，丙斜，丁斜の長さを「甲」，「乙」，「丙」，「丁」とする，また，対角線の長さを「己」，「庚」とする。

面積が最大となるのは，丁斜が四辺形が内接する円の直径の場合である。

トレミーの定理により，甲*丙 + 乙*丁 = 戊*己 である。

また，戊^2 = 丁^2 - 甲^2，己^2 = 丁^2 - 丙^2 なので，

(甲*丙 + 乙*丁)^2 = (丁^2 - 甲^2)*(丁^2 - 丙^2)
丁^4 - 丁^2*丙^2 - 丁^2*乙^2 - 丁^2*甲^2 - 2*丁*丙*乙*甲 = 0

よっt，甲乙丙丁が満たすべき式は以下のもの。

丁^3 - 丁*(丙^2 + 乙^2 + 甲^2) - 2*丙*乙*甲 = 0

using SymPy
@syms 甲, 乙, 丙, 丁, 戊, 己
eq = 丁^3 - 丁*(丙^2 + 乙^2 + 甲^2) - 2*丙*乙*甲;

eq(甲 => 12, 乙 => 9, 丙 => 2, 丁 => 16)

   0

eq(甲 => 2, 乙 => 7, 丙 => 11, 丁 => 14)

   0

 

算額（その1393）

算額（その1393）

四 武州川越町 八幡宮 寛政7年(1795)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円6個，共通接線，等脚台形
#Julia, #SymPy, #算額, #和算

共通接線2本を持つ 3 個の円が 2 セットある。甲円，乙円の直径が 81 寸，16 寸のとき，戊円の直径はいかほどか。

甲円，丙円，戊円と，乙円，丁円，戊円の 3 個の円が，共通接線 2 本で挟まれている。台形は必要な条件ではない。和算の心（その008）を 2 回適用し，連立方程式を解けばよい。

甲円，乙円，丙円，丁円，戊円の半径を r1, r2, r3, r4, r5 とする。
甲円，丙円，戊円のセットで，円の相似比を p
乙円，丁円，戊円のセットで，円の相似比を q
とする。
共通接線の交点から甲円との接点までの長さは 2r1*√p/(1 - p)
もう一組の共通接線の交点から乙円との接点までの長さは 2r2*√q/(1 - q)
一本の接線は共通なので，2 つの距離は等しい。
以上の関係式を eq0 としてまとめる。

using SymPy
@syms r1, r2, r5, p, q
p = sqrt(r5/r1)
q = sqrt(r5/r2)
eq0 = 2r1*√p/(1 - p) - 2r2*√q/(1 - q);

r1, r2 に実値を代入し，方程式を解き r5 を求める。

eq1 = eq0(r1 => 81//2, r2 => 16//2) |> simplify |> numerator
eq1 |> println

   2^(1/4)*r5^(1/4)*(-211*sqrt(2)*sqrt(r5) + 684)

直径は r5 の 2 倍で，分数になる。

2solve(eq1, r5)[2] |> println  # 2 of 2

   467856/44521

帯分数は小学校時代以来使わないが，467856/44521 = 10 + 22646/44521 である。

divrem(467856,44521)

   (10, 22646)

 

和算の心（その008）

和算の心（その008）

#Julia, #SymPy, #算額, #和算

直線 2 本に挟まれ，互いに外接する 3 個の円の位置関係

図のように，直線（接線） 2 本に甲円，乙円，丙円が接している。
甲円と丙円の直径（2r1, 2r3）が分かっているとき，2 直線の交点 O と甲円と直線の接点 A, B を結ぶ線分の長さ OA = OB を求める術を述べよ。

甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
とおく。

eq1, eq2 は，3 個の円が相似であることを表している。
eq3, eq4 は，互いに接する直線上の 2 円の間の水平距離に関する式である。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive, x2::positive, r3::positive, x3::positive
eq1 = r3/x3 - r2/x2
eq2 = r2/x2 - r1/x1
eq3 = (x2 - x3) - 2sqrt(r2*r3)
eq4 = (x1 - x2) - 2sqrt(r1*r2);

連立方程式を解き，r2, x1, x2, x3 を一度に求めることもできるが，逐次的に x1 だけを求めてみよう。

相似比を p とすれば，円の半径には以下の関係が成り立つ。
r2/r1 = r3/r2 = p
r2 = p*r1, r3 = p*r2 より r3 = p^2 * r1
甲円と丙円の半径がわかれば，相似比は「丙円と甲円の比の平方根」である。
p = √(r3/r1) 

直線との接点の座標にも
x2/x1 = r2/r1 = p
(x1 - 2√(r1*r2))/x1 = p
の関係があるので，甲円と直線の接点を結ぶ線分の長さは
x1 = 2r1*√p/(1 - p)
である。

甲円と丙円の直径が 20 寸，5寸のとき，2 直線の交点と甲円と直線の接点を結ぶ線分の長さは，56.5685424949238 寸である。

p = sqrt(r3/r1)
ans_x1 = 2r1*√p/(1 - p) |> simplify
ans_x1 |> println

   2*r1^(5/4)*r3^(1/4)/(sqrt(r1) - sqrt(r3))

ans_x1(r1 => 20, r3 => 5).evalf() |>  println

   56.5685424949238

---

連立方程式を解くと，x1 の式はかなり長く，SymPy では直接は簡約化することができないが，以下のように式変換すると上述の式と同じものが得られる。

res = solve([eq1, eq2, eq3, eq4], (r2, x1, x2, x3))[4];  # 4 of 4

res[1](r1 => 20, r3 => 5).evalf() |> println
res[2](r1 => 20, r3 => 5).evalf() |> println
res[3](r1 => 20, r3 => 5).evalf() |> println
res[4](r1 => 20, r3 => 5).evalf() |> println

   10.0000000000000
   56.5685424949238
   28.2842712474619
   14.1421356237310

res[2] |> println

   2*r1*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/r3

ルートの中の式を再定義する（分数式にするため / を // に変換）。

t = (r1^(7//2)*r3^(5//2) - r1^(5//2)*r3^(7//2) - r1^(3//2)*r3^(9//2) + sqrt(r1)*r3^(11//2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)
t |> println

   (r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)

√r1 を u, √r3 を v に置換し，因数分解する。

@syms u, v
t2 = t(√r1 => u, √r3 => v) |> factor
t2 |> println

   u*v^5/(u - v)^2

元の式に戻す。

t3 = (2r1*sqrt(t2)/r3)(u => √r1, v => √r3)
t3 |> println

   2*r1^(5/4)*r3^(1/4)/Abs(sqrt(r1) - sqrt(r3))

t3(r1 => 20, r3 => 5).evalf()　|> println

   56.5685424949238

using LaTeXStrings
function draw(r1, r3, more=false)
    #pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    gr(size=(700, 500), grid=false, aspectratio=1, label="", fontfamily="Helvetica")
   (r2, x1, x2, x3) = (r1*(2*r1^(7/2)*r3^(5/2) - 2*r1^(5/2)*r3^(7/2) - 2*r1^(3/2)*r3^(9/2) + 2*sqrt(r1)*r3^(11/2) + 4*r1^3*r3^3 - 8*r1^2*r3^4 + 4*r1*r3^5 + r3*(r3 - sqrt((4*r1^(7/2)*r3^(5/2) - 4*r1^(5/2)*r3^(7/2) - 4*r1^(3/2)*r3^(9/2) + 4*sqrt(r1)*r3^(11/2) + 8*r1^3*r3^3 - 16*r1^2*r3^4 + 8*r1*r3^5 + r3^2*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)))*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(2*(r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)), 2*r1*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/r3, r1*(2*r1^(7/2)*r3^(5/2) - 2*r1^(5/2)*r3^(7/2) - 2*r1^(3/2)*r3^(9/2) + 2*sqrt(r1)*r3^(11/2) + 4*r1^3*r3^3 - 8*r1^2*r3^4 + 4*r1*r3^5 + r3*(r3 - sqrt((4*r1^(7/2)*r3^(5/2) - 4*r1^(5/2)*r3^(7/2) - 4*r1^(3/2)*r3^(9/2) + 4*sqrt(r1)*r3^(11/2) + 8*r1^3*r3^3 - 16*r1^2*r3^4 + 8*r1*r3^5 + r3^2*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)))*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))/(r3*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4))*(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)), 2*sqrt((r1^(7/2)*r3^(5/2) - r1^(5/2)*r3^(7/2) - r1^(3/2)*r3^(9/2) + sqrt(r1)*r3^(11/2) + 2*r1^3*r3^3 - 4*r1^2*r3^4 + 2*r1*r3^5)/(r1^4 - 4*r1^3*r3 + 6*r1^2*r3^2 - 4*r1*r3^3 + r3^4)))
   p = r2/r1
   println(2r1*√p/(1 - p))
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g\n", r1, x1, r2, x2, r3, x3)
   p1 = plot()  # showaxis=false)
   circle(x1, r1, r1)
   circle(x2, r2, r2, :blue)
   circle(x3, r3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x1, r1, L"r_1,(x_1,r_1)", :red, :center, delta=-delta)
       point(x2, r2, L"r_2,(x_2,r_2)", :blue, :center, delta=-delta)
       point(x3, r3, L"r_3,(x_3,r_3)", :green, :center, delta=-delta)
       θ = 2atand(r1/x1)
       abline(0, 0, tand(θ), 0, x1)
       segment(0, 0, 1.2x1, 0)
       point(x1 - r1*sind(θ), r1 + r1*cosd(θ), "B", :black, :right, :bottom, delta=delta)
       point(x2 - r2*sind(θ), r2 + r2*cosd(θ))
       point(x3 - r3*sind(θ), r3 + r3*cosd(θ))
       point(x1, 0, "A", :black, :center, delta=-delta)
       point(x2, 0)
       point(x3, 0)
       point(0, 0, "O", :black, :center, delta=-delta)
       ylims!(-7delta, 2r1 + delta)
   end
end;

draw(20, 5, true)

ぶっかけうどん はな庄

高松市香川町川東上 ぶっかけうどん はな庄
ぶっかけうどんがおすすめのようなので，温玉ぶっかけ小をお願いした😄

算額（その1392）

算額（その1392）

十九 羽生市中手子林 八幡神社 文政元年(1818)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，三角形，正六角形
#Julia, #SymPy, #算額, #和算

（不等辺）三角形の中に，天円，地円，人円と正六角形を容れる。天円，地円の直径が 4 寸，6 寸のとき，人円の直径はいかほどか。

三角形の頂点座標を (3a, 0), (x0, y0), (b, 0)
正六角形の一辺の長さと中心座標を 2a, (0, √3a)
天円の半径と中心座標を r1, (x1, 2√3a + r1)
地円の半径と中心座標を r2, (2a, r2)
人円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms a::positive, b::negative,
     r1::positive, x1::positive,
     r2::positive,
     r3::positive, x3::negative
# a = √Sym(3)r2

eq1 = dist(-a, 0, -2a, √Sym(3)a, x3, r3) - r3^2
eq2 = dist(b, 0, -a, 2√Sym(3)a, x3, r3) - r3^2
eq3 = dist(b, 0, -a, 2√Sym(3)a, x1, 2√Sym(3)a + r1) - r1^2
eq4 = dist(3a, 0, 2a, √Sym(3)a, x1, 2√Sym(3)a + r1) - r1^2;

一度には解けないので，eq1, eq3, eq4 を連立方程式として，b, x1, x3 を求める。
なお，a = √3r2 であるが，a はそのままにしておく。

res = solve([eq1, eq3, eq4], (b, x1, x3))[3]  # 3 of 4

   (a*(24*sqrt(3)*a^3 - 48*a^2*r1 - 4*sqrt(3)*a*r1^2 + 3*r1^3)/(r1*(-12*a^2 + 4*sqrt(3)*a*r1 + 3*r1^2)), a - sqrt(3)*r1, -a - sqrt(3)*r3)

b, x1, x3 を eq2 に代入し(eq12 とする），r3 を求める。

eq12 = eq2(b => res[1], x1 => res[2], x3 => res[3]) |> simplify;
res2 = solve(eq12, r3)[1];  # 1 of 2

ここで，a = √Sym(3)r2 を代入し，因数分解する。
結果は，長精度整数を含む長い式になる。

ans_r3 = res2(a => √Sym(3)r2) |> factor;
ans_r3 |> println

   -(r1^16 + 20*r1^15*r2 + 70*r1^14*r2^2 - 876*r1^13*r2^3 - 5148*r1^12*r2^4 + 16974*r1^11*r2^5 + 110601*r1^10*r2^6 - 293274*r1^9*r2^7 - 1113426*r1^8*r2^8 + 4030560*r1^7*r2^9 + 2711232*r1^6*r2^10 - 27387072*r1^5*r2^11 + 34572096*r1^4*r2^12 + 22394880*r1^3*r2^13 - 95738112*r1^2*r2^14 + 90699264*r1*r2^15 - 30233088*r2^16 + sqrt(r1^32 + 34*r1^31*r2 + 321*r1^30*r2^2 - 1406*r1^29*r2^3 - 37418*r1^28*r2^4 - 54378*r1^27*r2^5 + 1866549*r1^26*r2^6 + 5556726*r1^25*r2^7 - 62431686*r1^24*r2^8 - 193177494*r1^23*r2^9 + 1667661669*r1^22*r2^10 + 3566941110*r1^21*r2^11 - 36334402551*r1^20*r2^12 - 22380759756*r1^19*r2^13 + 596476955268*r1^18*r2^14 - 616148677152*r1^17*r2^15 - 6288496227900*r1^16*r2^16 + 18622534127616*r1^15*r2^17 + 24103256904576*r1^14*r2^18 - 214678876124160*r1^13*r2^19 + 294385908328704*r1^12*r2^20 + 800063891564544*r1^11*r2^21 - 3753687721439232*r1^10*r2^22 + 5063157982642176*r1^9*r2^23 + 4824832605788160*r1^8*r2^24 - 33676529335271424*r1^7*r2^25 + 73507309933658112*r1^6*r2^26 - 99895501823016960*r1^5*r2^27 + 93672133367169024*r1^4*r2^28 - 61376067146612736*r1^3*r2^29 + 27116508430467072*r1^2*r2^30 - 7312316880125952*r1*r2^31 + 914039610015744*r2^32))/(r1*(r1 - 2*r2)^3*(r1 + 6*r2)^6*(2*r1 - 3*r2)*(r1^2 - 3*r1*r2 + 3*r2^2)^2)

式を簡約化しなくても，r1, r2 に数値を代入すれば，数値解が得られる。

ans_r3(r1 => 4//2, r2 => 6//2) |> println

   11/2

天円，地円の直径が 6 寸，4 寸のとき，人円の直径は 11 寸になる。

簡約化

SymPy では結果の式を自動では簡約化できないので，手作業を交えて簡約化を試みる。
r3 = (項1 - sqrt(項2)) / 項3 の形をしている。

まず，分子の (項1 - sqrt(項2)) を取り出す

num = numerator(ans_r3);

   -r1^16 - 20*r1^15*r2 - 70*r1^14*r2^2 + 876*r1^13*r2^3 + 5148*r1^12*r2^4 - 16974*r1^11*r2^5 - 110601*r1^10*r2^6 + 293274*r1^9*r2^7 + 1113426*r1^8*r2^8 - 4030560*r1^7*r2^9 - 2711232*r1^6*r2^10 + 27387072*r1^5*r2^11 - 34572096*r1^4*r2^12 - 22394880*r1^3*r2^13 + 95738112*r1^2*r2^14 - 90699264*r1*r2^15 + 30233088*r2^16 - sqrt(r1^32 + 34*r1^31*r2 + 321*r1^30*r2^2 - 1406*r1^29*r2^3 - 37418*r1^28*r2^4 - 54378*r1^27*r2^5 + 1866549*r1^26*r2^6 + 5556726*r1^25*r2^7 - 62431686*r1^24*r2^8 - 193177494*r1^23*r2^9 + 1667661669*r1^22*r2^10 + 3566941110*r1^21*r2^11 - 36334402551*r1^20*r2^12 - 22380759756*r1^19*r2^13 + 596476955268*r1^18*r2^14 - 616148677152*r1^17*r2^15 - 6288496227900*r1^16*r2^16 + 18622534127616*r1^15*r2^17 + 24103256904576*r1^14*r2^18 - 214678876124160*r1^13*r2^19 + 294385908328704*r1^12*r2^20 + 800063891564544*r1^11*r2^21 - 3753687721439232*r1^10*r2^22 + 5063157982642176*r1^9*r2^23 + 4824832605788160*r1^8*r2^24 - 33676529335271424*r1^7*r2^25 + 73507309933658112*r1^6*r2^26 - 99895501823016960*r1^5*r2^27 + 93672133367169024*r1^4*r2^28 - 61376067146612736*r1^3*r2^29 + 27116508430467072*r1^2*r2^30 - 7312316880125952*r1*r2^31 + 914039610015744*r2^32)

項1，項2 を因数分解した結果は，かなり短い式になる。

term1 = -r1^16 - 20*r1^15*r2 - 70*r1^14*r2^2 + 876*r1^13*r2^3 + 5148*r1^12*r2^4 - 16974*r1^11*r2^5 - 110601*r1^10*r2^6 + 293274*r1^9*r2^7 + 1113426*r1^8*r2^8 - 4030560*r1^7*r2^9 - 2711232*r1^6*r2^10 + 27387072*r1^5*r2^11 - 34572096*r1^4*r2^12 - 22394880*r1^3*r2^13 + 95738112*r1^2*r2^14 - 90699264*r1*r2^15 + 30233088*r2^16;
term1 = term1 |> factor
term1 |> println

   -(r1 - 2*r2)^2*(r1 + 6*r2)^6*(r1^2 - 3*r2^2)*(r1^2 - 6*r1*r2 + 6*r2^2)*(r1^2 - 3*r1*r2 + 3*r2^2)^2

term2 = r1^32 + 34*r1^31*r2 + 321*r1^30*r2^2 - 1406*r1^29*r2^3 - 37418*r1^28*r2^4 - 54378*r1^27*r2^5 + 1866549*r1^26*r2^6 + 5556726*r1^25*r2^7 - 62431686*r1^24*r2^8 - 193177494*r1^23*r2^9 + 1667661669*r1^22*r2^10 + 3566941110*r1^21*r2^11 - 36334402551*r1^20*r2^12 - 22380759756*r1^19*r2^13 + 596476955268*r1^18*r2^14 - 616148677152*r1^17*r2^15 - 6288496227900*r1^16*r2^16 + 18622534127616*r1^15*r2^17 + 24103256904576*r1^14*r2^18 - 214678876124160*r1^13*r2^19 + 294385908328704*r1^12*r2^20 + 800063891564544*r1^11*r2^21 - 3753687721439232*r1^10*r2^22 + 5063157982642176*r1^9*r2^23 + 4824832605788160*r1^8*r2^24 - 33676529335271424*r1^7*r2^25 + 73507309933658112*r1^6*r2^26 - 99895501823016960*r1^5*r2^27 + 93672133367169024*r1^4*r2^28 - 61376067146612736*r1^3*r2^29 + 27116508430467072*r1^2*r2^30 - 7312316880125952*r1*r2^31 + 914039610015744*r2^32;
term2 = term2 |> factor
term2 |> println

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

分子を書き直し，因数分解する。

num = term1 -(r1 - 2*r2)^2*(r1 + 6*r2)^6*(r1^2 - 6*r1*r2 + 6*r2^2)*(r1^2 - 3*r1*r2 + 3*r2^2)^3 |> factor;
num |> println

   -r1*(r1 - 2*r2)^2*(r1 + 6*r2)^6*(2*r1 - 3*r2)*(r1^2 - 6*r1*r2 + 6*r2^2)*(r1^2 - 3*r1*r2 + 3*r2^2)^2

分子/分母 を計算し，簡約化する。

term3 = denominator(ans_r3);

ans__r3 = num/term3 |> simplify
ans__r3 |> println

   (-r1^2 + 6*r1*r2 - 6*r2^2)/(r1 - 2*r2)

得られた式 (-r1^2 + 6*r1*r2 - 6*r2^2)/(r1 - 2*r2) は術と同じである。

ans__r3(r1 => 4//2, r2 => 6//2) |> println

   11/2

術は，
人円径 = 地円径*(6地円径 - 4天円径)/(2地円径 - 天円径) - 天円径
であるが，簡約化すると
人円径 = (6*地円径^2 - 6*地円径*天円径 + 天円径^2)/(2*地円径 - 天円径)
になる。

@syms 天円径, 地円径
人円径 = 地円径*(6地円径 - 4天円径)/(2地円径 - 天円径) - 天円径;
人円径 |> simplify |> println
人円径(天円径 => 4, 地円径 => 6) |> println

   (6*地円径^2 - 6*地円径*天円径 + 天円径^2)/(2*地円径 - 天円径)
   11

# (ox, oy) を中心，r を半径とする円に内接する正 n 多角形。fcol を指定すると塗りつぶし。
function polygon(ox, oy, r, n; φ=90, color=:black, lt=:solid, lw=0.5, fillcolor="")
   θ = range(0, 2π, length=n+1) .+ pi*φ/180
   if fillcolor == ""
       plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
             linecolor=color, linestyle=lt, linewidth=lw)
   else
       plot!(r .* cos.(θ) .+ ox, r .* sin.(θ) .+ oy,
             linecolor=color, linestyle=lt, linewidth=lw,
             seriestype=:shape, fillcolor=fillcolor)
   end
end;

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = √3r2
   r3 = (r1^2 - 6r1*r2 + 6r2^2)/(2r2 - r1)
   (b, x1, x3) = (a*(24*sqrt(3)*a^3 - 48*a^2*r1 - 4*sqrt(3)*a*r1^2 + 3*r1^3)/(r1*(-12*a^2 + 4*sqrt(3)*a*r1 + 3*r1^2)), a - sqrt(3)*r1, -a - sqrt(3)*r3)
   (x0, y0) = float.(intersection(b, 0, -a, 2√3a, 3a, 0, 2a, √3a))
   @printf("天円，地円の直径が %g, %g のとき，人円の直径は %g である。\n", 2r1, 2r2, 2r3)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  a = %g;  b = %g;  x1 = %g;  x3 = %g\n", r1, r2, r3, a, b, x1, x3)
   plot([b, x0, 3a, b], [0, y0, 0, 0], color=:green, lw=0.5)
   polygon(0, √3a, 2a, 6, φ=60, color=:magenta)
   circle(x1, 2√3a + r1, r1, :orange)
   circle(2a, r2, r2)
   circle(x3, r3, r3, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a, 0, "a", :green, :center, delta=-2delta)
       point(2a, 0, "2a", :green, :center, delta=-2delta)
       point(3a, 0, "3a", :green, :center, delta=-2delta)
       point(0, 2√3a, "2√3a", :magenta, :center, delta=-2delta)
       point(x1, 2√3a + r1, " 天円:r2,(x1,2√3a+r1)", :black, :left, :vcenter)
       point(2a, r2, " 地円:r2,(2a,r2)", :red, :left, :vcenter)
       point(x3, r3, "人円:r3\n(x3,r3)", :red, :center, delta=-delta)
       point(x0, y0, " (x0,y0)", :green, :left, :vcenter)
       point(b, 0, "b", :green, :center, delta=-2delta)
       plot!(ylims=(-8delta, y0 + 3delta), xlims=(b - 3delta, 3a + 25delta))
   end
end;

draw(4/2, 6/2, true)

晴屋製麺所

高松市今里町　晴屋製麺所

他に類を見ないぐらいの、かなりの細麺😀

揚げ過ぎじゃないけど、天ぷらの衣が堅い🤔

算額（その1391）

算額（その1391）

十八 大里郡岡部村岡 稲荷社
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円5個，斜線2本
#Julia, #SymPy, #算額, #和算

5 個の円が，2 本の直線に挟まれている。大円の直径が 2 寸のとき，中円の直径はいかほどか。

元の図を反時計回りに 90 °回転させたもので考える。
小円の位置が元図と違うが，大円と中円の中心を結ぶ直線に関して対称なのでこの図のままでよい。

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

eq1, eq2, eq3 は，大円，中円，右側の小円が互いに接しているという中心間の水平距離に関する条件式である。
eq4 は，大円と中円の水平距離と中円と左側の小円の水平距離（灰色で描いた相似の直角三角形の底辺と高さ）に関する条件式である。

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

using SymPy

@syms r1::positive, x1::positive,
     r2::positive, r3::positive,
     x3::positive

eq1 = x1 - 2sqrt(r1*r2)
eq2 = (x1 - x3) - 2sqrt(r1*r3)
eq3 = x3 - 2sqrt(r2*r3)
eq4 = (r1 - r2)/x1 - (r2 - r3)/(x3 + 2r3);

一度に解けないので，まず eq1, eq2, eq3 を解いて x1, r3, x3 を求める。

res = solve([eq1, eq2, eq3], (x1, r3, x3))[1]

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

eq4 に x1, r3, x3 を代入し（eq14 とする），r2 について解く。

eq14 = eq4(x1 => res[1], r3 => res[2], x3 => res[3])
ans_r2 = solve(eq14, r2)[1]
ans_r2 |> println

   r1/2

中円の半径 r2 は，大円の半径の 1/2 である。
大円の直径が 2 寸のとき，中円の直径は 1 寸である。

いきなり solve() でよいが， どういうふうに解かれているかは，以下のようなものでもあろう。
まず，x1, r2, x3 を代入した方程式を簡約化し，分数式になったので，その分子を取り出して，分子 = 0 を解く。

eq14 = eq4(x1 => res[1], r3 => res[2], x3 => res[3]) |> simplify |> numerator
eq14 |> println

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

式が因数分解できれば，解がすぐわかる。ルートを含むそのままの式は SymPy では因数分解できないので，変数を置き換えて因数分解した後元の変数に戻す。

@syms s, t
eq14(√r1 => s, √r2 => t) |> factor |> x -> x(s => √r1, t => √r2) |> println

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

方程式の解は，r1 - 2*r2 = 0，つまり，r2 = r1/2 である。

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

  r1 = 1;  r2 = 1.41421;  x1 = 0.5;  r3 = 0.171573;  x3 = 0.585786

function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = r1/2
   (x1, r3, x3) = (2*sqrt(r1)*sqrt(r2), r1*r2/(sqrt(r1) + sqrt(r2))^2, 2*sqrt(r1)*r2/(sqrt(r1) + sqrt(r2)))
   @printf("大円の直径が %g のとき，中円の直径は %g である。\n", 2r1, 2r2)
   @printf("r1 = %g;  x1 = %g;  r2 = %g;  r3 = %g;  x3 = %g\n", r1, x1, r2, r3, x3)
   plot()
   circle
   circle(x1, r1, r1)
   circle(0, r2, r2, :green)
   circle2(x3, r3, r3, :blue)
   circle(-x3 - 2r3, r3, r3, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       plot!([-x3 - 2r3, 0, 0, -x3 - 2r3], [r3, r3, r2, r3], color=:gray60, lw=0.5)
       plot!([0, x1, x1, 0], [r2, r2, r1, r2], color=:gray60, lw=0.5)
       point(x1, r1, "大円:r1,(x1,r1)", :red, :center, :bottom, delta=delta)
       point(0, r2, "中円:r2,(0,r2)", :green, :center, :bottom, delta=4delta)
       point(x3, r3, "小円:r3,(x3,r3)", :blue, :center, delta=-7delta)
       point(-x3 - 2r3, r3, "小円:r3,(-x3-2r3,r3)", :blue, :center, delta=-7delta, deltax=12delta)
       θ = 2atand((r1 - r2)/x1)
       (x01, y01) = (x1 - r1*sind(θ), r1 + r1*cosd(θ))
       (x02, y02) = (-r2*sind(θ), r2 + r2*cosd(θ))
       println(θ)
       abline(x01, y01, tand(θ), -1.2, 2.2)
       segment(-1.2, 0, 2.2, 0)
       ylims!(-8delta, 2r1 + 3delta)
   end
end;

draw(2/2, true)

 

算額（その1390）

算額（その1390）

十七 大里郡岡部村岡 稲荷社
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円8個，直線上，斜線2本
#Julia, #SymPy, #算額, #和算

直線と斜線 2 本を引き，円 8 個を載せる。丁円の直径が 33 寸のとき，丙円の直径はいかほどか。

甲円の半径と中心座標を r1, (-r1, r1), (r1, r1)
乙円の半径と中心座標を r2, (0, 2r4 + r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (0, r4), (2r1, r4)
とおく。

甲円と丁円が接することから，関係式 eq1 が成り立つ。
それを解けば，r1 は r4 の 4 倍であることが簡単にわかる。

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

using SymPy

@syms r1::positive, r2::positive,
     r3::positive, x3::positive,
     r4::positive, a::positive

eq1 = r1 - 2sqrt(r1*r4)
solve(eq1, r1)[1] |> println

   4*r4

甲円と同じ大きさの補助円（灰色）を描くと，左側の甲円の中心 (-r1, r1) と右側の丁円の中心 (2r1, r4) と補助円と x 軸の接点 (3r1, 0) は一直線上にある。

直角三角形 ⊿ABC と ⊿DBE の相似関係から sqrt((r1 + a)^2 + 16r1^2) = 4a が成り立つので a = 68r4/15 がわかる。

eq2 = sqrt((r1 + a)^2 + 16r1^2) - 4a
solve([eq1, eq2], (a, r1))[1] |> println

   (68*r4/15, 4*r4)

乙円と丙円，丙円と丁円が外接するので，以下の関係式が成り立つ。

eq3 = x3^2 + (2r4 + r2 - r3)^2 - (r2 + r3)^2
eq4 = x3^2 + (r3 - r4)^2 - (r3 + r4)^2;

左側の甲円の中心，乙円の中心，右側の丁円の中心 から (-r1, r1 + a) と (3r1, 0) を結ぶ直線までの距離はそれぞれ r1, r2, r4 である。
この中から，これまでの eq1, eq2, eq3 ,eq4 と独立な関係式を用いる。

eq5 = dist2(-r1, r1 + a, 3r1, 0, 0, 2r4 + r2, r2);
# eq6 = dist2(-r1, r1 + a, 3r1, 0, -r1, r1, r1)
# eq7 = dist2(-r1, r1 + a, 3r1, 0, 2r1, 0, r4);

以上の 5 本の方程式を連立させて解き，a, r1, r2, r3, x3 を求める。

solve([eq1, eq2, eq3, eq4, eq5], (a, r1, r2, r3, x3))[1]

   (68*r4/15, 4*r4, 33*r4/16, 49*r4/33, 14*sqrt(33)*r4/33)

丙円の半径 r3 は，丁円 の半径の 49/33 倍である。
丁円の直径が 33 寸のとき，丙円の直径は 49 寸である。

全てのパラメータは以下の通りである。

  r4 = 16.5;  a = 74.8;  r1 = 66;  r2 = 34.0312;  r3 = 24.5;  x3 = 40.2119

function draw(r4, more)
    pyplot(size=(800, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (a, r1, r2, r3, x3) = r4 .* (68/15, 4, 33/16, 49/33, 14/√33)
   @printf("r4 = %g;  a = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g\n", r4, a, r1, r2, r3, x3)
   plot()
   circle2(r1, r1, r1)
   circle(3r1, r1, r1, :gray60)
   circle(0, r4, r4, :blue)
   circle2(2r1, r4, r4, :blue)
   circle(0, 2r4 + r2, r2, :orange)
   circle2(x3, r3, r3, :green)
   segment(-r1, r1 + a, 3r1, 0, :magenta)
   segment(r1, r1 + a, -3r1, 0, :magenta)
   segment(-r1, r1, 3r1, 0, :gray60)
   segment(-r1, 0, -r1, a + r1)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(-r1, r1 + a, "B:(-r1,r1+a) ", :black, :right, :vcenter)
       point(-r1, r1, "C:(-r1,r1) ", :red, :right, :vcenter)
       point(3r1, 0, "E:3r1", :red, :center, delta=-3delta)
       point(2r1, 0, "2r1", :red, :center, delta=-3delta)
       point(r1, r1, "甲円:r1,(r1,r1)", :red, :left, :bottom, delta=2delta)
       point(0, 2r4 + r2, "乙円:r2\n(0,2r4+r2)", :orange, :left, :vcenter, deltax=delta)
       point(0, r4, "丁円:r4\n(0,r4)", :blue, :center, :vcenter)
       point(2r1, r4, "丁円:r4\n(2r1,r4)", :blue, :center, :vcenter)
       point(x3, r3, " 丙円:r3,(x3,r3)", :green, :left, :vcenter)
       θ = atand((r1 + a)/4r1)
       segment(-r1, r1, -r1 + r1*sind(θ), r1 + r1*cosd(θ), :gray60)
       point(-r1 + r1*sind(θ), r1 + r1*cosd(θ), " A", :green, :left, :bottom)
       point(-r1, 0, "D:-r1", :green, :center, delta=-3delta)
       plot!(xlims=(-5delta - 2r1 - r4, 4r1 + 5delta), ylims=(-10delta, r1 + a + 3delta))
   end
end;

draw(33/2, true)

算額（その1389）

算額（その1389）

 

番外九 武州 慈恩寺

埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

キーワード：球5個，外球，3 次元，ソディ・ゴセットの n+1 球定理

#Julia, #SymPy, #算額, #和算

 

外球の中に小球を 4 個容れる。外球の直径が 1 寸のとき，小球の直径はいかほどか。

1.　ソディ・ゴセットの n+1 球定理による解答

 

算額（その1374）で用いた方法である。

 

外球の半径を R，小球の半径を r とおき，以下の方程式をとき r を求める。

 

using SymPy

 

@syms r::positive, R::positive

eq = 3(4(1/r^2) + 1/R^2) - (4/r - 1/R)^2;

 

solve(eq, r)[1] |> println

 

   R*(-2 + sqrt(6))

 

小球の半径 r は，外球の半径 R の (√6 - 2) 倍である。

外球の直径が 1 寸のとき，小球の直径は √6 - 2 = 0.4494897427831779 寸である。

 

答えは「4分4厘有奇」と正確とは言い難いが，術は「6 の平方根から 2 を引いて外球の直径を掛ける」とあるので正確である。誤差の主原因は √6 に不正確な近似値 2.44 を使用しためであろうか。2.44 - 2 =0.44

 

s6 = 2.44

s6 - 2

 

   0.43999999999999995

 

2.　基礎的なところから計算する解答

 

半径が r の 4 つの小球の中心は，一辺が a = 2r の正四面体の頂点である。

小球が内接する外球の半径 R は、各球の中心から正四面体の重心までの距離に加え、各球の半径を加えることで求められる。

 

1. 一辺の長さが a の四面体において，重心から各頂点までの距離 d は

    d = sqrt(3/8)*a

2. 内接する球の半径 R は

    R = d + r

      = sqrt(3/8)*a + r

      = sqrt(3/8)*2r + r

      = (sqrt(3/8)*2 + 1)*r

         = (1 + √6/2)*r

   3. また，小球の半径 r は

       r = R/(sqrt(3/8)*2 + 1)

         = (√6 - 2)*R

 

外球の半径が R = 0.5 のとき，小球の半径は r = √6 - 2 = 0.22474487139158905

直径はその 2 倍の 2*0.22474487139158905 = 0.4494897427831781 である。

 

算額（その1387）

算額（その1387）

七二 加須市大字外野 棘脱地蔵堂 明治7年(1874)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円2個，直角三角形，正方形2個
#Julia, #SymPy, #算額, #和算

直角三角形の中に，正方形 2 個，日円，月円を容れる。月円の直径が 159 寸のとき，日円の直径はいかほどか。

注1：この問題を解くにあたっては，鈎：股の比の情報が必要である。それが与えられれば，日円と月円を含む直角三角形の相似比が決まり，日円と月円の相似比も同じなので，日円の直径も簡単に！決まる。

注2：推方（傾いた正方形？）は意味のない情報である。なので，求めもしないし，図にも描かない。

鈎，股をそのまま変数名「鈎」，「股」
正方形の一辺の長さを a
日円の半径と中心座標を r1, (r1, a + r1)
月円の半径と中心座標を r2, (a + r2, r2)
とおき，以下の連立方程式を「月円を定数として」解く。

術は「月円の直径を三倍して四で割る」と，身も蓋もないことをいっている。
つまり，日円を含む直角三角形と月園を含む直角三角形の相似比，ひいては元の直角三角形の鈎と股の比が 3:4 であることを前提としている。上の図も 3:4 のときのものである。相似比が違えば，見かけの異なる図になる。

この時代（?）は，鈎股弦（直角三角形）といえば辺の比が 3:4:5 であるものと，暗黙の前提としていたのであろう。

ということで，eq2 は 鈎/股 = 3/4 を表している。

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

using SymPy

@syms 鈎::positive, 股::positive, a::positive,
     r1::positive, r2::positive
eq1 = (鈎 - a)/a - 鈎/股
eq2 = 鈎/股 - 3//4
eq5 = r1/r2 - a/(股 - a)
eq6 = sqrt((鈎 - a)^2 + a^2) + sqrt(a^2 + (股 - a)^2) - sqrt(鈎^2 + 股^2)
eq3 = (鈎 - a) + a - sqrt((鈎 - a)^2 + a^2) - 2r1
eq4 = a + (股 - a) - sqrt(a^2 + (股 - a)^2) - 2r2;

得られる解は非常にきれいな結果である。

日円の半径 r1 は，月円の半径 r2 の 3/4 倍である。

res = solve([eq2, eq3, eq4, eq5], (鈎, 股, a, r1))[1]

   (21*r2/4, 7*r2, 3*r2, 3*r2/4)

その結果，「月円の直径が 159 寸」のとき，「日円の直径は 159*3/4 = 119.25 寸」が得られる。

一般的に，算額の問は「きれいな数の答えが得られる」ように設定されることが多い。

もし，問が「月円の直径が 160 寸」ならば，「日円の直径は 160*3/4 = 120 寸」が得られるのになぜそうしなかったのか。

鈎股弦の比が 3:4:5 というのを隠そうとしたのか。

人心を惑わす，問題に無関係な推方（傾いた正方形？）を持ち出した理由も含め，謎は尽きない。

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, a, r1) = (21*r2/4, 7*r2, 3*r2, 3*r2/4)
   @printf("月円の直径が %g のとき，日円の直径は %g である。\n",  2r2, 2r1)
   @printf("r2 = %g;  鈎 = %g;  股 = %g;  a = %g;  r1 = %g\n", r2, 鈎, 股, a, r1)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:green, lw=0.5)
   rect(0, 0, a, a, :blue)
   circle(r1, a + r1, r1)
   circle(a + r2, r2, r2, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, 鈎, "鈎 ", :green, :right, :vcenter)
       point(股, 0, "股", :green, :center, delta=-delta)
       point(a, a, "(a,a)", :blue, :left, :bottom, delta=delta)
       point(r1, a + r1, "日円:r1\n(r1,a+r1)", :red, :center, delta=-delta)
       point(a + r2, r2, "月円:r2\n(a+r2,r2)", :magenta, :center, delta=-delta)
       plot!(xlims=(-6delta, 股 + 5delta), ylims=(-6delta, 鈎 + 5delta))
   end
end;

draw(159/2, true)

 

算額（その1386）

算額（その1386）

五九 北埼玉郡南河原村南河原 観福寺 万延二年(1861)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：計算
#Julia, #SymPy, #算額, #和算

甲, 乙, 丙の正方形がある。面積の合計は 1133 である。一辺の長さの差は 5 である（大きさは，甲，乙，丙の順）。
それぞれの一辺の長さを求めよ。

甲, 乙, 丙の一辺の長さをそのまま変数名として，以下の連立方程式を解く。

using SymPy

@syms 甲::positive, 乙::positive, 丙::positive

eq1 = (甲 - 乙)　- 5
eq2 = (乙 - 丙)　- 5
eq2 = 甲^2 + 乙^2 + 丙^2 - 1133
eq3 = sum([甲, 乙, 丙].^2) - 1133
solve([eq1, eq2, eq3], (甲, 乙, 丙))[1]

   (24, 19, 14)

甲 = 24, 乙 = 19, 丙 = 14 である。

検算

sum([24, 19, 14].^2)

   1133

算額（その1385）

算額（その1385）

四五 秩父郡吉田町大字久長字頼母沢 野栗神社 安政3年(1856)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：計算
#Julia, #SymPy, #算額, #和算

甲, 乙, 丙, 丁, 戊, 己の立方体がある。体積の合計は 33291 である。
乙の体積は甲の体積の 7/9
丙の体積は乙の体積の 5/7
丁の体積は丙の体積の 3/5
戊の体積は丁の体積の 2/3
己の体積は戊の体積の 1/2
のとき，それぞれの一辺の長さを求めよ。

甲, 乙, 丙, 丁, 戊, 己の一辺の長さをそのまま変数名として，以下の連立方程式を解く。

using SymPy

@syms 甲::positive, 乙::positive, 丙::positive,
     丁::positive, 戊::positive, 己::positive

eq1 = 乙 - 7//9 * 甲
eq2 = 丙 - 5//7 * 乙
eq3 = 丁 - 3//5 * 丙
eq4 = 戊 - 2//3 * 丁
eq5 = 己 - 1//2 * 戊
eq6 = 甲^3 + 乙^3 + 丙^3 + 丁^3 + 戊^3 + 己^3 - 33291
eq6 = sum([甲, 乙, 丙, 丁, 戊, 己].^3) - 33291
solve([eq1, eq2, eq3, eq4, eq5, eq6], (甲, 乙, 丙, 丁, 戊, 己))[1]

   (27, 21, 15, 9, 6, 3)

甲 = 27, 乙 = 21, 丙 = 15, 丁 = 9, 戊 = 6, 己 = 3 である。

検算

sum([27, 21, 15, 9, 6, 3].^3)

   33291

 

算額（その1384）

算額（その1384）

一八 大里郡岡部村岡 稲荷社 文化14年(1817)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円9個，外円，弦2本
#Julia, #SymPy, #算額, #和算

外円の中に 2 本の弦を引き，甲円 1 個，乙円 1 個，丙円 6 個を容れる。甲円，乙円の直径が 22 寸，11 寸のとき，丙円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, r1 - R)
乙円の半径と中心座標を r2, (0, R - r2)
丙円の半径と中心座標を r3, (x31, y31), (x32, y32), (x33, y33)
とおき，以下の連立方程式を解く。

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

using SymPy

@syms R::positive, x01::positive, x02::positive,
     r1::positive, r2::positive, r3::positive,
     x31::positive, y31::positive,
     x32::positive, y32::positive,
     x33::positive, y33::positive;

R = r1 + r2
eq1 = x31^2 + y31^2 - (R - r3)^2
eq2 = x32^2 + y32^2 - (R - r3)^2
eq3 = x33^2 + y33^2 - (R - r3)^2
eq4 = (x31 - x32)^2 + (y31 - y32)^2 - 4r3^2
eq5 = (x31 - x33)^2 + (y31 - y33)^2 - 4r3^2
eq6 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, r1 - R) - r1^2
eq7 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, R - r2) - r2^2
eq8 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x32, y32) - r3^2
eq9 = dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x33, y33) - r3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (x01, x02, r3, x31, y31, x32, y32, x33, y33))

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=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return Float64.(v), r.f_converged
end;

function H(u)
   (x01, x02, r3, x31, y31, x32, y32, x33, y33) = u
   return [
       x31^2 + y31^2 - (R - r3)^2,
       x32^2 + y32^2 - (R - r3)^2,
       x33^2 + y33^2 - (R - r3)^2,
       (x31 - x32)^2 + (y31 - y32)^2 - 4r3^2,
       (x31 - x33)^2 + (y31 - y33)^2 - 4r3^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, r1 - R) - r1^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), 0, R - r2) - r2^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x32, y32) - r3^2,
       dist(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2), x33, y33) - r3^2
   ]
end;
(r1, r2) = (22, 11) ./ 2
R = r1 + r2

iniv = BigFloat[4, 13, 3.1, 13, 4, 13, -1.5, 10, 10]
res = nls(H, ini=iniv)

   ([4.069279408445208, 13.215553020559287, 3.0, 12.727922061357855, 4.5, 13.420835425335575, -1.459854953773714, 9.520851253161299, 9.570966064884825], true)

甲円，乙円の直径が 22 寸, 11 寸のとき，丙円の直径は 6 寸である。

全てのパラメータは以下のとおりである。
   R = 16.5;  x01 = 4.06928;  x02 = 13.2156;  r1 = 11;  r2 = 5.5;  r3 = 3;  x31 = 12.7279;  y31 = 4.5;  x32 = 13.4208;  y32 = -1.45985;  x33 = 9.52085;  y33 = 9.57097

function draw(r1, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
  (x01, x02, r3, x31, y31, x32, y32, x33, y33) = res[1]
   R = r1 + r2
   @printf("甲円，乙円の直径が %g, %g のとき，丙円の直径は %g である。\n", 2r1, 2r2, 2r3)
   @printf("R = %g;  x01 = %g;  x02 = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x31 = %g;  y31 = %g;  x32 = %g;  y32 = %g;  x33 = %g;  y33 = %g\n",
       R, x01, x02, r1, r2, r3, x31, y31, x32, y32, x33, y33)
   plot()
   circle(0, 0, R, :magenta)
   circle(0, r1 - R, r1, :blue)
   circle(0, R - r2, r2, :green)
   circle2(x31, y31, r3)
   circle2(x32, y32, r3)
   circle2(x33, y33, r3)
   segment(x01, sqrt(R^2 - x01^2), x02, -sqrt(R^2 - x02^2))
   segment(-x01, sqrt(R^2 - x01^2), -x02, -sqrt(R^2 - x02^2))
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R - r2, "乙円:r2,(0,R-r2)", :black, :center, :bottom, delta=delta)
       point(0, r1 - R, "甲円:r1,(0,r1-R)", :blue, :center, delta=-delta)
       point(0, R, "R", :green, :center, :bottom, delta=delta)
       point(x01, sqrt(R^2 - x01^2), "(x01,y01)", :black, :left, :bottom, delta=delta)
       point(x02, -sqrt(R^2 - x02^2), "(x02,y02)", :black, :left, delta=-delta, deltax=-delta)
       point(x31, y31, "丙円:r3\n(x31,y31)", :red, :center, :bottom, delta=delta/2)
       point(x32, y32, "丙円:r3\n(x32,y32)", :red, :center, :bottom, delta=delta/2)
       point(x33, y33, "丙円:r3\n(x33,y33)", :red, :center, :bottom, delta=delta/2)
   end
end;

draw(22/2, 11/2, true)