【题解】Codeforces 193 E Fibonacci Number

题目链接

Codeforces 193 E Fibonacci Number

题目大意

输入一个数 \(f\) ，求最小的 \(x\) 使得 \(F_x \equiv f \ (mod \ 10^{13})\)，其中 \(F = {0, 1, 1, 2, 3 ...}\) 为斐波那契数列，下标从0开始。

思路

和CCPC 2025 济南站的 A 题几乎一模一样，这题似乎还难一点，可惜我没有在比赛前刷到这道题。

考虑大步小步(BSGS)算法，不难猜到对于固定模数，斐波那契会出现循环，需要前置知识矩阵乘法、快速幂和快速幂求解斐波那契数列。

记转移矩阵为

\[A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \]

构造第 \(1\) 项斐波那契数列列向量为

\[F_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \]

第 \(i\) 项斐波那契数列列向量为

\[F_i = A^{i-1}F_1 \]

其中，第一行为斐波那契数列的第 \(i\) 项，第二行为斐波那契数列的第 \(i-1\) 项，可以用矩阵乘法的快速幂 \(O(log \ n)\) 求一项斐波那契数列。

下文的所有运算符，包括 \(mod\) 如果发生在 矩阵/向量 与 数值之间，默认为元素操作。

记 \(F_x^k = F_x \ mod \ 10^{k}\) 从小到大枚举 \(x \in [0, +\infty)\) 显然，\(F_x^k\) 会进入一个长度不超过 \(10^{2k}\) 的循环，记循环长度为 \(T\) ，即 \(F_x^k = F_{nT+x}^k\)

例如，\(f = 998244353, \ k=1\) 时，这个循环的长度为 60（横向优先阅读，列向量转置为行向量）:

(1, 1)(2, 1)(3, 2)(5, 3)(8, 5)
(3, 8)(1, 3)(4, 1)(5, 4)(9, 5)
(4, 9)(3, 4)(7, 3)(0, 7)(7, 0)
(7, 7)(4, 7)(1, 4)(5, 1)(6, 5)
(1, 6)(7, 1)(8, 7)(5, 8)(3, 5)
(8, 3)(1, 8)(9, 1)(0, 9)(9, 0)
(9, 9)(8, 9)(7, 8)(5, 7)(2, 5)
(7, 2)(9, 7)(6, 9)(5, 6)(1, 5)
(6, 1)(7, 6)(3, 7)(0, 3)(3, 0)
(3, 3)(6, 3)(9, 6)(5, 9)(4, 5)
(9, 4)(3, 9)(2, 3)(5, 2)(7, 5)
(2, 7)(9, 2)(1, 9)(0, 1)(1, 0)

此时考虑从低往高第递推考虑 \(k\) 位，记此时的周期 \(T\) 为 \(T_k\)，从小到大枚举 \(n\) 容易发现，对于任意 \(x\)， \(F_{nT_k+x}^k\) 会进入长度不超过 \(10^2\) （记长度为 \(t_k\)） 的循环。而且会有若干个 \(x\) 满足 “后k位匹配”，令它们在集合 \(S_k\) 中，即

\[S_k = \{ F_x^k | F_x^k = \begin{bmatrix} f \\ * \end{bmatrix} \ mod 10^k , x < T\} \]

其中 \(*\) 为任意值，对于刚刚 \(k=1\) 的情况,也就是最后一位匹配上的情况，\(S_k\)：

(3, 2)(3, 8)(3, 4)(3, 5)(3, 7)(3, 0)(3, 3)

注意，由于斐波那契数列与两项有关，故上述虽然最后一位匹配上了即满足条件，仍然不能因为“通配”而合并，需要用 vector 存下来。

选定一个 \(F_{x_k} \in S_k\) 若存在一个最小的 \(n=n'\) 使得

\[\begin{bmatrix} f \\ * \end{bmatrix} \equiv F_{n'T_k+x_k} \ (mod \ 10^{k+1}) \]

则可以发现 \(F_{n'T_k+x_k} \in S_{k+1}\)，而且新的 \(x_{k+1} = n'T_k+x_k\)，我们可以通过从小到大枚举 \(n\) 递推求解 \(S_{k+1}\)，例如从 \(k=1\) 转移到 \(k=2\):

n=0
  Null
n=1
  (13, 8) -> (53, 88)
  (33, 37) -> (53, 77)
n=2
  (33, 44) -> (53, 84)
  (93, 25) -> (53, 45)
  (73, 3) -> (53, 63)
n=3
  (73, 99) -> (53, 59)
n=4
  Null

注意，这里发现(13, 8)虽然以3结尾，而且最小，但是很多状态都转移不到，这就是为什么之前的通配符不能合并。

如果出现了某个 \(S\) 为空集，说明答案不存在。

重复向前递推，直到模数达到 \(10^{13}\) ，就找到了答案。

AC代码

注意：记得开 __int128，记得特判0和1。

#include <iostream>
#include <vector>
using namespace std;
typedef long long int64;
const int N = 2;
const int LOGM = 14;
const __int128 MODER = 100000000000000;

class Matrix {
public:
    __int128 arr[N][N];

    Matrix operator*(const Matrix &o) const {
        Matrix res = {0};

        for (int64 i = 0; i < N; i++) {
            for (int64 j = 0; j < N; j++) {
                for (int64 k = 0; k < N; k++) {
                    res.arr[i][j] += arr[i][k] * o.arr[k][j];
                }
                res.arr[i][j] %= MODER;
            }
        }

        return res;
    }

    Matrix operator%(const int64 m) {
        Matrix res = {0};

        for (int64 i = 0; i < N; i++) {
            for (int64 j = 0; j < N; j++) {
                res.arr[i][j] = arr[i][j] % m;
            }
        }

        return res;
    }

    bool operator==(const Matrix &o) {
        for (int64 i = 0; i < N; i++) {
            for (int64 j = 0; j < N; j++) {
                if (arr[i][j] != o.arr[i][j]) {
                    return false;
                }
            }
        }
        return true;
    }

    bool operator!=(const Matrix &o) {
        return !((*this) == o);
    }
};

Matrix E = {
        1, 0,
        0, 1
};

Matrix quickPow(Matrix a, int64 x) {
    if (x == 0) {
        return E;
    } else if (x == 1) {
        return a;
    } else {
        Matrix t = quickPow(a, x >> 1);
        Matrix ans = t * t;
        if (x & 1) {
            ans = ans * a;
        }
        return ans;
    }
}

const Matrix A = {1, 1, 1, 0};
const Matrix st{1, 0, 0, 0};
vector<int64> ans[LOGM];

int main() {
    int64 f;
    cin >> f;
    if(f == 0) {
        cout<<0<<endl;
        return 0;
    }
    if(f == 1) {
        cout<<1<<endl;
        return 0;
    }
    int64 moder = 10;
    int64 prod = 1;
    ans[0].push_back(0);
    for (int i = 1; i < LOGM; i++, moder *= 10) {
        int j = 0;
        bool flag = true;
        do {
            for(int64 y: ans[i-1]) {
                Matrix Y = quickPow(A, prod * j + y) * st% (moder*10);
                if (Y.arr[0][0] % (moder) == f % (moder)) {
                    ans[i].push_back(j * prod + y);
                    flag = false;
                }
            }
        } while (quickPow(A, 1 + prod * (++j)) * st % moder != A * st);
        prod *= j;
        if(flag) {
            cout<<-1<<endl;
            return 0;
        }
    }
    cout << ans[LOGM-1][0] + 1<< endl;
    return 0;
}

时间复杂度我不会算……

第一层: \(ans_1\) 大小不超过 10，周期 \(t_1\) 不超过100。
第 \(i\) 层由第 \(i-1\) 层转移而来，周期 \(t_i\) 不超过 \(t_{i-1} \cdot 10\) ， while 循环执行不超过10次，由于这一层的所有答案后缀相同，有且仅有最高位不同且必须不同，故答案集合不会变大 ，即\(|ans_i| \leq |ans_{i-1}|\)。
故最终 \(ans\) 集合大小不超过10。
故时间复杂度为 \(O(b \cdot a^2 \ \cdot log b^a) = O(a^3b \ log \ b)\)，其中\(a^b\)为模数，在本题中 \(a=10,\ b=13\)。