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

// Octonion.hpp
// 八元数 (Octonion / Cayley 数) テンプレートクラス
//
// o = e0 + e1·i + e2·j + e3·k + e4·l + e5·li + e6·lj + e7·lk
//
// Cayley-Dickson 構成: O = H + Hl  (H = 四元数)
// 乗法: (a, b)(c, d) = (ac - d*b, da + bc*)  [* = 共役]
//
// 性質:
//   - 非可換, 非結合 (代替的 alternative)
//   - ノルムは積性: |pq| = |p|·|q|  (合成代数)
//   - 逆元が存在: q^(-1) = conj(q)/|q|^2
//
// 用途: 弦理論, M理論, 例外型リー群 (G₂), Fano 平面

#pragma once

#include <iosfwd>
#include <sstream>
#include <string>
#include <array>
#include <cmath>
#include <concepts>
#include <limits>

namespace calx {

template<typename T>
concept OctonionScalar = requires(T a, T b) {
    { a + b } -> std::convertible_to<T>;
    { a - b } -> std::convertible_to<T>;
    { a * b } -> std::convertible_to<T>;
    { a / b } -> std::convertible_to<T>;
    { -a } -> std::convertible_to<T>;
    { T(0) };
    { T(1) };
};

// ================================================================
// Octonion<T> クラス
// ================================================================
// メモリレイアウト: e[0]..e[7]
// e[0] = 実部 (スカラー部)
// e[1]..e[7] = 虚部 (i, j, k, l, li, lj, lk)

template<OctonionScalar T>
class Octonion {
public:
    using value_type = T;

    std::array<T, 8> e;  // e[0]=実部, e[1..7]=虚部

    // ============================================================
    // コンストラクタ
    // ============================================================

    constexpr Octonion() noexcept(noexcept(T(0)))
        : e{T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0)} {}

    constexpr explicit Octonion(const T& scalar) noexcept
        : e{scalar, T(0), T(0), T(0), T(0), T(0), T(0), T(0)} {}

    constexpr Octonion(const T& e0, const T& e1, const T& e2, const T& e3,
                       const T& e4, const T& e5, const T& e6, const T& e7) noexcept
        : e{e0, e1, e2, e3, e4, e5, e6, e7} {}

    // 配列から構築
    constexpr explicit Octonion(const std::array<T, 8>& arr) noexcept : e(arr) {}

    // ============================================================
    // アクセサ
    // ============================================================

    [[nodiscard]] constexpr T  real()   const { return e[0]; }
    [[nodiscard]] constexpr T& real()         { return e[0]; }

    [[nodiscard]] constexpr const T& operator[](size_t i) const { return e[i]; }
    [[nodiscard]] constexpr T&       operator[](size_t i)       { return e[i]; }

    // ============================================================
    // 比較演算子
    // ============================================================

    [[nodiscard]] constexpr bool operator==(const Octonion& rhs) const {
        return e == rhs.e;
    }
    [[nodiscard]] constexpr bool operator!=(const Octonion& rhs) const {
        return e != rhs.e;
    }

    // ============================================================
    // 単項演算子
    // ============================================================

    [[nodiscard]] constexpr Octonion operator+() const { return *this; }
    [[nodiscard]] constexpr Octonion operator-() const {
        return Octonion(-e[0], -e[1], -e[2], -e[3],
                        -e[4], -e[5], -e[6], -e[7]);
    }

    // ============================================================
    // 加減算
    // ============================================================

    [[nodiscard]] constexpr Octonion operator+(const Octonion& rhs) const {
        return Octonion(e[0]+rhs.e[0], e[1]+rhs.e[1], e[2]+rhs.e[2], e[3]+rhs.e[3],
                        e[4]+rhs.e[4], e[5]+rhs.e[5], e[6]+rhs.e[6], e[7]+rhs.e[7]);
    }
    [[nodiscard]] constexpr Octonion operator-(const Octonion& rhs) const {
        return Octonion(e[0]-rhs.e[0], e[1]-rhs.e[1], e[2]-rhs.e[2], e[3]-rhs.e[3],
                        e[4]-rhs.e[4], e[5]-rhs.e[5], e[6]-rhs.e[6], e[7]-rhs.e[7]);
    }

    // ============================================================
    // 乗算 — Cayley-Dickson 構成
    // ============================================================
    // (a, b)(c, d) = (ac - d*b, da + bc*)
    // ここで a,b,c,d は四元数、* は四元数の共役。
    // 展開すると以下の 8×8 乗算テーブルになる。
    //
    // 乗算テーブル (Fano 平面に基づく):
    //   e_i · e_j の符号は下記配列で定義。
    //   乗算は O(1) に展開して直接計算。

    [[nodiscard]] constexpr Octonion operator*(const Octonion& q) const {
        const auto& a = e;
        const auto& b = q.e;
        return Octonion(
            a[0]*b[0] - a[1]*b[1] - a[2]*b[2] - a[3]*b[3]
          - a[4]*b[4] - a[5]*b[5] - a[6]*b[6] - a[7]*b[7],

            a[0]*b[1] + a[1]*b[0] + a[2]*b[3] - a[3]*b[2]
          + a[4]*b[5] - a[5]*b[4] - a[6]*b[7] + a[7]*b[6],

            a[0]*b[2] - a[1]*b[3] + a[2]*b[0] + a[3]*b[1]
          + a[4]*b[6] + a[5]*b[7] - a[6]*b[4] - a[7]*b[5],

            a[0]*b[3] + a[1]*b[2] - a[2]*b[1] + a[3]*b[0]
          + a[4]*b[7] - a[5]*b[6] + a[6]*b[5] - a[7]*b[4],

            a[0]*b[4] - a[1]*b[5] - a[2]*b[6] - a[3]*b[7]
          + a[4]*b[0] + a[5]*b[1] + a[6]*b[2] + a[7]*b[3],

            a[0]*b[5] + a[1]*b[4] - a[2]*b[7] + a[3]*b[6]
          - a[4]*b[1] + a[5]*b[0] - a[6]*b[3] + a[7]*b[2],

            a[0]*b[6] + a[1]*b[7] + a[2]*b[4] - a[3]*b[5]
          - a[4]*b[2] + a[5]*b[3] + a[6]*b[0] - a[7]*b[1],

            a[0]*b[7] - a[1]*b[6] + a[2]*b[5] + a[3]*b[4]
          - a[4]*b[3] - a[5]*b[2] + a[6]*b[1] + a[7]*b[0]
        );
    }

    // 除算: p / q = p * q^(-1)
    [[nodiscard]] constexpr Octonion operator/(const Octonion& rhs) const {
        return *this * rhs.inverse();
    }

    // スカラー乗除算
    [[nodiscard]] constexpr Octonion operator*(const T& s) const {
        return Octonion(e[0]*s, e[1]*s, e[2]*s, e[3]*s,
                        e[4]*s, e[5]*s, e[6]*s, e[7]*s);
    }
    [[nodiscard]] constexpr Octonion operator/(const T& s) const {
        return Octonion(e[0]/s, e[1]/s, e[2]/s, e[3]/s,
                        e[4]/s, e[5]/s, e[6]/s, e[7]/s);
    }
    friend constexpr Octonion operator*(const T& s, const Octonion& q) {
        return Octonion(s*q.e[0], s*q.e[1], s*q.e[2], s*q.e[3],
                        s*q.e[4], s*q.e[5], s*q.e[6], s*q.e[7]);
    }

    // ============================================================
    // 複合代入演算子
    // ============================================================

    constexpr Octonion& operator+=(const Octonion& rhs) {
        for (int i = 0; i < 8; ++i) e[i] += rhs.e[i];
        return *this;
    }
    constexpr Octonion& operator-=(const Octonion& rhs) {
        for (int i = 0; i < 8; ++i) e[i] -= rhs.e[i];
        return *this;
    }
    constexpr Octonion& operator*=(const Octonion& rhs) {
        *this = *this * rhs; return *this;
    }
    constexpr Octonion& operator/=(const Octonion& rhs) {
        *this = *this / rhs; return *this;
    }
    constexpr Octonion& operator*=(const T& s) {
        for (int i = 0; i < 8; ++i) e[i] *= s;
        return *this;
    }
    constexpr Octonion& operator/=(const T& s) {
        for (int i = 0; i < 8; ++i) e[i] /= s;
        return *this;
    }

    // ============================================================
    // 八元数固有の操作
    // ============================================================

    // 共役: conj(o) = e0 - e1·i - e2·j - ... - e7·lk
    [[nodiscard]] constexpr Octonion conj() const {
        return Octonion(e[0], -e[1], -e[2], -e[3],
                        -e[4], -e[5], -e[6], -e[7]);
    }

    // ノルムの二乗: |o|^2 = Σ e_i^2
    [[nodiscard]] constexpr T normSq() const {
        T s = T(0);
        for (int i = 0; i < 8; ++i) s += e[i] * e[i];
        return s;
    }

    // ノルム: |o|
    [[nodiscard]] T norm() const {
        using std::sqrt;
        return sqrt(normSq());
    }

    // 逆元: o^(-1) = conj(o) / |o|^2
    [[nodiscard]] constexpr Octonion inverse() const {
        T ns = normSq();
        if (ns == T(0)) {
            T nan = std::numeric_limits<T>::quiet_NaN();
            return Octonion(nan, nan, nan, nan, nan, nan, nan, nan);
        }
        T d = T(1) / ns;
        return Octonion(e[0]*d, -e[1]*d, -e[2]*d, -e[3]*d,
                        -e[4]*d, -e[5]*d, -e[6]*d, -e[7]*d);
    }

    // 単位八元数に正規化
    [[nodiscard]] Octonion normalized() const {
        T n = norm();
        if (n == T(0)) {
            T nan = std::numeric_limits<T>::quiet_NaN();
            return Octonion(nan, nan, nan, nan, nan, nan, nan, nan);
        }
        return *this / n;
    }

    // 虚部のノルム: |v| = sqrt(e1^2 + ... + e7^2)
    [[nodiscard]] T imagNorm() const {
        using std::sqrt;
        T s = T(0);
        for (int i = 1; i < 8; ++i) s += e[i] * e[i];
        return sqrt(s);
    }

    // ============================================================
    // 指数・対数
    // ============================================================

    // exp(o) = e^(e0) * (cos|v| + v̂·sin|v|)
    // v = 虚部, v̂ = v/|v|
    [[nodiscard]] Octonion exp() const {
        using std::exp; using std::sin; using std::cos;
        T vn = imagNorm();
        T ew = exp(e[0]);
        if (vn == T(0)) {
            return Octonion(ew);
        }
        T s = ew * sin(vn) / vn;
        T c = ew * cos(vn);
        return Octonion(c, e[1]*s, e[2]*s, e[3]*s,
                           e[4]*s, e[5]*s, e[6]*s, e[7]*s);
    }

    // log(o) = log|o| + v̂·acos(e0/|o|)
    [[nodiscard]] Octonion log() const {
        using std::log; using std::acos;
        T n = norm();
        T vn = imagNorm();
        if (vn == T(0)) {
            return Octonion(log(n));
        }
        T angle = acos(e[0] / n) / vn;
        return Octonion(log(n), e[1]*angle, e[2]*angle, e[3]*angle,
                                e[4]*angle, e[5]*angle, e[6]*angle, e[7]*angle);
    }

    // pow(o, t) = exp(t * log(o))  (スカラー指数)
    [[nodiscard]] Octonion pow(const T& t) const {
        return (log() * t).exp();
    }

    // ============================================================
    // 三角関数・双曲線関数
    // ============================================================
    // sin(o) = sin(a)cosh|v| + v̂·cos(a)sinh|v|
    // cos(o) = cos(a)cosh|v| - v̂·sin(a)sinh|v|

    [[nodiscard]] Octonion sin() const {
        using std::sin; using std::cos; using std::sinh; using std::cosh;
        T vn = imagNorm();
        if (vn == T(0)) return Octonion(sin(e[0]));
        T s = cos(e[0]) * sinh(vn) / vn;
        T c = sin(e[0]) * cosh(vn);
        return Octonion(c, e[1]*s, e[2]*s, e[3]*s,
                           e[4]*s, e[5]*s, e[6]*s, e[7]*s);
    }

    [[nodiscard]] Octonion cos() const {
        using std::sin; using std::cos; using std::sinh; using std::cosh;
        T vn = imagNorm();
        if (vn == T(0)) return Octonion(cos(e[0]));
        T s = -sin(e[0]) * sinh(vn) / vn;
        T c = cos(e[0]) * cosh(vn);
        return Octonion(c, e[1]*s, e[2]*s, e[3]*s,
                           e[4]*s, e[5]*s, e[6]*s, e[7]*s);
    }

    [[nodiscard]] Octonion tan() const {
        // tan(o) = sin(o) * cos(o)^(-1)  (右除算)
        return sin() * cos().inverse();
    }

    [[nodiscard]] Octonion sinh() const {
        using std::sin; using std::cos; using std::sinh; using std::cosh;
        T vn = imagNorm();
        if (vn == T(0)) return Octonion(sinh(e[0]));
        T s = cosh(e[0]) * sin(vn) / vn;
        T c = sinh(e[0]) * cos(vn);
        return Octonion(c, e[1]*s, e[2]*s, e[3]*s,
                           e[4]*s, e[5]*s, e[6]*s, e[7]*s);
    }

    [[nodiscard]] Octonion cosh() const {
        using std::sin; using std::cos; using std::sinh; using std::cosh;
        T vn = imagNorm();
        if (vn == T(0)) return Octonion(cosh(e[0]));
        T s = sinh(e[0]) * sin(vn) / vn;
        T c = cosh(e[0]) * cos(vn);
        return Octonion(c, e[1]*s, e[2]*s, e[3]*s,
                           e[4]*s, e[5]*s, e[6]*s, e[7]*s);
    }

    [[nodiscard]] Octonion tanh() const {
        return sinh() * cosh().inverse();
    }

    // ============================================================
    // 非結合性の検証ユーティリティ
    // ============================================================

    // 結合子 (associator): [a,b,c] = (a*b)*c - a*(b*c)
    // 八元数は alternative: [a,a,b] = [a,b,b] = 0
    static Octonion associator(const Octonion& a, const Octonion& b,
                               const Octonion& c) {
        return (a * b) * c - a * (b * c);
    }

    // ============================================================
    // 文字列変換・出力
    // ============================================================

    [[nodiscard]] std::string toString() const {
        static const char* units[] = {"", "i", "j", "k", "l", "li", "lj", "lk"};
        auto fmt = [](const T& v) -> std::string {
            std::ostringstream os;
            os << ((v == T(0)) ? T(0) : v);
            return os.str();
        };

        std::string s;
        for (int i = 0; i < 8; ++i) {
            std::string val = fmt(e[i]);
            if (val == "0" || val == "0.") continue;

            if (i == 0) {
                s += val;
            } else {
                if (val == "1") { s += "+"; s += units[i]; }
                else if (val == "-1") { s += "-"; s += units[i]; }
                else if (!val.empty() && val[0] == '-') { s += val; s += units[i]; }
                else { s += "+"; s += val; s += units[i]; }
            }
        }

        if (s.empty()) return "0";
        if (!s.empty() && s[0] == '+') s.erase(0, 1);
        return s;
    }

    friend std::ostream& operator<<(std::ostream& os, const Octonion& o) {
        return os << o.toString();
    }
};

// ================================================================
// 自由関数
// ================================================================

template<OctonionScalar T>
[[nodiscard]] constexpr Octonion<T> conj(const Octonion<T>& o) { return o.conj(); }

template<OctonionScalar T>
[[nodiscard]] constexpr T normSq(const Octonion<T>& o) { return o.normSq(); }

template<OctonionScalar T>
[[nodiscard]] T norm(const Octonion<T>& o) { return o.norm(); }

template<OctonionScalar T>
[[nodiscard]] constexpr Octonion<T> inverse(const Octonion<T>& o) { return o.inverse(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> normalize(const Octonion<T>& o) { return o.normalized(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> exp(const Octonion<T>& o) { return o.exp(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> log(const Octonion<T>& o) { return o.log(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> sin(const Octonion<T>& o) { return o.sin(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> cos(const Octonion<T>& o) { return o.cos(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> tan(const Octonion<T>& o) { return o.tan(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> sinh(const Octonion<T>& o) { return o.sinh(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> cosh(const Octonion<T>& o) { return o.cosh(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> tanh(const Octonion<T>& o) { return o.tanh(); }

template<OctonionScalar T>
[[nodiscard]] Octonion<T> pow(const Octonion<T>& o, const T& t) { return o.pow(t); }

// 型特性
template<typename T>
struct is_octonion : std::false_type {};

template<OctonionScalar T>
struct is_octonion<Octonion<T>> : std::true_type {};

template<typename T>
inline constexpr bool is_octonion_v = is_octonion<T>::value;

} // namespace calx
