筆者はふだん Ruby ばかり書いていて Rust は初心者だが，自分のレベルに合った何か面白いプログラムが書きたい。
以前書いた「Rust でごった煮マシン」のような，視覚に訴えるような何か。Ruby だと計算量的に辛いようなもの。
よしっ，複素関数の可視化だ。一変数の複素関数を可視化しよう。
複素関数の可視化って？
一変数の実関数1の可視化は，中学校の数学で「関数のグラフ」として出てくる。
二つの変数の関係を 2 次元の紙面・黒板に表すのは容易だ。
しかし，複素関数2となるとちょっと工夫が要る。複素数自体が 2 次元的な広がりを持つため，実関数と同じような考え方でグラフを描こうとすると，4 次元空間上に図示しなければならないことになる。
基本アイデア
以下のように考えた。
$z$ を独立変数とし，表現したい複素関数を $f$ とする。
まず，関数 $f$ の値域3として複素平面をそのまま使う。
その複素平面上の各点にどうにかして $f(z)$ の値を表現する。
1 点つまり 1 ピクセルで一つの複素数を表すとすれば，色しかないだろう。
幸い，色は 3 次元的なので4，2 次元的な複素数を表すことができそうだ。
では点 $z$ における値 $f(z)$ にどのような色を割り当てるか。
偏角を色相に対応させ，絶対値を明度に対応させよう。
色相はそもそも円環状になっているので，都合がよい。
絶対値と明度の対応付けだが，$0$ 付近も $\infty$ 付近もいい感じに階調が変化するようにしたい。
そこで，複素数のリーマン球面表示を使うことにする。リーマン球面の垂直方向の座標値を明度に対応させるのだ。
すると，$0$ は真っ黒，$\infty$ は真っ白になる。絶対値が $1$ の複素数はそのちょうど中間の明度だ5。
おそらく，私が考案するまでもなく，きっと，Mathematica だの R だのといったソフトで昔から似たようなことをやっているのだろうと思う。車輪の再発明に違いない。
いずれにしてもこの方式の重大な欠陥は，色覚障害の無い方を対象にしていることである。障害というと特殊な印象を持つかもしれないが，典型的でない色覚を持つ方は人口の数 % いる6。いずれこの欠陥を克服する表示方法に挑戦したい。アイデアはある。
さて，これがどのような色付けになっているかは，$f(z) = z$ という恒等関数に適用してみればよい。以下のようになる。
中心（$z=0$）付近は絶対値が小さいので暗い。
中心の周りをぐるっと回ると色相が赤→橙→黄→緑→青緑→青→青紫→赤紫のように変化する。
軸に付いた目盛りは一目盛りが 1 である。
しかし，んー，なんというか，意図したとおりには表示されているのだが，これが分かりやすいかというと，どうもね。
量子化
なめらかに変化する階調は捉えにくいのだ。
そこで，量子化（quantization）ということを行う。いや，この用語が妥当かどうかはよく分からない。離散化（discretization）のほうがいいのかな。
要するに，複素数の絶対値と偏角それぞれについて，階段関数7を使って飛び飛びの値にしてやり（これが量子化），その値に基づいて色を決定するのだ。
このようにすると，色の境目のところがちょうど等値線8のようになって，値の変化が摑みやすくなる。
このことは，だいぶ前に書いた以下の記事で実証済みだ。
ミンコフスキー距離を絵にしてみた - Qiita
絶対値を 16 分割，偏角を 12 分割すると以下のようになる。
絶対値と偏角の分割数をいろいろに変えてさまざまな複素関数を表示してみたところ，下記の設定がまあよさそう，というところに落ち着いた。
- 絶対値は 12 分割
- 偏角は 4 分割
つまり以下のとおりである。
「いくら何でも偏角の分割が粗すぎる」と思われるだろう。しかし，これでいい。このくらい粗くないと $f(z) = z^3$ みたいな関数はワケがわからなくなってしまうのだ。
偏角の分割数は 4 の倍数がよい。そのほうが実部・虚部の大きさが摑みやすいのだ。8 分割も試したが，やや複雑な関数になると，やはり 4 分割が把握しやすい。
偏角 0°，90°，180°，270° の色を何と呼ぶか微妙だが，この記事では仮に「赤」「黄土色」「青緑」「群青」とでも呼んでおく。
赤の部分は，虚部が実部より小さく，実部が正と特徴づけられる。同じようなことが他の三色にも言える。これが関数を把握するうえで重要である。
コード
では Rust のコードを。
まず
cargo new complex_function
とでもやってプロジェクトを作り，Cargo.toml に
[dependencies]
image = "0.23.9"
num-complex = "0.3"
palette = "0.5"
を追記。
image は画像を生成する。num-complex は複素数の演算を行う。palette は LCh 色空間と RGB の変換に用いる。
src/lib.rs
ComplexCanvas は，画像のサイズや複素平面上の表示する範囲などを保持し，画像化するためのもの。
extern crate image;
extern crate num_complex;
extern crate palette;
use num_complex::Complex64;
use palette::{Lch, Srgb};
use std::f64::consts::PI;
pub struct ComplexCanvasParameter {
pub width: u32, // 画像の水平方向のピクセル数
pub height: u32, // 画像の垂直方向のピクセル数
pub origin: (u32, u32), // 原点（z=0）の位置
pub unit: u32, // z=0 から z=1 までのピクセル数
}
pub struct ComplexCanvas {
width: u32,
height: u32,
origin: (u32, u32),
unit: u32,
re_range: (f64, f64), // 画像の左端と右端の実部
re_span: f64, // その差
im_range: (f64, f64), // 画像の下端と上端の虚部
im_span: f64, // その差
}
impl ComplexCanvas {
pub fn new(params: ComplexCanvasParameter) -> ComplexCanvas {
let unit_f64 = params.unit as f64;
let re_range = (
-(params.origin.0 as f64 / unit_f64),
((params.width - params.origin.0) as f64) / unit_f64,
);
let im_range = (
-(((params.height - params.origin.1) as f64) / unit_f64),
params.origin.1 as f64 / unit_f64,
);
ComplexCanvas {
width: params.width,
height: params.height,
origin: params.origin,
unit: params.unit,
re_range: re_range,
re_span: re_range.1 - re_range.0,
im_range: im_range,
im_span: im_range.1 - im_range.0,
}
}
// ピクセルのインデックスに対応する複素数
fn pix2z(&self, px: u32, py: u32) -> Complex64 {
Complex64::new(
self.re_range.0 + (px as f64 + 0.5) * self.re_span / self.width as f64,
self.im_range.1 - (py as f64 + 0.5) * self.im_span / self.height as f64,
)
}
// 与えられた複素関数を描画し，与えられたパスの画像ファイルに書き出す
// 複素関数はクロージャーで与える
pub fn draw<F>(&self, f: F, pathname: &str)
where
F: Fn(Complex64) -> Complex64,
{
let mut imgbuf = image::ImageBuffer::new(self.width, self.height);
let mut z: Complex64;
let (ox, oy) = (self.origin.0, self.origin.1);
// 関数を描画
for (ix, iy, pixel) in imgbuf.enumerate_pixels_mut() {
z = self.pix2z(ix, iy);
*pixel = z2color(&f(z));
}
// 座標軸を描画
const COORDINATE_COLOR: image::Rgb<u8> = image::Rgb([150, 150, 150]);
for ix in 0..self.width {
imgbuf.put_pixel(ix, oy - 1, COORDINATE_COLOR);
imgbuf.put_pixel(ix, oy, COORDINATE_COLOR);
}
for iy in 0..self.height {
imgbuf.put_pixel(ox - 1, iy, COORDINATE_COLOR);
imgbuf.put_pixel(ox, iy, COORDINATE_COLOR);
}
let w = self.unit / 10;
let mut ix = ox % self.unit;
while ix < self.width {
for iy in oy - w..oy + w {
imgbuf.put_pixel(ix - 1, iy, COORDINATE_COLOR);
imgbuf.put_pixel(ix, iy, COORDINATE_COLOR);
}
ix += self.unit;
}
let mut iy = oy % self.unit;
while iy < self.height {
for ix in ox - w..ox + w {
imgbuf.put_pixel(ix, iy - 1, COORDINATE_COLOR);
imgbuf.put_pixel(ix, iy, COORDINATE_COLOR);
}
iy += self.unit;
}
// 画像ファイルに書き出す
imgbuf.save(pathname).unwrap();
}
}
// 階段関数
// 0, max/n, 2max/n, 3max/n, ..., max を中心とする
// 幅 max/n の区間の中央の値を取る
fn step(v: f64, max: f64, n: u32) -> f64 {
(v / max * n as f64 + 0.5).floor() / (n as f64) * max
}
// 正の剰余
fn modulo(a: f64, b: f64) -> f64 {
((a % b) + b) % b
}
# [test]
fn test_modulo() {
assert!(modulo(-0.1, 2.0) == 1.9);
assert!(modulo(-PI * 0.5, 2.0 * PI) - PI * 1.5 < 0.00001);
}
# [test]
fn test_complex_arg() {
let z: Complex64 = Complex64::new(0.0, -1.0);
assert!((z.arg() - PI * 0.5) < 0.0001);
}
// 無限大を含む任意の複素数に RGB 色を割り当てる
// 0 が黒，∞ が白
// 偏角が色相（0 のとき赤，2π/3 のとき緑，4π/3 のとき青）
fn z2color(&z: &Complex64) -> image::Rgb<u8> {
let abs2 = z.norm_sqr();
// 偏角を 0 以上 2π 未満の範囲に
let arg = modulo(z.arg(), 2.0 * PI);
// 量子化された色相
let quantized_hue = step(arg, 2.0 * PI, 4);
// 輝度（0 以上 1 以下）
// Riemann 球面の垂直座標に基づいて決定
let luminance: f64 = abs2 / (1.0 + abs2);
// 量子化された輝度
// 0〜1 よりも若干内側に入れる
let quantized_luminance = (step(luminance, 1.0, 12) - 0.5) * 0.8 + 0.5;
// 彩度（0 以上 1 以下）
// 輝度に応じて彩度を抑える
let quantized_chroma = if quantized_luminance < 0.72 {
quantized_luminance / 0.72 * 0.38
} else {
(1.0 - quantized_luminance) / 0.28 * 0.38
};
// LCh の色を作る
let lch = Lch::new(
(quantized_luminance * 100.0) as f32,
(quantized_chroma * 100.0) as f32,
(quantized_hue as f32) / (PI as f32) * 180.0,
);
// sRGB に変換する
// その値をそのまま RGB 値とする
let srgb: Srgb = lch.into();
let (r, g, b) = srgb.into_components();
image::Rgb([(r * 255.0) as u8, (g * 255.0) as u8, (b * 255.0) as u8])
}
なんか座標軸の描画に行数を喰っているような。もっと簡潔に書けませんかね。
src/main.rs
こんなふうに使う。
まずは正弦関数（sine）を描かせてみよう。
type C = num_complex::Complex64;
mod lib;
use lib::{ComplexCanvas, ComplexCanvasParameter};
fn main() {
let complex_canvas = ComplexCanvas::new(ComplexCanvasParameter {
width: 320,
height: 320,
origin: (160, 160),
unit: 50,
});
complex_canvas.draw(&|z: C| z.sin(), "images/sin.png");
}
冒頭の
type C = num_complex::Complex64;
は Complex64 の別名を C としている。こんなことしなくていいのだけど，後ででてくる関数を簡潔に書きたいがため。
ComplexCanvas を作って，その draw 関数を呼ぶだけ。
draw の第一引数に表示させたい複素関数をクロージャーとして与える。
第二引数は画像のパス。ここでは，images というディレクトリーが存在することを仮定しているので，プログラムを実行する前に手で作っておこう。
さて，問題のクロージャーは
|z: C| z.sin()
という形をしている。| | の中に引数を書く。引数の型も書く必要があるので，|z: C| としている。この C はさきほど作った Complex64 の別名である。
返り値の型は省略している。
今の場合，クロージャーの中身は単一の式なので，{ } で囲まず z.sin() を裸で書いている。
draw の第一引数はクロージャーへの参照なので，& を付けて渡している。
これで以下のような絵が描ける。
$f(z) = \sin z$
以下のようなことが容易に見て取れる。
- $k\pi ; (k \in \boldsymbol Z)$ に零点を持つ
- 実軸方向に周期的（周期 $2\pi$）
- 実軸から離れるほど絶対値が大きくなる
- $\mathrm{Re}, z = (2k+1)\pi ;(k\in \boldsymbol Z)$ では実軸に関して対称（$f(\bar z)=f(z)$）
- $\mathrm{Re}, z = 2k\pi ;(k\in \boldsymbol Z)$ では実軸に関して反対称（$f(\bar z)=-f(z)$）
ところで $f(z)=z$ の図を再掲するので見比べてほしい。
$f(z) = z$
原点（$z=0$）付近がよく似ている。
これは，$f(z)=\sin z$ のマクローリン展開（原点におけるテイラー展開）が
\sin z = z - \frac{1}{3!}z^3 + \frac{1}{5!} z^5 + \cdots
であって，原点付近の支配項が $z$ であることを反映している。
ギャラリー
ではさまざまな関数をご覧いただこう。
f(z) = 1/z
逆数はどうなるのか。
complex_canvas.draw(&|z: C| z.inv(), "images/inv.png");
のように書けばよい。
Complex64 には逆数を返すメソッド inv() があるので，C::new(1.0, 0.0) / z なんて書かなくていい。
$f(z) = 1/z$
原点で無限大 $\infty$ になっている
三たび $f(z) = z$ の図を掲載する。見比べてみよう。
$f(z) = z$
絶対値の等値線が完全に同じ形をしている。これはリーマン球面の垂直座標に基づいて輝度を決定したことに由来する。
ある意味で $0$ と $\infty$ が対等に扱えていると言えるのではないか。このことは次に掲げる正接関数でも見て取れる。
f(z) = tan(z)
正接関数（tangent）にいこう。
complex_canvas.draw(&|z: C| z.tan(), "images/tan.png");
とすればいいね。
$f(z) = \tan z$
- 実軸方向に周期 $\pi$ を持つ
- 実軸付近は派手だが，実軸から離れるとのっぺりしている
- $k\pi ;(k\in \boldsymbol Z)$ に零点を持つ
- $(k + \frac12)\pi ;(k\in \boldsymbol Z)$ に無限大点9を持つ
- 零点付近は $f(z)=z$ の原点付近に似ている
- 無限大点付近は $f(z)=-1/z$ の原点付近に似ている。
f(z) = z^n
では $f(z) = z^n$ の形の冪を見よう。
以下のように書く。
for n in 1..=4 {
complex_canvas.draw(&|z: C| z.powi(n), &format!("images/z^{}.png", n));
}
クロージャーには外の変数が見えるので，
|z: C| z.powi(n)
のように書いていいわけだ。
こうなる。
$f(z) = z$
これ↑は何度も再掲しているやつ。
$f(z) = z^2$
$z$ が原点周りを 1 周回れば，$f(z)$ は 2 周回ることがハッキリと分かる。
$f(z) = z^3$
3 周回っている。うむ。
$f(z) = z^4$
えっと，横尾忠則さんの作品 を連想するのは私だけだろうか。
f(z) = z^2 + 1
次はごく簡単な多項式関数 $f(z) = z^2 + 1$。
こう書く。
complex_canvas.draw(&|z: C| z * z + 1.0, "images/z^2+1.png");
んー，
z * z + C::new(1.0, 0.0)
でなくて
z * z + 1.0
でもエラーにならないのはなんでだろう？ 誰か教えて。
$f(z) = z^2 + 1$
この関数は $f(z) = (x+i)(x-i)$ と因数分解でき，$\pm i$ に零点を持つことが容易に分かるが，上の図からもそれが読み取れる。
零点の付近が $f(z)=z$（の偏角をずらしたもの）に似ていることも何となく分かる。
実軸上では正の実数値を取ることもうかがえる。
余談だが，実関数 $f(x) = x^2 + 1$ は正の値しか取らず，2 次方程式 $f(x) = 0$ は解を持たない。
しかし，これを複素関数に拡張すると $f$ は 2 個の零点を持ち，それゆえ方程式 $f(z)=0$ は 2 個の解を持つ。
もう一つ重要なことは，原点から遠く離れたところでは $f(z) = z^2$ と同じように振る舞うということ。
確かにそんなふうに見える。
おわりに
Rust で 1 変数の複素関数を可視化した。
関数の性質がまあ分かるような絵が描けたと思う（これぞまさに自画自賛！）。
ごく簡単なプログラムだが，何かしら実用性のあるものが Rust で書けた。楽しい。
同じことを C でやろうとは思わない。C が悪いのではない。あれは玄人しか使っちゃいけない道具だからだ。
他にもいろんな関数を描いてみよう。一般の代数関数とか（←あまり意味が分からずに書いている）。
三角関数と双曲線関数を比べてみたりとか。
関数を外部から与えることはできないだろうか。
時間とともに変化する複素関数のアニメーションとかもやってみたい。
-
実関数は独立変数・従属変数ともに実数である関数 ↩
-
複素関数は独立変数・従属変数ともに複素数である関数 ↩
-
$z$ が取りうる値全体の空間 ↩
-
色が 3 自由度を持つ量であることを意味する。しかし，あとでも触れるが色覚は多様であり，色が 3 自由度を持つのは色覚上の多数者にとっての真実に過ぎない。 ↩
-
当初 HSV 色空間に基づいて考えていたため「明度」という言葉を使ったが，試行錯誤の末，実際には LCh 色空間を採用したので，「輝度」と呼ぶべきかもしれないが私にはよく分からない。LCh の L は luminance（輝度）なので，コード中のコメントでは「輝度」としておいた。 ↩
-
話が逸れるが，これが色覚障害などと呼ばれるのは，社会が典型的な色覚を持つ者の都合で作られているから，そうでない者が不利益を被るということだ。 ↩
-
ここでいう「階段関数」は，機械学習方面で使う段が一段しか無いやつではなく，
floor関数のようなもの。 ↩ -
地図の等高線（同じ高さの地点を結んでできる曲線）と同じように，何かの値が等しい点を結んでできる曲線。 ↩
-
そういう用語があるかは知らん。 ↩