BZOJ1951 || 洛谷P2480 [SDOI2010]古代猪文【Lucas定理+CRT+欧拉定理】

niiick

于 2019-04-02 13:14:21 发布

阅读量181

点赞数

分类专栏：数学;数论

本文链接：https://blog.csdn.net/niiick/article/details/88969862

版权

数学;数论专栏收录该内容

36 篇文章 3 订阅

订阅专栏

时空限制 1000ms / 128MB

题目描述

猪王国的文明源远流长，博大精深。

iPig在大肥猪学校图书馆中查阅资料，得知远古时期猪文文字总个数为N。当然，一种语言如果字数很多，字典也相应会很大。当时的猪王国国王考虑到如果修一本字典，规模有可能远远超过康熙字典，花费的猪力、物力将难以估量。故考虑再三没有进行这一项劳猪伤财之举。当然，猪王国的文字后来随着历史变迁逐渐进行了简化，去掉了一些不常用的字。

iPig打算研究古时某个朝代的猪文文字。根据相关文献记载，那个朝代流传的猪文文字恰好为远古时期的k分之一，其中k是N的一个正约数（可以是1和N）。不过具体是哪k分之一，以及k是多少，由于历史过于久远，已经无从考证了。

iPig觉得只要符合文献，每一种能整除N的k都是有可能的。他打算考虑到所有可能的k。显然当k等于某个定值时，该朝的猪文文字个数为N / k。然而从N个文字中保留下N / k个的情况也是相当多的。iPig预计，如果所有可能的k的所有情况数加起来为P的话，那么他研究古代文字的代价将会是G的P次方。

现在他想知道猪王国研究古代文字的代价是多少。由于iPig觉得这个数字可能是天文数字，所以你只需要告诉他答案除以999911659的余数就可以了。

输入格式：

输入文件ancient.in有且仅有一行：两个数N、G，用一个空格分开。

输出格式：

输出文件ancient.out有且仅有一行：一个数，表示答案除以999911659的余数。

说明

10%的数据中，1 <= N <= 50；
20%的数据中，1 <= N <= 1000；
40%的数据中，1 <= N <= 100000；
100%的数据中，1 <= G <= 1000000000，1 <= N <= 1000000000

题目分析

神题预定

首先一句话题意 $G^{\sum_{k|d}C_n^k}\bmod 999911659$
因为p=999911659是素数，所以根据欧拉定理的推论有 $G^{\sum_{k|d}C_n^k\bmod \varphi(p)}\bmod p$
其中 $\varphi(p)=p-1=999911658$
所以难点就在于求 $\sum_{k|d}C_n^k\bmod 999911658$

接下来就是神思路
999911658分解质因数为 $2 * 3 * 4679 * 35617$ ，发现每个质因子指数都是1
于是可以分别求 $a_i=\sum_{k|d}C_n^k\bmod p_i$
与 $x=\sum_{k|d}C_n^k\bmod 999911658$ 组成同余方程 $x\equiv a_i(\mod p_i)$
用中国剩余定理合并四个同余方程得x再快速幂即可

注意特判G=999911659得情况

#include<iostream>
#include<cstdio>
#include<cmath>
#include<queue>
#include<algorithm>
#include<cstring>
#include<iomanip>
using namespace std;
typedef long long lt;
typedef double dd;

int read()
{
    int f=1,x=0;
    char ss=getchar();
    while(ss<'0'||ss>'9'){if(ss=='-')f=-1;ss=getchar();}
    while(ss>='0'&&ss<='9'){x=x*10+ss-'0';ss=getchar();}
    return f*x;
}

const lt mod=999911659;
const int maxn=100010;
lt b[]={0,2,3,4679,35617};
lt a[15],fac[maxn];

lt qpow(lt a,lt k,lt p)
{
    lt res=1;
    while(k>0){
        if(k&1) res=(res*a)%p;
        a=(a*a)%p; k>>=1;
    }
    return res;
}

lt C(lt n,lt m,lt p)
{
    if(n<m) return 0;
    return fac[n]*qpow(fac[m],p-2,p)%p*qpow(fac[n-m],p-2,p)%p;
}

lt lucas(lt n,lt m,lt p)
{
    if(n<m) return 0;
    if(m==0) return 1;
    return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}

void init(lt p)
{
    fac[0]=1;
    for(lt i=1;i<=p;i++)
    fac[i]=fac[i-1]*i%p;
}

lt exgcd(lt a,lt b,lt &x,lt &y)
{
    if(b==0){ x=1; y=0; return a;}
    lt g=exgcd(b,a%b,x,y);
    lt tp=x; x=y; y=tp-a/b*y;
    return g;
}

lt exCRT()
{
    lt x,y;
    lt M=b[1],res=a[1];
    for(int i=2;i<=4;i++)
    {
        lt ai=M,bi=b[i],c=(a[i]-res%bi+bi)%bi;
        lt gcd=exgcd(ai,bi,x,y),bg=bi/gcd;
        if(c%gcd!=0) return -1; 
        
        x=x*c/gcd%bg;
        res+=x*M; M*=bg;
        res=(res%M+M)%M;
    }
    return (res%M+M)%M;
}

int main()
{
    lt N=read(),G=read(); 
    if(G==mod){ printf("0"); return 0;}
    
    for(int i=1;i<=4;++i)
    {
        int k; init(b[i]);
        for(k=1;k*k<=N;++k)
        if(N%k==0){
            a[i]=(a[i]+lucas(N,k,b[i]))%b[i];
            if(k*k!=N)a[i]=(a[i]+lucas(N,N/k,b[i]))%b[i];
        }
    }
    
    lt ans=qpow(G,exCRT(),mod);
    printf("%lld",ans);
    return 0;
}