近年、多くのプログラミング言語で「複素数 (complex number)」がサポートされるようになりました。たとえば、Ｃ言語では 1999 年に発行された規格 C99 で複素数型が導入されました。Go 言語や Python でも複素数型をサポートしていますし、複素数用の標準ライブラリを用意している言語 (C++, C#, F#, Ruby, Haskell など) も多くあります。
他の言語では、FORTRAN や Common Lisp が昔から複素数型をサポートしています。Common Lisp の場合、基本的な数学関数でも複素数を適用できるのであれば、引数に複素数を渡して計算することができます。.NET (C#, F#) の場合、ライブラリ System.Numerics に複素数を表す構造体 Complex が用意されているので、簡単に複素数を扱うことができます。今回は F# の複素数についてまとめてみました。
.NET の数はプリミティブな型で大きく分けると、整数 (int, long など) と実数 (float32, float など) の 2 種類があります。System.Numerics に用意されている構造体 BigInteger を使うと、任意の桁の整数 (多倍長整数) を扱うことができます。F# では bigint という別名の型が用意されています。
実数は浮動小数点数 (floating point number) として表現されます。浮動小数点数には IEEE 754 という標準仕様があり、近代的なプログラミング言語のほとんどは、IEEE 754 に準拠した浮動小数点数をサポートしています。浮動小数点数はすべての小数を正確に表現することはできません。このため、実数は近似的 (不正確) な値になります。
IEEE 754 には通常の数値以外にも、負のゼロ (-0.0)、正負の無限大 (∞, -∞)、NaN (Not a Number, 非数) といった値が定義されています。これらの値は .NET でも取り扱うことができます。F# の REPL (dotnet fsi) では、負のゼロを -0.0、正負の無限大を infinity と -infinity 、NaN を nan と表示します。
一般に、無限大は値のオーバーフロー、ゼロ除算 (数値 / 0.0)、数学関数の計算結果 (たとえば log(0.0)) などで発生します。なお、浮動小数点数のゼロ除算でエラー (例外) を送出する処理系 (たとえば Python など) もあります。float 型の無限大は .NET のライブラリ System.Double の定数 PositiveInfinity, NegativeInfinity に定義されています。
簡単な実行例を示します。
> open System;; > 1e308;; val it: float = 1e+308 > 1e308 * 2.0;; val it: float = infinity > -1e308;; val it: float = -1e+308 > -1e308 * 2.0;; val it: float = -infinity > 1.0 / 0.0;; val it: float = infinity > -1.0 / 0.0;; val it: float = -infinity > log(0.0);; val it: float = -infinity > Double.PositiveInfinity;; val it: float = infinity > Double.NegativeInfinity;; val it: float = -infinity
.NET の場合、float 型の無限大の判定はライブラリ System.Double に用意されている関数を使います。
> Double.IsInfinity(1.0 / 0.0);; val it: bool = true > Double.IsInfinity(-1.0 / 0.0);; val it: bool = true > Double.IsInfinity(-1.0);; val it: bool = false > Double.IsPositiveInfinity(1.0 / 0.0);; val it: bool = true > Double.IsPositiveInfinity(-1.0 / 0.0);; val it: bool = false > Double.IsNegativeInfinity(-1.0 / 0.0);; val it: bool = true > Double.IsNegativeInfinity(1.0 / 0.0);; val it: bool = false
無限大は他の数値と比較したり演算することもできますが、結果が NaN になることもあります。
> let a = 1.0 / 0.0;; val a: float = infinity > let b = -1.0 /0.0;; val b: float = -infinity > a < b;; val it: bool = false > a > b;; val it: bool = true > a = b;; val it: bool = false > a = a;; val it: bool = true > a < 0;; val it: bool = false > a > 0;; val it: bool = true > b < 0;; val it: bool = true > b > 0;; val it: bool = false > a + 100.0;; val it: float = infinity > a - 100.0;; val it: float = infinity > a * 100.0;; val it: float = infinity > a / 100.0;; val it: float = infinity > a - a;; val it: float = nan > a / a;; val it: float = nan > a + b;; val it: float = nan > a * 0.0;; val it: float = nan
負のゼロ (-0.0) は、計算結果が負の極めて小さな値でアンダーフローになったとき発生します。また、正の値を負の無限大で除算する、負の値を正の無限大で除算する、負の値と 0.0 を乗算しても -0.0 が得られます。
> -1e-323;; val it: float = -9.881312917e-324 > -1e-323 / 2.0;; val it: float = -4.940656458e-324 > -1e-323 / 4.0;; val it: float = -0.0 > -1.0 / (1.0 / 0.0);; val it: float = -0.0 > 1.0 / (-1.0 / 0.0);; val it: float = -0.0 > -1.0 * 0.0;; val it: float = -0.0
標準 (IEEE 754) では、演算子 (= など) による 0.0 と -0.0 の比較は等しいと判定されます。.NET の場合、-0.0 は関数 Double.IsNegative で判別することができます。また、任意の正の数を -0.0 で除算すると負の無限大になるので、それを使って判別することもできます。
> let a = 0.0;; val a: float = 0.0 > let b = -0.0;; val b: float = -0.0 > a = b;; val it: bool = true > Double.IsNegative(a);; val it: bool = false > Double.IsNegative(b);; val it: bool = true > 1.0 / a;; val it: float = infinity > 1.0 / b;; val it: float = -infinity > (1.0 / a) < 0 ;; val it: bool = false > (1.0 / b) < 0 ;; val it: bool = true
-0.0 は数学関数 (たとえば atan2 など) や複素数の演算処理などで使われます。
> atan2;; val it: (float -> float -> float) = <fun:it@83> > atan2 0.0 (-1.0);; val it: float = 3.141592654 > atan2 (-0.0) (-1.0);; val it: float = -3.141592654 > sqrt(0.0);; val it: float = 0.0 > sqrt(-0.0);; val it: float = -0.0
NaN は数ではないことを表す特別な値 (非数) です。一般的には 0.0 / 0.0 といった不正な演算を行うと、その結果は NaN になります。たとえば、負数の平方根は虚数になりますが、関数 sqrt に負数を与えると浮動小数点数では表現できないので NaN になります。
> 0.0 / 0.0;; val it: float = nan > sqrt(-1.0);; val it: float = nan > let a = 1.0 / 0.0;; val a: float = infinity > a / a;; val it: float = nan > Double.IsNaN(a / a);; val it: bool = true > Double.IsNaN(-0.0);; val it: bool = false
.NET の場合、NaN は関数 Double.IsNaN で判別することができます。
数学では複素数 z を \(x + yi\) と表記します。x を実部、y を虚部、i を虚数単位といいます。虚数単位は 2 乗すると -1 になる数です。実部と虚部の 2 つの数値を格納するデータ構造を用意すれば、プログラミング言語でも複素数を表すことができます。
.NET で複素数を利用するときはライブラリ System.Numerics を使います。.NET の複素数は構造体で、データ型は Complex になります。複素数はコンストラクタ Complex で生成します。
let z = new Complex(double x, double y)
第 1 引数 x が実部、第 2 引数 y が虚部になります。ここで「型略称」という機能を使うと、複素数の操作が簡単になります。
type 別名 = 既存の型
型略称は既存の型に別の名前を付ける機能です。.NET の場合、多倍長整数を表す構造体は System.Numerics.BigInteger ですが、F# では型略称により bigint という別名が付けられています。本ページでは複素数に complex という別名を付けることにしましょう。
> type complex = System.Numerics.Complex;;
type complex = Numerics.Complex
> let a = complex(1.0, 2.0);;
val a: complex = <1; 2>
> a;;
val it: complex = <1; 2> {Imaginary = 2.0;
Magnitude = 2.236067977;
Phase = 1.107148718;
Real = 1.0;}
> let b = complex(3.0, 4.0);;
val b: complex = <3; 4>
> b;;
val it: complex = <3; 4> {Imaginary = 4.0;
Magnitude = 5.0;
Phase = 0.927295218;
Real = 3.0;}
C# の場合、クラスや構造体は new 型名(...) でインスタンスや値を生成します。本ページ (F# 入門) は、インスタンスを生成するとき new を使ってきましたが、F# は new を省略して 型名(...) だけでもクラスのインスタンスや構造体の値を生成することができます。型略称の場合も同様です。new を付けてもいいですが、付けなくても複素数を生成することができます。
複素数 z の実部はプロパティ Real で、虚部はプロパティ Imaginary で取得することができます。複素数の虚部の符号を反転することを「複素共役」といいます。複素共役は関数 Conjugate で求めることができます。
簡単な例を示しましょう。
> a.Real;;
val it: float = 1.0
> a.Imaginary;;
val it: float = 2.0
> complex.Conjugate(a);;
val it: Numerics.Complex = <1; -2> {Imaginary = -2.0;
Magnitude = 2.236067977;
Phase = -1.107148718;
Real = 1.0;}
このほかに、Complex には次の定数が定義されています。
> complex.ImaginaryOne;;
val it: Numerics.Complex = <0; 1> {Imaginary = 1.0;
Magnitude = 1.0;
Phase = 1.570796327;
Real = 0.0;}
> complex.One;;
val it: Numerics.Complex = <1; 0> {Imaginary = 0.0;
Magnitude = 1.0;
Phase = 0.0;
Real = 1.0;}
> complex.Zero;;
val it: Numerics.Complex = <0; 0> {Imaginary = 0.0;
Magnitude = 0.0;
Phase = 0.0;
Real = 0.0;}
> complex.Infinity;;
val it: Numerics.Complex = <∞; ∞> {Imaginary = infinity;
Magnitude = infinity;
Phase = 0.7853981634;
Real = infinity;}
> complex.NaN;;
val it: Numerics.Complex = <NaN; NaN> {Imaginary = nan;
Magnitude = nan;
Phase = nan;
Real = nan;}
複素数は極形式 \(z = r (\cos \theta + i \sin \theta)\) で表すことができます。r を絶対値、\(\theta\) を偏角といいます。絶対値はプロパティ Magnitude で求めることができます。偏角は複素平面において正の実軸とベクトル (x, y) との角度を表します。偏角はプロパティ Phase で求めます。Phase の返り値 \(\theta\) は \(-\pi \leq \theta \leq \pi\) (\(\pi\) : 円周率) です。
極形式 \(z = r (\cos \theta + i \sin \theta)\) で複素数を生成するときは関数 FromPolarCoordinates を使います。
FromPolarCoordinates(magnitude: float, phase: float) => complex
第 1 引数が絶対値、第 2 引数が偏角 θ を表します。
複素数 z はオイラーの公式 \(e^{i\theta} = \cos \theta + i \sin \theta\) を使って \(z = r e^{i\theta}\) と書くこともできます。ここで、\(\theta\) に \(pi\) を代入するとオイラーの等式 \(e^{i\pi} + 1 = 0\) になります。また、下記のようにド・モアヴルの公式も導出することができます。
簡単な例を示しましょう。
> a;;
val it: complex = <1; 2> {Imaginary = 2.0;
Magnitude = 2.236067977;
Phase = 1.107148718;
Real = 1.0;}
> a.Magnitude;;
val it: float = 2.236067977
> a.Phase;;
val it: float = 1.107148718
> complex.FromPolarCoordinates(a.Magnitude, a.Phase);;
val it: Numerics.Complex = <1.0000000000000002; 2> {Imaginary = 2.0;
Magnitude = 2.236067977;
Phase = 1.107148718;
Real = 1.0;}
> complex(-1.0, 0.0).Phase;;
val it: float = 3.141592654
> complex(-1.0, 1.0).Phase;;
val it: float = 2.35619449
> complex(0.0, 1.0).Phase;;
val it: float = 1.570796327
> complex(1.0, 1.0).Phase;;
val it: float = 0.7853981634
> complex(1.0, 0.0).Phase;;
val it: float = 0.0
> complex(1.0, -1.0).Phase;;
val it: float = -0.7853981634
> complex(0.0, -1.0).Phase;;
val it: float = -1.570796327
> complex(-1.0, -1.0).Phase;;
val it: float = -2.35619449
> complex(-1.0, -0.0).Phase;;
val it: float = -3.141592654
.NET の場合、\(-1.0+0.0i\) の偏角 \(\theta\) は \(\pi\) になり、\(-1.0-0.0i\) の偏角は \(-\pi\) になります。ゼロと負のゼロを区別しないプログラミング言語では、偏角 \(\theta\) の範囲を \(-\pi \lt \theta \leq \pi\) に制限して、\(-1.0 + 0.0i \ (= -1.0 - 0.0j)\) の偏角を \(\pi\) とします。
複素数の四則演算は次のようになります。
これらの演算は F# の演算子 +, -, * , / で行うことができます。複素数の場合、大小の比較演算は使えませんが、等値の判定は演算子 = や <> で行うことができます。簡単な例を示しましょう。
> let a = complex(1.0, 2.0);;
val a: complex = <1; 2>
> let b = complex(3.0, 4.0);;
val b: complex = <3; 4>
> a + b;;
val it: Numerics.Complex = <4; 6> {Imaginary = 6.0;
Magnitude = 7.211102551;
Phase = 0.9827937232;
Real = 4.0;}
> a - b;;
val it: Numerics.Complex = <-2; -2> {Imaginary = -2.0;
Magnitude = 2.828427125;
Phase = -2.35619449;
Real = -2.0;}
> a * b;;
val it: Numerics.Complex = <-5; 10> {Imaginary = 10.0;
Magnitude = 11.18033989;
Phase = 2.034443936;
Real = -5.0;}
> a / b;;
val it: Numerics.Complex = <0.44; 0.08> {Imaginary = 0.08;
Magnitude = 0.4472135955;
Phase = 0.1798534998;
Real = 0.44;}
> a = complex(1.0, 2.0);;
val it: bool = true
> a = b;;
val it: bool = false
> a <> b;;
val it: bool = true
> b <> b;;
val it: bool = false
実部と虚部の値は無限大や NaN になることもあります。
> let c = complex(1e300, 1e300);;
val c: complex = <1E+300; 1E+300>
> c * c;;
val it: Numerics.Complex = <NaN; ∞> {Imaginary = infinity;
Magnitude = nan;
Phase = nan;
Real = nan;}
最後に F# のお勉強ということで、実際に複素数を実装してみましょう。今回は構造体を使うことにします。次のリストを見てください。
リスト : 複素数の定義 (complex.fsx)
[<Struct; CustomEquality; NoComparison>]
type Complex =
val Real: float
val Imag: float
// コンストラクタ
new(a: float, b: float) = { Real = a; Imag = b }
new(a: int, b: int) = { Real = float a; Imag = float b }
// メソッド
override this.ToString () =
sprintf "complex(%g, %g)" this.Real this.Imag
// 等値演算子
override this.Equals (other: obj) =
let that = other :?> Complex
this.Real = that.Real && this.Real = that.Real
// とても適当なので実際に使ってはいけない
override this.GetHashCode () =
(this.Real / this.Imag).GetHashCode()
名前は Complex としました。実部をフィールド変数 Real に、虚部を Imag にセットします。格納する数値の型は実数 (float) とします。複素数はコンストラクタ Complex() で生成します。引数 a が実部で、b が虚部です。a, b をフィールド変数 Real, Imag にセットします。引数が int の場合は float に変換してからセットします。
次に、メソッド ToString(), Equals(), GetHashCode() をオーバーライドします。複素数の場合、大小の比較演算はできませんが、等値の判定は演算子 = や <> で行うことができます。等値演算子の実装はメソッド Equals() をオーバーライドするだけです。このとき、GetHashCode() もオーバーライドしないとワーニングが表示されます。
Equals() と GetHashCode() だけをオーバーライドする場合、属性 CustomEquality と NoComparison の指定が必要になります。比較演算子を実装すると、CustomEquality と NoComparison の指定は不要になります。ちなみに、比較演算子はインターフェース Sysmte.IComparable のメソッド CompareTo() を実装するだけです。
次に、複素共役、絶対値、偏角を求めるメソッドを定義します。
リスト : 複素共役, 絶対値, 偏角
// 複素共役
member this.Conjugate() = Complex(this.Real, -this.Imag)
// 偏角
member this.Phase() = atan2 this.Imag this.Real
// 絶対値
member this.Magnitude() =
if this.Real = 0.0 then abs this.Imag
else if this.Imag = 0.0 then abs this.Real
else if abs this.Imag > abs this.Real then
let temp = this.Real / this.Imag
(abs this.Imag) * sqrt(1.0 + temp * temp)
else
let temp = this.Imag / this.Real
(abs this.Real) * sqrt(1.0 + temp * temp)
複素共役を求める Conjugate() と偏角を求める Phase() は簡単ですね。atan2 y x は直交座標においてベクトル (x, y) と x 軸との角度を求める関数です。角度 \(\theta\) の範囲は \(-\pi \leq \theta \leq \pi\) になります。簡単な例を示しましょう。
> atan2 0.0 1.0;; val it: float = 0.0 > atan2 1.0 1.0;; val it: float = 0.7853981634 > atan2 1.0 0.0;; val it: float = 1.570796327 > atan2 0.0 -1.0;; val it: float = 3.141592654 > atan2 -1.0 1.0;; val it: float = -0.7853981634 > atan2 -1.0 0.0;; val it: float = -1.570796327 > atan2 -1.0 -1.0;; val it: float = -2.35619449 > atan2 -0.0 -1.0;; val it: float = -3.141592654
\(y \geq 0\) の場合、atan2 の返り値は \(0 \leq \theta \leq \pi\) になり、y が負 (-0.0 も含む) の場合は \(-\pi \leq \theta \lt 0\) になります。
絶対値を求めるメソッド Magnitude は、定義 \(\sqrt{x^2 + y^2}\) をそのままプログラムすると二乗の計算でオーバーフローすることがあります。たとえば、1e300 + i 1e300 の絶対値を求めてみましょう。このとき、1e300 の二乗で無限大となります。
> 1e300 * 1e300;; val it: float = infinity
参考文献『Ｃ言語による最新アルゴリズム事典』によると、次のように場合分けすることで、\(x^2\) や \(y^2\) で生じ得る上位桁あふれを回避することができるそうです。
今回は参考文献『Ｃ言語による最新アルゴリズム事典』のプログラムを F# に移植しました。
それでは実際に試してみましょう。
> #load "complex.fsx";; ... 略 ... > open Complex;; > let a = Complex(1, 1);; val a: Complex = complex(1, 1) > let b = a.Conjugate();; val b: Complex = complex(1, -1) > let c = a.Conjugate().Conjugate();; val c: Complex = complex(1, 1) > Complex(1, 1).Magnitude();; val it: float = 1.414213562 > Complex(1, -1).Magnitude();; val it: float = 1.414213562 > Complex(1e300, 1e300).Magnitude();; val it: float = 1.414213562e+300 > Complex(1e301, 1e300).Magnitude();; val it: float = 1.004987562e+301 > Complex(1e300, 1e301).Magnitude();; val it: float = 1.004987562e+301 > Complex(1, 0).Phase();; val it: float = 0.0 > Complex(1, 1).Phase();; val it: float = 0.7853981634 > Complex(0, 1).Phase();; val it: float = 1.570796327 > Complex(-1, 1).Phase();; val it: float = 2.35619449 > Complex(-1, 0).Phase();; val it: float = 3.141592654 > Complex(1, -1).Phase();; val it: float = -0.7853981634 > Complex(0, -1).Phase();; val it: float = -1.570796327 > Complex(-1, -1).Phase();; val it: float = -2.35619449 > Complex(-1.0, -0.0).Phase();; val it: float = -3.141592654
複素数の四則演算は次のようになります。
リスト : 複素数の四則演算
static member ( + ) (x: Complex, y: Complex) =
Complex(x.Real + y.Real, x.Imag + y.Imag)
static member ( - ) (x: Complex, y: Complex) =
Complex(x.Real - y.Real, x.Imag - y.Imag)
static member ( * ) (x: Complex, y: Complex) =
let a = x.Real
let b = x.Imag
let c = y.Real
let d = y.Imag
Complex(a * c - b * d, b * c + a * d)
static member ( / ) (x: Complex, y: Complex) =
let a = x.Real
let b = x.Imag
let c = y.Real
let d = y.Imag
if abs c >= abs d then
let u = d / c
let v = c + d * u
Complex((a + b * u) / v, (b - a * u) / v)
else
let u = c / d
let v = c * u + d
Complex((a * u + b) / v, (b * u - a) / v)
除算の場合、絶対値の計算と同様にオーバーフローの対策が必要になります。今回は参考文献『Ｃ言語による最新アルゴリズム事典』のプログラムを F# に移植しました。
それでは実際に試してみましょう。
> let a = Complex(1, 2);; val a: Complex = complex(1, 2) > let b = Complex(3, 4);; val b: Complex = complex(3, 4) > let c = a + b;; val c: Complex = complex(4, 6) > let d = a - b;; val d: Complex = complex(-2, -2) > let e = a * b;; val e: Complex = complex(-5, 10) > let f = a / b;; val f: Complex = complex(0.44, 0.08) > let one = Complex(1, 0);; val one: Complex = complex(1, 0) > let g = one / Complex(1e300, 1e300);; val g: Complex = complex(5e-301, -5e-301) > let h = one / Complex(1e301, 1e300);; val h: Complex = complex(9.90099e-302, -9.90099e-303) > let i = one / Complex(1e300, 1e301);; val i: Complex = complex(9.90099e-303, -9.90099e-302)
//
// complex.fsx : 複素数
//
// Copyright (C) 2022 Makoto Hiroi
//
[<Struct; CustomEquality; NoComparison>]
type Complex =
val Real: float
val Imag: float
// コンストラクタ
new(a: float, b: float) = { Real = a; Imag = b }
new(a: int, b: int) = { Real = float a; Imag = float b }
// メソッド
override this.ToString () =
sprintf "complex(%g, %g)" this.Real this.Imag
// 等値演算子
override this.Equals (other: obj) =
let that = other :?> Complex
this.Real = that.Real && this.Real = that.Real
// とても適当なので実際に使ってはいけない
override this.GetHashCode () =
(this.Real / this.Imag).GetHashCode()
// 複素共役
member this.Conjugate() = Complex(this.Real, -this.Imag)
// 偏角
member this.Phase() = atan2 this.Imag this.Real
// 絶対値
member this.Magnitude() =
if this.Real = 0.0 then abs this.Imag
else if this.Imag = 0.0 then abs this.Real
else if abs this.Imag > abs this.Real then
let temp = this.Real / this.Imag
(abs this.Imag) * sqrt(1.0 + temp * temp)
else
let temp = this.Imag / this.Real
(abs this.Real) * sqrt(1.0 + temp * temp)
// 四則演算
static member ( + ) (x: Complex, y: Complex) =
Complex(x.Real + y.Real, x.Imag + y.Imag)
static member ( - ) (x: Complex, y: Complex) =
Complex(x.Real - y.Real, x.Imag - y.Imag)
static member ( * ) (x: Complex, y: Complex) =
let a = x.Real
let b = x.Imag
let c = y.Real
let d = y.Imag
Complex(a * c - b * d, b * c + a * d)
static member ( / ) (x: Complex, y: Complex) =
let a = x.Real
let b = x.Imag
let c = y.Real
let d = y.Imag
if abs c >= abs d then
let u = d / c
let v = c + d * u
Complex((a + b * u) / v, (b - a * u) / v)
else
let u = c / d
let v = c * u + d
Complex((a * u + b) / v, (b * u - a) / v)