算額（その1551）
百四十六 群馬県新田郡新田町下田中 田中神社 大正6年(1917)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：連立方程式
#Julia, #SymPy, #算額, #和算，#数学
合計 244 坪の土地を，甲乙 2 つに分筆する。（面積を正方形に換算すると，）乙の正方形の一辺は甲のそれよりも 2 間短い。甲乙それぞれの一辺の長さはいかほどか。
using SymPy
@syms 甲::positive, 乙::positive
eq1 = 甲^2 + 乙^2 - 244
eq2 = (甲 - 2) - 乙
solve([eq1, eq2], (甲, 乙))[1]
(12, 10)
算額（その1550）
百二十四 群馬県桐生市天神町 天満宮 明治11年(1878)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：連立方程式，整数方程式
#Julia, #SymPy, #算額, #和算，#数学
桃と梨を買った。桃と梨の両方合わせた個数と代金（の数値）が同じになった。ただし，桃は 9 個で代金は 5 文，梨は 3 個で 25 文である。個数と代金を求めよ。
using SymPy
@syms 桃, 梨
eq = (5(桃/9) + 25(梨/3)) - (桃 + 梨)
数学の試験ならいざ知らず，ブルートフォースで一発即答。
for 桃 = 9:9:1000, 梨 = 3:3:1000
(5(桃/9) + 25(梨/3)) == (桃 + 梨) && (println("桃 = $桃 個，代金 = $(5(桃/9))\n梨 = $梨 個，代金 = $(25(梨/3))"); break)
end
桃 = 99 個，代金 = 55.0
梨 = 6 個，代金 = 50.0
算額（その1549）
百三 群馬県高崎市八幡町 八幡宮 安政7年(1860)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：数列，部分和
#Julia, #SymPy, #算額, #和算，#数学
逐増する「平方垜」がある，第 1 項（底子 1）では 1 個，第 2 項までの合計は 6 個（5 個増えたことになる。この増え方が平方垜になる），第 3 項までの合計は 20個（14 個増えたことになる）であるとき，第 n 項までの合計は何個か。
注：垜（だ）とは，球，土嚢などを一番下から順次規則的に積み上げた構造物である。
規則はいろいろあるが
1. 三角数の積み上げ
三角数は、上から順に**1、2、3…**と増加する層で構成される。Tn = n(n + 1)/2
2. 四角数の積み上げ
四角数は正方形の形を形成する。Qn = n^2
3. 五角数の積み上げ
五角数は星型や角形を意識した積み方Pn = n(3n - 1)/2
4. 六角数（蜂の巣型） Hn = n(2n - 1)
など
平方垜は四角数の積み上げと思いきや，合計が 1, 6, 20, ... となり，差分が（問題文にもあるが）5, 14 となる。例が少ないが前を補うと 0, 1, 5, 14, ... で，更に差分を取ると 1, 4, 9, ...となるこれは平方数の数列であろうと推測できる。
当時の人は，平方垜といえばすぐにこの数列を書き下せるほど，常識だったのであろう。
ひとまずその仮定に基づくと，平方垜の第 n 項までの合計の数列は，
0, 1, 6, 20, 50, 105, ... 合計
0, 1, 5, 14, 30, 55, ... 各段の個数
1, 4, 9, 16, 25, ... 差分
「オンライン整数列大辞典」を参照すれば以下のものが見つかる。
A002415 4-dimensional pyramidal numbers: a(n) = n^2*(n^2-1)/12.
0, 0, 1, 6, 20, 50, 105, 196, 336, ...
ただし，この数式は n = 0 から始まっているので，その調整をする必要がある。
using SymPy
@syms n
an = n^2*(n^2-1)/12 # A002415
an |> println
n^2*(n^2 - 1)/12
an(n => n + 1) |> factor|> println # 項の調整
n*(n + 1)^2*(n + 2)/12
各項を計算して，推定どおりであることを確認する。
for n = 1:10
println("n = $n, an = $(n*(n + 1)^2*(n + 2)/12)")
end
n = 1, an = 1.0
n = 2, an = 6.0
n = 3, an = 20.0
n = 4, an = 50.0
n = 5, an = 105.0
n = 6, an = 196.0
n = 7, an = 336.0
n = 8, an = 540.0
n = 9, an = 825.0
n = 10, an = 1210.0
「術」は n*(n^3 +4n^2 + 5n + 2)/12 であるが，因数分解すると同じであることがわかる。
n*(n^3 +4n^2 + 5n + 2)/12 |> factor |> println
n*(n + 1)^2*(n + 2)/12
整数列大辞典に見つからないような場合は，以下のように多項式の係数を推定するというアプローチを取る。
using SymPy
@syms a1, a2, a3, a4
f(n) = a1*n + a2*n^2 + a3*n^3 + a4*n^4;
eq1 = f(1) - 1
eq2 = f(2) - 6
eq3 = f(3) - 20
eq4 = f(4) - 50
res = solve([eq1, eq2, eq3, eq4], (a1, a2, a3, a4))
res |> println
Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(a3 => 1/3, a1 => 1/6, a2 => 5/12, a4 => 1/12)
一般項を表す，同じ多項式が得られた。
@syms n
res[a1]*n + res[a2]*n^2 + res[a3]*n^3 + res[a4]*n^4 |> factor |> println
n*(n + 1)^2*(n + 2)/12
一般項がわかっている場合の和の数列は summatin 関数で求めることができる。
using SymPy
@syms n, m
# 一般項
# an = n*(n + 1)*(2n + 1)/6
# 和の一般項
summation(n*(n + 1)*(2n + 1)/6, (n, 1, m)) |> factor |> println
m*(m + 1)^2*(m + 2)/12
算額（その1548）
神壁算法 皇都 安井聖天堂 文化2年()
藤田貞資門人 皇都粟田 町野左一郎好謙
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
キーワード：数列
#Julia, #SymPy, #算額, #和算，#数学
最上段に 1 個，その下に 3 個，更にその下に 6 個，次いで，10，15...のように玉が積まれている。玉の重さは一番上が最も重く，以下順に 2 ずつ軽くなっている。玉の重さは全部で 6860 である。玉の個数は全部で 56 個である。一番上の玉の重さはいかほどか。
一般項を求めるのが肝要である。
一番上の玉の重さを w, n 段あるとする。
段ごとの玉の数は，1, 3, 6, 10, ... で，n 段目にある玉の数は n*(n + 1)/2 である。
最上段から n 段目までにある玉の重さの合計は n*(n + 1)*(n + 2)/6 * w - (n - 1)*n*(n + 1)*(n + 2)/4 である。
最上段から各段までの玉の数の合計は，1, 4, 10, 20, ... で，n*(n + 1)*(n + 2)/6 である。
最上段から n 段目までにある玉の重さの合計を K1
最上段から各段までの玉の数の合計を K2
とおき，以下の連立方程式を解く。
K1, K2 を実地数とすれば容易に解ける。
問の場合は，最上段の重さは 130 である。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms w::positive, n::positive, K1::positive, K2::positive
eq1 =n*(n + 1)*(n + 2)/6 * w - (n - 1)*n*(n+1)*(n+2)/4 - 6860
eq2 = n*(n + 1)*(n + 2)/6 - 56;
res = solve([eq1, eq2], (n, w))[1]
(6, 130)
K1, K2 を記号のまま解くと虚数解になるが，虚部はほとんど 0 なので，実部だけをとればよい。
eq1 =n*(n + 1)*(n + 2)/6 * w - (n - 1)*n*(n+1)*(n+2)/4 - K1
eq2 = n*(n + 1)*(n + 2)/6 - K2;
res = solve([eq1, eq2], (n, w))[1] # 1 of 3
((-(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(2/3)/3 - (-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) - 3^(1/6)*I*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(2/3) + 3^(5/6)*I*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 2)/((3^(1/3) - 3^(5/6)*I)*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3)), (I*K1*(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 3^(5/6)*K1*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) - I*K2*(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(2/3)/2 - 3*I*K2*(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 3*3^(1/6)*K2*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(2/3)/2 - 3*3^(5/6)*K2*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 3*I*K2)/(K2*(3^(5/6) + 3^(1/3)*I)*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3)))
res[1](K1 => 6860, K2 => 56).evalf() |> println
res[2](K1 => 6860, K2 => 56).evalf() |> println
6.0 + 5.47542772767019e-26*I
130.0 + 3.28284985518336e-24*I
15 段目までの球の総数は 680 個，重さの合計は 74120 である。
res[1](K1 => 74120, K2 => 680).evalf() |> println
res[2](K1 => 74120, K2 => 680).evalf() |> println
15.0 - 2.03553346724214e-19*I
130.0 - 1.17148914145633e-18*I
算額（その1547）
神壁算法 播州垂井 住吉大明神社 寛政12年(1800)
藤田貞資門人 播州小野 高瀬恒右衛門信之
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
キーワード：円6個，面積
#Julia, #SymPy, #算額, #和算，#数学
大円 3 個が交わっており，小円 3 個を容れる。黒積が 3329.29のとき，小円の直径はいかほどか。
大円の半径と中心座標を r1, (0, r1), (√3r1, -r1/2)
小円の半径と中心座標を r2, (x2, x2/√3)
とおく。
まず，r2 が与えられたとき，r1 を求める。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms r1::positive, x1::positive, y1::negative, r2::positive, x2::positive, y2::negative
x1 = √Sym(3)r1/2
y1 = -r1/2
y2 = -x2/√Sym(3)
eq1 = x1^2 + (2r1 - r2 - y1)^2 - (r1 + r2)^2
eq2 = x2^2 + (r1 - y2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, x2))[1]
(7*r2/6, 2*sqrt(3)*r2/3)
大円の半径は小円の半径の 7/6 倍である。
黒積は，図の灰色部分の 6 倍である。
灰色部分の面積は，(大円の面積/2 - 小円の面積/2) - 2(大円の面積/6 - 正三角形OBCの面積) である。
@syms r1, r2, 黒積, d
r1 = 7r2/6
eq = 6*((PI*r1^2/2 - PI*r2^2/2) - 2(PI*r1^2/6 - r1^2*√Sym(3)/4)) - 黒積
eq |> println
-59*pi*r2^2/36 + 49*sqrt(3)*r2^2/12 - 黒積
方程式を解き，小円の半径 r2 = 6*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3)) である。
黒積が 3329.29 のとき，小円の直径は 12*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3)) = 83.20006146161401 である。
ans_r2 = solve(eq, r2)[2]
ans_r2 |> println
6*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3))
# r2
黒積 = 3329.29
12*sqrt(黒積)/sqrt(-59*pi + 147*sqrt(3))
83.20006146161401
「術」は「sqrt(3329.29/(sqrt(64827) - 59*π))*12 = 83.20006146161398」，上の結果と一致する。
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r2 = 83.20006146161401/2
(r1, x2) = (7*r2/6, 2*sqrt(3)*r2/3)
y2 = -x2/√3
x1 = √3r1
y1 = -r1/2
円周率 = 3.16
円周率 = π
S = 6(円周率*r1^2/2 - 円周率*r2^2/2 - 2(円周率*r1^2/6 - r1^2*√3/4))
println("r1 = $r1, r2 = $r2 のとき，黒積は $S")
plot()
rotate(0, r1, r1)
rotate(0, 2r1 - r2, r2, :blue)
θ = -30:0.1:90
x = r1.*cosd.(θ)
y = r1.*sind.(θ) .+ r1
θ = 90:-0.1:-90
append!(x, r2.*cosd.(θ))
append!(y, r2.*sind.(θ) .+ (2r1 - r2))
θ = 150:-0.1:90
append!(x, r1.*cosd.(θ) .+ √3r1/2)
append!(y, r1.*sind.(θ) .- r1/2)
plot!(x, y, color=:gray80, seriestype=:shape, fillcolor=:gray80, lw=0.5)
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, r1, "大円:r1,(0,r1)", :red, :center, delta=-delta/2)
point(0, 2r1 - r2, "小円:r2,(0,2r1-r2)", :blue, :center, :bottom, delta=delta)
point(0, 0, "O", :green, :center, delta=-1.5delta)
point(r1*sind(30), r1 - r1*cosd(30), "A", :green, :left, delta=-delta)
point(r1*sind(60), r1 - r1*cosd(60), "B", :green, :left, delta=-delta)
point(0, r1, "C ", :green, :right, :bottom, delta=-delta)
end
end;
draw(true)
算額（その1546）
神壁算法 越後州長岡 蒼柴大明神 寛政8年(1796)
關流藤田貞資門人 長岡 石垣作右衛門光隆
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf
キーワード：円2個，等脚台形，対角線
#Julia, #SymPy, #算額, #和算
外円の中に等脚台形と上円，下円の 2 個の円を容れる。内斜（対角線）は旁斜（斜辺）より 4 寸長い。上円の直径が 1 寸のとき，下円の直径はいかほどか。
外円の半径と中心座標を R, (0, 0)
上円の半径と中心座標を r1, (0, R - r1)
下円の半径と中心座標を r2,（0, r2 - R)
上頭（上底），下頭（下底） を 2b，2a
内斜（対角線）と 旁斜（斜辺）の差を K; 内斜 > 傍斜
とおき，以下の方程式を解く。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms R, a, b, r1, r2, K
a = sqrt(R^2 - (2r2 - R)^2)
b = sqrt(R^2 - (R - 2r1)^2)
eq1 = sqrt((a + b)^2 + (2R - 2r1 - 2r2)^2) - sqrt((a - b)^2 + (2R - 2r1 - 2r2)^2) - K
res = solve(eq1, r2)[1] # 1 of 2
res |> println
K^2/(16*r1)
不思議なことに（？），解の中に外円の半径 R が出てこない。つまり，内斜（対角線）と 旁斜（斜辺）の差 K と，上円の半径 r1 が決まれば r2 は自動的に決まり，R はなんでもよい（R により a, b が決まるから）ということだ。
それはさておき，下円の半径は『「内斜 - 旁斜」の二乗を上円の半径の16倍で割る』。
上円の半径が 1/2 のとき，下円の半径は 4^2/(16*(1/2)) = 2 である。術と一致する。
function draw(K, r1, R, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 1/2
r2 = K^2/(16r1)
@printf("「旁斜-内斜」が %g，上円の直径が %g のとき，下円の直径は %g である。外円の直径は %g である。\n", K, 2r1, 2r2, 2R)
b = sqrt(R^2 - (R - 2r1)^2)
a = sqrt(R^2 - (2r2 - R)^2)
println("a = $a, b = $b")
l1 = sqrt((a + b)^2 + (2R - 2r1 - 2r2)^2)
l2 = sqrt((a - b)^2 + (2R - 2r1 - 2r2)^2)
println((l1, l2, l1 - l2))
plot([a , b, -b, -a, a], [2r2 - R, R - 2r1, R - 2r1, 2r2 - R, 2r2 - R], color=:magenta, lw=0.5)
circle(0, 0, R, :blue)
circle(0, R - r1, r1)
circle(0, r2 - R, r2, :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(0, r2 - R, "下円:r1,(a,r2-R)", :green, :center, delta=-delta/2)
point(0, R - r1, "上円:r2,(b,R-r1)", :red, :center, delta=-delta/2)
point(a, 2r2 - R, "(a,2r2-R) ", :magenta, :right, delta=-delta/2)
point(b, R - 2r1, "(b,R-2r1) ", :magenta, :right, delta=-delta/2)
dimension_line(-b, R - 2r1, a, 2r2 - R, "内斜", :orange, :left, :vcenter, deltax=2delta)
dimension_line(b, R - 2r1, a, 2r2 - R, "旁斜", :tomato, :left, :vcenter, deltax=delta)
end
end;
draw(4, 1/2, 6, true)
算額（その1545）
埼玉県川口市三ツ和 氷川神社 享和4年(1804)
山口正義：やまぶき 第 39 号
https://yamabukiwasan.sakura.ne.jp/ymbk39.pdf
キーワード：円5個，直角三角形
#Julia, #SymPy, #算額, #和算
直角三角形の中に，甲，乙，丙，丁，戊の 5 円を容れる。直角三角形の股（底辺）と勾配が与えられたとき，各辺の長さと各円の直径を求めよ。
図形としては算額（その206）と同じである。何を与えて何を求めるかが違うだけである。
注：勾配といえば Δy/Δx であるが，当時はその 1/10 の値を使っていたのだろうか，また，勾配の有効桁が 4 桁しかないので，本来は整数値を意図していたはずのものに端数が出てしまう。そこで，鈎は問題文の「高倍（勾配）四寸弐分八厘五毛」を使わずに，3 とする。
鈎，股，弦をそのまま，「鈎」，「股」，「弦」
甲円の半径と中心座標を r1, (股 - r1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
戊円の半径と中心座標を r5, (x5, r5)
とおき，順次解く。
長いので，先に答えを書いておく。
鈎 = 3; 股 = 7; 弦 = 7.61577
直径： 甲円 = 2.38423; 乙円 = 1.58596; 丙円 = 1.05496; 丁円 = 0.701746; 戊円 = 0.466793
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms 鈎, 股, 弦, r1, r2, x2, r3, x3, r4, x4, r5, x5
弦 = sqrt(鈎^2 + 股^2);
# r1:甲円の半径
eq1 = (鈎 + 股 - 弦) - 2r1
ans_r1 = solve(eq1, r1)[1]
ans_r1 |> println
股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
ans_r1(鈎 => 3, 股 => 7).evalf() |> println
1.19211344706805
@syms x1, x2
r1 = ans_r1 # 1.19211344706805
x1 = 股 - r1
eq2 = (x1 - x2) - 2sqrt(r1*r2)
eq3 = r1/x1 - r2/x2;
res = solve([eq2, eq3], r2, x2)[1];
# r2:乙円の半径
ans_r2 = res[1]
ans_r2 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))^2
ans_r2(鈎 => 3, 股 => 7).evalf() |> println
0.792979085478491
# x2:乙円の中心の x 座標
ans_x2 = res[2]
ans_x2 |> println
(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
ans_x2(鈎 => 3, 股 => 7).evalf() |> println
3.86333413034970
累円の半径は，公比 p = r2/r1 = 0.665187602261171 の 等比数列をなす。
# p:公比
p = res[1]/ans_r1
p |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^2)
p(鈎 => 3, 股 => 7).evalf() |> println
0.665187602261170
# r3:丙円の半径
r3 = ans_r2*p
r3 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))^2*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^2/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^4)
r3(鈎 => 3, 股 => 7).evalf() |> println
0.527479856512693
# r4:丁円の半径
r4 = ans_r2*p^2
r4 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^2*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)
r4(鈎 => 3, 股 => 7).evalf() |> println
0.350873060994744
# r5:戊円の半径
r5 = ans_r2*p^3
r5 |> println
(股 + 鈎 - sqrt(股^2 + 鈎^2))^4*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^4/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^8)
r5(鈎 => 3, 股 => 7).evalf() |> println
0.233396410141131
互いに接する円の半径が既知ならば，中心間の水平距離は d12 = 2sqrt(r1*r2) であるが，算額の場合は，距離も等比級数に従う。
# d23:乙円と丙円の中心の水平距離
d23 = 2sqrt(ans_r2*r3)
d23 |> println
2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6))
d23(鈎 => 3, 股 => 7).evalf() |> println
1.29349216344864
# d34:丙円と丁円の中心の水平距離
d34 = 2sqrt(r3*r4)
d34 |> println
2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^5*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^5/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^10))
d34(鈎 => 3, 股 => 7).evalf() |> println
0.860414950748014
# d45:丁円と戊円の中心の水平距離
d45 = 2sqrt(r4*r5)
d45 |> println
2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^7*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^7/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^5*(股 - 鈎 + sqrt(股^2 + 鈎^2))^14))
d45(鈎 => 3, 股 => 7).evalf() |> println
0.572337358037734
# x3:丙円の中心の x 座標
x2 = res[2]
x3 = x2 - d23
x3 |> println
-2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)) + (3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
x3(鈎 => 3, 股 => 7).evalf() |> println
2.56984196690106
# x4:丁円の中心の x 座標
x4 = x3 - d34
x4 |> println
-2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^5*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^5/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^10)) - 2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)) + (3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
x4(鈎 => 3, 股 => 7).evalf() |> println
1.70942701615304
# x5:戊円の中心の x 座標
x5 = x4 - d45
x5 |> println
-2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^7*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^7/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^5*(股 - 鈎 + sqrt(股^2 + 鈎^2))^14)) - 2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^5*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^5/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^3*(股 - 鈎 + sqrt(股^2 + 鈎^2))^10)) - 2*sqrt((股 + 鈎 - sqrt(股^2 + 鈎^2))^3*(3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))^3/((股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)*(股 - 鈎 + sqrt(股^2 + 鈎^2))^6)) + (3*股^2 + 股*鈎 - 股*sqrt(股^2 + 鈎^2) + 3*鈎^2 - 3*鈎*sqrt(股^2 + 鈎^2) - 2*sqrt(2*股^4 + 4*股^3*鈎 - 2*股^3*sqrt(股^2 + 鈎^2) + 6*股^2*鈎^2 - 4*股^2*鈎*sqrt(股^2 + 鈎^2) + 4*股*鈎^3 - 4*股*鈎^2*sqrt(股^2 + 鈎^2) + 4*鈎^4 - 4*鈎^3*sqrt(股^2 + 鈎^2)))/(股 - 鈎 + sqrt(股^2 + 鈎^2))
x5(鈎 => 3, 股 => 7).evalf() |> println
1.13708965811531
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(鈎, 股) = (3, 7)
弦 = sqrt(鈎^2 + 股^2)
r1 = 1.19211344706805
(r2, x2) = (0.792979085478493, 3.86333413034969)
(r3, r4, r5) = (0.527479856512695, 0.350873060994746, 0.233396410141133)
(x3, x4, x5) = (2.56984196690104, 1.70942701615302, 1.13708965811529)
@printf("鈎 = %g; 股 = %g; 弦 = %g\n", 鈎, 股, 弦)
@printf("直径： 甲円 = %g; 乙円 = %g; 丙円 = %g; 丁円 = %g; 戊円 = %g\n", 2r1, 2r2, 2r3, 2r4, 2r5)
plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:red, lw=0.5)
circle(股 - r1, r1, r1, 1)
circle(x2, r2, r2, 2)
circle(x3, r3, r3, 3)
circle(x4, r4, r4, 4)
circle(x5, r5, r5, 5)
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)
end
end;
draw(true)
算額（その1544）
埼玉県川口市三ツ和 氷川神社 享和4年(1804)
山口正義：やまぶき 第 39 号
https://yamabukiwasan.sakura.ne.jp/ymbk39.pdf
キーワード：円3個，連立方程式
#Julia, #SymPy, #算額, #和算
大中小の 3 個の円がある。直径の差は大円と中円で 2 寸，大円と小円で 4 寸（等差数列），3 個の円の面積の和は 6557 である。各円の直径を求めよ。
注：算額の「問」では，6557 の単位は「分坪」つまり「平方分」なので半径を分で表さなければならないので注意。半径を寸で表すなら「65.57 平方寸」とする。
また，円周率は 3.16 を使っている。
大円の半径を r1 とおき，以下の方程式を解く。
@syms r1
円周率 = 3.16
r2 = r1 - 1
r3 = r1 - 2
eq1 = 円周率.*(r1^2 + r2^2 + r3^2) - 65.57;
# r1: 大円の半径
ans_r1 = solve(eq1, r1)[2] |> println # 2 of 2
3.50000000000000
直径はそれぞれ，7 寸, 5 寸, 3 寸である。
算額（その1543）
千葉県富津市寺尾 六所神社 明治4年(1871)
山口正義：やまぶき 第 57 号
https://yamabukiwasan.sakura.ne.jp/ymbk57.pdf
山口正義：やまぶき 第 59 号
https://yamabukiwasan.sakura.ne.jp/ymbk59.pdf
キーワード：球9個，3次元
#Julia, #SymPy, #算額, #和算
外球の中に甲球 1 個，乙球 7 個を容れる。外球の直径が 12 寸，甲球の直径が 8 寸のとき，乙球の直径はいかほどか。
外球の半径と中心座標を R, (0, 0, 0)
甲球の半径と中心座標を r1, (0, 0, R - r1)
乙球の半径と中心座標を r2, (x2, y2, z2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms R, r1, r2, x2, y2, z2
eq1 = x2^2 + (R - r1 - z2)^2 - (r1 + r2)^2
eq2 = x2^2 + z2^2 - (R - r2)^2
eq3 = x2*sind(Sym(360)/14) - r2
res = solve([eq1, eq2, eq3], (r2, x2, z2))[2] # 2 of 2
(2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2), 2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/((R^2 - 2*R*r1*cos(2*pi/7) + r1^2)*sin(pi/7)), R*(R^2 - 2*R*r1 - r1^2 + 2*r1^2*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2))
# 2r2
2res[1](R => 12/2, r1 => 8/2).evalf() |> println
3.27511575021090
外球，甲球の直径がそれぞれ 12 寸，8 寸のとき，乙球の直径は 3.27511575021090 である。
答・術は 4.7071 寸としているが，山口も 3.275 有奇として異を唱えている。
function draw(R, r1, more=false)
pyplot(size=(500, 250), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, x2, z2) = (2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2), 2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/((R^2 - 2*R*r1*cos(2*pi/7) + r1^2)*sin(pi/7)), R*(R^2 - 2*R*r1 - r1^2 + 2*r1^2*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2))
p1 = plot()
circle(0, 0, R, :cyan)
circle(0, 0, x2, :blue)
rotate(x2, 0, r2, angle=360/7)
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, "R", :black, :center, :bottom, delta=delta/2)
point(x2, 0, "乙球:r2,(x2,0,z2)", :red, :right, delta=-delta/2, deltax=15delta)
end
p2 = plot()
circle(0, 0, R, :cyan)
circle(0, R- r1, r1, :blue)
circle(x2, z2, r2)
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, "R", :black, :center, :bottom, delta=delta/2)
point(0, R - r1, "甲球:(0,0, R-r1)", :blue, :center, delta=-delta/2)
point(x2, z2, "乙球:r2,(x2,0,z2)", :red, :right, delta=-delta/2, deltax=15delta)
end
plot(p1, p2)
end;
draw(12/2, 8/2, true)
算額（その1542）
四十 埼玉県熊谷市玉井 玉井大神社 嘉永元年(1848)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
山口正義：やまぶき 第 4 号
https://yamabukiwasan.sakura.ne.jp/ymbk4.pdf
キーワード：正五角形，面積分割
#Julia, #SymPy, #算額, #和算
正五角形の一辺が 1 寸のとき，この面積を三等分するときの截面（一辺の分割した長さ；図参照）はいかほどか。
計算しなくても，図を見ればわかる。
四角形OAEB = 五角形OECDE' ならば，三角形OBE が三角形OCE の二倍，すなわち CE が BC の 1/3 である。
あえて計算すると，以下のようになる。
正五角形が内接する円の半径を r，一辺の長さを「五角面」，截面を「截面」とおき，以下の連立方程式を解く。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms r, x, y, 五角面
s18 = Sym(18)
s36 = Sym(36)
r = 五角面/2sind(s36)
S = r^2*sind(s36)*cosd(s36) # 正五角形の面積/5
s1 = area([0 0; r*cosd(s18) r*sind(s18); x y]) # 面積/5 の更に 2/3
s2 = area([0 0; x y; r*sind(s36) -r*cosd(s36)]) # 面積/5 の更に 1/3
eq1 = s1 - 2S/3 # S/5 + s1 = S/5 + s2
eq2 = s2 - S/3;
res = solve([eq1, eq2], (x, y));
@syms d
ans_x = apart(res[x]) |> simplify
ans_x |> println
五角面*(sqrt(5) + 5)/12
ans_y = apart(res[y]) |> simplify
ans_y|> println
-五角面*(5*sqrt(2) + 3*sqrt(10))*sqrt(sqrt(5) + 5)/120
截面 = sqrt((r*sind(s36) - ans_x)^2 + (-r*cosd(s36) - ans_y)^2) |> simplify
截面 |> println
截面(五角面 => 1).evalf() |> println
sqrt(10)*sqrt(-15*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) - 6*sqrt(5)*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 20*sqrt(5) + 110)*sqrt(五角面^2)/(30*sqrt(5 - sqrt(5)))
0.333333333333333
截面は五角面の 1/3 である。
function draw(五角面, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = 五角面/2sind(36)
x = 五角面*(√5 + 5)/12
y = -五角面*(5√2 + 3√10)*sqrt(√5 + 5)/120
(x1, y1) = r.*(sind(36), -cosd(36))
(x2, y2) = r.*(cosd(18), sind(18))
length = sqrt((x - x1)^2 + (y - y1)^2)
@printf("正五角形の一辺の長さが %g のとき，面積を 3 分割するときの\"截面\"は %g である。\n", 五角面, length)
plot([-x, 0, x], [y, 0, y], color=:blue, lw=0.5, showaxis=false)
polygon(0, 0, r, 5, color=:red)
segment(0, 0, 0, r, :blue)
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(x, y, " (x,y)", :blue, :left, :vcenter)
point(x1, y1, " (x1,y1)", :red, :left, :bottom, delta=delta/2)
dimension_line(-x1, y1 - delta, x1, y1 - delta, "五角面", :red, :center, delta=-2delta)
dimension_line(-x1 - 1.5delta, y1 - delta/2, -x - 1.5delta, y - delta/2, "截面", :magenta, :right)
segment(0, 0, -x2, y2, :gray80)
segment(0, 0, -x1, y1, :gray80)
segment(0, 0, x1, y1, :gray80)
point(0, 0, " O", :blue, :left, :bottom, delta=delta/2)
point(0, r, "A", :blue, :center, :bottom, delta=delta/2)
point(-x2, y2, "B ", :blue, :right, :vcenter)
point(-x1, y1, " C", :blue, :left, :bottom, delta=delta/2)
point(x1, y1, "D ", :blue, :right, :bottom, delta=delta/2)
point(-x, y, "E ", :blue, :right, :bottom, delta=delta/2)
point(x, y, "E' ", :blue, :right, delta=-delta/2)
end
end;
draw(1, true)
算額（その1541）
埼玉県本庄市都島 正観寺 戸塚盛政の算額 享保11年(1726)
山口正義：やまぶき 第 26 号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf
山口正義：やまぶき 第 40 号
https://yamabukiwasan.sakura.ne.jp/ymbk40.pdf
甲，乙，丙，丁，戊の立方体がある。甲と乙の体積の和が 189，丙と丁と戊の体積の和が 36 である。辺の長さの差は一定である。それぞれの一辺の長さはいかほどか。
立方体:甲の一辺の長さを「甲」，辺の差を「差」とおき，以下の連立方程式を解く。
問題文に出てくる 2 つの定数を記号のままとして解こうとすると，"RecursionError('maximum recursion depth exceeded')" が発生するので解けない。それほど複雑だということであろう。
using SymPy
@syms 甲::positive, 差::positive
eq1 = 甲^3 + (甲 - 差)^3 - 189
eq2 = (甲 -2差)^3 + (甲 -3差)^3 + (甲 -4差)^3 - 36
res = solve([eq1, eq2], (甲, 差))[1]
(5, 1)
甲の一辺の長さは 5，差が 1 なので，順次 4, 3, 2, 1 となる。
山口の解説は「（まともにやれば）9 次方程式が得られるが，これをとくのはやっかいである」と書いているが，SymPy では何の苦も無く解ける（当たり前だが）。
丙，丁，戊の体積の和が 36 なので，丙^3 + 丁^3 + 戊^3 = 36, 丙 > 丁 > 戊, 各辺は整数なので，3, 2, 1 が解であることはすぐにわかり，5^3 + 4^3 = 189 となるのを確認すれば解を得るのは簡単である。直感で解いていはいけないのかな？
直感で解かれるのを避けるなら，以下のようなテストデータを作ればよい。
甲^3 + 乙^3 = 3699397765
丙^3 + 丁^3 + 戊^3 = 5120681355
using SymPy
@syms 甲::positive, 差::positive
eq1 = 甲^3 + (甲 - 差)^3 - 3699397765
eq2 = (甲 -2差)^3 + (甲 -3差)^3 + (甲 -4差)^3 - 5120681355
res = solve([eq1, eq2], (甲, 差))[1]
(1234, 13)
甲の一辺の長さは 1234，差は 13 である。
1. ブルートフォースで解く
甲の一辺の長さと差が未知なので，探索範囲を誤ると答えが出ないが，実行時間を見ながら範囲を広げていけば解が見つかる可能性は高い。
実行時間を短くするために色々策を弄することもできるが，まずは，単純なプログラムでやってみる。
以下のプログラムは 2 秒以内で正解を出す。
@time begin
K1 = 3699397765
K2 = 5120681355
for 差 = 1:100
for 甲 = 5:999999
乙 = 甲 - 差
丙 = 甲 - 2差
丁 = 甲 - 3差
戊 = 甲 - 4差
戊 > 0 && 甲^3 + 乙^3 == K1 && 丙^3 + 丁^3 + 戊^3 == K2 && (println("甲 = $甲, 差 = $差"); break)
end
end
end
甲 = 1234, 差 = 13
1.396169 seconds (99.00 M allocations: 1.476 GiB, 5.47% gc time, 0.70% compilation time)
新年最初の満月，「Wolf Moon（ウルフムーン）」は，まるで大自然が秘める静かな力を感じさせるような神秘的な輝きを放っています。
冷たい冬の夜空に浮かぶ満月。その光は凍てついた地表をやさしく照らし，まるで狼の遠吠えを待ち望んでいるかのようです。写真の中の月の輪郭が鮮明に浮かび上がり，そのクレーターの影が物語るのは，地球からどれだけ遠い距離を旅しているかということ…
Wolf Moon の由来と伝説
Wolf Moon の名前は，ネイティブアメリカンの伝承に由来します。厳しい冬の間，食糧を求めて遠吠えする狼の姿を見て，人々はこの満月に特別な名前を与えたといわれています。そのため，この満月は「絆」や「忍耐」の象徴ともされています。
Wolf Moon を楽しむためのヒント
• 月明かりの下で散歩を楽しむ
• ホットドリンクを片手に，夜空を見上げる
• 詩や音楽でこの幻想的な瞬間を心に刻む
今夜，あなたもこの特別な満月を見上げてみませんか？
算額（その1540）
竹束（ちくそく）問題
源為憲：口遊（くちずさみ），天禄元年(970)
中村信弥：和算家 北明 寺島宗伴
http://www.wasan.jp/terasima/terasima.html
図のような竹の束がある。周員（外側の赤い竹の本数）が 21 本のとき，総本数はいかほどか。
SymPy では，以下の 2 つの関数を使う
実際の数列 a を作る
a = sympy.sequence(第 n 項の式, (n, 最初の項, 最後の項))
部分列の和の一般項
summation(a.formula, (n, 最初の項, 最後の項))
竹の本数は，3, 9, 15, ... の等差級数なので，「第n項の式」は 6n - 3
using SymPy
@syms n, m
b = sympy.sequence(6n - 3, (n, 1, 10))
[3, 9, 15, 21, ...]
1 から m までの部分列の和の式は以下で求まる。
Sm = summation(b.formula, (n, 1, m))
Sm |> println
3*m^2
外側の竹が 21 本は，何項目かは (21 + 初項)/公差で求まる。
(21 + 3)/6 |> println
4.0
Sm の m に代入すると，総和が求まる。
Sm(m => 4) |> println
48
まだるっこしいので，ブルートフォースで計算する方が簡単である。
target = 74067
s = 0
for n = 1:999999 # 適当な数
a = 6n - 3
s += a
a == target && (println(s); break)
end
457197075
算額（その1539）
（三）兵庫県 (5) 兵庫県伊丹市寺本堂山 昆陽寺 嘉永7年(1854)，昭和43年(1968)復元奉納
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：直角三角形，菱形
#Julia, #SymPy, #算額, #和算
直角三角形内に菱形を入れる。鈎，股（，弦）の長さが与えられたとき，菱形の一辺の長さはいかほどか。
算額（その924）で，鈎・股を入れ替えただけのものである。
鈎，股，弦の長さをそのまま変数名として使う。
菱形の下側の一辺が鈎，股と交差する座標を (0, b), (a, 0)
菱形の一辺の長さを diamond
とおき，以下の連立方程式を解く。
include("julia-source.txt");
using SymPy
@syms 鈎::positive, 股::positive, 弦::positive,
a::positive, b::positive, diamond::positive
eq1 = a^2 + b^2 - diamond^2
eq2 = diamond/(股 - a) - 鈎/股
eq3 = b/a - 鈎/股
(a, b, diamond) = solve([eq1, eq2, eq3], (a, b, diamond))[1];
# diamond: 菱形の一辺の長さ
diamond |> println
鈎*(股^2 + 鈎^2 - 鈎*sqrt(股^2 + 鈎^2))/股^2
菱形の一辺の長さは 鈎*(股^2 + 鈎^2 - 鈎*sqrt(股^2 + 鈎^2))/股^2 である。
鈎，股の長さがそれぞれ 3 寸，4 寸のとき，菱形の一辺の長さは 15/8 = 1.875 寸である。
diamond |> println
diamond(鈎 => 3, 股 => 4) |> println
diamond(鈎 => 3, 股 => 4).evalf() |> println
鈎*(股^2 + 鈎^2 - 鈎*sqrt(股^2 + 鈎^2))/股^2
15/8
1.87500000000000
# a
a |> println
a(鈎 => 3, 股 => 4) |> println
a(鈎 => 3, 股 => 4).evalf() |> println
鈎*(-鈎 + sqrt(股^2 + 鈎^2))/股
3/2
1.50000000000000
# b
b |> println
b(鈎 => 3, 股 => 4) |> println
b(鈎 => 3, 股 => 4).evalf() |> println
鈎^2*(-鈎 + sqrt(股^2 + 鈎^2))/股^2
9/8
1.12500000000000
考察：「問」において，鈎，股の他に冗長にも弦も与えたときとしているのは，術で「鈎*弦/(鈎 + 弦)」といいたいためである。
function draw(鈎, 股, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
弦 = sqrt(鈎^2 + 股^2)
(a, b, diamond) = (鈎*(-鈎 + sqrt(股^2 + 鈎^2))/股, 鈎^2*(-鈎 + sqrt(股^2 + 鈎^2))/股^2, 鈎*(股^2 + 鈎^2 - 鈎*sqrt(股^2 + 鈎^2))/股^2)
@printf("鈎 %g; 股 = %g; 弦 = %g のとき， 菱形の一辺の長さ = %g\n", 鈎, 股, 弦, diamond)
@printf("a = %g; b = %g; diamond = %g\n", a, b, diamond)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
plot!([0, a, a, 0, 0], [b, 0, diamond, 鈎, b], color=:red, lw=1)
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, " 股", :blue, :left, :bottom, delta=delta/2)
point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
point(0, 鈎, " 鈎", :blue, :left, :bottom, delta=delta/2)
point(0, b, " b", :blue, :left, :vcenter)
plot!(xlims=(-3delta, 股+5delta))
end
end;
draw(3, 4, true)
算額（その1538）
(6) 大阪府茨木市大字総持寺 総持寺 嘉永7年(1854)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円6個，外円
#Julia, #SymPy, #算額, #和算
問は「如図大円之内五箇円径問」とだけあり，答は「中円一尺二寸，脇円八寸，小円六寸」とある。
注：大円の直径は中円の直径の 2 倍ではあるが，問で大円の直径を述べて置かなければ問が成り立たない。また，算額の図がどうなっているかのわからないが横に 3 個並んでいる円は等円のようだし，脇円，小円なのだろうが，答えにある直径を加えると2尺2寸になり，大円の直径より小さくなる。とりあえず脇円，小円を区別して立式するが，脇円が中円，小円，大円に接しているということは譲れない。また，「只云大円径二尺四寸」があるものとする。
大円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1, (0, r1)
脇円の半径と中心座標を r2, (R - r2, 0)
小円の半径と中心座標を r3, (0, 0)
とおき，以下の方程式を解く。
include("julia-source.txt")
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms R, r1, r2, r3
eq1 = r1 - R/2
eq2 = r3 - (2r1 - 2r2)
eq3 = r1^2 + (r3 + r3)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))[1] # 1 of 2
(R/2, R/3, R/3)
中円の半径は大円の半径の 1/2 倍，脇円は同じく 1/3 倍，小円も 1/3 倍ということで，脇円と小円は同じ大きさであるということになった。
大円の直径が 2 尺 4 寸のとき，中円の直径は 1 尺 2 寸，脇円，小円の直径は 8 寸である。
function draw(R, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = R.*(1/2, 1/3, 1/3)
plot()
circle(0, 0, R)
circle22(0, r1, r1, :green)
circle2(R - r2, 0, r2, :blue)
circle(0, 0, r3, :magenta)
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, "R", :red, :center, :bottom, delta=delta/2)
point(0, r1, "中円:r1,(0,r1)", :green, :center, delta=-delta/2)
point(R - r2, 0, "脇円:r2,(R-r2,0)", :blue, :center, delta=-delta/2)
point(0, 0, "小円:r3,(0,0)", :magenta, :center, delta=-delta/2)
end
end;
draw(24/2, true)
