算額（その317）

算額（その317）

p02-03 特集「算額」（杉浦）.indd
https://www.s-coop.net/lifestage/backnumber/2011/pdf/1106_02-03.pdf
御香宮の算額
甲，乙，丙の 3 つの正方形は，以下の条件を満たす。
1. 各面積の和はある一定の値を取る
2. 甲の一辺の 3 乗根，乙の一辺の 5 乗根，丙の一辺の 7 乗根の和はある一定の値を取る
3. 甲と乙の一辺の長さの差と，乙と丙の一辺の長さの差は等しい
このとき，甲，乙，丙それぞれの辺の長さを求めよ。

御香宮神社の算額
http://tadahikostar.blog21.fc2.com/blog-entry-3847.html
には，鮮明な写真が掲載されている。

読んだところでは，
1. 甲乙丙の各平積の和が既知，
2. 甲の開立方，乙の開四乗方，丙の開六乗方の和が既知
3. 甲，乙，丙の一辺の長さは等差数列（甲>乙>丙）
で，条件 2. が違う。

また，
御香宮（京都・伏見）の算額の問題 - 高精度計算サイト - Keisan
https://keisan.casio.jp/exec/user/1465639202
では，
2. 甲の1辺の3乗と乙の1辺の5乗と丙の1辺の7乗がある値をとり
と書かれている。これは明らかに誤りである。
「開立」は立方根を表す。三乗は「再自乗」という。

以下では，「甲の開立方，乙の開四乗方，丙の開六乗方の和」を採用する。

既知である値を記号（変数）として解くのは solve() でも無理である。
また，既知の値を具体的に与えても solve() では解けないので，nlsolve() で数値解を求める。

具体的に，既知の値を求めておく。甲，乙，丙の辺の長さをそれぞれ 6, 4, 2 として，条件 1, 2 で表される数値を計算する。

(甲, 乙, 丙) = (6, 4, 2)
甲^2 + 乙^2 + 丙^2 |> println  # a とする
甲^(1/3) + 乙^(1/4) + 丙^(1/6) |> println  # b とする

   56
   4.353796203514608

   丙^2 + 乙^2 + 甲^2 - 56,  # eq1
   丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461,  # eq2
   丙 - 2*乙 + 甲,  # eq3

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)
   (甲, 乙, 丙) = u
   return [
       丙^2 + 乙^2 + 甲^2 - 56,  # eq1
       丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461,  # eq2
       丙 - 2*乙 + 甲,  # eq3
   ]
end;

初期値により，二通りの解が示される。両方ともに，条件を満たしている。

iniv = [big"7.0", 4.1, 1.1]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[5.999999999999978421814505360998991236330686683731818877458530801081713663484236, 4.000000000000010789092747335389143058164701865210416514719281049866373768067708, 2.000000000000043156370989309779294879998717046689014151980031298601986612507766], true)

(甲, 乙, 丙) = res[1]  # 6, 4, 2
丙^2 + 乙^2 + 甲^2 - 56,  # eq1
丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461,  # eq2
丙 - 2*乙 + 甲  # eq3

   (2.566627137897993161354744300647618625279842621939623293996006529315262445625456e-25, -3.867992033135278266328028707026632666657885527869404966426262217141461630526278e-27, -4.904726014337897321140092172924707876804539997277609351100695166069929253914544e-65)

iniv = [big"3.0", 2.1, 1.2]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[4.417532066579428278706241618180857546089192311716748510275041555425975958939494, 4.319756156028245623996284797567623698490811039184724945910874439881535311042085, 4.221980245477062969286327976954389850892429766652701381546707324388269228906724], true)

(甲, 乙, 丙) = res[1]  # 4.417532066579429, 4.319756156028245, 4.221980245477063
丙^2 + 乙^2 + 甲^2 - 56,  # eq1
丙^0.166666666666667 + 乙^0.25 + 甲^0.333333333333333 - 4.35379620351461,  # eq2
丙 - 2*乙 + 甲  # eq3

   (7.928474500472361481718927810447226817906506437873289558404530451020621884955364e-25, -5.071109598580144750966960369570845052151000656128305908528664840899463931562674e-27, 5.117456576204827560281226639245275222567489206961363157351246079614929159341669e-65)