２変数関数の等高線図描画プログラム

ダウンロード

関数 Func(x,y) の返り値を高さとして、その等高線をウィンドウに表示する。

ウィンドウの隅をマウスの左ボタンを押しながら動かして
ウィンドウの大きさを変更したり、最大化ボタンを押すなどすると、
すぐに計算し直して表示する。

高さが 0.0 の等高線は黒色、
1.0 以上の等高線は水色、
負の等高線は青色で表示する。
等高線の間隔は高さが 1.0 毎である。
等高線の間隔の高さを例えば 0.2 毎にしたい場合は関数の戻り値を 5.0 倍にする。

/* 以下の三つの定数で平面上のどの領域が表示されるかが決定される */
#define CENTER_X 0.0 /* ウィンドウ中心点の平面上での位置 */
#define CENTER_Y 0.0
#define SCALE    2.0 /* ウィンドウにこれを半径とした円が入る */

Func(CENTER_X-SCALE,CENTER_Y+SCALE) が Func(CENTER_X,CENTER_Y) よりも
大きければ計算時に ceil() を使用するモードになり、
そうでなければ計算時に floor() を使用するモードになる。
CHFC を定義すると ceil() を使用するか floor() 使用するかを逆転する。
この時、表示される等高線の位置が少し変化する。

simple.cpp、cont.cpp、cont2.cpp では処理時間の表示を１ミリ秒単位に
するため timeBeginPeriod(1) と timeEndPeriod(1) を呼び出している。
このため、TDM-GCC-64 ではコンパイル時に -lwinmm をコマンドラインの
最後に加えておく。timeBeginPeriod(1) と timeEndPerido(1) を呼び出さない
ようにするにはコンパイル時に NO_PERIOD を定義する。

シンプルで遅くてメモリを節約するバージョン

【simple.h simple.c simple.rc】

バッファを経由せず SetPixel() で直接画面に描画し表示する。

simple.c に
double Func(const double x,const double y) { return(2-y)*(2+y)*y*y-x*x; }
のように関数を書く。

高速バージョン

【cont.h cont.cpp cont.rc】

cont.cpp に
real Func(const real&x,const real&y) { return(2-y)*(2+y)*y*y-x*x; }
のように関数を書く。real は通常 double(６４ビット) である。

描画高速化

描画はバッファに対して行い、画面表示は BitBlt() で一気に行う。
さらに３２ビットカラーの場合は SetPixel() を使わずにビットマップイメージに
直接３２ビットカラー値を書き込んで描画して高速化している。
CreateBitmap はデバイス依存なので、色が正しく表示されない可能性がある。
この場合 PXL を定義して３２ビットカラーでも強制的に SetPixel() で描画すると
正しい色で表示される。この時メモリが節約されるが遅くなる。

マルチスレッド

cont.cpp はマルチスレッドに対応していて計算を高速化している。
FETCH_INC() マクロが排他制御を行っている。
CRIT を定義すると FETCH_INC() をクリティカルセクションで実現する。
これは 80386 のマルチＣＰＵで実行する場合に定義する。
VisualC++2005 以降や TDM-GCC 等では CRIT を定義しない場合 FETCH_INC() を
組み込み関数 _InterlockedExchangeAdd で実現する。
VisualC++2003 以前や BorlandC++ で CRIT を定義しない場合 FETCH_INC() は
以下の i486 以降のマシン語ルーチン
    "\xB8\x01\x00\x00\x00" /* mov eax,1   */
    "\x0F\xC1\x01"         /* lock xadd dword ptr[ecx],eax */
    "\xC3";                /* ret */
を呼び出す。
MCH を定義すると VisualC++2005 以降や TDM-GCC 等でも
FETCH_INC() はマシン語ルーチンを呼び出す。
g++ のバージョンが古い場合は MCH を定義する。

STH を定義すると全てシングルスレッドで計算する。
メモリが節約されるが遅くなる。

ベクトル化

AVX 等を定義することによりベクトル化による高速化も選択できる。

AVX/AVXF を定義すると AVX で計算する。   (real は２５６ビットになる)
SSE4/SSE4F を定義すると SSE4 で計算する。(real は１２８ビットになる)
SSE2/SSE2F を定義すると SSE2 で計算する。(real は１２８ビットになる)

AVXF を定義すると計算精度は単精度で８個の数値を同時に計算する。
SSE4F/SSE2F を定義すると計算精度は単精度で４個の数値を同時に計算する。
AVX を定義すると計算精度は倍精度で４個の数値を同時に計算する。
SSE4/SSE2 を定義すると計算精度は倍精度で２個の数値を同時に計算する。


６４ビット TDM-GCC で path C:\TDM-GCC-64\bin として
g++ -Wall -m64 -DAVX cont.cpp -lgdi32 -lwinmm
または
g++ -Wall -m64 -DAVXF cont.cpp -lgdi32 -lwinmm
とするとコンパイルはできるが実行はできない。

しかし
g++ -Wall -O4 -m64 -DAVX cont.cpp -lgdi32 -lwinmm
または
g++ -Wall -O4 -m64 -DAVXF cont.cpp -lgdi32 -lwinmm
のように -O4 または -O3 を指定すると実行できる。
-m64 の代わりに -m32 のように３２ビットでコンパイルしたり
SSE2 SSE2F SSE4 SSE4F の場合は -O4 や -O3 を加えなくても実行できる。


VisualC++2005 以降や TDM-GCC 等ではＣＰＵが対応しているかどうかを
組み込み関数 __cpuid() を使って確認する。
VisualC++2003 以前ではマシン語ルーチンを呼び出して確認する。
MCH を定義すると VisualC++2005 以降や TDM-GCC 等でも
マシン語ルーチンを呼び出す。
g++ のバージョンが古い場合は MCH を定義する。


(1)Func(x,y) 内でベクトル演算が行われる演算／関数
+ - * / += -= *= /= ++ -- 等の四則演算
ceil fabs floor fmax fmin ldexp modf

(2)Func(x,y) 内でベクトル演算が行われないが対応している関数
acos asin atan cos cosh exp log log10 sin sinh sqrt tan tanh
aton2 hypot pow

Func(x,y) 内に(1)が多ければベクトル化の効果があるが、
(1)に対して(2)が多ければベクトル化しない方が高速になる場合がある。


ベクトル化した場合 ldexp(x,i) はオーバーフローやアンダーフローの確認は
しない。倍精度の場合 i は±1020、単精度の場合±120 程度までしか正しく
計算できない。

また modf(x,p) は x が±２０億程度までしか正しく計算できない。


SSE2 では SSE4 と違って
 _mm_floor_ps, _mm_floor_pd, _mm_ceil_ps, _mm_ceil_pd 
が使えない。

SSE2 が定義されているときは、
floor の代わりに高さに 65536 を足してから整数化し 65536 を引いて戻す。
ceil の場合は 65536 を引いてから整数化し65536 を足して戻す。

従って高さがちょうど 65536 か -65536 の等高線は表示されない場合がある。
また、-2.4e-324 〜 -3.6e-12 の数値は
floor すれば -1 になるはずなのに 0 となってしまう。
同様に、1.0-5.5e-17 〜 1.0-7.2e-12 の数値は
floor すれば 0 になるはずなのに 1 となってしまう。

SSE2F が定義されていれば前述の 65536 が 256 になる。
従って高さがちょうど 256 か -256 の等高線は表示されない場合がある。
また、-7.0e-46 〜 -7.6e-6 の数値は
floor すれば -1 になるはずなのに 0 となってしまう。
同様に、1.0-2.9e-8 〜 1.0-1.52e-6 の数値は
floor すれば 0 になるはずなのに 1 となってしまう。

SSE4/SSE4F/AVX/AVXF を定義すれば
この不具合は起きない。また SSE2/SSE2F を定義せずにベクトル化を諦めても
この不具合は起きない。

開発用のマクロ定義

TEST を定義すると高さ計算のみをシングルスレッドで
１００回繰り返して性能を計る。

BUG_TEST を定義するとマルチスレッドで排他制御しない悪影響の確認ができる。

更に高速なバージョン

【cont.h cont2.cpp cont2.rc】

Func(x,y) の代わりに Func(Fx(x),Fy(y)) を高さとする。
cont2.cpp では実行時に Func() が呼び出される回数は cont.cpp と同じだが、
Fx() と Fy() の呼び出し回数は必要最低限である。
従って Func() の演算が cont.cpp に比べて少なくなった場合は効果がある。

(例1) cont.cpp で
real Func(const real&x,const real&y) { return(2-y)*(2+y)*y*y-x*x; }
としていたなら、cont2.cpp では、
real Fx(const real&x) { return x*x; }
double Fy(const double y) { return(2-y)*(2+y)*y*y; }
real Func(const real&fx,const real&fy) { return fy-fx; }
のように関数を書く。
この場合 Func() 内の演算数が少なくなるので cont2.cpp の方が高速になる。

(例2) cont.cpp で
real Func(const real&x,const real&y) { return y*(x-1)+(1-y)*x+(y+1)*(x+1); }
のようにしていた場合、cont2.cpp では
real Fx(const real&x) { return x; }
double Fy(const double y) { return y; }
real Func(const real&fx,const real&fy)
{ return fy*(fx-1)+(1-fy)*fx+(fy+1)*(fx+1); }
と関数を書くしかないので cont2.cpp にしても全く効果はない。

関数の書き方以外の仕様は cont.cpp と同じである。