「円周率をそんなに求めて何の役に立つの？」と思う人は多いでしょう。 はっきり言って、ほぼ自己満足の世界です。 コンピュータの耐久試験くらいにしか役に立ちません。
Z80 というのは 1980 年頃によく使われた 8bit CPU なのですが、「なぜ今さらそんな CPU で？」と思うでしょう。 はい、自己満足です。
そういったことは気にせず、とにかく高速に計算したいのです。意地です。
1985～1986 年あたりの「I/O」というパソコン雑誌に、円周率 1 万桁（正確には「小数第 1 万位まで」ですが「桁」と表現することにします） を求めるプログラムがいくつか載っていました。 主に 6809 という CPU が対象でしたが、Z80 を対象にしたプログラムもあり（I/O 1985 年 10 月号、機種は FP-1100）、計算時間はたしか 1 時間 48 分 55 秒でした。
そのプログラムを PC-8001mkIISR というパソコンに移植し、6809 用のプログラムも参考にしながら高速化していました。1 時間 2 分程度まで縮まったものの、1 時間を切ることができませんでした。
そのことが長年ずっとモヤモヤしていたのですが、 最近になってふと再チャレンジしてみて、最終的には 1607 秒（26 分 47 秒）と、 半分以下の時間で求めることができるようになりました。 当時の自分が知ったら驚愕したことでしょう。
円周率にまつわる式は非常に多いので、詳しく知りたい方は Wikipedia の 円周率 などを読んでください。
ある程度以上のマシンパワーで多くの桁数を求めるなら ガウス＝ルジャンドルのアルゴリズム や Chudnovsky の公式 などを使用するのですが、非力な 8bit CPU で 1 万桁程度を求める場合は、arctan の組み合わせで求めるのが一番速いと思います。
その系統の式は無数にありますが（詳しくは
arctan 関係式一覧
を参照）、その中でも有名なものの一つとしてマチンの公式、
$$ \frac{\pi}{4} = 4 \arctan \left( \frac{1}{5} \right)
- \arctan \left( \frac{1}{239} \right) $$
があり、今回の 1 万桁計算ではこの公式が一番高速でした。
【'19/05/25 追記】
その後の改良により、マチンの公式よりもシュテルマーの公式（その 1）の方が速くなりました。
arctan は、
$$ \arctan(x) = \sum_{n=0}^{\infty} \frac{(-1)^n}{2n+1}
x^{2n+1} $$
で求めることができます。マチンの公式の両辺を 4 倍した式、
$$ \pi = 16 \arctan \left( \frac{1}{5} \right)
- 4 \arctan \left(\frac{1}{239} \right) $$
を上記の arctan の式に当てはめると、
$$ \begin{align}
\pi &= 16 \left( \frac{1}{5} - \frac{1}{3 \times 5^3}
+ \frac{1}{5 \times 5^5} - \frac{1}{7 \times 5^7} + \cdots \right) \\
&- 4 \left( \frac{1}{239} - \frac{1}{3 \times 239^3}
+ \frac{1}{5 \times 239^5} - \frac{1}{7 \times 239^7} + \cdots \right)
\end{align} $$
となります。第 1 項を変形すると、
$$ \frac{16}{5}
- \frac{16}{5 \times 5^2 \times 3}
+ \frac{16}{5 \times 5^2 \times 5^2 \times 5}
- \frac{16}{5 \times 5^2 \times 5^2 \times 5^2 \times 7}
+ \cdots $$
となるので、16÷5 から始めて 52
で割っていき、それを奇数で割った値を加減算していくことになります。
第 2 項も同様に計算しますが、係数が -4 で負号が付いており、 展開すると加算と減算が逆になるので注意してください。
少し長くなりますが、計算手順は以下のようになります。
WORK1～3 は多倍長整数（任意長の固定小数点数）で、 整数部に数桁と、小数部に求める桁数が最低限必要です。
求める桁数を $ n $ とすると、必要なサイズは $ n \div \log_{10} 256 $ Byte（＋整数部 1Byte）程度です。1 万桁なら約 4153+1 Byte になります。 ただし加減算で計算誤差が蓄積するので、数バイト多めに確保しておきます。
d は 2 ずつ増えていきますが、arctan の分母を $ x $、求める桁数を $ n $ とすると、d は最大で $ n \div \log_{10} x $ 程度になります。 マチンの公式で 1 万桁なら約 14307 なので、14bit あれば足ります（214=16384）。
最後の「10 進数に変換」ですが、多倍長整数は 2 進数表現なので、人間用に 10 進数へ変換する必要があります。 今回は 1Byte のうち上位 4bit と下位 4bit にそれぞれ 10 進数 1 桁を格納します（「packed BCD（Binary Coded Decimal）」と呼ばれます）。
BCD_RESULT が 10 進数の格納先です。初めの 1Byte が整数部、その後に小数部が続きます。
これらを踏まえると、主に以下の計算ルーチンが必要になります。
この計算手順の確認ということで、C
言語で簡易的に（多倍長整数の代わりに倍精度実数を使用して）15
桁ほどを求めるサンプル machin.c
を置いておきます。手順自体はそれほど難しくないかと思います。
Z80 を積んだ実機を用意するのは大変ですし、 仮に用意しても計算に時間がかかりすぎるので、Z80 エミュレータを作成しました。このエミュレータだと普通の PC で 1 万桁の計算が 10 秒程度で終わります。
さらに Z80 アセンブラも作成し、アーカイブに C/C++ ソースと Windows 用実行ファイルを収めてあります。
ちなみにエミュレータは昔作ったものをベースに簡略化したものです。 アセンブラは一から作りました。バグってませんよーに。
実装の解説に入る前に、今回決めた条件を書いておきます。
ここからは実装の解説に入りますが、Z80 で多少のプログラミングをできるという前提で解説します。 そうでない方はこちらで勉強してください。
……というのも厳しすぎるので、C 言語で記述してから Z80 に置き換える、という流れで解説していきます。
まずは比較的簡単な多倍長整数同士の加減算を C 言語で書いてみます。 筆算と同じように、下のほうの位から計算していきます。
typedef unsigned char u8;
/* size バイトの dest[] に src[] を足す */
void MPAdd(u8 dest[], u8 src[], int size)
{
int i;
int carry = 0;
for(i = size - 1; i >= 0; i--) {
int n = dest[i] + src[i] + carry;
carry = (n >> 8) & 1;
dest[i] = (u8)n;
}
return;
}
/* size バイトの dest[] から src[] を引く */
void MPSub(u8 dest[], u8 src[], int size)
{
int i;
int carry = 0;
for(i = size - 1; i >= 0; i--) {
int n = dest[i] - src[i] - carry;
carry = (n >> 8) & 1;
dest[i] = (u8)n;
}
return;
}
加算も減算も処理はほとんど同じです（n の計算の 1 行が違うだけです）。
ちなみに繰り下がりは英語で（carry ではなく）borrow と言い、CPU によっては 0/1 が逆になりますが、今回は Z80 の動作に合わせます。
さて、加算ルーチンを Z80 に置き換えてみます。たとえば WORK1 に WORK3 を足す場合は、DE レジスタに WORK1、HL レジスタに WORK3、BC レジスタにワークサイズをセットして呼び出すことにします。
ちなみに今回、ループには Z80 講座の「16bit ループの技」を多用しています。
; [DE..DE+BC-1] += [HL..HL+BC-1] MPADD: DEC BC ADD HL,BC EX DE,HL ADD HL,BC EX DE,HL INC C ; 16bit ループの技（ただし DEC BC は初めに実行済み） INC B LD A,C LD C,B LD B,A ; Cフラグは最後に実行した ADD HL,BC で 0 になっている MPADD_L1: LD A,(DE) ADC A,(HL) ; 減算なら SBC A,(HL) LD (DE),A DEC DE DEC HL DJNZ MPADD_L1 DEC C JR NZ,MPADD_L1 RET
早速 C フラグのクリアを省略するというケチなことをしていますが、 こういうのはバグの元になるのでお勧めはできません。でも C フラグの変化は体に染み付いているので、反射的にやってしまいます。
減算は ADC A,(HL) を SBC A,(HL)
に替えるだけです。これで多倍長整数同士の加減算が完成しました。
ソースは mpadd1.z80s です。
次に多倍長整数を小さな整数（16bit 以下）で割る処理を実装します。 円周率計算の大部分を占める計算で非常に重要なので、 何度も改良していきます。
筆算と同じように、上のほうの位から計算していきます。 C 言語ではこんな感じです。
typedef unsigned char u8;
/* size バイトの src[] を d で割って dest[] に格納する */
void MPDiv(u8 dest[], u8 src[], int d, int size)
{
int i;
int r = 0;
for(i = 0; i < size; i++) {
r = (r << 8) + src[i];
dest[i] = (u8)(r / d);
r %= d;
}
return;
}
しかし除算（/）と剰余算（%）が入っていて、このままでは Z80 に置き換えるのが面倒です。そこで 1bit ずつコツコツ求めることにします。
ちなみに 10 進数の筆算では「除数の何倍まで引けるか」というのを 0～9 倍まで考えなければいけませんが、2 進数なら 0 と 1 しかありませんから、「除数が引けるか」だけで済みます。
typedef unsigned char u8;
/* size バイトの src[] を d で割って dest[] に格納する */
void MPDiv(u8 dest[], u8 src[], int d, int size)
{
int i, j;
int r = 0;
for(i = 0; i < size; i++) {
u8 q = 0;
u8 n = src[i];
for(j = 0; j < 8; j++) {
r <<= 1;
if ((n & 0x80) != 0) {
r++;
}
q <<= 1;
if (r >= d) {
r -= d;
q++;
}
n <<= 1;
}
dest[i] = q;
}
return;
}
ちょっと難しいかもしれませんが、これでかなり Z80 に置き換えやすくなったかと思います。
さて、これを Z80 に置き換えるわけですが、まず内側のループ内だけ考えてみます。
; HL=余り(r), DE=除数(d), A=非除数(n), C=商(q) SLA C ADD A,A ADC HL,HL SBC HL,DE JR NC,L1 ADD HL,DE JP L2 L1: INC C L2:
Z80 には 16bit 比較命令が無いので、とりあえず引いてみて C フラグが立ったら足して戻しています。
これを 8 回ループするか 8 個並べれば 1Byte 求まりますが、その前に高速化してみます。まず商を C レジスタに 1bit ずつ格納していますが、A レジスタも 1bit ずつ空いていくので、ここに格納してしまいましょう。
; HL=余り(r), DE=除数(d), A=非除数(n):商(q) ADD A,A ADC HL,HL SBC HL,DE JR NC,L1 ADD HL,DE JP L2 L1: INC A L2: ; ここから2bit目 ADD A,A
2bit 目を求める時点でも C
フラグが残っているので、ADC A,A とすれば INC A
を省けそうです。しかしよく考えると 0 と 1 が逆になるので、最後に
CPL でまとめて反転してしまいましょう。
; HL=余り(r), DE=除数(d), A=非除数(n):商(q) ADD A,A ADC HL,HL SBC HL,DE JR NC,L1 ADD HL,DE L1: ; ここから2bit目 ADC A,A ADC HL,HL : L7: ; ここから8bit目 ADC A,A ADC HL,HL SBC HL,DE JR NC,L8 ADD HL,DE L8: ADC A,A CPL
これで内側のループは完成です……と言いたいところですが、
実は一つ問題があります。除数が 15bit
に収まっていればいいのですが、それを超えると
ADC HL,HL
で 16bit からはみ出てしまうことがあるのです。
マチンの公式では第 2 項で 2392=57121 での除算が必要で、これは 15bit の最大値 32767 を超えていますから、除数 16bit 用の除算ルーチンも用意することにします。
ようは ADC HL,HL で C
フラグが立ったら、必ず DE が引けるものとして処理します。
; HL=余り(r), DE=除数(d), A=非除数(n):商(q) ADD A,A ADC HL,HL JR NC,L1 OR A SBC HL,DE OR A JP L2 L1: SBC HL,DE JR NC,L2 ADD HL,DE L2: ; ここから2bit目 ADC A,A
C フラグの変化が都合悪く、クリアするために入れている
OR A をなんとかしたい感じです。
SUB 命令に 16bit 版が無いのも痛いですね。
そこで「DE を引く」代わりに「-DE を足す」ことで C フラグの変化が逆になる（ただし DE≠0）ことを利用します。BC レジスタに -DE をセットしておくと、
; HL=余り(r), DE=除数(d), BC=-除数(-d), A=非除数(n):商(q) ADD A,A ADC HL,HL JR NC,L1 ADD HL,BC JP L2 L1: SBC HL,DE JR NC,L2 ADD HL,DE L2: ; ここから2bit目 ADC A,A
と書くことができます。
とりあえず高速化はここまでにしておきましょう。
外側のループも付け加えたものが mpdiv1.z80s です。
整数を 2 進数から 10 進数に変換したことはあっても、 小数を変換したことのある人は少ないのではないでしょうか。
整数では「10 で割った余り」を繰り返し求めることで変換できますが、小数では「小数部を 10 倍した整数部」を繰り返し求めることで変換できます。 たとえば 3.1415 なら、0.1415 を 10 倍して 1.415、0.415 を 10 倍して 4.15、といった要領です。
これを C 言語で書いてみます。ただし packed BCD では 10 進数 2 桁ずつ格納しますから、10 倍ではなく 100 倍して 2 桁ずつ求めることにします。
100 倍するのは筆算と同じように、下のほうの位から計算していきます。
typedef unsigned char u8;
/* binSize バイトの bin[] を packed BCD に変換して
bcdSize バイトの bcd[] に格納する
ただし bin[], bcd[] は小数部とする */
void Bin2Bcd(u8 bcd[], int bcdSize, u8 bin[], int binSize)
{
int i, j;
for(i = 0; i < bcdSize; i++) {
int n = 0;
for(j = binSize - 1; j >= 0; j--) {
int m = bin[j] * 100 + n;
bin[j] = (u8)m;
n = m >> 8;
}
bcd[i] = (u8)(((n / 10) << 4) + (n % 10));
}
return;
}
これを Z80 に置き換えるのですが、まず 100 倍する部分をどう置き換えるかが問題になります。Z80 には乗算命令が無いので、他の方法で乗算しなければなりません。
たとえば「HL=A*100」は、100 の 2 進数表現「1100100」を踏まえて、
LD H,0 LD L,A ; HL=A LD D,H LD E,L ; DE=A ADD HL,HL ; HL=A*2 ADD HL,DE ; HL=A*3 ADD HL,HL ; HL=A*6 ADD HL,HL ; HL=A*12 ADD HL,HL ; HL=A*24 ADD HL,DE ; HL=A*25 ADD HL,HL ; HL=A*50 ADD HL,HL ; HL=A*100
と書けます。これはこれで普通の方法なのですが、 高速化のためにデータテーブルを使用することにします。
8bit 値（256 通り）を 100 倍して 16bit 値を得るので、16bit×256=512Byte のテーブルになります。
; MUL100TBL_H=データテーブル先頭の上位バイト LD HL,MUL100TBL_H*100H LD D,L ; L=0を利用して節約 LD E,L LD BC,100 L1: LD (HL),E INC H LD (HL),D DEC H EX DE,HL ADD HL,BC EX DE,HL INC L JR NZ,L1
0 から始めて 100 ずつ足しているだけです。 なお、テーブルを参照しやすいように 100H の倍数アドレスから格納しています。
あと DE レジスタを 0
クリアするところで、ちょっとセコい技を使っています。LD DE,0
よりも 1Byte と 2clk の節約になっています。
ではこのテーブルを使って、多倍長整数の小数部分を 100 倍してみます。整数部は 1Byte とします。
; [HL+1..HL+BC-1] *= 100 DEC BC ADD HL,BC DEC BC ; 整数部1Byteを除外 INC C ; 16bit ループの技（ただし DEC BC は初めに実行済み） INC B LD A,C LD C,B LD B,A XOR A LD D,MUL100TBL_H L1: LD E,(HL) EX DE,HL ADC A,(HL) EX DE,HL LD (HL),A DEC HL INC D LD A,(DE) DEC D DJNZ L1 DEC C JR NZ,L1 ADC A,B ; Cフラグ（繰り上がり）を足す（Bレジスタは0）
これで A レジスタに 0～99 の値が求まりますが、欲しいのは packed BCD 形式の 00H～99H です。また除算するのもちょっと面倒ですね。
そこでこれもテーブルで処理することにしましょう。 100Byte のテーブルになります。
; BCDTBL_H=データテーブル先頭の上位バイト LD HL,BCDTBL_H*100H XOR A L1: LD (HL),A INC L ADD A,1 DAA JR NC,L1
1 足して DAA 命令で 10 進補正、を繰り返すだけです。 DAA 命令のおかげで簡単ですね。
ここまでをまとめて、外側のループも付け加えたものが
tobcd1.z80s です。
ここまで作ったルーチンにメインルーチンを加えて、円周率 1
万桁を求めてみます。まずは C 言語版
pi1.c です。
少し補足すると、main() 関数の work2Skip という変数に、work2 の先頭から何バイトが 0 になったかを格納しています。 これがワークサイズ（binSize）になったらループ終了です。
コンパイルして実行すると、1～2 秒ほどで計算結果を 10010 桁表示して終了します。でもこれでは結果が正しいかどうかわかりませんね。
アーカイブの中の
pi100k.txt に円周率 10 万桁が入っています。
比較すると 10006 桁まで正しいことがわかります。
手作業で比較するのは面倒なので、Z80 版では自動化することにします。
次に Z80 版です。メインルーチンを加えたのが
pi1.z80s
です。プログラムについての補足は、
といったところです。
メモリ配置については、プログラムを 0000H から、その後に除算テーブルなどのワーク、続いて WORK1～3 という配置にしました。計算結果（BCD_RESULT）は WORK1 を避けて WORK2～WORK3 にまたがるような配置にすれば大丈夫です。
それでは Z80 ソースをアセンブルし、エミュレータで実行してみましょう。
アーカイブの中にある runpi.bat
を実行するだけです。
ソースがアセンブルされると pi.z80
というバイナリファイルが作成され、
エミュレータがそのファイルを読み込んで実行します。
しばらく（数十秒程度）待つと、
40460522295 clk (10132 sec., 2:48:52) Match : 10006 digit(s) 続行するには何かキーを押してください . . .
と出ると思います。10132 秒で 10006 桁まで正確な値が求まりました。
なお、コマンドプロンプトを開いて runpi.bat
を実行する場合、第 1 引数に Z80
ソースファイル名（省略時 pi1.z80s）、第 2 引数に BCD
格納アドレス（16 進数、省略時 6000）を指定することもできます。
他のプログラムを実行する場合に利用してください。
また、エミュレータは 64KB のメモリを全て
mem.bin
に出力して終了します。デバッグに利用できます。
除算のたびに上位の桁から 0 になっていきますが、 今はその桁も毎回除算していて非常に無駄です。
PWORK2_SKIP にこのバイト数が入っているので、 そこを飛ばして除算することで無駄が省けます。 その後の加減算のときにも利用します。
同様に 10 進数への変換で、100 倍するたびに下位の桁が少しずつ意味の無い値になっていきます。 具体的には $ \log 100 \div \log 256 $（約 0.83）Byte ずつ計算を減らすことができます。
これらの改良をしたのが
pi2.z80s
です。計算時間は 4766
秒と、一気に半分以下になりました。
これは I/O 1986 年 9 月号に載っていた 6809 版の手法です。6809 版は 52=25 での除算専用でしたが、今回は任意の 8bit 値での除算に対応しました。
まずバイト単位で除算する手順を確認します。被除数の 1Byte 目を p1、除数を d とし、商 q1 と余り r1 を求めます。q1 はそのまま商の 1Byte 目になります。
次に被除数の 2Byte 目を p2 とし、r1×256+p2 を d で割って商 q2 と余り r2 を求めます。q2 が商の 2Byte 目になります。以下これを繰り返します。
さて、r1×256+p2 を d で割る際に、r1×256 と p2 に分けて割ることにします。 まず r1×256 を d で割った商を s1、余りを t1 とします。そして p2+t1 を d で割った商を Q2 とすると、q2=s1+Q2 になります（余りは r2 になります）。
3Byte 目以降も同様に考えて整理すると、p2+t1 から Q2、s2、t2 の 3 つの値を求めるテーブルを作成することで、 後は簡単な計算をするだけで済みそうです（ちなみに r は必要無くなります）。 ただし p2+t1 は 9bit になりますが、テーブルのサイズを節約するために、256 以上の場合は補正処理を行うことにします（後述します）。
こうして 8bit 値から 3Byte の値を求めるテーブルになり、サイズは 768Byte になります。
この処理を C 言語で書いてみます。 ただしテーブル参照はせず、それぞれを毎回計算します。
typedef unsigned char u8;
/* size バイトの dest[] を d（1～0xff）で割って dest[] に格納する */
void MPDiv8(u8 dest[], u8 src[], int d, int size)
{
int i;
u8 s = 0, t = 0;
u8 qq = 256 / d;
u8 rr = 256 % d;
for(i = 0; i < size; i++) {
u8 p = src[i];
int n = p + t;
u8 q = 0;
if (n >= 256) {
n = n - 256 + rr;
q = qq;
}
u8 tq = n / d;
u8 r = n % d;
u8 ts = (r << 8) / d;
t = (r << 8) % d;
q += s + tq;
dest[i] = q;
s = ts;
}
return;
}
初めのほうに qq, rr という変数がありますが、これが p2+t1 が 256 以上になった場合の補正値です。 それぞれ 256÷d の商と余りなので、この値はテーブルの中に存在しています （p2+t1=1 のときの s2 と t2）。
次に Z80 のソース
mpdiv2.z80s
です。MAKE_DIVTBL がテーブル作成ルーチン、MPDIV8T
が除算ルーチンです。結構ややこしくて難しいかもしれませんが、
コメントを入れてあるので参考にしてください。
この除算ルーチンを組み込んだのが
pi3.z80s
です。計算時間は 3482 秒（58 分 02
秒）で、あっさり 1 時間を切りました。
ここまでが、ほぼ昔作ったプログラムです。 当時のプログラムは発掘していないので、記憶を頼りに作り直したのですが、 なぜか当時よりも速くなってしまいました。
たしか当時は、奇数での除算を直接 WORK3=WORK2÷d とせず、WORK2 から WORK3 へコピーしてから除算していたと思います。 コピーにそんなに時間がかかっているとは思わなかったのです。 これに気付いたときは拍子抜けでした。
気を取り直して、さらなる高速化を目指します。 まずは各計算にどれくらいの時間がかかっているか調べてみましょう。
実はエミュレータにはちょっと細工がしてあります。 通常のクロックカウンタとは別にカウンタが 3 つあり、
LD B,B でカウンタ 1 起動、LD C,C で停止LD D,D でカウンタ 2 起動、LD E,E で停止LD H,H でカウンタ 3 起動、LD L,L で停止という動作をします（それぞれの命令自体は 0clk に変えてあります）。計測したい部分の前後をこれらの命令で挟むことで、 部分的なクロック数（時間）を計測することができます。
主要な部分の計測結果は以下の通りです。
| CALC1 | 2223 | ÷52 | 383 |
|---|---|---|---|
| ÷(2n+1) | 1667 | ||
| +/- | 172 | ||
| CALC2 | 1095 | ÷2392 | 554 |
| ÷(2n+1) | 490 | ||
| +/- | 51 | ||
| CONVERT | 165 | ||
| total | 3482 |
奇数での除算がかなりの時間を占めています。 それ以外で少し気になるのが 2392 での除算ですが、試しに 1÷239（テーブル利用）と 1÷2392 の計算クロック数を計測してみましょう。
| 1÷239 | 424397 |
|---|---|
| 1÷2392 | 1871155 |
4 倍以上の差が付いていますから、テーブルを利用した 239 での除算 2 回にしたほうが圧倒的に速くなります。
さらに奇数での除算も、除数が 8bit に収まっているうちはテーブルを利用したほうが速そうです （ただしテーブル作成時間のロスはあります）。 毎回除算テーブルを書き換えるので、 テーブルをもう一組用意することにします。
あとついでに、加減算をループ展開しておきます。1 ループで 8Byte
処理し、端数はループの途中にジャンプすることで調整します
（ジャンプテーブルを使用）。ソースは
mpadd2.z80s です。
ループ展開は確実に速くはなるのですが、 プログラムが汚くなるのが難点ですね。
除算と加減算を改良したのが
pi4.z80s
です。計算時間は 3024 秒
と、そこそこ速くなりました。
除算で 1bit ずつ求める場合、初めに使用した「被除数から除数を引き、0 以上なら商を 1 とし、負になったら除数を足して戻し商を 0 とする」 方法を「回復法」と言います。
これに対し「被除数が 0 以上なら除数を引いて商を 1 とし、被除数が負なら除数を足して商を -1 とする」方法があり、 「非回復法」（または「引きっ放し法」）と言います。
それほど大きな違いではありませんが、16bit 加減算が 1 つ減るので、うまくいけば少し速くなりそうです。
では C 言語で書いてみます。
typedef unsigned char u8;
/* size バイトの src[] を d で割って dest[] に格納する */
void MPDiv(u8 dest[], u8 src[], int d, int size)
{
int i, j;
int r = 0;
for(i = 0; i < size; i++) {
u8 q = 0;
u8 n = src[i];
for(j = 0; j < 8; j++) {
r <<= 1;
if ((n & 0x80) != 0) {
r++;
}
q <<= 1;
if (r >= 0) {
r -= d;
q++;
} else {
r += d;
q--;
}
n <<= 1;
}
/* 1Byte 単位で完結させるための処理 */
if (r < 0) {
r += d;
q--;
}
dest[i] = q;
}
return;
}
8bit 求めた後に被除数が負になっていたら補正処理をしています。 こうしないとバイト単位で処理が完結しません。 その代わり、内側のループの初め（j=0）では必ず r が 0 以上になっています。Z80 版ではこの点を踏まえて処理を省きます。
さて、内側のループ＋後処理を Z80 に置き換えてみます。回復法の 16bit 版と同様に BC レジスタに -d（－除数）をセットしておきます。商の 0 と 1 が逆になるのも同じで、最後に CPL でまとめて反転します。
; HL=余り(r), DE=除数(d), BC=-除数(-d), A=非除数(n):商(q) ADD A,A ADC HL,HL ADD HL,BC ; ここから2bit目 ADD A,A ADC HL,HL ; 符号がCフラグに反映 JR C,L2_1 ADD HL,BC JP L2_2 L2_1: INC A ; 反転しているので、商から1引いていることに相当 ADD HL,DE L2_2: ; ここから3bit目 ADD A,A : ; ここから8bit目とループ後の処理 ADD A,A ADC HL,HL JR C,L8_1 ADD A,A SBC HL,DE ; ADD HL,BCではSフラグが変化しない JP L8_2 L8_1: INC A ADD A,A ADC HL,DE ; ADD HL,DEではSフラグが変化しない L8_2: JP P,L9 ADD HL,DE INC A L9: CPL
商が反転しているせいでわかりにくいですが、
こんな処理になりました。外側のループも付け加えたものが
mpdiv3.z80s
ですが、途中の JP L2_2 などは無駄なので、被除数が 0
以上のときと負のときで処理を分けています。
なお、除数は 15bit までです。 符号付きで処理しているため、16bit に対応するのは難しいと思います。
この除算に置き換えたのが
pi5.z80s
です。計算時間は 2819 秒
と、少しだけ速くなりました。
なんとか除数 9bit 以上の除算もテーブルを利用して高速化できないかと考え、 「だいたいの商を求めて補正」という方法にたどり着きました。 まず除数が 15bit（4000H～7FFFH）の場合で解説します。
1Byte ずつ処理すると、上の桁からの余りが 15bit、これから処理する被除数が 8bit の計 23bit を 15bit で割ることになります。この 23bit を p、15bit（除数）を d とします。
テーブルはせいぜい 8bit で参照したいので、p の上位 8bit で参照することにします。「p の上位 8bit 以外を 0 にした値」を P とし、P を d で割った商を Q とします。そしてテーブルには Q と Q×d を格納しておきます。この Q×d を M とします。
除算処理では、まず p の上位 8bit でテーブル参照して Q と M を得ます。 次に p から M を引き、これを R とします。 Q は商に近い値ですが小さめの値が出るので、補正のために「R が d 以上なら d を引いて Q に 1 を足す」を繰り返します。これで Q が正しい商、R が正しい余りになります。
さて、テーブルには Q と M を格納するわけですが、Q は 8bit でいいとして、M は何ビット必要なのでしょうか。
そこで、被除数と除数の組み合わせを総当たりで調べるプログラムを C 言語で書いて、M が下位 16bit で足りることと、R から d を引く操作は最大で 2 回ということがわかりました。これならテーブルのサイズは除数 8bit と同じ 768Byte で済みます。
除数が 9～14bit の場合は「上位 8bit」のビット位置がずれるだけで、 それ以外は同様です。なお除数 16bit では M が 16bit では足りないので、ここでは除数は 15bit までとします。
それでは C 言語で書いてみます。テーブル参照はせず、毎回計算します。
typedef unsigned char u8;
/* size バイトの src[] を d（0x100～0x7fff）で割って dest[] に格納する */
void MPDiv9_15(u8 dest[], u8 src[], int d, int size)
{
int i;
int s = 0;
int t = d;
while((t & 0x4000) == 0) {
t <<= 1;
s++;
}
int mask = ~(0x7fff >> s);
int r = 0;
for(i = 0; i < size; i++) {
int p = (r << 8) + src[i];
u8 q = (p & mask) / d;
r = (p - q * d) & 0xffff;
if (r >= d) {
r -= d;
q++;
if (r >= d) {
r -= d;
q++;
}
}
dest[i] = q;
}
return;
}
わりとシンプルではないでしょうか。 考え方自体は、除数 8bit のものより単純だと思います。
そして Z80 版が
mpdiv4.z80s
です。MPDIV9_15 が除算ルーチンですが、
除数のビット数によってビット位置調整処理が変わるため、
速度優先ということでビット数ごとに処理を分けています。
テーブル作成ルーチン MPDIV9_15_MAKETBL
は、各ビット数ごとの処理ルーチン内から呼び出しています。
この除算を組み込んだのが
pi6.z80s
です。ただし計算桁数が減ってくるとテーブル作成時間が無視できないので、
有効部分が DIV9_15_MINSIZE Byte
未満になったら、従来の除算ルーチンで処理しています。
そして計算時間は 1933 秒 と、大幅に速くなりました。
そろそろ高速化のネタも尽きてきたので、プログラムを整理します。
まず有効桁数をチェックする CHECK_WORK ルーチンに、呼び出し側で共通している処理を押し込んで、 プログラムサイズを削減します。
それと MPDIV8T ルーチンは WORK2
に対してしか呼び出さないため（WORK3=WORK2÷d 用は
MPDIV8T2）、src=dest を前提にします。速度的にはループ内の
INC DE の 6clk が減るだけですが、
ループ回数が多いのでかなり影響があります。
整理したのが
pi7.z80s
です。計算時間は 1897 秒（31 分 37
秒）で、これが今の自己ベストです。
【'19/05/11 追記】
その後いろいろと改良し、計算時間は
1795 秒（29 分 55 秒）と、ついに 30
分を切ることができました。ソースは
pi7_2.z80s
です。主な改良点は以下の通りです。
さらに 16bit 除算のテーブル版を大幅に高速化できたため、239 で 2 回割っていたところを 2392 で直接割るように戻すことによって 1781 秒（29 分 41 秒）まで縮まりました。 改良版ソースは「もっと桁数を」の追記を参照してください。
【'19/05/25 追記】
さらに改良し、計算時間は 1607 秒（26 分 47 秒）になりました。 改良点と改良版ソースは「もっと桁数を」の追記を参照してください。
9～15bit 除算をテーブル化する前に考えた手法です。 結局必要無くなってソースも残っていませんが、概要だけ解説します。
除数が 8bit×8bit の積に分解できた場合は、8bit 除算を 2 回行います。そのたびにテーブルを作り直しますが、 それでもビット単位で除算するより速いからです。
問題は 8bit×8bit の組み合わせを求める方法です。 まともに素因数分解などしていたら時間がかかりすぎるので、 もっと軽い方法を考えました。
対象の除数は 2n+1 なので、まず 3～255（実際はもっと少なくてよい） の各奇数用のカウンタを用意して 1 で初期化し、n が増えるたびに全てに 2 を足します。各奇数値以上になったらその値を引きます。
ようするに「(2n+1) % 各奇数値」（剰余）を記録します。 この値が 0 なら割り切れるというわけです。
あとはその中の最大値を探し、2n+1 をその値で割って 8bit に収まっていたら分解成功です。
当然ながら素数だと分解できませんし、分解できても 8bit に収まらない場合もありますが、それなりに速くはなりました。
しかしその後 9～15bit 除算もテーブル使用で高速化できたため、この手法は丸ごと消滅しました。 合掌。
15bit に収まらない 16bit 値での除算も、頑張ってテーブル化してみました。
mpdiv5.z80s
です。1 万桁では使用しませんでしたが、この後使用します。
基本的には 9～15bit 除算と同様ですが、16bit に収まらない部分を頑張って補正しています。 かなりややこしくなってしまったので、詳しい解説は省略します。 もっと良い方法があるかもしれません。
【'19/05/11 追記】
除算テーブルに M（Q×d）を 24bit
全て登録することで、テーブルサイズは 768Byte から 1024Byte
に増えましたが、かなり素直で高速な処理になりました。これを
2392 での除算に使用することで、1
万桁の計算も速くなりました。ソースは
mpdiv6.z80s です。
【'19/05/25 追記】
除算テーブルの商に 1 加算することで、テーブルサイズは再び 768Byte で済むようになって速度も速くなりました。
今回は 1 万桁の計算が主な目的でしたが、桁数が増えると
2n+1（奇数）による除算で問題が出てきます。除数が 16bit
を超えると計算が非常に面倒になるので避けたいところです。
計算時間は基本的には桁数の 2 乗に比例しますが、 除数が大きくなりすぎるとさらに遅くなってしまうわけです。
arctan の中の分母で一番小さいものを $ x $、奇数の最大値を $ d $ とすると、求められる桁数は $ d \log_{10} x $ くらいです。マチンの公式（$x=$5）だと奇数 16bit（$d=$65535）で 45806 桁ほどになります。
他の公式では、シュテルマーの公式の一つに、
$$ \frac{\pi}{4} = 6 \arctan \left( \frac{1}{8} \right)
+ 2 \arctan \left(\frac{1}{57} \right)
+ \arctan \left(\frac{1}{239} \right) $$
というのがあります（以後「シュテルマー 1」と表記します）。
この公式だと奇数 16bit（$d=$65535）で
59184 桁ほど求められます。
と言っても、メモリが 64KB だと 5 万桁くらいが限界ですが。
この公式にはコンピュータに都合の良い部分があります。
両辺を 4 倍して第 1 項を展開すると、
$$ 3
- \frac{3}{8^2 \times 3}
+ \frac{3}{8^2 \times 8^2 \times 5}
- \frac{3}{8^2 \times 8^2 \times 8^2 \times 7}
+ \cdots $$
となり、3 から始めて 82
で割っていくわけですが、これは（82=26 なので）右 6bit
シフトするだけで済むため、大幅な時間短縮になります。
もっと分母の大きな公式では、
再びシュテルマーの公式の一つになりますが、
$$ \frac{\pi}{4} = 44 \arctan \left( \frac{1}{57} \right)
+ 7 \arctan \left(\frac{1}{239} \right)
- 12 \arctan \left(\frac{1}{682} \right)
+ 24 \arctan \left(\frac{1}{12943} \right) $$
というのがあります（以後「シュテルマー 2」と表記します）。
これなら奇数 15bit（$d=$32767）でも 57534
桁ほど求められます。
ちなみに 682=2×11×31、12943=7×432 と素因数分解できます。除算を複数回に分けるときの参考までに。
それぞれの公式で、ほぼ限界まで計算してみました。
pi8.z80s です。pi9.z80s です。pi10.z80s です。それぞれの所要時間は以下の通りです。
| 公式＼桁数 | 10000 | 45000 | 50000 |
|---|---|---|---|
| マチン | 1897 | 39102 | - |
| シュテルマー 1 | 2036 | 40060 | 49885 |
| シュテルマー 2 | 2136 | 41771 | 51429 |
4 万 5 千桁まではマチンの公式が、5 万桁ではシュテルマー 1 が速いという結果になりました。シュテルマー 2 は残念ながら一歩及びませんでしたが、計算式が少し複雑な分、 まだ改良の余地があるかもしれません。
【'19/05/11 追記】
その後いろいろと改良しました。 主な改良点は「プログラムの整理」と「他の手法」に追記してあります。
pi8_2.z80spi9_2.z80spi10_2.z80s| 公式＼桁数 | 10000 | 45000 | 50000 |
|---|---|---|---|
| マチン | 1781 | 35146 | - |
| シュテルマー 1 | 1826 | 35742 | 44124 |
| シュテルマー 2 | 1998 | 39702 | 48959 |
シュテルマー 1 とシュテルマー 2 の差が開きました。どうやら 5 万桁ではシュテルマー 2 の出番は無さそうです。
【'19/05/25 追記】
さらに改良しました。主な改良点は以下の通りです。
以下、改良版ソースです。シュテルマー 2 については割愛します。
pi8_3.z80spi9_3.z80s| 公式＼桁数 | 10000 | 45000 | 50000 |
|---|---|---|---|
| マチン | 1618 | 31293 | - |
| シュテルマー 1 | 1607 | 31311 | 38540 |
8bit 除算と 9～16bit 除算の速度差が縮まった影響で、1 万桁ではシュテルマー 1 のほうがマチンより速くなりました。
高速化のためにいろんな手法を考えてみましたが、 最終的にはわりと普通の手法の組み合わせになったかもしれません。まぁ 9bit 以上での除算のテーブル化は、わりと上手くいったとは思います。
昔はプログラムを紙に書いて、 クロック数を数えながら最適化し、ハンドアセンブルして実行して確認、 なんてことをやっていました。
今回は思い付いたことをすぐ反映して実行＆確認できたので、 効率が格段に違いました。そこが昔とは違うところですね。
仕事でこういった高速化の必要が出てくると、 時間が限られているので結構あせりますが、 趣味だと気楽なもので楽しいです。パズルを解く楽しさと似ています。特に Z80 は工夫の余地が多くて飽きません。
とはいえ、Z80 でプログラミングするのは、さすがにこれが最後かもしれませんね。