// power.hpp
// 汎用冪乗関数（右向きバイナリ法）
//
// 移植元: lib++20 基本関数.h Pow()
// 任意の乗法を持つ型（整数, 浮動小数点, ModularInt, 多項式, 行列, 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 calx {

    // T(1) で乗法単位元を取得できる型向け
    // 整数, 浮動小数点, ModularInt, 多項式 等
    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;

        // 右向きバイナリ法（MSB から LSB へ走査）
        // 左向きバイナリ法より約3倍高速（lib++20 実測値）
        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;
    }

    // 乗法単位元を明示指定するバージョン
    // T(1) が使えない型（Matrix 等）向け
    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 calx