(Continues from Number Theory #2)

0.1 Basic Identities

$$\bullet\;[n=1]=\sum_{d\mid n}\mu(d)$$

$$\bullet\;n=\sum_{d\mid n}\varphi(d)\;\implies\;\varphi(n)=\sum_{d\mid n}\mu(d)\cdot\frac{n}{d}$$

$$\bullet\;\sigma(n)=\sum_{d\mid n}d$$

Indicator / sieve identities

• Coprimality indicator: $\displaystyle[\gcd(x,n)=1]=\sum_{d\mid\gcd(x,n)}\mu(d)$
• Congruence indicator: $[x\equiv y\pmod{d}]=[d\mid(x-y)]$

Classical divisor sums

• $\displaystyle\sum_{d\mid n}\varphi(d)=n$
• $\displaystyle\sum_{d\mid n}\mu(d)=[n=1]$
• $\displaystyle\sum_{d\mid n}\tau(d)=\sum_{d\mid n}\tau(n/d)\quad(\tau=1*1,\ \tau\text{ is the divisor-count function})$
• $\sigma=\mathrm{id}*\mathbf{1},\quad\tau=\mathbf{1}*\mathbf{1},\quad\varphi=\mu*\mathrm{id}$

GCD variants

• $\displaystyle\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)$
• $\displaystyle\sum_{d\mid n}\mu(d)\Bigl\lfloor\frac{m}{d}\Bigr\rfloor$ = “count of $x\in[1,m]$ coprime to $n$” (inclusion-exclusion over prime factors of $n$)
• $\displaystyle d(xy)=\sum_{a\mid x}\sum_{b\mid y}[\gcd(a,b)=1]$

GCD-layered substitution (very common)

(Group $k$ by $\gcd(k,n)=d$.)

1 Trick 1 — Loop Shortening via Divisibility

When iterating $i=1,\ldots,n$, if the expression is non-zero only when $k\mid i$, shorten the loop:

$$\sum_{i=1}^{n}[k\mid i]=\sum_{i=1}^{\lfloor n/k\rfloor}1=\Bigl\lfloor\frac{n}{k}\Bigr\rfloor$$

Example: $\displaystyle\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=k]$

Substitute $i=i'k$, $j=j'k$; the condition $\gcd(i,j)=k$ becomes $\gcd(i',j')=1$:

Key point: non-zero entries appear only every $k$ steps in both loops, so the range shrinks by factor $k$.

2 Trick 2 — Convert GCD-Divisors to Individual Divisors

Core idea: $\gcd(i,j)$ varies with $i,j$ so its divisors are hard to enumerate directly. Use $d\mid\gcd(i,j)\iff d\mid i\text{ and }d\mid j$ to decouple the two indices.

Example: $\displaystyle\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1]$

The last sum is evaluated in $O\!\left(\sqrt{\min(n,m)}\right)$ via divisibility blocking.

Example 1: Luogu P3455 [POI2007] ZAP-Queries

• Trick 1: reduce to $\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor}[\gcd(i,j)=1]$.
• Trick 2: $=\displaystyle\sum_{d=1}^{\min(\lfloor n/k\rfloor,\lfloor m/k\rfloor)}\mu(d)\Bigl\lfloor\frac{\lfloor n/k\rfloor}{d}\Bigr\rfloor\Bigl\lfloor\frac{\lfloor m/k\rfloor}{d}\Bigr\rfloor$.

Supplementary identity: $\displaystyle\Bigl\lfloor\frac{\lfloor n/k\rfloor}{d}\Bigr\rfloor=\Bigl\lfloor\frac{n}{kd}\Bigr\rfloor$. The two forms are equivalent; the right-hand form is cleaner in derivations, the left-hand form is more natural in code.

int main(){
    long long t;
    read(t);
    get_mu(50000);
    while(t--){
        static long long a, b, d;
        read(a); read(b); read(d);
        static long long max_rep, ans;
        max_rep = min(a, b); ans = 0;
        for(long long l=1, r; l<=max_rep; l=r+1){
            r = min(a/(a/l), b/(b/l));
            ans += (a/(l*d)) * (b/(l*d)) * (sum[r] - sum[l-1]);
        }
        printf("%lld\n", ans);
    }
    return 0;
}

Example 2: Luogu P2522 [HAOI2011] Problem b

Reduce to four $(1,\cdot)$-based queries via 2D prefix-sum inclusion-exclusion:

3 Trick 3 — Reduce a General $f(\gcd)$ Sum to $[\gcd=1]$ Form

Replace the inner expression $f(\gcd(i,j))$ by $f(d)\cdot[\gcd(i',j')=1]$ by factoring out $d=\gcd(i,j)$.

The most common Möbius pattern: $\displaystyle\sum_{i=1}^n\sum_{j=1}^mf(\gcd(i,j))$.

Set $d=\gcd(i,j)$, $i=i'd$, $j=j'd$; then $\gcd(i',j')=1$ (otherwise $d$ is not maximal).

Add an outer loop over $d$ and shrink the inner ranges (Trick 1):

4 Trick 4 — Substitute $T=dk$ to Identify a Dirichlet Convolution

When outer loop is $d$ and inner variable is $k$, and only their product $T=dk$ appears in floor expressions, swap to outer loop over $T$ and enumerate divisors $d\mid T$.

Starting from Trick 3:

Step 1. Expand $[\gcd(i,j)=1]$ with Möbius, swap summation order:

Step 2. Replace divisibility counts:

Step 3. Set $T=dk$. For fixed $T$, all pairs $(d,k)$ with $dk=T$ correspond to divisors $d\mid T$ with $k=T/d$. Since $dk\leq\min(n,m)$, we have $T\leq\min(n,m)$:

Step 4. Recognise the Dirichlet convolution $g=f*\mu$:

General Formula 1. Combining Tricks 3 and 4 (assuming $n\leq m$):

$$\sum_{i=1}^n\sum_{j=1}^m f(\gcd(i,j))=\sum_{T=1}^{\min(n,m)}g(T)\Bigl\lfloor\frac{n}{T}\Bigr\rfloor\Bigl\lfloor\frac{m}{T}\Bigr\rfloor,\quad g=f*\mu$$

Limitation: the inner expression must depend only on $\gcd(i,j)$ (a univariate function). This formula does not generalise to expressions involving $i$ or $j$ independently.

4.1 Example: ZAP-Queries via General Formula 1

With $f(x)=[x=k]$:

Apply Trick 1 (set $T=ik$):

Evaluate with divisibility blocking.

5 Trick 5 — Divisor Symmetry: $d\leftrightarrow n/d$

When $d$ runs over all divisors of $n$, $d$ and $n/d$ are in bijection.

Example:

(Replace $d\to n/d$ throughout the divisor sum.)

6 Trick 6 — Diagonal Symmetry: Upper Triangle vs Full Matrix

When $f(i,j)=f(j,i)$, convert between the full double sum and the upper-triangle sum.

6.1 Upper-triangle and lower-triangle

Both sides sum the upper-triangle region (including the diagonal); only the summation order differs.

6.2 Upper-triangle sum and full-matrix sum

“Upper triangle including diagonal” = (full matrix + diagonal) / 2.

Equivalently:

7 Trick 7 — Value-Domain Möbius with $\mathrm{cnt}[\,]$

Overview: standard Möbius problems operate on $1,2,\ldots,n$. When the input is an array $a[1\cdots n]$ and summation involves number-theoretic functions of array values, introduce:

• $w[y]$ = number of times value $y$ appears in $a[\cdot]$
• $\mathrm{cnt}[d]$ = number of elements in $a[\cdot]$ divisible by $d$

Example. Given $a[1\cdots n]$ with $1\leq a_i\leq n$, compute:

Split into two parts:

Computing $P_2$ (non-coprime pairs):

Expand $[\gcd=1]=\sum_{d\mid\gcd}\mu(d)$, decouple with Trick 2:

Computing $P_1$ (sum of GCDs):

Use Euler inversion $\mathrm{id}=\varphi*\mathbf{1}$, i.e.\ $\gcd(a,b)=\sum_{d\mid\gcd(a,b)}\varphi(d)$, then Trick 2:

Final answer:

Computing $\mathrm{cnt}[d]$: build $w[y]$ (frequency array), then for each $d$ accumulate $w[d]+w[2d]+\cdots$ — harmonic-series cost $O(V\log V)$.

Complexity: $O(V\log V)$ overall, $V=\max\{a_i\}$.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int N = 100000 + 10;   // set to max(a[i]) + 5 as needed
const int P = 998244353;

int n, v;
int a[N], w[N];

bool st[N];
int primes[N], tot;
int mu[N]; ll phi[N];
ll cnt[N];
vector<int> g[N];  // g[x][0] = smallest prime factor of x

// sieve + phi + mu
void get_primes(int nmax) {
    phi[1] = 1;
    for (int i = 1; i <= nmax; i++) mu[i] = 1;
    for (int i = 2; i <= nmax; i++) {
        if (!st[i]) {
            primes[tot++] = i;
            mu[i] = -1;
            phi[i] = i - 1;
            for (int j = i + i; j <= nmax; j += i) {
                st[j] = true;
                g[j].push_back(i);
                if ((j / i) % i == 0) mu[j] = 0;
                else if (mu[j] != 0) mu[j] = -mu[j];
            }
        } else {
            int p = g[i][0];        // smallest prime factor of i
            int now = i / p;
            if (now % p == 0) phi[i] = 1LL * p * phi[now];
            else              phi[i] = 1LL * (p - 1) * phi[now];
        }
    }
}

// cnt[d] = number of elements in a[] divisible by d
void get_cnt() {
    for (int i = 1; i <= n; i++) w[a[i]]++;
    cnt[1] = n;
    for (int d = 2; d <= v; d++) {
        for (int j = d; j <= v; j += d) {
            cnt[d] += w[j];   // no modular reduction: cnt[d] <= n
        }
    }
}

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

    cin >> n;
    v = 0;
    for (int i = 1; i <= n; i++) {
        cin >> a[i];
        v = max(v, a[i]);
    }

    get_primes(v);
    get_cnt();

    ll ans = 1LL * n * n % P;
    for (int d = 1; d <= v; d++) {
        ll term = (((phi[d] - mu[d]) % P) + P) % P;
        ll cd   = cnt[d] % P;
        ans = (ans + term * cd % P * cd) % P;
    }
    cout << ans << '\n';
    return 0;
}