【GDOI2018模拟8.12】求和

Description:

1<=n<=10^10,k = 1、40

题解：

一看就是反演题，立马进行简单变形：
原式=$\sum_{x = 1}^k \sum_{d = 1}^n f_x(d) *\sum_{i = 1}^n \sum_{j = 1}^n [gcd(i, j) = d]$
$=\sum_{x = 1}^k \sum_{d = 1}^n f_x(d) *\sum_{i = 1}^{\lfloor {n\over i} \rfloor} \sum_{j = 1}^{\lfloor {n\over i} \rfloor} [gcd(i, j) = 1]$
$=\sum_{x = 1}^k \sum_{d = 1}^n f_x(d) *((\sum_{i = 1}^{\lfloor {n\over i} \rfloor} φ(i) *2) - 1)$

当k=1时，$f_1(x) = μ(x)$,直接杜教筛μ和φ的前缀和即可。
当k=40时，就不好做了，首先需要注意到2^40>10^10。

设$λ(n) = f_{max}(n)$
若$n = \prod p_{i}^{q_i}$，则$λ(n) = (-1)^{\sum{q_i}}$。

设$F_d(n) = \sum_{i = 1}^n f_d(i)$。
实际上可以用λ来容斥出F_d，思路是这样的：
用总和，减去指数至少为d+1的，就可以枚举一个i，那么$i_{d+1}*j(j<=n/i^{d+1})$一定是指数大于i+1的，这样求会有重复，直接乘上一个μ即可。
所以：
$F_d(n) = \sum_{i = 1}^n μ(i) * \sum_{j = 1}^{n \over{i^{d +1}}}λ(i^{d +1} *j)$
λ显然是一个完全积性函数，所以：
$F_d(n) = \sum_{i = 1}^n μ(i) * λ(i^{d +1}) * \sum_{j = 1}^{n \over{i^{d +1}}}λ( j)$
缩小范围：
$F_d(n) = \sum_{i = 1}^{n^{1\over{d+1}}} μ(i) * λ(i^{d +1}) * \sum_{j = 1}^{n \over{i^{d +1}}}λ( j)$

注意式子后面是λ的前缀和，杜教筛求出这个，对每个$F_d(n)$就能快速求了，复杂度证明不会，但是感性认识一下，这个是指数级减小的，所以应该是很小很小的。

最后一个问题是如何杜教筛λ。

$\sum_{d|n}λ(d)$
$=\sum_{d|n}λ(d)*I({n \over d})$
这是个裸的狄利克雷卷积式，假设$n=p^q$,随便算一下，易得：
$\sum_{d|n}λ(d)= [n是完全平方数]$

$\sum_{i = 1}^n \sum_{j |i} λ(j)$
$=\sqrt n$
$=\sum_{i = 1}^n \sum_{j = 1}^{\lfloor {n\over j }\rfloor}λ(j)$
设$S(n)= \sum_{i = 1}^n λ(i)$
则$S(n) =\sqrt n - \sum_{i = 2}^n S(\lfloor {n\over i}\rfloor)$

#include<cmath>
#include<cstdio>
#define ll long long
#define ld long double
#define fo(i, x, y) for(ll i = x; i <= y; i ++)
#define max(a, b) ((a) > (b) ? (a) : (b))
using namespace std;

const ll N = 5000000, M = 16516513;
const ll mo = 1 << 30, ni_2 = 5e8 + 4;
ll n, m;
bool bz[N + 5];
ll p[N], s[N + 5], mu[N + 5], v[N + 5], c[N + 5], phi[N + 5];
ll la[N + 5], cc[N + 5], sm[N + 5], ss[N + 5];
ll ans;


ll h[M][2];
ll hash(ll x) {
    ll y = x % M;
    while(h[y][0] != 0 && h[y][0] != x) y = (y + 1) % M;
    return y;
}
ll dg(ll n) {
    if(n <= N) return s[n];
    ll y = hash(n);
    if(h[y][0] == n) return h[y][1];
    h[y][0] = n;
    ll ans = sqrt(n * 1.0);
    fo(i, 2, n) {
        ll j = n / (n / i);
        ans -= (j - i + 1) * dg(n / i) % mo;
        i = j;
    }
    ans = (ans % mo + mo) % mo;
    h[y][1] = ans;
    return ans;
}

ll h2[M][2];
ll hash2(ll x) {
    ll y = x % M;
    while(h2[y][0] != 0 && h2[y][0] != x) y = (y + 1) % M;
    return y;
}
ll dg2(ll n) {
    if(n <= N) return phi[n];
    ll y = hash2(n);
    if(h2[y][0] == n) return h2[y][1];
    h2[y][0] = n;
    ll n_ = n % mo;
    ll ans = n_ * (n_ + 1) / 2 % mo;
    fo(i, 2, n) {
        ll j = n / (n / i);
        ans -= (j - i + 1) * dg2(n / i); ans %= mo;
        i = j;
    }
    ans = (ans % mo + mo) % mo;
    h2[y][1] = ans;
    return ans;
}

void Build() {
    mu[1] = 1; c[1] = 1; phi[1] = 1;
    fo(i, 2, N) {
        if(!bz[i]) p[++ p[0]] = i, v[i] = 1, c[i] = 1, s[i] = 1, mu[i] = -1, phi[i] = i - 1;
        fo(j, 1, p[0]) {
            int k = i * p[j];
            if(k > N) break;
            bz[k] = 1;
            s[k] = s[i] ^ 1;
            if(i % p[j] == 0) {
                v[k] = v[i] + 1;
                c[k] = max(c[i], v[k]);
                mu[k] = 0;
                phi[k] = phi[i] * p[j];
                break;
            }
            v[k] = 1;
            c[k] = c[i];
            mu[k] = -mu[i];
            phi[k] = phi[i] * phi[p[j]];
        }
    }
}

ll SS(ll p) {
    if(p <= N) return ss[p];
    ll gg = sqrt(p * 1.0);
    fo(k, 1, gg) cc[k] = k;
    ll sum = 0;
    fo(d, 1, m) {
        ll mm = pow(p, (ld) 1 / (d + 1));
        fo(i, 1, mm) {
            cc[i] *= i;
            ll lb = (((d + 1) % 2) ? la[i] : 1) * mu[i];
            sum += lb * dg(p / cc[i]); sum %= mo;
        }
    }
    return sum;
}

int main() {
    Build();
    scanf("%lld %lld", &n, &m);
    fo(i, 1, N) phi[i] = (phi[i - 1] + phi[i]) % mo;
    fo(i, 1, N) {
        ss[i] = ss[i - 1];
        if(m == 1)
            ss[i] = (ss[i] + (c[i] <= 1) * (s[i] ? -1 : 1)) % mo; else
            ss[i] = (ss[i] +  (s[i] ? -1 : 1) * (m - c[i] + 1)) % mo;
    }
    fo(i, 1, N) la[i] = s[i] ? -1 : 1, s[i] = s[i - 1] + la[i];
    fo(i, 1, N) sm[i] = sm[i - 1] + mu[i];
    fo(i, 1, n) {
        ll c = n / i;
        ll j = n / c;
        ans += (SS(j) - SS(i - 1)) * (2 * dg2(n / i) - 1); ans %= mo;
        i = j;
    }
    ans = (ans % mo + mo) % mo;
    printf("%lld", ans);
}