Processingの複素数ライブラリ

https://qiita.com/hextomino/items/13d74872625e52a1329a

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

https://processing.org/reference/libraries/

https://github.com/vmathmachine/Complex-Numbers/tree/latest

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

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

# 複素数の計算

二乗してみます。

```math
i^2=-1
```

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

```
-1
```

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

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

```
-1
```

次に、sqrt命令。

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

```java:processing
import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.sqrt());
```
完璧
```
0.7071067811865476+0.7071067811865475i
```

指数部に実数0.5を使う。

```java:processing
import complexnumbers.*;
Complex i = new Complex(0,1);
print(i.pow(0.5));
```
同じ結果
```
0.7071067811865476+0.7071067811865475i
```

虚数の虚数乗は？実数になるか？

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

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

```
0.20787957635076193
```
完璧である。

参考までに、オイラーの公式より

```math
e^{i\theta}=\cos\theta+i\sin\theta
```
```math
\theta=\frac{\pi}{2}とすれば、
```
```math
e^{\frac{i\pi}{2}}=\cos\frac{\pi}{2}+i\sin\frac{\pi}{2}=0+i=i
```
ここで、両辺をi乗する。
```math
e^{\frac{i\pi}{2}i}=i^i
```
これを、まとめて
```math
e^{-\frac{\pi}{2}}=i^i
```


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

```java: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("");
}
```

```
.................................................................
.................................................................
.................................................................
.................................................................
................................................*................
..............................................*****..............
..........................................*.*********.*..........
.......................................*********************.....
...................................*************************.....
..................................****************************...
......................********...*****************************...
....................************.****************************....
......*........*******************************************.......
....................************.****************************....
......................********...*****************************...
..................................****************************...
...................................*************************.....
.......................................*********************.....
..........................................*.*********.*..........
..............................................*****..............
................................................*................
.................................................................
.................................................................
.................................................................
.................................................................

```

```
.................................................................
.................................................................
.................................................................
.................................................................
......................................................*..........
.............................................*.*.....***.........
.........................*******.........***************.........
.....................*************.*.**********************......
........................*************************************....
........................*.************************************...
................................******************************...
..............................********************************...
............................*******************************......
..............................********************************...
................................******************************...
........................*.************************************...
........................*************************************....
.....................*************.*.**********************......
.........................*******.........***************.........
.............................................*.*.....***.........
......................................................*..........
.................................................................
.................................................................
.................................................................
.................................................................

```

```java
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秒ぐらいです。
マルチスレッド化したくなる。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/16add766-ef02-d871-9615-00f9ad134f4b.png)

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

上の式について、2乗のところを7まで0.01刻みでアニメーション
![mandelFP.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/454a4107-7663-31dd-a073-f6cd3f51ebef.gif)
同じく2から0まで、0.01刻みでアニメーション
![mandelFP2_0.png.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/510257/496c3cd0-c9bf-ca1d-bc94-4a86aa60daee.gif)

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

ド・モアブルの定理より

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


```java: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);
}
```

# 参考

ド・モアブルの公式

https://en.wikipedia.org/wiki/De_Moivre%27s_formula

https://qiita.com/hextomino/items/062f19e628a15774c252