bzoj 2313 分组 2012-01-15

http://www.zybbs.org/JudgeOnline/problem.php?id=2313

_______________________________________________

对于n要求方案数：

对于n的约数s，即将n个数分为s组，每组有n/s=p个数共有：

ans (p)=(C(p,n)*C(p,n-p)*C(p,n-2p)*...*C(p,p))/(s!)

=(n!*（n-p)!*...*p!)/((p!)^s*(n-p)!*(n-2p)!..1!*s!)

=(n!)/((p!)^s*s!)

种方案

所以总共的方案数就有 total=sigma( ans(si))

题目要求的是  m^total mod (10^9-401)

设prime=10^9-401

prime-1 分解质因数为 di[i]=｛2，13，5281，7283｝

首先需要知道：①费尔马小定理： a^(p-1)=1(mod p)

②  （k-1)! =-1 (mod k)

所以 m^(prime-1)=1 （1 mod prime);

m^ans=m^(ans mod (prime-1)) (mod prime);

所以现在需要求 ans mod (prime-1)

等价于求满足 ans =xi (mod di[i])  (1<=i<=4)  最后将每一个yi 用中国剩余定理合并即求出 x=ans( mod (prime-1))。

最后即可算出每一个x的和即为 total mod (prime-1)

对于 ans=(n!)/((p!)^s*s!) mod di

怎么求 n! mod di 呢？

因为 n!=1*2..(di-1)*di*(di+1)*..*(q*di)*(q*di+1)...*n （mod di)

=((di-1)! mod di)^q*(n-q*di)!*q!*(di)^q ( mod di)

=(-1)^q*(n-q*di)!*q!*di^q ( mod di)

因为（n-q*di)不是特别大（n-q*di)! mod di的值可以预处理出来

所以可以用递归求一个数的阶乘模di的值   求出来的值可以化简为 fac*di^pow

将分数ans 的分子存成  di^pow1*fac1

分母存成  di^pow2*fac2

易证pow1>=pow2 所以当 pow1>pow2时  ans= (fac1/fac2) (mod di);

当 pow1=pow2时  ans= 0 (mod di)

！！！程序合并ans =xi (mod di[i])   用的是求线性模方程组的方法，具体百度

_______________________________________________

  1 Program Stone;
  2 const di:array[1..4]of longint=(2,13,5281,7283);
  3       prime=999999599;
  4 type data=record pow,fac:int64;end;
  5 var i,j,t:longint;
  6     n,m,ans:int64;
  7     aj:array[1..4,0..8000]of int64;
  8     s_up,s_low:array[1..4]of data;
  9     ci:array[1..4]of int64;
 10   function pow(c,x,p:int64):int64;  //求 c^x mod p 的值
 11   var i,k:int64;
 12    begin
 13      pow:=1;
 14      i:=x;
 15      k:=c;
 16      while i>0 do
 17       begin
 18         if i mod 2=1 then pow:=k*pow mod p;
 19         i:=i shr 1;
 20         k:=k*k mod p;
 21       end;
 22    end;
 23
 24   function factorial(p,a:int64):data;  //求 p! mod di[a] 的值，存储成 fac*di[a]^pow
 25   var i,j:int64;
 26       s:data;
 27    begin
 28       i:=p div di[a];j:=p mod di[a];
 29       s.fac:=1;s.pow:=0;
 30       if i>0 then s:=factorial(i,a);
 31       if i mod 2=1 then s.fac:=-s.fac;
 32       s.fac:=(s.fac*aj[a,j]) mod di[a];
 33       s.pow:=(s.pow+i);
 34       factorial:=s;
 35    end;
 36   function exgcd(a,b:int64;var x,y:int64):int64;   //扩展欧几里德
 37
 38   var t:int64;
 39    begin
 40      if b=0 then begin x:=1;y:=0;exit(a);end;
 41      exgcd:=exgcd(b,a mod b,x,y);
 42      t:=x;
 43      x:=y;
 44      y:=t-(a div b)*y;
 45    end;
 46   function work(a,b,n:int64):int64;   //求解  ax=n (mod b)的x 值
 47   var y,g:int64;
 48    begin
 49      g:=exgcd(a,b,work,y);
 50      work:=((work*(n div g) mod (b div g)+(b div g))mod (b div g));  //求出的解应该为最小整数解
 51    end;
 52   function lcm(a,b:int64):int64;   //求最小公倍数
 53   var i,j,k,gcd:int64;
 54    begin
 55      i:=a;j:=b;
 56      while j<>0 do
 57       begin
 58         k:=i;i:=j;j:=k mod j;
 59       end;
 60      lcm:=a*b div i;
 61    end;
 62   Procedure doans(i:int64);   //对于将n分成每组有i 个的方案数
 63   var j:longint;
 64       ai,yi,xi,k:int64;
 65       s:data;
 66    begin
 67      for j:=1 to 4 do    //s_up表示分子，s_low表示分母
 68         begin
 69          s_low[j]:=factorial(n div i,j);
 70          s:=factorial(i,j);
 71          s_low[j].pow:=s_low[j].pow+s.pow*(n div i);
 72          s_low[j].fac:=s_low[j].fac*pow(s.fac,n div i,di[j]);
 73          if s_up[j].pow<>s_low[j].pow then ci[j]:=0
 74                                       else ci[j]:=work(s_low[j].fac,di[j],s_up[j].fac);
 75         end;
 76         ai:=di[1];yi:=ci[1];    //合并ans =xi (mod di[i])
 77         for j:=2 to 4 do
 78         begin
 79           xi:=work(ai,di[j],(ci[j]-yi));
 80           yi:=ai*xi+yi;
 81           ai:=lcm(ai,di[j]);
 82         end;
 83         ans:=(ans+yi)mod (prime-1);
 84    end;
 85   Procedure main;
 86   var i,j:longint;
 87       s:data;
 88    begin
 89       if n>1 then ans:=2 else ans:=1;
 90       for i:=1 to 4 do
 91         s_up[i]:=factorial(n,i);
 92       for i:=2 to round(sqrt(n)) do    //枚举n的约数
 93        if n mod i=0 then
 94          begin
 95            doans(i);
 96            if n div i<>i then doans(n div i);
 97          end;
 98    end;
 99 Begin
100  assign(input,‘input.in‘);reset(input);
101   readln(t);
102   for i:=1 to 4 do
103    begin
104      aj[i,0]:=1;
105      for j:=1 to di[i] do
106        aj[i,j]:=aj[i,j-1]*j mod di[i];    //预处理 i! mod di[i]的值 (i<=di)
107    end;
108   for i:=1 to t do
109    begin
110      readln(n,m);
111      main;
112      writeln(pow(m,ans,prime));
113    end;
114 end.
115
116