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