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

// IntSequence.hpp
// 整数列（Fibonacci 数、Lucas 数）
// アルゴリズム: square-only fast doubling, O(S(n) log n)
//   S(n) = 自乗コスト。一般乗算を完全に除去。

#pragma once

#include "IntBase.hpp"
#include "IntOps.hpp"
#include "MpnOps.hpp"

#include <cstdint>
#include <cstring>
#include <utility>

namespace sangi {

/**
 * @brief Fibonacci 数・Lucas 数
 *
 * Square-only fast doubling: (F(k), F(k-1)) ペアを追跡し、
 * 各ステップ 2 回の自乗 + 加減算のみで計算 (一般乗算なし)。
 *
 * 公式:
 *   F(2k-1) = F(k)^2 + F(k-1)^2
 *   F(2k+1) = 4*F(k)^2 - F(k-1)^2 + 2*(-1)^k
 *   F(2k)   = F(2k+1) - F(2k-1)
 *
 * GMP 対応: mpz_fib_ui, mpz_fib2_ui, mpz_lucnum_ui, mpz_lucnum2_ui
 */
class IntSequence {
public:
    /**
     * @brief Fibonacci 数 F(n)
     *
     * F(0)=0, F(1)=1, F(n)=F(n-1)+F(n-2)
     *
     * @param n 非負整数
     * @return F(n)
     */
    static Int fibonacci(uint64_t n) {
        return fibonacci2(n).first;
    }

    /**
     * @brief Fibonacci 数のペア {F(n), F(n-1)}
     *
     * @param n 非負整数 (n=0 のとき {F(0), F(-1)} = {0, 1})
     * @return {F(n), F(n-1)}
     */
    static std::pair<Int, Int> fibonacci2(uint64_t n) {
        if (n == 0) return { Int(0), Int(1) };
        if (n == 1) return { Int(1), Int(0) };

        // F(93) は uint64_t に収まる最大の Fibonacci 数
        static constexpr uint64_t fib_table[] = {
            0,1,1,2,3,5,8,13,21,34,55,89,144,233,377,610,987,
            1597,2584,4181,6765,10946,17711,28657,46368,75025,
            121393,196418,317811,514229,832040,1346269,2178309,
            3524578,5702887,9227465,14930352,24157817,39088169,
            63245986,102334155,165580141,267914296,433494437,
            701408733,1134903170,1836311903,2971215073ULL,
            4807526976ULL,7778742049ULL,12586269025ULL,
            20365011074ULL,32951280099ULL,53316291173ULL,
            86267571272ULL,139583862445ULL,225851433717ULL,
            365435296162ULL,591286729879ULL,956722026041ULL,
            1548008755920ULL,2504730781961ULL,4052739537881ULL,
            6557470319842ULL,10610209857723ULL,17167680177565ULL,
            27777890035288ULL,44945570212853ULL,72723460248141ULL,
            117669030460994ULL,190392490709135ULL,308061521170129ULL,
            498454011879264ULL,806515533049393ULL,1304969544928657ULL,
            2111485077978050ULL,3416454622906707ULL,5527939700884757ULL,
            8944394323791464ULL,14472334024676221ULL,23416728348467685ULL,
            37889062373143906ULL,61305790721611591ULL,99194853094755497ULL,
            160500643816367088ULL,259695496911122585ULL,420196140727489673ULL,
            679891637638612258ULL,1100087778366101931ULL,
            1779979416004714189ULL,2880067194370816120ULL,
            4660046610375530309ULL,7540113804746346429ULL,
            12200160415121876738ULL,
        };

        // 小さい n: テーブルルックアップ
        if (n <= 93) {
            return { Int(fib_table[n]), Int(fib_table[n > 0 ? n - 1 : 0]) };
        }

        // n <= 50000: mpn 直接版 (GMP 方式 3 バッファ + テーブル開始)
        // n > 50000: IntOps 版 (大サイズでは square が支配的で差が小さい)
        if (n > 50000) {
            int bits = 0;
            for (uint64_t tmp = n; tmp > 0; tmp >>= 1) ++bits;
            return fibonacci2_intops(n, bits);
        }

        // テーブル開始: n をテーブル範囲まで半減し、反復回数を削減
        // GMP 方式: nfirst = n >> steps, そこから steps 回 doubling
        uint64_t nfirst = n;
        int steps = 0;
        while (nfirst > 93) { nfirst /= 2; steps++; }

        // F(n) のサイズ推定: ~0.694*n bits = ceil(0.694*n/64) limbs
        size_t max_limbs = static_cast<size_t>(static_cast<double>(n) * 0.694 / 64.0) + 8;
        if (max_limbs < 8) max_limbs = 8;

        // 3 バッファ (fp, f1p, xp) + scratch
        // GMP 方式: fp=F[k], f1p=F[k-1], xp=scratch (ポインタスワップ不要)
        size_t alloc = 2 * max_limbs + 4;
        size_t scratch_sz = mpn::square_scratch_size(max_limbs);
        size_t total = 3 * alloc + scratch_sz;

        // 小サイズはスタックバッファ (ヒープ確保+ゼロ初期化を回避)
        static constexpr size_t STACK_LIMIT = 2048;  // 16KB
        uint64_t stack_mem[STACK_LIMIT];
        std::vector<uint64_t> heap_mem;
        uint64_t* mem_ptr;
        if (total <= STACK_LIMIT) {
            mem_ptr = stack_mem;
        } else {
            heap_mem.resize(total);
            mem_ptr = heap_mem.data();
        }

        uint64_t* fp  = mem_ptr;                // F[k]
        uint64_t* f1p = mem_ptr + alloc;         // F[k-1]
        uint64_t* xp  = mem_ptr + 2 * alloc;    // scratch (F[k]² を保持)
        uint64_t* scratch = mem_ptr + 3 * alloc;

        // テーブルから初期値を設定
        fp[0] = fib_table[nfirst];
        f1p[0] = fib_table[nfirst > 0 ? nfirst - 1 : 0];
        size_t size = 1;

        // GMP 方式 doubling ループ:
        //   xp ← fp²          (F[k]²)
        //   fp ← f1p²         (F[k-1]², fp を上書き — 入力は f1p なので安全)
        //   f1p ← xp + fp     (F[2k-1] = F[k]² + F[k-1]²)
        //   fp ← 4*xp - fp ± 2 (F[2k+1] = 4·F[k]² - F[k-1]² + 2·(-1)^k)
        //   bit=1: f1p ← fp - f1p  → (fp, f1p) = (F[2k+1], F[2k])
        //   bit=0: fp ← fp - f1p   → (fp, f1p) = (F[2k], F[2k-1])
        for (int i = steps - 1; i >= 0; --i) {
            bool k_odd = (n >> (i + 1)) & 1;  // nfirst の対応ビット = 現在の k の LSB
            // 注: k の値は n >> i (上位ビット群) なので、k の LSB は n >> (i+1) の LSB ではなく
            //     n の bit (i) の上のビット列の最下位。GMP と同様に n & (1 << (i+1)) で判定。
            //     ただし i+1 が 64 以上になることはない (n <= 50000)。

            // xp = fp² (F[k]²)
            mpn::square(xp, fp, size, scratch);
            // fp = f1p² (F[k-1]²) — fp ≠ f1p なので安全
            mpn::square(fp, f1p, size, scratch);
            size *= 2;
            size -= (xp[size - 1] == 0);

            // f1p = xp + fp (F[2k-1])
            f1p[size] = mpn::add(f1p, xp, size, fp, size);

            // F[2k+1] = 4·F[k]² - F[k-1]² + 2·(-1)^k
            // GMP のビットOR トリック: 平方数の bit1 は常に 0 なので |= 2 で ±2 を吸収
            if (k_odd) {
                // (-1)^k = -1 → -2: fp (= F[k-1]²) に 2 を OR して被減数を増やす
                fp[0] |= 2;
            }
            uint64_t cy = mpn::lshift(xp, xp, size, 2);  // 4·F[k]²
            if (!k_odd) {
                // (-1)^k = +1 → +2: xp (= 4·F[k]²) に 2 を OR
                xp[0] |= 2;
            }
            // fp = xp - fp (= 4·F[k]² ± 2 - F[k-1]² ∓ 2 = F[2k+1])
            cy -= mpn::sub(fp, xp, size, fp, size);
            fp[size] = cy;
            size += (fp[size] != 0);

            // F[2k] = F[2k+1] - F[2k-1]: 不要な方を置き換え
            if ((n >> i) & 1) {
                // bit=1: (fp, f1p) → (F[2k+1], F[2k])
                mpn::sub(f1p, fp, size, f1p, size);
            } else {
                // bit=0: (fp, f1p) → (F[2k], F[2k-1])
                mpn::sub(fp, fp, size, f1p, size);
                // F[2k] < F[2k+1] なので先頭ゼロが出る可能性
                while (size > 0 && fp[size - 1] == 0) --size;
            }
        }

        // Int 構築 (最終結果のみ): fp=F(n), f1p=F(n-1)
        // span + PreNormalized で一時 vector 確保を回避
        while (size > 0 && fp[size - 1] == 0) --size;
        size_t f1size = size;
        while (f1size > 0 && f1p[f1size - 1] == 0) --f1size;

        Int result_fn, result_fn1;
        if (size > 0) {
            result_fn = Int::fromRawWordsPreNormalized(
                std::span<const uint64_t>(fp, size), 1);
        }
        if (f1size > 0) {
            result_fn1 = Int::fromRawWordsPreNormalized(
                std::span<const uint64_t>(f1p, f1size), 1);
        }
        return { std::move(result_fn), std::move(result_fn1) };
    }

private:
    // 大サイズ用: IntOps square-only fast doubling
    static std::pair<Int, Int> fibonacci2_intops(uint64_t n, int bits) {
        Int a(1), b(0);  // F(1), F(0)
        bool neg = true;
        Int a2, b2, f2km1, f2kp1, f2k;
        static const Int TWO(2);

        for (int i = bits - 2; i >= 0; --i) {
            IntOps::square(a, a2);
            IntOps::square(b, b2);

            // F(2k-1) = a² + b²
            IntOps::addUnchecked(a2, b2, f2km1);

            // F(2k+1) = 4*a² - b² + 2*(-1)^k
            IntOps::leftShift(a2, 2, f2kp1);
            IntOps::subUnchecked(f2kp1, b2, f2kp1);
            if (neg) {
                IntOps::subUnchecked(f2kp1, TWO, f2kp1);
            } else {
                IntOps::addUnchecked(f2kp1, TWO, f2kp1);
            }

            // F(2k) = F(2k+1) - F(2k-1)
            IntOps::subUnchecked(f2kp1, f2km1, f2k);

            if ((n >> i) & 1) {
                a = std::move(f2kp1);
                b = std::move(f2k);
                neg = true;
            } else {
                a = std::move(f2k);
                b = std::move(f2km1);
                neg = false;
            }
        }
        return { std::move(a), std::move(b) };
    }

public:
    /**
     * @brief Lucas 数 L(n)
     *
     * L(0)=2, L(1)=1, L(n)=L(n-1)+L(n-2)
     * L(n) = F(n-1) + F(n+1) = 2*F(n+1) - F(n)
     *
     * @param n 非負整数
     * @return L(n)
     */
    static Int lucas(uint64_t n) {
        if (n == 0) return Int(2);
        auto [fn, fn_1] = fibonacci2(n);
        // L(n) = F(n) + 2*F(n-1)
        return fn + fn_1 + fn_1;
    }

    /**
     * @brief Lucas 数のペア {L(n), L(n-1)}
     *
     * @param n 非負整数 (n=0 のとき {L(0), L(-1)} = {2, -1})
     * @return {L(n), L(n-1)}
     */
    static std::pair<Int, Int> lucas2(uint64_t n) {
        if (n == 0) return { Int(2), Int(-1) };
        auto [fn, fn_1] = fibonacci2(n);
        // L(n) = F(n) + 2*F(n-1)
        // L(n-1) = 2*F(n) - F(n-1)
        Int ln = fn + fn_1 + fn_1;
        Int ln_1 = fn + fn - fn_1;
        return { std::move(ln), std::move(ln_1) };
    }
};

} // namespace sangi