bzoj 3328: PYXFIB

题意：

给定$n, k, p$，求

$\sum_{i = 0}^{\lfloor \frac{n}{k} \rfloor} \binom{n}{ik} F_{ik} \mod p$

其中$F$表示斐波那契数列。

题解：

$\begin{aligned} & \sum_{i = 0}^{\lfloor \frac{n}{k} \rfloor} \binom{n}{ik} F_{ik} \\\ =& \sum_{i = 0}^n \binom{n}{i} F_i [k | i] \\\ =& \sum_{i = 0}^n \binom{n}{i} (A^i)_{0, 0} \frac{1}{k} \sum_{j = 0}^{k - 1} \omega_k^{ij} \\\ =& \frac{1}{k} \sum_{j = 0}^{k - 1} \sum_{i = 0}^n \binom{n}{i} (A^i)_{0, 0} \omega_k^{ij} \\\ =& \frac{1}{k} \sum_{j = 0}^{k - 1} ((A \omega_k^j + I)^n)_{0, 0} \end{aligned}$

其中$A$是斐波那契数列的转移矩阵。

代码：

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

int T, mod, K, ans;
long long n;
struct matrix
{
    int x, y;
    long long a[2][2];
    matrix() : a() {}
    long long *operator[](int x) { return a[x]; }
} I, a;

matrix operator*(matrix x, matrix y)
{
    matrix ans;
    ans.x = x.x;
    ans.y = y.y;
    for (int i = 0; i < x.x; i++)
        for (int j = 0; j < y.y; j++)
            for (int k = 0; k < y.x; k++)
                (ans[i][j] += x[i][k] * y[k][j] % mod) %= mod;
    return ans;
}
matrix operator*(matrix x, int y)
{
    for (int i = 0; i < x.x; i++)
        for (int j = 0; j < x.y; j++)
            x[i][j] = x[i][j] * y % mod;
    return x;
}
matrix operator+(matrix x, matrix y)
{
    for (int i = 0; i < x.x; i++)
        for (int j = 0; j < x.y; j++)
            (x[i][j] += y[i][j]) %= mod;
    return x;
}
matrix Pow(matrix x, long long y)
{
    matrix ans = I;
    for (; y; y >>= 1, x = x * x)
        if (y & 1) ans = ans * x;
    return ans;
}
int Pow(int x, int y)
{
    int ans = 1;
    for (; y; y >>= 1, x = 1LL * x * x % mod)
        if (y & 1) ans = 1LL * ans * x % mod;
    return ans;
}
int get_rt()
{
    int hh = mod - 1, num = 0;
    static int p[20];
    for (int i = 2; i * i <= hh; i++)
        if (hh % i == 0)
        {
            p[num++] = i;
            while (hh % i == 0) hh /= i;
        }
    if (hh > 1) p[num++] = hh;
    int g = 2;
    while (1)
    {
        bool tf = false;
        for (int i = 0; i < num; i++)
            if (Pow(g, (mod - 1) / p[i]) == 1)
            {
                tf = true;
                break;
            }
        if (!tf) return g;
        g++;
    }
}
int main()
{
    scanf("%d", &T);
    I.x = I.y = a.x = a.y = 2;
    I[0][0] = I[1][1] = a[0][0] = a[0][1] = a[1][0] = 1;
    while (T--)
    {
        scanf("%lld%d%d", &n, &K, &mod);
        int g = Pow(get_rt(), (mod - 1) / K), now = 1;
        // printf("%d\n", g);
        ans = 0;
        for (int i = 0; i < K; i++, now = 1LL * now * g % mod)
            (ans += Pow(a * now + I, n).a[0][0]) %= mod;
        printf("%lld\n", 1LL * ans * Pow(K, mod - 2) % mod);
    }
}