Processingの複素数ライブラリ

　以前、上記の投稿をしたときには、processingの複素数ライブラリを見つけることができなかったんですが、見直したら、ちゃんとありました。javaのライブラリを探していたのが行けなかったのかもしれません。

ライブラリ名称「Complex Numbers」
制作者「Math Machine」
マニュアルは、GitHubからZipでダウンロード。HTMLで提供。

以下、過去の投稿同様に、確認していきます。

複素数の計算

二乗してみます。

i^2=-1

processing

import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.mul(i));

-1

C#のpowを使うと精度誤差が出たのに、processingだと問題なし。

processing

import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.pow(2));

-1

次に、sqrt命令。

\sqrt{i}=\frac{1}{\sqrt{2}}+\frac{1}{\sqrt{2}}i

processing

import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.sqrt());

完璧

0.7071067811865476+0.7071067811865475i

指数部に実数0.5を使う。

processing

import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.pow(0.5));

同じ結果

0.7071067811865476+0.7071067811865475i

虚数の虚数乗は？

i^i=e^{-\pi/2}=\frac{1}{\sqrt{e^\pi}}=0.20787957\cdots

processing

import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.pow(i));

0.20787957635076193

完璧である

計算速度は?

整数pow(2)は十分速いが、実数pow(2.2)になると、100倍ぐらい遅くなる。内部で場合分けして、最適化しているようである。

processing

import complexnumbers.*;

for (double y = -1.25d; y <= 1.25d; y += 2.5d/24.0d)
{
    for (double x = -2.0d; x <= 0.5d; x += 2.5d/64.0d)
    {
        Complex z = new Complex(0, 0);
        Complex c = new Complex(x, y);
        for (int n = 0; n < 256; n++)
        {
            z = z.pow(2);
            z = z.add(c);

            if (z.abs() > 2)
            {
                print(".");
                break;
            }
            if (n == 255)
            {
                print("*");
            }
        }
    }
    println("");
}

.................................................................
.................................................................
.................................................................
.................................................................
................................................*................
..............................................*****..............
..........................................*.*********.*..........
.......................................*********************.....
...................................*************************.....
..................................****************************...
......................********...*****************************...
....................************.****************************....
......*........*******************************************.......
....................************.****************************....
......................********...*****************************...
..................................****************************...
...................................*************************.....
.......................................*********************.....
..........................................*.*********.*..........
..............................................*****..............
................................................*................
.................................................................
.................................................................
.................................................................
.................................................................

.................................................................
.................................................................
.................................................................
.................................................................
......................................................*..........
.............................................*.*.....***.........
.........................*******.........***************.........
.....................*************.*.**********************......
........................*************************************....
........................*.************************************...
................................******************************...
..............................********************************...
............................*******************************......
..............................********************************...
................................******************************...
........................*.************************************...
........................*************************************....
.....................*************.*.**********************......
.........................*******.........***************.........
.............................................*.*.....***.........
......................................................*..........
.................................................................
.................................................................
.................................................................
.................................................................

import complexnumbers.*;

size(512, 512);
loadPixels();
int px, py;
py=0;
for (double y = -1.25d; y < 1.25d; y += 2.5d/512.0d)
{
  px=0;
  for (double x = -2.0d; x < 0.5d; x += 2.5d/512.0d)
  {
    Complex z = new Complex(0, 0);
    Complex c = new Complex(x, y);
    for (int n = 0; n < 256; n++)
    {
      z = z.pow(2.2);
      z = z.add(c);

      if (z.abs() > 2)
      {
        pixels[(int)(py*width+px)] = color(n*10%256);
        break;
      }
      if (n == 255)
      {
        pixels[(int)(py*width+px)] = color(0, 0, 255);
      }
    }
    px++;
  }
  py++;
}
updatePixels();

512x512の繰り返し最大256で、
時間は8秒ぐらいです。
マルチスレッド化したくなる。

z_{n+1}=z_n^2+c

上の式について、2乗のところを7まで0.01刻みでアニメーション

同じく2から0まで、0.01刻みでアニメーション

複素数のn乗(nは実数)

このライブラリは複素数の複素数乗が計算できる。複素数の実数乗に限定して、関数を書いたら、体感２倍ぐらいに速くなった。10倍ぐらいにならないかなと思ったが、大した効果はなかった。

ド・モアブルの定理より

z^n=|z|^n(cos(n\theta)+isin(n\theta))

processing

import complexnumbers.*;

float f=2.0;

void setup() {
  size(256, 256);
}

void draw() {
  loadPixels();
  int py=0;
  for (double y = -1.25d; y < 1.25d; y += 2.5d/256.0d)
  {
    int px=0;
    for (double x = -2.0d; x < 0.5d; x += 2.5d/256.0d)
    {
      Complex z = new Complex(0, 0);
      Complex c = new Complex(x, y);
      for (int n = 0; n < 256; n++)
      {
        //z = z.pow(f);
        z = myPow(z, f);
        z = z.add(c);

        if (z.abs() > 2)
        {
          pixels[(int)(py*width+px)] = color(n*10%256);
          n=0;
          break;
        }
        if (n == 255)
        {
          pixels[(int)(py*width+px)] = color(0, 0, 255);
        }
      }
      px++;
    }
    py++;
  }
  updatePixels();
  text(str(f), 5, 10);
  //saveFrame("mandelFP#####.png");
  f+=0.01;
}

Complex myPow(Complex c, double n) {
  double r = c.abs();
  double th = c.arg();
  double prn = Math.pow(r, n);
  double re = prn*(Math.cos(n*th));
  double im = prn*(Math.sin(n*th));
  return new Complex(re, im);
}

参考

ド・モアブルの公式