Complex デモ — 複素数演算

sangi の Complex<T> は double から多倍長 Float まで統一的に扱える複素数クラスである。このページでは 3 つのデモで代表的なユースケースを紹介する。

Demo 1: 複素数の基本演算

Complex<double> で四則演算、絶対値 (abs)、偏角 (arg)、共役 (conj) を計算する。 using sangi::literals を宣言すると 3.0_i のようなユーザー定義リテラルで虚数を直感的に記述できる。

=== Complex<double> Basic Operations === a = (2+3i) b = (1-1i) a + b = (3+2i) a - b = (1+4i) a * b = (5+1i) a / b = (-0.5+2.5i) abs(a) = 3.60555 arg(a) = 0.982794 (rad) conj(a) = (2-3i)

Demo 1 ソースコード

// demo1_basic.cpp — Complex<double> の基本演算
#include <math/core/Complex.hpp>
#include <iostream>
#include <iomanip>
using namespace sangi;
using namespace sangi::literals;  // 3.0_i リテラルを有効化

int main() {
    std::cout << std::setprecision(6);
    std::cout << "=== Complex<double> Basic Operations ===\n";

    // ユーザー定義リテラルで複素数を生成
    auto a = 2.0 + 3.0_i;   // (2, 3)
    auto b = 1.0 - 1.0_i;   // (1, -1)

    std::cout << "a = " << a << "\n";
    std::cout << "b = " << b << "\n\n";

    // 四則演算
    std::cout << "a + b = " << (a + b) << "\n";
    std::cout << "a - b = " << (a - b) << "\n";
    std::cout << "a * b = " << (a * b) << "\n";
    std::cout << "a / b = " << (a / b) << "\n\n";

    // 絶対値・偏角・共役
    std::cout << "abs(a) = " << abs(a) << "\n";
    std::cout << "arg(a) = " << arg(a) << "  (rad)\n";
    std::cout << "conj(a) = " << conj(a) << "\n";

    return 0;
}

sangi::literals 名前空間を使うと 3.0_i が Complex<double>(0, 3.0) に変換される。通常の算術演算子との組み合わせで、数学的な記法に近いコードが書ける。

Demo 2: マンデルブロ集合の判定

複素数 $c$ がマンデルブロ集合に属するかどうかを、漸化式 $z_{n+1} = z_n^2 + c$ ($z_0 = 0$) の反復で判定する。 100 回の反復で $|z| > 2$ にならなければ、その点は集合に属すると見なす。格子上の数点を判定し、結果を表示する。

=== Mandelbrot Set Membership === c = (0+0i) : IN (|z| = 0.000000, iter = 100) c = (0.25+0i) : IN (|z| = 0.500000, iter = 100) c = (-1+0i) : IN (|z| = 0.000000, iter = 100) c = (1+0i) : OUT (|z| = 2.396970, iter = 3) c = (-0.5+0.5i) : IN (|z| = 0.627930, iter = 100) c = (0.3+0.5i) : OUT (|z| = 2.301270, iter = 10) c = (-1.5+0i) : OUT (|z| = 2.261458, iter = 24) c = (-0.75+0.1i) : IN (|z| = 1.324718, iter = 100)

Demo 2 ソースコード

// demo2_mandelbrot.cpp — マンデルブロ集合の判定
#include <math/core/Complex.hpp>
#include <iostream>
#include <iomanip>
#include <vector>
using namespace sangi;

bool isMandelbrot(Complex<double> c, int maxIter, int& iterOut,
                  double& absOut) {
    Complex<double> z(0.0, 0.0);
    for (int i = 0; i < maxIter; ++i) {
        z = z * z + c;
        if (abs(z) > 2.0) {
            iterOut = i + 1;
            absOut = abs(z);
            return false;  // 発散 → 集合に属さない
        }
    }
    iterOut = maxIter;
    absOut = abs(z);
    return true;  // 収束 → 集合に属する
}

int main() {
    std::cout << std::setprecision(6) << std::fixed;
    std::cout << "=== Mandelbrot Set Membership ===\n";

    std::vector<Complex<double>> points = {
        {0.0, 0.0}, {0.25, 0.0}, {-1.0, 0.0}, {1.0, 0.0},
        {-0.5, 0.5}, {0.3, 0.5}, {-1.5, 0.0}, {-0.75, 0.1}
    };

    constexpr int maxIter = 100;
    for (auto& c : points) {
        int iter;
        double absZ;
        bool inSet = isMandelbrot(c, maxIter, iter, absZ);
        std::cout << "c = " << std::setw(16) << std::left << c
                  << ": " << (inSet ? "IN " : "OUT")
                  << "  (|z| = " << absZ
                  << ", iter = " << iter << ")\n";
    }

    return 0;
}

マンデルブロ集合は「$z_{n+1} = z_n^2 + c$ が発散しない $c$ の集合」である。 $|z| > 2$ になった時点で発散が確定するため、閾値 2 で早期打ち切りが可能である。原点 $(0, 0)$ や $(-1, 0)$ は集合に属し、$(1, 0)$ はわずか 3 反復で発散する。

Demo 3: 複素数の数学関数

オイラーの公式 $e^{i\pi} = -1$ を Complex<double> で検証し、さらに複素数の sin, cos, exp, log, sqrt を計算する。最後に多倍長 Complex<Float> で $e^{i\pi} + 1 = 0$ を 100 桁精度で検証する。

=== Euler's Formula: exp(i*pi) === exp(i*pi) = (-1+1.22465e-16i) exp(i*pi) + 1 = (0+1.22465e-16i) (imag error is double rounding noise) === Complex Math Functions === z = (1+1i) sin(z) = (1.29846+0.634964i) cos(z) = (0.83373-0.988898i) exp(z) = (1.46869+2.28736i) log(z) = (0.346574+0.785398i) sqrt(z) = (1.09868+0.45509i) === High-Precision: Complex<Float> (100 digits) === exp(i*pi) + 1 = real: 0.0000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000000000 imag: 0.0000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000000000

Demo 3 ソースコード

// demo3_math_functions.cpp — 複素数の数学関数と高精度計算
#include <math/core/Complex.hpp>
#include <math/core/Float.hpp>
#include <iostream>
#include <iomanip>
#include <numbers>
using namespace sangi;
using namespace sangi::literals;

int main() {
    // --- オイラーの公式 ---
    std::cout << "=== Euler's Formula: exp(i*pi) ===\n";
    auto ipi = Complex<double>(0.0, std::numbers::pi);
    auto euler = exp(ipi);
    std::cout << "exp(i*pi) = " << euler << "\n";
    std::cout << "exp(i*pi) + 1 = " << (euler + 1.0) << "\n";
    std::cout << "  (虚部の誤差は double の丸め誤差)\n";

    // --- 複素数の数学関数 ---
    std::cout << "\n=== Complex Math Functions ===\n";
    auto z = 1.0 + 1.0_i;
    std::cout << std::setprecision(6);
    std::cout << "z = " << z << "\n";
    std::cout << "sin(z)  = " << sin(z) << "\n";
    std::cout << "cos(z)  = " << cos(z) << "\n";
    std::cout << "exp(z)  = " << exp(z) << "\n";
    std::cout << "log(z)  = " << log(z) << "\n";
    std::cout << "sqrt(z) = " << sqrt(z) << "\n";

    // --- 多倍長 Complex<Float> ---
    std::cout << "\n=== High-Precision: Complex<Float> (100 digits) ===\n";
    constexpr int digits = 100;
    Float::setDefaultPrecision(digits);

    Float pi = Float::pi(digits);
    Complex<Float> ipi_hp(Float(0), pi);
    Complex<Float> result = exp(ipi_hp) + Complex<Float>(Float(1));

    std::cout << "exp(i*pi) + 1 =\n";
    std::cout << "  real: " << result.real().toString(digits) << "\n";
    std::cout << "  imag: " << result.imag().toString(digits) << "\n";

    return 0;
}

Complex<double> ではオイラーの公式に $\approx 10^{-16}$ の丸め誤差が残るが、 Complex<Float> を使えば任意精度で $e^{i\pi} + 1 = 0$ を検証できる。 sangi の Complex<T> はテンプレートパラメータを切り替えるだけで精度を変更でき、 API は完全に共通である。

ソースコードと実行方法

example_complex.cpp (完全なソースコード)

// example_complex.cpp
// Complex<T> demo: basic operations, Mandelbrot membership, math functions,
// high-precision Complex<Float>

#include <math/core/Complex.hpp>
#include <math/core/mp/Float.hpp>
#include <iostream>
#include <iomanip>
#include <vector>
#include <numbers>

using namespace sangi;
using namespace sangi::literals;

// --- Demo 1: Basic arithmetic ---
void demo1_basic() {
    std::cout << std::setprecision(6);
    std::cout << "=== Complex<double> Basic Operations ===\n";

    auto a = 2.0 + 3.0_i;
    auto b = 1.0 - 1.0_i;

    std::cout << "a = " << a << "\n";
    std::cout << "b = " << b << "\n\n";

    std::cout << "a + b = " << (a + b) << "\n";
    std::cout << "a - b = " << (a - b) << "\n";
    std::cout << "a * b = " << (a * b) << "\n";
    std::cout << "a / b = " << (a / b) << "\n\n";

    std::cout << "abs(a) = " << abs(a) << "\n";
    std::cout << "arg(a) = " << arg(a) << "  (rad)\n";
    std::cout << "conj(a) = " << conj(a) << "\n";
}

// --- Demo 2: Mandelbrot set membership ---
bool isMandelbrot(Complex<double> c, int maxIter, int& iterOut,
                  double& absOut) {
    Complex<double> z(0.0, 0.0);
    for (int i = 0; i < maxIter; ++i) {
        z = z * z + c;
        if (abs(z) > 2.0) {
            iterOut = i + 1;
            absOut = abs(z);
            return false;
        }
    }
    iterOut = maxIter;
    absOut = abs(z);
    return true;
}

void demo2_mandelbrot() {
    std::cout << "\n=== Mandelbrot Set Membership ===\n";
    std::cout << std::setprecision(6) << std::fixed;

    std::vector<Complex<double>> points = {
        {0.0, 0.0}, {0.25, 0.0}, {-1.0, 0.0}, {1.0, 0.0},
        {-0.5, 0.5}, {0.3, 0.5}, {-1.5, 0.0}, {-0.75, 0.1}
    };

    constexpr int maxIter = 100;
    for (auto& c : points) {
        int iter;
        double absZ;
        bool inSet = isMandelbrot(c, maxIter, iter, absZ);
        std::cout << "c = " << std::setw(16) << std::left << c
                  << ": " << (inSet ? "IN " : "OUT")
                  << "  (|z| = " << absZ
                  << ", iter = " << iter << ")\n";
    }
}

// --- Demo 3: Math functions and high-precision computation ---
void demo3_math() {
    std::cout << "\n=== Euler's Formula: exp(i*pi) ===\n";
    std::cout << std::setprecision(6) << std::defaultfloat;
    auto ipi = Complex<double>(0.0, std::numbers::pi);
    auto euler = exp(ipi);
    std::cout << "exp(i*pi) = " << euler << "\n";
    std::cout << "exp(i*pi) + 1 = " << (euler + 1.0) << "\n";
    std::cout << "  (imag error is double rounding noise)\n";

    std::cout << "\n=== Complex Math Functions ===\n";
    auto z = 1.0 + 1.0_i;
    std::cout << "z = " << z << "\n";
    std::cout << "sin(z)  = " << sin(z) << "\n";
    std::cout << "cos(z)  = " << cos(z) << "\n";
    std::cout << "exp(z)  = " << exp(z) << "\n";
    std::cout << "log(z)  = " << log(z) << "\n";
    std::cout << "sqrt(z) = " << sqrt(z) << "\n";

    std::cout << "\n=== High-Precision: Complex<Float> (100 digits) ===\n";
    constexpr int digits = 100;
    Float::setDefaultPrecision(digits);

    Float pi = Float::pi(digits);
    Complex<Float> ipi_hp(Float(0), pi);
    Complex<Float> result = exp(ipi_hp) + Complex<Float>(Float(1));

    std::cout << "exp(i*pi) + 1 =\n";
    std::cout << "  real: " << result.real().toString(digits) << "\n";
    std::cout << "  imag: " << result.imag().toString(digits) << "\n";
}

int main() {
    demo1_basic();
    demo2_mandelbrot();
    demo3_math();
    return 0;
}

API の詳細は Complex API リファレンスを参照のこと。

ビルドと実行

cd sangi
mkdir build && cd build
cmake .. -G "Visual Studio 17 2022" -A x64
cmake --build . --config Release --target example-complex
examples\Release\example-complex.exe