组合数取模

D - Combinations
Given n different objects, you want to take k of them. How many ways to can do it?

For example, say there are 4 items; you want to take 2 of them. So, you can do it 6 ways.

Take 1, 2

Take 1, 3

Take 1, 4

Take 2, 3

Take 2, 4

Take 3, 4

Input
Input starts with an integer T (≤ 2000), denoting the number of test cases.

Each test case contains two integers n (1 ≤ n ≤ 106), k (0 ≤ k ≤ n).

Output
For each case, output the case number and the desired value. Since the result can be very large, you have to print the result modulo 1000003.

Sample Input
3

4 2

5 0

6 4

Sample Output
Case 1: 6

Case 2: 1

Case 3: 15
这题目有两种方法

第一种是通俗易懂但要复杂一点的方法，就是你要分解质因数，用桶排记录它们的个数，组合的公式 C（n，m）= n！/（m！*（n-m）！）将n! 的所有质因数累加，m！和（n-m)！的质因数减去，

最后算出这些质因数相乘的值。这就要用到线性筛啊，快速幂等等

附上代码（很长)

  1 var
  2   n,m,i,j,k,t,sum,p:longint;
  3   ans:int64;
  4   mark:array[1..200000] of boolean;
  5   prim,num,markx:array[1..200000] of int64;
  6 function prime(x:longint):longint;
  7   var i,j,k:longint;
  8   begin
  9     k:=0;
 10     for i:=2 to x do
 11       begin
 12         if mark[i]
 13           then
 14             begin
 15               k:=k+1;
 16               prim[k]:=i;
 17               markx[i]:=k;
 18             end;
 19         for j:=1 to k do
 20           begin
 21             if i*prim[j]>m
 22               then break;
 23             mark[i*prim[j]]:=false;
 24             if i mod prim[j]=0
 25               then break;
 26           end;
 27       end;
 28     prime:=k;
 29   end;
 30 function f(a,b:longint):int64;
 31   var y:int64;
 32   begin
 33     t:=1;
 34     y:=a;
 35     while b<>0 do
 36       begin
 37         if (b and 1)=1
 38           then t:=t*y mod p;
 39         y:=y*y mod p;
 40         b:=b shr 1;
 41       end;
 42     f:=t;
 43   end;
 44  
 45 begin
 46   readln(m,n,p);
 47   for i:=1 to m do
 48     mark[i]:=true;
 49   sum:=prime(m);
 50   for i:=2 to m do
 51     begin
 52       if mark[i]
 53         then
 54           begin
 55             inc(num[markx[i]]);
 56             continue;
 57           end;
 58       k:=i;
 59       j:=1;
 60       while (prim[j]<=trunc(sqrt(i)))and(k>1) do
 61         begin
 62           if k mod prim[j]=0
 63             then
 64               while k mod prim[j]=0 do
 65                 begin
 66                   k:=k div prim[j];
 67                   num[j]:=num[j]+1;
 68                 end;
 69           j:=j+1;
 70         end;
 71       if k>1
 72         then inc(num[markx[k]]);
 73     end;
 74   for i:=2 to n do
 75     begin
 76       if mark[i]
 77         then
 78           begin
 79             dec(num[markx[i]]);
 80             continue;
 81           end;
 82       k:=i;
 83       j:=1;
 84       while (prim[j]<=trunc(sqrt(i)))and(k>1) do
 85         begin
 86           if k mod prim[j]=0
 87             then
 88               while k mod prim[j]=0 do
 89                 begin
 90                   k:=k div prim[j];
 91                   num[j]:=num[j]-1;
 92                 end;
 93           j:=j+1;
 94         end;
 95       if k>1
 96         then dec(num[markx[k]]);
 97     end;
 98       for i:=2 to m-n do
 99     begin
100       if mark[i]
101         then
102           begin
103             dec(num[markx[i]]);
104             continue;
105           end;
106       k:=i;
107       j:=1;
108       while (prim[j]<=trunc(sqrt(i)))and(k>1) do
109         begin
110           if k mod prim[j]=0
111             then
112               while k mod prim[j]=0 do
113                 begin
114                   k:=k div prim[j];
115                   num[j]:=num[j]-1;
116                 end;
117           j:=j+1;
118         end;
119       if k>1
120         then dec(num[markx[k]]);
121     end;
122  
123   ans:=1;
124   for i:=1 to sum do
125      if num[i]>0
126        then ans:=(ans*f(prim[i],num[i])) mod p;
127   writeln(ans);
128 end.

View Code

第二种就要简单多了，直接算，但是因为要取模，模数又很大，不能直接取模，要用到逆元，就是乘上逆元再模。求逆元呢，要用到扩展欧几里得。再就是一点组合数的运算了

 1 var
 2 l,t,i,modd:longint;
 3 f,inv:array[0..1100005]of int64;
 4 n,m:int64;
 5 function exgcd(a,b:int64;var x,y:int64):int64;
 6  var r,t:int64;
 7  begin
 8    if b=0 then
 9     begin
10       x:=1;
11       y:=0;
12       exit(a);
13     end;
14     r:=exgcd(b,a mod b,x,y);
15     t:=x;
16     x:=y;
17     y:=t-a div b*y;
18     exit(r);
19  end;
20 function cal(d,m:int64):int64;
21  var x,y:int64;
22   begin
23    exgcd(d,m,x,y);
24    exit((x mod m+m) mod m);
25  end;
26  function c(a,b:int64):int64;
27   begin
28    exit(f[a]*inv[b]mod modd *inv[a-b] mod modd);
29   end;
30 begin
31 readln(t);
32 modd:=1000003;
33 f[0]:=1;inv[0]:=1;
34 for i:=1 to 1010000 do
35  begin
36   f[i]:=f[i-1]*i mod modd;
37   inv[i]:=cal(f[i],modd);
38  end;
39 
40 for l:=1 to t do
41  begin
42   readln(n,m);
43 //  if m=0 then begin writeln(1); continue; end; 不要想当然乱加条件，先试试程序对于特殊样例能不能过，能过就不要加初始条件了
44   writeln('Case ',l,': ',c(n,m));
45  end;
46 end.

NOIP2018 rp++