XCPC wiki - 数列

斐波那契数列的 20 种求法

矩阵快速幂

题意 设$a_{n}=\begin{cases}a_{n-1}+a_{n-2},&n \geqslant 3\\1,&n=1\lor n=2\end{cases}$，求$a_n\bmod m$。

解 设$\mathbf{M}=\begin{bmatrix}0 &1\\ 1 & 1\end{bmatrix}$​，则

$$\begin{bmatrix}F_n & 0 \\ F_{n+1} & 0\end{bmatrix} =\mathbf{M}^{n-1}\begin{bmatrix}1 & 0 \\ 1 & 0 \end{bmatrix}$$

于是$F_n=\left(\mathbf{M}^{n-1}\right)_{11}+\left(\mathbf{M}^{n-1}\right)_{12}$​。

推广 设$a_n=pa_{n-1}+qa_{n-2}$，给定$p,~q,~a_1,~a_2,~n,~m$，求$a_n\bmod m$。（VJwin11I - 这是签到题？要不先去看看其他的？（广义斐波那契数列） ）

解 本题修改了系数和首项，若设$\mathbf{M}=\begin{bmatrix}0 &1\\ q & p\end{bmatrix}$​，类似地，有

$$\begin{bmatrix}F_n & 0 \\ F_{n+1} & 0\end{bmatrix}=\mathbf{M}^{n-1}\begin{bmatrix}a_1 & 0 \\ a_2 & 0 \end{bmatrix}$$

于是$F_n=a_1\cdot\left(\mathbf{M}^{n-1}\right)_{11}+a_2\cdot\left(\mathbf{M}^{n-1}\right)_{12}$。

#include <cstdio>
#include <cmath>
using ll = long long;
inline ll read() {
    ll S = 0, C = 0; char Y = getchar();
    while (Y > 57 || Y < 48) C |= Y == 45, Y = getchar();
    while (Y <= 57 && Y >= 48) S = (S << 3) + (S << 1) + Y - 48, Y = getchar();
    return C ? -S : S;
}

ll MOD;

struct matrix {
    ll a1, a2, b1, b2;
    matrix(ll a1, ll a2, ll b1, ll b2) : a1(a1), a2(a2), b1(b1), b2(b2) {}
    matrix operator*(const matrix& y) {
        matrix ans((a1 * y.a1 + a2 * y.b1) % MOD,
            (a1 * y.a2 + a2 * y.b2) % MOD,
            (b1 * y.a1 + b2 * y.b1) % MOD,
            (b1 * y.a2 + b2 * y.b2) % MOD);
        return ans;
    }
};

matrix qpow(matrix a, ll n) {
    matrix ans(1, 0, 0, 1);
    while (n) {
        if (n & 1) ans = ans * a;
        a = a * a;
        n >>= 1;
    }
    return ans;
}

int main() {
    ll p = read(), q = read(), a1 = read(), a2 = read(), n = read();
    MOD = read();

    matrix M(0, 1, q, p);
    matrix ans = qpow(M, n - 1);
    matrix M0(a1, 0, 0, a2);
    printf("%lld\n", (a1 * ans.a1 + a2 * ans.a2) % MOD);

    return 0;
}

根号分治

题意 设$a_n=pa_{n-1}+qa_{n-2}$，给定$p,~q,~a_1,~a_2,~n,~m$，求$a_n\bmod m$。

思路 由于$n$ 巨大，从头开始一一推算数列的每一项是不可能的；又由于$m$ 巨大，利用线性递推数列具有模周期性也是不可能的。我们将原本的递推公式，即数列中某项与前两项的关系，推导得出数列中某项和与其遥遥相隔的连续两项之间的关系（比如$a_{100}$ 与$a_1,~a_2$ 之间的关系）。

解 设$a_n=c_k \cdot a_{n-k} + d_k \cdot a_{n-k-1}$，易知$c_1=p，d_1=q$，且$\{c_n\},~\{d_n\}$ 的递推公式是

$$\begin{cases}c_n = pc_{n-1}+d_{n-1}\\d_n = qc_{n-1}\end{cases}$$

如此，我们由最初两项向后推算时每步的步长可以非常大，忽略中间的许多项，从而节省时间。

具体实现方法是：取定$N$，首先根据$\{c_n\},~\{d_n\}$ 的递推公式计算出$c_N,~d_N$。然后根据已知的$a_1,~a_2$，计算出$a_3$，以及

$$\begin{cases} a_{N+2}=c_N a_2 + d_N a_1 \\ a_{N+3}=c_N a_3 + d_N a_2 \end{cases}$$

然后计算出$a_{N+4}$，再根据$a_{N+2},~a_{N+3},~a_{N+4}$ 计算出$a_{2N+3},~a_{2N+4},~a_{2N+5}$，以此类推，这样就由起初的连续三项一下子推出$N+1$ 项后的连续三项。

当$n$ 巨大时，用这样的方法不断向前推接近$n$，到与$n$ 的距离不足$N$ 时再一项项推算过去，直到求出$a_n$。

时间复杂度为$\Theta(\frac{n}{N}+N)$，因此我们取$N=\sqrt{n}$ 能最大程度上减小时间，最终时间复杂度为$\Theta(\sqrt{n})$，可以在给定时间内完成。

#include <cstdio>
#include <cmath>
using ll = long long;
inline ll read() {
    ll S = 0, C = 0; char Y = getchar();
    while (Y > 57 || Y < 48) C |= Y == 45, Y = getchar();
    while (Y <= 57 && Y >= 48) S = (S << 3) + (S << 1) + Y - 48, Y = getchar();
    return C ? -S : S;
}

int main() {
    ll p = read(), q = read(), a1 = read(), a2 = read(), a3, n = read(), MOD = read();
    p %= MOD, q %= MOD;
    a1 %= MOD, a2 %= MOD, a3 = (p * a2 + q * a1) % MOD;

    if (n == 1) return 0 * printf("%lld\n", a1);
    if (n == 2) return 0 * printf("%lld\n", a2);
    if (n == 3) return 0 * printf("%lld\n", a3);

    ll N = (ll)sqrt(n) + 7;
    if (n > N) {
        ll c = p, d = q;
        for (int i = 2; i <= N; ++i) {
            ll c1 = (p * c + d) % MOD, d1 = (q * c) % MOD;
            c = c1, d = d1;
        }
        while (n > N) {
            a1 = (c * a2 + d * a1) % MOD;
            a2 = (c * a3 + d * a2) % MOD;
            a3 = (p * a2 + q * a1) % MOD;
            n -= N + 1;
        }
    }
    for (int i = 4; i <= n; ++i) {
        a1 = a2, a2 = a3;
        a3 = (p * a2 + q * a1) % MOD;
    }
    printf("%lld\n", a3);

    return 0;
}

等差数列

CF928E - Vlad and an Odd Ordering

题意 将小于等于$n$ 的正整数依以下规则排成一列：最前面是从小到大的所有奇数，接着是将这些奇数的$2$ 倍中未使用的数从小到大排在后面，然后再将这些奇数的$3$ 倍中未使用的数从小到大排在后面数，直到所有数都排完，求这个数列中第$k$ 个数。

思路 注意到只会放置奇数的$2$ 的幂次倍，且每次放置的数构成等差数列，因此只需知道第$k$ 个数所在的等差数列即可。考虑每个等差数列的首项$a$，公差$d$，依

$$num=\Big\lfloor\frac{n-a}{d}\Big\rfloor+1$$

计算项数$num$​，如果$k>num$​，令$k:=k-num$​，考虑下一个等差数列，最终依

$$a_k=a+(k-1)d$$

计算答案。

int a = 1, d = a << 1, num = (n - a) / d + 1;
while (k > num) {
    k -= num;
    a <<= 1;
    d <<= 1;
    num = (n - a) / d + 1;
}
printf("%d\n", a + (k - 1) * d);

或

int n = read(), k = read();
int pos = (n + 1) >> 1, pow = 2, ans = 2 * k - 1;
while (ans > n) {
    ans = (pow << 1) * (k - pos) - pow;
    pos += (n - pow) / (pow << 1) + 1;
    pow <<= 1;
}
printf("%d\n", ans);

2024 杭电多校 7008 - 循环图

[[…/contests/2024杭电多校/2024-08-05：2024“钉耙编程”中国大学生算法设计超级联赛（7）|2024-08-05：2024“钉耙编程”中国大学生算法设计超级联赛（7）]]

坎格鲁斯普雷被困在了一张循环图里，这张循环图有无数个节点，初始时坎格鲁斯普雷在$1$ 号节点。

循环图是边存在着一定循环关系的图，循环图里面的边可以用循环周期$n$ 和$m$ 对三元组$(u_i,v_i,w_i)$$(1 \leq u_i \leq n , u_i + 1 \leq v_i \leq 2 \times n,1 \leq w_i \leq 10^9)$ 表示。每对三元组$(u_i,v_i,w_i)$ 表示，对于循环图内所有的满足$s = u_i + k \times n$ ，$t = v_i + k \times n$$( k \in N )$ 的点对$(s,t)$ ，都存在有$w_i$ 条从点$s$ 通往点$t$ 的边。

现在，坎格鲁斯普雷知道了这张循环图的第$L$ 个节点到第$R$ 个节点各存在着一个出口，坎格鲁斯普雷需要到达这些节点中的任意一个才能逃出循环图（到达有出口存在的节点后不一定要立刻逃出）。坎格鲁斯普雷想请你帮他算算，他有多少种逃出这张循环图的方式。由于答案可能很大，你需要输出答案对$10^9+7$ 取模后的结果。

#include <iostream>
#include <algorithm>
#include <vector>
using ll = long long;
using Z = ModInt<1000000007>;
using namespace std;

struct Matrix {
    int n;
    vector<vector<Z>> a;

    Matrix(int _n) : n(_n) {
        a.resize(n, vector<Z>(n, 0));
    }

    void unit() {
        for (int i = 0; i < n; ++i) {
            a[i][i] = 1;
        }
    }

    Matrix operator + (Matrix b) {
        Matrix res(n);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                res.a[i][j] = a[i][j] + b.a[i][j];
            }
        }
        return res;
    }

    Matrix operator * (Matrix b) {
        Matrix res(n);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                for (int k = 0; k < n; ++k) {
                    res.a[i][j] += a[i][k] * b.a[k][j];
                }
            }
        }
        return res;
    }

    Matrix operator ^ (ll d) {
        Matrix M(n);
        M.unit();
        for (int bit = 60; bit >= 0; --bit) {
            M = M * M;
            if (d & (1ll << bit)) {
                M = M * *this;
            }
        }
        return M;
    }

    Matrix operator >> (ll d) {
        Matrix M(n), res(n);
        M.unit();
        for (int bit = 60; bit >= 0; --bit) {
            res = res * M + res;
            M = M * M;
            if (d & (1ll << bit)) {
                M = M * *this;
                res = res + M;
            }
        }
        return res;
    }
};

void eachT() {
    int n, m;
    cin >> n >> m;

    ll L, R;
    cin >> L >> R;
    L--, R--;    // 0-index

    Matrix E(n), M(n);
    while (m--) {
        int x, y, w;
        cin >> x >> y >> w;
        x--, y--;
        if (y < n) E.a[x][y] = w;
        else M.a[x][y - n] = w;
    }

    // 初始值
    vector<Z> f(n);
    f[0] = 1;
    for (int x = 0; x < n; ++x) {
        for (int y = 0; y < n; ++y) {
            f[y] += f[x] * E.a[x][y];
        }
    }

    // 广义 Floyd
    for (int k = 0; k < n; ++k) {
        for (int x = 0; x < n; ++x) {
            for (int y = 0; y < n; ++y) {
                M.a[x][y] += M.a[x][k] * E.a[k][y];
            }
        }
    }

    ll Lq = L / n, Lr = L % n;
    ll Rq = R / n, Rr = R % n;

    Z ans = 0;
    if (Lq < Rq) {
        Matrix MM = M ^ Lq;
        vector<Z> g(n);
        for (int x = 0; x < n; ++x) {
            for (int y = 0; y < n; ++y) {
                g[y] += f[x] * MM.a[x][y];
            }
        }

        // from Lq*(n-1)+Lr to Lq*n
        for (int y = Lr; y < n; ++y) {
            ans += g[y];
        }

        // from Lq*n to Rq*n
        MM = M >> (Rq - Lq - 1);
        for (int x = 0; x < n; ++x) {
            for (int y = 0; y < n; ++y) {
                ans += g[x] * MM.a[x][y];
            }
        }

        // from Rq*n to Rq*n+Rr
        MM = M ^ Rq;
        for (int x = 0; x < n; ++x) {
            for (int y = 0; y <= Rr; ++y) {
                ans += f[x] * MM.a[x][y];
            }
        }
    }
    else {
        Matrix MM = M ^ Lq;
        for (int x = 0; x < n; ++x) {
            for (int y = Lr; y <= Rr; ++y) {
                ans += f[x] * MM.a[x][y];
            }
        }
    }

    cout << ans << '\n';
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout.tie(nullptr);

    int T;
    cin >> T;

    while (T--) {
        eachT();
    }

    return 0;
}