// power.hpp
// Generic exponentiation function (right-to-left binary method)
//
// Supports any type with multiplication (integers, floating point, ModularInt, polynomials, matrices, etc.)
//
// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

#pragma once

#include <concepts>
#include <bit>
#include <cassert>
#include <type_traits>

namespace sangi {

    // For types where the multiplicative identity can be obtained via T(1)
    // (integers, floating point, ModularInt, polynomials, etc.)
    template<typename T>
        requires requires(T a, const T& b) { a = a * b; }
    [[nodiscard]] T power(const T& base, unsigned int exponent) {
        if (exponent == 0) return T(1);
        if (exponent == 1) return base;

        // Right-to-left binary method (scan from MSB to LSB)
        // About 3x faster than left-to-right binary method (sangi measurement)
        T y(1);
        unsigned int bit = std::bit_width(exponent) - 1;
        while (true) {
            y = y * y;
            if ((exponent >> bit) & 1u) {
                y = y * base;
            }
            if (bit == 0) break;
            --bit;
        }
        return y;
    }

    // Version that takes the multiplicative identity explicitly
    // For types where T(1) is unavailable (Matrix, etc.)
    template<typename T>
        requires requires(T a, const T& b) { a = a * b; }
    [[nodiscard]] T power(const T& base, unsigned int exponent, const T& identity) {
        if (exponent == 0) return identity;
        if (exponent == 1) return base;

        T y = identity;
        unsigned int bit = std::bit_width(exponent) - 1;
        while (true) {
            y = y * y;
            if ((exponent >> bit) & 1u) {
                y = y * base;
            }
            if (bit == 0) break;
            --bit;
        }
        return y;
    }

} // namespace sangi