Problem: tower_rush (2023 round2)

Alice and Bob are playing *Tower Rush,* a two-phase game involving \(N\) types of blocks, numbered from \(1\) to \(N\). Blocks of type \(i\) have height \(H_i\), and no two types of blocks have the same height.

Phase 1 of the game consists of \(K \ge 2\) alternating turns, with Alice going first. On each player's turn, they will choose a block type that has not yet been chosen. Note that if \(K\) is odd, this phase would end with Alice having chosen one more type of block than Bob. From this point on, players have access to infinitely many blocks of each of the types they chose.

Phase 2 consists of an indefinite number of alternating turns, with Bob going first. Each player starts with a tower of height \(0\). On each player's turn, they may pick a block of any type \(i\) available to them, and extend their tower height by \(H_i\). They may also choose to skip their turn, leaving their tower unchanged.

In a break from tradition, Alice and Bob want to work *together* to see if it's possible for Alice to build a tower that's exactly \(D\) units taller than Bob's. In how many different ways can phase 1 be played out such that it will be possible for Alice to get her tower to be exactly \(D\) units taller than Bob's in phase 2? Two ways are considered different if there exists an \(i\) such that different block types are chosen on turn \(i\) between the two ways.

As this value may be large, output it modulo \(1{,}000{,}000{,}007\).


# Constraints

\(1 \leq T \leq 24\)
\(2 \leq N \leq 1{,}000{,}000\)
\(2 \leq K \leq \min(N, 20)\)
\(1 \leq D \leq 1{,}000{,}000\)
\(1 \leq H_i \leq 1{,}000{,}000\)

The sum of \(N\) across all test cases is at most \(7{,}500{,}000\).


# Input Format

Input begins with an integer \(T\), the number of test cases. For each case, the first line contains the integers \(N\), \(K\), and \(D\). Then, there is a line containing the \(N\) integers \(H_1\), ..., \(H_N\).


# Output Format

For the \(i\)th test case, print "`Case #i:` " followed by a single integer, the number of unique ways to accomplish Bob's goal, modulo \(1{,}000{,}000{,}007\).


# Sample Explanation

In the first sample case, there are \(54\) valid choices of blocks in phase 1. One such example is \(\{H_2 = 4,\, H_3 = 5\}\) for Alice and \(\{H_4 = 7\}\) for Bob. Four turns in phase 2 can go as follows:
* Bob adds \(H_4 = 7\). His tower now has height \(7\).
* Alice adds \(H_3 = 5\). Her tower now has height \(5\).
* Bob chooses not to add a block. His tower still has height \(7\).
* Alice adds \(H_3 = 5\). Her tower now has height \(10\).

Alice’s tower is now exactly \(3\) units taller than Bob’s, so this is a valid choice.

An example of an invalid choice in the first test case would be \(\{H_1 = 2,\, H_2=4\}\) for Alice and \(\{H_5 = 8\}\) for Bob, since no matter how they build their towers, it is impossible for the difference in the height of their towers to ever be \(3\).


Sample input:
2
5 3 3
2 4 5 7 8
3 2 4
3 6 9


Sample output:
Case #1: 54
Case #2: 0


Solution: First note that, for a set of blocks, it does not matter how they are distributed between the players. That is, if two blocks get swapped this will not change the outcome (possible or not to have difference of d). This comes from the fact that for a given distribution of blocks, we can always form two towers of equal height (with height being the LCM of all elements in the set).

With that, we can effectively either remove or add a block with height \(H_i\) to a tower, which corresponds to either adding or subtracting \(H_i\) to the current difference between the towers. It is a well known result that given integers \((a,b)\), with the equation \(a*x + b*y = k\), where \(x\) and \(y\) can assume any integer values, \(k\) can assume any value that is a multiple of \(\text{gcd}(a,b)\). This results can is also valid for multiple variables (see Bézout's identity). For a given set of blocks, we can make it so that the difference assumes any value that is a multiple of the gcd of the chosen blocks.

Now the problem simplifies to, for each divisor \(d\) of \(D\), consider how many subsets of size \(K\) are there such that the gcd is \(d\). This can be solved in \(\mathcal{O}(MAXV/i)\). It is easy to convert this problem to, how many subsets of size \(k\) are there such that GCD of values is \(1\). This final problem can be solved with Inclusion-Exclusion principle: for each integer \(x\), it is easy to count how many subsets of size \(k\) are there such that the gcd is a multiple of \(x\) (let's call this \(f(x)\)), this is simply \(\binom{\text{freq}[x]}{k}\) where \(\text{freq}[x]\) is how many integers are divisible by \(x\). The solution is

\[\sum_{x=1}^{MAXV} [f(x)*\mu(x)]\]

where \(\mu(x)\) is the Mobious function (that is, for square free integers, \(\mu(x)\) is \(1\) if \(x\) has an even number of prime factors, or \(-1\) if it has an odd number of prime factors, which is related to inclusion exclusion). This works in \(\mathcal{O}(N \log N)\) for any value of \(K\), not necessarily \(K \le 20\).


Code: #include <cassert>
#include <iostream>
#include <vector>
using namespace std;

const int MAX = 1000010;
const int MOD = 1000000007;
using LL = long long;

class mint {
  LL val;

  mint& normalize(LL x) {
    val = x % MOD;
    if (val < 0) {
      val += MOD;
    }
    return *this;
  }

  static LL inverse(LL a) {
    LL u = 0, v = 1, m = MOD;
    while (a != 0) {
      LL t = m / a;
      m -= t * a;
      swap(a, m);
      u -= t * v;
      swap(u, v);
    }
    assert(m == 1);
    return u;
  }

 public:
  mint(LL x = 0) { normalize(x); }

  LL operator()() const { return val; }

  template <typename U>
  explicit operator U() const { return static_cast<U>(val); }

  mint& operator=(LL x) { return normalize(x); }
  mint& operator=(const mint& x) { val = x.val; return *this; }

  mint& operator+=(const mint& x) { return normalize(val + x.val); }
  mint& operator-=(const mint& x) { return normalize(val - x.val); }
  mint& operator*=(const mint& x) { return normalize(val * x.val); }
  mint& operator/=(const mint& x) { return *this *= mint(inverse(x.val)); }

  mint& operator+=(LL x) { return *this += mint(x); }
  mint& operator-=(LL x) { return *this -= mint(x); }
  mint& operator*=(LL x) { return *this *= mint(x); }
  mint& operator/=(LL x) { return *this /= mint(x); }

  mint& operator++() { return *this += 1; }
  mint& operator--() { return *this -= 1; }
  mint operator++(int) { mint z(*this); ++*this; return z; }
  mint operator--(int) { mint z(*this); --*this; return z; }

  friend mint operator+(mint x, const mint& y) { return x += y; }
  friend mint operator*(mint x, const mint& y) { return x *= y; }
  friend mint operator-(mint x, const mint& y) { return x -= y; }
  friend mint operator/(mint x, const mint& y) { return x /= y; }

  friend mint operator+(mint x, LL y) { return x += y; }
  friend mint operator*(mint x, LL y) { return x *= y; }
  friend mint operator-(mint x, LL y) { return x -= y; }
  friend mint operator/(mint x, LL y) { return x /= y; }

  friend mint operator+(LL x, mint y) { return y += x; }
  friend mint operator*(LL x, mint y) { return y *= x; }
  friend mint operator-(LL x, const mint& y) { mint z(x); return z -= y; }
  friend mint operator/(LL x, const mint& y) { mint z(x); return z /= y; }

  bool operator <(const mint& x) const { return val < x.val; }
  bool operator==(const mint& x) const { return val == x.val; }
  bool operator >(const mint& x) const { return val > x.val; }
  bool operator!=(const mint& x) const { return val != x.val; }
  bool operator<=(const mint& x) const { return val <= x.val; }
  bool operator>=(const mint& x) const { return val >= x.val; }

  bool operator <(LL x) const { return val < x; }
  bool operator==(LL x) const { return val == x; }
  bool operator >(LL x) const { return val > x; }
  bool operator!=(LL x) const { return val != x; }
  bool operator<=(LL x) const { return val <= x; }
  bool operator>=(LL x) const { return val >= x; }
};

class factorial : vector<mint> {
  void lazy_eval(int n) {
    for (LL p = size(); n >= p; ++p) {
      push_back(back() * p);
    }
  }

 public:
  factorial() { push_back(1); }

  mint choose(LL n, LL k) {
    if (n < 0 || k < 0 || n < k) {
      return 0;
    }
    if (k == 0 || k == n) {
      return 1;
    }
    if (n >= MOD && n <= k + MOD - 1) {
      return 0;
    }
    lazy_eval(n);
    return at(n) / (at(k) * at(n - k));
  }

  mint operator[](LL n) {
    lazy_eval(n);
    return at(n);
  }
};

factorial F;
vector<vector<LL>> divisors(MAX);
vector<LL> mu(MAX);

void build_mu() {
  vector<bool> B(MAX);
  vector<LL> primes;
  mu[1] = 1LL;
  for (LL i = 2; i < MAX; i++) {
    if (!B[i]) {
      primes.push_back(i);
      mu[i] = -1;
    }
    for (auto& p : primes) {
      LL k = (LL)i * p;
      if (k > MAX) {
        break;
      }
      B[k] = 1;
      if (i % p != 0) {
        mu[k] = -mu[i];
      } else {
        mu[k] = 0;
        break;
      }
    }
  }
}

LL solve() {
  int N, K, D;
  cin >> N >> K >> D;
  vector<int> H(N), freq(MAX);
  for (int i = 0; i < N; i++) {
    cin >> H[i];
    for (LL d : divisors[H[i]]) {
      freq[d]++;
    }
  }
  mint ans = 0;
  for (LL d : divisors[D]) {
    for (LL j = 0; j < MAX; j++) {
      LL i = (LL)d * j;
      if (i > MAX) {
        break;
      }
      ans += F.choose(freq[i], K) * mu[j];
    }
  }
  return (LL)(ans * F[K]);
}

int main() {
  ios_base::sync_with_stdio(false);
  cin.tie(nullptr);

  for (int i = 1; i < MAX; i++) {
    for (int j = i; j < MAX; j += i) {
      divisors[j].push_back(i);
    }
  }
  build_mu();

  int T;
  cin >> T;
  for (int t = 1; t <= T; t++) {
    cout << "Case #" << t << ": " << solve() << endl;
  }
  return 0;
}