Pythonでガンマ関数とゼータ関数をプロットするためのメモ書き（複素数対応版）

複素数対応版を探すのに苦労したのでメモ。
精度については調査していない。

# 複素数の使い方

虚数単位として`j`を付けるだけ。`.real`と`.imag`で実部、虚部を取り出せる。

```python
1+2j
```

https://qiita.com/shuhei/items/f5cf6c83fcfb5dd24c2d

# ガンマ関数

scipy.specialの[`gamma`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gamma.html#scipy.special.gamma)を使う。ガンマ関数の逆数[`rgamma`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.rgamma.html#scipy.special.rgamma)も用意されている。
[array_like型](https://numpy.org/devdocs/glossary.html#term-array_like)を引数に取れる。

https://qiita.com/sci_Haru/items/a699fc6b86fa3ef65a28

[`mpmath.gamma`](https://mpmath.org/doc/current/functions/gamma.html#gamma-function)もある。

# (リーマン)ゼータ関数

[`mpmath.zeta`](https://mpmath.org/doc/current/functions/zeta.html)を使えば複素数領域の計算が可能。

なお、scipy.specialの[`zeta`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.zeta.html), [`zetac`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.zetac.html)関数は記事作成時点(2023/1/8)では複素数非対応。

https://python.atelierkobato.com/zeta/

## インストール
```
pip install mpmath
```

# プロット

`matplotlib.pyplot`を使う。

https://qiita.com/skotaro/items/08dc0b8c5704c94eafb9

https://python.atelierkobato.com/matplotlib/

https://jakevdp.github.io/PythonDataScienceHandbook/04.12-three-dimensional-plotting.html

# ソースコード

## scipyのガンマ関数

表示の都合で値`6.0`で上限を設定しています。

```py
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import *

Z_LIMIT = 6.0

def f(x, y):
    return np.clip(abs(gamma(x + y*1j)),0,Z_LIMIT)

x = np.linspace(-4, 4, 160)
y = np.linspace(-2, 2, 80)

X, Y = np.meshgrid(x, y)
Z = f(X, Y)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(elev=20, azim=225)
ax.set_zlim3d(top=Z_LIMIT)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis', edgecolor='none')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()
```

## mpmathのガンマ関数

mpmathのgamma関数の引数は、array_like型を取れないので、`for`文を使って配列要素ごとに計算して設定してます。

```py
from mpl_toolkits import mplot3d
import numpy as np
import matplotlib.pyplot as plt

import mpmath
mpmath.mp.dps = 20

Z_LIMIT = 6.0

def f(x, y):
    z = np.zeros((len(x), len(x[0])))
    for i in range(len(x)):
        for k in range(len(x[0])):
            z[i][k] = abs(mpmath.gamma(x[i][k] + y[i][k]*1j))
            if z[i][k] > Z_LIMIT :
                z[i][k] = Z_LIMIT
    return z

x = np.linspace(-4, 4, 160)
y = np.linspace(-2, 2, 80)

X, Y = np.meshgrid(x, y)
Z = f(X, Y)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(elev=20, azim=225)
ax.set_zlim3d(top=Z_LIMIT)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='viridis', edgecolor='none')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()
```



# ガンマ関数（の絶対値）のプロット結果 

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/438183/9222b50f-0145-918f-56c5-a02121c3673c.png)

scipyのgammaとmpmath(dps=20設定時)のgammaの差（厳密には`abs`で絶対値取ったあとの差なので、参考程度に。）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/438183/fa9f14a0-2899-50c5-9053-d7987a1fe4f8.png)


# 参考サイト（数学）

https://en.wikipedia.org/wiki/Gamma_function

https://en.wikipedia.org/wiki/Riemann_zeta_function