标签归档:数论

[2018 Multi-University Training Contest 10] C. Calculate

Link: http://acm.hdu.edu.cn/showproblem.php?pid=6428

Statement:

Given \(1 \leq A, B, C \leq 10^7\), compute

\[ \sum_{i=1}^A \sum_{j=1}^B \sum_{k=1}^C \phi(\gcd(i, j^2, k^3)) \]

no more than 10 test cases.

Solution:

Consider the number of triples (i, j, k) such that \( \gcd(i, j^2, k^3) = d\), denoted as \(f(d)\). By Mobius inversion,

\[ f(d) = \sum_{i} \mu(i) g(di) \]

where \(g(d)\) is the number of triples satisfying \( d | \gcd(i, j^2, k^3) \), and the answer should be

\[ \sum_{d=1}^A f(d) \phi(d) \]

Now we try to compute \(g(d)\):

\[ g(d) = \sum_{i = 1}^A [d | i] \sum_{j = 1}^B [d | j^2] \sum_{k = 1}^C [d | k^3] \]

Let \(d = \prod {p_i}^{a_i} \) denote the unique factorization of d. Note that \(d | k^n \) whenever \( \prod {p_i}^{\lceil a_i / n \rceil} | k \).

Let \( h_n(d) = \prod {p_i}^{\lceil a_i / n \rceil} \), which denotes the minimum positive integer whose nth power is divisible by d. Now \(g(d)\) could be rewritten as

\[ g(d) = \lfloor \frac{A}{d} \rfloor \lfloor \frac{B}{h_2(d)} \rfloor \lfloor \frac{C}{h_3(d)} \rfloor \]

Plug in answer:

\[ \sum_{d=1}^A f(d) \phi(d) = \sum_{d=1}^A \phi(d) \sum_{i} \mu(i) g(di) = \sum_{q = 1}^A g(q) \sum_{ij = q} \mu(i) \phi(j) \]

where the Dirichlet convolution \( (\mu * \phi)(q) = \sum_{ij = q} \mu(i) \phi(j) \) is multiplicative and thus could be effectively precomputed (along with \(h_2\), \(h_3\)) by the sieve of Euler in linear time. Hence each test case could be done in \(O(A)\) time.

Source Code:

#include <bits/stdc++.h>
using namespace std;

#ifdef __LOCAL_DEBUG__
# define _debug(fmt, ...) fprintf(stderr, "\033[94m%s: " fmt "\n\033[0m", \
    __func__, ##__VA_ARGS__)
#else 
# define _debug(...) ((void) 0)
#endif
#define rep(i, n) for (int i=0; i<(n); i++)
#define Rep(i, n) for (int i=1; i<=(n); i++)
#define range(x) (x).begin(), (x).end()
typedef long long LL;
typedef unsigned long long ULL;

namespace sieve {
  constexpr int MAXN = 10000007;
  bool p[MAXN];
  int prime[MAXN], sz;
  int pval[MAXN], pcnt[MAXN];
  unsigned h2[MAXN], h3[MAXN], phi[MAXN], muphi[MAXN];

  void exec(int N = 10000007) {
    p[0] = p[1] = 1;
    
    pval[1] = 1;
    pcnt[1] = 0;
    h2[1] = h3[1] = phi[1] = muphi[1] = 1;

    for (int i = 2; i < N; i++) {
      if (!p[i]) {
        prime[sz++] = i;
        for (LL j = i; j < N; j *= i) {
          int b = j / i;
          pval[j] = i * pval[b];
          pcnt[j] = pcnt[b] + 1;
          h2[j] = h2[b] * (pcnt[j] % 2 == 1 ? i : 1);
          h3[j] = h3[b] * (pcnt[j] % 3 == 1 ? i : 1);
          phi[j] = (i == j ? i-1 : phi[b] * i);
          muphi[j] = phi[j] - phi[b];
        }
      }
      for (int j = 0; i * prime[j] < N; j++) {
        int x = i * prime[j]; p[x] = 1;
        if (i % prime[j] == 0) {
          pval[x] = pval[i] * prime[j];
          pcnt[x] = pcnt[i] + 1;
        } else {
          pval[x] = prime[j];
          pcnt[x] = 1;
        }
        if (x != pval[x]) {
          int b = x / pval[x];
          h2[x] = h2[b] * h2[pval[x]];
          h3[x] = h3[b] * h3[pval[x]];
          muphi[x] = muphi[b] * muphi[pval[x]];
        }
        if (i % prime[j] == 0) break;
      }
    }
  }
}

int main() {
  using namespace sieve;
  exec();
  int T; scanf("%d", &T);
  while (T--) {
    unsigned a, b, c; scanf("%d%d%d", &a, &b, &c);
    unsigned res = 0;
    for (unsigned i = 1; i <= a; i++) 
      res += muphi[i] * (a / i) * (b / h2[i]) * (c / h3[i]); 
    printf("%d\n", res & 0x3fffffff);
  }
  return 0;
}