【BZOJ3328】PYXFIB（数学）

最新推荐文章于 2020-01-25 21:22:31 发布

weixin_30394981

最新推荐文章于 2020-01-25 21:22:31 发布

阅读量75

点赞数

原文链接：http://www.cnblogs.com/zyt1253679098/p/10184001.html

版权

什么都不会的数学蒻菜瑟瑟发抖……Orz橙子（和兔子）

题目：

BZOJ3328

分析：

橙子给我安利的数学题……（然后我就看着他因为矩阵乘法多模了一次卡了一天常数qwq表示同情）

先考虑一个子问题：求下面这个式子的值：

\[\sum_{i=0}^n C_n^i F_i\]

首先我们知道\(F_i=\left[\begin{matrix}1 & 1\\1 & 0\end{matrix}\right]^i\)。设那个矩阵为\(A\)，即\(F_i=A^i\)。（注意这题斐波那契数列下标从\(0\)开始，所以\(F_2=2\)。）

（不知道？你把\(\left[\begin{matrix}F_i & F_{i-1}\\F_{i-2} & 0\end{matrix}\right]\)乘一下\(A\)试试。一开始左下方的值并不影响计算结果。）

然后\(\sum\limits_{i=0}^n C_n^i A^i\)这个东西好像长得挺像二项式定理的。事实上，如果两个矩阵相乘可以交换（即\(AB=BA\)），则这两个矩阵之和的幂可以用二项式定理展开。因此，上面那个式子就是\((A+I)^n\)，其中\(I\)是单位矩阵（任何矩阵与单位矩阵相乘都是可交换的）。

~~（于是我推到这里忘了原题是什么兴高采烈地去给橙子说我会了~~

那么问题来了，如何“筛”掉当\(i\)不是\(k\)的倍数的时候对答案的贡献呢？介绍一个神奇的东西：单位根。

首先先介绍原根。若在模\(p\)（这里只讨论\(p\)是质数的情况）意义下，\(a(0\leq a<p)\)是一个满足对于任意\(i(1\leq i < p - 1)\)都满足\(a^i\neq 1\)的整数，则\(a\)是\(p\)的原根。（根据欧拉定理，\(a^i=a^{i \mod \phi(p)}\)，所以不讨论指数不小于\(\phi(p)\)即\(p-1\)的情况。）原根满足对于任意小于\(\phi(p)\)的非负整数\(i\)，\(a^i\)在模\(p\)意义下互不相同。以下的“原根”都特指最小原根，记为\(g\)。求法？从\(2\)开始暴力枚举\(a\)，把\(i=p'\)（其中\(p'\)是\(p-1\)的质因数）挨个试一遍看有没有\(a^i\equiv 1\mod p\)的情况就完了。（并不知道为什么是对的并且跑得飞快）

模\(p\)的\(k\)次单位根\(w_k=g^{\frac{p-1}{k}}\)，所以\(w_k^k \equiv 1\mod p\)。

在模\(p\)的意义下，\(p\)的\(k\)次单位根\(w_k\)具有如下性质：

\[\sum_{i=0}^{k-1}w_k^{ij}=\begin{cases}k\ (j|k)\\0\ (otherwise)\end{cases}\]

证明？当\(j|k\)，\(ij\)全是\(k\)的倍数，所以\(w_k^{ij}=1\)，全部加起来就是\(k\)。

否则，大力等比数列求和公式\(w_k^0\cdot \frac{1-w_k^{jk}}{1-w_k^j}\)。上面\(w_k^{jk}=1\)，所以和是\(0\)。这个有点像FFT中的单位根，只是一个在单位圆上，一个在模域上，都是循环的（注意模域中指数的循环节是\(\phi(p)=p-1\)）。

基于这个性质，我们达到了筛掉\(i\)不是\(k\)的倍数的情况的目的。那么要求的那个式子就变成了：

\[\sum_{i=0}^n C_n^i A^i\cdot\frac{1}{k}\sum_{j=0}^{k-1}w_k^{ij}\]

发现后面那一堆和\(i\)没什么关系，提到前面来（相当于把\(w_k^{ij}\)分配进去）：

\[\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n C_n^i A^iw_k^{ij}\]

好像挺像二项式定理的，只是\(w_k^j\)的指数比较丑。但我们可以这样……

\[\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n C_n^i A^iw_k^{(-i)\cdot(-j)}\]

然后设\(x_j=w_k^{-j}\)，得到：

\[\frac{1}{k}\sum_{j=0}^{k-1}x_j^{-n}\sum_{i=0}^n C_n^i A^ix_j^{n-i}\]

于是可以上矩阵二项式定理了（注意因为\(x_j\)是一个整数，所以要乘上单位矩阵\(I\)才能与\(A\)相加）。

\[\frac{1}{k}\sum_{j=0}^{k-1}x_j^{-n}(x_jI+A)^n\]

照着这个就可以直接算了qwq

代码：

矩阵乘法的时候一定不要模两次！一定不要模两次！一定不要模两次！否则你会像橙子一样卡一天常数。

以及记着用上面提到过的欧拉定理把所有指数搞成非负的。

#include <cstdio>
#include <cstring>
#include <cctype>
#include <algorithm>
#define _ 0
using namespace std;

namespace zyt
{
    template<typename T>
    inline bool read(T &x)
    {
        char c;
        bool f = false;
        x = 0;
        do
            c = getchar();
        while (c != EOF && c != '-' && !isdigit(c));
        if (c == EOF)
            return false;
        if (c == '-')
            f = true, c = getchar();
        do
            x = x * 10 + c - '0', c = getchar();
        while(isdigit(c));
        if (f)
            x = -x;
        return true;
    }
    template<typename T>
    inline void write(T x)
    {
        static char buf[20];
        char *pos = buf;
        if (x < 0)
            putchar('-'), x = -x;
        do
            *pos++ = x % 10 + '0';
        while (x /= 10);
        while (pos > buf)
            putchar(*--pos);
    }
    typedef long long ll;
    ll n;
    int k, p;
    struct Matrix
    {
        int n, m, data[2][2];
        Matrix(const int _n = 0, const int _m = 0)
            : n(_n), m(_m)
        {
            for (int i = 0; i < n; i++)
                memset(data[i], 0, sizeof(int[m])); 
        }
        ~Matrix() {}
        Matrix operator + (const Matrix &b) const
        {
            Matrix ans(n, m);
            for (int i = 0; i < n; i++)
                for (int j = 0; j < m; j++)
                    ans[i][j] = (data[i][j] + b[i][j]) % p;
            return ans;
        }
        Matrix operator * (const Matrix &b) const
        {
            Matrix ans(n, b.m);
            for (int i = 0; i < n; i++)
                for (int k = 0; k < m; k++)
                    for (int j = 0; j < b.m; j++)
                        ans[i][j] = (ans[i][j] + 
                            (ll)data[i][k] * b[k][j] % p) % p;
            return ans;
        }
        Matrix operator * (const ll &x) const
        {
            Matrix ans(n, m);
            for (int i = 0; i < n; i++)
                for (int j = 0; j < m; j++)
                    ans[i][j] = ((ll)data[i][j] * x) % p;   
            return ans;
        }
        const int *operator [] (const int a) const
        {
            return data[a]; 
        }
        int *operator [] (const int a)
        {
            return data[a]; 
        }
    };
    inline Matrix get_identity(const int n)
    {
        Matrix ans(n, n);
        for (int i = 0; i < n; i++)
            ans[i][i] = 1;
        return ans; 
    }
    inline Matrix power(Matrix a, ll b)
    {
        Matrix ans = get_identity(a.n);
        while (b)
        {
            if (b & 1)
                ans = ans * a;
            a = a * a;
            b >>= 1;
        }
        return ans;
    }
    inline int power(int a, ll b)
    {
        int ans = 1;
        while (b)
        {
            if (b & 1)
                ans = (ll)ans * a % p;
            a = (ll)a * a % p;
            b >>= 1;
        }
        return ans;
    }
    inline int inv(const int a)
    {
        return power(a, p - 2); 
    }
    pair<int, int> prime[20];
    int cnt;
    inline void get_prime(int n)
    {
        cnt = 0;
        for (int i = 2; i * i <= n; i++)
        {
            if (n % i == 0)
                prime[cnt++] = make_pair(i, 0);
            while (n % i == 0)
                ++prime[cnt - 1].second, n /= i;
        }
        if (n != 1)
            prime[cnt++] = make_pair(n, 1);
    }
    inline int get_g(const int p)
    {
        get_prime(p - 1);
        for (int i = 2; i < p; i++)
        {
            bool flag = true;
            for (int j = 0; j < cnt && flag; j++)
                if (power(i, (p - 1) / prime[j].first) == 1)
                    flag = false;
            if (flag)
                return i;
        }
    }
    int work()
    {
        int T;
        read(T);
        Matrix I = get_identity(2);
        Matrix A(2, 2);
        A[0][0] = A[0][1] = A[1][0] = 1;
        while (T--)
        {
            read(n), read(k), read(p);
            int omega = power(get_g(p), (p - 1) / k), ans = 0;
            int tmp = power(omega, p - k - 1);
            for (int i = 0; i > -k; i--)
            {
                tmp = (ll)tmp * omega % p;
                ans = (ans + (ll)power(tmp, -n) * power(I * tmp + A, n)[0][0] % p) % p;
            }
            write((ll)ans * inv(k) % p);
            putchar('\n');
        }
        return (0^_^0); 
    }
}
int main()
{
    freopen("3328.in", "r", stdin);
    freopen("3328.out", "w", stdout);
    return zyt::work(); 
}

转载于:https://www.cnblogs.com/zyt1253679098/p/10184001.html

weixin_30394981

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
【BZOJ3328】PYXFIB（数学）

什么都不会的数学蒻菜瑟瑟发抖……Orz橙子（和兔子）题目：BZOJ3328分析：橙子给我安利的数学题……（然后我就看着他因为矩阵乘法多模了一次卡了一天常数qwq表示同情）先考虑一个子问题：求下面这个式子的值：\[\sum_{i=0}^n C_n^i F_i\]首先我们知道\(F_i=\left[\begin{matrix}1 & 1\\1 & 0\end{matrix...
复制链接

扫一扫