0%

组合数取模

Long long ago, there was a beautiful boy, he was so lonely in the summer night. At his most lonely moment, he decided to invite his boyfriends to his bedroom, and they will be happy together. The beautiful boy owned n boyfriends,and he would be happy only with exactly m boyfriends in his bedroom. So how many ways he could be happy? As the beautiful boy had too many boyfriends, and he was so lonely, the answer could be very large. You just need give the answer module 10^18+7. Input There are multiple test cases(No more than 66). For each test case: There are two integers n,m (1<=n,m<=10^5) The last line contains two integers: 0 0 Output For each cases,output a line with a integer ,which is the answer module 1e18+7. If the beautiful could never be happy,you can just output his voice from heart: "How lonely night~"(without quotes) Sample Input 3 2 10 5 3 5 0 0 Sample Output 3 252 How lonely night~ Hint For the first case: The beautiful boy had 3 friends :A,B,C He could invite A+B,A+C,B+C. So the answer should be 3 For the third case: The beautiful boy had 3 friends, but he wanted 5 friends,so he would never be happy.(So pity!) '~' is shift+`

之前一直不会乘法逆元,不过现在会了... 如果a和n互质,ax+ny=1,用扩展欧几里得算法求出的x就是a对于n的乘法逆元,即ax同余1(mod n) 证明比较显然... 然后计算组合数过程中的除法就变成了乘法... 但是这道题的中间结果爆long long,用一个快速乘的技巧即可

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#include<bits/stdc++.h>
typedef long long LL;
using namespace std;
LL t[100010]={};
LL cheng(LL a,LL b,LL p)
{
    if(b==1) return a%p;
    LL e;
    e=cheng(a,b/2,p);
    e=e+e;
    e%=p;
    if(b%2==1) e+=a;
    e%=p;
    return e;
}
void ex_gcd(LL a,LL b,LL &x,LL &y,LL &d)
{
if(!b)
{
d=a;
x=1;
y=0;
}
else
{
ex_gcd(b,a%b,y,x,d);
y-=x*(a/b);
}
}
LL inv(LL t,LL p)
{
LL d,x,y;
ex_gcd(t,p,x,y,d);
return d==1?(x%p+p)%p:-1;
}
LL C(LL n,LL m,LL p)
{
m=min(m,n-m);
if(m==0) return 1;
LL ans=1;
for(LL i=n;i>=n-m+1;i--)
{
ans=cheng(ans,i,p);
}
ans=cheng(ans,t[m],p);
return ans;
}
int main()
{
LL p;
p=1e18;
p+=7;
LL a,b;
t[0]=1;
for(int i=1;i<=100000;i++)
{
t[i]=cheng(t[i-1],inv(i,p),p);
}
while(1)
{
scanf("%lld%lld",&a,&b);
if(a==b&&a==0)
{
return 0;
}
if(a<b)
{
printf("How lonely night~\n");
}
else
{
printf("%lld\n",C(a,b,p));
}
}
return 0;
}

当然,逆元也可以不使用扩展欧几里得,也可以使用递推法 inv(a) = (p - p / a) * inv(p % a) % p 这个方法理论上是O(n)的,时间复杂度比扩展欧几里得好一些

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
LL inv[100010]={};
LL cheng(LL a,LL b,LL p)
{
    if(b==1) return a%p;
    LL e;
    e=cheng(a,b/2,p);
    e=e+e;
    e%=p;
    if(b%2==1) e+=a;
    e%=p;
    return e;
}
LL C(LL n,LL m,LL p)
{
m=min(m,n-m);
if(m==0) return 1;
LL ans=1;
for(LL i=n;i>=n-m+1;i--)
{
ans=cheng(ans,i,p);
}
ans=cheng(ans,inv[m],p);
return ans;
}
int main()
{
LL p;
p=1e18;
p+=7;
LL a,b;
inv[0]=inv[1]=1;
for(int i=2;i<=100000;i++) inv[i]=cheng(p-p/i,inv[p%i],p);
    for(int i=2;i<=100000;i++) inv[i]=cheng(inv[i],inv[i-1],p);
while(1)
{
scanf("%lld%lld",&a,&b);
if(a==b&&a==0)
{
return 0;
}
if(a<b)
{
printf("How lonely night~\n");
}
else
{
printf("%lld\n",C(a,b,p));
}
}
}

http://www.lydsy.com/JudgeOnline/problem.php?id=2982

Description

LMZ有n个不同的基友,他每天晚上要选m个进行河蟹,而且要求每天晚上的选择都不一样。那么LMZ能够持续多少个这样的夜晚呢?当然,LMZ的一年有10007天,所以他想知道答案mod 10007的值。(1<=m<=n<=200,000,000) Input

第一行一个整数t,表示有t组数据。(t<=200) 接下来t行每行两个整数n, m,如题意。 Output

T行,每行一个数,为C(n, m) mod 10007的答案。

这道题使用的方法是Lucas定理,然后逆元是用费马小定理求的...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
LL n,m,p;
LL quick_mod(LL a,LL b)
{
LL ans=1;
a%=p;
while(b)
{
if(b&1)
{
ans=ans*a%p;
b--;
}
b>>=1;
a=a*a%p;
}
return ans;
}
LL C(LL n,LL m)
{
if(m>n) return 0;
LL ans=1;
for(int i=1;i<=m;i++)
{
LL a=(n+i-m)%p;
LL b=i%p;
ans=ans*(a*quick_mod(b,p-2)%p)%p;
}
return ans;
}
LL Lucas(LL n,LL m)
{
if(m==0) return 1;
return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}
int main()
{
int T;
scanf("%d",&T);
p=10007;
while(T--)
{
scanf("%lld%lld",&n,&m);
printf("%lld\n",Lucas(n,m));
}
return 0;
}

组合数范围较小时,还有传统的杨辉三角求法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
#include<cstdio>
const int N = 2000 + 5;
const int MOD = (int)1e9 + 7;
int comb[N][N];//comb[n][m]就是C(n,m)
void init(){
    for(int i = 0; i < N; i ++){
        comb[i][0] = comb[i][i] = 1;
        for(int j = 1; j < i; j ++){
            comb[i][j] = comb[i-1][j] + comb[i-1][j-1];
            comb[i][j] %= MOD;
        }
    }
}
int main(){
    init();
}

以上就是一些基本的组合数取模的方法,但是如果我们取模的数和其他数不互质怎么办? 在组合数范围比较小时(小于1000000?) 我们可以用一个数组a[1000000] a[i]代表我们最后结果需要乘a[i]个i,因为最后组合数的结果为整数,所以我们把所有的i进行质因数分解,然后把a[i]分配到对应的质因子中,最后我们就可以保证a数组中所有的值都为非负数,即我们进行的还是乘法运算,不会有精度问题,就是时间复杂度和空间复杂度比较大,但是还可以接受

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#include<bits/stdc++.h>
#define LL unsigned long long
using namespace std;
int n,m;
const LL mod=1000000000000000007ll;
LL a[100011];
LL qadd(LL a,LL t)
{
if (t==1) return a;
LL tem;
tem=qadd(a,t/2);
tem=tem+tem;
tem%=mod;
if (t%2==1)
{
tem+=a;
tem%=mod;
}
return tem;
}
LL qpow(LL a,LL t)
{
if (t==1) return a;
LL tem;
tem=qpow(a,t/2);
tem=qadd(tem,tem);
if (t%2==1) tem=qadd(tem,a);
return tem;
}
void work()
{
memset(a,0ll,sizeof(a));
for (int i=2;i<=n;i++) a[i]++;
for (int i=2;i<=m;i++) a[i]--;
for (int i=2;i<=n-m;i++) a[i]--;

for (int i=n;i>1;i--)
{
if (a[i]!=0)
for (int j=2;j<=sqrt(i);j++)
{
if (i%j==0)
{
a[j]+=a[i];
a[i/j]+=a[i];
a[i]=0;
break;
}
}
}
//for (int i=1;i<=n;i++) cout<<a[i]<<" ";
//cout<<endl;
LL ans=1ll;
for (int i=2;i<=n;i++)
{
if (a[i]!=0)
ans=qadd(ans,qpow(i,a[i]));
ans%=mod;
}
cout<<ans<<endl;
}
int main()
{
int T;
while (cin>>n>>m,n!=0&&m!=0)
{
if (n<m)
{
cout<<"How lonely night~"<<endl;
}
else
{
work();
}
}

}