bzoj1951 [Sdoi2010]古代猪文(Lucas+CRT+exgcd)

Description

“在那山的那边海的那边有一群小肥猪。他们活泼又聪明,他们调皮又灵敏。他们自由自在生活在那绿色的大草坪,他们善良勇敢相互都关心……” ——选自猪王国民歌 很久很久以前,在山的那边海的那边的某片风水宝地曾经存在过一个猪王国。猪王国地理位置偏僻,实施的是适应当时社会的自给自足的庄园经济,很少与外界联系,商贸活动就更少了。因此也很少有其他动物知道这样一个王国。 猪王国虽然不大,但是土地肥沃,屋舍俨然。如果一定要拿什么与之相比的话,那就只能是东晋陶渊明笔下的大家想象中的桃花源了。猪王勤政爱民,猪民安居乐业,邻里和睦相处,国家秩序井然,经济欣欣向荣,社会和谐稳定。和谐的社会带给猪民们对工作火红的热情和对未来的粉色的憧憬。 小猪iPig是猪王国的一个很普通的公民。小猪今年10岁了,在大肥猪学校上小学三年级。和大多数猪一样,他不是很聪明,因此经常遇到很多或者稀奇古怪或者旁人看来轻而易举的事情令他大伤脑筋。小猪后来参加了全猪信息学奥林匹克竞赛(Pig Olympiad in Informatics, POI),取得了不错的名次,最终保送进入了猪王国大学(Pig Kingdom University, PKU)深造。 现在的小猪已经能用计算机解决简单的问题了,比如能用P++语言编写程序计算出A + B的值。这个“成就”已经成为了他津津乐道的话题。当然,不明真相的同学们也开始对他刮目相看啦~ 小猪的故事就将从此展开,伴随大家两天时间,希望大家能够喜欢小猪。 题目描述 猪王国的文明源远流长,博大精深。 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的余数就可以了。

Input
有且仅有一行:两个数N、G,用一个空格分开。

Output
有且仅有一行:一个数,表示答案除以999911659的余数。

Sample Input
4 2

Sample Output
2048

HINT
10%的数据中,1 <= N <= 50;
20%的数据中,1 <= N <= 1000;
40%的数据中,1 <= N <= 100000;
100%的数据中,1 <= G <= 1000000000,1 <= N <= 1000000000。

分析:
我都懒得喷题面了

很显然答案就是
G^ΣC(n,i)%p

1)有这样一个神奇的公式

费马定理
A^x = A^(x % φ(p) + φ(p)) (mod p) x>=p

简单来说就是如果模数是一质数,在计算快速幂的时候,可以直接把指数%(p-1)

由于p是质数,φ(p)=p-1=999911658

2)所以现在要求的是 ΣC(n,i)%(p-1) {实际上就是指数%(p-1)}

3)φ(p)=p-1=999911658=2*3*4679*35617
每个质因子只有一个
所以我们可以求出ΣC(n,i)%2,ΣC(n,i)%3,ΣC(n,i)%4679,ΣC(n,i)%35617
然后用中国剩余定理合并起来

4)复习中国剩余定理(不如说是重新学习)
图:
这里写图片描述
中国剩余定理,简称CRT
有若干个同余方程,形如:
x%mi=ai 且 ai互质的方程组
M=π(mi)
Mi=M/mi
Mi*ki≡1 (mod mi)
ki和Mi是mod mi意义下的逆元
mod M意义下的唯一解就是:

ans=Σ(Mi*ki*ai)

5)对于n,m>=1e5 可以用Lucas定理
C(n,m)=Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)

特别的

n < m时C(n,m)=0

这也算是一道挺复杂的题了
下面就说一说每一部分中需要注意的问题

首先就是4个计算之前的预处理
这几个计算都是求组合数,显然,我们需要预处理阶乘和阶乘逆元
这里写图片描述
这里需要注意,我的循环只到p-1,因为在Lucas中,每次我都要进行mod p操作
所以真正进行组合数计算的数字一定是小于p的
逆元就是喜闻乐见的线性求法了

之后主程序中枚举n的所有约数

函数doit就是进行了4个单独的计算

注意

每一部分的模数都不同,分别是2,3,4679,35617

喜闻乐见的Lucas,不用多说
这里写图片描述
可见在预处理之后,组合数就可以O(1)出解了

比较困难的就是中国剩余定理的实现
M就是p-1
实际上仔细看的话,其实程序中的语句把分析中的所有操作都模拟了

Mi*ki≡1 (mod mi)
等价于:Mi*ki+y*mi=1 (mod mi)
已知Mi=phimod/mi,和mi
exgcd求解ki就可以了

注意Σ答案一定是对phimod取模

这里写图片描述
exgcd用的是pair存储的,其实就是一个华丽的外表

这里我们用了一个高级的表达:typedef
实际上就是给pair<>起了一个简单的名字

tip

在计算每一部分的时候,模数都是该部分特有的
这些函数包括:
doit(),Lucas(),C(),makepreCRT()
我又无意中发现,mmm也在写有exgcd的题,提高警惕!!!不颓废fighting

第一遍写完之后,又是狂WA
发现还是主程序的问题,在枚举约数的时候写崩了

这里写代码片
#include<cstdio>
#include<iostream>
#include<iostream>
#include<cmath>
#define ll long long

using namespace std;

const ll mod=999911659;
const ll phimod=999911658;   //999911658=2*3*4679*35617
ll n,g,ans[5];
ll num[4]={2,3,4679,35617},jc[4][35618],inv[4][35618];
typedef pair<ll,ll> abcd;

void makepreCRT(ll p,ll *jc,ll *inv)   //每一部分的模数都不一样 
{
    ll i;
    jc[0]=1;
    for (i=1;i<p;i++)    //这里为什么只到p-1,因为Lucas的%操作,使得每个数都小于p 
        jc[i]=(jc[i-1]%p*i%p)%p;
    inv[1]=1;
    for (i=2;i<p;i++) 
        inv[i]=((p-p/i)%p*inv[p%i]%p)%p;  
    inv[0]=1; 
    for (i=1;i<=p;i++)
        inv[i]=(inv[i-1]%p*inv[i]%p)%p;   //阶乘的逆元 
}

ll KSM(ll a,ll b,ll p)
{
    ll t=1;
    a%=p;
    while (b)
    {
        if (b&1)
           t=(t%p*a%p)%p;
        b>>=1;
        a=(a%p*a%p)%p;
    }
    return t%p;
}

ll C(ll n,ll m,ll p,ll *jc,ll *inv)
{
    if (m>n) return 0;
    return (jc[n]%p*inv[m]%p*inv[n-m]%p)%p; 
}

ll Lucas(ll n,ll m,ll p,ll *jc,ll *inv)
{
    if (m>n) return 0;
    ll ans=1;
    while (m)
    {
        ans=(ans%p*C(n%p,m%p,p,jc,inv)%p)%p;
        n/=p; m/=p;
    }
    return ans%p;
}

void doit(ll x)
{
    int i;
    for (i=0;i<4;i++)
    {
        ans[i]+=Lucas(n,x,num[i],jc[i],inv[i]);
        ans[i]%=num[i];   //在计算每一部分的时候,模数都是该部分特有的 
    }
}

abcd exgcd(ll a,ll b)
{
    if (!b) return abcd(1,0);
    abcd temp=exgcd(b,a%b);
    return abcd(temp.second,temp.first-(a/b)*temp.second);
}

ll CRT()
{
    int i;
    ll re=0;
    for(i=0;i<4;i++)
    {
        abcd temp=exgcd(phimod/num[i],num[i]);    //Mi mi
        ll x=(temp.first%phimod*(phimod/num[i])%phimod+phimod)%phimod;   //Mi*ki  mod M
        re+=x*ans[i]%phimod;
        re%=phimod;   //Σ(Mi*ki*ai) mod M
    }
    return re;
}

int main()
{
    for (int i=0;i<4;i++)
        makepreCRT(num[i],jc[i],inv[i]);   //为CRT做准备 
    scanf("%lld%lld",&n,&g);
    for (ll i=1;i*i<=n;i++)
    {
        if (n%i==0)
        {
            doit(i);
            if (n/i!=i) doit(n/i); 
        }
    }  
    ans[4]=CRT();
    ans[4]=KSM(g,ans[4]+phimod,mod);   //A^x = A^(x % φ(p) + φ(p))  
    printf("%lld",ans[4]);
    return 0;
}
原文地址:https://www.cnblogs.com/wutongtong3117/p/7673236.html