BASIC 《 プログラミング言語 《 大系 《 彷徨の神殿
『スピログラフ』

第3回 トロコイド

Part 3. Trochoid
2006.12.25

スピログラフという名称はHasbro社の商標で、実は一般には「内トロコイド」あるいは「外トロコイド」という名称になります。

内トロコイド・外トロコイドを描きつつ、その性質についても考えていきます。

  1. スピログラフ
    1. 内トロコイドの式
    2. スピログラフを描く
    3. 外トロコイドも描ける
      1. 試しに描く
      2. 式で確認
    4. 経過を省略して描く
  2. スピログラフの性質
    1. 軌跡が原点を通る条件
    2. 周回数の法則
    3. 周回数の実装
    4. 経過を省略して描く
  3. おわりに

1. スピログラフ

スピログラフ（spirograph）、一般的に言えば内トロコイド（hypotrochoid）は、内サイクロイドをちょっと変えただけのものです。

ちなみに、spiro-は「らせん」とかの意味だそうです。

trochoidの語源に関しては、力学で出てくる「滑車」が「trochlea」で、その辺が関係しているようです。

1.1. 内トロコイドの式

図1*1. スピログラフ

図1のように、スピログラフは内サイクロイドとほとんど同じです。

r1、r2、θ、φは同じ。

ただ違うのは、点Pが動円の円周上ではなく、動円の中心Cからr3のところにあるということです。

つまり、第2回の式(12)のr2がr3に変わるだけです。

そうすると、第2回の式(10)は次のようになります。

(1)
CP=(r3 cos(φ-θ),-r3 sin(φ-θ))=(r3 cos((r1-r2)/r2)θ,-r3 sin((r1-r2)/r2)θ)

第2回の式(8)と式(1)から、点Pの座標は、

(2)
P=C+CP=((r1-r2)cosθ+r3 cos((r1-r2)/r2)θ,(r1-r2)sinθ-r3 sin((r1-r2)/r2)θ)

と表せます。

第2回の式(11)との微妙な違いを見てください。

r3=r2にすれば、内サイクロイドになります。

1.2. スピログラフを描く

式(2)を元に、第2回のリスト7をベースにして、スピログラフを描きます。

リスト1. spirograph0.bas
!
! スピログラフを描く
!      Copyright(c) カイン Cain 2006
!
DECLARE EXTERNAL PICTURE circle
!
! 定数設定
!  半径
LET r1 = 7 ! 定円半径
LET r2 = 3 ! 動円半径
LET r3 = 2 ! 動円内描画点の半径
LET prompt$="半径r1,r2,r3を入力（デフォルト値 " &
&            & STR$(r1) & "," &STR$(r2) & "," & STR$(r3) & ")"
WHEN EXCEPTION IN
   INPUT PROMPT prompt$ : r1,r2,r3
USE
END WHEN
LET ar1 = ABS( r1 )
LET ar2 = ABS( r2 )
LET ar3 = ABS( r3 )
!  周回数
…
LET rounds = ABS( rounds )
!
LET win_size = ar1    !定円
IF r1 * r2 < 0 THEN   ! r1,r2のどちらか一方が負の場合は外トロコイド
   LET win_size = win_size + ar2 * 2
ELSEIF ar2 > ar1 THEN ! 内トロコイドだが、動円の方が大きい場合
   LET win_size = win_size + ( ar2 - ar1 ) * 2
END IF
IF ar3 > ar2 THEN     ! 動円の外の軌跡を描く場合
   LET win_size = win_size + ( ar3 - ar2 )
END IF
SET WINDOW -(win_size+1), win_size+1, -(win_size+1), win_size+1
DRAW grid
!
! 定円描画
DRAW circle( 0, 0, r1 )
!
! スピログラフ描画
! 関数定義
DEF xp(t) = ( r1 - r2 ) * COS( t ) + r3 * COS( ( r1 - r2 ) / r2 * t )
DEF yp(t) = ( r1 - r2 ) * SIN( t ) - r3 * SIN( ( r1 - r2 ) / r2 * t )
!
LET div = 100 ! 360°を分割する数
LET round_time = 1 ! t=0～2π を描くのにウェイトをかける秒数
…

リスト1では、主に第2回のリスト7からの変更箇所を示しています。

リスト1を実行すると、次のような図形が描けます。

IF ～ ELSEIF ～ ELSE ～ END IF 文
複数の条件で条件付けをする時は、 ELSEIF を使います。 
IF 条件1 THEN
   条件1が真の時に実行されるプログラム
ELSEIF 条件2 THEN
   条件1が偽で条件2が真の時に実行されるプログラム
ELSE
   条件1も条件2も偽の時に実行されるプログラム
END IF
という形になります。
更に ELSEIF を連ねて条件を増やすこともできます。
絶対値 ar1 、 ar2 、 ar3
これまでは r1 、 r2 、 r3 の絶対値をそのまま r1 、 r2 、 r3 に取っていましたが、リスト1では、別の変数 ar1 、 ar2 、 ar3 にとっています。
これは、 r1 、 r2 、 r3 が負の値の場合でも、その値をそのまま使って図形を描くためです。
具体的な話は次項で書きます。
ウィンドウサイズ
リスト1では、ウィンドウサイズを次のように計算しています。 
LET win_size = ar1    !定円
まず、初めに定円の大きさを基準とします。 
IF r1 * r2 < 0 THEN   ! r1,r2のどちらか一方が負の場合は外トロコイド
   LET win_size = win_size + ar2 * 2
この「負の場合」については、次項で。 
ELSEIF ar2 > ar1 THEN ! 内トロコイドだが、動円の方が大きい場合
   LET win_size = win_size + ( ar2 - ar1 ) * 2
END IF
この「動円の方が大きい場合」は、前回の内サイクロイドでやりました。 
IF ar3 > ar2 THEN     ! 動円の外の軌跡を描く場合
   LET win_size = win_size + ( ar3 - ar2 )
END IF
|r3|>|r2|の場合は、点Pが動円をはみ出すので、その分を追加しています。
例えば、 r1,r2,r3 を 3,1,2 としてスピログラフを描くと、次のようになります。
実際の定規ならこういう描き方は難しいのですが、これはプログラムなので、こういったこともできます。

やっとスピログラフを描くところまで辿り着きました。

1.3. 外トロコイドも描ける

リスト1で外トロコイド（epitrochoid）を描くこともできます。

1.3.1. 試しに描く

ここまでは、r1、r2が正の場合のみ考えてきました。

では、r1、r2の片方を負の値にしてみたら、どうなるでしょうか。

リスト1を、 r1,r2,r3 を 3,-1,1 と入力して実行してみましょう。

どこかで見た感じ。r1,r2 を 3,1 とした外サイクロイドと同じ形です。

このように、r1とr2の正負を異なるようにすると、外トロコイドが描けます。

一石二鳥、面白い性質ですね。

リスト1の

IF r1 * r2 < 0 THEN   ! r1,r2のどちらか一方が負の場合は外トロコイド
   LET win_size = win_size + ar2 * 2

の部分は、この場合を考慮して、動円の直径を加えています。

条件の r1 * r2 < 0 は、掛け合わせて正になるか負になるかで、 r1 と r2 の符号が同じか異なるかを判定しています。

1.3.2. 式で確認

このことを式で確認してみます。

式(2)にr2=-r2'を代入します（r2'>0）。

(3)
P=((r1-r2)cosθ+r3 cos((r1-r2)/r2)θ,(r1-r2)sinθ-r3 sin((r1-r2)/r2)θ)=((r1+r2')cosθ+r3 cos((r1+r2')/(-r2'))θ,(r1+r2)sinθ-r3 sin((r1+r2')/(-r2'))θ)=((r1+r2')cosθ+r3 cos((r1+r2')/r2')θ,(r1+r2)sinθ+r3 sin((r1+r2')/r2')θ)
図2*2. 外トロコイド

第2回の式(7)と似た形の式が得られました。

二項目の符号が逆になっているのは、点Pが、動円内で180度反対にあることを意味します。

図2のような感じです。

なので、リスト1で描いた外サイクロイドは、前回のプログラムで描いた外サイクロイドを少し回転させたような描かれ方になっています。

r2だけでなくr3も負にしてしまえば、同じになります。

1.4. 経過を省略して描く

描かれた結果だけを素早く見ることができるように、描画経過のアニメーションを省いたプログラムも作っておきます。

リスト2. spirograph1.bas
…
LET div = 100 ! 360°を分割する数
!
SET LINE COLOR 4
!
! 描画
FOR t = 0 TO 2 * PI * rounds STEP 2 * PI / div
   PLOT LINES : xp(t), yp(t) ;
NEXT t
END
…

リスト2では、主にリスト1からの変更箇所を示しています。

これで例えば、 r1,r2,r3 を 31,16,11 、 rounds を 16 と入力して実行すると、次のような図形が描かれます。

アニメーションしていたら時間のかかるこのようなスピログラフも、一瞬で描けます。

2. スピログラフの性質

ここでは、スピログラフの持つ性質を調べていきます。

2.1. 軌跡が原点を通る条件

点Pの軌跡が原点を通る場合と通らない場合とがあります。

その条件は、どのように求まるでしょうか？

点Pが原点を通っている時は、次の図のような状態になります。

この時、 |r1| = |r2| + |r3| です。

さらに外トロコイドの場合も含めて考えれば、一般に点Pが原点を通る条件は、 |r3| = |r1 - r2| となります。

  1. r1=5,r2=3の場合
    2. r3=1 r3=2 r3=3
  2. r1=3,r2=-1の場合
    3. r3=3 r3=4 r3=5

真中の図だけ原点を通っていることが確認できます。

2.2. 周回数の法則

いろんなパラメータでスピログラフを描いていると、周回数 — つまり、何周したら最初の点に戻るか — には、それなりの法則性が感じられてきます。

r1とr2の関係で決まるようです。

まずは、いくつかのパターンに分けて考えてみます。

なお、周回数を確認する時は、経過を描いてくれるリスト1の方が適しているでしょう。

  1. r1がr2の倍数
    1周で終わります。
    ex. r1=12,r2=6 ⇒ rounds=1
  2. r2がr1の倍数
    r2/r1周で終わります。
    ex. r1=6,r2=12 ⇒ rounds=12/6=2
  3. r1とr2が互いに素（つまり、公約数が1だけ）
    r2周で終わります。
    ex. r1=13,r2=8 ⇒ rounds=8
  4. それ以外
    r2/(r1とr2の最大公約数)。
    ex. r1=12,r2=9 ⇒ rounds=3

…で、結局まとめると、周回数は r2/GCD(r1,r2) になります（ただし、GCDは最大公約数（Greatest Common Divisor））。

Aの場合、GCD(r1,r2)はr2ですし、BならGCD(r1,r2)=r1、CならGCD(r1,r2)=1というわけで、結局Dに集約されます。

なお、 r2/GCD(r1,r2) は、 LCM(r1,r2)/r1 でも同じです（ただし、LCMは最小公倍数（Lowest Common Multiple））。あてはめてみれば分かります。

スピログラフについて調べてたら、「スピログラフは初等整数論の入門になる」みたいなことが書いてあって（参考文献[3]）、「円を転がすことと整数論と何の関係があるんだろうな～」と思ったのですが、こういうことだったんですね。

いろいろなr1とr2の組み合わせを試すことで、最大公約数や最小公倍数が感覚的に身につくかもしれません。

2.3. 周回数の実装

周回数の法則がわかったので、それをプログラムに組み込みます。

リスト3. spirograph2.bas
…
DECLARE EXTERNAL PICTURE circle
DECLARE EXTERNAL FUNCTION gcd
!
…
LET ar3 = ABS( r3 )
!  周回数
IF INT(r1) <> r1 OR INT(r2) <> r2 THEN ! 非整数なら、以下の条件は適用しない
   LET rounds = 3
ELSE                                   ! 整数なら r2 ÷ (r1とr2の最大公約数)
   LET rounds = r2 / gcd( r1, r2 )
END IF
LET prompt$="周回数を入力（デフォルト値 " & STR$(rounds) & ")"
…
! 描画
FOR t = 0 TO 2 * PI * rounds STEP 2 * PI / div
…
NEXT t
! 最後に軌跡描画
SET LINE COLOR 4
PLOT LINES: prev_xp, prev_yp; xp(t), yp(t)
END
!
! （擬似）円描画ルーチン
!
EXTERNAL PICTURE circle( xc, yc, r )
…
END PICTURE
!
! 最大公約数関数：ユークリッドの互除法
!  参考：（仮称）十進BASICチュートリアル
!
EXTERNAL FUNCTION gcd( a, b )
   DO
      LET r = MOD( a, b )
      IF r=0 THEN EXIT DO
      LET a=b
      LET b=r
   LOOP
   LET gcd = b
END FUNCTION

リスト3では、主にリスト1からの変更箇所を示しています。

入力された r1,r2 を基に、rounds のデフォルト値として、適切な周回数を計算しています。

外部関数定義
外部関数定義は、外部絵定義と同じようにサブルーチンの一種で、絵を描く代わりに、何らかの値を返します。
これまで、 DEF 文で定義した関数が出てきましたが、あれをもっと複雑にできるようにしたものです。
まずは、プログラムの先頭で宣言をします。 
DECLARE EXTERNAL FUNCTION 関数名
外部関数本体は、次のように定義します。 
EXTERNAL FUNCTION 関数名(引数1,引数2,…)
   処理
   LET 関数名 = 戻り値
END FUNCTION
サブルーチン内で何らかの計算をし、その計算結果（戻り値（return value））を、関数名と同じ名前の変数に代入して返します。
関数の使い方は、 組み込み関数や DEF 文で定義した関数と同じで、式中で使います。
DO ～ LOOP 構文
FOR ～ NEXT 構文は、決まった回数の繰り返しでしたが、 DO ～ LOOP 構文 は、無限に繰り返すことができます。
もちろん、本当に無限に繰り返してはプログラムの実行が終わらないので、 EXIT DO 文 で終わらせます。
リスト3の gcd 関数では、「 r=0 になったらループを抜ける」という条件付けをしています。
余り
a を b で割った余りは、 MOD 関数を用いて、MOD(a,b) として求めます。
ユークリッドの互除法
「ユークリッドの互除法」は、アルゴリズムの教科書には必ず出てくるような、最大公約数を求めるための手順です。
（仮称）十進BASICについてくるPDF形式のチュートリアルにも載っています。
基本的な考え方は、「aもbもgcd(a,b)の倍数なのだから、aをbで割った余りもgcd(a,b)の倍数」といった感じでしょうか。
INT 関数
INT 関数は、「引数を超えない最大の整数」を返します。
INT(3.2) は3、 INT(-3.2) は-4です。
最大公約数を求めたりするのは、 r1,r2 が整数の時のみできることなので、リスト3では、「 INT(r1) が r1 自体と変わらなければ整数」という判定に使っています。
最後に軌跡描画
リスト3では、 END 文の手前で、ちょっと線を書き足してます。
なぜ書き足しているのかというと、その最後の PLOT LINES 文をコメントアウト(*3)して実行してみれば分かります。 
!PLOT LINES: prev_xp, prev_yp; xp(t), yp(t)
すると、最後のところにちょっと空白ができてしまいます。

なぜこのような空白ができるのかは、コラムで。

これで、適切な周回数を一々入力する必要がなくなりました。

計算精度と誤差
誤差の発生

以下のプログラムを実行してみましょう。

リスト4. verification-error.bas
LET n = 10
LET i = 0
FOR t=0 TO PI STEP PI/n
   PRINT "t(" ; STR$(i) ; ") =" ; t
   LET i = i + 1
NEXT t
PRINT "tx=";t
PRINT "PI=";PI
END

PRINT 文は、パラメータをテキストウィンドウに出力します。

パラメータは、間隔を空けずに表示したい時はセミコロン（ ; ）、でなければコンマ（ , ）で区切ります。

このプログラムは、 0 から PI までを n 等分し、それを t(i) として PRINT しています。

この場合、 0 から PI を10等分しているので、 t(0) ～ t(10) の11個の値が表示され、 t(10) は PI と同じになるものと期待されます。

しかし実行してみると、テキストウィンドウには次のように表示されます。

t(0) = 0
t(1) = .314159265358979
t(2) = .628318530717958 
t(3) = .942477796076937 
t(4) = 1.25663706143592 
t(5) = 1.5707963267949 
t(6) = 1.88495559215388 
t(7) = 2.19911485751286 
t(8) = 2.51327412287184
t(9) = 2.82743338823082
tx= 3.1415926535898
PI= 3.14159265358979

t(9) までしか表示されず、その値も PI にはだいぶ遠い値になっています。

なぜこのような実行結果になってしまったのかは、ループを抜けた後の t の値 tx と PI を比較すれば分かります。

tx の方が、 PI よりほんのちょっとだけ、大きくなっています。

FOR t=0 TO PI という文は、「 t が PI 以下の間だけ繰り返す」という意味なので、 t が計算誤差によってほんのちょっと PI を超えてしまったために、一足早くループが終わってしまったのです。

そして、この誤差がどのように発生するのかは、 t(1) ～ t(4) の末尾の桁に現れています。

t(1) から9、8、7と誤差無く来ているのに、 t(4) では小数点以下の桁数が1桁減り、そこで誤差が発生しています。

（仮称）十進BASICの数値変数の桁数は、標準では十進数15桁です。 t(4) では、 t(3) から一桁繰り上がって、それまでなかった一の位が現れてしまっているために、その分小数点以下の桁数が1桁減って誤差が生じた、ということです。

解決方法

一つの解決方法は、リスト3のように、ループを抜けた後で一つ余分に処理する方法です。

もう一つ、「あえて精度を下げる」という解決方法もあります。

以下のプログラムを実行してみましょう。

リスト5. verification-low-precision.bas
LET p = 3.14159265
LET n = 10
LET i = 0
FOR t=0 TO p STEP p/n
   PRINT "t(" ; STR$(i) ; ") =" ; t
   LET i = i + 1
NEXT t
PRINT "tx=";t
PRINT "p =";p
END

リスト4との違いは、円周率として PI ではなく、より精度を落とした p を定義して使っているということです。

出力は次のようになります。

t(0) = 0 
t(1) = .314159265 
t(2) = .62831853 
t(3) = .942477795 
t(4) = 1.25663706 
t(5) = 1.570796325 
t(6) = 1.88495559 
t(7) = 2.199114855 
t(8) = 2.51327412 
t(9) = 2.827433385 
t(10) = 3.14159265 
tx= 3.455751915 
p = 3.14159265

しっかり t(10) まで出力され、 t(10) と p は一致しています。

これは、 p の精度を落としたことで、10で割ったものを足していっても15桁を超えることが無くなり、途中で誤差が生じなくなったためです（ p を n で割ってしっかり割り切れ、加算していっても15桁を超えないことが条件）。

普通に考えれば、「精度が高い方が正確な結果が得られるだろう」ということになります。しかしこのように、精度が高いことでかえって誤差が生じてしまうこともあります。

大事なのは、前にも書いた通り「目的に応じた方法」です。

2.4. 経過を省略して描く

最後に、描画経過を省略したものを示します。

リスト6. spirograph3.bas
…
LET div = 100 ! 360°を分割する数
!
SET LINE COLOR 4
!
! 描画
FOR t = 0 TO 2 * PI * rounds STEP 2 * PI / div
   PLOT LINES : xp(t), yp(t) ;
NEXT t
! 最後に軌跡描画
PLOT LINES: xp(t), yp(t)
END
…

リスト6では、主にリスト3からの変更箇所を示しています。

3. おわりに

ここで示したプログラムも、改良の余地はいろいろあります。

例えば、参考文献[2]にあるようなイメージを、複数のスピログラフを重ねることで描いてみるのもいいかもしれません。

このイメージはパブリックドメインということなので、転載させてもらうことにします（ただし、60%に縮小しています）。

このイメージが掲載されているページには、描くのに使った設定が記されています。

例えば、左上の三角っぽいのは、

輪 歯車 穴
144/96 64 1, 3, 5, 7, ... , 25

となっています。

輪の「144/96」というのは、「外側の歯車の歯の数が144、内側の歯の数が96の輪」、歯車の「64」というのは「歯の数が64の歯車」を表すようです。

各歯は等間隔で並んでいるので、これはそのまま円周の比、つまりは半径r1とr2の比ということです。

また、「穴」欄の数値は、歯車の中の点ということで、歯の数との正確な比率は分かりませんが（歯車によって比率は違うかもしれません）、つまりは半径r3に対応します。

というわけで、上のスピログラフは、 r1,r2 を 96,64 にし、 r3 を 小さい値から 48 ぐらいまで重ねて描けば描けそうです。

そんな感じで、改造して遊んでいただければ幸いです。

たとえば…
トロコイド

第2回のリスト4をトロコイドに変えたらどうなるでしょうか。

スピログラフのサブルーチン化

リスト6のスピログラフ描画部分を spirograph( r1, r2, r3, rounds, c ) というサブルーチンにしたらどうなるでしょうか。ただし、 c はカラーインデックスを指定します。

『スピログラフ』 第3回 トロコイド (完)
【まとめ】
  1. ファイル庫/Spiro/figure/fig-spirograph.odg（OpenOffice.org Draw）にて出力。
  2. ファイル庫/Spiro/figure/fig-epitrochoid.odg（OpenOffice.org Draw）にて出力。
  3. プログラムの一部をコメントにして実行されないようにすること。コメントの記号を外せば、また実行できる。
参考文献
[1]. サイクロイドとスピログラフ（http://www.geocities.jp/sgwr0/sp/sp.html）
[2]. Wikipedia英語版「Spirograph」項目（http://en.wikipedia.org/wiki/Spirograph）
[3]. 2kspiro（http://www.zusaku.com/2kspiro.html）
BASIC 《 プログラミング言語 《 大系 《 彷徨の神殿
2006.12.25
Copyright © KHA-IN カイン 2006
inserted by FC2 system