Mandelbrot Set (1)

マンデルブロ集合は z_k+1 = z _k² + c（ k → ∞）で k を無限大にしたときの極限 {z_k} が無限大に発散しない複素数 c 全体の集合です。普通によく見るマンデルブロ集合は z0 = 0 としています。またマンデルブロ集合の彩色の基本は有限回数内の計算で無限大に発散しない点の集合を黒で描くのが基本とされているようです。

形式的な堂々巡りがこんな形を描くなんて。描かずにはいられない。

どうするの？
簡単なテストを。複素数の計算には complex.h の _Complex型 と Accelerate.framework の 構造体 DSPSplitComplex, DSPDoubleSplitComplex をそれぞれ使ってみました。おおざっぱな時間がわかればよいので計算も適当です。Xcode で Command Line Tool, Type を Foundation で新規プロジェクトを作成、Accelerate.framework を追加し、実装ファイルに <Accelerate/Accelerate.h> をインポートします。

#import <Foundation/Foundation.h>
#import <Accelerate/Accelerate.h>


int main (int argc, const char * argv[]) {
    NSAutoreleasePool * pool = [[NSAutoreleasePool alloc] init];
 NSLog(@"---start---");
 
 ///----------- 時間 -----------///
 clock_t start, end;
 
 ///----------- z = z * z + c の計算をする。 -----------///
 double width, height, countLimit;
 width = 600.0;
 height = 600.0;
 countLimit = 10000.0;
 
 /////////// Normal ///////////
 
 double _Complex z, c;
 
 
 c = (-1.5) + 0.0 * I;
 
 sleep(1);
 start = clock();
 for (NSInteger x = 0; x < width; ++x) {
  for (NSInteger y = 0;y < height ; ++y) {
   
   z = c;
   
   for (NSInteger count = 0; count < countLimit; ++count) {
  
    z = z * z + c;
  
   }
   
  }
 }
 end = clock();
 NSLog(@"complex:%fs, %e",(double)(end - start) / CLOCKS_PER_SEC, creal(z));
 
 ////////// vDSP 基本設定 ///////////
 vDSP_Stride stride = 1;
 vDSP_Length length = 1 << (int)ceil(log2(height));
 
 /////////// vDSP DSPSplitComplex ///////////
 DSPSplitComplex zF, cF;
 
 zF.realp = (float *)calloc(length,sizeof(float));
 zF.imagp = (float *)calloc(length,sizeof(float));
 cF.realp = (float *)calloc(length,sizeof(float));
 cF.imagp = (float *)calloc(length,sizeof(float));
 
 for (NSInteger i = 0; i < height; ++i) {
  cF.realp[i] = (-1.5); 
 }
 
 sleep(1);
 start = clock();
 for (NSInteger x = 0; x < width; ++x) {
  vDSP_zvmov(&cF, stride, &zF, stride, length);
  for (NSInteger count = 0; count < countLimit; ++count) {
    
   vDSP_zvmul(&zF, stride, &zF, stride, &zF, stride, length, 1);
   
   vDSP_zvadd(&zF, stride, &cF, stride, &zF, stride, length);
    
  }
 }
 end = clock();
 NSLog(@"DSPSplitComplex:%fs, %e",(double)(end - start) / CLOCKS_PER_SEC, zF.realp[1]);
 
 free(zF.realp);
 free(zF.imagp);
 free(cF.realp);
 free(cF.imagp);
 
 
 ////////// vDSP DSPDoubleSplitComplex ///////////
 DSPDoubleSplitComplex zD, cD;
 
 zD.realp = (double *)calloc(length,sizeof(double));
 zD.imagp = (double *)calloc(length,sizeof(double));
 cD.realp = (double *)calloc(length,sizeof(double));
 cD.imagp = (double *)calloc(length,sizeof(double));
 

 for (NSInteger i = 0;i < height; ++i) {
  cD.realp[i] = (-1.5);
 }
 

 
 sleep(1);
 start = clock();
 for (NSInteger x = 0; x < width; ++x) {
  vDSP_zvmovD(&cD, stride, &zD, stride, length);
  for (NSInteger count = 0; count < countLimit; ++count) {
   
   vDSP_zvmulD(&zD, stride, &zD, stride, &zD, stride, length, 1);
   
   vDSP_zvaddD(&zD, stride, &cD, stride, &zD, stride, length);
   
  }
 }
 end = clock();
 NSLog(@"DSPDoubleSplitComplex:%fs, %e",(double)(end - start) / CLOCKS_PER_SEC, zD.realp[1]);
 
 free(zD.realp);
 free(zD.imagp);
 free(cD.realp);
 free(cD.imagp);
    
 NSLog(@"---end---");
 [pool drain];
    return 0;
}

結果は？
vDSPvsNormal[789:903] ---start---
vDSPvsNormal[789:903] complex:13.345002s, -1.473274e+00
vDSPvsNormal[789:903] DSPSplitComplex:8.198444s, -1.282221e+00
vDSPvsNormal[789:903] DSPDoubleSplitComplex:16.541649s, -1.473274e+00
vDSPvsNormal[789:903] ---end---

ガウス平面上の（-1.5, 0）は振動し続け、計算回数が多くなる点なので適当に選んでみました。float と double で誤差があります。ただ単に集合を描くのであれば、座標を各点に置き換えて計算結果を書き出してやればよいですが、もやもやした感じを出したい場合には個別の点について発散条件に至ったときの計算回数を保存してやる必要があります。試していませんが塗り重ねてもできるかも知れません。