BZOJ 3328: PYXFIB 解题报告

BZOJ 3328: PYXFIB

题意

给定\(n,p,k(1\le n\le 10^{18},1\le k\le 20000,1\le p\le 10^9,p \ is \ prime,k|(p-1))\),求
\[ \sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}\binom{n}{ik}F_{ik} \]
其中\(F_{i}\)代表斐波那契数列的第\(i\)项。


首先分析一波题目条件,有个很奇怪的条件\(k|\varphi(p)\)

如果对NTT或者单位根有所了解,我们可以猜到这个条件可以取出单位根来。
\[ w_k=g^{\frac{p-1}{k}} \]
然而如果没发现,推一波式子之后也可以发现需要使用单位根。

不过可以先考虑如何处理fib,一个简单的想法是构造矩阵乘法,比如
\[ A=\begin{bmatrix}1&1\\1& 0\end{bmatrix} \]
可以得到
\[ A^i=\begin{bmatrix}F_i&F_{i-1}\\\dots&\dots\end{bmatrix} \]
然后我们直接把式子转换成求矩阵
\[ \begin{aligned} &\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}\binom{n}{ik}A^{ik}\\ =&\sum_{i=0}^n\binom{n}{i}A^i[k|i]\\ =&\sum_{i=0}^n\binom{n}{i}A^i\frac{1}{k}\sum_{j=0}^{k-1}w_k^{ij}\\ =&\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^n\binom{n}{i}A^i(w_k^j)^i\\ =&\frac{1}{k}\sum_{j=0}^{k-1}(w_k^jA+I)^n \end{aligned} \]
复杂度\(O(Tk\log n)\)


Code:

#include <cstdio>
#include <cctype>
#define ll long long
const int SIZE=1<<21;
char ibuf[SIZE],*iS,*iT;
//#define gc() (iS==iT?(iT=(iS=ibuf)+fread(ibuf,1,SIZE,stdin),iS==iT?EOF:*iS++):*iS++)
#define gc() getchar()
template <class T>
void read(T &x)
{
    x=0;char c=gc();
    while(!isdigit(c)) c=gc();
    while(isdigit(c)) x=x*10+c-'0',c=gc();
}
ll n;
int k,p;
inline int add(int a,int b){return a+b>=p?a+b-p:a+b;}
#define mul(a,b) (1ll*(a)*(b)%p)
inline int qp(int d,int k)
{
    int f=1;
    while(k)
    {
        if(k&1) f=mul(f,d);
        d=mul(d,d);
        k>>=1;
    }
    return f;
}
struct Matrix
{
    int a,b,c,d;
    Matrix friend operator *(Matrix a,Matrix b)
    {
        Matrix c;
        c.a=add(mul(a.a,b.a),mul(a.b,b.c));
        c.b=add(mul(a.a,b.b),mul(a.b,b.d));
        c.c=add(mul(a.c,b.a),mul(a.d,b.c));
        c.d=add(mul(a.c,b.b),mul(a.d,b.d));
        return c;
    }
}A;
inline Matrix mqp(Matrix d,ll k)
{
    Matrix f=d;--k;
    while(k)
    {
        if(k&1) f=f*d;
        d=d*d;
        k>>=1;
    }
    return f;
}
int getg(int p)
{
    int s[30]={},phi=p-1;
    for(int i=2;i*i<=phi;i++)
    {
        if(phi%i==0)
        {
            s[++s[0]]=i;
            while(phi%i==0) phi/=i;
        }
    }
    if(phi!=1) s[++s[0]]=phi;
    for(int i=2;;i++)
    {
        int flag=1;
        for(int j=1;j<=s[0];j++)
            if(qp(i,(p-1)/s[j])==1)
            {
                flag=0;
                break;
            }
        if(flag) return i;
    }
}
int main()
{
    int T;read(T);
    while(T--)
    {
        read(n),read(k),read(p);
        int w=qp(getg(p),(p-1)/k),ans=0;
        for(int j=0;j<k;j++)
        {
            A.a=A.b=A.c=qp(w,j),A.d=1;
            A.a=add(A.a,1);
            A=mqp(A,n);
            ans=add(ans,A.a);
        }
        printf("%lld\n",mul(ans,qp(k,p-2)));    
    }
    return 0;
}

2019.5.8

转载于:https://www.cnblogs.com/butterflydew/p/10830216.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值