コンビネーション数モデリングのまとめ
本稿では,C(N,M)modPについて,異なるデータ範囲に対する3つのアルゴリズムを紹介する.
例一:FZU 2020
データ範囲n,m,p(1<=m<=n<=10^9,m<=10^4,mC(N,M)=N*(N-1)*(N-2)*...*(N-M+1)/M!
mは非常に小さく,分子分母にはm項があることに気づき,O(m)暴力で分子分母mod pの値を求めることができる.
またm時間複雑度O(m+log 2 p),空間複雑度O(1).
code:
例二:BZOJ 2982
データ範囲1<=m<=n<=2,000,000,000,p=10007.
本題の割合は1 mの範囲を2*10^8に増やした.pは10007しかなく、質量数である.
そこでlucasの定理を直接運用することができる.
Lucasの定理:
where
and
これがあれば、m,nを直接p進法に分割することができます.
C(Nk,Mk)mod pの求め方は,例1のO(P)時間で解くことに変換できる.
時間複雑度O(p),空間複雑度O(p).
code:
例3:HDU 3439
nの全配列の中でn-k転位列の方案数mod pだけあることを求めます.
データ範囲(0<=k<=n<=10^9,1<=p<=10^5,n!=0)
ターゲットの答えはC(n,k)*f(n-k)mod pで、f(i)はiです!全配列誤りスキーム数.
f(n−k)mod pは直接ループノードを探すことができ,ループノード長はm*2である.
次にC(n,k)mod pを要求する.本題ではpはまだ小さいが,pは素数を保証しないため,例2の方法は使用できない.
どうしようかな、まずp分解素因数p=p 1^k 1*p 2^k 2.*.*pt^kt.そしてc(n,k)mod pi^kiの値をそれぞれ求め,中国の残りの定理を用いて方程式群と同列に統合すればよい.
肝心なのはc(n,k)mod pi^kiを求めることだ.
c(n,k)=n!/(k!(n-k)!).分子分母中のすべての素因数piを蹴り出すと、分母はpi^kiと互いに質を交換し、拡張ユークリッドでmod pi^kiの逆元を求めればよい.
O(logPi)の時間内に分子分母に含まれる素因数piの個数を求め,それらの差をTとすることができる.
そしてすべての問題はnを求めることに転化します!素因数pを蹴ったmod p^cの値にpi^Tを乗じます.
前処理F[i]は、1からiまでpで割り切れない数の積を表す.P=p^cとする.1からnを求めてP個ごとに節を分けます.各グループの前p-1個の数mod Pは同余であるため、素因数pを含まない積はF[n%P]*F[P-1]^(n/P)であり、残りの素因数pの数に対してpを除去した後に積は(n/p)!に変換される.このようにして同じサブ問題に変換します!
1つのインスタンスでシミュレーション:20!Mod 9.
20!=(1*2*4*5*7*8*…*16*17*19)*(3*6*9*12*15*18)=(1*2*4*5*7*8)^2*(19*20)*3^6(1*2*3*4*5*6)=F[9-1]^(20/9)*F[20%9]*((20/3)! Mod 9).
ここで、F[0]からF[9−1]までの値は1,1,2,2,8,4,4,1,8である.
このように問題は完璧に解決された.
時間的複雑度はO(Σpi^ki).空間複雑度O(m).
code:
例一:FZU 2020
データ範囲n,m,p(1<=m<=n<=10^9,m<=10^4,m
mは非常に小さく,分子分母にはm項があることに気づき,O(m)暴力で分子分母mod pの値を求めることができる.
またm
code:
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <string>
using namespace std;
typedef long long LL;
LL n,k,p,a,b;
int get()
{
int f=0,v=0;char ch;
while(!isdigit(ch=getchar()))if(ch=='-')break;
if(ch=='-')f=1;else v=ch-48;
while(isdigit(ch=getchar()))v=v*10+ch-48;
if(f==1)return -v;else return v;
}
LL Pow(LL a,LL b,LL p)
{
LL res=1;
for(;b>0;b/=2)
{
if(b%2==1)res=res*a%p;
a=a*a%p;
}
return res;
}
int main()
{
int t=get();
while(t--)
{
n=get(),k=get(),p=get();
a=1,b=1;
for(int i=1;i<=k;i++)a=a*(n-i+1)%p,b=b*i%p;
b=Pow(b,p-2,p);
printf("%d
",a*b%p);
}
return 0;
}
例二:BZOJ 2982
データ範囲1<=m<=n<=2,000,000,000,p=10007.
本題の割合は1 mの範囲を2*10^8に増やした.pは10007しかなく、質量数である.
そこでlucasの定理を直接運用することができる.
Lucasの定理:
where
and
これがあれば、m,nを直接p進法に分割することができます.
C(Nk,Mk)mod pの求め方は,例1のO(P)時間で解くことに変換できる.
時間複雑度O(p),空間複雑度O(p).
code:
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <string>
using namespace std;
const int p=10007;
int n,k,inv[p],fac[p],a[10],b[10],len,ans;
int get()
{
int f=0,v=0;char ch;
while(!isdigit(ch=getchar()))if(ch=='-')break;
if(ch=='-')f=1;else v=ch-48;
while(isdigit(ch=getchar()))v=v*10+ch-48;
if(f==1)return -v;else return v;
}
int Pow(int a,int b)
{
int res=1;
for(;b>0;b/=2)
{
if(b%2==1)res=res*a%p;
a=a*a%p;
}
return res;
}
int C(int n,int k)
{
if(k>n)return 0;
if(k==0||n==0||n==k)return 1;
return fac[n]*inv[k]%p*inv[n-k]%p;
}
int main()
{
int T=1;
for(int i=1;i<p;i++)T=T*i%p,inv[i]=Pow(T,p-2),fac[i]=T;
T=get();
while(T--)
{
n=get(),k=get(); ans=1;
memset(a,0,sizeof(a));
memset(b,0,sizeof(b));
for(len=0;k>0;k/=p)a[++len]=k%p;
for(len=0;n>0;n/=p)b[++len]=n%p;
for(int i=1;i<=len;i++)ans=ans*C(b[i],a[i])%p;
printf("%d
",ans);
}
return 0;
}
例3:HDU 3439
nの全配列の中でn-k転位列の方案数mod pだけあることを求めます.
データ範囲(0<=k<=n<=10^9,1<=p<=10^5,n!=0)
ターゲットの答えはC(n,k)*f(n-k)mod pで、f(i)はiです!全配列誤りスキーム数.
f(n−k)mod pは直接ループノードを探すことができ,ループノード長はm*2である.
次にC(n,k)mod pを要求する.本題ではpはまだ小さいが,pは素数を保証しないため,例2の方法は使用できない.
どうしようかな、まずp分解素因数p=p 1^k 1*p 2^k 2.*.*pt^kt.そしてc(n,k)mod pi^kiの値をそれぞれ求め,中国の残りの定理を用いて方程式群と同列に統合すればよい.
肝心なのはc(n,k)mod pi^kiを求めることだ.
c(n,k)=n!/(k!(n-k)!).分子分母中のすべての素因数piを蹴り出すと、分母はpi^kiと互いに質を交換し、拡張ユークリッドでmod pi^kiの逆元を求めればよい.
O(logPi)の時間内に分子分母に含まれる素因数piの個数を求め,それらの差をTとすることができる.
そしてすべての問題はnを求めることに転化します!素因数pを蹴ったmod p^cの値にpi^Tを乗じます.
前処理F[i]は、1からiまでpで割り切れない数の積を表す.P=p^cとする.1からnを求めてP個ごとに節を分けます.各グループの前p-1個の数mod Pは同余であるため、素因数pを含まない積はF[n%P]*F[P-1]^(n/P)であり、残りの素因数pの数に対してpを除去した後に積は(n/p)!に変換される.このようにして同じサブ問題に変換します!
1つのインスタンスでシミュレーション:20!Mod 9.
20!=(1*2*4*5*7*8*…*16*17*19)*(3*6*9*12*15*18)=(1*2*4*5*7*8)^2*(19*20)*3^6(1*2*3*4*5*6)=F[9-1]^(20/9)*F[20%9]*((20/3)! Mod 9).
ここで、F[0]からF[9−1]までの値は1,1,2,2,8,4,4,1,8である.
このように問題は完璧に解決された.
時間的複雑度はO(Σpi^ki).空間複雑度O(m).
code:
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cmath>
#include <string>
using namespace std;
typedef long long LL;
const int maxn=200003;
int ans,n,k,p,T,F[maxn],c[maxn],a[maxn],t[maxn],tot;
int get()
{
int f=0,v=0;char ch;
while(!isdigit(ch=getchar()))if(ch=='-')break;
if(ch=='-')f=1;else v=ch-48;
while(isdigit(ch=getchar()))v=v*10+ch-48;
if(f==1)return -v;else return v;
}
int calc()
{
F[0]=1,F[1]=0;
for(int i=2;i<2*p;i++)F[i]=LL(i-1)*(F[i-1]+F[i-2])%p;
return F[(n-k)%(2*p)];
}
void init(int p)
{
tot=0;
int pp=p;
for(int i=2;i*i<=pp;i++)
{
if(p%i!=0)continue;
a[++tot]=i,c[tot]=1,t[tot]=0;
while(p%i==0)c[tot]*=i,p/=i,t[tot]++;
}
if(p!=1)a[++tot]=p,c[tot]=p,t[tot]=1;
}
void ext_gcd(LL a,LL b,LL &x,LL &y)
{
if(b==0)
{
x=1,y=0;
return;
}
ext_gcd(b,a%b,x,y);
LL tp=x;
x=y;
y=tp-a/b*y;
}
LL Pow(LL a,LL b,int P)
{
LL res=1;
for(;b>0;b/=2)
{
if(b%2==1)res=res*a%p;
a=a*a%p;
}
return res;
}
LL calc(int n,int p,int _p)
{
if(n<_p)return F[n];
return (LL(F[n%p])*calc(n/_p,p,_p)%p*Pow(F[p-1],n/p,p)%p+p)%p;
}
int count(int n,int p)
{
if(n<p)return 0;
return n/p+count(n/p,p);
}
int C(int n,int k,int p,int _p,int t)
{
F[0]=1;
for(int i=1;i<=p;i++)
if(i%_p!=0)F[i]=LL(F[i-1])*i%p;else F[i]=F[i-1];
int t3=count(n,_p)-count(k,_p)-count(n-k,_p);
if(t3>=t)return 0;
LL t1=calc(n,p,_p),t2=calc(k,p,_p)*calc(n-k,p,_p)%p;
LL x,y;
ext_gcd(t2,p,x,y);
return (t1*x%p*Pow(_p,t3,p)+p)%p;
}
int work()
{
int res=0,tp,inv;
LL x,y;
for(int i=1;i<=tot;i++)
{
ext_gcd(p/c[i],c[i],x,y);
tp=C(n,k,c[i],a[i],t[i]);
res+=x*tp*(p/c[i])%p,res%=p;
}
return (res+p)%p;
}
int main()
{
T=get();
for(int i=1;i<=T;i++)
{
n=get(),k=get(),p=get();
ans=calc();
init(p);
ans=LL(ans)*work()%p;
printf("Case %d: %d
",i,ans);
}
return 0;
}