F# Programming

構造体

「構造体 (structer)」はクラスのようにユーザーが定義するデータ型の一つです。クラスは「参照型」ですが、構造体は「値型」になります。値型は変数に代入するときや関数の引数に渡すとき、その値がコピーされます。基本的な数値型は値型です。構造体も同様に、内部の値がコピーされます。

●構造体の定義

構造体は type を使って宣言します。

1. type 名前 = struct
     // フィールド変数
     val [ mutable ] [ access-modifier ] 変数名: 型名
     ...
     // コンストラクタ
     new(...) = 式
     ...
     // メソッド
     member 自己識別子.名前(引数: 型, ...) = 式
     ...
   end

2. [<Struct>]
   type 名前 =
     // フィールド変数
     ...
     // コンストラクタ
     ...
     // メソッド
     ...

構造体は 1. のように type 名前 = struct ... end で定義します。属性 [<Struct>] を指定すると、struct と end を省略した軽量構文を使うことができます。フィールド変数は val を使って宣言します。クラスとは異なり、let 束縛や do 束縛はありません。構造体はクラスや他の構造体を継承することはできません。構造体を継承してサブクラスを定義することもできません。ただし、インターフェースは継承することができます。

構造体の値はクラスと同様にコンストラクタで生成します。

1. new 型名(...)
2. 型名(...)

C# の場合、クラスや構造体は 1. のように new を使ってインスタンスや値を生成します。本ページ (F# 入門) は、インスタンスを生成するとき new を使ってきましたが、F# は new を省略して 2. のように 型名(...) だけでもクラスのインスタンスや構造体の値を生成することができます。

フィールド変数はコンストラクタで初期化します。構造体は引数なしのコンストラクタを定義することはできません。引数なしでコンストラクタを呼び出すと、フィールド変数はゼロ初期化 (数値であれば 0 や 0.0、参照型であれば null) されます。val のアクセス制御は access-modifier で指定することができます。省略した場合は public になります。

●簡単な使用例

簡単な使用例として、平面座標を表す型 Point を構造体で定義してみましょう。

リスト : 構造体の使用例

type Point = struct
  val x: float
  val y: float
  new(a: float, b: float) = { x = a; y = b }
  member this.distance(p: Point) =
    let dx = this.x - p.x
    let dy = this.y - p.y
    sqrt(dx * dx + dy * dy)
end

プログラムは簡単なので説明は割愛させていただきます。以下に実行例を示します。

> type Point = struct
-   val x: float
-   val y: float
-   new(a: float, b: float) = { x = a; y = b }
-   member this.distance(p: Point) =
-     let dx = this.x - p.x
-     let dy = this.y - p.y
-     sqrt(dx * dx + dy * dy)
- end;;
[<Struct>]
type Point =
  new: a: float * b: float -> Point
  val x: float
  val y: float
  member distance: p: Point -> float

> let a = Point();;
val a: Point = FSI_0002+Point

> a.x;;
val it: float = 0.0

> a.y;;
val it: float = 0.0

> let b = Point(1.0, 1.0);;
val b: Point = FSI_0002+Point

> b.x;;
val it: float = 1.0

> b.y;;
val it: float = 1.0

> a.distance(b);;
val it: float = 1.414213562

●ジェネリック

構造体はクラスと同様にジェネリックを使用することができます。簡単な例として、2 つの値を格納する構造体 Pair を定義してみましょう。

リスト : ジェネリック構造体

type Pair<'a, 'b> = struct
  val mutable fst: 'a
  val mutable snd: 'b
  new(x: 'a, y: 'b) = {fst = x; snd = y}
end

プログラムは簡単なので説明は割愛させていただきます。以下に実行例を示します。

> type Pair<'a, 'b> = struct
-   val mutable fst: 'a
-   val mutable snd: 'b
-   new(x: 'a, y: 'b) = {fst = x; snd = y}
- end;;
[<Struct>]
type Pair<'a,'b> =
  new: x: 'a * y: 'b -> Pair<'a,'b>
  val mutable fst: 'a
  val mutable snd: 'b

> let a = Pair<int,float>();;
val a: Pair<int,float> = FSI_0011+Pair`2[System.Int32,System.Double]

> a.fst;;
val it: int = 0

> a.snd;;
val it: float = 0.0

> let b = Pair<int,float>(1, 1.2345);;
val b: Pair<int,float> = FSI_0011+Pair`2[System.Int32,System.Double]

> b.fst;;
val it: int = 1

> b.snd;;
val it: float = 1.2345

> let mutable c = Pair<int,float>();;
val mutable c: Pair<int,float> = FSI_0011+Pair`2[System.Int32,System.Double]

> c.fst;;
val it: int = 0

> c.fst <- 10;;
val it: unit = ()

> c.fst;;
val it: int = 10

構造体でフィールド変数を書き換える場合、フィールド変数を mutable で宣言するだけではなく、構造体の値を格納する変数も mutable で宣言する必要があります。ご注意くださいませ。

●値型のレコード

F# のレコードはデフォルトだと「参照型」ですが、属性 [<Struct>] を指定すると「値型」に変更することができます。簡単な例として、構造体 Pair をレコードで書き直してみましょう。最初は参照型のレコードです。

> type PairR<'a, 'b> = { mutable Fst: 'a; mutable Snd: 'b };;
type PairR<'a,'b> =
  {
    mutable Fst: 'a
    mutable Snd: 'b
  }

> let a = {Fst = 1; Snd = 1.234};;
val a: PairR<int,float> = { Fst = 1
                            Snd = 1.234 }

> a.Fst;;
val it: int = 1

> a.Fst <- 10;;
val it: unit = ()

> a;;
val it: PairR<int,float> = { Fst = 10
                             Snd = 1.234 }

レコード PairR のフィールド変数 Fst, Snd は mutable なので値を書き換えることができます。この場合、PairR は参照型なので、レコードを格納している変数 a が mutable でなくても値を書き換えることができます。

次は値型のレコードです。

> [<Struct>] type PairR<'a, 'b> = { mutable Fst: 'a; mutable Snd: 'b };;
[<Struct>]
type PairR<'a,'b> =
  {
    mutable Fst: 'a
    mutable Snd: 'b
  }

> let a = {Fst = 1; Snd = 1.234};;
val a: PairR<int,float> = { Fst = 1
                            Snd = 1.234 }

> a.Fst;;
val it: int = 1

> a.Fst <- 10;;
=> エラー

> let mutable b = {Fst = 123; Snd = 123.456};;
val mutable b: PairR<int,float> = { Fst = 123
                                    Snd = 123.456 }

> b.Fst;;
val it: int = 123

> b.Fst <- 456;;
val it: unit = ()

> b.Fst;;
val it: int = 456

属性 [<Struct>] を指定すると、PairR は値型のレコードになります。変数 a に格納されたレコードのフィールド変数を書き換えようとすると、変数 a が mutable でないためエラーが送出されます。変数 b のように mutable として宣言すると、値型のレコードのフィールド変数を書き換えることができます。

複素数

近年、多くのプログラミング言語で「複素数 (complex number)」がサポートされるようになりました。たとえば、Ｃ言語では 1999 年に発行された規格 C99 で複素数型が導入されました。Go 言語や Python でも複素数型をサポートしていますし、複素数用の標準ライブラリを用意している言語 (C++, C#, F#, Ruby, Haskell など) も多くあります。

他の言語では、FORTRAN や Common Lisp が昔から複素数型をサポートしています。Common Lisp の場合、基本的な数学関数でも複素数を適用できるのであれば、引数に複素数を渡して計算することができます。.NET (C#, F#) の場合、ライブラリ System.Numerics に複素数を表す構造体 Complex が用意されているので、簡単に複素数を扱うことができます。今回は F# の複素数についてまとめてみました。

●.NET の数

.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 に用意されている関数を使います。

• bool IsInfinity (d: float)
  引数 d が正または負の無限大のとき真を返す
• bool IsPositiveInfinity (d: float);
  引数 d が正の無限大のとき真を返す
• bool IsNegativeInfinity (d: float);
  引数 d が負の無限大のとき真を返す

> 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;;
[<Struct>]
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 には次の定数が定義されています。

• ImaginaryOne, 0 + 1i
• One, 1 + 0i
• Zero, 0 + 0i
• Infinity, ∞ + ∞ i
• NaN, NaN + NaN i

> 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

参考文献 1 によると、次のように場合分けすることで、\(x^2\) や \(y^2\) で生じ得る上位桁あふれを回避することができるそうです。

今回は 参考文献 1 のプログラムを F# に移植しました。

それでは実際に試してみましょう。

> #load "complex.fsx";;
... 略 ...

> open Complex;;
> Complex(1, 1);;
val it: Complex = complex(1, 1) {Imag = 1.0;
                                 Real = 1.0;}

> Complex(1, 1).Conjugate();;
val it: Complex = complex(1, -1) {Imag = -1.0;
                                  Real = 1.0;}

> Complex(1, 1).Conjugate().Conjugate();;
val it: Complex = complex(1, 1) {Imag = 1.0;
                                 Real = 1.0;}

> 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)

除算の場合、絶対値の計算と同様にオーバーフローの対策が必要になります。今回は 参考文献 1 のプログラムを F# に移植しました。

それでは実際に試してみましょう。

> let a = Complex(1, 2);;
val a: Complex = complex(1, 2)

> let b = Complex(3, 4);;
val b: Complex = complex(3, 4)

> a + b;;
val it: Complex = complex(4, 6) {Imag = 6.0;
                                 Real = 4.0;}

> a - b;;
val it: Complex = complex(-2, -2) {Imag = -2.0;
                                   Real = -2.0;}

> a * b;;
val it: Complex = complex(-5, 10) {Imag = 10.0;
                                   Real = -5.0;}

> a / b;;
val it: Complex = complex(0.44, 0.08) {Imag = 0.08;
                                       Real = 0.44;}
> let one = Complex(1, 0);;
val one: Complex = complex(1, 0)

> one / Complex(1e300, 1e300);;
val it: Complex = complex(5e-301, -5e-301) {Imag = -5e-301;
                                            Real = 5e-301;}

> one / Complex(1e301, 1e300);;
val it: Complex =
  complex(9.90099e-302, -9.90099e-303) {Imag = -9.900990099e-303;
                                        Real = 9.900990099e-302;}

> one / Complex(1e300, 1e301);;
val it: Complex =
  complex(9.90099e-303, -9.90099e-302) {Imag = -9.900990099e-302;
                                        Real = 9.900990099e-303;}

●参考文献, URL

1. 奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
2. IEEE 754 -- Wikipedia
3. IEEE 754における負のゼロ - Wikipedia
4. NaN - Wikipedia

●プログラムリスト

//
// 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)