// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// example_complex.cpp
// Complex<T> demo: basic operations, Mandelbrot membership, math functions,
// high-precision Complex<Float>
//
// Build:
//   cd build && cmake --build . --config Release --target example-complex

#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;
}