计算系数(NOIP2011 day2)

给定一个多项式 (ax + by)^k 。

输入a,b,k,n,m。

请求出多项式展开后 xⁿy^m 项的系数。

【数据范围】
对于 30%的数据，有 0≤k≤10；
对于 50%的数据，有 a = 1，b = 1；
对于 100%的数据，有 0≤k≤1,000，0≤n, m≤k，且 n + m = k，0≤a，b≤1,000,000。

分析：

根据二项式定理可知：xⁿy^m的系数是C(k,n)*aⁿ*b^m。

关键是求C(k,n)的值：

方法一：

本题中k(<=1000)的值比较小，可以利用递推式C(k,n)=C(k-1,n-1)+C(k-1,n)直接求得。

const
  maxn=1000;
  md=10007;
var
  f:array[0..maxn,0..maxn] of longint;
  a,b,n,k,m,s,i,j:longint;
begin
  assign(input,'factor.in'); reset(input);
  assign(output,'factor.out'); rewrite(output);
  readln(a,b,k,n,m);
  a:=a mod md;
  b:=b mod md;
  f[0,0]:=1;
  for i:=1 to k do begin f[i,0]:=1; f[i,i]:=1; end;
  for i:=2 to k do
    for j:=1 to i-1 do f[i,j]:=(f[i-1,j]+f[i-1,j-1])mod md;
  s:=1;
  for i:=1 to n do s:=(s*a)mod md;
  for i:=1 to m do s:=(s*b)mod md;
  writeln((f[k,n]*s)mod md);
end.

方法二：

如果k的范围再大的话，如k<=50000，就无法直接使用方法一中的递推式。

另一种解决问题的方法是：根据组合公式C(k,n)=(k*(k-1)*...*(k-n+1))/(1*2*...*n)

分子和分母分解质因数后约分即可（当然要先求出k以内的素数）。

const
  maxk=100000;
  md=10007;
  fin='factor.in';
  fout='factor.out';
var
  p,e:array[1.. maxk] of longint;
  f:array[1..maxk] of boolean;
  a,b,k,n,m,cnt,s:longint;
  procedure init;
    begin
      assign(input,fin); reset(input);
      readln(a,b,k,n,m);
      a:=a mod md;
      b:=b mod md;
      close(input);
    end;
  procedure make_p;//生成素数表
    var i,j:longint;
    begin
      fillchar(f,sizeof(f),true);
      cnt:=0;
      for i:=2 to k do
        if f[i] then
          for j:=2 to k div i do f[i*j]:=false;
      for i:=2 to k do
        if f[i] then
          begin inc(cnt); p[cnt]:=i; end;
    end;
  procedure make_e; //生成指数
    var i,j,now:longint;
    begin
      fillchar(e,sizeof(e),0);
      for i:=k-n+1 to k do  //处理分子指数
        begin
          j:=i; now:=1;
          while j>1 do
            begin
              while j mod p[now]=0 do
                begin inc(e[now]); j:=j div p[now]; end;
              inc(now);
            end;
        end;
      for i:=1 to n do     //处理分母指数
        begin
          j:=i; now:=1;
          while j>1 do
            begin
              while j mod p[now]=0 do
                begin dec(e[now]); j:=j div p[now]; end;
              inc(now);
            end;
        end;
      end;
  procedure calc;
    var i,j:longint;
    begin
      s:=1;
      for i:=1 to cnt do
        if e[i]>0 then
          for j:=1 to e[i] do s:=(s*p[i])mod md;
      for i:=1 to n do s:=(s*a)mod md;
      for i:=1 to m do s:=(s*b)mod md;
    end;
  procedure print;
    begin
      assign(output,fout); rewrite(output);
      writeln(s);
      close(output);
    end;
begin
  init;
  make_p;
  make_e;
  calc;
  print;
end.