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

Description

求

∑i=0⌊nk⌋(nik)Fik

1≤n≤1018,1≤k≤2×104,p为质数且pmodk=1

Solution

设g为p的原根，ω=gp−1k。

F=(1110)

∑i=0⌊nk⌋(nik)Fik====∑i=1n(ni)Fi[k∣i]1k∑i=0n(ni)Fi∑j=0k−1ωij1k∑j=0k−1ωjn∑i=0n(ni)Fi(ω−j)n−i1k∑j=0k−1ωjn(F+ω−jI)n

剩下的东西矩阵快速幂搞一搞就好啦。

#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;}