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

// ModularInt.hpp
#ifndef CALX_MODULAR_INT_HPP
#define CALX_MODULAR_INT_HPP
// モジュラー整数クラス（C++20/23対応、コンセプト対応版）
//
// このファイルは、効率的なモジュラー整数演算を提供するクラステンプレートを定義します。
// ModularInt<P>は、素数Pを法とする剰余体Z/PZの要素を表現します。
// 基本的な算術演算（加算、減算、乗算、除算）、冪乗演算、比較演算をサポートします。
//
// 主な特徴:
// - テンプレートパラメータで法（モジュラス）を指定
// - C++20のコンセプトとの統合
// - 効率的な演算アルゴリズム
// - 代数的構造との統合（モノイド、群、リング、体）
// - ストリーム入出力サポート
//
// 使用例:
// ```
// ModularInt<1000000007> a = 42;  // 42 mod 10^9+7
// ModularInt<1000000007> b = 987654321;
// auto c = a * b;  // 41477417 mod 10^9+7
// auto d = a.pow(12345);  // a^12345 mod 10^9+7
// ```

#include <math/concepts/algebraic_concepts.hpp>
#include <math/core/common.hpp>
#include <math/core/traits.hpp>

#include <iosfwd>
#include <type_traits>
#include <cstdint>
#include <concepts>

namespace calx {

// モジュラー整数クラス（C++20コンセプト対応版）
template <int P>
class ModularInt {
private:
    int value_;

    // 乗算時のオーバーフロー回避のためのバレット削減法の代わりに
    // 64ビット整数を使用した簡易実装
    static inline constexpr int mul_mod(int64_t a, int64_t b) {
        return static_cast<int>((a * b) % P);
    }

    // 64ビット整数の高速なmod計算
    static inline constexpr int fast_mod(int64_t x) {
        return static_cast<int>(x % P);
    }

public:
    // デフォルトコンストラクタ（0で初期化）
    constexpr ModularInt() noexcept : value_(0) {}

    // 整数値からの変換コンストラクタ
    template <std::integral T>
    constexpr ModularInt(T v) noexcept : value_(fast_mod(static_cast<int64_t>(v) % P + P)) {}

    // コピーコンストラクタ
    constexpr ModularInt(const ModularInt& other) noexcept = default;

    // ムーブコンストラクタ
    constexpr ModularInt(ModularInt&& other) noexcept = default;

    // 代入演算子
    constexpr ModularInt& operator=(const ModularInt& other) noexcept = default;

    // ムーブ代入演算子
    constexpr ModularInt& operator=(ModularInt&& other) noexcept = default;

    // 整数型への代入演算子
    template <std::integral T>
    constexpr ModularInt& operator=(T v) noexcept {
        value_ = fast_mod(static_cast<int64_t>(v) % P + P);
        return *this;
    }

    // デストラクタ
    ~ModularInt() = default;

    // アクセサ
    constexpr int value() const noexcept { return value_; }

    // 加算
    constexpr ModularInt operator+(const ModularInt& other) const noexcept {
        int sum = value_ + other.value_;
        return ModularInt(sum >= P ? sum - P : sum);
    }

    // 整数との加算
    template <std::integral T>
    constexpr ModularInt operator+(T other) const noexcept {
        return *this + ModularInt(other);
    }

    // 減算
    constexpr ModularInt operator-(const ModularInt& other) const noexcept {
        int diff = value_ - other.value_;
        return ModularInt(diff < 0 ? diff + P : diff);
    }

    // 整数との減算
    template <std::integral T>
    constexpr ModularInt operator-(T other) const noexcept {
        return *this - ModularInt(other);
    }

    // 乗算
    constexpr ModularInt operator*(const ModularInt& other) const noexcept {
        return ModularInt(mul_mod(value_, other.value_));
    }

    // 整数との乗算
    template <std::integral T>
    constexpr ModularInt operator*(T other) const noexcept {
        return *this * ModularInt(other);
    }

    // 除算（逆元を利用）
    ModularInt operator/(const ModularInt& other) const {
        if (other.value_ == 0) {
            throw MathError("Division by zero in ModularInt");
        }
        return (*this) * other.inverse();
    }

    // 整数との除算
    template <std::integral T>
    ModularInt operator/(T other) const {
        return *this / ModularInt(other);
    }

    // 単項マイナス演算子
    constexpr ModularInt operator-() const noexcept {
        return ModularInt(value_ == 0 ? 0 : P - value_);
    }

    // 累積加算
    constexpr ModularInt& operator+=(const ModularInt& other) noexcept {
        value_ += other.value_;
        if (value_ >= P) value_ -= P;
        return *this;
    }

    // 整数との累積加算
    template <std::integral T>
    constexpr ModularInt& operator+=(T other) noexcept {
        return *this += ModularInt(other);
    }

    // 累積減算
    constexpr ModularInt& operator-=(const ModularInt& other) noexcept {
        value_ -= other.value_;
        if (value_ < 0) value_ += P;
        return *this;
    }

    // 整数との累積減算
    template <std::integral T>
    constexpr ModularInt& operator-=(T other) noexcept {
        return *this -= ModularInt(other);
    }

    // 累積乗算
    constexpr ModularInt& operator*=(const ModularInt& other) noexcept {
        value_ = mul_mod(value_, other.value_);
        return *this;
    }

    // 整数との累積乗算
    template <std::integral T>
    constexpr ModularInt& operator*=(T other) noexcept {
        return *this *= ModularInt(other);
    }

    // 累積除算
    ModularInt& operator/=(const ModularInt& other) {
        if (other.value_ == 0) {
            throw MathError("Division by zero in ModularInt");
        }
        return (*this) *= other.inverse();
    }

    // 整数との累積除算
    template <std::integral T>
    ModularInt& operator/=(T other) {
        return *this /= ModularInt(other);
    }

    // 前置インクリメント
    constexpr ModularInt& operator++() noexcept {
        ++value_;
        if (value_ == P) value_ = 0;
        return *this;
    }

    // 後置インクリメント
    constexpr ModularInt operator++(int) noexcept {
        ModularInt old = *this;
        ++(*this);
        return old;
    }

    // 前置デクリメント
    constexpr ModularInt& operator--() noexcept {
        if (value_ == 0) value_ = P;
        --value_;
        return *this;
    }

    // 後置デクリメント
    constexpr ModularInt operator--(int) noexcept {
        ModularInt old = *this;
        --(*this);
        return old;
    }

    // 等価比較演算子
    constexpr bool operator==(const ModularInt& other) const noexcept {
        return value_ == other.value_;
    }

    // 整数との等価比較演算子
    template <std::integral T>
    constexpr bool operator==(T other) const noexcept {
        return *this == ModularInt(other);
    }

    // 不等価比較演算子
    constexpr bool operator!=(const ModularInt& other) const noexcept {
        return value_ != other.value_;
    }

    // 整数との不等価比較演算子
    template <std::integral T>
    constexpr bool operator!=(T other) const noexcept {
        return *this != ModularInt(other);
    }

    // 大小比較演算子
    constexpr bool operator<(const ModularInt& other) const noexcept {
        return value_ < other.value_;
    }

    // 整数との大小比較演算子
    template <std::integral T>
    constexpr bool operator<(T other) const noexcept {
        return *this < ModularInt(other);
    }

    constexpr bool operator<=(const ModularInt& other) const noexcept {
        return value_ <= other.value_;
    }

    template <std::integral T>
    constexpr bool operator<=(T other) const noexcept {
        return *this <= ModularInt(other);
    }

    constexpr bool operator>(const ModularInt& other) const noexcept {
        return value_ > other.value_;
    }

    template <std::integral T>
    constexpr bool operator>(T other) const noexcept {
        return *this > ModularInt(other);
    }

    constexpr bool operator>=(const ModularInt& other) const noexcept {
        return value_ >= other.value_;
    }

    template <std::integral T>
    constexpr bool operator>=(T other) const noexcept {
        return *this >= ModularInt(other);
    }

    // 整数型への暗黙的変換
    constexpr operator int() const noexcept { return value_; }

    // 逆元の計算（拡張ユークリッドアルゴリズム）
    constexpr ModularInt inverse() const {
        if (value_ == 0) {
            throw MathError("Division by zero in ModularInt::inverse");
        }

        int a = value_;
        int b = P;
        int u = 1, v = 0;

        while (b > 0) {
            int q = a / b;
            std::swap(a, b);
            b = b - q * a;
            std::swap(u, v);
            v = v - q * u;
        }

        // 結果が負の場合は正の値に変換
        return ModularInt(u < 0 ? u + P : u);
    }

    // 冪乗演算（繰り返し二乗法）
    constexpr ModularInt pow(int64_t n) const {
        if (n < 0) {
            return inverse().pow(-n);
        }

        ModularInt res(1);
        ModularInt base(*this);

        while (n > 0) {
            if (n & 1) res *= base;
            base *= base;
            n >>= 1;
        }

        return res;
    }

    // ストリーム出力
    friend std::ostream& operator<<(std::ostream& os, const ModularInt& m) {
        return os << m.value_;
    }

    // ストリーム入力
    friend std::istream& operator>>(std::istream& is, ModularInt& m) {
        int val;
        is >> val;
        m = ModularInt(val);
        return is;
    }
};

// フリー関数

// 加算（右辺が整数）
template <int P, std::integral T>
constexpr ModularInt<P> operator+(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) + right;
}

// 減算（右辺が整数）
template <int P, std::integral T>
constexpr ModularInt<P> operator-(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) - right;
}

// 乗算（右辺が整数）
template <int P, std::integral T>
constexpr ModularInt<P> operator*(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) * right;
}

// 除算（右辺が整数）
template <int P, std::integral T>
ModularInt<P> operator/(T left, const ModularInt<P>& right) {
    return ModularInt<P>(left) / right;
}

// 等価比較（左辺が整数）
template <int P, std::integral T>
constexpr bool operator==(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) == right;
}

// 不等価比較（左辺が整数）
template <int P, std::integral T>
constexpr bool operator!=(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) != right;
}

// 大小比較（左辺が整数）
template <int P, std::integral T>
constexpr bool operator<(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) < right;
}

template <int P, std::integral T>
constexpr bool operator<=(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) <= right;
}

template <int P, std::integral T>
constexpr bool operator>(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) > right;
}

template <int P, std::integral T>
constexpr bool operator>=(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) >= right;
}

// 絶対値関数
template <int P>
constexpr int abs(const ModularInt<P>& m) noexcept {
    return m.value();  // モジュラー整数は常に非負なので、値をそのまま返す
}

// 高速冪乗計算
template <int P>
constexpr ModularInt<P> pow(const ModularInt<P>& base, int64_t exponent) {
    return base.pow(exponent);
}

// モジュラー整数のコンセプト対応を確認するための特殊化
template <int P>
struct numeric_traits<ModularInt<P>> {
    using value_type = ModularInt<P>;
    using category = floating_point_tag;  // 実数風の振る舞い
    
    static constexpr bool is_supported = true;
    static constexpr bool is_complex = false;
    static constexpr bool is_integer = false;  // 体として扱うため
    static constexpr bool is_floating_point = false;
    
    static constexpr ModularInt<P> zero() { return ModularInt<P>(0); }
    static constexpr ModularInt<P> one() { return ModularInt<P>(1); }
    static constexpr ModularInt<P> epsilon() { return ModularInt<P>(1); }
    
    static constexpr ModularInt<P> abs(const ModularInt<P>& value) { return value; }
    static constexpr ModularInt<P> conj(const ModularInt<P>& value) { return value; }
    static constexpr int norm(const ModularInt<P>& value) { return static_cast<int>(static_cast<int64_t>(value.value()) * value.value() % P); }
    static bool pivotBetter(const ModularInt<P>& a, const ModularInt<P>&) {
        return a.value() != 0; // 有限体では全非ゼロ要素が同等
    }
};

} // namespace calx

namespace std {
    // std::hash特殊化
    template <int P>
    struct hash<calx::ModularInt<P>> {
        std::size_t operator()(const calx::ModularInt<P>& m) const noexcept {
            return std::hash<int>{}(m.value());
        }
    };
}

#endif // CALX_MODULAR_INT_HPP