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

// IntFactorTable.cpp
// 素因数テーブルの実装

#include "math/core/mp/Int/IntFactorTable.hpp"
#include "math/core/mp/Int/IntOps.hpp"
#include "math/core/mp/Int/IntPrime.hpp"
#include "math/core/mp/Int/IntGCD.hpp"
#include <cmath>
#include <fstream>
#include <algorithm>

namespace calx {

// コンストラクタ
IntFactorTable::IntFactorTable(int max_value)
    : max_value_(max_value)
{
    if (max_value_ <= 0) {
        max_value_ = 0;
        return;
    }

    // 奇数のみ格納（偶数は素因数 2 で固定）
    // 奇数 n のインデックスは n/2
    size_t table_size = (max_value_ + 1) / 2;
    table_.resize(table_size, 0);

    // ファイルから読み込みを試みる
    if (max_value_ >= 100000) {  // 大きいテーブルの場合のみファイル使用
        if (loadFromFile()) {
            return;  // 読み込み成功
        }
    }

    // 読み込み失敗 or ファイル未使用 → テーブルを構築
    buildTable();

    // ファイルに保存（次回起動時の高速化）
    if (max_value_ >= 100000) {
        saveToFile();
    }
}

// エラトステネスの篩でテーブルを構築（最大素因数版）
void IntFactorTable::buildTable() {
    // table[n/2] に奇数 n の最大素因数を格納
    // 0 は「素数」または「未処理」を意味する

    if (max_value_ < 1) return;

    // 1 の素因数は 1（特殊ケース）
    if (max_value_ >= 1) {
        table_[oddIndex(1)] = 1;
    }

    // sqrt(max_value) まで篩にかける
    int sqrt_max = static_cast<int>(std::sqrt(static_cast<double>(max_value_)));

    // すべての奇数素数 p について、その倍数に p を書き込む
    // 小さい素数から順に処理するため、最後に書き込まれた値が最大素因数になる
    for (int p = 3; p <= sqrt_max; p += 2) {
        // p が素数かチェック（table[p] == 0 なら素数）
        if (table_[oddIndex(p)] == 0) {
            // p の倍数に対して最大素因数 p を記録
            // p の最初の奇数合成数倍から開始（p 自身は素数なのでスキップし、p*2 から）
            int start = 2 * p;  // p の最初の合成数倍（偶数）
            if ((start & 1) == 0) {
                start += p;  // 奇数にする（p*2 が偶数なら p*3 にする）
            }

            // p の倍数を 2*p ずつ進む（奇数のみ）
            for (int m = start; m <= max_value_; m += 2 * p) {
                // ★ LPF の特徴: 条件判定なしで無条件に上書き ★
                table_[oddIndex(m)] = static_cast<uint32_t>(p);
            }
        }
    }

    // sqrt_max より大きい素数 p について、その倍数を処理
    // これらは p^2 > max_value なので、p の小さい倍数のみが範囲内
    for (int p = sqrt_max + (sqrt_max % 2 == 0 ? 1 : 2); p <= max_value_; p += 2) {
        // p が素数かチェック
        if (table_[oddIndex(p)] == 0) {
            // p の最初の奇数合成数倍から開始
            // int64_t を使用: 2*p が INT_MAX を超える場合のオーバーフローを回避
            int64_t start = 2LL * p;
            if ((start & 1) == 0) {
                start += p;  // 奇数にする
            }
            if (start > max_value_) continue;  // 範囲外なら skip

            // p の倍数を 2*p ずつ進む（奇数のみ）
            for (int64_t m = start; m <= max_value_; m += 2LL * p) {
                table_[oddIndex(static_cast<int>(m))] = static_cast<uint32_t>(p);
            }
        }
    }
}

// 最大素因数を取得
int IntFactorTable::operator[](int n) const {
    // 範囲外
    if (n < 0 || n > max_value_) {
        return 0;
    }

    // 特殊ケース
    if (n == 0 || n == 1) {
        return 0;  // 素因数なし
    }

    // 2 は素数
    if (n == 2) {
        return 0;
    }

    // 偶数（2 以外）
    if ((n & 1) == 0) {
        // 偶数の場合、奇数部分の最大素因数と 2 を比較して大きい方を返す
        int odd_part = n;
        while ((odd_part & 1) == 0) {
            odd_part >>= 1;
        }
        if (odd_part == 1) {
            return 2;  // 2 の冪乗
        }
        // 奇数部分の最大素因数を取得
        int odd_lpf = table_[oddIndex(odd_part)];
        if (odd_lpf == 0) {
            // odd_part は素数
            // 奇数部分の素因数と 2 を比較
            return (odd_part > 2) ? odd_part : 2;
        }
        // table に記録されている最大素因数と 2 を比較
        return (odd_lpf > 2) ? odd_lpf : 2;
    }

    // 奇数の場合
    int lpf = table_[oddIndex(n)];
    if (lpf == 0) {
        // n は素数
        return 0;
    }
    return lpf;
}

// 素因数分解
std::vector<std::pair<Int, int>> IntFactorTable::factorize(const Int& n) const {
    std::vector<std::pair<Int, int>> factors;


    // 特殊状態のチェック
    if (n.isNaN() || n.isInfinite()) {
        return factors;  // 空
    }

    // 0 の場合
    if (n.isZero()) {
        factors.push_back({Int::Zero(), 1});
        return factors;
    }

    // 負の数の場合は絶対値を分解
    Int abs_n = abs(n);

    // 1 の場合
    if (abs_n.isOne()) {
        return factors;  // 空
    }

    Int remaining = abs_n;

    // 2 で割り切れるだけ割る（LSBの0ビット数 = 2の指数）
    size_t exp_2 = remaining.countTrailingZeros();
    if (exp_2 > 0) {
        factors.push_back({Int(2), static_cast<int>(exp_2)});
        remaining = remaining >> static_cast<int>(exp_2);
    }

    // remaining が 1 になったら終了
    if (remaining.isOne()) {
        return factors;
    }

    // テーブルを使って分解（remaining が max_value 以下の間）
    bool finished_in_table = false;
    while (remaining <= Int(max_value_) && !remaining.isOne()) {
        int n_int = remaining.toInt();
        int lpf = (*this)[n_int];

        if (lpf == 0) {
            // remaining は素数
            factors.push_back({remaining, 1});
            finished_in_table = true;
            break;
        }

        // lpf で割り切れるだけ割る
        Int p(lpf);
        int exp = 0;
        while (true) {
            Int quotient = remaining / p;
            Int remainder_div = remaining % p;

            if (!remainder_div.isZero()) {
                break;
            }

            remaining = quotient;
            exp++;
        }

        if (exp > 0) {
            factors.push_back({p, exp});
        }
    }

    // remaining > max_value の場合は 1リム・バッチ GCD で小素数因子を抽出
    if (!remaining.isOne() && !finished_in_table) {
        IntPrime::extractSmallPrimeFactors(remaining, factors);

        // 小素数で割り切った後も remaining > 1 なら、素数または大きな合成数
        if (!remaining.isOne()) {
            factors.push_back({remaining, 1});
        }
    }

    // 素因数を昇順にソート
    std::sort(factors.begin(), factors.end(),
              [](const std::pair<Int, int>& a, const std::pair<Int, int>& b) {
                  return a.first < b.first;
              });

    return factors;
}

// ファイルに保存
bool IntFactorTable::saveToFile(const std::string& filename) const {
    std::ofstream ofs(filename, std::ios::binary);
    if (!ofs) {
        return false;
    }

    // max_value を書き込み
    ofs.write(reinterpret_cast<const char*>(&max_value_), sizeof(max_value_));

    // table のサイズを書き込み
    uint32_t table_size = static_cast<uint32_t>(table_.size());
    ofs.write(reinterpret_cast<const char*>(&table_size), sizeof(table_size));

    // table の内容を書き込み
    ofs.write(reinterpret_cast<const char*>(table_.data()), table_size * sizeof(uint32_t));

    return ofs.good();
}

// ファイルから読み込み
bool IntFactorTable::loadFromFile(const std::string& filename) {
    std::ifstream ifs(filename, std::ios::binary);
    if (!ifs) {
        return false;
    }

    // max_value を読み込み
    int file_max_value;
    ifs.read(reinterpret_cast<char*>(&file_max_value), sizeof(file_max_value));
    if (!ifs || file_max_value < max_value_) {
        return false;  // ファイルのテーブルが小さい
    }

    // table のサイズを読み込み
    uint32_t file_table_size;
    ifs.read(reinterpret_cast<char*>(&file_table_size), sizeof(file_table_size));
    if (!ifs || file_table_size < static_cast<uint32_t>(table_.size())) {
        return false;  // サイズ不足
    }

    // table の内容を読み込み
    ifs.read(reinterpret_cast<char*>(table_.data()), table_.size() * sizeof(uint32_t));
    if (!ifs) {
        return false;
    }

    return true;
}

} // namespace calx