Linux Programming / お気楽 Rust プログラミング超入門

簡単なプログラム

●多倍長整数

Rust で多倍長整数の演算を行う場合、クレート num を使うのがお手軽 (他のクレートもある)
num のインストールは拙作のページ「外部クレートの使い方」を参照
モジュール bigint にある構造体 BigInt が多倍長整数を表す
使うときは use num::bigint::BigInt; とする
Copy トレイトは実装されていないので、借用チェッカーには注意すること
算術演算子は多重定義によりそのまま使用できる
プリミティブ型 (i32 など) との演算も可能
BigInt は参照でも演算可能

●BigInt の生成

配列から生成する場合

fn new(sign: Sign, digits: Vec<u32>) -> BigInt
fn from_slice(sign: Sign, slice: &[u32]) -> BigInt
fn from_bytes_be(sign: Sign, bytes: &[u8]) -> BigInt
fn from_bytes_le(sign: Sign, bytes: &[u8]) -> BigInt
fn parse_bytes(buf: &[u8], radix: u32) -> Option<BigInt>
fn from_radix_be(sign: Sign, buf: &[u8], radix: u32) -> Option<BigInt>
fn from_radix_le(sign: Sign, buf: &[u8], radix: u32) -> Option<BigInt>

上記の関数は配列 [u32], [u8] を BigInt に変換する
引数に基数 radix の指定が無い関数は、配列の要素の型が基数になる
Sign は符号を表す列挙型

enum Sign { Minus, NoSign, Plus }

使用するときは use num::bigint::Sign とする
_be はビッグエンディアン、_le はリトルエンディアンを表す
引数に基数 radix を指定する関数は Option<BigInt> を返す
radix の範囲は 2 .. 256
parse_byte の引数 buf の指定にはバイトリテラル b"..." を使うと便利
文字列のかわりに [u8] を生成する

文字列から BigInt を生成する場合、num::Num トレイトの関数 from_str_radix を使うと便利

fn from_str_radix(s: &str, radix: u32) -> Result<BigInt, Self::FromStrRadixErrr>

引数 radix の範囲は 2 ..= 36
整数値を BigInt に変換する場合は関数 from を使うと便利

fn from(n: T) -> Option<BigInt>

T は i32, u32 などのプリミティブ型
num::FromPrimitive トレイトにはプリミティブを BigInt に変換する関数 from_XXX (XXX はデータ型) もある

fn from_i32(n: i32) -> Option<Self>
fn from_u32(n: u32) -> Option<Self>
fn from_i64(n: i64) -> Option<Self>
fn from_u64(n: u64) -> Option<Self>
... その他いろいろ ...

num::bigint::ToBigInt トレイトにもプリミティブを BigInt に変換するメソッドがある

fn to_bigint(&self) -> Option<BigInt>

●BigInt の変換

配列に変換する場合

fn to_bytes_be(&self) -> (Sign, Vec<u8>)
fn to_bytes_le(&self) -> (Sign, Vec<u8>)
fn to_u32_digits(&self) -> (Sign, Vec<u32>)
fn to_u64_digits(&self) -> (Sign, Vec<u64>)
fn to_radix_be(&self, radix: u32) -> (Sign, Vec<u8>)
fn to_radix_le(&self, radix: u32) -> (Sign, Vec<u8>)

引数 radix の範囲は 2 .. 256
文字列に変換する場合

fn to_str_radix(&self, radix: u32) -> String

引数 radix の範囲は 2 ..= 36
BigInt をプリミティブに変換するにはトレイト num::ToPrimitive あるメソッド to_XXX (XXX はデータ型) を使う

fn to_i32(&self) -> Option<i32>
fn to_u32(&self) -> Option<u32>
fn to_i64(&self) -> Option<i64>
fn to_u64(&self) -> Option<u64>
... その他いろいろ ...

返り値は Option<プリミティブ型>

●その他のトレイト

num::One

fn one() -> Self       // 1 を返す
fn set_one(&mut self)  // 1 をセット
fn is_one() -> bool    // 1 ならば true を返す

num::Zero

fn zero() -> Self       // 0 を返す
fn set_zero(&mut self)  // 0 をセット
fn is_zero() -> bool    // 10ならば true を返す

num::Integer

fn gcd(&self, other: &Self) -> Self   // 最大公約数
fn lcm(&self, other: &Self) -> Self   // 最小公倍数
fn is_even(&self) -> bool             // 偶数か？
fn is_odd(&self) -> bool              // 奇数か？
fn extended_gcd(&self, other: &Self) -> ExtendedGcd<Self>  // 拡張ユークリッドの互除法

num::Signed

fn abs(&self) -> Self          // 絶対値
fn is_positive(&self) -> bool  // 正か？
fn is_negative(&self) -> bool  // 負か？

●簡単な使用例

リスト : 多倍長整数の使用例 (samplenum/src/main.rs)

use num::bigint::{BigInt, Sign, ToBigInt};
use num::Num;
use num::One;
use num::Zero;

// 階乗
fn fact(n: i32) -> BigInt {
    if n == 0 {
        BigInt::one()
    } else {
        n * fact(n - 1)
    }
}

// フィボナッチ数
fn fibo(n: i32) -> BigInt {
    fn iter(n: i32, a: BigInt, b: BigInt) -> BigInt {
        if n == 0 {
            a
        } else {
            let c = a + &b;
            iter(n - 1, b, c)
        }
    }
    iter(n, BigInt::zero(), BigInt::one())
}

// 繰り返し版
fn fiboi(n: i32) -> BigInt {
    let mut a = BigInt::zero();
    let mut b = BigInt::one();
    let mut i = 0;
    while i < n {
        let c = a + &b;
        a = b;
        b = c;
        i += 1;
    }
    a
}

// 組み合わせの数
fn combination(n: i32, r: i32) -> BigInt {
    if n == 0 || r == 0 {
        BigInt::one()
    } else {
        combination(n, r - 1) * (n - r + 1) / r
    }
}

// カタラン数
fn catalan(n: i32) -> BigInt {
    combination(2 * n, n) / (n + 1)
}

// モンモール数
fn montmort(n: i32) -> BigInt {
    let mut a = BigInt::zero();
    let mut b = BigInt::one();
    let mut i = 1;
    while i < n {
        let c = (a + &b) * (i + 1);
        a = b;
        b = c;
        i += 1;
    }
    a
}

// 分割数 (二重再帰なのでとても遅い)
fn partition_number(n: i32) -> BigInt {
    fn part_num(n: i32, k: i32) -> BigInt {
        if n == 1 || k == 1 {
            BigInt::one()
        } else if n < 0 || k < 1 {
            BigInt::zero()
        } else {
            part_num(n - k, k) + part_num(n, k - 1)
        }
    }
    part_num(n, n)
}

// 動的計画法による高速化
fn partition_number_fast(n: usize) -> BigInt {
    let mut table = vec![BigInt::one(); n + 1];
    for k in 2 .. n + 1 {
        for m in k .. n + 1 {
            let a = &table[m] + &table[m - k];
            table[m] = a;
        }
    }
    table[n].clone()
}

fn main() {
    let a: Vec<u32> = vec![0,0,0,1];
    println!("Plus {:?} => {}", &a, BigInt::new(Sign::Plus, a.clone()));
    println!("Minus {:?} => {}", &a, BigInt::from_slice(Sign::Minus, &a));

    let b: Vec<u8> = vec![1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9];
    println!("b = {:?}", b);
    println!("b => {}", BigInt::from_bytes_be(Sign::Plus, &b));
    println!("b => {}", BigInt::from_bytes_le(Sign::Plus, &b));
    println!("b => {:?}", BigInt::from_radix_be(Sign::Plus, &b, 10));
    println!("b => {:?}", BigInt::from_radix_le(Sign::Plus, &b, 10));

    println!("1234567890123456789 => {:?}", BigInt::parse_bytes(b"1234567890123456789", 10));

    println!("1234567890123456789 => {:?}", BigInt::from_str_radix("1234567890123456789", 10));
    println!("0xffffffffu32 => {:?}", BigInt::from(0xffffffffu32));
    println!("0xffffffffffffffffu64 => {:?}", BigInt::from(0xffffffffffffffffu64));
    println!("-123456789 => {:?}", BigInt::from(-123456789));
    println!("123456789 => {:?}", 123456789.to_bigint());
    println!("0x7fffffffffffffffu64 => {:?}", 0x7fffffffffffffffu64.to_bigint());

    println!("----- fact -----");
    for n in 20 .. 31 {
        println!("fact({}) = {}", n, fact(n));
    }
    println!("----- fibo -----");
    for n in 50 .. 61 {
        println!("fibo({}) = {}", n, fibo(n));
        println!("fiboi({}) = {}", n, fiboi(n));
    }
    println!("----- combination -----");
    for n in 30 .. 51 {
        println!("({}, {}) = {}", n * 2, n, combination(n * 2, n));
    }
    println!("----- catalan -----");
    println!("catalan(50) = {}", catalan(50));
    println!("catalan(100) = {}", catalan(100));
    println!("----- montmort -----");
    for n in 1 .. 21 {
        println!("montmort({}) = {}", n, montmort(n));
    }
    println!("montmort(50) = {}", montmort(50));
    println!("----- partition_number -----");
    for n in 0 .. 11 {
        println!("partition_number({}) = {}", n, partition_number(n));
    }
    println!("partition_number(40) = {}", partition_number(40));
    println!("partition_number(60) = {}", partition_number(60));
    println!("partition_number_fast(80) = {}", partition_number_fast(80));
    println!("partition_number_fast(100) = {}", partition_number_fast(100));
    println!("partition_number_fast(200) = {}", partition_number_fast(200));
}

●実行例

$ cargo run
   Compiling samplenum v0.1.0 (/home/mhiroi/rust/samplenum)
    Finished dev [unoptimized + debuginfo] target(s) in 3.62s
     Running `target/debug/samplenum`
Plus [0, 0, 0, 1] => 79228162514264337593543950336
Minus [0, 0, 0, 1] => -79228162514264337593543950336
b = [1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
b => 22475995960490625424076925365263601926342665
b => 201405995052208963761489336373750075956068865
b => Some(1234567890123456789)
b => Some(9876543210987654321)
1234567890123456789 => Some(1234567890123456789)
1234567890123456789 => Ok(1234567890123456789)
0xffffffffu32 => 4294967295
0xffffffffffffffffu64 => 18446744073709551615
-123456789 => -123456789
123456789 => Some(123456789)
0x7fffffffffffffffu64 => Some(9223372036854775807)
----- fact -----
fact(20) = 2432902008176640000
fact(21) = 51090942171709440000
fact(22) = 1124000727777607680000
fact(23) = 25852016738884976640000
fact(24) = 620448401733239439360000
fact(25) = 15511210043330985984000000
fact(26) = 403291461126605635584000000
fact(27) = 10888869450418352160768000000
fact(28) = 304888344611713860501504000000
fact(29) = 8841761993739701954543616000000
fact(30) = 265252859812191058636308480000000
----- fibo -----
fibo(50) = 12586269025
fiboi(50) = 12586269025
fibo(51) = 20365011074
fiboi(51) = 20365011074
fibo(52) = 32951280099
fiboi(52) = 32951280099
fibo(53) = 53316291173
fiboi(53) = 53316291173
fibo(54) = 86267571272
fiboi(54) = 86267571272
fibo(55) = 139583862445
fiboi(55) = 139583862445
fibo(56) = 225851433717
fiboi(56) = 225851433717
fibo(57) = 365435296162
fiboi(57) = 365435296162
fibo(58) = 591286729879
fiboi(58) = 591286729879
fibo(59) = 956722026041
fiboi(59) = 956722026041
fibo(60) = 1548008755920
fiboi(60) = 1548008755920
----- combination -----
(60, 30) = 118264581564861424
(62, 31) = 465428353255261088
(64, 32) = 1832624140942590534
(66, 33) = 7219428434016265740
(68, 34) = 28453041475240576740
(70, 35) = 112186277816662845432
(72, 36) = 442512540276836779204
(74, 37) = 1746130564335626209832
(76, 38) = 6892620648693261354600
(78, 39) = 27217014869199032015600
(80, 40) = 107507208733336176461620
(82, 41) = 424784580848791721628840
(84, 42) = 1678910486211891090247320
(86, 43) = 6637553085023755473070800
(88, 44) = 26248505381684851188961800
(90, 45) = 103827421287553411369671120
(92, 46) = 410795449442059149332177040
(94, 47) = 1625701140345170250548615520
(96, 48) = 6435067013866298908421603100
(98, 49) = 25477612258980856902730428600
(100, 50) = 100891344545564193334812497256
----- catalan -----
catalan(50) = 1978261657756160653623774456
catalan(100) = 896519947090131496687170070074100632420837521538745909320
----- montmort -----
montmort(1) = 0
montmort(2) = 1
montmort(3) = 2
montmort(4) = 9
montmort(5) = 44
montmort(6) = 265
montmort(7) = 1854
montmort(8) = 14833
montmort(9) = 133496
montmort(10) = 1334961
montmort(11) = 14684570
montmort(12) = 176214841
montmort(13) = 2290792932
montmort(14) = 32071101049
montmort(15) = 481066515734
montmort(16) = 7697064251745
montmort(17) = 130850092279664
montmort(18) = 2355301661033953
montmort(19) = 44750731559645106
montmort(20) = 895014631192902121
montmort(50) = 11188719610782480504630258070757734324011354208865721592720336801
----- partition_number -----
partition_number(0) = 0
partition_number(1) = 1
partition_number(2) = 2
partition_number(3) = 3
partition_number(4) = 5
partition_number(5) = 7
partition_number(6) = 11
partition_number(7) = 15
partition_number(8) = 22
partition_number(9) = 30
partition_number(10) = 42
partition_number(40) = 37338
partition_number(60) = 966467
partition_number_fast(80) = 15796476
partition_number_fast(100) = 190569292
partition_number_fast(200) = 3972999029388

●複素数

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

●Rust の数

Rust の数はプリミティブな型として、整数 (int32, 1nt64 など) と実数 (f32, f64) の 2 種類があります。外部クレート num をインストールすると、多倍長整数、有理数 (分数)、複素数を使用することができます。

実数は浮動小数点数 (floating point number) として表現されます。浮動小数点数には IEEE 754 という標準仕様があり、近代的なプログラミング言語のほとんどは、IEEE 754 に準拠した浮動小数点数をサポートしています。浮動小数点数はすべての小数を正確に表現することはできません。このため、実数は近似的な値になります。

IEEE 754 には通常の数値以外にも、負のゼロ (-0.0)、正負の無限大 (∞, -∞)、NaN (Not a Number, 非数) といった値が定義されています。これらの値は Rust でも取り扱うことができます。負のゼロは -0、正負の無限大は inf と -inf、NaN は NaN と表示されます。

●無限大

一般に、無限大は値のオーバーフロー、ゼロ除算 (数値 / 0.0)、数学関数の計算結果 (たとえば log(0.0)) などで発生します。また、標準ライブラリのモジュール f64 には、正負の無限大を表す定数 INFINITY, NEG_INFINITY が定義されています。

簡単な実行例を示します。

リスト : 無限大

fn main() {
    let a: f64 = 1e308;
    let b: f64 = f64::INFINITY;
    let c: f64 = f64::NEG_INFINITY;
    println!("a = {:e}", a);
    println!("b = {}", b);
    println!("c = {}", c);
    println!("a * 2 = {}", a * 2.0);
    println!("-a * 2 = {}", -a * 2.0);
    println!("1.0 / 0.0 = {}", 1.0 / 0.0);
    println!("-1.0 / 0.0 = {}", -1.0 / 0.0);
    println!("0.0_f64.ln() = {}", 0.0_f64.ln());

    println!("a.is_infinite() = {}", a.is_infinite());
    println!("b.is_infinite() = {}", b.is_infinite());
    println!("c.is_infinite() = {}", b.is_infinite());

    println!("b == b => {}", b == b);
    println!("b == c => {}", b == c);
    println!("b > 0.0 => {}", b > 0.0);
    println!("b < 0.0 => {}", b < 0.0);
    println!("c > 0.0 => {}", c > 0.0);
    println!("c < 0.0 => {}", c < 0.0);
    println!("b > c => {}", b > c);
    println!("b < c => {}", b < b);

    println!("b + 10.0 = {}", b + 10.0);
    println!("b - 10.0 = {}", b - 10.0);
    println!("b * 10.0 = {}", b * 10.0);
    println!("b / 10.0 = {}", b / 10.0);

    println!("b + b = {}", b + b);
    println!("b * b = {}", b * b);
    println!("b - b = {}", b - b);
    println!("b / b = {}", b / b);
    println!("b + c = {}", b + c);
    println!("b * 0.0 = {}", b * 0.0);
}

a = 1e308
b = inf
c = -inf
a * 2 = inf
-a * 2 = -inf
1.0 / 0.0 = inf
-1.0 / 0.0 = -inf
0.0_f64.ln() = -inf
a.is_infinite() = false
b.is_infinite() = true
c.is_infinite() = true
b == b => true
b == c => false
b > 0.0 => true
b < 0.0 => false
c > 0.0 => false
c < 0.0 => true
b > c => true
b < c => false
b + 10.0 = inf
b - 10.0 = inf
b * 10.0 = inf
b / 10.0 = inf
b + b = inf
b * b = inf
b - b = NaN
b / b = NaN
b + c = NaN
b * 0.0 = NaN

Rust の場合、無限大の判定はモジュール f64 のメソッド is_infinite() を使います。無限大は他の数値と比較したり演算することもできますが、結果が NaN になることもあります。

●負のゼロ

負のゼロ (-0.0) は、計算結果が負の極めて小さな値でアンダーフローになったとき発生します。また、正の値を負の無限大で除算する、負の値を正の無限大で除算する、負の値と 0.0 を乗算しても -0.0 が得られます。

簡単な実行例を示します。

リスト : 負のゼロ

fn main() {
    let a: f64 = -1e-323;
    let b: f64 = f64::INFINITY;

    println!("a = {:e}", a);
    println!("b = {}", b);
    println!("a / 4.0 = {}", a / 4.0);
    println!("-1.0 / b = {}", -1.0 / b);
    println!("1.0 / -b = {}", 1.0 / -b);
    println!("-1.0 * 0.0 = {}", -1.0 * 0.0);

    println!("0.0 == -0.0 => {}", 0.0 == -0.0);
    println!("0.0 != -0.0 => {}", 0.0 != -0.0);
    println!("0.0_f64.is_sign_negative() = {}", 0.0_f64.is_sign_negative());
    println!("(-0.0_f64).is_sign_negative() = {}", (-0.0_f64).is_sign_negative());

    println!("0.0_f64.sqrt() = {}", 0.0_f64.sqrt());
    println!("-0.0_f64.sqrt() = {}", -0.0_f64.sqrt());
    println!("0.0_f64.atan2(-1.0) = {}", 0.0_f64.atan2(-1.0));
    println!("-0.0_f64.atan2(-1.0) = {}", -0.0_f64.atan2(-1.0));
}

a = -1e-323
b = inf
a / 4.0 = -0
-1.0 / b = -0
1.0 / -b = -0
-1.0 * 0.0 = -0
0.0 == -0.0 => true
0.0 != -0.0 => false
0.0_f64.is_sign_negative() = false
(-0.0_f64).is_sign_negative() = true
0.0_f64.sqrt() = 0
-0.0_f64.sqrt() = -0
0.0_f64.atan2(-1.0) = 3.141592653589793
-0.0_f64.atan2(-1.0) = -3.141592653589793

標準では、演算子 (Rust では ==) による 0.0 と -0.0 の比較は等しいと判定されます。Rust の場合、-0.0 はモジュール f64 のメソッド is_sign_negative で判別することができます。また、任意の正の数を -0.0 で除算すると負の無限大になるので、それを使って判別することもできます。なお、-0.0 は数学関数 (モジュール f64 のメソッド sqrt や atan2 など) や複素数の演算処理などで使われます。

●非数

NaN は数ではないことを表す特別な値 (非数) です。一般的には 0.0 / 0.0 といった不正な演算を行うと、その結果は NaN になります。NaN はモジュール f64 に定数 NAN として定義されています。

リスト : 非数

fn main() {
    let a: f64 = f64::INFINITY;
    let nan: f64 = f64::NAN;
    println!("a = {}", a);
    println!("nan = {}", nan);
    println!("nan / nan = {}", nan / nan);
    println!("a / a = {}", a / a);
    println!("a.is_nan() = {}", a.is_nan());
    println!("nan.is_nan() = {}", nan.is_nan());
    println!("a.is_finite() = {}", a.is_finite());
    println!("nan.is_finite() = {}", nan.is_finite());
    println!("1.0_f64.is_finite() = {}", 1.0_f64.is_finite());
}

a = inf
nan = NaN
nan / nan = NaN
a / a = NaN
a.is_nan() = false
nan.is_nan() = true
a.is_finite() = false
nan.is_finite() = false
1.0_f64.is_finite() = true

Rust の場合、NaN はモジュール f64 のメソッド is_nan() で判別することができます。メソッド is_finite() は引数が無限大でも NaN でもなければ真を返します。

●Rust (クレート num) の複素数

数学では複素数 \(z\) を \(x + yi\) と表記します。x を実部、y を虚部、i を虚数単位といいます。虚数単位は 2 乗すると -1 になる数です。実部と虚部の 2 つの数値を格納するデータ構造を用意すれば、プログラミング言語でも複素数を表すことができます。Rust (外部クレート num) では複素数を構造体 Complex で表します。

struct Complex<T> { re: T, im: T }
type Complex32 = Complex<f32>;
type Complex64 = Complex<f64>;

Compelx はジェネリクスですが、本ページでは Complex64 を使って説明することにします。フィールド re が実部を、フィールド im が虚部を表します。使用するときは use num::complex::Complex; としてください。(use num::Complex; でも大丈夫なようです。)

複素数はコンストラクタ Complex::new(x, y) で生成します。

fn new(re: T, im: T) -> Complex<T>

複素数の虚部の符号を反転することを「複素共役」といいます。複素共役はメソッド conj() で求めることができます。

fn conj(&self) -> Complex<T>

複素数は極形式 \(z = |z|(\cos \theta + i \sin \theta)\) で表すことができます。このとき、\(|z|\) を絶対値、\(\theta\) を偏角といいます。絶対値はメソッド norm() で求めることができます。偏角は複素平面において正の実軸とベクトル (x, y) との角度を表します。偏角はメソッド arg() で求めることができます。arg() の返り値 \(\theta\) は \(-\pi \leq \theta \leq \pi \) (\(\pi\) : 円周率) です。

fn norm(self) -> T
fn arg(self) -> T

このほかに、複素数の絶対値と偏角をタプルに格納して返すメソッド to_polar()、絶対値と偏角から複素数を生成する関数 from_polar(絶対値, 偏角) が用意されています。

fn to_polar(self) -> (T, T)
fn from_polar(r: T, theta: T) -> Complex<T>

簡単な例を示しましょう。

リスト : クレート num の複素数

use num::complex::Complex;

fn main() {
    let a = Complex::new(1.234, 5.6789);
    println!("a = {}", a);
    println!("a.re = {}", a.re);
    println!("a.im = {}", a.im);
    println!("a.conj() = {}", a.conj());
    println!("i = {}", Complex::<f64>::i());
    println!("-i = {}", -Complex::<f64>::i());
    let b = Complex::new(1.0, 1.0);
    println!("b = {}", b);
    println!("b.norm() = {}", b.norm());
    println!("b.arg() = {}", b.arg());
    println!("b.to_polar() = {:?}", b.to_polar());
    println!("from_polar(b.norm(), b.arg()) = {}", Complex::from_polar(b.norm(), b.arg()));
    for (x, y) in [(1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0),
                   (1.0, -1.0), (0.0, -1.0), (-1.0, -1.0), (-1.0, -0.0)] {
        println!("Complex::new({}, {}).arg() = {}", x, y, Complex::new(x, y).arg());
    }
}

a = 1.234+5.6789i
a.re = 1.234
a.im = 5.6789
a.conj() = 1.234-5.6789i
i = 0+1i
-i = -0-1i
b = 1+1i
b.norm() = 1.4142135623730951
b.arg() = 0.7853981633974483
b.to_polar() = (1.4142135623730951, 0.7853981633974483)
from_polar(b.norm(), b.arg()) = 1.0000000000000002+1i
Complex::new(1, 0).arg() = 0
Complex::new(1, 1).arg() = 0.7853981633974483
Complex::new(0, 1).arg() = 1.5707963267948966
Complex::new(-1, 1).arg() = 2.356194490192345
Complex::new(-1, 0).arg() = 3.141592653589793
Complex::new(1, -1).arg() = -0.7853981633974483
Complex::new(0, -1).arg() = -1.5707963267948966
Complex::new(-1, -1).arg() = -2.356194490192345
Complex::new(-1, -0).arg() = -3.141592653589793

関数 i() は虚数単位 i を返します。Rust (クレート num) の場合、-1.0+0.0i の偏角 θ は pi になり、-1.0-0.0i の偏角は -pi になります。ゼロと負のゼロを区別しないプログラミング言語では、偏角 θ の範囲を -pi < θ <= pi に制限して、-1.0 + 0.0j (== -1.0 - 0.0j) の偏角を pi とします。

●複素数の四則演算

複素数の四則演算は次のようになります。

\( \begin{array}{l} (a + bi) + (c + di) = (a + c) + (b + d)i \\ (a + bi) - (c + di) = (a - c) + (b - d)i \\ (a + bi) \times (c + di) = (ac - bd) + (bc + ad)i \\ (a + bi) \div (c + di) = \dfrac{ac + bd + (bc - ad)i}{c^2 + d^2} \end{array} \)

これらの計算は演算子の多重定義により +, -, * , / で行うことができます。他のプリミティブ型 (f64 など) との演算も可能です。複素数の場合、大小の比較演算は使えませんが、等値の判定は演算子 == や != で行うことができます。簡単な例を示しましょう。

リスト : 複素数の四則演算 (main() に追加)

fn main() {
    //
    // ・・・省略・・・
    //
    let b = Complex::new(1.0, 2.0);
    let c = Complex::new(3.0, 4.0);
    println!("b = {}", b);
    println!("c = {}", c);
    println!("b + c = {}", b + c);
    println!("b - c = {}", b - c);
    println!("b * c = {}", b * c);
    println!("b / c = {}", b / c);
    println!("b == b => {}", b == b);
    println!("b == c => {}", b == c);
    println!("b != c => {}", b != c);
    let d = Complex::new(1e300, 1e300);
    println!("d = {:?}", d);
    println!("d.norm() = {:e}", d.norm());
    let e = d * d;
    println!("e = d * d = {}", e);
    println!("e.is_infinite() = {}", e.is_infinite());
    println!("e.is_nan() = {}", e.is_nan());
    println!("1.0 / d = {}", 1.0 / d);
}

... 略 ...
b = 1+2i
c = 3+4i
b + c = 4+6i
b - c = -2-2i
b * c = -5+10i
b / c = 0.44+0.08i
b == b => true
b == c => false
b != c => true
d = Complex { re: 1e300, im: 1e300 }
d.norm() = 1.4142135623730952e300
e = d * d = NaN+infi
e.is_infinite() = false
e.is_nan() = true
1.0 / d = 0+0i

Rust (クレート num) の場合、Complex を書式 {} で表示すると、実部と虚部の実数は小数点形式 (Ｃ言語の書式 %f) で表示されます。{:?} で表示すると、実部と虚部の実数は指数形式 (Ｃ言語の書式 %e) で表示されます。f64 は書式 {} で表示すると小数点形式、{:e} で表示すると指数形式で表示されます。Ｃ言語には %f と %e のどちらか短いほうで表示する書式 %g がありますが、Rust にはそのような書式はないようです。ご注意くださいませ。

また、実部と虚部の値は無限大や NaN になることもあります。Complex にも無限大や NaN を判定するメソッド is_infinite(), is_nan() が用意されています。

ところで、複素数 \(x + yi\) の絶対値は \(\sqrt {x^2 + y^2}\) で求めることができますが、絶対値を求める場合、定義をそのままプログラムすると二乗の計算でオーバーフローすることがあります。たとえば、1e300 + 1e300i の絶対値を求めてみましょう。このとき、1e300 の二乗でオーバーフローします。Rust の場合は無限大になります。

参考文献『Ｃ言語による最新アルゴリズム事典』によると、次のように場合分けすることで、\(x^2\) や \(y^2\) で生じ得る上位桁あふれを回避することができるそうです。

\( \sqrt{x^2 + y^2} = \begin{cases} |x| \sqrt{1 + \left(\dfrac{y}{x}\right)^2} \quad & \mathrm{if} \ |x| \geq |y| \\ |y| \sqrt{1 + \left(\dfrac{x}{y}\right)^2} \quad & \mathrm{other} \end{cases} \)

同様に、割り算でもオーバーフローの対策が必要になります。クレート num のメソッド norm() はオーバーフローの対策が行われているようですが、割り算は対策されていないようです。1.0 / d を計算すると 0 + 0i になってしまいました。正しい値は 5e-301 - 5e-301 i になります。num の Complex で大きな値を計算するときには注意してください。

●参考文献, URL

奥村晴彦,『Ｃ言語による最新アルゴリズム事典』, 技術評論社, 1991
プログラミング言語における複素数の実装格付けチェック, (雑記帳 | 人生やっていき)

●プログラムリスト

//
// クレート num の使用例 (samplenum/src/main.rs)
//
// Copyright (C) 2022-2024 Makoto Hiroi
//
use num::bigint::{BigInt, Sign, ToBigInt};
use num::Num;
use num::One;
use num::Zero;
use num::complex::Complex;

// 階乗
fn fact(n: i32) -> BigInt {
    if n == 0 {
        BigInt::one()
    } else {
        n * fact(n - 1)
    }
}

// フィボナッチ数
fn fibo(n: i32) -> BigInt {
    fn iter(n: i32, a: BigInt, b: BigInt) -> BigInt {
        if n == 0 {
            a
        } else {
            let c = a + &b;
            iter(n - 1, b, c)
        }
    }
    iter(n, BigInt::zero(), BigInt::one())
}

// 繰り返し版
fn fiboi(n: i32) -> BigInt {
    let mut a = BigInt::zero();
    let mut b = BigInt::one();
    let mut i = 0;
    while i < n {
        let c = a + &b;
        a = b;
        b = c;
        i += 1;
    }
    a
}

// 組み合わせの数
fn combination(n: i32, r: i32) -> BigInt {
    if n == 0 || r == 0 {
        BigInt::one()
    } else {
        combination(n, r - 1) * (n - r + 1) / r
    }
}

// カタラン数
fn catalan(n: i32) -> BigInt {
    combination(2 * n, n) / (n + 1)
}

// モンモール数
fn montmort(n: i32) -> BigInt {
    let mut a = BigInt::zero();
    let mut b = BigInt::one();
    let mut i = 1;
    while i < n {
        let c = (a + &b) * (i + 1);
        a = b;
        b = c;
        i += 1;
    }
    a
}

// 分割数 (二重再帰なのでとても遅い)
fn partition_number(n: i32) -> BigInt {
    fn part_num(n: i32, k: i32) -> BigInt {
        if n == 1 || k == 1 {
            BigInt::one()
        } else if n < 0 || k < 1 {
            BigInt::zero()
        } else {
            part_num(n - k, k) + part_num(n, k - 1)
        }
    }
    part_num(n, n)
}

// 動的計画法による高速化
fn partition_number_fast(n: usize) -> BigInt {
    let mut table = vec![BigInt::one(); n + 1];
    for k in 2 .. n + 1 {
        for m in k .. n + 1 {
            let a = &table[m] + &table[m - k];
            table[m] = a;
        }
    }
    table[n].clone()
}

// 多倍長整数のテスト
fn test_bigint() {
    let a: Vec<u32> = vec![0,0,0,1];
    println!("Plus {:?} => {}", &a, BigInt::new(Sign::Plus, a.clone()));
    println!("Minus {:?} => {}", &a, BigInt::from_slice(Sign::Minus, &a));

    let b: Vec<u8> = vec![1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9];
    println!("b = {:?}", b);
    println!("b => {}", BigInt::from_bytes_be(Sign::Plus, &b));
    println!("b => {}", BigInt::from_bytes_le(Sign::Plus, &b));
    println!("b => {:?}", BigInt::from_radix_be(Sign::Plus, &b, 10));
    println!("b => {:?}", BigInt::from_radix_le(Sign::Plus, &b, 10));

    println!("1234567890123456789 => {:?}", BigInt::parse_bytes(b"1234567890123456789", 10));

    println!("1234567890123456789 => {:?}", BigInt::from_str_radix("1234567890123456789", 10));
    println!("0xffffffffu32 => {:?}", BigInt::from(0xffffffffu32));
    println!("0xffffffffffffffffu64 => {:?}", BigInt::from(0xffffffffffffffffu64));
    println!("-123456789 => {:?}", BigInt::from(-123456789));
    println!("123456789 => {:?}", 123456789.to_bigint());
    println!("0x7fffffffffffffffu64 => {:?}", 0x7fffffffffffffffu64.to_bigint());

    println!("----- fact -----");
    for n in 20 .. 31 {
        println!("fact({}) = {}", n, fact(n));
    }
    println!("----- fibo -----");
    for n in 50 .. 61 {
        println!("fibo({}) = {}", n, fibo(n));
        println!("fiboi({}) = {}", n, fiboi(n));
    }
    println!("----- combination -----");
    for n in 30 .. 51 {
        println!("({}, {}) = {}", n * 2, n, combination(n * 2, n));
    }
    println!("----- catalan -----");
    println!("catalan(50) = {}", catalan(50));
    println!("catalan(100) = {}", catalan(100));
    println!("----- montmort -----");
    for n in 1 .. 21 {
        println!("montmort({}) = {}", n, montmort(n));
    }
    println!("montmort(50) = {}", montmort(50));
    println!("----- partition_number -----");
    for n in 0 .. 11 {
        println!("partition_number({}) = {}", n, partition_number(n));
    }
    println!("partition_number(40) = {}", partition_number(40));
    println!("partition_number(60) = {}", partition_number(60));
    println!("partition_number_fast(80) = {}", partition_number_fast(80));
    println!("partition_number_fast(100) = {}", partition_number_fast(100));
    println!("partition_number_fast(200) = {}", partition_number_fast(200));
}

fn main() {
    test_bigint();
    // 複素数のテスト
    println!("----- Complex -----");
    let a = Complex::new(1.234, 5.6789);
    println!("a = {}", a);
    println!("a.re = {}", a.re);
    println!("a.im = {}", a.im);
    println!("a.conj() = {}", a.conj());
    println!("i = {}", Complex::<f64>::i());
    println!("-i = {}", -Complex::<f64>::i());
    let b = Complex::new(1.0, 1.0);
    println!("b = {}", b);
    println!("b.norm() = {}", b.norm());
    println!("b.arg() = {}", b.arg());
    println!("b.to_polar() = {:?}", b.to_polar());
    println!("from_polar(b.norm(), b.arg()) = {}", Complex::from_polar(b.norm(), b.arg()));
    for (x, y) in [(1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 0.0),
                   (1.0, -1.0), (0.0, -1.0), (-1.0, -1.0), (-1.0, -0.0)] {
        println!("Complex::new({}, {}).arg() = {}", x, y, Complex::new(x, y).arg());
    }
    let b = Complex::new(1.0, 2.0);
    let c = Complex::new(3.0, 4.0);
    println!("b = {}", b);
    println!("c = {}", c);
    println!("b + c = {}", b + c);
    println!("b - c = {}", b - c);
    println!("b * c = {}", b * c);
    println!("b / c = {}", b / c);
    println!("b == b => {}", b == b);
    println!("b == c => {}", b == c);
    println!("b != c => {}", b != c);
    let d = Complex::new(1e300, 1e300);
    println!("d = {:?}", d);
    println!("d.norm() = {:e}", d.norm());
    let e = d * d;
    println!("e = d * d = {}", e);
    println!("e.is_infinite() = {}", e.is_infinite());
    println!("e.is_nan() = {}", e.is_nan());
    println!("1.0 / d = {}", 1.0 / d);
    println!("1e300 * 1e300 = {}", 1e300 * 1e300);
}