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])   用的是求线性模方程组的方法，具体百度

_______________________________________________

Program Stone;
const di:array[1..4]of longint=(2,13,5281,7283);
      prime=999999599;
type data=record pow,fac:int64;end;
var i,j,t:longint;
    n,m,ans:int64;
    aj:array[1..4,0..8000]of int64;
    s_up,s_low:array[1..4]of data;
    ci:array[1..4]of int64;
  function pow(c,x,p:int64):int64;  //求 c^x mod p 的值
  var i,k:int64;
   begin
     pow:=1;
     i:=x;
     k:=c;
     while i>0 do
      begin
        if i mod 2=1 then pow:=k*pow mod p;
        i:=i shr 1;
        k:=k*k mod p;
      end;
   end;

  function factorial(p,a:int64):data;  //求 p! mod di[a] 的值，存储成 fac*di[a]^pow
  var i,j:int64;
      s:data;
   begin
      i:=p div di[a];j:=p mod di[a];
      s.fac:=1;s.pow:=0;
      if i>0 then s:=factorial(i,a);
      if i mod 2=1 then s.fac:=-s.fac;
      s.fac:=(s.fac*aj[a,j]) mod di[a];
      s.pow:=(s.pow+i);
      factorial:=s;
   end;
  function exgcd(a,b:int64;var x,y:int64):int64;   //扩展欧几里德

  var t:int64;
   begin
     if b=0 then begin x:=1;y:=0;exit(a);end;
     exgcd:=exgcd(b,a mod b,x,y);
     t:=x;
     x:=y;
     y:=t-(a div b)*y;
   end;
  function work(a,b,n:int64):int64;   //求解  ax=n (mod b)的x 值
  var y,g:int64;
   begin
     g:=exgcd(a,b,work,y);
     work:=((work*(n div g) mod (b div g)+(b div g))mod (b div g));  //求出的解应该为最小整数解
   end;
  function lcm(a,b:int64):int64;   //求最小公倍数
  var i,j,k,gcd:int64;
   begin
     i:=a;j:=b;
     while j<>0 do
      begin
        k:=i;i:=j;j:=k mod j;
      end;
     lcm:=a*b div i;
   end;
  Procedure doans(i:int64);   //对于将n分成每组有i 个的方案数
  var j:longint;
      ai,yi,xi,k:int64;
      s:data;
   begin
     for j:=1 to 4 do    //s_up表示分子，s_low表示分母
        begin
         s_low[j]:=factorial(n div i,j);
         s:=factorial(i,j);
         s_low[j].pow:=s_low[j].pow+s.pow*(n div i);
         s_low[j].fac:=s_low[j].fac*pow(s.fac,n div i,di[j]);
         if s_up[j].pow<>s_low[j].pow then ci[j]:=0
                                      else ci[j]:=work(s_low[j].fac,di[j],s_up[j].fac);
        end;
        ai:=di[1];yi:=ci[1];    //合并ans =xi (mod di[i])   
        for j:=2 to 4 do
        begin
          xi:=work(ai,di[j],(ci[j]-yi));
          yi:=ai*xi+yi;
          ai:=lcm(ai,di[j]);
        end;
        ans:=(ans+yi)mod (prime-1);
   end;
  Procedure main;
  var i,j:longint;
      s:data;
   begin
      if n>1 then ans:=2 else ans:=1;
      for i:=1 to 4 do
        s_up[i]:=factorial(n,i);
      for i:=2 to round(sqrt(n)) do    //枚举n的约数
       if n mod i=0 then
         begin
           doans(i);
           if n div i<>i then doans(n div i);
         end;
   end;
Begin
 assign(input,'input.in');reset(input);
  readln(t);
  for i:=1 to 4 do
   begin
     aj[i,0]:=1;
     for j:=1 to di[i] do
       aj[i,j]:=aj[i,j-1]*j mod di[i];    //预处理 i! mod di[i]的值 (i<=di)
   end;
  for i:=1 to t do
   begin
     readln(n,m);
     main;
     writeln(pow(m,ans,prime));
   end;
end.