bzoj 3328: PYXFIB unit root inversion

Posted by jblack on Sun, 24 Nov 2019 20:49:04 +0100

Description

Input
The first line is a positive integer, representing the data group data, and the next T line
Three positive integers per line N,K,P

Output
T line, output an integer for each line, representing the result

Sample Input
1
1 2 3

Sample Output
1

HINT

Source
By Wcmg

Analysis:
The number of combinations is very large, so the number of combinations should be taken into account.
Simplify the appeal formula to get
=∑i=0n[i mod k==0](ni)∗F(i)=\sum_{i=0}^{n}[i\ mod\ k==0]\binom{n}{i}*F(i)=i=0∑n​[i mod k==0](in​)∗F(i)
Considering unit root inversion, we get
=∑i=0n1k∗∑j=0k−1(wkj)i∗(ni)∗F(i)=\sum_{i=0}^{n}\frac{1}{k}*\sum_{j=0}^{k-1}(w_{k}^j)^i*\binom{n}{i}*F(i)=i=0∑n​k1​∗j=0∑k−1​(wkj​)i∗(in​)∗F(i)
Exchange cycle,
=1k∗∑j=0k−1∑i=0n(wkj)i∗(ni)∗F(i)=\frac{1}{k}*\sum_{j=0}^{k-1}\sum_{i=0}^{n}(w_{k}^j)^i*\binom{n}{i}*F(i)=k1​∗j=0∑k−1​i=0∑n​(wkj​)i∗(in​)∗F(i)
Where, F(i)=AiF(i)=A^iF(i)=Ai, AAA is the recurrence matrix of Fibonacci sequence. Then there is a binomial expansion.
=1k∗∑j=0k−1(wkjA+I)n=\frac{1}{k}*\sum_{j=0}^{k-1}(w_{k}^jA+I)^n=k1​∗j=0∑k−1​(wkj​A+I)n
The direct O(klogn)O(klogn)O(klogn) is good.

Code:

/**************************************************************
    Problem: 3328
    User: liangzihao
    Language: C++
    Result: Accepted
    Time:8932 ms
    Memory:1300 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long
 
using namespace std;
 
LL G,n,w,wn,ans,mod;
int T,k,cnt;
int p[101];
 
struct matrix{
    LL a[3][3];
}A,B;
 
matrix operator *(matrix a,matrix b)
{
    matrix c;
    for (int i=0;i<=2;i++)
    {
        for (int j=0;j<=2;j++) c.a[i][j]=0;
    }
    for (int k=1;k<=2;k++)
    {
        for (int i=1;i<=2;i++)
        {
            for (int j=1;j<=2;j++) c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
        }
    }
    return c;
}
 
LL ksm(LL x,LL y)
{
    if (y==1) return x;
    LL c=ksm(x,y/2);
    c=c*c%mod;
    if (y&1) c=c*x%mod;
    return c;
}
 
void divide(int x)
{
    for (int i=2;i<=trunc(sqrt(x));i++)
    {
        if (x%i==0)
        {
            p[++cnt]=i;
            while (x%i==0) x/=i;
        }
    }
    if (x>1) p[++cnt]=x;
}
 
void findroot(int x)
{
    for (int i=2;i<x;i++)
    {
        int flag=0;
        for (int j=1;j<=cnt;j++)
        {
            if (ksm(i,(x-1)/p[j])==1)
            {
                flag=1;
                break;
            }
        }
        if (!flag)
        {
            G=i;
            return;
        }
    }
}
 
void power(LL p)
{
    if (p==1)
    {
        B=A;
        return;
    }
    power(p/2);
    B=B*B;
    if (p&1) B=B*A;
}
 
int main()
{
    scanf("%d",&T);
    while (T--)
    {
        scanf("%lld%d%lld",&n,&k,&mod);
        cnt=0;      
        divide(mod-1);
        findroot(mod);
        wn=ksm(G,(mod-1)/k);
        w=1;
        ans=0;
        for (int j=0;j<k;j++)
        {
            A.a[1][1]=1,A.a[1][2]=w;
            A.a[2][1]=w,A.a[2][2]=(w+1)%mod;
            power(n);
            ans=(ans+B.a[2][2])%mod;
            w=(w*wn)%mod;
        }
        ans=ans*ksm(k,mod-2)%mod;
        printf("%lld\n",ans);
    }
}