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

### 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$

$\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;
}


# k次幂前缀和的线性做法

$\sum_{i=1}^{x} i^k \bmod p$

### 线性预处理逆元

$i^{-1} = -\left\lfloor \frac{p}{i} \right\rfloor \cdot (p \bmod i)^{-1} \bmod p$

$-\frac{\left\lfloor \frac{p}{i} \right\rfloor}{p-\left\lfloor \frac{p}{i} \right\rfloor i} = \frac{1}{i}$

### 拉格朗日插值法

$f(x) = \sum_{i=0}^{k+1} y_i l_i(x)$

$l_i(x) = \prod_{0 \leq m \leq k+1, m \neq i} \frac{x – x_m}{x_j – x_m}$

$f(n) = \sum_{i = 0}^{k+1} (-1)^{k+i+1} f(i) \frac{n(n-1) \cdots (n-i+1)(n-i-1) \cdots (n-d)}{(d-i)!i!}$

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

#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;

#define pow owahgrauhgso

ULL n, k, m;

ULL powmod(ULL b, ULL e) {
ULL r = 1;
while (e) {
if (e & 1) r = r * b % m;
b = b * b % m;
e >>= 1;
}
return r;
}

const int MAXN = 1000006;
ULL inv[MAXN];

void init_inv() {
inv[1] = 1;
for (int i = 2; i <= k+1; i++) {
inv[i] = (m - m / i) * inv[m % i] % m;
assert(inv[i] * i % m == 1);
}
}

ULL pow[MAXN];
ULL prime[MAXN], cnt;

void sieve() {
pow[1] = 1;
for (int i = 2; i <= k+1; i++) {
if (!pow[i]) {
pow[i] = powmod(i, k);
prime[cnt++] = i;
for (int j = 0; j < cnt && i*prime[j] <= k+1; j++) {
pow[i * prime[j]] = pow[i] * pow[prime[j]] % m;
if (i % prime[j] == 0) break;
}
}
}
}

ULL sum[MAXN];
ULL ans[MAXN];

auto addmod = [](ULL a, ULL b) -> ULL {return (a+b)%m;};

ULL lagrange() {
ULL p;
p = 1;
for (int i=0; i<=k+1; i++) {
if (i) p = p * inv[i] % m;
ans[i] = (k+1-i)&1 ? m-sum[i] : sum[i];
ans[i] = ans[i] * p % m;
p = p * (m + n - i) % m;
}
p = 1;
for (int i=k+1; i>=0; i--) {
if (k+1-i) p = p * inv[k+1-i] % m;
ans[i] = ans[i] * p % m;
p = p * (m + n - i) % m;
}
}

int main() {
cin >> n >> k >> m;
init_inv(); sieve();
cout << lagrange() << endl;
return 0;
}

# 2017-2018年问题求解（三）期末上机考试 Day1 题解

### A. 智能机器

$$n = q_1^{a_1} q_2^{a_2} \cdots q_m^{a_m}$$

$$d(n) = (1+a_1)(1+a_2)\cdots(1+a_m)$$

$$d(n^k) = (1+ka_1)(1+ka_2)\cdots(1+ka_m)$$