杜教篩

問題描述

求 \(S(n)=\sum\limits_{i=1}^{n}f(i)\)，其中 \(f(i)\) 為積性函式。

\(1\le n\le2^{31}\)。

演算法解決

顯然，線性篩是 \(\mathcal O(n)\) 的，過不了。

考慮推數學式子。

構造積性函式 \(h,g\)，使得 \(h=f*g\)。

有：

\[\sum\limits_{i=1}^{n}h(i)=\sum\limits_{i=1}^{n}\sum\limits_{d|i}f(\frac{i}{d})g(d)\\ =\sum\limits_{d=1}^{n}g(d)\sum\limits_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}f(i)\\ =\sum\limits_{d=1}^{n}g(d)S(\left\lfloor\frac{n}{d}\right\rfloor) \]

因為我們要算的是 \(S(n)\)

，那把這項單獨拎出來不就得了。

所以，有：

\[S(n)=\frac{\sum\limits_{i=1}^{n}h(i)-\sum\limits_{d=2}^{n}g(d)S(\left\lfloor\frac{n}{d}\right\rfloor)}{g(1)} \]

於是我們只要能快速算出 \(h\) 的字首和就可以了。

什麼叫快速算出 \(h\) 的字首和呢？

比如說 \(\mu * I = \epsilon\) 的字首和就非常的好算。

\(\varphi * I=id\) 的字首和也非常的好算。

實戰演練

\(\mu\) 的字首和

因為 \(\mu * I=\epsilon\)，那麼取 \(f=\mu,g=I,f*g=\epsilon\)

。

inline ll GetSumu(int n)
{
    if (n <= N)
        return sumu[n]; // sumu是提前篩好的字首和
    if (Smu[n])
        return Smu[n]; // 記憶化
    ll ret = 1ll;      // 單位元的字首和就是 1
    for (int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        ret -= (r - l + 1) * GetSumu(n / l);
        // (r - l + 1) 就是 I 在 [l, r] 的和
    }
    return Smu[n] = ret; // 記憶化
}
 

\(\varphi\) 的字首和

因為 \(\varphi * I =id\)，那麼取 \(f=\varphi,g=I,f*g=id\)。

inline ll GetSphi(int n)
{
    if (n <= N)
        return sump[n]; // 提前篩好的
    if (Sphi[n])
        return Sphi[n];             // 記憶化
    ll ret = 1ll * n * (n + 1) / 2; // f * g = id 的字首和
    for (int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        ret -= (r - l + 1) * GetSphi(n / l);
        // 同上，因為兩個的 g 都是 I
    }
    return Sphi[n] = ret; // 記憶化
}

P3768 簡單的數學題

求：

\[\sum\limits_{i=1}^{n}\varphi(i)\cdot i \]

令 \(f= \varphi \cdot id,g=id\)，考慮迪利克雷卷積的形式得到 \((f*g)(n)=\sum\limits_{d|n}(\varphi(d)\cdot d)\cdot (\frac{n}{d})=n\sum\limits_{d|n}\varphi(d)=n^2\)，即 \((f * g)(i)=i^2\)。

這樣就可以快速求得 \((f*g)(i)\) 的字首和，為 \(\frac{n(n+1)(2n+1)}{6}\)。

杜教篩即可。

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

#define int unsigned long long

int read()
{
    char c = getchar();
    int x = 0, f = 1;
    while (c < '0' || c > '9')
    {
        if (c == '-')
            f = -f;
        c = getchar();
    }
    while (c >= '0' && c <= '9')
        x = x * 10 + c - '0', c = getchar();
    return x * f;
}

const int N = 5e6 + 5;

#define g(x) (((((x % mod) * ((x + 1) % mod) / 2) % mod) * ((x % mod) * ((x + 1) % mod) / 2 % mod)) % mod)

#define g2(x) (((x % mod) * ((x + 1) % mod) % mod * ((2 * x + 1) % mod) % mod * inv6 % mod) % mod)

int n, phi[N], pr[N], f[N], top, mod, inv6, inv2;

bool vis[N];

map<int, int> mp;

void init()
{
    phi[1] = f[1] = 1;
    for (int i = 2; i <= n; ++i)
    {
        if (!vis[i])
            pr[++top] = i, phi[i] = i - 1;
        for (int j = 1; j <= top; ++j)
        {
            if (i * pr[j] > n)
                break;
            vis[i * pr[j]] = 1;
            if (i % pr[j] == 0)
            {
                phi[i * pr[j]] = phi[i] * pr[j];
                break;
            }
            phi[i * pr[j]] = phi[i] * phi[pr[j]];
        }
        f[i] = f[i - 1] + ((1ll * phi[i] * i) % mod * i) % mod;
        f[i] %= mod;
    }
}

int fpow(int x, int k)
{
    int base = x, ans = 1;
    while (k)
    {
        if (k & 1)
            ans *= base, ans %= mod;
        base *= base, base %= mod;
        k >>= 1;
    }
    return ans;
}

int phi_sum(int x)
{
    if (x <= n)
        return f[x];
    if (mp[x])
        return mp[x];
    int sumId = g(x % mod) % mod, ans = 0;
    int l, r;
    for (l = 2; l <= x; l = r + 1)
    {
        r = (x / (x / l));
        ans += ((1ll * ((g2(r) - g2((l - 1)) + mod) % mod) * phi_sum(x / l)) % mod);
        ans %= mod;
    }
    return mp[x] = ((sumId - ans + mod) % mod);
}

int solve1(int x)
{
    int l, r, ans = 0;
    for (l = 1; l <= x; l = r + 1)
    {
        r = (x / (x / l));
        int ret = g(x / l) % mod;
        ret = ret * ((phi_sum(r) - phi_sum(l - 1) + mod) % mod), ret %= mod;
        ans += ret, ans %= mod;
    }
    return (ans + mod) % mod;
}

signed main()
{
    mod = read();
    inv6 = fpow(6, mod - 2);
    int x = read();
    n = 5000000;
    init();
    printf("%lld", solve1(x));
    return 0;
}

本文來自部落格園，作者：蒟蒻orz，轉載請註明原文連結：https://www.cnblogs.com/orzz/p/15727115.html