Heronの式

前回のブログ「四面体の体積」に名称だけ出てきたHeronの式というのがある. 三角形の3辺の長さをそれぞれa, b, cとし, (a+b+c)/2をSとすると, その三角形の 面積は√S(S-a)(S-b)(S-c)というのである.

たとえば, 辺の長さが3,4,5のPythagoras三角形ではS=6だから, 平方根号の中は 6×3×2×1=36で, 面積は6だ. 右の図, 1辺が2の正三角形では, S=3で, 根号内は3×1×1×1=3で, 面積は√3である.

Heronの式を覚えた時, これが正しいという証明を見たかどうかは分らない. 勿論, Sの次元は長さ, S-a,S-b,S-cも長さなので, 根号内の次元は長さの4乗; 開平すると 2乗になり, つまり面積だ. そういえば, 前回のブログにあった, 四面体の体積を6本の 辺の長さから求める式も体積の次元になっている. あたりまえだ.

でも定量的に合っているかを調べるのが今回の頭の体操である. 方針としては, 根号の あるのがウンシャンだから, 面積の二乗の値を考えることにする. 例題の三角形は こういう頂点名と辺名を持つ. 計算の途中で時々数値的にチェックするのにもこれを 使う. この三角形は∠Aが直角なので, 面積は√8×√18/2=6だ.

式を使わざるを得ないから, Texの出力のpdfにも助けてもらおう.

式(0)はHeronの式で, 2乗すると(1). Sを定義で戻し, 展開したのが式(2)である. この程度 の計算はWolframAlphaにやって貰う.

三角形ABCの各頂点の対辺をa,b,cとすると, その各々の2乗もpdfにある通り. これを 式(2)の分子に入れて展開するわけだが, これは中々面倒である. WolframAlphaも 時間がかかるとか文句をいうし, 我が家のRapsberry PiのMathematicaは今は ネットに繋っていないので, 使い勝手が悪い. 結局は多項式の乗算のプログラムを Schemeでちょこっと作り, それで展開した. 時々例題の値を入れて験算する.

そして展開した結果は項数21のpdfにある式である. これも目がくらくらするから, x_ix_jy_ky_lについて, ijと klの対の項の係数を表にしてみた. さらに例題の座標で験算すると576. 16で 割ると36だ. しめしめ.

さてこの式は本当に三角形の面積であろうか. 頂点の座標(x₀y₀), (x₁y₁),(x₂y₂)の三角形の 面積(頂点の回り方によってはその負数)の式は, 幸いにも私は高校生の頃から知っている. 次のpdfにある式だ. 行列式を計算式にしてその分子の2乗を展開すると, またもや項数21の多項式が得られた. 今回もその係数の表を作るとブラボー! 先の表の1/4の値である. という訳けで Heronの式は面積になるのであった.

だが待てよ. この行列式はなぜ正しいのか. こちらも気になったので, 計算してみた.

図の三角形の, 座標軸に平行な線での外接長方形を作る. それから長方形の辺と 三角形の辺で作る3個の三角形の面積を引く. そういう式を作ってみる.

WoflramAlpha様に

Expand[(x2-x1)(y0-y2)-(x0-x1)(y0-y1)/2-(x2-x1)(y1-y2)/2-(x2-x0)(y0-y2)/2]

と伺いを立てたら

(x0y1)/2-(x0y2)/2-(x1y0)/2+(x1y2)/2+(x2 y0)/2-(x2y1)/2

のお告げが出た. 整理すると

(x0y1+x1y2+x2y0-x0y2-x1y0-x2y1)/2

となり, 行列式の値と一致している. もう一度ブラボー! ちゃんとした証明でもないが, 私にとっては十分納得出来る頭の体操であった.

四面体の体積

筑波大の三谷さんのツィッターに四面体の体積を計算する話があった.

角封筒の上の辺の中央と下の辺の両端を結ぶ折目をつけ, それを山折りにすると, 合同な 二等辺三角形4枚による四面体が出来る. 元の封筒の横と縦が1と√2の時の 体積を計算せよという問題である. ちょっとやってみたところ, 存外簡単であった.

下の図で, 赤い線は元の封筒を展開したもの. 下の辺と上の辺は繋っている. この斜めの 線を谷折りにすると, ABCを底面, ACD'を手前の面, ABDを向うの面, 辺DEとD'E'は 繋っていて斜面だが上の面になる.

錐体の体積を計算するには, 底面積と高さが必要だが, 底面積は簡単だ. 辺BCが1, 辺OAは辺ABが√2だから, √7/2である.

高さについては, 三角形ABDは, 辺ABを折目として頂点Dが起き上がるから, Dに対応する 平面上の点はABと直交する破線で示す線に沿って下がり, AO上のGまで来ると下から上って来たD’と 出会って頂点になる. GAの距離はADの距離にtan αを掛けたものだが, tan αはOB/OAだしADは1だから, GA=1/√7である. よって 高さは√6/√7.

求める体積は (√7/2)*(√6)/((√7)*2*3)=1/2√6

三谷さんのツィッターにはヒントと称して, 四面体の6辺の長さからその体積を 求める式があった. 平面三角形のHeronの公式の三次元版である. そこにある式の 添字の付け方を, 私流に修正したのが次の図である.

このように添字を付けると, 四面体の体積Vの2乗は

の一番上の式で得られる. しかしこれも目がくらくらするので, 始めの3項と終りの4項を それぞれ6変数のtfと3変数のtgに分け, 2行目以下のように定義する. 最初のV²を計算する式を書いたLatexは次の様だ. 計算式と同様に分解してある.

\newcommand{\sq}[1]{l_#1^2}
\newcommand{\tf}[6]{\sq#1\sq#2(\sq#3+
\sq#4+\sq#5+\sq#6-\sq#1-\sq#2)}
\newcommand{\tg}[3]{-\sq#1\sq#2\sq#3}
$V^2=\frac{1}{144}[\tf012345+\tf234501\\
+\tf450123\tg025\tg034\tg124\tg135]$

これを使ってpythonで計算したのが次だ.

import math
def sq(x):
  return x*x
def tf(a,b,c,d,e,f):
  return sq(a)*sq(b)
    *(sq(c)+sq(d)+sq(e)+sq(f)-sq(a)-sq(b))
def tg(a,b,c):
  return -sq(a)*sq(b)*sq(c)
def v2(l0,l1,l2,l3,l4,l5):
  return (tf(l0,l1,l2,l3,l4,l5)+tf(l2,l3,l4,l5,l0,l1)
    +tf(l4,l5,l0,l1,l2,l3)
    +tg(l0,l2,l5)+tg(l0,l3,l4)+tg(l1,l2,l4)
    +tg(l1,l3,l5))/144
l0=1;l1=1;l2=math.sqrt(2);l3=math.sqrt(2)
l4=math.sqrt(2);l5=math.sqrt(2)
print(math.sqrt(v2(l0,l1,l2,l3,l4,l5)))

=> 0.20412414523193143

この値は確かに1/2√6 である.

ところで, 最初の封筒の四面体の体積は, 形状が特殊であったせいか, 簡単であった. 底面の3点の座標と, その各点から頂点までの距離が与えられている時にも体積は 得られるだろうか. 以下はその計算である. 式が沢山あるので, Latexでpdfに したもので示す. 底辺の3点, A, B, Cの座標をそれぞれ(x₀,y₀,0), (x₁,y₁,0),(x₂,y₂,0); 頂点Dの座標を(x,y,z)とする. A,B,Cから頂点までの距離をそれぞれ l₀,l₁,l₂とする.

今回もpythonで計算すると

x0=0;y0=-1/2;x1=math.sqrt(7)/2;y1=0;x2=0;y2=1/2
l0=math.sqrt(2);l1=1;l2=math.sqrt(2)
a0=2*(x0-x1); b0=2*(y0-y1)
c0=sq(l1)-sq(l0)+sq(x0)-sq(x1)+sq(y0)-sq(y1)
a1=2*(x1-x2); b1=2*(y1-y2)
c1=sq(l2)-sq(l1)+sq(x1)-sq(x2)+sq(y1)-sq(y2)

x=(c0*b1-c1*b0)/(a0*b1-a1*b0)
y=(a0*c1-a1*c0)/(a0*b1-a1*b0)

z=math.sqrt(sq(l0)-sq(x0-x)-sq(y0-y))
print(x,y,z)

#=> 0.9449111825230684 0.0 0.9258200997725514

期待通りの値が得られた. でもこんな実数の値より, 1/2√6 の方が嬉しい.

Fixed Day Number

前回のブログのアルゴリズムに, 年初から前月末までの総日数を計算する式:

(367 * month - 362) // 12

があった. monthを1から12まで変えて値を計算すると.

[0, 31, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336]

が得られる. 演算子「//」で小数点以下を落しているが, その落された値を調べると, ( (4777 - 367 * month) // 12 の図も右に示す.)

のような図になる. どちらも2月は30日あるとする.

まず左の図に注目する. 1月から2月は下がり, 2月から3月は上がり, 4月へは下がり, ... 8月と9月へは下がり, ... つまり大の月の後は下がり, 小の月の後は上がるらしい.

計算して見ると下がる量は0.4167, 上がる量は0.5833である.

もう一度式を見ると, monthが1増えると, 367 / 12 つまり 30.5833 増える. これと30及び31の差であった. これで納得できるが, Calendrical Calculation にはこういう説明があった.

(7 * m - 2) // 12 + 30 * (m - 1)

この第2項 30 * (m - 1) はすべて小の月として前月末までの総日数.

第1項はm = 1, 2, ..., 12について,

0, 1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6
  31 30 31 30 31 30 31 31 30 31 30

である(上の行), 下の行は1月から11月までの日数で, こうして見ると, 大の月の 上で上の行の値は1増える. つまり小の月とした時の誤差を得ている.

そこで上の第1項と第2項を足すと, (367 * m - 362) // 12になるのであった.

2月を30日にしたという発想がよかったのだ. 一方, 年末から計算する方の式はこう導く.

欲しい値は1月, 2月, ..., 12月で, 2月も30日とするから

[367, 336, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]

である. 30の倍数の列

[360, 330, 300, 270, 240, 210, 180, 150, 120, 90, 60, 30]

に誤差の列

  [7, 6, 6, 5, 5, 4, 4, 3, 2, 2, 1, 1]
  

を足したものだ. これは7から年初からの計算の誤差を引いたのもだ. つまり,
7-(7 * m - 2) // 12
でこれを(13 - m) * 30 に足す. これを一纏めにしたい.

まず下の計算の1行目のように, 床関数を天井関数にする. そして1を足して床関数に戻す. しかしmの値によっては, 整数になることがあり, 1を足すと失敗する. 従って, 2行目のように, 1の代りに 11/12を足してみる. 3行目はそれに30日の倍数を足したところ. 最後に分数の外に ある部分を分子に足す. すると(4777-367m)//12が得られる.

どうだろうか.

Fixed Day Number

私の好きな本に, Edward M. ReingoldとNachum DershowitzのCalendrical Caluculation がある. 随分昔に購入したもので, 手元にあるのは2004年のThird Printing.

要するに暦に関するアルゴリズムが満載で面白い. この本では暦の変換の基準は, Fixed Day Numberといい, Gregorian暦が過去に遡って使われたとして, その1年1月1日(の正子)を 第1日とするものである. いちいちfixed day numberといわず, R.D.と いうらしい. (Rata Die)

似たようのものに, -4712年1月1日正午から起算するユリウス日があるが, これは途中でJulian暦 からGregolian暦になるので変換が面倒である.

Calendrical Caluculationにある, Gregorian暦の年月日からそのR.Dを計算する アルゴリズムを見てみよう. 本書のアルゴリズムはCommon Lispで記述されているが, 今日はPythonに翻訳する.

def fixeddaynumber(year, month, day):
    return (365 * (year - 1) + (year - 1) // 4
            - ((year - 1) // 100) + (year - 1) // 400
            + (367 * month - 362) // 12
            + (0 if (month <= 2) else
               (-1 if leap_year (year) else -2))
            + day)
     

最初の365 * (year - 1)は前年までがすべて平年であったとしての, 前年の 最後の日までの総日数である. しかし, 4年に一度閏年があるから, その分を足す. + (year - 1) // 4. Gregorian暦では, 100で整除される年は平年に するのであったからそれを引く. - ((year - 1) // 100). しかし400で 整除されるなら閏年にするからそれを足す. + (year - 1) // 400. 前年末までの総日数はこれで計算出来た.

次は前月末までの日数を計算する. 目標としては, 1月なら0. 2月なら31, 3月は 平年なら59, 閏年なら60, ... 閏年の補正は後でやることにし, 本書にあるアルゴリズム は(367 * month - 362) // 12というものだ. monthを 1から12まで変えてこの値を計算すると:

print(list(map (lambda m:(367*m-362)//12, range(1,13))))
[0, 31, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336]
    

3月からは, 平年なら2, 閏年なら1多いことがわかる. 従ってこの関数の month番目を使い, monthが1, 2ならこの値そのまま, 3以降はこの値から 平年は2, 閏年は1を引いて使う.

これで前年末までと前月末までの総日数が分ったから, 後はdayを足せばよい.

閏年の補正に本書ではこういう式を使う.

def leap_year (year):
    return (year % 4) == 0 & ((year % 400)
                              in [100, 200, 300])
    

これで1年1月1日を計算すると:

print (fixeddaynumber (1, 1, 1))
1
 

なるほど. また, おととし, 開業150年と話題になった鉄道開通の日, 1872年10月14日(明治5年には日本は まだ旧暦を使っていたが, これは太陽暦に換算してある)のR.D.を計算してみる.

各々の項の値は682915, 467, -18, 4, 275, -1, 14で, R.D.は683656 である.

ところで私は何度も暦の計算のプログラムを書いたが, その年の終までの日数を計算し, その月の前までの日数を引いて日を足すのが好きだ. 上のアルゴリズムに year - 1が4回もあるのが気に食わないのである.

それで普通は年末から今月0日までの日数を引くことになるが, 大体はそこは定数の 表を利用する. しかし, Calendricalの本にあるような式が使えないかと 考えた. Calendricalの式と同様, 1,2月は特別扱いにしてもよいとして, 定数など を変えてテストした. 欲しい値は:

[[365, 366], [334, 335], 306, 275, 245, 214,
 184, 153, 122, 92, 61, 31]
    

である. 左端が1月と2月で, かっこ内は左が平年, 右が閏年. 右端が12月. その結果, (4777 - 367 * m) // 12がよいらしかった. つまり

print(list(map (lambda m:(4777 - 367 * m) // 12,
                range(1,13))))
[367, 336, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]
    

1月, 2月は平年は2を引き, 閏年は1を引く.

ここで使う定数4777は非常にクリティカルで, 4776や 4778に変えてみると, 前者は7月が小さく, 後者は2月が大きい.

print(list(map (lambda m:(4776 - 367 * m) // 12,
                range(1,13))))
[367, 336, 306, 275, 245, 214, 183, 153, 122, 92, 61, 31]

print(list(map (lambda m:(4778 - 367 * m) // 12, 
                range(1,13))))
[367, 337, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]

これを利用して書いたfixed day number関数は次の通り:

def fixeddaynumber2(year, month, day):
    return (365 * year + year // 4 - (year // 100)
            + year // 400 - (4777 - 367 * month) // 12
            + (0 if (month > 2) else
               (-1 (year % 400 == 0)
                if (year % 100 == 0) else
                (year % 4 == 0) -2))
            + day)

こちらで鉄道記念日を計算したのと, 本年6月30日の値の7で割った剰余とを見ると:

print (fixeddaynumber2 (1872, 10, 14))
683656

print ((fixeddaynumber2 (2024, 6, 30)) % 7)
0       

だから, R.D.を7で割った剰余が曜日になるわけだ. つまり, R.D.1の日, 1年1月1日は月曜であった.

年末から月始めまでの日数を今回のように式で計算すると, 実引数に0月とか１３月が与え られても計算してしまう恐れがある. 表なら範囲外としてエラーに出来る. それを 心配したとしても, 日の値に範囲外が与えられる可能性はあるのだから, まぁ我慢することに しよう.