はじめに
ちょい前ツイッターで
Wozさん、数学もできるのね。自然対数の底e（『博士の愛した公式』に出てくるヤツ）を116,000桁演算するという、しかも、Apple ][ で。アセンブリ言語で書いてあるし。ほとんど読んでないけど、そのうちな（笑）https://t.co/Qt6AX1EaMF
— savasio (@savasio) September 15, 2020
という呟きを見掛けた。
10進数で 116,000桁は 2進数に換算すると
$\frac{116000}{\log_{10}256} = 48167.95737...$
約 47kBということになるので、初期の APPLE][ の全 RAM 48kB の大半を使って 2進数で $e$ を求めた後 10進数に変換して出力するのかなと思ったのだけど、$e$ の計算を単純に
$e=\displaystyle １+\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots$
でやったとしても階乗の逆数とそれを足してく領域が要るからメモリは倍必要だよなあ、どうやってんだ? と不思議に思ったのでした。
以下は個人的なメモなので余程の暇人以外読まなくて良い多分読んだところで現代人の為になるものではない
件の記事は英語版 Wikipedia の『e (mathematical constant)』の注釈にもリンクが挙げられており無名なものではないようです。この記事『スティーブ・ウォズニアック氏 ～ 読んでみた』は件の記事の翻訳ではなく、件の記事自体が解説であり解説の解説に意味はないと思うのでこの記事で解説は行いません。興味が湧いた方は原文に当たって下さい。約40年前の 6502アセンブラ + BASIC による方法は多倍長計算がお手軽に扱える言語がフツーに使える現代では参考になる箇所はあまりないと思いますが、個人的には
- 2つの領域を使わずに $e$ を多倍長計算する方法
- 6502 の命令のオペランドを書き換える高速化テクニック
- 近頃あまり見なくなった 80桁 66行の連続帳票用紙を想定した出力形式
辺りが興味深く面白かったところです。
Cで書き直してみた
件の記事のプログラムの多くの部分は 6502 のアセンブラで書かれています。6502 のプログラムでは通常はスクラッチパッドとしてゼロページを使用するところを件の記事のプログラムでは高速化のために命令のオペランドを直接書き換えるテクニックが多用されていて読み易いものではなかったため、プログラムを理解する補助として C に書き直してみました。
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#define N 9720
#define NUMBYT 0x3800
uint8_t e[NUMBYT];
void ecalc1(void)
{
for (uint16_t n = N; n > 1; n--) {
uint32_t a = 1;
for (int i = 0; i < NUMBYT; i++) {
a = 0x100 * a + e[i];
e[i] = a / n;
a %= n;
}
}
}
uint8_t prodlo[0x100];
uint8_t prodhi[0x100];
void init(void)
{
for (int i = 0; i < 0x100; i++) {
prodlo[i] = 100 * i % 0x100;
prodhi[i] = 100 * i / 0x100;
}
}
int mult(void)
{
int a = 0;
int x = 0;
int c = 0;
for (int i = NUMBYT - 1; i >= 0; i--) {
a = prodhi[x];
x = e[i];
a += prodlo[x] + c;
c = a / 0x100;
e[i] = a % 0x100;
}
return prodhi[x] + c;
}
int main(void)
{
ecalc1();
init();
bool odd = false;
for (int page = 1; page <= 10; page++) {
printf("\n%6sE%63sPAGE %02d\n\n", "", "", page);
for (int line = 1; line <= 60; line++) {
printf("%6s", page == 1 && line == 1 ? "E=2." : "");
for (int group = 1; group <= 12; group++) {
for (int dig = 1; dig <= 5; dig++) {
static int num;
if (!odd) {
num = mult();
}
printf("%d", odd ? num % 10 : num / 10);
odd = !odd;
}
printf(" ");
}
printf("\n");
if (page == 10 && line == 35) {
goto end;
}
}
printf("\n\n\n");
}
end:;
}
件の記事によると APPLE][ 実機では計算に 10時間掛かるとのことですが現代の PC では 1秒やそこらの時間で計算してしまいます。ハードウェアの進歩を実感ぜずにはいられません。
APPLE][がないならMSXを使えばいいじゃない
逆に、$e$ の計算に10時間とか掛ける当時の感覚を体感したいと思い実機で動かしてみたかったのですが、APPLE][ とか持っていないので似たようなものとしてどこのご家庭にもありふれた MSX で動作させてみることとしました。
MAKE.BAT
M80 ECALC1,=ECALC1/Z/R/L
L80 ECALC1,ECALC1.BIN/N/Y/E
M80 EPRNT,=EPRNT/Z/R/L
L80 EPRNT,EPRNT.BIN/N/Y/E
ECALC1.MAC
N EQU 9720
NUMBYT EQU 3800H
EBUF EQU 9000H
ASEG
ORG 0100H
DB 0FEH
DW BEGIN, END, 0
.PHASE 8E00H
BEGIN:
PUBLIC ECALC1
ECALC1: LD DE, -N
NXTDVSR:LD HL, 1
EXX
LD HL, EBUF
LD BC, HIGH NUMBYT
NXTBYTE:LD A, (HL) ;(7+1)
EXX ;(4+1)
LD B, 8 ;(7+1)
RLA ;(4+1)
NXTBIT: ADC HL, HL ;(15+2)
ADD HL, DE ;(11+1)
JR C, $+4 ;(12+1/7+1)
SBC HL, DE ;(15+2)
RLA ;(4+1)
DJNZ NXTBIT ;(13+1/8+1)
EXX ;(4+1)
LD (HL), A ;(7+1)
INC L ;(4+1)
DJNZ NXTBYTE ;(13+1/8+1)
INC H ;(4+1)
DEC C ;(4+1)
JR NZ, NXTBYTE ;(12+1/7+1)
EXX
INC DE
LD A, E
AND D
INC A
JR NZ, NXTDVSR
RET
END EQU $-1
.DEPHASE
END
EPRNT.MAC
RESULT EQU 8FFFH
PRODLO EQU 0D000H
PRODHI EQU 0D100H
N EQU 9720
NUMBYT EQU 3800H
EBUF EQU 9000H
ASEG
ORG 0100H
DB 0FEH
DW BEGIN, END, 0
.PHASE 8F00H
BEGIN:
PUBLIC INIT
INIT: XOR A
LD B, A
LD DE, PRODLO
LD HL, PRODHI
PRODGEN:LD (DE), A
LD (HL), B
ADD A, 100
JR NC, $+3
INC B
INC E
INC L
JR NZ, PRODGEN
RET
PUBLIC MULT
MULT: LD H, HIGH PRODLO
XOR A
LD L, A
LD BC, HIGH NUMBYT
LD DE, EBUF+NUMBYT-1
MULBYTE:EX DE, HL
LD E, (HL)
EX DE, HL
ADC A, (HL)
LD (DE), A
INC H
LD A, (HL)
DEC H
DEC E
DJNZ MULBYTE
DEC D
DEC C
JR NZ, MULBYTE
ADC A, 0
LD (RESULT), A
RET
END EQU $-1
.DEPHASE
END
ビルドは MSX-DOS 上で M80 と L80 を使用しています。MAKE.BAT を実行すると BLOAD 形式の ECALC1.BIN と EPRNT.BIN、及び両プリントファイルが生成されます。
ECALC1.PRN
MSX.M-80 2.00 04-Oct-20 PAGE 1
25F8 N EQU 9720
3800 NUMBYT EQU 3800H
9000 EBUF EQU 9000H
0000' ASEG
ORG 0100H
0100 FE DB 0FEH
0101 8E00 8E2C DW BEGIN, END, 0
0105 0000
.PHASE 8E00H
8E00 BEGIN:
PUBLIC ECALC1
8E00 11 DA08 ECALC1: LD DE, -N
8E03 21 0001 NXTDVSR:LD HL, 1
8E06 D9 EXX
8E07 21 9000 LD HL, EBUF
8E0A 01 0038 LD BC, HIGH NUMBYT
8E0D 7E NXTBYTE:LD A, (HL) ;(7+1)
8E0E D9 EXX ;(4+1)
8E0F 06 08 LD B, 8 ;(7+1)
8E11 17 RLA ;(4+1)
8E12 ED 6A NXTBIT: ADC HL, HL ;(15+2)
8E14 19 ADD HL, DE ;(11+1)
8E15 38 02 JR C, $+4 ;(12+1/7+1)
8E17 ED 52 SBC HL, DE ;(15+2)
8E19 17 RLA ;(4+1)
8E1A 10 F6 DJNZ NXTBIT ;(13+1/8+1)
8E1C D9 EXX ;(4+1)
8E1D 77 LD (HL), A ;(7+1)
8E1E 2C INC L ;(4+1)
8E1F 10 EC DJNZ NXTBYTE ;(13+1/8+1)
8E21 24 INC H ;(4+1)
8E22 0D DEC C ;(4+1)
8E23 20 E8 JR NZ, NXTBYTE ;(12+1/7+1)
8E25 D9 EXX
8E26 13 INC DE
8E27 7B LD A, E
8E28 A2 AND D
8E29 3C INC A
8E2A 20 D7 JR NZ, NXTDVSR
8E2C C9 RET
8E2C END EQU $-1
.DEPHASE
END
MSX.M-80 2.00 04-Oct-20 PAGE S
Macros:
Symbols:
8E00 BEGIN 9000 EBUF 8E00I ECALC1
8E2C END 25F8 N 3800 NUMBYT
8E12 NXTBIT 8E0D NXTBYTE 8E03 NXTDVSR
No Fatal error(s)
EPRNT.PRN
MSX.M-80 2.00 04-Oct-20 PAGE 1
8FFF RESULT EQU 8FFFH
D000 PRODLO EQU 0D000H
D100 PRODHI EQU 0D100H
25F8 N EQU 9720
3800 NUMBYT EQU 3800H
9000 EBUF EQU 9000H
0000' ASEG
ORG 0100H
0100 FE DB 0FEH
0101 8F00 8F32 DW BEGIN, END, 0
0105 0000
.PHASE 8F00H
8F00 BEGIN:
PUBLIC INIT
8F00 AF INIT: XOR A
8F01 47 LD B, A
8F02 11 D000 LD DE, PRODLO
8F05 21 D100 LD HL, PRODHI
8F08 12 PRODGEN:LD (DE), A
8F09 70 LD (HL), B
8F0A C6 64 ADD A, 100
8F0C 30 01 JR NC, $+3
8F0E 04 INC B
8F0F 1C INC E
8F10 2C INC L
8F11 20 F5 JR NZ, PRODGEN
8F13 C9 RET
PUBLIC MULT
8F14 26 D0 MULT: LD H, HIGH PRODLO
8F16 AF XOR A
8F17 6F LD L, A
8F18 01 0038 LD BC, HIGH NUMBYT
8F1B 11 C7FF LD DE, EBUF+NUMBYT-1
8F1E EB MULBYTE:EX DE, HL
8F1F 5E LD E, (HL)
8F20 EB EX DE, HL
8F21 8E ADC A, (HL)
8F22 12 LD (DE), A
8F23 24 INC H
8F24 7E LD A, (HL)
8F25 25 DEC H
8F26 1D DEC E
8F27 10 F5 DJNZ MULBYTE
8F29 15 DEC D
8F2A 0D DEC C
8F2B 20 F1 JR NZ, MULBYTE
8F2D CE 00 ADC A, 0
8F2F 32 8FFF LD (RESULT), A
8F32 C9 RET
8F32 END EQU $-1
.DEPHASE
END
MSX.M-80 2.00 04-Oct-20 PAGE S
Macros:
Symbols:
8F00 BEGIN 9000 EBUF 8F32 END
8F00I INIT 8F1E MULBYTE 8F14I MULT
25F8 N 3800 NUMBYT 8F08 PRODGEN
D100 PRODHI D000 PRODLO 8FFF RESULT
No Fatal error(s)
DISK-BASIC 上で
ECALC.BAS
10 CLEAR 300,&H8DFF
20 BLOAD "ECALC1.BIN"
30 DEFUSR=&H8E00
40 A=USR(0)
50 BSAVE "E1.BIN",&H9000,&HC7FF
60 END
を実行すると $e$ の小数点以下を 2進数で計算し、結果が E1.BIN のファイル名で保存されます。
計算の前後に開始と終了の時刻の出力を加えて計算時間を計ったところでは MSXPLAYer の MSX2+ モードで動作させて 7時間7分14秒ほどとなりました。
計算時間のほとんどを占めているであろう ECALC1.MAC の除算部分のサイクル数の合計を MSX2+ の CPU 速度(Z80 3.579545MHz M1サイクル 1ウェイト)に換算すると少し遅い感じですが、割り込み処理の負荷を含めるとこんなものかもしれません。
E1.BIN が保存されてる状態で
FORMAT.BAS
10 CLEAR 300, &H8EFF: DEFINT A-Z: BLOAD "E1.BIN": BLOAD "EPRNT.BIN": OPEN "CRT:" FOR OUTPUT AS #1
20 DEFUSR0=&H8F00: DEFUSR1=&H8F14: A=USR0(0): ODDEVEN=0
30 FOR PAGE=1 TO 10: PRINT #1, : PRINT #1, " E"; SPC(63); "PAGE "; HEX$(PAGE/10); HEX$(PAGE MOD 10): PRINT #1,
40 FOR LI=1 TO 60: IF PAGE>1 OR LI>1 THEN 50 ELSE PRINT #1, " E=2.";: GOTO 60
50 PRINT #1, " ";
60 FOR GROUP=1 TO 12
70 FOR DIG=1 TO 5: GOSUB 200: NEXT DIG
80 PRINT #1, " ";: NEXT GROUP
90 PRINT #1, : IF PAGE=10 AND LI=35 THEN 110 ELSE NEXT LI: REM QUIT AFTER 34500 DIGITS
100 PRINT #1, : PRINT #1, : PRINT #1, : NEXT PAGE
110 CLOSE #1: END
190 REM
192 REM SUBROUTINE 200 PRINTS NEXT DIG
194 REM
200 IF ODDEVEN=1 THEN 220 ELSE A=USR1(0)
210 PRINT #1, HEX$(PEEK(&H8FFF)/10); : GOTO 230
220 PRINT #1, HEX$(PEEK(&H8FFF) MOD 10);
230 ODDEVEN=1-ODDEVEN: RETURN
を実行すると MSX の画面に $e$ が 10進数で出力されます。10行目の "CRT:" を "LPT:" に変更することでプリンタへの出力となることを期待していますが動作を確認しておりません。
おわりに
おわりです。
2021/04/22 追記
「Cで書き直してみた」版が原作に倣い 1バイト単位でメモリを読み込み、32bit÷16bit の除算を行っていたものを、現代風に 4バイト単位でメモリを読み込み 64bit÷32bit の除算を行うよう改変してみたので記録として残しておく。
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#define N 9720
#define NUMBYT 0x3800
uint32_t e[NUMBYT / sizeof(uint32_t)];
void ecalc1(void)
{
for (uint32_t n = N; n > 1; n--) {
uint64_t a = 1;
for (int i = 0; i < NUMBYT / sizeof(uint32_t); i++) {
a = (a << 32) + e[i];
e[i] = a / n;
a %= n;
}
}
}
int mult(void)
{
uint64_t a = 0;
for (int i = NUMBYT / sizeof(uint32_t) - 1; i >= 0; i--) {
a += 100000 * (uint64_t)e[i];
e[i] = a;
a >>= 32;
}
return a;
}
int main(void)
{
ecalc1();
for (int page = 1; page <= 10; page++) {
printf("\n%6sE%63sPAGE %02d\n\n", "", "", page);
for (int line = 1; line <= 60; line++) {
printf("%6s", page == 1 && line == 1 ? "E=2." : "");
for (int group = 1; group <= 12; group++) {
printf("%05d ", mult());
}
printf("\n");
if (page == 10 && line == 35) {
goto end;
}
}
printf("\n\n\n");
}
end:;
}
除算の回数が減った分高速化され、処理時間が 30% 程度となった。
更に、インラインアセンブラを試してみた。4バイト読み込み 64bit÷32bit の除算を繰り返している点は変わらない。
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#define N 9720
#define NUMBYT 0x3800
uint32_t e[NUMBYT / sizeof(uint32_t)];
void ecalc1(void)
{
__asm __volatile(
" lea %0, %%rcx \n"
" mov %1, %%esi \n"
"0: xor %%edi, %%edi \n"
" mov $1, %%edx \n"
"1: mov (%%rcx, %%rdi), %%eax \n"
" add $4, %%rdi \n"
" div %%esi \n"
" mov %%eax, -4(%%rcx, %%rdi) \n"
" cmp %2, %%rdi \n"
" jne 1b \n"
" sub $1, %%esi \n"
" cmp $1, %%esi \n"
" jne 0b \n"
:: "m"(e), "i"(N), "i"(NUMBYT): "rax", "rcx", "rdx", "rdi", "rsi", "memory", "cc"
);
}
int mult(void)
{
int ret;
__asm(
" xor %%ecx, %%ecx \n"
" mov $100000, %%r8d \n"
" lea %1, %%rdi \n"
" mov %2/4-1, %%esi \n"
"1: mov (%%rdi, %%rsi, 4), %%eax \n"
" mul %%r8d \n"
" add %%ecx, %%eax \n"
" mov %%eax, (%%rdi, %%rsi, 4) \n"
" adc $0, %%edx \n"
" mov %%edx, %%ecx \n"
" sub $1, %%esi \n"
" jns 1b \n"
: "=c"(ret): "m"(e), "i"(NUMBYT): "rax", "rdx", "rdi", "rsi", "r8", "memory", "cc"
);
return ret;
}
int main(void)
{
ecalc1();
for (int page = 1; page <= 10; page++) {
printf("\n%6sE%63sPAGE %02d\n\n", "", "", page);
for (int line = 1; line <= 60; line++) {
printf("%6s", page == 1 && line == 1 ? "E=2." : "");
for (int group = 1; group <= 12; group++) {
printf("%05d ", mult());
}
printf("\n");
if (page == 10 && line == 35) {
goto end;
}
}
printf("\n\n\n");
}
end:;
}
更に 3割程度? 速くなった(気がする)。
読み込む単位を 8バイトに変更し 128bit÷64bit の除算を繰り返すよう変更することで除算の回数を更に 1/2 に減らしてみた。
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#define N 9720
#define NUMBYT 0x3800
uint64_t e[NUMBYT / sizeof(uint64_t)];
void ecalc1(void)
{
__asm __volatile(
" lea %0, %%rcx \n"
" mov %1, %%esi \n"
"0: xor %%edi, %%edi \n"
" mov $1, %%edx \n"
"1: mov (%%rcx, %%rdi), %%rax \n"
" add $8, %%edi \n"
" div %%rsi \n"
" mov %%rax, -8(%%rcx, %%rdi) \n"
" cmp %2, %%edi \n"
" jne 1b \n"
" sub $1, %%esi \n"
" cmp $1, %%esi \n"
" jne 0b \n"
:: "m"(e), "i"(N), "i"(NUMBYT): "rax", "rcx", "rdx", "rdi", "rsi", "memory", "cc" );
}
int mult(void)
{
int ret;
__asm(
" xor %%ecx, %%ecx \n"
" mov $100000, %%r8d \n"
" lea %1, %%rdi \n"
" mov %2/8-1, %%esi \n"
"1: mov (%%rdi, %%rsi, 8), %%rax \n"
" mul %%r8 \n"
" add %%rcx, %%rax \n"
" mov %%rax, (%%rdi, %%rsi, 8) \n"
" adc $0, %%rdx \n"
" mov %%rdx, %%rcx \n"
" sub $1, %%esi \n"
" jns 1b \n"
: "=c"(ret): "m"(e), "i"(NUMBYT): "rax", "rdx", "rdi", "rsi", "r8", "memory", "cc"
);
return ret;
}
int main(void)
{
ecalc1();
for (int page = 1; page <= 10; page++) {
printf("\n%6sE%63sPAGE %02d\n\n", "", "", page);
for (int line = 1; line <= 60; line++) {
printf("%6s", page == 1 && line == 1 ? "E=2." : "");
for (int group = 1; group <= 12; group++) {
printf("%05d ", mult());
}
printf("\n");
if (page == 10 && line == 35) {
goto end;
}
}
printf("\n\n\n");
}
end:;
}
なんだか倍程度遅くなった。128bit÷64bit の除算が重いのか、原因は未確認(その内気が向いたら確認する)。