鸣谢
http://www.myexception.cn/program/1879933.html
http://blog.csdn.net/popoqqq/article/details/39891263
题意:链接
方法:各种数论
解析:我希望我以后不要再看到数论题TAT
好,又是看了各种题解后写的一道题,大爷的题解说的好厉害但是好像并不是那么显而易见?所以不如自己来写一发,随便理解就好了。
首先呢要明确这道题要你干什么。
就是有n个东西,然后m个人,给每个人w[i]个物品,问你方案数。
随便写写就会归纳成
c(n,sum)∗c(sum,w[1])∗c(sum−w[1],w[2])...
一直乘下去就好,最终取余p。
我靠这不就是lucas定理么,这么(哔)的一道题还有啥做的。
等等样例?这个p居然不是质数!(谁告诉你是质数了!)
噢难道我可以把这道题理解为c(n,m)终极版么- -
好,步入正题,既然知道了如此坑爹的条件后我们怎么做这道题呢?
首先,我们可以把p分解质因数,变成
p=p1a1∗p2a2...pnan
其中,pi均为质数。
之后怎么搞呢?对的!中国剩余定理!
将每个 piai 单独计算最终再用天朝剩余定理合并就好了。
可是组合数怎么计算呢?
c(n,m)我们把他看做三个阶乘,分别为n!,m!,(n-m)!
但是如果直接用逆元的话并不可以,因为有可能在这些阶乘中出现取模的数,这样的话就直接变为0了。所以我们要把这些数提取出来。
引用一名大牛的方法
n!=t1∗piu1
m!=t2∗piu2
(n−m)!=t3∗piu3
所以t1,t2,t3我们怎么求呢?
假设19!,pi=3(popoqqq大爷举的栗子)
19!=1*2…*19
19!=1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19∗36∗(1∗2∗3∗4∗5∗6)
所以后面的123456又是个子问题可以递归求解,而这个前面的乱七八糟一堆东西我们只需要把它乘到一起就好了!并且我们提取出来了6个pi,现在继续递归下去的话,乱七八糟的东西再次乘到一起,pi继续提取,然后继续递归就搞定n!了!
现在我们只需要 t1∗niyuan(t2)∗niyuan(t3)∗piu1−u2−u3
这样我们就可以搞定所有的c(n,m)最终乘到一起就好了!
代码:
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define N 100010
#define mp make_pair
#define pa pair<ll,ll>
#define fi first
#define se second
using namespace std ;
typedef long long ll;
ll p,n,m,w[N];
ll prime[N],mod[N],cnt[N],a[N];
int tot;
ll quick_my(ll a,ll b,ll M)
{
ll ret=1;
while(b)
{
if(b&1)ret=(ret*a)%M;
a=(a*a)%M;
b>>=1;
}
return ret;
}
void exgcd(ll a,ll b,ll &x,ll &y,ll &gcd)
{
if(!b)
{
x=1,y=0,gcd=a;
return ;
}
exgcd(b,a%b,y,x,gcd);
y=y-a/b*x;
}
void get_factor(ll x)
{
for(int i=2;i*i<=x;i++)
{
if(x%i==0)
{
prime[++tot]=i,cnt[tot]=0,mod[tot]=1;
while(x%i==0){x/=i,cnt[tot]++,mod[tot]*=i;}
}
}
if(x>1)prime[++tot]=x,cnt[tot]=1,mod[tot]=x;
}
ll inverse(ll a,ll b)
{
ll xx,yy,d;
exgcd(a,b,xx,yy,d);
return (xx%b+b)%b;
}
pa fact(ll k,ll n)
{
if(n==0)return mp(0,1);
int x=n/prime[k],y=n/mod[k];
ll ans=1;
if(y)
{
for(int i=2;i<mod[k];i++)
{
if(i%prime[k]!=0)ans=(ans*1ll*i)%mod[k];
}
ans=quick_my(ans,y,mod[k]);
}
for(int i=y*mod[k]+1;i<=n;i++)
{
if(i%prime[k]!=0)ans=ans*1ll*i%p;
}
pa tmp=fact(k,x);
return mp(x+tmp.fi,ans*tmp.se%p);
}
ll calc(int k,ll n,ll m)
{
if(n<m)return 0;
pa a=fact(k,n),b=fact(k,m),c=fact(k,n-m);
return quick_my(prime[k],a.fi-b.fi-c.fi,mod[k])*a.se%mod[k]*inverse(b.se,mod[k])%mod[k]*inverse(c.se,mod[k])%mod[k];
}
ll china()
{
ll gcd,y,x=0;
for(int i=1;i<=tot;i++)
{
ll r=p/mod[i];
exgcd(mod[i],r,gcd,y,gcd);
x=(x+r*y*a[i])%p;
}
return (x+p)%p;
}
ll work(ll n,ll m)
{
for(int i=1;i<=tot;i++)a[i]=calc(i,n,m);
return china();
}
int main()
{
scanf("%lld%lld%lld",&p,&n,&m);
get_factor(p);
int sum=0;
for(int i=1;i<=m;i++)scanf("%lld",&w[i]),sum+=w[i];
if(sum>n){printf("Impossible\n");return 0;}
ll ans=work(n,sum)%p;
for(int i=1;i<=m;i++)
{
ans=ans*work(sum,w[i])%p;
sum-=w[i];
}
printf("%lld\n",ans);
}