はじめに
この記事は @amidoazuki さんの投稿『36年前の沢口靖子さんのプロフィールを確認した話』のなかで紹介されていた 36年前の雑誌記事にインスパイアされたものです。
記事中の記述
を目にすれば、「それじゃ整数型 BASIC で動かしてみようか」とは誰もが考えることですね。
整数型 BASICとは
1970年代後半のマイコンブームの時代に一時流行った Tiny BASIC に端を発する小型の BASIC言語の処理系で、フットプリントの小ささと実装の簡便さを重視して BASIC言語の仕様から機能を絞り実装がされていて扱える数値も浮動小数点型をサポートせず整数演算のみとなっているものの総称として俗に「整数型 BASIC」という呼び方がされていました。この場合の整数演算とは現代の Ruby や Python のように演算結果に応じて勝手にサイズが伸長されるようなものではなく 16bit 固定が一般的でした。
あまり実用的なものではなかったことと、ある程度のメモリ容量を持ったメーカー製の高機能な BASIC言語が最初から使えるパソコンが当たり前となった以降では整数型 BASIC の利点も薄くなったので、件の雑誌記事が掲載された 1984年の時点では安価なホビー用パソコンや一部のホビー向けの処理系を除けば整数型 BASIC は既に時代遅れなものとなっていました。
動作環境
今回は整数型 BASIC の動作環境として Ichigojam を使用しています。
PALO ALTO TINY BASIC や FAMILY BASIC、APPLEII 6K BASIC、m5 BASIC-I、JR-100 等も候補として考えたのですが、現代となっては動作環境を用意するのも面倒であり、対して web 上の実行環境が用意されている Ichigojam はお手軽で良いのでこれを使用することとしました。
コード
Ichigojam 用に実装したコードは以下となっております。Ichigojam は BASIC プログラムの最大容量が 1KB と控えめな仕様なのでそれに則り 3分割しています。
100 P=#C02:O=2
110 PRINT "ｾｲﾈﾝｶﾞｯﾋﾟ : ";:GOSUB 1000
120 PRINT "ｼｭｯｼﾝﾁ : ";:GOSUB 1000
130 PRINT "ｼﾝﾁｮｳ,ﾀｲｼﾞｭｳ,B･W･H ﾉ ｻｲｽﾞ : ";:GOSUB 1000
140 PRINT "ｹﾂｴｷｶﾞﾀ,ｾｲｻﾞ : ";:GOSUB 1000
150 PRINT "ｶｿﾞｸｺｳｾｲ : ";:GOSUB 1000
160 PRINT "ｼｭﾐ : ";:GOSUB 1000
170 PRINT "ｽｷﾅ ｽﾎﾟｰﾂ : ";:GOSUB 1000
180 PRINT "ﾄｸｲﾅ ｶﾓｸ : ";:GOSUB 1000
270 END
1000 A=0:Q=P+PEEK(P)+4:IFPEEK(P+1)-39P=Q:GOTO1000
1005 C=PEEK(P+O):O=O+1:IFC>47A=10*A+C-48:GOTO1005ELSEIF!CP=Q:O=2
1010 B=(56*A+7*A/10-A/14*A/307*A/8)/108
1030 IF B=35 THEN PRINT:RETURN
1040 PRINT CHR$(B);
1050 GOTO 1000
2000 '161,100,93,89,104,89,95,95,67,372,372,386,374,62,429,62,386,374,364
2010 '62,389,67,95,102,110,193,214,62,100,102,209,201,62,62,108,93,89,102
2020 '108,89,108,106,67,126,62,374,478,398,62,429,398,384,478,62,386,478,67
2030 '401,401,62,424,424,62,362,413,62,421,62,100,413,475,62,374,396,478
2040 '379,67,372,475,374,478,379,374,475,389,355,367,62,369,364,374,478,374
2050 '475,389,355,367,67,406,413,391,62,391,364,369,364,67,369,364,384,478
2060 '62,374,478,62,391,376,62,406,478,389,398,67
100 P=#C02:O=2
190 PRINT "ｽｷﾅ ﾀﾍﾞﾓﾉ : ";:GOSUB 1000
200 PRINT "ｲｯﾃﾐﾀｲ ｸﾆ : ";:GOSUB 1000
210 PRINT "ｽｷﾅ ﾀﾚﾝﾄ ﾐｭｰｼﾞｼｬﾝ : ";:GOSUB 1000
220 PRINT "ﾘｿｳﾉ ﾀｲﾌﾟﾉ ﾀﾞﾝｾｲ : ";:GOSUB 1000
230 PRINT "ｷﾗｲﾅ ﾀｲﾌﾟﾉ ﾀﾞﾝｾｲ : ";:GOSUB 1000
270 END
1000 A=0:Q=P+PEEK(P)+4:IFPEEK(P+1)-39P=Q:GOTO1000
1005 C=PEEK(P+O):O=O+1:IFC>47A=10*A+C-48:GOTO1005ELSEIF!CP=Q:O=2
1010 B=(56*A+7*A/10-A/14*A/307*A/8)/108
1030 IF B=35 THEN PRINT:RETURN
1040 PRINT CHR$(B);
1050 GOTO 1000
2060 '386,458,398,478,62,464
2070 '364,62,411,458,64,62,456,360,379,478,464,408,62,447,67,391,364,391,62
2080 '77,424,364,389,478,62,374,478,62,364,398,456,367,411,62,408,384,469
2090 '79,67,453,360,439,475,62,362,475,461,62,372,429,384,360,391,67,401
2100 '355,357,426,481,461,62,426,355,367,376,475,62,406,478,62,364,357,389
2110 '355,413,364,406,62,398,421,389,364,62,426,408,62,374,478,62,364,364
2120 '411,67,389,403,384,379,406,62,389,478,389,475,62,374,389,478,355,367
2130 '62,421,62,426,408,62,374,478,62,364,450,67
100 P=#C02:O=2
240 PRINT "ﾃﾞﾋﾞｭｰﾉ ｷｯｶｹ : ";:GOSUB 1000
250 PRINT "ﾄﾞﾝﾅ ｵﾝﾅﾉﾋﾄﾆ ﾅｯﾃｲｷﾀｲｶ : ";:GOSUB 1000
260 PRINT "ｽｷﾅ ｲﾛ : ";:GOSUB 1000
270 END
1000 A=0:Q=P+PEEK(P)+4:IFPEEK(P+1)-39P=Q:GOTO1000
1005 C=PEEK(P+O):O=O+1:IFC>47A=10*A+C-48:GOTO1005ELSEIF!CP=Q:O=2
1010 B=(56*A+7*A/10-A/14*A/307*A/8)/108
1030 IF B=35 THEN PRINT:RETURN
1040 PRINT CHR$(B);
1050 GOTO 1000
2130 '398,478,364,95,374,364,62
2140 '408,367,434,367,62,389,475,406,478,467,458,62,421,62,372,360,406,478
2150 '341,389,355,475,62,406,478,62,453,367,389,355,367,62,389,437,389,398
2160 '67,393,364,374,379,62,374,458,62,391,432,478,406,62,374,472,364,364
2170 '62,408,62,364,472,467,464,62,372,475,411,62,421,62,426,408,62,413,62
2180 '411,461,398,364,67,389,469,62,426,481,475,379,67,133,151,126,159,141
解説
元の記事のプログラムのプロフィール表示のサブルーチンは
1000 READ A
1010 B=INT(300*SIN(A/10*3.1416/180))
1020 C$=CHR$(B)
1030 IF C$="#" THEN PRINT:RETURN
1040 PRINT C$;
1050 GOTO 1000
となっています。1000行で変数 A に数値を読み込み、1010行で変数 A の値をラジアンに変換して SIN 関数を呼び、SIN 関数が返した結果を 300倍し、INT 関数で小数点以下を切り捨てた結果を文字コードとして変数 B に格納し、1020行以降で表示する、という流れです。変数 A に読み込まれる数値は 10 で割って円周率の近似値を掛け、180 で割って SIN 関数の引数としており 0.1度を単位とした値となっています。
整数型 BASIC にこのプログラムを移植するとして問題となるのが SIN 関数です。多くの整数型 BASIC は小数点以下の値を扱う算術関数をサポートしておりません。Ichigojam は例外的に SIN 関数をサポートしていますが、1度を単位とした入力に対し半径 256 の円に対するサイン値を返す仕様であり今回の目的には精度が足りません。先の 1010行の内容を Ichigojam BASIC の仕様に合わせると
1010 B=300*SIN(A/10)/256
と書き換えられそうですが、最大 256 を返す SIN 関数の結果を 300倍すると 16bit 整数の範囲を超えてオーバーフローしてしまうため
1010 B=(300/4)*SIN(A/10)/(256/4)
と書き換えます。これを実行した結果が以下となります。
Ichigojam webで実行
Ichigojam で今回の目的に必要な精度を得るには別のアプローチが必要となります。
INT(300*SIN(A/10*3.1416/180)) の近似式の実装
ラジアンの値 x を引数に持つ正弦関数 sin(x) はテイラー展開すると
sin(x)=\sum_{n=0}^∞ \frac{(-1)^n}{(2n+1)!}x^{2n+1}
となります。総和記号を使わないで表記すると
sin(x)=x-\frac{x^3}{6}+\frac{x^5}{120}-\frac{x^7}{5040}+\frac{x^9}{362880}-\frac{x^{11}}{39916800}+\frac{x^{13}}{6227020800} ・・・
こんな感じです。
無限に繰り返し計算は行えないので今回の目的について第何項まで求めれば十分な精度が得られそうかを次の C プログラムで求めてみました。
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
unsigned long fact(unsigned long x)
{
return x <= 1 ? 1 : x * fact(x - 1);
}
double mysin(double x, int m)
{
double s = 0.0;
for (int n = 0; n < m; n++) {
s += pow(-1, n) / fact(2 * n + 1) * pow(x, 2 * n + 1);
}
return s;
}
bool test(int n)
{
static const int t[] = {
62, 64, 67, 77, 79, 89, 93, 95, 100, 102,
104, 106, 108, 110, 126, 133, 141, 151, 159, 161,
193, 201, 209, 214, 341, 355, 357, 360, 362, 364,
367, 369, 372, 374, 376, 379, 384, 386, 389, 391,
393, 396, 398, 401, 403, 406, 408, 411, 413, 421,
424, 426, 429, 432, 434, 437, 439, 447, 450, 453,
456, 458, 461, 464, 467, 469, 472, 475, 478, 481,
};
for (int i = 0; i < (int)(sizeof(t) / sizeof(t[0])); i++) {
int a = t[i];
double r = a / 10.0 * 3.1416 / 180.0;
int b = (int)(300.0 * sin(r));
int c = (int)(300.0 * mysin(r, n));
if (b != c) {
return false;
}
}
return true;
}
int main(void)
{
for (int i = 1; i <= 9; i++) {
printf("%2d: %s\n", i, test(i) ? "OK" : "Err");
}
}
オリジナルのプログラムで使用されている 62 から 481 の範囲の暗号化された数値 70種類すべてについてライブラリの sin() 関数と独自実装した mysin() 関数とで復調した結果が同じとなるかを判定し、mysin() 関数の実装に第何項までの計算が必要かを求めています。
このプログラムを実行すると
1: Err
2: Err
3: OK
4: OK
5: OK
6: OK
7: OK
8: OK
9: OK
とりあえずのところ浮動小数点演算では第 3項まで求めれば良いらしいことがわかりました。要するに
sin(x)≒x-\frac{x^3}{6}+\frac{x^5}{120}
で良いということです。
あと、オリジナルのプログラムで使用している円周率の近似値 3.1416 はあまり正確な値ではありませんが、今回の目的には円周率の近似値はどの程度の範囲で不正確でも良いかを確認します。
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
double mysin(double x)
{
return x - pow(x, 3) / 6.0 + pow(x, 5) / 120.0;
}
bool test(double pi)
{
static const int t[] = {
62, 64, 67, 77, 79, 89, 93, 95, 100, 102,
104, 106, 108, 110, 126, 133, 141, 151, 159, 161,
193, 201, 209, 214, 341, 355, 357, 360, 362, 364,
367, 369, 372, 374, 376, 379, 384, 386, 389, 391,
393, 396, 398, 401, 403, 406, 408, 411, 413, 421,
424, 426, 429, 432, 434, 437, 439, 447, 450, 453,
456, 458, 461, 464, 467, 469, 472, 475, 478, 481,
};
for (int i = 0; i < (int)(sizeof(t) / sizeof(t[0])); i++) {
int a = t[i];
int b = (int)(300.0 * sin(a / 10.0 * 3.1416 / 180.0));
double r = a / 10.0 * pi / 180.0;
int c = (int)(300.0 * mysin(r));
if (b != c) {
return false;
}
}
return true;
}
int main(void)
{
for (double pi = 3.1000; pi <= 3.2000; pi += 0.0001) {
printf("%f: %s\n", pi, test(pi) ? "OK" : "Err");
}
}
暗号化した文字がすべて正しく復調される円周率の近似値を 3.1000 から 3.2000 の範囲で 0.0001刻みで調べます。
このプログラムを実行すると
～ 略 ～
3.141400: Err
3.141500: OK
3.141600: OK
3.141700: OK
～ 略 ～
3.153300: OK
3.153400: OK
3.153500: OK
3.153600: Err
～ 略 ～
3.1415 から 3.1535 の範囲であれば今回の目的には使用できそうな感じです。案外範囲が広いので、計算に都合の良い値を選べそうです。今回は円周率の近似値として 10進数で切りの良い数字である 3.15 を使用することにします。
さて、先ほど求めた正弦関数の近似式
sin(x)≒x-\frac{x^3}{6}+\frac{x^5}{120}
をオリジナルのプログラムの
1010 B=INT(300*SIN(A/10*3.1416/180))
に当てはめると
1010 X=A/10*3.1416/180:B=INT(300*(X-X^3/6+X^5/120))
となります。これを整数演算化できれば Ichigojam で 36年前の沢口靖子さんのプロフィールが確認できそうです。
円周率の近似値を PI に置き換え X を展開すると
1010 B=INT(300*(A*PI/1800-(A*PI/1800)^3/6+(A*PI/1800)^5/120))
となります。
第1～3項を加えた後で 300倍するのを各項を加えた後で 300倍するよう式を変更します。
1010 B=INT(300*A*PI/1800-300*(A*PI/1800)^3/6+300*(A*PI/1800)^5/120)
整数演算でも小数点以下の数を扱いたいのでゲタを履かせます。具体的には各項を 108倍し、加算した後で 108 で割る様変更します。108 は元の 300倍に掛けても 16bit の整数の上限である 32767 を超えない範囲の値で適当に切りの良い数字なので採用しています。
1010 B=(INT(32400*A*PI/1800)-INT(32400*(A*PI/1800)^3/6)+INT(32400*(A*PI/1800)^5/120))/108
各項を約分します。
1010 B=(INT(18*A*PI)-INT(A^3/(1080000/PI^3))+INT(A^5/(69984000000000/PI^5)))/108
PI は 3.15 を使用すると決めたのでこの値で展開すると
1010 B=(INT(56.7*A)-INT(A^3/34553.50)+INT(A^5/225655535942))/108
大体こんな感じになります。
一応 C プログラムで、暗号化された 70種類の数値について複合結果が一致するか確認しておきましょう。
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
bool test(void)
{
static const int t[] = {
62, 64, 67, 77, 79, 89, 93, 95, 100, 102,
104, 106, 108, 110, 126, 133, 141, 151, 159, 161,
193, 201, 209, 214, 341, 355, 357, 360, 362, 364,
367, 369, 372, 374, 376, 379, 384, 386, 389, 391,
393, 396, 398, 401, 403, 406, 408, 411, 413, 421,
424, 426, 429, 432, 434, 437, 439, 447, 450, 453,
456, 458, 461, 464, 467, 469, 472, 475, 478, 481,
};
for (int i = 0; i < (int)(sizeof(t) / sizeof(t[0])); i++) {
int a = t[i];
int b = (int)(300.0 * sin(a / 10.0 * 3.1416 / 180.0));
int c = ((int)(56.7 * a) - (int)(pow(a, 3) / 34553.50) + (int)(pow(a, 5) / 225655535942.0)) / 108;
if (b != c) {
return false;
}
}
return true;
}
int main(void)
{
printf("%s\n", test() ? "OK" : "Err");
}
OK
問題ないようです。
さて、準備が揃ったところで整数演算化なのですが、
1010 B=(INT(56.7*A)-INT(A^3/34553.50)+INT(A^5/225655535942))/108
の第1項、INT(56.7*A)は56*A+7*A/10とすることで簡単に整数演算化できます。除算の際に切り捨てが発生しますが 108倍のゲタを履かせているので問題はないものとします。
第2項のINT(A^3/34553.50)は、A^3 だけで簡単に 16bit 整数の範囲を超えてオーバーフローしてしまうため少し工夫が必要となります。
A/除数1*A/除数2*A/INT(34553.50/除数1/除数2)
としてオーバーフローしない様少しづつに分けて割ることを考えてみました。除数1と除数2はそれぞれ整数の適切な値を選ぶ必要があります。
第3項のINT(A^5/225655535942)は、A の値が取り得る最大の 471 だとしてもさほど大きな値とならないため、多少計算に誤差が出ても許されそうです。A^5 はすぐにオーバーフローとなってしまうため第2項と同じく少しづつ割る必要はありますが、取り敢えずのところ 225655535942 の 1/5乗に近い値である 186 で割ってみることとします。
A/186*A/186*A/186*A/186*A/186
第2項の除数1と除数2にどのような値を使用すれば良いのか見当がつかないので頭が悪い方法ですが総当たりで調べてみました。
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
const double pi = 3.15;
bool test(int d1, int d2, int d3)
{
static const int t[] = {
62, 64, 67, 77, 79, 89, 93, 95, 100, 102,
104, 106, 108, 110, 126, 133, 141, 151, 159, 161,
193, 201, 209, 214, 341, 355, 357, 360, 362, 364,
367, 369, 372, 374, 376, 379, 384, 386, 389, 391,
393, 396, 398, 401, 403, 406, 408, 411, 413, 421,
424, 426, 429, 432, 434, 437, 439, 447, 450, 453,
456, 458, 461, 464, 467, 469, 472, 475, 478, 481,
};
if (d3 <= 0) return false;
for (int i = 0; i < (int)(sizeof(t) / sizeof(t[0])); i++) {
int a = t[i];
int b = (int)(300.0 * sin(a / 10.0 * 3.1416 / 180.0));
if (a / d1 * a > 32767) return false;
if (a / d1 * a / d2 * a > 32767) return false;
int c = ((56 * a + 7 * a / 10) - (a / d1 * a / d2 * a / d3) + (a / 186 * a / 186 * a / 186 * a / 186 * a / 186)) / 108;
if (b != c) {
return false;
}
}
return true;
}
int main(void)
{
for (int d1 = 1; d1 <= 482; d1++) {
for (int d2 = 1; d2 <= 482; d2++) {
int d3 = (int)(1080000.0 / pow(pi, 3) / d1 / d2);
if (test(d1, d2, d3)) {
printf("%d %d %d\n", d1, d2, d3);
}
}
}
}
8 473 9
8 474 9
8 475 9
9 376 10
9 377 10
9 378 10
9 421 9
9 422 9
9 471 8
9 472 8
9 473 8
9 474 8
10 379 9
浮動小数点演算と sin()関数を使用した場合とで同じ結果が得られる組み合わせを複数確認することができました。この内の最初の結果を整数型BASIC のプログラムに実装すると
1010 B=(56*A+7*A/10-A/8*A/473*A/9+A/186*A/186*A/186*A/186*A/186)/108
ということになります。
Ichigojam webで実行
さて、以上で当初の目的である整数型BASICを使用して 36年前の沢口靖子さんのプロフィールを確認することは達成できたのですが、折角なのでもう少し詰めてみましょう。
先ほど第3項をA/186*A/186*A/186*A/186*A/186としたのは「多少計算に誤差が出ても許されそう」という判断からですが、逆に「誤差が許されるのならばもっと簡単な計算でも良いのでは?」という気がします。いっそのことA/除数4くらいに簡略化できないものでしょうか? 折角なので先の C プログラムを小変更し総当たりで確かめることにしました。
#include <stdio.h>
#include <math.h>
#include <stdbool.h>
const double pi = 3.15;
bool test(int d1, int d2, int d3, int d4)
{
static const int t[] = {
62, 64, 67, 77, 79, 89, 93, 95, 100, 102,
104, 106, 108, 110, 126, 133, 141, 151, 159, 161,
193, 201, 209, 214, 341, 355, 357, 360, 362, 364,
367, 369, 372, 374, 376, 379, 384, 386, 389, 391,
393, 396, 398, 401, 403, 406, 408, 411, 413, 421,
424, 426, 429, 432, 434, 437, 439, 447, 450, 453,
456, 458, 461, 464, 467, 469, 472, 475, 478, 481,
};
if (d3 <= 0) return false;
for (int i = 0; i < (int)(sizeof(t) / sizeof(t[0])); i++) {
int a = t[i];
int b = (int)(300.0 * sin(a / 10.0 * 3.1416 / 180.0));
if (a / d1 * a > 32767) return false;
if (a / d1 * a / d2 * a > 32767) return false;
int c = ((56 * a + 7 * a / 10) - (a / d1 * a / d2 * a / d3) + (a / d4)) / 108;
if (b != c) {
return false;
}
}
return true;
}
int main(void)
{
for (int d1 = 1; d1 <= 482; d1++) {
for (int d2 = 1; d2 <= 482; d2++) {
int d3 = (int)(1080000.0 / pow(pi, 3) / d1 / d2);
for (int d4 = 1; d4 <= 482; d4++) {
if (test(d1, d2, d3, d4)) {
printf("%d %d %d %d\n", d1, d2, d3, d4);
}
}
}
}
}
～ 略 ～
14 307 8 482
～ 略 ～
14 308 8 482
除数4 として暗号化された数値の最大値である 481 より大きい値が使える組み合わせが 2つ見つかりました。
これは、第3項が常に 0 となっても正しい答えが求められるということなので、要は第3項は不要となります。整数型BASIC のプログラムに実装すると
1010 B=(56*A+7*A/10-A/14*A/307*A/8)/108
ということとなります。
浮動小数点演算では第3項まで求める必要があったものが、より計算精度の低い 16bit の整数演算で第2項までの計算で求まってしまうのは理不尽な気はしますが、除算の際の切り捨てやら 108 倍のゲタ履きやらの結果として「たまたま」良い感じのパタンが見つかったのだと思います。
READ～DATA の代用
オリジナルのプログラムでは暗号化した数値を BASIC の DATA文で記述し、READ 命令でひとつづつ読み出していますが、
1000 READ A
～ 略 ～
2000 DATA 161,100,93,89,104,89,95,95,67,372,372,386,374,62,429,62,386,374,364
～ 略 ～
今回使用した Ichigojam の BASIC には READ～DATA の機能がありません。仕方がないのでコメント文に記述した数値を 1文字づつ PEEK関数で読み出し数値化することで代用しています。
100 P=#C02:O=2
～ 略 ～
1000 A=0:Q=P+PEEK(P)+4:IFPEEK(P+1)-39P=Q:GOTO1000
1005 C=PEEK(P+O):O=O+1:IFC>47A=10*A+C-48:GOTO1005ELSEIF!CP=Q:O=2
～ 略 ～
2000 '161,100,93,89,104,89,95,95,67,372,372,386,374,62,429,62,386,374,364
～ 略 ～
やっていることは単純なのですが、Ichigojam の BASIC のテキストがどのようにメモリに格納されているかの知識が必要となります。簡単な解説を以下の記事に書いたので興味がある方は参照してみて下さい。
『Ichigojam BASICでQuine』
他
Ichigojam の BASIC は文字列型の変数をサポートしていないためオリジナルのプログラムの
1020 C$=CHR$(B)
1030 IF C$="#" THEN PRINT:RETURN
1040 PRINT C$;
の部分は下記に書き換えています。
1030 IF B=35 THEN PRINT:RETURN
1040 PRINT CHR$(B);
1030行は
1030 IF B=ASC("#") THEN PRINT:RETURN
と書いたほうが読み易いのですが、今回はメモリ使用量をケチるために "#" のアスキーコードを直書きしております。
おわりに
おわりです。