裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その1031)

2024年06月04日 | Julia

算額(その1031)

三十一 一関市舞川 観福寺内地蔵堂前額 明治34年(1901)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市
.

http://www.wasan.jp/yamamura/yamamura.html

2 個の大円が外接し,それぞれの内部に中円と 3 個の小円,更に外部に 2 個の小円を配置する。小円の直径が与えられたときに中円の直径を求めよ。

注:うっかり手を付けると条件不足ではないかと気づく。問をよく読めば「接大小円各五所」と書いてある。最初はそれぞれの大円が小円と 5 箇所の接点で接していると読んでいたのだが,共通接点だとすれば,接点は全体で 5 箇所になる。図では共通接点で接していないように見える端の 2 小円(下記で中心が (x3, y31), (0, y32) の2円)は大円と共通の接点(図の緑色の点)で接しているということだ。

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

include("julia-source.txt");

using SymPy
@syms r1::positive,
     r2::positive, r3::positive, x3::positive,
     y31::positive, y32::positive
eq1 = (r1 - x3)^2 + y31^2 - (r1 - r3)^2
eq2 = r1^2 + y32^2 - (r1 + r3)^2
eq3 = (2r1 - r2 - x3)^2 + y31^2 - (r2 + r3)^2
eq4 = (x3 - r3)^2 + y31^2 - 4r3^2
eq5 = y31/(r1 - x3) - y32/r1  # 端にある2個の小円と大円が一点で接する
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x3, y31, y32))[1]

   (r3*(2 + sqrt(5)), r3*(2*sqrt(5) + 5)/3, r3*(1 + sqrt(5))/2, r3*(-1 + sqrt(5))*sqrt(2*sqrt(5) + 5)/2, r3*sqrt(2*sqrt(5) + 5))

中円の半径 r2 は 小円の半径 r3 の (2√5 + 5)/3 倍である。
術の「二十個開平方加五個以三個除之乗小円径」は「20 の平方根に5 を加え3で割ったものを小円径に掛ける」であるから,術と一致する。
小円の直径が 1 寸のとき,中円の直径は 3.1573786516665265 寸である。

その他のパラメータは以下のとおりである。
r3 = 0.5;  r1 = 2.11803;  r2 = 1.57869;  x3 = 0.809017;  y31 = 0.951057;  y32 = 1.53884

(2√5 + 5)/3

   3.1573786516665265

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r3 = 1/2
   (r1, r2, x3, y31, y32) = (r3*(2 + sqrt(5)), r3*(2*sqrt(5) + 5)/3, r3*(1 + sqrt(5))/2, r3*(-1 + sqrt(5))*sqrt(2*sqrt(5) + 5)/2, r3*sqrt(2*sqrt(5) + 5))
   @printf("小円の直径が %g のとき,中円の直径は %g である。\n", 2r3, 2r2)
   @printf("r3 = %g;  r1 = %g;  r2 = %g;  x3 = %g;  y31 = %g;  y32 = %g\n", r3, r1, r2, x3, y31, y32)
   plot()
   circle2(r1, 0, r1, :blue)
   circle22(0, y32, r3)
   circle4(x3, y31, r3)
   circle2(r3, 0, r3)
   circle2(2r1 - r2, 0, 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(r1, 0, "大円:r1,(r1,0)", :blue, :center, delta=-delta/2)
       point(2r1 - r2, 0, "中円:r2,(2r1-r2,0)", :green, :center, :bottom, delta=delta/2)
       point(r3, 0, "小円:r3,(r3,0) ", :black, :right, :vcenter)
       point(x3, y31, "小円:r3,(x3,y31)", :black, :left, :bottom, delta=delta/2)
       point(0, y32, "小円:r3,(0,y32)", :black, :left, :bottom, delta=delta/2)
       point(x3/2, (y32 + y31)/2)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1030)

2024年06月04日 | Julia

算額(その1030)

三十 一関市山ノ目 配志和神社 嘉永5年(1852)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

外円の中に菱形を 2 個,等円を 2 個容れる。等円の直径が 1 寸のとき,外円の直径はいかほどか。

注:菱形についての記述がない。2 つの菱形の頂点が外円の中心になると仮定して解を求める。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x, y)
菱形の中心を (0, R/2), (0, -R/2)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive,
     r::positive, x::positive, y::negative
eq1 = x^2 + y^2 - (R - r)^2
eq2 = dist2(0, 0, √Sym(3)R/2, R/2, x, y, r)
eq3 = dist2(0, 0, √Sym(3)R/6, -R/2, x, y, r)
solve([eq1, eq2, eq3], (R, x, y))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r*(1 + sqrt(2)), r*(1 + sqrt(3))/2, r*(1/2 - sqrt(3)/2))

外円の半径 R は 等円の半径 r の (1 + √2) 倍である。
等円の直径が 1 寸のとき,外円の直径は 2.414213562373095 寸である。

答えは「大円径五寸一分九厘有奇」となっている。条件の差によるのだろうか?

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (R, x, y) = (r*(1 + sqrt(2)), r*(1 + sqrt(3))/2, r*(1/2 - sqrt(3)/2))
   @printf("小円の直径が %g のとき,外円の直径は %g\n", 2r, 2R)
   @printf("r = %g;  R = %g;  x = %g;  y = %g\n", r, R, x, y)
   plot(R .* [0, √3/2, 0, -√3/2, 0], R .* [0, 1/2, 1, 1/2, 0], color=:green, lw=0.5)
   plot!(R .* [0, √3/6, 0, -√3/6, 0], R .* [-1, -1/2, 0, -1/2, -1], color=:green, lw=0.5)
   circle(0, 0, R)
   circle2(x, y, 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, "小円:r,(x, y)", :red, :center, delta=-delta/2)
       point(R, 0, "R", :blue, :left, delta=-delta/2)
       point(√3R/2, R/2, " (√3R/2,R/2)  ", :green, :right, :vcenter)
       point(√3R/6, -R/2, "(√3R/6,-R/2) ", :green, :right, :vcenter)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1029)

2024年06月04日 | Julia

算額(その1029)

三十 一関市山ノ目 配志和神社 嘉永5年(1852)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

正方形の中に大円(半円)と小円を 2 個ずつ容れる。小円の直径が 1 寸のとき,大円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive,
     y2::positive
eq1 = r2^2 + (y2 - r1)^2 - (r1 - r2)^2
eq2 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
solve([eq1, eq2], (r1, y2))[1]

   (9*r2/4, 3*r2)

大円の半径 r1 は小円の半径 r2 の 9/4倍である。
小円の直径が 1 寸のとき,大円の直径は 2.25 寸である。

術では大円の直径は小円の直径の √392 + 3 = 22.79898987322333 倍としている。流石に,山村もこれはおかしいと「答」に合わすように √4.8 + 3 = 5.19 倍とした。
しかし,略図を描いてみるだけで,これもおかしいと気づかないものか。算額の図は往々にして与えられた条件に基づくものでないことがあるが,この場合は条件を変えようがない。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 1/2
   (r1, y2) = (9*r2/4, 3*r2)
   @printf("小円の直径が %g のとき,大円の直径は %g\n", 2r2, 2r1)
   @printf("r2 = %g;  r1 = %g;  y2 = %g\n", r1, r2, y2)
   plot([0, 2r1, 2r1, 0, 0], [0, 0, 2r1, 2r1, 0], color=:green, lw=0.5)
   circle(2r1 - y2, r2, r2)
   circle(2r1 - r2, y2, r2)
   circle(r1, 0, r1, :blue, beginangle=0, endangle=180)
   circle(2r1, r1, r1, :blue, beginangle=90, endangle=270)
   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(2r1 - y2, r2, "小円:r2,(2r1-y2,r2)", :red, :center, delta=-delta/2)
       point(2r1 - r2, y2, "小円:r2,(2r1-r2,y2)", :red, :center, delta=-delta/2)
       point(r1, 0, "r1", :blue, :center, delta=-delta/2)
       point(2r1, 0, "2r1", :blue, :center, delta=-delta/2)
       point(2r1, 2r1, "(2r1,2r1)", :blue, :right, :bottom, delta=delta/2)
       point(2r1, r1, "", :blue)
       plot!(ylims=(-4delta, 2r1+3delta))
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1028)

2024年06月04日 | Julia

算額(その1028)

二十三 一関市萩荘字八幡 達小袋八幡神社 弘化3年(1840)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

大円に中心を通る線を 4 本引き,等円を 8 個描く。等円の直径が与えられたとき,大円の直径を求めよ。

大円の半径と中心座標を R, (0, 0)
小円の半径と中心座標を r, ((R - r)/√2, (R - r)/√2), (R + r, 0)
斜線の端点座標を (x0, y0)
斜線の勾配を θ°
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r::positive,
     θ::positive, x0::positive, y0::positive
@syms R, r, θ, x0, y0
x0 = R + 2r
θ = asind(y0/sqrt(x0^2 + y0^2))
eq1 = (R - r)*sind(Sym(45) - θ) - (R + r)*sind(θ)
eq2 = dist2(0, 0, x0, y0, (R - r)/√Sym(2), (R - r)/√Sym(2), r);

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)
   (R, y0) = u
   return [
       -y0*(R + r)/sqrt(y0^2 + (R + 2*r)^2) + (R - r)*sin(pi*(-asin(y0/sqrt(y0^2 + (R + 2*r)^2))/pi + 1/4)),  # eq1
       R^4 + 2*R^3*r - 2*R^3*y0 - 5*R^2*r^2 + R^2*y0^2 - 12*R*r^3 + 6*R*r^2*y0 - 2*R*r*y0^2 - 4*r^4 - 4*r^3*y0 - r^2*y0^2,  # eq2
   ]
end;
r = 1/2
iniv = BigFloat[2, 1]
res = nls(H, ini=iniv)

   ([1.4872641552444894, 0.6466017927282794], true)

小円の直径が 1 のとき,大円の直径は 2*1.4872641552444894 = 2.974528310488979 である。

山村の術の解説は途中で複数の間違いがあり,最終的な答えも間違っている(ほかにも多くの誤りがあるが,校正などしていないのだろうか)。正しくは,以下のとおり。
小円径 = 1
天 = sqrt(2) + 3
大円径 = 小円径 * (sqrt(sqrt(4天 + 2) + 天))

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = 1/2
   (R, y0) = res[1]
   x0 = R + 2r
   θ = asind(y0/sqrt(x0^2 + y0^2))
   println("小円の直径が $(2r) のとき,大円の直径 = $(2R)")
   @printf("r = %g;  R = %g;  θ = %g;  x0 = %g;  y0 = %g\n", r, R, θ, x0, y0)
   plot()
   circle(0, 0, R, :blue)
   circle4((R - r)/√2, (R - r)/√2, r, :magenta)
   circle42(0, (R + r), r, :magenta)
   segment(x0, y0, -x0, -y0)
   segment(-x0, y0, x0, -y0)
   segment(y0, x0, -y0, -x0)
   segment(-y0, x0, y0, -x0)
   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(x0, y0, "(x0,y0) ", :green, :right, :bottom, delta=delta/2)
       circle(0, 0, 0.4R, beginangle=0, endangle=θ)
       circle(0, 0, 0.42R, beginangle=0, endangle=θ)
       point(0.4R, 0, " θ", :red, :left, :bottom, delta=delta/2, mark=false)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(R + r, 0, "等円:r\n(R+r,0)", :magenta, :center, delta=-delta/2)
       point((R - r)/√2, (R - r)/√2, "等円:r\n((R-r)/√2,(R-r)/√2)", :magenta, :center, :bottom, delta=delta/2)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1027)

2024年06月04日 | Julia

算額(その1027)

二十三 一関市萩荘字八幡 達小袋八幡神社 弘化3年(1840)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

大円の中に中円,乙円,小円を 3 個ずつと中央に甲円を 1 個入れる。中円の直径が与えられたとき,大円の直径を求めよ。

注:「問」では「小円及び甲乙円各三個容れる」とあるが,図では「甲」はなく,「小」が 4 個描かれている。中央にあるのが甲で(当然1個しかない),それ以外の小円が「小」である。また,図では中円は接していないように描かれているが,実際には接している。

大円の半径と中心座標を R, (0, 0)
中円の半径と中心座標を r1, (x1, y1); x1 = (R - r1)*√3/2; y1= (R - r1)/2
甲円の半径と中心座標を R - 2r1, (0, 0)
乙円の半径と中心座標を r2, (0, R - r2)
小円の半径と中心座標を r3, (0, R - 2r2 - r3)
とおき,以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive,
     r2::positive, x2::positive, y2::positive,
     r3::positive, y3::positive
x1 = (R - r1)*√Sym(3)/2
y1 = (R - r1)/2
eq1 = x1 - r1
eq2 = x1^2 + (R - r2 - y1)^2 - (r1 + r2)^2
eq3 = x1^2 + (R - 2r2 - r3 - y1)^2 - (r1 + r3)^2
solve([eq1, eq2, eq3], (R, r2, r3))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r1*(3 + 2*sqrt(3))/3, 4*sqrt(3)*r1/33 + 3*r1/11, 6*sqrt(3)*r1/253 + 19*r1/253)

大円の半径は,中円の半径の (3 + 2√3)/3 倍である。
中円の直径が 181 のとき,大円の半径は 390.0007974466445 である。

その他のパラメータは以下のとおりである。
r1 = 90.5;  R = 195;  r2 = 43.6819;  r3 = 10.5138

181*(3 + 2√3)/3

   390.0007974466445

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 181/2
   (R, r2, r3) = r1 .* ((3 + 2√3)/3, (4√3 + 9)/33, (6√3 + 19)/253)
   x1 = (R - r1)*√3/2
   y1 = (R - r1)/2
   @printf("中円の直径が %g のとき,大円の直径は %g である。\n", 2r1, 2R)
   @printf("r1 = %g;  R = %g;  r2 = %g;  r3 = %g\n", r1, R, r2, r3)
   plot()
   circle(0, 0, R, :blue)
   circle(0, 0, R - 2r1, :green)
   rotate(0, r1 - R, r1)
   rotate(0, R - r2, r2, :orange)
   rotate(0, R - 2r2 - r3, r3, :purple)
   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, y1, "中円:r1,(x1,y1)", :red, :center, delta=-delta/2)
       point(0, R - 2r2 - r3, "  小円:r3,(0,R-2r2-r3)", :purple, :left, :vcenter)
       point(0, 0, "   甲円:R-2r1,(0,0)", :green, :left, :vcenter)
       point(0, R - r2, "  乙円:r2,(0,R-r2)", :orange, :left, :vcenter)
       point(0, R, " R", :blue, :left, :bottom, delta=delta/2)
   end
end;

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

(2 * 3 * 5 * 7 * 11 * 13)**10 は何桁の数か?

2024年06月04日 | ブログラミング

こちらにも投稿しておこう。

受験生は log10(2), log10(3) は知っているとの前提で

import math

もともとの数値

(2 * 3 * 5 * 7 * 11 * 13)**10

    596421543032951827553785377955490490000000000

それは知らない(わからない)ことにしておいて,

10 乗する数は

2 * 3 * 5 * 7 * 11 * 13

    30030

まず,30030 よりほんの少し小さい数を 2, 3, 5 のべき乗の積で作る

2 * 4 * 4 * 5 * 10 * 16, 2**8 * 100

    (25600, 25600)

25600**10 の常用対数を取る

10*(8*math.log10(2) + 2)

    44.0823996531185

次に,30030 よりほんの少し大きい数を 2, 3, 5 のべき乗の積で作る

2 * 4 * 4 * 6 * 10 * 16, 2**10 * 3 * 10

    (30720, 30720)

30720**10 の常用対数を取る

10*(10*math.log10(2) + math.log10(3) + 1)

    44.87421211359475

3 つの数値は,同じ桁数になったようだ...

(2 * 3 * 5 * 7 * 11 * 13)**10 の常用対数は 25600 の常用対数より大きく, 30720 の常用対数より小さい

10*(8*math.log10(2) + 2) < math.log10((2 * 3 * 5 * 7 * 11 * 13)**10) < math.log10((2**10 * 3 * 10)**10)

    True

さらに,桁数は同じ 45 桁である。

10*(8*math.log10(2) + 2), math.log10((2 * 3 * 5 * 7 * 11 * 13)**10), math.log10((2**10 * 3 * 10)**10)

    (44.0823996531185, 44.77555332198981, 44.874212113594744)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村