2.675 を小数点以下3桁目を四捨五入せよ。
2.68 となることを予想しているのだが，
★ Python
>>> round(2.675, 2)
2.67
で，decimal パッケージの Decimal を使いましょう...ってさ
★ R
> round(2.675, 2)
[1] 2.67
★ Julia
julia> round(2.675, digits=2)
2.68
★ AWK には round 関数はないのだが int(2.675*100+0.5)/100 というような関数を定義すればよい。
$ awk 'BEGIN{print int(2.675*100+0.5)/100}'
2.68
まあ，完全かどうかは虱潰しにでも調べるしかないか?
算額（その24）
大きな円の中に，半径の異なる 4 種類の円が，図のように配置されている。各円の半径を求めよ。
岩手県東磐井郡藤沢町 藤勢寺 元治2年
http://www.wasan.jp/iwate/fujise.html
外側の円の半径を 2 とし，内部の円の中心座標と半径を小さい順に (0, y1, r1), (0, y2, r2), (x3, y3, r3), (1, 0, 1) と定める。
using SymPy
@syms r1::positive, r2::positive, r3::positive,
y1::positive, y2::positive, x3::positive, y3::positive;
未知数は r1, r2, r3, x3, y1, y2, y3 の 7 個で，以下の 7 つの方程式を立てて解く。
eq1 = x3^2 + (y1 - y3)^2 - (r1 + r3)^2;
eq2 = x3^2 + (y3 - y2)^2 - (r2 + r3)^2;
eq3 = 1 + y2^2 - (1 + r2)^2 |> expand
eq4 = (1 - x3)^2 + y3^2 - (1 + r3)^2;
eq5 = y1 - y2 - r1 - r2;
eq6 = y1 - 2 + r1;
eq7 = x3^2 + y3^2 - (2 - r3)^2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r2, r3, x3, y1, y2, y3));
解は以下のごとくである。
name = ["r1", "r2", "r3", "x3", "y1", "y2", "y3"]
for i = 1:7
println("$(name[i]) = $(res[1][i]) = $(res[1][i].evalf())")
end
r1 = 8/11 - 2*sqrt(5)/11 = 0.320714913181856
r2 = 16*sqrt(5)/209 + 46/209 = 0.391277931291850
r3 = 7/11 - sqrt(5)/11 = 0.433084729318201
x3 = 1/11 + 3*sqrt(5)/11 = 0.700745812045397
y1 = 2*sqrt(5)/11 + 14/11 = 1.67928508681814
y2 = 68/209 + 60*sqrt(5)/209 = 0.967292242344437
y3 = 2/11 + 6*sqrt(5)/11 = 1.40149162409079
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=:top)
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")
(r1, r2, r3, x3, y1, y2, y3) = (8/11 - 2*sqrt(5)/11, 16*sqrt(5)/209 + 46/209, 7/11 - sqrt(5)/11, 1/11 + 3*sqrt(5)/11, 2*sqrt(5)/11 + 14/11, 68/209 + 60*sqrt(5)/209, 2/11 + 6*sqrt(5)/11)
println("r1 = $r1, r2 = $r2, r3 = $r3")
println("y1 = $y1, y2 = $y2, y3 = $y3")
plot()
circle(0, y1, r1, :green)
circle(0, -y1, r1, :green)
circle(0, y2, r2, :blue)
circle(0, -y2, r2, :blue)
circle( x3, y3, r3, :magenta)
circle(-x3, y3, r3, :magenta)
circle( x3, -y3, r3, :magenta)
circle(-x3, -y3, r3, :magenta)
circle( 1, 0, 1)
circle(-1, 0, 1)
if more
point(0, y1, "(0,y1,r1)", :brown, :center)
point(0, y2, "(0,y2,r2)", :blue, :center)
point(x3, y3, "(x3,y3,r3)", :magenta, :center)
point(1, 0, "(1,0,1)", :magenta, :center)
vline!([0], color=:black, linewidth=0.25)
hline!([0], color=:black, linewidth=0.5)
end
circle(0, 0, 2, :black)
end;
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], ftol=1e-14)
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
v = r.zero
end
return v, r.f_converged
end;
function h(u)
r1, r2, r3, x3, y1, y2, y3 = u
return [
x3^2 + (y1 - y3)^2 - (r1 + r3)^2,
x3^2 + (y3 - y2)^2 - (r2 + r3)^2,
1^2 + y2^2 - (1 +r2)^2,
(1 - x3)^2 + y3^2 - (1 + r3)^2,
(y1 - y2) - (r1+ r2),
y1 - 2 + r1,
x3^2 + y3^2 - (2 - r3)^2
]
end;
iniv = [0.354644124017785, 0.33761490117551013, 0.4639474199979425, 0.6081577400061725, 1.3714562740209328, 0.8883769604436443, 1.4105324143047882] # 初期値
res = nls(h, ini=iniv);
name = ["r1", "r2", "r3", "x3", "y1", "y2", "y3"]
for i = 1:7
println("$(name[i]) = $(res[1][i])")
end
r1 = 0.32071491318185646
r2 = 0.3912779312918499
r3 = 0.4330847293182009
x3 = 0.7007458120453972
y1 = 1.6792850868181435
y2 = 0.9672922423444372
y3 = 1.4014916240907944
