* ES-BASIC #12
-(by [[K]], 2019.11.21)
** はじめに
-このページは[[esbasic0007]]の内容が古くなってしまったので、全面的に書き直したものです。
** (19-1) ベンチマーク ver.1.1
-ES-BASICが他の言語と比べてどのくらい速いのか、もしくはどのくらい遅いのかを比べてみようと思って、どんなプログラムで比較したらいいのかを考えてみました。
-それで世間の人はどんなプログラムで比較しているのかなと思って調べてみたら、以下のようなものがありました。
--[1]再帰を使いまくってfibo(42)を計算する。
--[2]あえて収束の悪い級数を使って（arctan(1)の級数展開）、円周率を計算する。
--[3]単純な10億回ループ。
--[4]自然数を与えたときに、その平方根までの奇数で割り算して余りを出す方法で、素数かどうかを判定する関数を書き、その関数を使ってかなり大きな数まで素数かどうかを判定する。
-これらのベンチマークは以下の''「問題」''があると思いました。
--[1]フィボナッチ数を計算するのなら、普通は再帰なんか使わずに漸化式を使って計算する。そのほうが何億倍も速い。あえて極端に遅いアルゴリズムを選んで、それで優劣を決めるのは実際の性能を反映しているだろうか？フィボナッチ数を再帰で計算しようとすれば、それは相当な回数の関数呼び出しが発生するのだが、それはもはや関数呼び出し命令のベンチマークにしかなっていない気がする。普通のプログラムはそこまで関数呼び出しばかりするわけではなく、他の演算の性能のほうがずっと重要だと思う。
--[2]普通、円周率を計算するのには、もっと収束の速い級数を使うので、このベンチマークも意味がない。プログラミング言語の作者は、アホなアルゴリズムで性能が出るかどうかは正直どうでもいいだろう。それよりも実際的なアルゴリズムでの性能を考えて言語設計をしたに違いない。力士を集めて優劣をつける時に、なぜ将棋で競わなければいけないのか。相撲でやればいいではないか。だから、あえてダメなアルゴリズムを持ってきて、それでベンチマークを取るのは有意義ではないと思う。
--[3]単純10億回ループは、まあ最速なアルゴリズムではあると思うが、しかしループの中で何もしないのなら、そもそもループする必要がないと思うので、やはり不自然だ。10億回やらなければいけない処理があって、それで10億回ループしているのなら何の問題もないのだが・・・。
--[4]単に素数探しをしたいだけなら、エラトステネスのふるいのアルゴリズムを使うのが普通だし、そのほうが何倍も速い。だからこれもベンチマークとしては適切ではないと思う。
-ということで、''「まともなアルゴリズムを使っているのにそれなりに処理時間がかかるような処理」を探さなければいけない''わけです。
----
-・・・それで思いついたのが「マンデルブロ集合」です。これは具体的には
--http://k.osask.jp/files/esb20190827b.png
-この絵を計算して描画する処理です。
-プログラムは以下のようになります（blaライブラリを使っています）。
// "19-1a.c"
#include <stdio.h>
#include <time.h>
#define STATIC static
#include "bla.c"
#define MAX 447
#define STEP (14.0 / 0x100000)
#define DN 4
int bla_main(int argc, const char **argv)
{
int x, y, sx, sy, n, sn, c;
float cx, cy, zx, zy, xx, yy, tx, ty;
bla_Window *win = bla_openWin(512, 384, "mandel", 1);
for (y = 0; y < 384; y++) {
for (x = 0; x < 512; x++) {
sn = 0;
for (sx = 0; sx < DN; sx++) {
cx = (float) 0x4750 / 0x10000 + (x * DN + sx) * (STEP / DN);
for (sy = 0; sy < DN; sy++) {
cy = (float) -0x1e8 / 0x10000 - (y * DN + sy) * (STEP / DN);
zx = zy = 0.0;
for (n = 0; n < MAX; n++) {
xx = zx * zx;
yy = zy * zy;
if (xx + yy > 4.0) break;
tx = xx - yy;
ty = zx * zy * 2.0;
zx = tx + cx;
zy = ty + cy;
}
sn += n;
}
}
n = sn / (DN * DN);
c = bla_rgb(0, 0, 0);
if (n < 256)
c = bla_rgb(n, 0, 0);
else if (n < MAX)
c = bla_rgb(255, n - 255, 0);
bla_setPix(win, x, y, c);
}
bla_leapFlushAll(win, 300);
}
printf("time=%.3f[sec]\n", (float) clock() / CLOCKS_PER_SEC);
bla_end();
return 0; // dummy
}
-このプログラムを手元の環境で実行したら、1.815秒でした。
-しかしこのプログラムには問題があります。
--[1]グラフィック描画をサポートしていない言語には移植できない。
--[2]グラフィック描画処理が重い環境では、言語の純粋な計算力の比較ではなくなるのではないか？
-ということで、上記のプログラムからグラフィック処理を全部抜き取って、各ピクセルごとのnの値の合計を計算することにしました。
--本当は合計を計算する必然性はなかったのですが、しかし各ピクセルごとに計算したあたりを使わないと、C言語などでは最適化によってnの計算そのものを省略することがあり得て、それだとさすがに困ってしまうので、合計を計算することにしたのです。
--本当は合計を計算する必然性はなかったのですが、しかし各ピクセルごとに計算した値を使わないと、C言語などでは最適化によってnの計算そのものを省略することがあり得て、それだとさすがに困ってしまうので、合計を計算することにしたのです。
----
-ということでできあがったのが以下のプログラムです。
// "19-1b.c"
#include <stdio.h>
#include <time.h>
#define MAX 447
#define STEP (14.0 / 0x100000)
#define DN 4
int main()
{
int x, y, sx, sy, n, sn, h = 0;
float cx, cy, zx, zy, xx, yy, tx, ty;
for (y = 0; y < 384; y++) {
for (x = 0; x < 512; x++) {
sn = 0;
for (sx = 0; sx < DN; sx++) {
cx = (float) 0x4750 / 0x10000 + (x * DN + sx) * (STEP / DN);
for (sy = 0; sy < DN; sy++) {
cy = (float) -0x1e8 / 0x10000 - (y * DN + sy) * (STEP / DN);
zx = zy = 0.0;
for (n = 0; n < MAX; n++) {
xx = zx * zx;
yy = zy * zy;
if (xx + yy > 4.0) break;
tx = xx - yy;
ty = zx * zy * 2;
zx = tx + cx;
zy = ty + cy;
}
sn += n;
}
}
n = sn / (DN * DN);
h += n;
}
}
printf("h=%d\n", h);
printf("time=%.3f[sec]\n", (float) clock() / CLOCKS_PER_SEC);
return 0; // dummy
}
-これを実行すると、手元の環境では以下のようになりました。
h=28656157
time=1.807[sec]
----
-私は当初、この「19-1b.c」で満足していました。しかし、このプログラムを実行するには、float型（つまり浮動小数点演算）を言語がサポートしている必要があって、まあたいていの言語はサポートしているわけですが、しかし自作言語では最初は整数型だけになるのが普通で、そうなると言語の完成度がかなり上がらないとこのベンチマークは試せないということになります。
-それは私には不満でした。できるだけ早期の段階で性能テストできる方が、基本部分の設計を早い段階で練ることができるわけで、後になったらもう直したり試行錯誤したりするのが大変になるかもしれないからです。
-ということで、floatの利用をやめて、cx～tyは小数部が24ビットの固定小数点にしてみました。
-これで足し算や引き算はうまく行くのですが、掛け算は少し問題があります。というのは乗算してから24ビットシフトする必要があるのですが、乗算結果は一時的に32ビットを超えてしまうので、64bit整数にキャストしなければいけません。しかしC言語では、掛け算ごとにキャストして頑張るよりも、そもそもcx～tyを最初から64bit整数にしておく方がずっと高速でした。ということで、こんなプログラムになりました。
// "19-1c.c"
#include <stdio.h>
#include <time.h>
int main()
{
int x, y, sx, sy, n, sn, h = 0;
long long cx, cy, zx, zy, xx, yy, tx, ty;
for (y = 0; y < 384; y++) {
for (x = 0; x < 512; x++) {
sn = 0;
for (sx = 0; sx < 4; sx++) {
cx = (x * 4 + sx) * 56 + 4673536;
for (sy = 0; sy < 4; sy++) {
cy = (y * 4 + sy) * -56 - 124928;
zx = zy = 0;
for (n = 0; n < 447; n++) {
xx = (zx * zx) >> 24;
yy = (zy * zy) >> 24;
if (xx + yy > 0x4000000) break;
tx = xx - yy;
ty = (zx * zy) >> 23;
zx = tx + cx;
zy = ty + cy;
}
sn += n;
}
}
n = sn / 16;
h += n;
}
}
printf("h=%d\n", (int) h);
printf("time=%.3f[sec]\n", (float) clock() / CLOCKS_PER_SEC);
return 0; // dummy
}
-これを実行すると、手元の環境では以下のようになりました。
h=28656617
time=1.177[sec]
-「19-1b.c」の結果と比べると、hが少しだけ増えています。これはfloatと24ビット固定小数点での精度の差に由来するもので、本質的なものではありません。実際この計算結果で画面描画してみても、違いを感じることはできませんでした。
-そして1.807秒と比較すると1.5倍以上も速くなっていたのです！・・・つまり、C言語のように浮動小数点演算を十分にサポートできている場合であっても、整数演算だけで代用できるのであれば、そのほうがずっと高速だということです。なるほど勉強になります。
-そういう代用を検討する場合は、64ビット整数が使えると非常に重宝するので、64ビット化というのはこういう風に活用すればいいんだなと悟りました（私はそれまで、整数なんて32ビットもあれば十分だ、64ビット整数なんてあっても使いきれない、メモリの無駄だ！とか思っていたのです）。
----
-この「19-1c.c」での大幅な高速化が気に入って、私はもっと高速化できる方法はないかなと考えてみることにしました。すると、まず、txやtyがいらないと思いました。
zy = ((zy * zx) >> 23) + cy;
zx = xx - yy + cx;
-とでもすればそれで済む話だだからです。
-さらに、xやyを工夫すれば、sxやsyもなくせると思いました。
-さらに加えて、zx=zy=0で回り始めるのなら、n=0のときのif文でbreakすることなんてありえないわけで、すると zx=cx; zy=cy; になってn=1になるので、この初期条件で回してもいいはずです。その方が速いはずです。
-さらに、cxが変化しなければ、最初のzx*zxの値は変化しないはずです。だからcyのループに入る前にxx0として計算して置いてこれを流用するようにします。
-さらに、 n=1～447 でループするのではなく、n= -446～0 でループするようにします。これならループ中は0との比較をするだけですむのでわずかに高速化できるかもしれません。このつじつまを合わせるために、snの初期値は0ではなく、447*16=7152にします。
-ええとなぜこのような努力をするのかというと、「まともなアルゴリズムを使う」ことがベンチマークでは何よりも大事だと私は思っているからです。良いアルゴリズムが使えるのに、それを使わずに計算速度の優劣を競っても、有意ではないと思っているのです。
-これらすべての変更を盛り込んだ究極形が以下の「19-1d.c」です。
// "19-1d.c"
#include <stdio.h>
#include <time.h>
int main()
{
int n, sn, h = 0;
long long x, y, cx, cy, zx, zy, xx, yy, xx0;
for (y = -125152; y > -211168; y -= 224) {
for (x = 4673760; x < 4788448; x += 224) {
sn = 7152;
for (cx = x - 224; cx < x; cx += 56) {
xx0 = (cx * cx) >> 24;
for (cy = y + 224; cy > y; cy -= 56) {
zy = cy;
zx = cx;
xx = xx0;
yy = (cy * cy) >> 24;
n = -446;
do {
zy = (zy * zx) >> 23;
zx = xx + cx - yy;
zy = zy + cy;
n++;
if (n >= 0) goto skip;
xx = (zx * zx) >> 24;
yy = (zy * zy) >> 24;
} while (xx + yy <= 0x4000000);
sn += n;
skip:
;
}
}
h += sn >> 4;
}
}
printf("h=%d\n", h);
printf("time=%.3f[sec]\n", (float) clock() / CLOCKS_PER_SEC);
return 0; // dummy
}
-これを実行すると、手元の環境では以下のようになりました。
h=28656617
time=1.115[sec]
-hの値が全く変わっていないので、正しく計算できていることがよくわかります。「19-1c」と比較すると1.055倍くらい高速化できたことになります。
----
-参考までに、floatを使いつつ精一杯高速化するとこうなります。これはC言語においては全然速くはありませんが、一時的に64ビットになるような乗算が不得意なら、こちらのアルゴリズムの方が速いかもしれません。
// "19-1e.c"
#include <stdio.h>
#include <time.h>
int main()
{
int x, y, ix, iy, n, sn, h = 0;
float cx, cy, zx, zy, xx, yy, xx0;
for (y = -125152; y > -211168; y -= 224) {
for (x = 4673760; x < 4788448; x += 224) {
sn = 7152;
for (ix = x - 224; ix < x; ix += 56) {
cx = ix * (1.0 / 16777216.0);
xx0 = cx * cx;
for (iy = y + 224; iy > y; iy -= 56) {
cy = iy * (1.0 / 16777216.0);
zy = cy;
zx = cx;
xx = xx0;
yy = cy * cy;
n = -446;
do {
zy = zy * zx * 2.0;
zx = xx + cx - yy;
zy = zy + cy;
n++;
if (n >= 0) goto skip;
xx = zx * zx;
yy = zy * zy;
} while (xx + yy <= 4.0);
sn += n;
skip:
;
}
}
h += sn >> 4;
}
}
printf("h=%d\n", h);
printf("time=%.3f[sec]\n", (float) clock() / CLOCKS_PER_SEC);
return 0; // dummy
}
-これを実行すると、手元の環境では以下のようになりました。
h=28655896
time=1.761[sec]
-hの値がまた少し変わっていますが、これもおそらく精度に由来する誤差であり、気にしなくてよいと思われます。
** (19-2) gccの最適化について
-今回、19-1b → 19-1c → 19-1d と改良してきたわけですが、この改良のつど、処理時間は確実に短くなっており、つまり人間のセンスによる最適化は、コンパイラによる自動的な最適化を十分に上回っていたことになります。もし上回ることができていなければ、コンパイラは最適化の結果として同じバイナリにたどり着いて同じ処理時間になっていたはずなのです。そうならなかったのは、コンパイラの最適化は、まだまだ人間には及ばないという証拠なのです。
-だから私はES-BASICにおいて言語が最適化を頑張ることを放棄して、人間による最適化を支援する方向で言語を設計しているわけです。
** つづく
-[[esbasic0013]]へつづく
* こめんと欄
#comment
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
