* ES-BASIC #7
-(by [[K]], 2019.09.21)
** (9-1) ベンチマーク#1
-ES-BASICが他の言語と比べてどのくらい速いのか、もしくはどのくらい遅いのかを比べてみようと思って、どんなプログラムで比較したらいいのかを考えてみました。
-それで世間の人はどんなプログラムで比較しているのかなと思って調べてみたら、以下のようなものがありました。
--[1]再帰を使いまくってfibo(42)を計算する。
--[2]あえて収束の悪い級数を使って（arctan(1)の級数展開）、円周率を計算する。
--[3]単純な10億回ループ。
--[4]自然数を与えたときに、その平方根までの奇数で割り算して余りを出す方法で、素数かどうかを判定する関数を書き、その関数を使ってかなり大きな数まで素数かどうかを判定する。
-これらのベンチマークは以下の問題があると思いました。
--[1]フィボナッチ数を計算するのなら、普通は再帰なんか使わずに漸化式を使って計算する。そのほうが何億倍も速い。あえて極端に遅いアルゴリズムを選んで、それで優劣を決めるのは実際の性能を反映しているだろうか？フィボナッチ数を再帰で計算しようとすれば、それは相当な回数の関数呼び出しが発生するのだが、それはもはや関数呼び出し命令のベンチマークにしかなっていない気がする。普通のプログラムはそこまで関数呼び出しばかりするわけではなく、他の演算の性能のほうがずっと重要だと思う。
--[2]普通、円周率を計算するのには、もっと収束の速い級数を使うので、このベンチマークも意味がない。プログラミング言語の作者は、アホなアルゴリズムで性能が出るかどうかは正直どうでもいいだろう。それよりも実際的なアルゴリズムでの性能を考えて言語設計をしたに違いない。力士を集めて優劣をつける時に、なぜ将棋で競わなければいけないのか。相撲でやればいいではないか。だから、あえてダメなアルゴリズムを持ってきて、それでベンチマークを取るのは有意義ではないと思う。
--[3]単純10億回ループは、まあ最速なアルゴリズムではあると思うが、しかしループの中で何もしないのなら、そもそもループする必要がないと思うので、やはり不自然だ。10億回やらなければいけない処理があって、それで10億回ループしているのなら何の問題もないのだが・・・。
--[4]単に素数探しをしたいだけなら、エラトステネスのふるいのアルゴリズムを使うのが普通だし、そのほうが何倍も速い。だからこれもベンチマークとしては適切ではないと思う。
-ということで、まともなアルゴリズムを使っているのにそれなりに処理時間がかかるような処理を探さなければいけないわけです。・・・ということで、マンデルブロ集合の計算かなと思いました。せっかく計算してもグラフィック描画ができないと表示できないのですが、表示部分があるとグラフィックができない言語ではテストできなくなってしまうので、プログラムを改造して、画像の表示はしないことにして、でも演算結果の集計値を計算してprintすることにします。そうでないと、ちゃんと計算できているのかどうかがわからないので。
-私の手元でベンチマークしたらこんな感じです：
|言語|処理時間|備考|
|Python 3.7.4|RIGHT:324.276秒|コードは(9-5)|
|Ruby 2.6.4|RIGHT:98.800秒|コードは(9-6)|
|ES-BASIC (32bit)|RIGHT:2.906秒|コードは(9-3)、JITコンパイルに要した時間も含む|
|gcc(C言語)(32bit)|RIGHT:1.954秒|コードは(9-2)、もちろん最適化レベル最大、コンパイル時間含まず|
|ES-BASIC (32bit、レジスタ変数利用)|RIGHT:1.820秒|コードは(9-7)、JITコンパイルに要した時間も含む|
|gcc(C言語)(64bit)|RIGHT:1.195秒|コードは(9-4)、もちろん最適化レベル最大、コンパイル時間含まず|
-ちなみにこのベンチマークプログラムは、以下の画像の各ピクセルの色を決定するための計算をしているだけです。
--http://k.osask.jp/files/esb20190827b.png
** (9-2) C言語版(32bit向け)
-もちろんC++もこれで行けます。
#include <stdio.h>
#include <time.h>
#define MAX 447
#define STEP (14.0 / 0x100000)
#define DN 4
int main()
{
int x, y, n, sx, sy, sn, h = 0;
double zx, zy, cx, cy, xx, yy;
for (y = 0; y < 384; y++) {
for (x = 0; x < 512; x++) {
sn = 0;
for (sx = 0; sx < DN; sx++) {
cx = (double) 0x4750 / 0x10000 + (x * DN + sx) * (STEP / DN);
for (sy = 0; sy < DN; sy++) {
cy = (double) 0x1e8 / 0x10000 + (y * DN + sy) * (STEP / DN);
zx = zy = 0;
for (n = 0; n < MAX; n++) {
xx = zx * zx;
yy = zy * zy;
if (xx + yy > 4) break;
zy = zx * zy * 2 - cy;
zx = xx - yy + cx;
}
sn += n;
}
}
// n = sn / (DN * DN);
n = sn >> 4;
h += n; // 表示の代わりに集計
}
}
printf("%d\n", h);
printf("(%.3f[sec])\n", clock() / (double) CLOCKS_PER_SEC);
return 0;
}
-うまく実行できると 28654738 が表示されると思いますが、浮動小数点の実装の微妙な違いで、もしかすると多少は違う値になるかもしれません。でもそれでもこの計算結果を画像化した場合は大差ないと思われますので問題ないと私は考えます。1%以上違う値になるのであれば、それはこのベンチマークを正しく再現できていない可能性が高いです。
** (9-3) ES-BASIC版(32bit向け)
-ES-BASICは現状のバージョンでは浮動小数点演算をサポートしていません（情けない）。整数演算のみです。・・・だから符号部1ビット、整数部7ビット、小数部24ビットの固定小数点を使って演算しています。
-ES-BASICでも定数に名前を付ける機能（ALIAS命令）を持っているので、もう少しきれいに書くこともできるのですが、このプログラムはその機能を作る前から作っていたこともあり、面倒だったので清書していません。
1 H=0
2 FOR Y=0,383
3 FOR X=0,511
4 SN=0
5 FOR SX=0,3
6 CX=(X*4+SX)*56+4673536
7 FOR SY=0,3
8 CY=(Y*4+SY)*56+124928
9 ZX=0
10 ZY=0
11 FOR N=0,446
12 XX=ZX*>>ZX,24; // これはZXとZXをIMULで計算して64bitの結果を生成した後で、SHRD命令で24ビット右シフトして下位32bitを結果とする命令
13 YY=ZY*>>ZY,24
14 T=XX+YY
15 IF T>0X4000000 GOTO SKIP
16 ZY=ZY*>>ZX,23
17 ZX=XX-YY+CX
18 ZY=ZY-CY
19 NEXT
20 LABEL SKIP
21 SN=SN+N
22 NEXT
23 NEXT
24 N=SN>>4; // N=SN/16 (N=0...447)
25 H=H+N; // 表示の代わりに集計
26 NEXT
27 NEXT
28 PRINT H
-これは実行すると 28656617 という結果になります。オリジナルの(9-2)と比べると0.007%ほど大きいことになります。これは固定小数点演算化したことによる誤差だと思われます。問題ないという認識です。
** (9-4) C言語版(64bit向け)
-もちろんC++もこれで行けます。
-ES-BASICでは「浮動小数点演算が使えないから」という理由で整数演算だけで頑張りましたが、x64ならSHRD命令が使えなくても普通のシフト命令で代用できるので、C言語でも整数命令だけでこのベンチマークをやることができます。そうすると、(9-2)より速くなります。
#include <stdio.h>
#include <time.h>
#define MAX 447
#define DN 4
int main()
{
long long x, y, sx, sy, zx, zy, cx, cy, xx, yy;
int n, sn, h = 0;
for (y = 0; y < 384; y++) {
for (x = 0; x < 512; x++) {
sn = 0;
for (sx = 0; sx < DN; sx++) {
cx = (x * 4 + sx) * 56 + 4673536;
for (sy = 0; sy < DN; sy++) {
cy = (y * 4 + sy) * 56 + 124928;
zx = zy = 0;
for (n = 0; n < MAX; n++) {
xx = (zx * zx) >> 24;
yy = (zy * zy) >> 24;
if (xx + yy > 0x4000000) break;
zy = (zx * zy) >> 23;
zx = xx - yy + cx;
zy = zy - cy;
}
sn += n;
}
}
// n = sn / (DN * DN);
n = sn >> 4;
h += n; // 表示の代わりに集計
}
}
printf("%d\n", h);
printf("(%.3f[sec])\n", clock() / (double) CLOCKS_PER_SEC);
return 0;
}
-結果としては当然のことながら(9-3)と同じく28656617 が表示されます。
-32bitのC/C++では、いろいろ工夫してもdoubleを使った場合よりも速くすることはできませんでした（インラインアセンブラでIMUL+SHRDを使えば速くできますが）。
** (9-5) Python版
-私はPython初心者なので、もっといい書き方があるかもしれません。もっといい書き方を知っている人がいたら、ぜひ教えてほしいです。
import time
start = time.time()
h = 0
CX0 = 0x4750 / 0x10000
CY0 = 0x1e8 / 0x10000
STP = 14.0 / 0x100000 / 4
for y in range(384):
for x in range(512):
sn = 0
for sx in range(4):
cx = CX0 + (x * 4 + sx) * STP
for sy in range(4):
cy = CY0 + (y * 4 + sy) * STP
zx = zy = 0
n = 0
while n < 447:
xx = zx * zx
yy = zy * zy
if xx + yy > 4: break
zy = zx * zy * 2 - cy
zx = xx - yy + cx
n = n + 1
sn += n
n = sn >> 4
h += n
print(h)
print(time.time() - start)
-これを Python 3.7.4 で実行したところ、28654739になりました。(9-2)とはわずかに違いますが、まあ誤差の範囲だと思います。
** (9-6) Ruby版
-私はRuby初心者なので、もっといい書き方があるかもしれません。もっといい書き方を知っている人がいたら、ぜひ教えてほしいです。
start_time = Time.now
h = 0
CX0 = 0x4750 / 65536.0
CY0 = 0x1e8 / 65536.0
STP = 14.0 / 0x100000 / 4
y = 0
while y < 384 do
x = 0
while x < 512 do
sn = 0
sx = 0
while sx < 4 do
cx = CX0 + (x * 4 + sx) * STP
sy = 0
while sy < 4 do
cy = CY0 + (y * 4 + sy) * STP
zx = zy = 0
n = 0
while n < 447 do
xx = zx * zx
yy = zy * zy
if xx + yy > 4 then
break
end
zy = zx * zy * 2 - cy
zx = xx - yy + cx
n = n + 1
end
sn += n
sy = sy + 1
end
sx = sx + 1
end
n = sn >> 4
h += n
x = x + 1
end
y = y + 1
end
p h
p Time.now - start_time
-これを ruby 2.6.4 で実行したところ、28654739になりました。(9-5)と同じになりました。
** (9-7) ES-BASIC版(32bit向け、レジスタ変数利用)
-ES-BASICはいくつかの変数をレジスタに割り当てることで、大幅に高速化することができます。これを使うとどのくらい速くなるのか、どのくらい可読性が落ちるのかを確認したかったので、書いてみました。
1 H=0
2 $EDI = EDI; // EDIレジスタを自由に使いたいときの前処理
3 FOR Y=0,383
4 FOR X=0,511
5 SN=0
6 FOR SX=0,3
7 CX=(X*4+SX)*56+4673536
8 FOR SY=0,3
9 CY=(Y*4+SY)*56+124928
10 ALIAS ZX:ECX, ZY:EBX, XX:ESI, YY:EDI
11 ZX=0
12 ZY=0
13 FOR N=0,446
14 XX=ZX*>>ZX,24
15 YY=ZY*>>ZY,24; // 結果はEAXにも入っている.
16 EAX=EAX+XX
17 IF EAX>0X4000000 GOTO SKIP
18 ZY=ZX*>>ZY,23
19 XX=XX-YY+CX
20 ZX=XX
21 ZY=ZY-CY
22 NEXT
23 LABEL SKIP
24 SN=SN+N
25 NEXT
26 NEXT
27 N=SN>>4; // N=SN/16 (N=0...447)
28 H=H+N; // 表示の代わりに集計
29 NEXT
30 NEXT
31 PRINT H
-これは実行すると 28656617 という結果になります。(9-3)と同じです。
-3行ほど増えていますが、これだけ速くなってくれるのであれば、私なら十分に許容範囲です。
** (9-8) ES-BASIC版(64bit向け、レジスタ変数利用)
-ES-BASICの64bit版ができたので試してみました。
1 $DI = RDI; // RDIレジスタを自由に使いたいときの前処理
2 ALIAS ZX:RCX, ZY:RDX, XX:RBX, YY:RSI, N:RDI, CX:R08, CY:R09
3 ALIAS SX:R10, SY:R11, SN:R12, H:R13, X:R14, Y:R15
4 H=0
5 FOR Y=0,383
6 FOR X=0,511
7 SN=0
8 FOR SX=0,3
9 CX=(X*4+SX)*56+4673536
10 FOR SY=0,3
11 CY=(Y*4+SY)*56+124928
12 ZX=0
13 ZY=0
14 FOR N=0,446
15 XX=(ZX*ZX)>>24
16 YY=(ZY*ZY)>>24
17 RAX=XX+YY
18 IF RAX>0X4000000 GOTO SKIP
19 ZY=(ZX*ZY)>>23
19 ZY=(ZY*ZX)>>23
20 XX=XX-YY
21 ZX=XX+CX
22 ZY=ZY-CY
23 NEXT
24 LABEL SKIP
25 SN=SN+N
26 NEXT
27 NEXT
28 N=SN>>4; // N=SN/16 (N=0...447)
29 H=H+N; // 表示の代わりに集計
30 NEXT
31 NEXT
32 PRINT H
-これは実行すると 28656617 という結果になります。(9-3)と同じです。
-3行ほど増えていますが、これだけ速くなってくれるのであれば、私なら十分に許容範囲です。
* こめんと欄
#comment
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
