HDU 6706 huntian oy [杜教筛]

Problem

求解$\sum ^{n}{i=1}\sum ^{i}{j=1}\gcd ( i^{a}-j^{a} , i^b-j^b)[\gcd(i,j)=1][\gcd(a,b)=1]$

$1\le n\le {10}^9$

Solution

因为$gcd(i,j)=1$

原式等价于$\sum ^{n}{i=1}\sum ^{i}{j=1}i^{\gcd (a,b)}-j^{\gcd (a,b)}[\gcd(i,j)=1][\gcd(a,b)=1]$

又$gcd(a,b)=1$

式子化简为$\sum ^{n}{i=1}\sum ^{i}{j=1}i-j[\gcd(i,j)=1]$

对称变换得$(\sum ^{n}{i=1}\sum ^{i}{j=1}j[\gcd(i,j)=1])-1$

$\sum _{i=1}^{n}i[gcd(i,n)=1]=\frac{n\ast \phi (n)+[n=1]}{2}$

所以原式等价于$\frac{(\sum _{i=1}^{n}i\ast \phi (i))-1}{2}$

记$S(n)=\sum _{i=1}^{n}i\ast \phi (i)，f(n)=n\ast \phi (n)$

考虑到狄利克雷卷积$(\phi \ast I)(n)=id$

我们对$f$卷上单位函数$id$来消除系数

$(f\ast id)(n)=\sum _{d|n}f(d)\ast \frac{n}{d}=\sum _{d|n}d\ast \phi (d)\ast \frac{n}{d}=n\ast \sum _{d|n}\phi (d)=n^2$

$\sum _{i=1}^n{i}^2=\sum _{i=1}^n\sum _{d|i}d\ast \phi (d)\ast {\displaystyle \frac{i}{d}}$

枚举约数转枚举倍数有 $\sum _{d=1}^n\sum _{d|i}\frac{i}{d}\ast \phi (\frac{i}{d})\ast d$

提取$d$之后发现我们可以得到原函数

$\sum _{d=1}^{n}d\sum _{d|i}{\displaystyle \frac{i}{d}}\ast \phi (\frac{i}{d})=\sum _{d=1}^{n}dS(\lfloor {\displaystyle \frac{n}{d}}\rfloor )$

则有$S(n)=\sum _{i=1}^n{i}^2-\sum _{i=2}^{n}iS(\lfloor \frac{n}{i}\rfloor )$

预处理$S$的$n^{\frac{2}{3}}$项，递归数值分块计算即可

#include <bits/stdc++.h>
typedef long long ll;
const int MOD = 1e9 + 7, N = 2e6 + 5;
bool notp[N];
int prime[N], pnum;
ll phi[N];
void sieve() {
    memset(notp, 0, sizeof(notp));
    notp[0] = notp[1] = phi[1] = 1;
    pnum = 0;
    for (int i = 2; i < N; i++) {
        if (!notp[i]) {
            prime[++pnum] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= pnum && prime[j] * i < N; j++) {
            notp[prime[j] * i] = 1;
            if (i % prime[j] == 0) {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
}
ll inv(ll a, ll m) {
    if (a == 1) return 1;
    return inv(m % a, m) * (m - m / a) % m;
}
ll inv6, inv2;
void init() {
    sieve();
    for (int i = 1; i < N; i++) phi[i] = (1ll * phi[i] * i + phi[i - 1]) % MOD;
}
ll mp[2][N], n;
ll f(ll x) { return (1ll * x * (x + 1) / 2) % MOD; }
std::unordered_map<ll, ll> hs;
ll S(ll n) {
    if (n < N) return phi[n];
    if (hs[n]) return hs[n];
    ll res = (n % MOD) * ((n + 1) % MOD) % MOD * ((2 * n + 1) % MOD) % MOD *
             inv6 % MOD;
    for (ll i = 2, last; i <= n; i = last + 1) {
        last = n / (n / i);
        res += MOD - (f(last) - f(i - 1) + MOD) % MOD * S(n / i) % MOD;
        res %= MOD;
    }
    return hs[n] = res;
}
int T;
int main() {
    scanf("%d", &T);
    memset(mp, -1, sizeof(mp));
    init();
    inv6 = inv(6, MOD);
    inv2 = inv(2, MOD);
    while (T--) {
        int a, b;
        scanf("%lld%d%d", &n, &a, &b);
        printf("%lld\n", (S(n) - 1 + MOD) % MOD * inv2 % MOD);
    }
    return 0;
}