[BZOJ2313]分组

把问题拆成两问：

  1. 求出当m=1的时候的方案数ans
  2. 计算m^ans mod (10^9-401)

先考虑第二问。

显然第一问的ans很大很大，而答案只是要 mod (109-401)的值，而109-401刚好是一个质数。由费马小定理知，a^n mod p == a^(n mod (p-1)) mod p。因此我们要求出ans mod (p-1)，之后快速幂处理就好了。麻烦的就是如何求出ans mod (p-1)。

接下来考虑第一问。

考虑每一个盒子放i个球，那么就总共需要n/i个盒子，这种情况下的方案数就是： 
C(n, i) * C(n-i, i) * C(n-2i, i) * ... * C(i, i) / P(n/i, n/i)
根据排列组合的定义，可以化简为： 
n! / (i!^(n/i) * (n/i)!)

这就是ans的精确值。

现在我们要求的是ans mod (10^9-401-1)的值，就必须用中国剩余定理进行拆分了。 
x ≡ count[0] (mod A[0])
x ≡ count[1] (mod A[1])
x ≡ count[2] (mod A[2])
x ≡ count[3] (mod A[3])

其中count[i]是ans mod A[i]的值，A[] = {2, 13, 5281, 7283}。问题变成求满足上述同余方程组的x，得到count[]后用中国剩余定理处理（见主程序的中国剩余定理部分）。

棘手的是如何求出ans%A[i]

观察ans，我们把问题转成求x/y mod p（p为A[]中的某个质数）。 不妨令 
x = a * p^b
y = c * p^d

那么x/y mod p就要看b和d的关系了。由计算ans的公式可以发现，b >= d。 当b > d时，显然ans ≡ 0 (mod p) 当b == d时，ans ≡ a/c (mod p)。计算模上的除法一般转成乘上其乘法逆元，因此扩展欧几里德解之（见程序AdBModp）。

如何将n!转成a*p^b形式

注意到 
n! mod p =  1      * 2     * ... * (p-1)  * p
          * (p+1)  * (p+2) * ... * (2p-1) * 2p
          * ...
          * (kp+1) * (kp+2) * ... * (kp+m)  mod p
k就是[n/p]
m就是n-kp

         =  1 * 2 * ... * (p-1) * p
          * 1 * 2 * ... * (p-1) * 2p
          * ...
          * 1 * 2 * ... * m    mod p

         =  [(n%p)! * (p-1)!^(n/p)] * [(n/p)! * p^(n/p)]    mod p

假设n! = a * p^b 初始a = 1, b = 0 因为(n-n/p)! mod p、(p-1)!(n/p) mod p很小，可以预处理算出，为x。 那么a = ax, b = b+(n/p)。 然后递归处理(n/p)!最后就得到了n! = a * p^b

end

好了，基本所有的都讲完了，我觉得讲的有点像递归……为了得到那个同余方程组费了好大劲 ……到最后还要记得用中国剩余定理求解（具体在程序中有注释）

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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

/**************************************************************
    Problem: 2313
    User: abcdabcd987
    Language: C++
    Result: Accepted
    Time:144 ms
    Memory:1504 kb
****************************************************************/
 
#include <cmath>
#include <iostream>
using namespace std;
 
const int MODER = 1000000000-401;
const long long A[4] = {2, 13, 5281, 7283};
typedef pair<long long, long long> Pair;
int T, n, m;
long long X, Y, fac[4][7284];
inline long long sqr(const long long a) { return a*a; }
long long Power(const long long a, const long long x, const long long MODER)
{ /* a^x mod MODER */
  if (!x)
    return 1;
  long long res = sqr(Power(a, x/2, MODER))%MODER;
  if (x&1)
    res *= a;
  return res%MODER;
}
long long ExtendedGcd(const long long a, const long long b)
{
  if (!b)
  {
    X = 1, Y = 0;
    return a;
  }
  const long long gcd = ExtendedGcd(b, a%b);
  const long long t = X;
  X = Y;
  Y = t-a/b*Y;
  return gcd;
}
Pair FnModP(const long long n, const long long p)
{ /* n! mod p = p^ret.first * ret.second */
  if (n < A[p])
    return Pair(0, fac[p][n]);
  Pair x(FnModP(n/A[p], p));
  const long long ndp = n/A[p];
  x.first = x.first+ndp;
  x.second = (x.second*fac[p][n%A[p]])%A[p]*Power(fac[p][A[p]-1], ndp, A[p])%A[p];
  return x;
}
long long AdBModp(const Pair& fa, Pair tmp1, Pair tmp2, const int pown, const int p)
{ /* fa/(tmp1^pown+tmp2) mod p */
  if (fa.first > tmp1.first*pown+tmp2.first)
    return 0;
  tmp1.second = Power(tmp1.second, pown, A[p]);
  long long gcd = ExtendedGcd(tmp1.second*tmp2.second%A[p], A[p]);
  return X*fa.second/gcd%A[p];
}
int main()
{
  ios::sync_with_stdio(false);
 
  for (int i = 0; i < 4; ++i)
  {
    fac[i][0] = 1;
    for (int j = 1; j <= A[i]; ++j) fac[i][j] = fac[i][j-1]*j%A[i];
  }
 
  cin >> T;
  for (; T; --T)
  {
    cin >> n >> m;
    long long count[4] = {0};
    for (int p = 0; p < 4; ++p)
    {
      Pair fa(FnModP(n, p));
      const int sqrtn = static_cast<int>(sqrt(static_cast<double>(n)));
      /* To get the following remainder equations
        d === count[0] (mod A[0])
        d === count[1] (mod A[1])
        d === count[2] (mod A[2])
        d === count[3] (mod A[3])
      */
      for (int i = 1; i <= sqrtn; ++i)
        if (!(n%i))
        {
          const int ndi = n/i;
          Pair tmp1(FnModP(i, p));
          Pair tmp2(FnModP(ndi, p));
          count[p] = (count[p]+AdBModp(fa, tmp2, tmp1, i, p))%A[p];
          if (n-i*i)
            count[p] = (count[p]+AdBModp(fa, tmp1, tmp2, ndi, p))%A[p];
        }
      count[p] = (count[p]+A[p])%A[p];
    }
    /*  Chinese remainder theorem
                                     d === last.first               (mod last.second)
                                     d === count[i]                 (mod A[i]       )
Let                          d = last.second*k+last.first
Then I have   last.second*k+last.first === count[i]                 (mod A[i]       )
Let   last.second*k+last.first = A[i]*Y+count[i]
Use Extended GCD to solve (k, Y)
The new remainder equation is        d === last.second*k+last.first (mod LCM(last.second, A[i]))
    */
    Pair last(count[0], A[0]);
    for (int i = 1; i < 4; ++i)
    {
      long long gcd = ExtendedGcd(last.second, A[i]);
      long long tmp = last.second*A[i]/ExtendedGcd(last.second, A[i]);
      long long k = X*(count[i]-last.first)/gcd%tmp;
      last.first = (k*last.second+last.first)%tmp;
      last.second = tmp;
    }
    /* (a^b) mod p = (a^(b mod (p-1))) mod p*/
    cout << Power(m, last.first, MODER) << endl;
  }
}

Comments