2011年8月28日日曜日

多面体描画道楽

Archimedes多面体13種のうち, rhomb(斜方とか菱形)という形容詞のつくものが2つある. 斜方立方八面体と斜方二十十二面体である. 何ゆえに斜方といわれるのか.





これらの図に示すように, この2種類の多面体には, それぞれ3種類の面がある.

斜方立方八面体では, 青は元々の立方体の面, 緑は正八面体の面である. 斜方二十十二面体では, 青は元々の正十二面体の面, 緑は正二十面体の面である.

では, それぞれの赤の面はなにか.

赤の面を延長して, 緑の面の上方での交点を決め, 赤の正方形に外接する菱形を考えることが出来る. 赤は元々そういう面の多面体の面である.

この多面体がどちらも菱形になるので, 斜方といわれる所以である. つまり, 斜方(赤)立方(青)八(緑)面体, 斜方(赤)二十(緑)十二(青)面体という命名であったわけだ.

ではその斜方多面体を描いてみよう.


上は斜方立方八面体である. x, y, zの対称軸も描いてある. z軸(青)のまわりに少し回転し, y軸(緑)のまわりにも僅かに回転したようになっている.


そのz軸, y軸まわりの回転をやめ, x軸方向から眺めたのが, この上の図だ. そこで, 最初の赤い正方形を含む面を考えてみると, この図の赤線で示す菱形が得られる. その座標が分かれば, 次の図が描けるわけだ.



これは, 対角線の比が1:√2の菱形が12個で構成されているので, 斜方十二面体(rhombic dodecahedron)という.


斜方二十十二面体の方は, 面が多いせいか, 多少手ごわい.



これが元の斜方二十十二面体で, 上と同様にz軸, y軸について回転してある.



回転角を0とし, x方向から眺めると, このように見えるはずだ. 赤線はそこにかぶせた菱形である. この方は, 菱形30枚で構成され, 斜方三十面体(triacontahedron)という.



菱形の対角線の比は1:φ(黄金比)である.

同一の菱形だけで出来る多面体はこの2つだけらしい. またKeplerはこの2つの多面体のあることを知っていたといわれる.

2011年8月27日土曜日

3シリンダ機関車

ずいぶん前のこのブログに3シリンダ機関車のことを書いた. そこで引用しているURLが, 最近なくなっているのに気づいたので, もう一度3シリンダ機関車の弁装置について書くことにする.

今でも蒸気機関車は人気があり, 各地で復活されているが, その構造はあまり理解されていないに違いない.

簡単にいってしまえば, 機関車の左右の前方につけた2個のシリンダ(気筒)の中にピストンを入れ, ピストンの両側から交互に蒸気を補給してピストンを前後に動かす. ピストンの往復運動を主連棒で繋ぎ, 動輪を回転して機関車を進める. ピストンがどちらかの端に寄っていると, 蒸気の補給が出来ず, 死点といって力が出せないが, 反対側のピストンが最も力の出る中央にいるように, 左右のピストンの位相を90度づらせて配置してある. このことは, 鉄道博物館で実物の機関車を見ると確認出来る.

ピストンの前後に交互に蒸気の補給をしているのが, シリンダと並んだ蒸気室で, その中に滑り弁があり, これが蒸気室の中を前後に滑りながら, 蒸気をピストンの一方に送り込む.

Walschaertの走り装置の図はここにあるので参照されたい.

滑り弁は, 動輪に動力を伝える主連棒と90度づれたエキセン棒から, 逆転装置を経由した弁棒で押し引きされて運動する. この辺の装置は大変に込み入っているから, 説明は省略!

ところで, 表題の3シリンダ機関車は, 2つのシリンダでは力不足の時, 第3のシリンダを両シリンダの中央に搭載し, 合わせて3本の主連棒で動輪を駆動する方式である. 動輪軸の中央はクランク状になっていて, そこに中央の主連棒が接続されている.

問題は, 左右の動輪のエキセン棒のようなものが, 中央のクランクからも出ているかということだが, そういう複雑な走り装置にはなっておらず, 中央のシリンダの弁の制御は両側の弁の制御から構成していたというのが前回のブログであった.

下の図を見て欲しい. Cylinder AとCylinder Bと書いたのが, 両端のシリンダの蒸気室で, 左(機関車の後方)から弁棒が入り, 途中に滑り弁があり, 弁棒の延長が右方に突き出している. 機関車の先頭にAOD, BDCのような梃子があり, 支点AはCylinder Aの弁棒の先に, 支点BはCylinder Bの弁棒の先に固定されている. 梃子だから, 図で左右に揺れると支点の上下移動もあるはずだが, 今は各支点は左右にだけ動くとしている.

支点Oは固定点である. またOD=OA/2, つまりAが動くとDはAの半分の振幅で逆位相に動く. そのDを支点として, DB=DCの梃子があり, CがCylinder Cの弁棒になっている. そうすると, CはAとBが120度の位相差で動くと, さらに120度の位相差で動き, 3番目の滑り弁が望み通りに動くのである.

このアニメーションがここにある.

アニメの下でくるくる回る3芒星は, 120度の位相で回転することを示す.

2011年8月26日金曜日

ベクトルの和

「とっておきの数学パズル」に, 多面体の各面の, 直交外向きで, 面積に比例する大きさのベクトルの和が0になることを証明せよというのがあった. Q5.13

多分そうなりそうだが, 証明となると....

「いかにして問題を解くか」には, 「一様な四面体の重心を求めよ」という問題が解けなければ, 平面の類似の問題「一様な三角形の重心を求めよ」を考えてみるという教訓がある.

その伝でいくなら, 多角形の各辺の, 直交外向きで, 辺の長さに比例する大きさのベクトルの和が0になることを証明せよとなる.

たとえば, 下の図の左は, 最も簡単な多角形の三角形の場合で, 辺の長さをa,b,cとする. (5:3:4になっている.) 従って, 外向きのベクトルは, 黒い矢印のように出来る. ベクトルの和を計算するのだから, aの先端までbの付け根を移動し(青矢印), その先端までcの付け根を移動する(青矢印)と, 3つのベクトルで構成される三角形は, もとの三角形と相似になり, cの先端はaの付け根に一致し, 和が0なことが分かる.



しかし, こういう方法はどうだろうか. (図の右)

元の三角形の左右方向をx軸. 上下方向をy軸とする. 各ベクトルのx, y成分を考える. bはx成分しか, cはy成分しか持たない. aは, 図にax, ayと書いて示したx成分とy成分がある. ベクトルaの大きさが, 斜面aの長さに比例するなら, axはaの斜めの面の, x方向からみた断面(青線で示す)の長さに比例するわけで, この断面の長さはbの辺の長さに等しく, 従ってaxの大きさはbの大きさを同じで方向が反対である.従って, ベクトルのx成分の和は0.

同様にして, y方向の成分の和も0. 従って, ベクトル全体の和も0になる.

これを3次元の問題でも使ってみよう. ある面の, その面積に比例した直交外向きベクトルaのx成分は, ベクトルとx軸のなす角をθとすると, a cosθ, この面のx方向から見た断面積もa cosθのはずである.



したがって, この多面体の直交外向きのベクトルのx+方向の和は, x+方向から見た断面積になり, x-方向の和は同じ断面積で反対向きになり, x方向全体は打ち消す.

y,z方向も同じわけだから, 結局和は0というので証明出来たのではないか.

ところで, ある時気づいた, 多面体を水の中に置くという考えも棄て難い. ある面に立てた, 面積に比例する直交外向きのベクトルは, その面に働く水圧に相当する. 水中の多面体が流れないのは, その水圧の総和が0だからだと思うと, この証明は終りらしいが.

2011年8月12日金曜日

3での整除

TAOCPに, レジスタ上の64ビットの数が3で整除出来るかどうかの判定プログラムをMMIXで書けという問題7.1.3-215がある.

十進の場合, 各桁の和が3で整除出来るなら, 元の数も3で整除出来るのは知られている.

理由は1も10も100も..., つまり10nが3を法として1であるからだ.

二進の場合も同様に考えられる. 二進の1, 10, 100,...は, 交互に3を法として1と-1だから, ビットに右から0,1,2,...と番号をつければ, 偶数番目の1の数と奇数番目の1の数の差が3の倍数の時に3で整除出来る.

これは十進法で11で整除出来る判定と同じだ. 3003や1023や5060は11で整除出来る.

だから, (nu x)をxの横方向の和をとる関数, (band x y)をxとyとのビットごとのandをとる関数とし, 32ビットでお許し頂くと

(= (modulo (- (nu (band x #x55555555)) (nu (band x #xaaaaaaaa))) 3) 0)

ということになる. 最後にまた3での剰余をとるのは違反っぽいが, 十進の場合の各桁の和の3の整除と同じで, 大きい数をいきなり割るのよりは簡単であり, とりあえずはこういう考えでよかろう.

この辺でTAOCPの解答をみると, 一見妙なプログラムである. MMIXのアセンブリコードに慣れていない向きには申し訳ないが, コメントをつけてある.



最後の3で割った剰余をみるのに, 001001...を左シフトし, 符号ビットが1になるのを利用している. これはうまい.

ではその前は何をしたか. レジスタ0は偶数ビットと奇数ビットの横方向の和(これをr0とする). レジスタ1は偶数ビットの横方向の和(r1とする). ところで, 奇数ビットの横方向の和(r2とする)も欲しい. r2=r0-r1のはず. レジスタ0はその後32-r0になり, 2ADDUでレジスタ1の2倍に足される.

つまりレジスタ1はr1+r1+(32-r0)=32+r1-r2になり, #x249...をこれだけ左シフトすることになる. bの定数は32ビット目, 8文字目の次が#x9なので, r1=r2のとき符号ビットが1, 負数になり, 3で整除出来ることが分かる.

巧妙なプログラムであった. 斯くの如きプログラムの解読は楽しい.

2011年8月11日木曜日

多面体描画道楽

立方体は正方形6枚で出来ている. 正方形は1つの内角が90度. 従って正方形全体で内角の和は360度. それがF=6枚あるから, 2160度, つまり12π(パイ)になる.



一方頂点の数Vは8. それに対して内角の総和は(2V-4)π=12πという関係をEulerは知っていたらしい.

別の例. 正十二面体は正五角形がF=12枚. 頂点はHamilton世界巡礼でお馴染のV=20. 正五角形の内角は108度. 従って内角の総和は, 108*5*12=6480. 6480/180=36(π). 2V-4も36だ.

さらに2種類の面のある, 切頭二十面体(サッカーボール)は, 正五角形が12, 正六角形が20, 頂点は60. 内角の総和は, 108*5*12+120*6*20=20880, 20880/180=116. 2*60-4=116となる.

こう考えてみた.

Eulerの式に, 面の数F, 辺の数E, 頂点の数Vとすると, F+V=E+2という有名なのがある.

いま, ある多面体がn1角形, n2角形,...,nF角形のFの面で出来ていたとする.

n角形の内角の和は(n-2)πであったから, 計算したい内角の総和Sは
S=(n1-2)π+(n2-2)π+...+(nF-2)π=(n1+n2+...+nF)π-2Fπ.

一方, この多面体の辺の数を計算すると, 各多角形の辺は2回ずつ寄与しているから, E=(n1+n2+...+nF)/2.

従って, n1+n2+...+nF=2E. これをSの式に入れると,

S=2Eπ-2Fπ. Eulerの式から, E-F=V-2.

ゆえにS=2(V-2)π.

なるほど.

また, Descartesの定理というのもある.



左下は, 正四面体の頂点は, 正三角形が3個集まっていて, その内角の和が, 360度より足りない度合い(deficiency)が180度であることを示す. 正四面体には頂点が4個あるから,deficiencyの総和は720度. 4πである.

これを他の正多面体でやって見ると, やはり, どれも720度になるのである.

サッカーボールではどうか.

この図のように, 120度, 120度, 108度に挟まれた角は12度. それが60個あるから, やはり720度である.




Eulerの方は, サッカーボールの図でいうと, 348度の方を総和している. 60倍すれば先ほどの20880度になる.

V個のradianで表わしたdeficiencyをd1, d2,...,dVとする.

Eulerの方の内角の和は, (2π-d1)+(2π-d2)+...(2π-dV)
=2πV-(d1+d2+...+dV) と足す.

これが2(V-2)πだったから, d1+d2+...+dV=2πV-2(V-2)π=4π.

これもなるほどであった.

正多面体では, 面の数の順に正四面体, 立方体(正六面体), 正八面体, 正十二面体, 正二十面体と並べる. 正四面体の頂点は尖っているから, あとに行くほど丸くなるかと思うのは誤解であって, 正八面体も正二十面体も結構尖っている. それは正八面体は立方体より, 正二十面体は正十二面体よりも, deficiencyが大きく, 頂点の数が少ないためであろう.

サッカーボールがよく転がるのは, deficiencyが12度で, V=60もあるからである.

2011年8月4日木曜日

3つの円

「とっておきの数学パズル」に, 「3つの円」という問題があった. 5.14



平面上で離れている, 半径の違う2つの円の, 円の間で交わらない(外側の)共通接線の交点を, それらの円の焦点という. 3つの円があると, 焦点は3個出来るわけだが, その3点は一直線上にあることを証明せよというのである.

こういう図で考えてみる.


円O0は中心の座標をx0,y0, 半径をr0, 円O1は座標をx1,y1, 半径をr1とする. また共通接線の1本と, 円との接点をT0, T1, 共通接線の交点P01の座標を, x01, y01とする.

三角形O0P01T0と三角形O1P01T1とは相似で, 大きさの比はr0 : r1である. 従って
(x01-x0) : (x01-x1)=r0 : r1,   (y01-y0) : (y01-y1)=r0 : r1.

これから, x01=(r0x1-r1x0)/(r0-r1),   y01=(r0y1-r1y0)/(r0-r1).

x12, y12, x20, y20についても, 同様な式が得られる.

ここで, 3点(x0,y0), (x1,y1), (x2,y2)が一直線上にあることの検査は, このブログの2008年4月24日の「同一直線上の3点」にあるように, 3点で作る三角形の面積が0かどうかで判定出来る.

今回もP01, P12, P20の三角形の面積を求めれば良い.

まず座標は
x01=(r0x1-r1x0)/(r0-r1)
y01=(r0y1-r1y0)/(r0-r1)
x12=(r1x2-r2x1)/(r1-r2)
y12=(r1y2-r2y1)/(r1-r2)
x20=(r2x0-r0x2)/(r2-r0)
y20=(r2y0-r0y2)/(r2-r0)
だから, 三角形の面積の2倍は
x01y12+x12y20+x20y01-x12y01-x20y12-x01y20

通分して分母を払い, 計算を続けると

  (r0x1-r1x0)(r1y2-r2y1)(r2-r0)+(r1x2-r2x1)(r2y0-r0y2)(r0-r1)
+(r2x0-r0x2)(r0y1-r1y0)(r1-r2)-(r1x2-r2x1)(r0y1-r1y0)(r2-r0)
-(r2x0-r0x2)(r1y2-r2y1)(r0-r1)-(r0x1-r1x0)(r2y0-r0y2)(r1-r2)

展開すると,
 (r0r1x1y2-r0r2x1y1-r1r1x0y2+r1r2x0y1)(r2-r0)
+(r1r2x2y0-r1r0x2y2-r2r2x1y0+r2r0x1y2)(r0-r1)
+(r2r0x0y1-r2r1x0y0-r0r0x2y1+r0r1x2y0)(r1-r2)
-(r1r0x2y1-r1r1x2y0-r2r0x1y1+r2r1x1y0)(r2-r0)
-(r2r1x0y2-r2r2x0y1-r0r1x2y2+r0r2x2y1)(r0-r1)
-(r0r2x1y0-r0r0x1y2-r1r2x0y0+r1r0x0y2)(r1-r2)

さらに展開する.
+r0r1r2x1y2-r0r2r2x1y1-r1r1r2x0y2+r1r2r2x0y1
-r0r0r1x1y2+r0r0r2x1y1+r0r1r1x0y2-r0r1r2x0y1
+r0r1r2x2y0-r0r0r1x2y2-r0r2r2x1y0+r0r0r2x1y2
-r1r1r2x2y0+r0r1r1x2y2+r1r2r2x1y0-r0r1r2x1y2
+r0r1r2x0y1-r1r1r2x0y0-r0r0r1x2y1+r0r1r1x2y0
-r0r2r2x0y1+r1r2r2x0y0+r0r0r2x2y1-r0r1r2x2y0
-r0r1r2x2y1+r1r1r2x2y0+r0r2r2x1y1-r1r2r2x1y0
+r0r0r1x2y1-r0r1r1x2y0-r0r0r2x1y1+r0r1r2x1y0
-r0r1r2x0y2+r0r2r2x0y1+r0r0r1x2y2-r0r0r2x2y1
+r1r1r2x0y2-r1r2r2x0y1-r0r1r1x2y2+r0r1r2x2y1
-r0r1r2x1y0+r0r0r1x1y2+r1r1r2x0y0-r0r1r1x0y2
+r0r2r2x1y0-r0r0r2x1y2-r1r2r2x0y0+r0r1r2x0y2

同じrの項をまとめると,
r0r0r1(-x1y2-x2y2-x2y1+x2y1+x2y2+x1y2)
r0r0r2(+x1y1+x1y2+x2y1-x1y1-x2y1-x1y2)
r0r1r1(+x0y2+x2y2+x2y0-x2y0-x2y2-x0y2)
r0r2r2(-x1y1-x1y0-x0y1+x1y1+x0y1+x1y0)
r1r1r2(-x0y2-x2y0-x0y0+x2y0+x0y2+x0y0)
r1r2r2(+x0y1+x1y0+x0y0-x1y0-x0y1-x0y0)
r0r1r2(+x1y2-x0y1+x2y0-x1y2+x0y1-x2y0
  -x2y1+x1y0-x0y2+x2y1-x1y0+x0y2)

となり, かっこ内はプラスマイナス打ち消し, めでたく0になった.

実際に上のように書いて計算したわけではない. 変数の順を覚えておき, 引数だけを使って計算した. 最後に, ブログ用に変数名を追加するには, utilispのプログラムを書いて実行した.

ところで, 数式処理システムrisa/asirを使ってみたら,

X01=(r0*x1-r1*x0)/(r0-r1);
X12=(r1*x2-r2*x1)/(r1-r2);
X20=(r2*x0-r0*x2)/(r2-r0);

Y01=(r0*y1-r1*y0)/(r0-r1);
Y12=(r1*y2-r2*y1)/(r1-r2);
Y20=(r2*y0-r0*y2)/(r2-r0);

X01*Y12+X12*Y20+X20*Y01-X12*Y01-X20*Y12-X01*Y20;

を入力すると, 瞬時に0が返り, 安堵する.

さて, 例によって, 幾何学的に証明できないかと考える. 幾何の証明で大切なのは, 図をなるべく正確に描くことだ. 幸い, このブログの最初の図を描いたPostScriptのプログラムがあるから, それを多少修正すると以下の図が得られる.

今度は3つの円O0, O1, O2の中心をA, B, Cとしている. また, 焦点をE, F, Gとした. BからEFに平行線を引き, AFとの交点をDとする.

円の中心O0, O1と焦点P01の距離の比は, 半径r0, r1の比であった.

EはA,Bの焦点だから, BE/AE=r1/r0
FはA,Cの焦点だから, CF/AF=r2/r0

BD//EFにしたから, DF/AF=r1/r0
従って, CF/DF==r2/r1

また, GがB,Cの焦点だから, CF/DR=r2/r1=CG/GB
従って, FG//BD//EF

QED.

この方が楽だなぁ.

2011年8月2日火曜日

条件付きの倍数

三角関数表, 対数表をはじめとして, 数表の存在感の薄い昨今である. なにしろ手元の計算機で関数値が一発で計算出来る世の中だからだ. だが今でも, 理科年表には, 10.0から100.9までの四桁の対数表があり, なんとなく空しい.

計算機が手元にない, 古き良き時代には, 数表は重要であった. Charles Babbageが階差機関を考案したのも, 数表を自動作成したいからであった.

江戸時代に歩いて日本地図を作成した伊能忠敬に関連した展示会にいってみると, 漢数字で書かれた各種の数表も並べられていて, そんな昔から数表が使われていたことに驚く.

ところで対数表の話だ. 高木先生の解析概論には気になる脚注がある. 「Wolframノ表ニハ1000以下ノ素数ノ自然対数ノ50桁ノ表ガ掲ゲラレテヰル. コノ表ハ既ニ少年がうすガ愛用シタモノデアル. 」だ.

なるほど, ある数の対数は, その素因数の対数の和なので, 素数の対数表だけあれば充分なのであった. しかし50桁も足すのは大変であったろう. また常用対数に変換するには, ln 10で割らなければならず, これも一仕事であるが, 昔の人は計算は得意であったからだとまずは納得する.

常用対数表は, スコットランド人のJohn Napierが1614年に作成, 出版したといわれている. そこで, 2014年に, 出版400年を記念してエディンバラで会合を開こうという話が進んでいる.

Napierの後, 何人もが常用対数表を計算し, 刊行したが, 先頃, Edward Sangの対数表の話をウェブで見つけた. Sangは1875年に, 素数について28桁の常用対数表を完成させた. その計算方法は次のようだ.

小さい素数の方から順に計算していく. さて, 素数1619の常用対数を計算するには, その1619を何倍かし, 下の方の桁を...999999のようにする. 9は何桁でもよいが, 多い方がよい. 例えば

714021*1619=1155999999

を作る. この計算が, このブログのタイトルの, ...9999で終るという「条件付きの倍数」である.

この倍数の作り方は次の通り.

まず1619の1から9までの倍数表を用意する.

(map (lambda (x) (display (list x (* x 1619))) (newline))
(a2b 1 10))
(1 1619)
(2 3238)
(3 4857)
(4 6476)
(5 8095)
(6 9714)
(7 11333)
(8 12952)
(9 14571)

被乗数が2と5以外の素数だから, 最後の桁には1から9が出揃う.

次にこの表を見ながら, 下の桁から積の各桁が9になるように1619の倍数を足していく. まず1倍の1619が9で終るから, 1619と書く. 下から2桁目は, 1だから, 最後が8で終る倍数, 3238をそこに足す. すると33999になる. 4桁目は3だから, 最後が6で終る倍数, 6476を足す. 今度は6509になる. 0のところを9にしなければならにから, 9で終る1619をもう一度足す. 2269になる. 6のところを9にするので, 11333を足す. こうしてここまでで1155999999が得られる.

1619 1
3238 2
----------
33999
6476 40
----
6509
1619 1
-----
2269
11333 7
-----
11559

この数を素因数分解すると, 1155999999=3×7×11×11×281×1619 である. いま常用対数を求めようとしている1619以外の素因数は, すべてこれより小さいから, それらの常用対数は分かっているはずだ.

従って, 1155999999の常用対数が得られれば, 1619の常用対数は,
log 1619=log 1155999999 - log 3 - log 7 - log 11 - log 11 - log 281
で得られるわけである.

ところで, log 1155999999=log(1156×106-1).
1156=17×17×2×2だから, この常用対数は既知である.

一方, x=1/(1156×106)とおいて,

log(1156×106-1)=log(1156×106×(1-x))=log(1156×106)+log(1-x).

自然対数をlnと書くと, ln(1+x)のTaylor展開は

ln(1+x)=x-x2/2+x3/3-x4/4+...    (|x|≤1, x≠-1)

だったから,

log(1-x)=ln(1-x)/ln 10=(-x-x2/2-x3/3-x4/4-...)/ln 10

Schemeで計算してみよう. Schemeのlog関数は自然対数lnのことである.

(define x (/ 1. (* 1156 1000000)))

x=>8.650519031141868e-10

(/ (* x x) 2)=>3.741573975407383e-19

(/ (* x x x) 3)=>2.157770458712447e-28

(log 10)=>2.302585092994046

常用対数のlog 1156×106

(/ (log (* 1156 (expt 10 6))) (log 10))=>9.06295783408451

従って log 11559999として

(- 9.06295783408451
(/ (+ x (/ (* x x) 2) (/ (* x x x) 3)) 2.302585092994046))
=> 9.062957833708824

が得られる. そこでlog 1619は

(- 9.062957833708824
(/ (+ (log 3) (log 7) (log 11) (log 11) (log 281)) (log 10)))
=>3.2092468487533754

log 1619をSchemeで直接計算すると

(/ (log 1619) (log 10)) => 3.2092468487533736

最後の2桁が違うがまぁよかろう.

こんな数値計算をするもの久しぶりだ.