[FZU 1833 Counting Problem]数论、组合数、扩展欧几里得算法解模线性方程

【题目地址】http://acm.fzu.edu.cn/problem.php?pid=1833
【题目大意】
求P(sum(Ai),n*n) / (A1! * A2! * A3! …… * Ak!)
根据出题人福大核武的指示，本题考的是数论，并不是组合数学，而由于这个公式已经水得连我都能推出来，就直接写算了。。。
【算法分析】
首先做之前应该清楚几点：
1、一个数n的素因子是O(lg n)级别的。
2、没有任何辅助工具(素数表之类)的时候，将一个数分解质因数的复杂度是O(sqrt(n))级别的。
3、AC大神的 n! 质因数分解模板的复杂度是好吧，你也许会觉得这非常脑残，但是，我就是脑残地没搞清。。。

做法：
1、对M分解质因数
2、对公式的分母中红字部分分解质因数，并且分开自己拥有的质因数和M所拥有的质因数。
3、由于Sum<=100*5000，所以可以暴力求P(sum,n*n)。就是直接从n*n乘到n*n-sum+1。
乘的过程中把M所含的质因子也分离开来，和上面的那个蓝色的进行幂减。幂减以后，把各部分再分别合起来，弄成分子s1和分母s2，我们就可以惊讶地发现。。。gcd(s2，M)=1。
4、由于gcd(s2，M)=1，那么对应的 s1 = s2*x (mod M) 的解就变成了唯一的。（这个不知道对不对=_=，不会证明）
5、对于那个模线性方程，变形得 s1= s2*x + M*y。这个用欧几里得拓展算法搞定。

【其他】
Orz AC大神~
据其指示，时间和第一名居然是一样的~
1WA、1TLE。。。
WA：
最后一步搞了一会儿，一开始以为直接 s1 * s2^(M-2) % M，但是原来M必须为质数才能搞。
上面的第4步得到的gcd(s2,M)=1并不能让这个公式成立。。。
TLE：
我脑残地不听AC大神的劝告，直接枚举答案求逆元去了。。。

再Orz 一下我的各种脑残。。。果然数论我已经水到一种境界了。。。
【CODE】
#include #include #include #include using namespace std;
typedef long long lld;
int n,mod,cn,mtot,tot2;
int a[105];
int pl[50],cl[50],pl2[5005],cl2[5005];
bool ss[5005],hash[5005];
lld s1,s2;

void init(){
int i,j;
tot2=0;
memset(ss,true,sizeof(ss));
for (i=2;i*i<=5000;i++)
for (j=2;i*j<=5000;j++)
ss[i*j]=false;
for (i=2;i<=5000;i++)
if (!hash[i] && ss[i])
pl2[++tot2]=i;
}

void dealm(){
int i,j,x=mod;
bool flag;
for (i=2;i*i<=mod;i++){
flag=true;
while (x%i==0){
if (i<=5000) hash[i]=true;
if (flag) pl[++mtot]=i;
x/=i;
flag=false;
}
}
if (x!=1){
pl[++mtot]=x;
if (x<=5000) hash[x]=true;
}
}

void Div(int *list,int n){
if (!n) return;
int i,tmp;
for (i=1;i<=mtot && pl[i]<=n;i++){
tmp=n;
while (tmp/pl[i]){
list[i]-=tmp/pl[i];
tmp/=pl[i];
}
}
}

void Div2(int *list,int n){
if (!n) return;
int i,tmp;
for (i=1;i<=tot2 && pl2[i]<=n;i++){
tmp=n;
while (tmp/pl2[i]){
list[i]+=tmp/pl2[i];
tmp/=pl2[i];
}
}
}

void Fen(int x){
if (!x) return;
int i;
for (i=1;i<=mtot && x>=pl[i];i++)
while (x%pl[i]==0){
cl[i]++;
x/=pl[i];
}
s1=s1*x%mod;
}

lld pow_mod(lld x,lld cf){
if (cf==0) return 1;
lld res=pow_mod(x,cf/2);
res=res*res%mod;
if (cf&1) res=res*x%mod;
return res;
}

lld exgcd(lld a,lld b,lld &x,lld &y){
if (b==0){
x=1; y=0;
return a;
}
lld res=exgcd(b,a%b,x,y),tx=x,ty=y;
x=ty; y=tx-a/b*ty;
return res;
}

void solve(){
lld x,y,rate=s1/exgcd(s2,mod,x,y);
x=x%mod*rate%mod;
x=(x+mod)%mod;
cout << x << endl;
}

int main(){
while (scanf("%d%d%d",&n,&mod,&cn)!=EOF){
mtot=0; tot2=0;
s1=1; s2=1;
memset(cl,0,sizeof(cl));
memset(hash,false,sizeof(hash));
memset(cl2,0,sizeof(cl2));
int sum=0;
for (int i=1;i<=cn;i++){
scanf("%d",&a[i]);
sum+=a[i];
}
if (mod==1){
cout << 0 << endl;
continue;
}
dealm();
init();
for (int i=1;i<=cn;i++){
Div(cl,a[i]);
Div2(cl2,a[i]);
}
for (int i=n*n;i>n*n-sum;i–) Fen(i);
for (int i=1;i<=mtot;i++)
s1=s1*pow_mod(pl[i],cl[i])%mod;
for (int i=1;i<=tot2;i++)
s2=s2*pow_mod(pl2[i],cl2[i])%mod;
solve();
}
}