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

// ModularInt.hpp
#ifndef SANGI_MODULAR_INT_HPP
#define SANGI_MODULAR_INT_HPP
// Modular integer class (C++20/23 ready, concept-aware version)
//
// This file defines a class template that provides efficient modular integer
// arithmetic.
// ModularInt<P> represents an element of the residue field Z/PZ for the
// prime modulus P.
// It supports basic arithmetic operations (addition, subtraction,
// multiplication, division), exponentiation, and comparison operations.
//
// Main features:
// - Modulus is specified through a template parameter
// - Integrates with C++20 concepts
// - Efficient algorithms
// - Integrates with algebraic structures (monoid, group, ring, field)
// - Stream I/O support
//
// Usage example:
// ```
// 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 sangi {

// Modular integer class (C++20 concepts version)
template <int P>
class ModularInt {
private:
    int value_;

    // Simple implementation using a 64-bit integer in place of Barrett
    // reduction to avoid overflow during multiplication
    static inline constexpr int mul_mod(int64_t a, int64_t b) {
        return static_cast<int>((a * b) % P);
    }

    // Fast modulo of a 64-bit integer
    static inline constexpr int fast_mod(int64_t x) {
        return static_cast<int>(x % P);
    }

public:
    // Default constructor (initializes to 0)
    constexpr ModularInt() noexcept : value_(0) {}

    // Converting constructor from an integer value
    template <std::integral T>
    constexpr ModularInt(T v) noexcept : value_(fast_mod(static_cast<int64_t>(v) % P + P)) {}

    // Copy constructor
    constexpr ModularInt(const ModularInt& other) noexcept = default;

    // Move constructor
    constexpr ModularInt(ModularInt&& other) noexcept = default;

    // Assignment operator
    constexpr ModularInt& operator=(const ModularInt& other) noexcept = default;

    // Move assignment operator
    constexpr ModularInt& operator=(ModularInt&& other) noexcept = default;

    // Assignment from an integer type
    template <std::integral T>
    constexpr ModularInt& operator=(T v) noexcept {
        value_ = fast_mod(static_cast<int64_t>(v) % P + P);
        return *this;
    }

    // Destructor
    ~ModularInt() = default;

    // Accessor
    constexpr int value() const noexcept { return value_; }

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

    // Addition with an integer
    template <std::integral T>
    constexpr ModularInt operator+(T other) const noexcept {
        return *this + ModularInt(other);
    }

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

    // Subtraction with an integer
    template <std::integral T>
    constexpr ModularInt operator-(T other) const noexcept {
        return *this - ModularInt(other);
    }

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

    // Multiplication with an integer
    template <std::integral T>
    constexpr ModularInt operator*(T other) const noexcept {
        return *this * ModularInt(other);
    }

    // Division (uses the modular inverse)
    ModularInt operator/(const ModularInt& other) const {
        if (other.value_ == 0) {
            throw MathError("Division by zero in ModularInt");
        }
        return (*this) * other.inverse();
    }

    // Division by an integer
    template <std::integral T>
    ModularInt operator/(T other) const {
        return *this / ModularInt(other);
    }

    // Unary minus operator
    constexpr ModularInt operator-() const noexcept {
        return ModularInt(value_ == 0 ? 0 : P - value_);
    }

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

    // Compound addition with an integer
    template <std::integral T>
    constexpr ModularInt& operator+=(T other) noexcept {
        return *this += ModularInt(other);
    }

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

    // Compound subtraction with an integer
    template <std::integral T>
    constexpr ModularInt& operator-=(T other) noexcept {
        return *this -= ModularInt(other);
    }

    // Compound multiplication
    constexpr ModularInt& operator*=(const ModularInt& other) noexcept {
        value_ = mul_mod(value_, other.value_);
        return *this;
    }

    // Compound multiplication with an integer
    template <std::integral T>
    constexpr ModularInt& operator*=(T other) noexcept {
        return *this *= ModularInt(other);
    }

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

    // Compound division by an integer
    template <std::integral T>
    ModularInt& operator/=(T other) {
        return *this /= ModularInt(other);
    }

    // Pre-increment
    constexpr ModularInt& operator++() noexcept {
        ++value_;
        if (value_ == P) value_ = 0;
        return *this;
    }

    // Post-increment
    constexpr ModularInt operator++(int) noexcept {
        ModularInt old = *this;
        ++(*this);
        return old;
    }

    // Pre-decrement
    constexpr ModularInt& operator--() noexcept {
        if (value_ == 0) value_ = P;
        --value_;
        return *this;
    }

    // Post-decrement
    constexpr ModularInt operator--(int) noexcept {
        ModularInt old = *this;
        --(*this);
        return old;
    }

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

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

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

    // Inequality with an integer
    template <std::integral T>
    constexpr bool operator!=(T other) const noexcept {
        return *this != ModularInt(other);
    }

    // Less-than operator
    constexpr bool operator<(const ModularInt& other) const noexcept {
        return value_ < other.value_;
    }

    // Less-than with an integer
    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);
    }

    // Implicit conversion to an integer type
    constexpr operator int() const noexcept { return value_; }

    // Inverse computation (extended Euclidean algorithm)
    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;
        }

        // If the result is negative, convert to a positive value
        return ModularInt(u < 0 ? u + P : u);
    }

    // Exponentiation (repeated squaring)
    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;
    }

    // Stream output
    friend std::ostream& operator<<(std::ostream& os, const ModularInt& m) {
        return os << m.value_;
    }

    // Stream input
    friend std::istream& operator>>(std::istream& is, ModularInt& m) {
        int val;
        is >> val;
        m = ModularInt(val);
        return is;
    }
};

// Free functions

// Addition (right-hand side is an integer)
template <int P, std::integral T>
constexpr ModularInt<P> operator+(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) + right;
}

// Subtraction (right-hand side is an integer)
template <int P, std::integral T>
constexpr ModularInt<P> operator-(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) - right;
}

// Multiplication (right-hand side is an integer)
template <int P, std::integral T>
constexpr ModularInt<P> operator*(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) * right;
}

// Division (right-hand side is an integer)
template <int P, std::integral T>
ModularInt<P> operator/(T left, const ModularInt<P>& right) {
    return ModularInt<P>(left) / right;
}

// Equality (left-hand side is an integer)
template <int P, std::integral T>
constexpr bool operator==(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) == right;
}

// Inequality (left-hand side is an integer)
template <int P, std::integral T>
constexpr bool operator!=(T left, const ModularInt<P>& right) noexcept {
    return ModularInt<P>(left) != right;
}

// Ordering (left-hand side is an integer)
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;
}

// Absolute value
template <int P>
constexpr int abs(const ModularInt<P>& m) noexcept {
    return m.value();  // Modular integers are always non-negative, so just return the value
}

// Fast exponentiation
template <int P>
constexpr ModularInt<P> pow(const ModularInt<P>& base, int64_t exponent) {
    return base.pow(exponent);
}

// Specialization that confirms ModularInt's concept compatibility
template <int P>
struct numeric_traits<ModularInt<P>> {
    using value_type = ModularInt<P>;
    using category = floating_point_tag;  // Behaves like a real-number-like type

    static constexpr bool is_supported = true;
    static constexpr bool is_complex = false;
    static constexpr bool is_integer = false;  // Treated as a field
    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; // In a finite field, all non-zero elements are equally good
    }
};

} // namespace sangi

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

#endif // SANGI_MODULAR_INT_HPP