Minecraft三角関数計算機報告

みなさん久しぶりです!(・ω・)ノ

あかげふです。

最近は専ら三角関数計算機にリソース全振りして制作してましたが、最新のupdateによるRender Dragonへの仕様変更により計算機開発環境を潰されたのと製作中の端末が限界を迎えつつあるので究極の三角関数計算機を作る野望は一旦損切りで打ち切りを決断しました。リソースを2ヶ月全振りしたにもかかわらずプロジェクトをお蔵入りに飛ばすのは流石に私の精神がズタズタになるので進捗報告して退去します(解決策とモチベが見つかったら帰ってくると思います)

信号強度をカラフルにするリソースパックを新versionで作成してくれる方がいたならば再開を考えようと思います。

去年2021年4月に私はCarry Cancel Adderを搭載したBEでは極小かつ最速の十進数12桁出力の加減乗除算機を制作した。そして中3からの野望であった究極的な三角関数計算機の制作に舵を切ったが、完璧主義的思考に囚われサイズの制約を厳しくしすぎた結果8月の長期休暇では完成に至らずモチベも消えたため暫く放置していた。今年2022年2月に再開を決めて、GoodNotesで設計図を詳細に書きながら計算機の構築を進めた。去年の8月に直面していた計算の技術面での困難を突破して無事に完成にいつくという三段であったが、ここからの道のりも長かった。2月末に数学物理好きの学部生が集団で宿泊して学問を語り合う数物セミナーというイベントに参加して準備等で忙しかったので数学にリソースを振っていたが、3月は本格的にMinecraftにリソースを振った。友達と3月にスキー旅行に行ったが、そのときに何度か転倒して半月板を損傷し、結局歩行に少し障害が残ることとなった。病院に診察に行くと手術をする見込みとなり、4月の中旬に全身麻酔を伴う手術をした。初めての体験であった。入院は1週間であったが、その間も痛みで夜寝られず、ストレス等で苦しみながら計算機の制作に励んだ。モチベはとうに枯渇して、数学や物理を学ぶほうが楽しいことも承知しておきながら中3からの意地を押し通すことだけで精神を突き動かして計算機の制作にのめり込んだのであった。計算機制作は設計ミスの連続で、メモリの数が1列不足したり、(パイプライン処理には回路全体の遅延を同期させておく必要があるわけだが)遅延量が合わなかったりが多発して、当初3月で終わるはずの計算機制作は4月までのびてしまったのであった。無事に退院して記事を書く今は筆者は3週間の松葉杖の生活なのであるが、ようやく慣れて健常者より歩行速度が早くなった。そんなさかな、悲劇が訪れた。1.18.30のupdateである。私は8年間BEのMinecraftをやり続けているが、Java勢と比べて計算機開発環境は決して良いものではなかった。何より時間を調整するツールが無いため、random tick rateを1000とか4000とかに引き上げることで処理を遅くして、1rstickの長さを伸ばして描画が追いつくようにして画面を録画して動作を確認するという手法を取っていた。そして私の計算機の大きな特徴である信号強度演算には、各点のワイヤーの信号強度を厳密に把握しておく必要があるため、瞬時に判断できるよう信号強度に応じて色がレインボーに変わる仕様にしてあった。このリソースパックは昔私にCarry Cancel Adderの技術を紹介してくれたYakumi氏が作成したものである。私はゲーム関連のプログラミングに疎いのでこのリソパが失効するのを恐れていた。最新のupdateにより、描画仕様が大きく変更されて現行のリソースパックは完全に使えなくなってしまった。完璧主義でかなり心を刳りながら進捗を生んでいたが、致命的な自体が起きた場合は即撤退すると決めていたので、updateの件が発覚した直後に撤退を決めた。

とまぁ経緯はこのくらいにして、技術について説明しましょうか。

Carry Cancel Adder(CCA)の仕組みについては僕の過去のblogを参照すればよい。今回は10進数演算でどのようにして三角関数を計算しようとしたかについて語る。

私の開発した10進CCA連続加算器は13桁である。三角関数を計算するにはCORDICというアルゴリズムが強力である。これは2進数の場合、事前計算したテーブルのROMとbit-shiftと加減算だけでsin,cosの計算を実現する。tanはsin/cosを除算器で求める。まず2進数でこのアルゴリズムを説明しておこう。

数3の知識を用いる。複素平面に $原点O(0),点P(x+yi),点Q(u+vi)$ があったとする。点Pを反時計に原点中心で $\theta$ だけ回転すると、その位置は

$(\cos\theta+i\sin\theta)(x+yi)$

となる。点Qを $u+vi=(1+i\tan\theta)(x+yi)$ とすると、 $△OPQ$ は $∠OPQ=\pi/2$ かつ $∠POQ=\theta$ の直角三角形である。このように点Pを $\theta$ だけ回転して、半直線OP方向に $1/\cos\theta$ 倍にしてQの位置に移動する操作を便宜的に独自に「 $\theta$ だけ擬回転する」と呼称しよう。 $|OQ|= |1+i\tan\theta| |x+yi|=\frac{1}{\cos\theta} |OP|$

であるので、「 $\theta$ だけ擬回転する」と「 $-\theta$ だけ擬回転する」場合で|OQ|の長さに違いがないことがわかる。「半径が擬回転の符号によらないこと」が1つめのポイントである。

$\tan\theta=T$ とおくと擬回転で成分は

$u+vi=(x-Ty)+(y+Tx)i、(u,v)=(x-Ty,y+Tx)$ というように変化する。２進数で計算する場合、 $T=2^{-n}$ であるとき、 $Ty$ は $y$ を下位にn bit だけshiftすると求められるので乗算する必要がない。そこで、 $\tan\theta_n=2^{-n}(n=0,1,...)$ となる $\theta_n$ を用意しておいて、 $\varepsilon_n\in \{-1,1\}$ を使って

\begin{align}\theta=\sum_{n=0}^N \varepsilon_n \theta_n\end{align}

と表せるとしよう。点 $P_0(R+0i)$ (R>0)を「 $\varepsilon_0 \theta_0$ 擬回転」したものを点 $P_1$ 。そして点 $P_1$ を $\varepsilon_1\theta_1$ 擬回転」したものを点 $P_2$ 。...「 $\varepsilon_N \theta_N$ 擬回転」したものを点 $P_N$ 。という操作を続けざまに行うと、最終的にx軸となす角が $\theta$ の点 $P_N$ に行き着く。この点の座標は複素数で書くと

\begin{align}\prod_{n=0}^N(1+\tan\varepsilon_n \theta_n)(R+0i)\\ =(1+\tan\varepsilon_0 \theta_0)(1+\tan\varepsilon_1 \theta_1)\cdots(1+\tan\varepsilon_N \theta_N)(R+0i)\end{align}

となる。擬回転1回はbit shiftして足すという操作だけであったので、これをN+1回合成すれば $P_N$ の座標がわかる。 $|OP_N|=1$ となるように $R$ を規格化する。

\begin{align}1=|OP_N| &=|1+\tan\varepsilon_0 \theta_0| |1+\tan\varepsilon_1 \theta_1|\cdots|1+\tan\varepsilon_N \theta_N| |R+0i|\\ &=\frac{R}{\cos\theta_0\cos\theta_1\cdots \cos\theta_N}\end{align}

こうすることで、 $P_N$ のx座標は $\cos\theta$ , $P_N$ のy座標は $\sin\theta$ となる。ここまでをまとめると、 $P_n(x_n,y_n),\tan\theta_n=2^{-n},0\lt \theta_n\lt\pi/2,\varepsilon_n\in \{-1,1\},(n=0,1,...,N)$ $\theta=\varepsilon_0 \theta_0+\varepsilon_1 \theta_1+\cdots+\varepsilon_N \theta_N$ の形で $\theta$ を分解できたとき、 $(x_0,y_0)=(R,0)$

\begin{align}(x_{n+1},y_{n+1})=(x_n-\varepsilon_n2^{-n}y_n,y_n+\varepsilon_n2^{-n}x_n)\end{align}

という漸化式に従うと $P_N(x_0,y_0)=(\cos\theta,\sin\theta)$ を得る。

あとは $\theta=\varepsilon_0 \theta_0+\varepsilon_1 \theta_1+\cdots+\varepsilon_N \theta_N$ という分解の仕方を考えればよい。 $\theta_n\fallingdotseq\tan\theta_n=2^{-n}$ という近似をすれば、 $\theta\fallingdotseq\varepsilon_0 2^{-0}+\varepsilon_12^{-1}+\cdots+\varepsilon_N 2^{-N}$ となる。

つまり2進数と同様の分解の仕方で考えることができる(この事の正当性については後で挙げる記事を参照してほしい)

$\theta$ の符号を $\varepsilon_0$ として、

$\theta-\varepsilon_0\theta_0$ の符号を $\varepsilon_1$ として、

...、

$\theta-\varepsilon_0\theta_0-\cdots-\varepsilon_{N-1}\theta_{N-1}$ の符号を $\varepsilon_N$ とすればよい。

なお一般の $\theta$ に対し $\theta-\varepsilon_0\theta_0-\cdots-\varepsilon_{N}\theta_{N}\neq 0$ であるが、

打ち切り誤差を発生させて $\theta-\varepsilon_0\theta_0-\cdots-\varepsilon_{N}\theta_{N}\fallingdotseq 0$ とする。

この話を10進数化して考える。

擬回転で変化する成分は $u+vi=(x-(\tan\theta)y)+(y+(\tan\theta)x)i$ となって、これが計算しやすいように2進数の場合は「 $\tan\theta=2^{-n}$ 」となるように選んだのであった。安直に $\tan\theta_n=10^{-n}$ として考えてみる。CORDICのポイントである「擬回転したあとの半径は符号(回転が時計か反時計か)によらない」ことから、

$\theta=\varepsilon_0 \theta_0+\varepsilon_1 \theta_1+\cdots+\varepsilon_N \theta_N$

という形の分解を考えることになる。 $\theta\fallingdotseq 10^{-n}$ という近似で考えてみるとわかるが、つねに $\theta\fallingdotseq\varepsilon_0 10^{-0}+\varepsilon_1 10^{-1}+\cdots+\varepsilon_N 10^{-N}$ となる $\varepsilon_n\in \{-1,1\}$ が存在するとは限らない。結論を言うと、分解の仕方は次のようになる。

$\ N=4M,\ n=4s+t,\ t\in\{1,2,3,4\},(c_1,c_2,c_3,c_4)=(4,2,2,1)$
$\tan\theta_0= 1,\tan\theta_n= \dfrac {c_t}{10^{s+1}}$

$\theta=\varepsilon_0 \theta_0+\varepsilon_1 \theta_1+\cdots+\varepsilon_N \theta_N$

このような分解の仕方に至った考え方や、定義域内の $\theta$ に対して $\theta=\varepsilon_0 \theta_0+\varepsilon_1 \theta_1+\cdots+\varepsilon_N \theta_N$ のかたちの分解 $(\varepsilon_0 ,\varepsilon_1 ,\ldots,\varepsilon_N)$ が存在すること、そして再帰的に分解できることは自明ではないが、私のmathlogの記事に書いた。

mathlog.info

\begin{align}R&=\cos\theta_0\cos\theta_1\cdots \cos\theta_N\\ &=\frac{1}{\sqrt{1+\tan\theta_0^2}}\frac{1}{\sqrt{1+\tan\theta_1^2}}\cdots\frac{1}{\sqrt{1+\tan\theta_N^2}}\\ &=\frac{1}{\sqrt{1+10^{-2×0}}}\frac{1}{\sqrt{1+1^2×10^{-2×1}}}\frac{1}{\sqrt{1+2^2×10^{-2×1}}}\frac{1}{\sqrt{1+2^2×10^{-2×1}}}\frac{1}{\sqrt{1+4^2×10^{-2×1}}}\cdots\frac{1}{\sqrt{1+4^2×10^{-2M}}}\end{align}

擬回転の成分変化は

$(x_0,y_0)=(R,0)$

\begin{align}(x_{n+1},y_{n+1})=(x_n-\varepsilon_n\dfrac {c_t}{10^{s+1}}y_n,y_n+\varepsilon_n\frac {c_t}{10^{s+1}}x_n)\end{align}

となる。 $c_t$ は掛け算を行うのではなく、 $c_t回\frac {\varepsilon_nx_n}{10^{s+1}},\frac {\varepsilon_ny_n}{10^{s+1}}$ を足し引きすると考える。

こうすることで10進数でCORDICを行うことができ、加減乗除とシフト( $10^{-n}$ で計算を実行できる。ROMに記憶させる必要があるのは事前に計算した $\theta_0,\theta_1,...,\theta_N,R$ の値である。

以上が僕が独自に10進数用に改変したCORDICアルゴリズムであるが、Minecraftで究極的な計算機を作るためにはさらに工夫を要した。その結果制御配線が複雑になった。そのゴミのような回路群を解説する気はないが、本質的なアイデアについては説明しようと思う。そもそもMinecraftでは計算機に使いやすいROM、特に10進数計算なので信号強度を所定のタイミングで引き出せるROMは構造が限られてくるし、10進数1桁あたりの大きさは3.75×2×2ブロックである。なので記憶する数を増やすと列数が増えて制御配線が大きくなるし、サイズも大きくなる。制御配線を増やしてでも記憶素子の数を減らすことが重要である。繰り上がり無遅延でとどく最大のCCAは13桁で、下位2,3桁は誤差として棄却するとした。度数法での計算で、入力は上位2桁が整数、下位11桁が小数部、出力は上位1桁が整数部とした。なので結果の打切り誤差や累積する加算の誤差は $10^{-12}$ 、 $\theta$ の誤差は $10^{-11}$ 程度になるように調整した。

私が製作中だった計算機では

$\theta=\varepsilon_0 \theta_0+\varepsilon_1 \theta_1+\cdots+\varepsilon_N \theta_N$

という分解の過程を3段階に分割して極力記憶素子の量をへらすようにした。

最初は通常のCORDICどおりのアルゴリズムで計算する。

次にn>N/3のときは段階が変わる。

$\theta_k=\arctan z$ をテイラー近似すると $\theta_k=z-z^3/3$ となる。 $z^3/3\fallingdotseq 10^{-12}$ 程度であれば $\theta_k=z$ という線形近似で十分である。だからN/3に近い4で割ると1あまる整数Kをもってくると $n\gt K$ のときは $\theta_nは\theta_K$ を下位にシフトして $c_t$ 倍すると得られる。

さらにn>N/2のときは段階が変わる。

擬回転の成分変化ではx,y成分のシフトを行うが、シフトではみ出した分は加算器の桁数からはみでるので捨てられる。n>N/2のときシフトする数は加算器の桁数の半分以上なので、2度同じ桁数をシフトすると加算器上では丸め誤差で0になる。このことから擬回転による成分変化の漸化式を簡約化することができる。L=N/2=2Mとするとn>Lのとき

\begin{align}(x_{n+1},y_{n+1})=(x_n-\varepsilon_n\dfrac {c_t}{10^{s+1}}y_n,y_n+\varepsilon_n\frac {c_t}{10^{s+1}}x_n)\end{align}

という漸化式は

\begin{align}(x_{n+1},y_{n+1})=(x_n-\varepsilon_n\dfrac {c_t}{10^{s+1}}y_L,y_n+\varepsilon_n\frac {c_t}{10^{s+1}}x_L)\end{align}

となる。計算機の内部構造に関する話をすると、中にはCCAの連続加算器が搭載されているので加算を次々行うことはできるが、シフトはできないので一旦値を取り出してシフトレジスタでシフトするステップがあるのでシフト数が増えると遅延量が増える。この近似を行うことにより、シフトさせるのは $(x_L,y_L)$ という量なのでいちいち加算器から値を取り出す必要がなく、高速化につながる。さらに私の前回の乗除算機の仕組みを見て類推してもらうと、第３段階は乗除算の構造と全く同じとなる。適切に近似を行った結果、第3段階はCORDICのアルゴリズムの手続きからは外れたものとなった。 $\theta$ の部分剰余を $\theta_L$ で回復型除算で割った結果で得られる符号をダイレクトに $(x_L,y_L)$ (をシフトしたもの)の足し引きの符号につなげる。