CONTENT: 线性筛与积性函数 III
DETAIL: Min_25筛



BB

今天因为BZOJ上多了个测试数据0，调了两个钟 QAQQQQQ

Introduction


Note



【模板】Min_25筛 - luogu P5325

Description
定义积性函数f(x)，且f(p^k)=p^k(p^k−1)（p是一个质数），求
sum(f(x))
对1e9 + 7取模。
Sample Input

10
1000000000

Sample Output

263
710164413

Solution

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = (int)1e6 + 15;
const ll MOD = (ll)1e9 + 7;

const ll inv6 = 166666668;
bool isNotPrime[N];
ll prime[N], tot;
ll sum1[N], sum2[N], g1[N], g2[N], w[N];
ll idx1[N], idx2[N];
int gTot;
ll sqrtN;

inline void init(int n) {
    for(int i = 2; i <= n; i++) {
        if(!isNotPrime[i]) {
            prime[++tot] = i;
            sum1[tot] = (sum1[tot - 1] + i) % MOD;
            sum2[tot] = (sum2[tot - 1] + (ll)i * i) % MOD;
        }
        for(int j = 1; i * prime[j] <= n; j++) {
            isNotPrime[i * prime[j]] = true;
            if(i % prime[j] == 0) {
                break;
            }
        }
    }
}

inline void initG(ll n) {
    // calculate g1[gTot] = 1 + 2 + 3 + ...
    // calculate g2[gTot] = 1^2 + 2^2 + 3^3 + ...
    for(ll i = 1, r; i <= n; i = r + 1) {
        w[++gTot] = n / i;
        r = n / (n / i);
        if(w[gTot] <= sqrtN) {
            idx1[n / i] = gTot;
        } else {
            idx2[n / (n / i)] = gTot;
        }

        ll x = w[gTot] % MOD;
        g1[gTot] = (x * (x + 1) / 2 + MOD - 1) % MOD;
        g2[gTot] = (x * (x + 1) % MOD * (2 * x + 1) % MOD * inv6 + MOD - 1) % MOD;
    }
    // enumerate j, then enumerate n
    // g(i, j - 1) can be calculated before g(n, j)
    for(int i = 1; i <= tot; i++) {
        for(int j = 1; j <= gTot && (ll)prime[i] * prime[i] <= w[j]; j++) {
            ll k = w[j] / prime[i] <= sqrtN ? idx1[w[j] / prime[i]] : idx2[n / (w[j] / prime[i])];
            g1[j] = (g1[j] - (ll)prime[i] * (g1[k] - sum1[i - 1] + MOD) % MOD + MOD) % MOD;
            g2[j] = (g2[j] - (ll)prime[i] * prime[i] % MOD * (g2[k] - sum2[i - 1] + MOD) % MOD + MOD) % MOD;
        }
    }
}

inline ll calcS(ll x, int j, ll n) {
    if(prime[j] >= x) {
        return 0;
    }
    ll k = x <= sqrtN ? idx1[x] : idx2[n / x];
    ll ans = ((ll)g2[k] - g1[k] - (sum2[j] - sum1[j]) + 2LL * MOD) % MOD;
    for(int i = j + 1; i <= tot && (ll)prime[i] * prime[i] <= x; i++) {
        for(ll e = 1, pie = prime[i]; pie <= x; pie *= prime[i], e++) {
            ll xx = pie % MOD;
            ans = (ans + xx * (xx - 1) % MOD * (calcS(x / pie, i, n) + (e != 1))) % MOD;
        }
    }
    return ans;
}

int main() {
    ll n;
    while(~scanf("%lld", &n)) {
        sqrtN = sqrt(n);
        init(sqrtN);
        initG(n);
        printf("%lld\n", (calcS(n, 0, n) + 1) % MOD);
    }
}


【模板】杜教筛（Sum） - luogu P4213

Description
求sum(mu(i))与sum(phi(i))
Sample Input

6
1
2
8
13
30
2333

Sample Output

1 1
2 0
22 -2
58 -3
278 -3
1655470 2

Solution
对于莫比乌斯函数，直接直接构造f(x)=-1，这是显然的
对于欧拉函数，直接构造f(x)=x-1，这也是显然的= =
那么就拆分为求g1(x)=sum(1)，g2(x)=sum(x)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = (int)1e5 + 15;

const int LIM = (1LL << 31) - 1;
bool isNotPrime[N];
int prime[N];
int tot;
int sum1[N], g1[N];
ll sum2[N], g2[N];
int w[N];
int idx1[N], idx2[N];
int gTot;
int sqrtN;

inline void init(int n) {
    tot = 0;
    for(int i = 2; i <= n; i++) {
        if(!isNotPrime[i]) {
            prime[++tot] = i;
            sum1[tot] = sum1[tot - 1] + 1;
            sum2[tot] = sum2[tot - 1] + i;
        }
        for(int j = 1; i * prime[j] <= n; j++) {
            isNotPrime[i * prime[j]] = true;
            if(i % prime[j] == 0) {
                break;
            }
        }
    }
}

inline void initG(int n) {
    // calculate g1[gTot] = \sum 1
    // calculate g2[gTot] = \sum i
    gTot = 0;
    for(int i = 1, r; i <= n; i = r + 1) {
        w[++gTot] = n / i;
        r = n / (n / i);
        if(n / i <= sqrtN) {
            idx1[n / i] = gTot;
        } else {
            idx2[n / (n / i)] = gTot;
        }

        int x = w[gTot];
        g1[gTot] = x - 1;
        g2[gTot] = (ll)x * ((ll)x + 1) / 2 - 1;
        if(r == LIM) {
            break;
        }
    }
    for(int i = 1; i <= tot; i++) {
        for(int j = 1; j <= gTot && prime[i] <= w[j] / prime[i]; j++) {
            int k = w[j] / prime[i] <= sqrtN ? idx1[w[j] / prime[i]] : idx2[n / (w[j] / prime[i])];
            g1[j] -= (g1[k] - sum1[i - 1]);
            g2[j] -= (ll)prime[i] * (g2[k] - sum2[i - 1]);
        }
    }
}

inline int calcSMu(int x, int j, int n) {
    if(prime[j] >= x) {
        return 0;
    }
    int k = x <= sqrtN ? idx1[x] : idx2[n / x];
    int ansMu = -(g1[k] - sum1[j]);
    for(int i = j + 1; i <= tot && prime[i] <= x / prime[i]; i++) {
        ansMu -= calcSMu(x / prime[i], i, n);
    }
    return ansMu;
}

inline ll calcSPhi(int x, int j, int n) {
    if(prime[j] >= x) {
        return 0;
    }
    int k = x <= sqrtN ? idx1[x] : idx2[n / x];
    ll ansPhi = g2[k] - g1[k] - (sum2[j] - sum1[j]);
    for(int i = j + 1; i <= tot && prime[i] <= x / prime[i]; i++) {
        int pie = prime[i];
        for(int e = 1; pie <= x; pie *= prime[i], e++) {
            ansPhi += 1LL * (pie - pie / prime[i]) * (calcSPhi(x / pie, i, n) + (e != 1));
            if(pie > x / prime[i]) {
                break;
            }
        }
    }
    return ansPhi;
}

int main() {
    init(N - 2);

    int t;
    scanf("%d", &t);
    while(t--) {
        int n;
        scanf("%d", &n);
        sqrtN = sqrt(n);
        initG(n);

        if(n == 0) {
            printf("0 0\n");
            continue;
        }

        printf("%lld %d\n", calcSPhi(n, 0, n) + 1, calcSMu(n, 0, n) + 1);
    }
}