題目描述

小 X 駕駛著他的飛船準備穿梭過一個 \(n\) 維空間，這個空間裏每個點的坐標可以用 \(n\) 個實數表示，即 \((x_1,x_2,\dots,x_n)\)。

為了穿過這個空間，小 X 需要在這個空間中選取 \(c\)（\(c\geq 2\)）個點作為飛船停留的地方，而這些點需要滿足以下三個條件：

每個點的每一維坐標均為正整數，且第 \(i\) 維坐標不超過 \(m_i\)。
第 \(i+1\)（\(1\leq i<c\)）個點的第 \(j\)（\(1\leq j\leq n\)）維坐標必須嚴格大於第 \(i\) 個點的第 \(j\) 維坐標。
存在一條直線經過所選的所有點。在這個 \(n\)

輸入格式

第一行包含一個正整數 \(T\)，表示有 \(T\) 組數據需要求解。

每組數據包含兩行，第一行包含兩個正整數 \(n,c\)（\(c\geq 2\)

第二行包含 \(n\) 個正整數，依次表示 \(m_1,m_2,…,m_n\)。

輸出格式

共 \(T\) 行，每行包含一個非負整數，依次對應每組數據的答案。

樣例

input

3
2 3
3 4
3 3
3 4 4
4 4
5 9 7 8

output

2
4
846

explanation

對於第一組數據，共有兩種可行的方案：一種選擇 \((1,1),(2,2),(3,3)\)，另一種選擇 \((1,2),(2,3),(3,4)\)。

限制與約定

\[T \leq 1000, n \leq 11, c \leq 20, m_i \leq 10^5\]

時間限制: 1s, 空間限制: 512MB.

題解

首先，我們發現只要滿足\(x_i < y_i, i = 1,2,\dots,n\) ，（其中\((x_1, x_2,\dots,x_n),(y_1, y_2,\dots,y_n)\)分別是第一個點和最後一個點的坐標）且點的不同排列算同一種即可。

那麽，我們枚舉這兩個點之後，以其為端點的線段上共有（除去兩端）\(\gcd(y_1-x_1,y_2-x_2,\dots.y_n-x_n)\)個點。其中選出\(c-2\)個點的方式有\(C_{\gcd(y_1-x_1,y_2-x_2,\dots.y_n-x_n)}^{c-2}\)種，所以總方案數是

\[ans=\sum_{\substack{1\leq x_i < y_i\leq m_i \\ i = 1,2,\dots,n}}C_{\gcd(y_1-x_1,y_2-x_2,\dots,y_n-x_n)}^{c-2}\]

我們發現，所有\(p_i=y_i-x_i\)固定之後，方案數恰有

\[C_{\gcd(p_1, p_2,\dots,p_n)}^{c-2}\prod_{i=1}^n(m_i-p_i)\]

種，這是因為\(p_i\)固定後\(x_i\)只能取\(1,2,\dots,m_i-p_i\).

於是我們有

\[ans=\sum_{\substack{1\leq q_i < m_i \\ i = 1,2,\dots,n}}C_{\gcd(p_1, p_2,\dots,p_n)}^{c-2}\prod_{i=1}^n(m_i-p_i)\]

優先枚舉\(d=\gcd(p_1, p_2,\dots,p_n)\)（下式中\(m_{min}\)表示\(\min(m_1,m_2,\dots,m_n)\):

\[ans=\sum_{d=1}^{m_{min}-1}C_d^{c-2}\sum_{\substack{1\leq q_i < m_i \\ i = 1,2,\dots,n\\ \gcd(p_1, p_2,\dots,p_n)=d}}\prod_{i=1}^n(m_i-p_i)\]

令

\[f(d)=\sum_{\substack{1\leq q_i < m_i \\ i = 1,2,\dots,n\\ \gcd(p_1, p_2,\dots,p_n)=d}}\prod_{i=1}^n(m_i-p_i)\]

則有

\[\begin{aligned}\sum_{d\mid l}f(l)&=\sum_{\substack{1\leq q_i < m_i \\ i = 1,2,\dots,n\\ d\mid\gcd(p_1, p_2,\dots,p_n)}}\prod_{i=1}^n(m_i-p_i)\\&=\sum_{\substack{1\leq q_i‘d < m_i \\ i = 1,2,\dots,n}}\prod_{i=1}^n(m_i-p_i‘d)\\&=\prod_{i=1}^n\sum_{p‘=1}^{\left\lfloor\frac{m_i-1}d\right\rfloor}(m_i-p‘d)\\&=\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}d\right\rfloor-\frac{d\left\lfloor\frac{m_i-1}d\right\rfloor\left(1+\left\lfloor\frac{m_i-1}d\right\rfloor\right)}2\right]\end{aligned}\]

最後一步是等差數列求和。

由莫比烏斯反演可知：

\[f(d)=\sum_{d\mid l}\mu\left(\frac ld \right)\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}d\right\rfloor-\frac{d\left\lfloor\frac{m_i-1}d\right\rfloor\left(1+\left\lfloor\frac{m_i-1}d\right\rfloor\right)}2\right]\]

於是有

\[\begin{aligned}ans&=\sum_{d=1}^{m_{min}-1}C_d^{c-2}\sum_{d\mid l}\mu\left(\frac ld \right)\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]\\&=\sum_{l=1}^{m_{min}-1}\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]\sum_{d\mid l}\mu\left(\frac ld \right)C_d^{c-2}\end{aligned}\]

易知\(\prod_{i=1}^n\left[m_i\left\lfloor\frac{m_i-1}l\right\rfloor-\frac{l\left\lfloor\frac{m_i-1}l\right\rfloor\left(1+\left\lfloor\frac{m_i-1}l\right\rfloor\right)}2\right]\)在將所有\(\left\lfloor\frac{m_i-1}l\right\rfloor\)作為常數看待的情況下是一個關於\(l\)的不超過\(n\)次（是因為\(mod10007\)之後有可能次數更低）的多項式，而且引起至少一個\(\left\lfloor\frac{m_i-1}l\right\rfloor\)改變的\(l\)有\(O(n\sqrt m)\)個，於是可以對每一個所有\(\left\lfloor\frac{m_i-1}l\right\rfloor\)都不改變的段分別求和。由於此時它是關於\(l\)的多項式，我們只需預處理出所有的\(S_{c,p,t}=\sum_{l=1}^tl^p\sum_{d\mid l}\mu\left(\frac ld \right)C_d^{c-2}\)即可。至於如何在只改變一個因式時快速維護成績中各項的系數，可以用線段樹實現（代碼中為ckw線段樹）。

完。

Code

#include <algorithm>
#include <cstdio>
const int N = 12;
const int C = 21;
const int M = 100001;
const int mod = 10007;
int fac[mod], inv[mod];
int mu[M], f[C][M], S[C][M][N];
int p[N * 1000], m[N];
inline int CC(int a, int b) {
  if (!b) return 1;
  if (a % mod < b % mod) return 0;
  return CC(a / mod, b / mod) * fac[a % mod] % mod * inv[fac[b % mod]] % mod * inv[fac[(a - b) % mod]] % mod;
}
int a[64][N], po[64];
int mdp[N];
int n, mm, c;
void upd(int x) {
  int l = x * 2, r = x * 2 + 1;
  po[x] = po[l] + po[r];
  for (int i = 0; i <= po[x]; ++i) a[x][i] = 0;
  for (int i = 0; i <= po[l]; ++i)
    for (int j = 0; j <= po[r]; ++j)
      a[x][i + j] = (a[x][i + j] + a[l][i] * a[r][j] % mod) % mod;
}
int sum(int l, int r) {
  int ans = 0;
  for (int i = 0; i <= po[1]; ++i)
    ans = (ans + a[1][i] * (S[c][r][i] - S[c][l - 1][i]) % mod) % mod;
  return ans;
}
int main() {
  fac[0] = 1;
  for (int i = 1; i < mod; ++i) fac[i] = fac[i - 1] * i % mod;
  inv[1] = 1;
  for (int i = 2; i < mod; ++i) inv[i] = mod - (mod / i) * inv[mod % i] % mod;
  mu[1] = 1;
  for (int i = 1; i < M; ++i)
    for (int j = i * 2; j < M; j += i)
      mu[j] -= mu[i];
  for (int c = 2; c < C; ++c) {
    for (int i = 1; i < M; ++i)
      f[c][i] = 0;
    for (int d = c - 1; d < M; ++d) {
      int cc = CC(d - 1, c - 2);
      for (int i = 1; i * d < M; ++i)
        f[c][i * d] = (f[c][i * d] + cc * mu[i]) % mod;
    }
    for (int j = 0; j < N; ++j)
      S[c][0][j] = 0;
    for (int i = 1; i < M; ++i)
      for (int j = 0, ij = 1; j < N; ++j, ij = ij * (i % mod) % mod)
        S[c][i][j] = (S[c][i - 1][j] + f[c][i] * ij) % mod;
  }
  int T;
  scanf("%d", &T);
  while (T--) {
    mm = 10000000;
    scanf("%d%d", &n, &c);
    for (int i = 0; i < n; ++i)
      scanf("%d", &m[i]), mm = std::min(mm, m[i]);
    if (mm < c) {
      printf("0\n");
      continue;
    } if (n == 1) {
      printf("%d\n", CC(m[0], c));
      continue;
    }
    int t = 0;
    p[t++] = 1;
    for (int i = 0; i < n; ++i) {
      int last = 1;
      while (last < m[i] - 1) {
        last = (m[i] - 1) / ((m[i] - 1) / last) + 1;
        p[t++] = last;
      }
    }
    std::sort(p, p + t);
    t = (int)(std::unique(p, p + t) - p);
    p[t] = mm;
    int d = 1;
    while (d < n) d <<= 1;
    for (int i = 0; i < n; ++i) mdp[i] = m[i];
    for (int i = 1; i <= d * 2; ++i) a[i][po[i] = 0] = 1;
    int ans = 0;
    for (int i = 0; p[i] < mm && i < t; ++i) {
      int l = p[i], r = p[i + 1] - 1; //[l, r]
      for (int j = 0; j < n; ++j)
        if ((m[j] - 1) / l != mdp[j]) {
          mdp[j] = (m[j] - 1) / l;
          int k = d + j;
          po[k] = 1;
          a[k][0] = m[j] % mod * (mdp[j] % mod) % mod;
          a[k][1] = -mdp[j] % mod * ((mdp[j] + 1) % mod) % mod * inv[2] % mod;
          while (k /= 2)
            upd(k);
        }
      ans = (ans + sum(l, r)) % mod;
    }
    printf("%d\n", (ans + mod) % mod);
  }
  return 0;
}

UOJ#54 [WC2010]重建計劃