[数论][二项式定理][矩阵乘法] BZOJ 3328: PYXFIB

$Description$

求

$$\sum_{i=0}^{\lfloor{n\over k}\rfloor}{{n\choose {ik}}F_{ik}}$$ $1\le n\le10^{18},1\le k\le2\times10^4,p$为质数且 $p\bmod k = 1$

$Solution$

设$g$为$p$的原根，$\omega=g^{{p-1}\over k}$。

$$F=\left( \begin{array}{ccc} 1 & 1 \\ 1 & 0 \end{array} \right)$$
$$\begin{eqnarray} \sum_{i=0}^{\lfloor{n\over k}\rfloor}{{n\choose {ik}}F^{ik}}&=&\sum_{i=1}^n{n\choose i}F^i[k\mid i]\\ &=&{1\over k}\sum_{i=0}^n{n\choose i}F^i\sum_{j=0}^{k-1}\omega^{ij}\\ &=&{1\over k}\sum_{j=0}^{k-1}\omega^{jn}\sum_{i=0}^n{n\choose i}F^i(\omega^{-j})^{n-i}\\ &=&{1\over k}\sum_{j=0}^{k-1}\omega^{jn}(F+\omega^{-j}I)^n \end{eqnarray}$$
剩下的东西矩阵快速幂搞一搞就好啦。

#include <bits/stdc++.h>
using namespace std;

const int N = 202020;
typedef long long ll;

int test;
ll n;
int k, P, g, w, Iw, Ppcnt, Ik;
struct Matrix {
    int a[3][3];
    Matrix(void) {
    }
    inline int *operator [](int x) {
        return a[x];
    }
    inline friend Matrix operator *(Matrix a, Matrix b) {
        Matrix res;
        for (int i = 1; i <= 2; i++)
            for (int j = 1; j <= 2; j++) {
                res[i][j] = 0;
                for (int k = 1; k <= 2; k++) {
                    res[i][j] += (ll)a[i][k] * b[k][j] % P;
                    while (res[i][j] >= P) res[i][j] -= P;
                }
            }
        return res;
    }
    inline friend Matrix operator *(int a, Matrix b) {
        for (int i = 1; i <= 2; i++)
            for (int j = 1; j <= 2; j++)
                b[i][j] = (ll)b[i][j] * a % P;
        return b;
    }
    inline friend Matrix operator *(Matrix a, int b) {
        return b * a;
    }
    inline friend Matrix operator +(Matrix a, Matrix b) {
        for (int i = 1; i <= 2; i++)
            for (int j = 1; j <= 2; j++) {
                a[i][j] += b[i][j];
                while (a[i][j] >= P) a[i][j] -= P;
            }
        return a;
    }
};
Matrix I, A, Ans, O;
int pp[N];

inline void Add(int &x, int a) {
    x += a; while (x >= P) x -= P;
}
inline int Pow(int a, ll b) {
    int c = 1;
    while (b) {
        if (b & 1) c = (ll)c * a % P;
        b >>= 1; a = (ll)a * a % P;
    }
    return c;
}
inline int Inv(int x) {
    return Pow(x, P - 2);
}
inline Matrix Pow(Matrix a, ll b) {
    Matrix c = I;
    while (b) {
        if (b & 1) c = c * a;
        b >>= 1; a = a * a;
    }
    return c;
}
inline int PRoot(int P) {
    if (P == 3) return 2;
    int phi = P - 1;
    double m = sqrt(phi);
    Ppcnt = 0;
    for (int i = 2; i <= m; i++)
        if (phi % i == 0) {
            pp[++Ppcnt] = i;
            if (phi / i > i) pp[++Ppcnt] = phi / i;
        }
    for (int g = 2; ; g++) {
        for (int i = 1; i <= Ppcnt; i++)
            if (Pow(g, pp[i]) == 1) break;
            else if (i == Ppcnt) return g;
    }
}
inline Matrix F(int x) {
    return Pow(Inv(x), n) * Pow(x * I + A, n);
}

int main(void) {
    freopen("1.in", "r", stdin);
    scanf("%d", &test);
    I[1][1] = I[2][2] = 1;
    A[1][1] = A[1][2] = A[2][1] = 1;
    while (test--) {
        scanf("%lld%d%d", &n, &k, &P);
        Ans = O; g = PRoot(P);
        w = Pow(g, (P - 1) / k); Iw = Inv(w);
        for (int i = 0; i < k; i++)
            Ans = Ans + F(Pow(Iw, i));
        Ans = Ans * Inv(k);
        printf("%d\n", Ans[1][1]);
    }
    return 0;
}