GigglingVivian: 点の回転から複素数_Complexの計算速度を検証してみる

C言語の話題です。

はじめに

C言語では、複素数（_Complex）が仕様に加わりました。C99の拡張です。
それによって複素数が計算に使えるようになり、プログラミングの幅が広がった…
はずなんですが…

「…それ何に使うんだ？」

というのが正直な感想です。

複素数を日常生活で使うことはまずありません。
プログラムの世界でも、本来複素数を使う式は複素数を使わない方法で書き直されています。コンパイラではdouble型（実数）の高速化や、SIMDを使った並列処理が自動で行われるようになり、floatよりもdoubleの方が早いと言われる時代です。

もし複素数の計算速度が遅ければ、それを利用する必要はあまりありません。

具体的な方法

複素数を利用するもっとも身近な式として、２次元の点の回転を考えます。

点Aが与えられ、それをangleだけ回転させます。
angleの回転のx成分はMx（Mx = cos(angle);）
angleの回転のy成分はMy（My = sin(angle);）
Rx,Ryが求めたい座標だとします。

点を回転させる３つの方法を比較します。

１、複素数の乗算

ガウス平面において、乗算は点の回転を意味しますので、

p = Ax + Ay * i;
q = Mx + My * i;

result = p * q;

Rx = resultの実部;
Ry = resultの虚部;

２、加法定理

cos(a + b) = cos(a) * cos(b) - sin(a) * sin(b);
sin(a + b) = sin(a) * cos(b) + cos(a) * sin(b);

角度同士であればこの公式の通りですが、cosはx座標、sinはy座標だと考えると、次のように書き直すことができます。

Rx = Ax * Mx - Ay * My;
Ry = Ay * Mx + Ax * My;

３、アフィン変換

加法定理を行列の形に整理すると一次変換になりますが、アフィン変換では拡大縮小などのほか、平行移動を加えて３×３の行列に拡張されます。加法定理であげた式よりも冗長な計算が必要になりますが、プログラムでは良く利用されるため、随所で高速化されています。
今回はCoreGraphicsライブラリを利用しますので、

CGAffineTransform transform = CGAffineTransformRotation(CGAffineTransformIdentity, angle);
R = CGPointApplyAffineTransform(A, transform);

で計算します。

予想

仮説１：コンパイラの高速化が全ての方法で、同じようにかかる

だとすれば、加法定理が最も早く、次に複素数、最も低速なのがアフィンだと考えられます。
複素数の内部的な計算が加法定理と同じであれば、簡潔な加法定理が早く済みます。
アフィンは計算回数が多いため、最も低速です。

仮説２：GPUがアフィン変換を計算する

だとすれば、アフィン変換が最も早く、次に加法定理、最も低速なのが複素数だと考えられます。
アフィン変換はアニメーションの中核になっています。
高速なGPUが計算を行い、徹底的に最適化されている可能性は十分にあります。

仮説３：複素数が専用の方法で計算されている

だとすれば、複素数が最も早く、次に加法定理、最も低速なのがアフィンだと考えられます。

やってみた

読み出しと代入、三角関数の実行は計測時間に含めない。
「必要なデータが与えられてから、回転後の数値が確定するまで」を1回とし、
10000回の計算にかかる時間を計測する。

検証のためのコードは次の通りです。

int main(int argc, const char * argv[])
{
    @autoreleasepool {
        
        CGPoint point;
        point.x = 3.2f;
        point.y = 4.f;
        
        const int count = 10000;
        
        CGPoint xy_result;
        
        NSTimeInterval transform_end, mult_end, xy_end, simd_end;
        
        for (int ang = 0; ang < 36; ang++)
        {
            CGFloat angle = 10.f * ang - 180.f;
            angle = (angle / 180.f) * M_PI;
            
            {
                //アフィン変換
                
                CGAffineTransform transform;
                
                transform = CGAffineTransformRotate(CGAffineTransformIdentity, angle);
                
                NSDate *begin = [NSDate date];
                
                for (int i = 0; i < count; i++) {
                    
                    xy_result = CGPointApplyAffineTransform(point, transform);
                }
                
                transform_end = [[NSDate date] timeIntervalSinceDate:begin];
            }
            
            {
                //複素数のかけ算
                
                _Complex double value, mul_value, comp_result;
                comp_result = 0;
                
                value = point.x + (I * point.y);
                mul_value = cos(angle) + (I * sin(angle));
                
                NSDate *begin = [NSDate date];
                
                for (int i = 0; i < count; i++) {
                    
                    comp_result = value * mul_value;
                }
                
                xy_result.x = creal(comp_result);
                xy_result.y = cimag(comp_result);
                
                mult_end = [[NSDate date] timeIntervalSinceDate:begin];
            }
            
            {
                //加法定理
                
                CGPoint transform_xy;
                transform_xy.x = cos(angle);
                transform_xy.y = sin(angle);
                
                NSDate *begin = [NSDate date];
                
                for (int i = 0; i < count; i++) {
                    
                    xy_result.x = point.x * transform_xy.x - point.y * transform_xy.y;
                    xy_result.y = point.y * transform_xy.x + point.x * transform_xy.y;
                }
                
                xy_end = [[NSDate date] timeIntervalSinceDate:begin];
            }
            
            {
                //加法定理をSIMDにする
                //floatのため、計算精度は落ちる
                
                CGPoint transform_xy;
                transform_xy.x = cosf(angle);
                transform_xy.y = sinf(angle);
                
                __m128 simd_point, simd_transform;
                float temp[4];
                
                simd_point = _mm_set_ps((float)point.x, (float)point.y, (float)point.y, (float)point.x);
                simd_transform = _mm_set_ps((float)transform_xy.x, (float)transform_xy.y, (float)transform_xy.x, (float)transform_xy.y);
                
                NSDate *begin = [NSDate date];
                
                for (int i = 0; i < count; i++) {
                    
                    _mm_stream_ps(temp, _mm_mul_ps(simd_point, simd_transform));

                    xy_result.x = (CGFloat)(temp[3] - temp[2]);
                    xy_result.y = (CGFloat)(temp[1] + temp[0]);
                }
                
                simd_end = [[NSDate date] timeIntervalSinceDate:begin];
            }
            
            printf("%.0f,%f,%f,%f,%f\n", angle, transform_end, mult_end, xy_end, simd_end);
        }
    }
    return 0;
}

実行結果

角度	アフィン変換	複素数の乗算	加法定理	SIMDを使った加法定理
-180	0.000161	0.000034	0.000046	0.000066
-170	0.000155	0.000034	0.000049	0.000068
-160	0.000155	0.000034	0.000053	0.000067
-150	0.000158	0.000034	0.000047	0.000068
-140	0.000157	0.000034	0.000047	0.000067
-130	0.000157	0.000034	0.000046	0.000066
-120	0.000158	0.000034	0.000047	0.000068
-110	0.000159	0.000034	0.000046	0.000068
-100	0.000159	0.000034	0.000047	0.000069
-90	0.000159	0.000034	0.000047	0.000068
-80	0.000159	0.000034	0.000046	0.000067
-70	0.000160	0.000034	0.000047	0.000069
-60	0.000160	0.000033	0.000046	0.000069
-50	0.000160	0.000034	0.000046	0.000068
-40	0.000160	0.000033	0.000046	0.000069
-30	0.000157	0.000033	0.000046	0.000070
-20	0.000169	0.000034	0.000047	0.000069
-10	0.000160	0.000034	0.000046	0.000069
0	0.000159	0.000034	0.000046	0.000068
10	0.000160	0.000034	0.000046	0.000068
20	0.000159	0.000034	0.000047	0.000067
30	0.000158	0.000034	0.000047	0.000068
40	0.000159	0.000061	0.000125	0.000067
50	0.000165	0.000034	0.000046	0.000068
60	0.000163	0.000034	0.000047	0.000068
70	0.000163	0.000034	0.000046	0.000068
80	0.000164	0.000033	0.000046	0.000068
90	0.000163	0.000034	0.000046	0.000092
100	0.000166	0.000034	0.000047	0.000068
110	0.000251	0.000034	0.000042	0.000066
120	0.000166	0.000034	0.000043	0.000068
130	0.000164	0.000034	0.000047	0.000068
140	0.000163	0.000034	0.000047	0.000069
150	0.000164	0.000034	0.000047	0.000068
160	0.000163	0.000033	0.000046	0.000068
170	0.000164	0.000033	0.000046	0.000068

やるじゃないか複素数
目に見えて早いです。

SIMDが遅いのは、キャストに時間がかかっていることと、floatであること、式がSIMDに向いていないことが原因だと思います。コンパイラの最適化が強力なため、SIMDの利点が活かされていないようです。

しかし、数値の代入や三角関数の呼び出しもforループの中に含めると、次のようになります。

角度	アフィン変換	複素数の乗算	加法定理	SIMDを使った加法定理
-180	0.000886	0.000697	0.000579	0.000546
-170	0.000956	0.000689	0.000585	0.000539
-160	0.000893	0.000752	0.000562	0.000539
-150	0.001062	0.000704	0.000574	0.000613
-140	0.000952	0.000713	0.000592	0.000534
-130	0.000872	0.000693	0.000579	0.000532
-120	0.000865	0.000674	0.000574	0.000567
-110	0.000878	0.000687	0.000570	0.000558
-100	0.000858	0.000886	0.000574	0.000546
-90	0.000948	0.000713	0.000614	0.000534
-80	0.000872	0.000693	0.000579	0.000532
-70	0.000880	0.000712	0.000586	0.000532
-60	0.000901	0.000684	0.000566	0.000538
-50	0.000861	0.000674	0.000574	0.000538
-40	0.000879	0.000712	0.000587	0.000532
-30	0.000878	0.000693	0.000622	0.000535
-20	0.000977	0.000712	0.000586	0.000534
-10	0.000877	0.000685	0.000566	0.000538
0	0.000490	0.000175	0.000113	0.000263
10	0.001221	0.000739	0.000569	0.000653
20	0.000881	0.000828	0.000605	0.000535
30	0.000876	0.000712	0.000676	0.000535
40	0.000867	0.000693	0.000579	0.000534
50	0.000857	0.000678	0.000574	0.000538
60	0.000892	0.000687	0.000562	0.000538
70	0.000895	0.000694	0.000580	0.000535
80	0.000885	0.000713	0.000587	0.000532
90	0.000875	0.000693	0.000833	0.000678
100	0.000872	0.000731	0.001071	0.000539
110	0.000898	0.000687	0.000567	0.000552
120	0.000928	0.000683	0.000579	0.000544
130	0.001024	0.000715	0.000587	0.000535
140	0.000869	0.000693	0.000579	0.000532
150	0.000879	0.000712	0.000586	0.000531
160	0.000892	0.000697	0.000560	0.000539
170	0.000879	0.000675	0.000574	0.000539

複素数は明らかに遅くなります。先の表でも見た通り計算自体は遅くないのですが、creal、cimagなど、データを読み出す部分がインライン定義されていないため、ボトルネックになります。

考察

複素数には専用の計算方法が使われているらしい。

回転を繰り返す場面であれば複素数による高速化も可能性があるのではないかと思うのですが、そのほかの場面では、doubleを使う方が早いのだと考えるべきだと思います。

GigglingVivian

2014年6月3日火曜日

点の回転から複素数_Complexの計算速度を検証してみる