マンデルブロ集合の面積

はじめに

　この記事は、マンデルブロ集合の面積を求めることをテーマに、計算速度を測ったり、piとeとiについて、調べることを目的とした実験の一部です。マンデルブロ集合の周囲の長さは無限大ですが、面積は約1.5の無理数と決まっているので、題材として扱いやすいのです。
　アルゴリズムの選定をしたのち、マルチスレッド化等の高速化をしました。

面積の解

　サイトによると√(6π-1)-e

\sqrt{6\pi-1}-e

で求まるようだが、解析解なのか、近似式なのか？
エクセル式は

=SQRT(6*PI()-1)-EXP(1)

1.506591651485500000000000000000 

wolframによる計算結果は1.5065916514855032852705345234769657413403164236938572609274770450

mrob.comによるとGPUでカウントして
1.5065918561
小数点以下６桁が一致

自作プログラムでは、
1.5066597
小数点以下3桁が一致

processingサンプルのプログラムは遅い？

　processingにはマンデルブロ集合を描くサンプルが入っている。しかし、平方根の計算が入っていて、なんだか遅そう。そこでマンデルブロ集合の面積（ピクセル数カウント）を計算しながら検証。

条件

　M1Macで実行。
　サイズ 512x512 で 65536回の設定で計算。
　結果が63277ピクセルであることを確認。

結果

①サンプルをいじって、必要な部分だけを残して計算。
25.7秒
②平方根のないコード
14.7秒

結論

　やっぱり標準サンプルのコードは遅く、1.75倍の差が出ました。平方根関数が遅いだけでなく、floatしか扱えないので精度問題も残ります。②のコードでは、floatで計算してもdoubleで計算しても、速度に違いが出ませんでした。なぜか？（追記：コード中の数字にdを明記した）

検証コード

mandelArea1.pde

size(512, 512);
noLoop();
background(255);

float w = 2.5;
float h = 2.5;

float xmin = -2.0;
float ymin = -1.25;

int maxiterations = 65536;
int count = 0;

float xmax = xmin + w;
float ymax = ymin + h;

float dx = (xmax - xmin) / (width);
float dy = (ymax - ymin) / (height);

float y = ymin;
for (int j = 0; j < height; j++) {
  float x = xmin;
  for (int i = 0; i < width; i++) {
    float a = x;
    float b = y;
    int n = 0;
    float max = 4.0;
    while (n < maxiterations) {
      float aa = a * a;
      float bb = b * b;
      float abs = sqrt(aa + bb);
      if (abs > max) {
        break;
      }
      float twoab = 2.0 * a * b;
      a = aa - bb + x;
      b = twoab + y;
      n++;
    }
    if (n == maxiterations) {
      count++;
    }
    x += dx;
  }
  y += dy;
}
println(count);
println(millis());

mandelArea2.pde

double w = 2.5;
double h = 2.5;

double xmin = -2.0;
double ymin = -1.25;

double xmax = xmin + w;
double ymax = ymin + h;

double dx = (xmax - xmin) / 512;
double dy = (ymax - ymin) / 512;

int maxiteration = 65536;
int count = 0;

void setup() {
  size(512, 512);
  noLoop();

  double y = ymin;
  for (int j = 0; j < 512; j++) {
    double x = xmin;
    for (int i = 0; i < 512; i++) {
      count += mandelbrotJulia(0, 0, x, y, maxiteration);
      x += dx;
    }
    y += dy;
  }
  println(count);
  println(millis());
}

void draw() {
}

int mandelbrotJulia(double zr, double zi, 
  double cr, double ci, 
  int iteration) {

  for (int n=0; n<iteration; n++) {
    double tmp;
    tmp = zr*zr - zi*zi + cr;
    zi = 2*zr*zi + ci;
    zr = tmp;
    if (zr*zr + zi*zi>4) {
      return 0;
    }
  }
  return 1;
}

スレッドで高速化

　速い方のコードで、マルチスレッド化した。processingではthreadの機能があるが、値を渡せないので使わず。過去にjavaのRunnableを使ったが、値が若干おかしくなって、直せなかったので、今回は使わず。
　Runnableではなく、値を渡して返すことも可能なCallableを使った。計算したい範囲などの値を渡して、発散・非発散かを返す。とりあえず、64スレッドに分けて、同時に実行。無駄が生じるが、無視。
　また収束することが明らかな、カージオイドといくつかの円を計算から除外した。

結果

　M1Macで65536x65536pixelで、反復65536で40分ぐらい。面積は1.5066597

コード

mandelArea2021_2.pde

import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.concurrent.Callable;

void setup() {
  size(512, 512);
  noLoop();

  println("div,iter,pixels,area,millis");

  double[] w = new double[64];
  double[] h = new double[64];
  double[] xmin = new double[64];
  double[] ymin = new double[64];

  for (int i=0; i<64; i++) {
    w[i] = 2.5d/8.0d;
    h[i] = 2.5d/8.0d;
    xmin[i] = -2.0d+2.5d/8.0d*(i%8);
    ymin[i] = -1.25d+2.5d/8.0d*(i/8);
  }


  for (int div=16384; div<=65536/4; div*=2) {
    for (int iter=16384; iter<=65536/4; iter*=2) {
      //スレッドの実行
      ExecutorService[] ex = new ExecutorService[64];
      Future<Integer>[] futureResult = new Future[64];
      for (int i=0; i<64; i++) {
        ex[i] = Executors.newSingleThreadExecutor();
        futureResult[i] = ex[i].submit(new calcArea(w[i], h[i], xmin[i], ymin[i], div/8, iter));
      }

      //結果の取得
      try {
        int result = 0;
        for (int i=0; i<64; i++) {
          result += futureResult[i].get();
        }
        //println(result);
        print(div+",");
        print(iter+",");
        print(result+",");
        print(result*2.5d*2.5d/div/div+",");
        println(millis()+"msec");
      } 
      catch (Exception e) {
        System.out.println(e);
      }
    }
  }
}

void draw() {
}

public class calcArea implements Callable<Integer> {
  double xmax, ymax, dx, dy;
  double w, h, xmin, ymin, size;
  int maxIteration;
  int areaPixels = 0;

  public Integer call() {
    double y = ymin;
    for (int j = 0; j < size; j++) {
      double x = xmin;
      for (int i = 0; i < size; i++) {
        if (outOfCardioidCircles(x, y) ) {
          areaPixels += mandelbrotJulia(0d, 0d, x, y, maxIteration);
        } else {
          areaPixels++;
        }
        x += dx;
      }
      y += dy;
    }
    return areaPixels;
  }

  public calcArea(double _w, double _h, double _xmin, double _ymin, int _size, int _maxIteration) {
    w = _w;
    h = _h;
    xmin = _xmin;
    ymin = _ymin;
    size = _size;
    maxIteration = _maxIteration;

    xmax = xmin + w;
    ymax = ymin + h;
    dx = (xmax - xmin) / size;
    dy = (ymax - ymin) / size;
  }

  int mandelbrotJulia(double zr, double zi, 
    double cr, double ci, 
    int iteration) {

    for (int n=0; n<iteration; n++) {
      double tmp;
      tmp = zr*zr - zi*zi + cr;
      zi = 2d*zr*zi + ci;
      zr = tmp;
      if (zr*zr + zi*zi>4d) {
        return 0;
      }
    }
    return 1;
  }
}

boolean outOfCardioidCircles(double ix, double iy) {
  //カージオイド
  double cx = -ix+0.25d;
  double cy = iy;
  double a =0.5d;
  //円
  double cx0 = -1.d-ix;
  double cy0 = 0d-iy;
  double r0 = 0.25d;
  double cx1 = -1.309d-ix;
  double cy1 = 0d-iy;
  double r1 = 0.055d;
  double cx2 = -0.125d-ix;
  double cy2 = -0.744d-iy;
  double r2 = 0.09d;
  double cx3 = -0.125d-ix;
  double cy3 = 0.744d-iy;
  double r3 = 0.09d;
  if (
    (cx*cx+cy*cy)*(cx*cx+cy*cy-2.*a*cx)-a*a*cy*cy < 0d ||
    cx0*cx0+cy0*cy0 < r0*r0 ||
    cx1*cx1+cy1*cy1 < r1*r1 ||
    cx2*cx2+cy2*cy2 < r2*r2 ||
    cx3*cx3+cy3*cy3 < r3*r3
    ) {
    return false;
  } else {
    return true;
  }
}