bzoj 3328: PYXFIB 数论

3328: PYXFIB

Time Limit: 10 Sec  Memory Limit: 256 MB
Submit: 130  Solved: 41
[Submit][Status][Discuss]

Description

Input

第一行一个正整数，表示数据组数据 ，接下来T行
每行三个正整数N,K,P

Output

T行，每行输出一个整数，表示结果

Sample Input

1
1 2 3

Sample Output

1

HINT

Source

　　思路与莫比乌斯反演相似，通过二项式巧妙地解决组合数的问题，总结一个技巧：对于同余系P，g为原根，w=g^((p-1)/k)，那么sigma(w^(ij))==[i mod k==0]，其他一些非人类的插值，代换这里就不多说了。太科幻了。。。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long qword;
int mod;
qword pow_mod(qword x,qword y,qword mod=::mod)
{
        qword ret=1;
        while (y)
        {
                if (y&1)
                        ret=ret*x%mod;
                x=x*x%mod;
                y>>=1;
        }
        return ret;
}
struct matrix
{
        qword mat[2][2];
        matrix()
        {
                memset(mat,0,sizeof(mat));
        }
        void Init0()
        {
                memset(mat,0,sizeof(mat));
                mat[0][0]=mat[1][1]=1;
        }
        void Init1()
        {
                memset(mat,0,sizeof(mat));
                mat[1][0]=1;
                mat[0][0]=mat[0][1]=1;
        }
        void Print()
        {
                for (int i=0;i<2;i++)
                {
                        for (int j=0;j<2;j++)
                        {
                                printf("%lld ",mat[i][j]);
                        }
                        printf("
");
                }
                printf("
");
        }
};
matrix operator *(matrix m1,matrix m2)
{
        matrix ret;
        for (int k=0;k<2;k++)
        {
                for (int i=0;i<2;i++)
                {
                        for (int j=0;j<2;j++)
                        {
                                ret.mat[i][j]=(ret.mat[i][j]+m1.mat[i][k]*m2.mat[k][j]%mod)%mod;
                        }
                }
        }
        return ret;
}
matrix operator *(matrix m1,qword k)
{
        for (int i=0;i<2;i++)
                for (int j=0;j<2;j++)
                        m1.mat[i][j]=m1.mat[i][j]*k%mod;
        return m1;
}
matrix operator +(matrix m1,matrix m2)
{
        for (int i=0;i<2;i++)
                for (int j=0;j<2;j++)
                        m1.mat[i][j]=(m1.mat[i][j]+m2.mat[i][j])%mod;
        return m1;
}

matrix matrix_pow(matrix m1,qword y)
{
        matrix res;
        res.Init0();
        while (y)
        {
                if (y&1)
                        res=res*m1;
                m1=m1*m1;
                y>>=1;
        }
        return res;
}
void Analyse(int x,vector<int> &vec)
{
        for (int i=1;i*i<=x;i++)
        {
                if (i*i==x){
                        vec.push_back(i);
                        sort(vec.begin(),vec.end());
                        return ;
                }
                if (x%i==0)
                {
                        vec.push_back(i);
                        vec.push_back(x/i);
                }
        }
        sort(vec.begin(),vec.end());
}


int main()
{
    //    freopen("input.txt","r",stdin);
        int nn;
        scanf("%d",&nn);
        while (nn--)
        {
                qword n;
                int t,p;
                scanf("%lld%d%d",&n,&t,&p);
                mod=p;
                int g=-1;//原根
                vector<int> fpm1;
                Analyse(p-1,fpm1);
        //        for (int i=0;i<fpm1.size();i++)
        //                printf("%d ",fpm1[i]);
        //        printf("
");
                for (int i=1;i<p;i++)
                {
                        bool flag=true;
                        for (int j=0;j<(int)fpm1.size()-1;j++)
                        {
                                if (pow_mod(i,fpm1[j])==1)
                                {
                                        flag=false;
                                        break;
                                }
                        }
                        if (flag)
                        {
                                g=i;
                                break;
                        }
                }
                /*
                qword x0=1;
                for (int i=1;i<p-2;i++)
                {
                        x0=x0*g;
                        if (x0==1)throw 1;
                }*/
                qword w=pow_mod(g,(p-1)/t);
                qword invw=pow_mod(w,mod-2);
                qword xnow=1;
                qword ixnow=1;
                matrix res;
                for (int i=0;i<t;i++)
                {
                        matrix matx;
                        matx.Init1();
                        matx.mat[0][0]=(matx.mat[0][0]+xnow)%mod;
                        matx.mat[1][1]=(matx.mat[1][1]+xnow)%mod;
                        matx=matrix_pow(matx,n);
                        matx=matx*pow_mod(ixnow,n);
                        res=res+matx;
                        xnow=xnow*invw%mod;
                        ixnow=ixnow*w%mod;
                }
                res=res*pow_mod(t,mod-2);
                printf("%lld
",res.mat[0][0]);
        }
}