隐藏页面特效

数论分块（二）

上部分链接：数论分块（一）

建议没看过数论分块（一）的优先去看第一部分，第二部分题目难度较高。

P1829 [国家集训队] Crash的数字表格 / JZPTAB

给定 n,m，求：

n∑i=1m∑j=1lcm(i,j)(mod20101009)

  • 对于 100% 的数据，保证 1≤n,m≤107。

将式子变形：

n∑i=1m∑j=1lcm(i,j)=n∑i=1m∑j=1ijgcd(i,j)=n∑d=1n∑i=1m∑j=1ijd[gcd(i,j)=d]=n∑d=1⌊nd⌋∑i=1⌊md⌋∑j=1ijd[gcd(i,j)=1]=n∑d=1d⌊nd⌋∑i=1⌊md⌋∑j=1ij∑k∣gcd(i,j)μ(k)=n∑d=1d∑k=1μ(k)⌊nd⌋∑k∣i⌊md⌋∑k∣jij=n∑d=1d∑k=1μ(k)⌊nkd⌋∑i=1⌊mkd⌋∑j=1ijk2=n∑d=1d∑k=1μ(k)k2⌊nkd⌋∑i=1i⌊mkd⌋∑j=1j=n∑d=1d∑k=1μ(k)k2(⌊nkd⌋+1)⌊nkd⌋2·(⌊mkd⌋+1)⌊mkd⌋2

令 T=kd，按 k 枚举改为按 T 枚举。

=n∑d=1d∑k=1μ(k)k2(⌊nkd⌋+1)⌊nkd⌋(⌊mkd⌋+1)⌊mkd⌋4=n∑d=1dn∑T=1,d∣Tμ(Td)(Td)2(⌊nT⌋+1)⌊nT⌋(⌊mT⌋+1)⌊mT⌋4=n∑T=1∑d∣TTμ(Td)(Td)(⌊nT⌋+1)⌊nT⌋(⌊mT⌋+1)⌊mT⌋4=n∑T=1T∑d∣Tμ(d)d(⌊nT⌋+1)⌊nT⌋(⌊mT⌋+1)⌊mT⌋4

记 f(T)=∑d∣Tdμ(d)，根据直觉 f(T) 为积性函数，我们取几个值来看看。

f(1)=1

f(p)=1−p

f(pk)=1−p

f(pq)=1−p−q+pq=(1−p)(1−q)

那么我们可以通过线性筛筛出 T×f(T) 的前缀和。

后面一坨做数论分块，总体时间复杂度 O(n+√n)。

PS: 本题还有杜教筛的做法，可以做到时间复杂度 O(n23)

参考代码

复制#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7 + 10;
const int mod = 20101009;
int n, m, cnt, primes[N];
ll f[N];
bool st[N];

void init()
{
    f[1] = 1;
    for (int i = 2; i <= n; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, f[i] = 1 - i + mod;
        for (int j = 0; primes[j] * i <= n; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0)
            {
                f[primes[j] * i] = f[i];
                break;
            }
            f[primes[j] * i] = f[i] * (1 - primes[j] + mod) % mod;
        }
    }
    for (int i = 1; i <= n; i ++ ) f[i] = (f[i] * i % mod + f[i - 1]) % mod;
}

ll calc(ll n, ll m)
{
    return ((n + 1) * n / 2 % mod) * ((m + 1) * m / 2 % mod) % mod ;
}

int main()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    init();
    ll res = 0;
    for (int l = 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (res += (f[r] - f[l - 1]) * calc(n / l, m / l) % mod) %= mod;
    }
    cout << (res + mod) % mod;
    return 0;
}

P3327 [SDOI2015] 约数个数和

设 d(x) 为 x 的约数个数，给定 n,m，求

n∑i=1m∑j=1d(ij)

对于 100% 的数据，1≤T,n,m≤50000。

n∑i=1m∑j=1d(ij)=n∑i=1m∑j=1∑d∣ij1

设 xy=d，考虑把 d 分解有什么情况，由于所有的 d 值是唯一的，所以每一个 x 和 y 的枚举需要唯一，不妨设 x∣i，y∣j，考虑 xy 什么时候唯一对应 ij 的因数。

可以发现，如果 xy 不能唯一对应 ij 的因数时，说明 xy 有多种拆解方式，而 x∣i，y∣j，因此，如果 x,y 含有共同的因子，这 xy 是可以再分的，因此所以它的充要条件是 [gcd(x,y)=1]，原始可继续改写。

n∑i=1m∑j=1∑d∣ij1=n∑i=1m∑j=1∑x∣i∑y∣j[gcd(x,y)=1]

又有 x,y 的求和是独立的，可以改写。

n∑i=1m∑j=1∑x∣i∑y∣j[gcd(x,y)=1]=n∑x=1m∑y=1n∑x∣im∑y∣j[gcd(x,y)=1]]=n∑x=1m∑y=1⌊nx⌋⌊my⌋[gcd(x,y)=1]

考虑莫比乌斯反演 ∑d∣nμ(d)=[n=1]。

n∑x=1m∑y=1⌊nx⌋⌊my⌋[gcd(x,y)=1]=n∑x=1m∑y=1⌊nx⌋⌊my⌋∑d∣gcd(x,y)μ(d)=n∑d=1μ(d)n∑d∣xm∑d∣y⌊nx⌋⌊my⌋=n∑d=1μ(d)⌊nd⌋∑i=1⌊md⌋∑j=1⌊ndi⌋⌊mdj⌋=n∑d=1μ(d)(⌊nd⌋∑i=1⌊ndi⌋)(⌊md⌋∑j=1⌊mdj⌋)

设 S(n)=n∑i=1⌊ni⌋，则原式可替换为：

n∑d=1μ(d)S(⌊nd⌋)S(⌊md⌋)

预处理一下 μ(n) 的前缀和与 S(n) 的值，进行数论分块即可。

参考代码

复制#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e4 + 10;
int cnt, st[N], primes[N];
int n, m;
ll S[N], mu[N];

void init()
{
    mu[1] = 1;
    for (ll i = 2; i < N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
        for (int j = 0; primes[j] * i < N; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0) break;
            mu[primes[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
    {
        mu[i] += mu[i - 1];
        for (int l = 1, r; l <= i; l = r + 1)
        {
            r = i / (i / l);
            S[i] += 1ll * (r - l + 1) * (i / l);
        }
    }
}

ll calc(int n, int m)
{
    ll res = 0;
    for (int l = 1, r; l <= min(n, m); l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        res += (mu[r] - mu[l - 1]) * S[n / l] * S[m / l];
    }
    return res;
}

void solve()
{
    scanf("%d%d", &n, &m);
    printf("%lld\n", calc(n, m));
}

int main()
{
    init();
    int T;
    scanf("%d", &T);
    while (T -- ) solve();
    return 0;
}

P4466 [国家集训队] 和与积

考虑 a+b∣ab 有怎样的性质，设 d=gcd(a,b),a=id,b=jd，则原式可化为 i+j∣ijd，由于 gcd(i,j)=gcd(i+j,j)=gcd(i,i+j)=1，故可知 i+j∣d。

因此我们的目标是对这样的 d 的个数求和。

n∑i=1i−1∑j=1⌊di+j⌋[gcd(i,j)=1][id≤n]

合法的 d 的上界为 ⌊ni⌋，所以

n∑i=1i−1∑j=1⌊⌊ni⌋i+j⌋[gcd(i,j)=1]=n∑i=1i−1∑j=1⌊ni(i+j)⌋[gcd(i,j)=1]

可以发现 i 的上界不超过 √n，可以缩小枚举范围，并对后面来一发莫反。

√n∑i=1i−1∑j=1⌊ni(i+j)⌋[gcd(i,j)=1]=√n∑i=1i−1∑j=1⌊ni(i+j)⌋∑d∣gcd(i,j)μ(d)=√n∑i=1i−1∑j=1⌊ni(i+j)⌋∑d∣i,d∣jμ(d)=√n∑d=1μ(d)⌊√nd⌋∑i=1i−1∑j=1⌊⌊nid2⌋i+j⌋=√n∑d=1μ(d)⌊√nd⌋∑i=1i−1∑j=1⌊⌊nid2⌋i+j⌋=√n∑d=1μ(d)⌊√nd⌋∑i=12i−1∑j=i+1⌊⌊nid2⌋j⌋

后面可以数论分块，对复杂度积分（舍去常数）有：

∫√n1dx∫x1√nx2ydy→∫√n1k√nxdx→k√n(√n)12→n34

因此复杂度的上界为 O(n34)，可以通过。

加了一些剪枝，398ms 的最优解。

参考代码

复制#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 60000;
int n;
int primes[N], mu[N], st[N], cnt;

void init(int maxn)
{
    mu[1] = 1;
    for (int i = 2; i <= maxn; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
        for (int j = 0; i * primes[j] <= maxn; j ++ )
        {
            st[i * primes[j]] = 1;
            if (i % primes[j] == 0) break;
            mu[i * primes[j]] = -mu[i];
        }
    }
}

int calc(int st, int ed, int mul)
{
    int res = 0;
    ed = min(ed, mul);
    for (int l = st, r; l <= ed; l = r + 1)
    {
        r = min(ed, mul / (mul / l));
        res += (r - l + 1) * (mul / l);
    }
    return res;
}

int main()
{
    cin >> n;
    int upper = sqrt(n);
    init(upper);
    ll res = 0;
    for (int d = 1; d <= upper; d ++ )
    {
        if (!mu[d]) continue;
        for (int i = 1; i * d <= upper; i ++ )
            res += mu[d] * calc(i + 1, (i << 1) - 1, n / (d * d * i));
    }
    cout << res;
    return 0;
}

P5572 [CmdOI2019] 简单的数论题

给出 n,m 求下列式子的值 ：

n∑i=1m∑j=1φ(lcm(i,j)gcd(i,j))mod23333

对于所有测试点， T≤3×104, m≤n≤5×104。

先看看这个式子的真面目。

n∑i=1m∑j=1φ(lcm(i,j)gcd(i,j))=n∑i=1m∑j=1φ(ijgcd2(i,j))=n∑d=1n∑i=1m∑j=1φ(ijd2)[gcd(i,j)=d]=n∑d=1n∑i=1m∑j=1φ(id)φ(id)[gcd(i,j)=d]=n∑d=1⌊nd⌋∑i=1⌊md⌋∑j=1φ(i)φ(j)[gcd(i,j)=1]=n∑d=1⌊nd⌋∑i=1⌊md⌋∑j=1φ(i)φ(j)∑k∣gcd(i,j)μ(k)=n∑k=1μ(k)n∑d=1⌊nd⌋∑k∣i⌊md⌋∑k∣jφ(i)φ(j)=n∑k=1μ(k)n∑d=1⌊ndk⌋∑i=1⌊mdk⌋∑j=1φ(ki)φ(kj)=n∑k=1μ(k)n∑d=1(⌊ndk⌋∑i=1φ(ki))(⌊mdk⌋∑j=1φ(kj))=n∑k=1n∑d=1μ(d)(⌊ndk⌋∑i=1φ(di))(⌊mdk⌋∑j=1φ(dj))

直到这里，步骤和上一题几乎完全一致，我们考虑用 T=dk 换元。

n∑k=1n∑d=1μ(d)(⌊ndk⌋∑i=1φ(di))(⌊mdk⌋∑j=1φ(dj))=n∑T=1∑d∣Tμ(d)(⌊nT⌋∑i=1φ(di))(⌊mT⌋∑j=1φ(dj))

观察式子后两坨满足 f(d,s)=s∑i=1φ(di)，预处理的复杂度满足调和级数，总复杂度 O(nlnn)。

记 g(x,y,z)=z∑T=1∑d∣Tμ(d)(x∑i=1φ(di))(y∑j=1φ(dj))，把答案用端点 ⌊ni⌋，⌊mi⌋ 拆分成若干线段 [l,r]，一个线段的答案贡献为 g(⌊ni⌋,⌊mi⌋,r)−g(⌊ni⌋,⌊mi⌋,l−1)，但是令人悲伤的是，预处理复杂度是 O(n2mlogn) 的，无法接受，没办法直接上数论分块。

如果我们暴力计算式子，则对于每次询问我们的复杂度是 O(nlogn) 的，也无法通过，于是我们想到用根号分治来解决这个问题。

设置阈值 B，对 ⌊nB⌋ 以内的 z 暴力遍历，时间复杂度是 O(nBlognB)，其余的 g(x,y,z) 使用数论分块，将预处理时间分摊为 O(B2nlogn)，总的时间复杂度应为 O(TnBlognB+T√Blogn+B2nlogn)，取 B 取 20 到 100 中任意一个数都可以，取 50 比较快。

参考代码

复制#include<bits/stdc++.h>
using namespace std;
const int N = 5e4 + 10, B = 50, mod = 23333;
int n, m;
int st[N], prime[N], phi[N], mu[N], cnt;
vector<int> factor[N], f[N], g[B + 5][B + 5];

int calc(int x, int y, int z)
{
    int res = 0;
    for (auto t : factor[z])
        (res += mu[t] * f[t][x] % mod * f[t][y] % mod + mod) %= mod;
    return res;
}

void init()
{
    phi[1] = mu[1] = 1;
    for (int i = 2; i < N; i ++ )
    {
        if (!st[i]) prime[cnt ++ ] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 0; prime[j] * i < N; j ++ )
        {
            st[prime[j] * i] = 1;
            if (i % prime[j] == 0)
            {
                phi[prime[j] * i] = phi[i] * prime[j];
                break;
            }
            phi[prime[j] * i] = phi[i] * (prime[j] - 1);
            mu[prime[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
        for (int j = i; j < N; j += i)
            factor[j].push_back(i);
    for (int i = 1; i < N; i ++ )
    {
        f[i].push_back(0);
        for (int j = 1; i * j < N; j ++ )
            f[i].push_back((f[i][j - 1] + phi[i * j]) % mod);
    }
    for (int i = 1; i <= B; i ++ )
    {
        for (int j = 1; j <= B; j ++ )
        {
            g[i][j].push_back(0);
            for (int k = 1; k * max(i, j) < N; k ++ )
                g[i][j].push_back((g[i][j][k - 1] + calc(i, j, k)) % mod);
        }
    }
}

void solve()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1; i * B <= m; i ++ ) (ans += calc(n / i, m / i, i)) %= mod;
    for (int l = m / B + 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (ans += g[n / l][m / l][r] - g[n / l][m / l][l - 1] + mod) %= mod;
    }
    cout << ans << endl;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL), cout.tie(NULL);
    init();
    int T;
    cin >> T;
    while (T -- ) solve();
    return 0;
}

P4240 毒瘤之神的考验

T 次询问，每次给定 n,m，求：

(n∑i=1m∑j=1φ(ij))mod998244353

对于 100% 的数据，1≤T≤104，1≤n,m≤105。

对于 φ(ij) 有（p 为质数）：

φ(ij)=ij∏p∣ijp−1p=ij∏p∣ip−1p∏p∣jp−1p∏p∣i,p∣jp−1p=i∏p∣ip−1pj∏p∣jp−1p∏p∣gcd(i,j)p−1p=φ(i)φ(j)gcd(i,j)φ(gcd(i,j))

一如既往的化式子：

n∑i=1m∑j=1φ(ij)=n∑i=1m∑j=1φ(i)φ(j)gcd(i,j)φ(gcd(i,j))=n∑d=1n∑d∣im∑d∣jφ(i)φ(j)dφ(d)[gcd(i,j)=d]=n∑d=1dφ(d)⌊nd⌋∑i=1⌊md⌋∑j=1φ(id)φ(jd)[gcd(i,j)=1]=n∑d=1dφ(d)⌊nd⌋∑i=1⌊md⌋∑j=1φ(id)φ(jd)∑k∣gcd(i,j)μ(k)=n∑d=1dφ(d)n∑k=1μ(k)⌊nd⌋∑k∣i⌊md⌋∑k∣jφ(id)φ(jd)=n∑d=1dφ(d)n∑k=1μ(k)⌊nkd⌋∑i=1⌊mkd⌋∑j=1φ(idk)φ(jdk)

记 T=kd 有：

n∑T=1∑d∣Tdφ(d)μ(Td)⌊nT⌋∑i=1⌊mT⌋∑j=1φ(Ti)φ(Tj)

这和上一个题也是类似的，记 f(d,s)=s∑i=1φ(di)，用 O(nlnn) 的时间复杂度预处理出。

记 g(x,y,z)=z∑T=1∑d∣Tdφ(d)μ(Td)x∑i=1y∑j=1φ(Ti)φ(Tj)

考虑根号分治，设阈值为 B，爆算 ⌊nB⌋ 以内的 z，预处理出剩下的部分，复杂度和上一题完全一致，甚至式子都是这么的形似，时间复杂度 O(TnBlognB+T√Blogn+B2nlogn)，取 B 为 45 左右即可。

参考代码

复制#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5 + 10, B = 45, mod = 998244353;
int n, m, st[N], primes[N], phi[N], mu[N], inv[N], cnt;
vector<int> factor[N], f[N], g[B + 5][B + 5];

int calc(int x, int y, int z)
{
    int res = 0;
    for (auto t : factor[z])
        (res += 1ll * t * inv[phi[t]] * mu[z / t] % mod) %= mod;
    res = (1ll * res * f[z][x] % mod * f[z][y] % mod + mod) % mod;
    return res;
}

void init()
{
    inv[1] = phi[1] = mu[1] = 1;
    for (int i = 2; i < N; i ++ ) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
    for (int i = 2; i < N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 0; 1ll * primes[j] * i < N; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0)
            {
                phi[primes[j] * i] = phi[i] * primes[j];
                break;
            }
            phi[primes[j] * i] = phi[i] * (primes[j] - 1);
            mu[primes[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
        for (int j = i; j < N; j += i)
            factor[j].push_back(i);
    for (int i = 1; i < N; i ++ )
    {
        f[i].push_back(0);
        for (int j = 1; i * j < N; j ++ )
            f[i].push_back((f[i][j - 1] + phi[i * j]) % mod);
    }
    for (int i = 1; i <= B; i ++ )
    {
        for (int j = 1; j <= B; j ++ )
        {
            g[i][j].push_back(0);
            for (int k = 1; k * max(i, j) < N; k ++ )
                g[i][j].push_back((g[i][j][k - 1] + calc(i, j, k)) % mod);
        }
    }
}

void solve()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1; i * B <= m; i ++ ) (ans += calc(n / i, m / i, i)) %= mod;
    for (int l = m / B + 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (ans += (g[n / l][m / l][r] - g[n / l][m / l][l - 1]) % mod) %= mod;
    }
    cout << (ans + mod) % mod << '\n';
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);cout.tie(NULL);
    init();
    int T;
    cin >> T;
    while (T -- ) solve();
    return 0;
}

我们发现其实不用每次都去枚举质因数，存储进一个系数数组即可，这样时间复杂度少一个 log(n) 级别。

复制#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e5 + 10, B = 45, mod = 998244353;
int n, m, st[N], primes[N], phi[N], mu[N], inv[N], coff[N], cnt;
vector<int> f[N], g[B + 5][B + 5];

void init()
{
    inv[1] = phi[1] = mu[1] = 1;
    for (int i = 2; i < N; i ++ ) inv[i] = 1ll * (mod - mod / i) * inv[mod % i] % mod;
    for (int i = 2; i < N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i, phi[i] = i - 1, mu[i] = -1;
        for (int j = 0; primes[j] * i < N; j ++ )
        {
            st[primes[j] * i] = 1;
            if (i % primes[j] == 0)
            {
                phi[primes[j] * i] = phi[i] * primes[j];
                break;
            }
            phi[primes[j] * i] = phi[i] * (primes[j] - 1);
            mu[primes[j] * i] = -mu[i];
        }
    }
    for (int i = 1; i < N; i ++ )
        for (int j = i; j < N; j += i)
            (coff[j] += (1ll * mu[j / i] * i * inv[phi[i]] % mod + mod) % mod) %= mod;
    for (int i = 1; i < N; i ++ )
    {
        f[i].push_back(0);
        for (int j = 1; i * j < N; j ++ )
            f[i].push_back((f[i][j - 1] + phi[i * j]) % mod);
    }
    for (int i = 1; i <= B; i ++ )
    {
        for (int j = 1; j <= B; j ++ )
        {
            g[i][j].push_back(0);
            for (int k = 1; k * max(i, j) < N; k ++ )
                g[i][j].push_back((g[i][j][k - 1] + 1ll * coff[k] * f[k][i] % mod * f[k][j] % mod) % mod);
        }
    }
}

void solve()
{
    cin >> n >> m;
    if (n > m) swap(n, m);
    int ans = 0;
    for (int i = 1; i * B <= m; i ++ ) (ans += 1ll * coff[i] * f[i][n / i] % mod * f[i][m / i] % mod) %= mod;
    for (int l = m / B + 1, r; l <= n; l = r + 1)
    {
        r = min(n / (n / l), m / (m / l));
        (ans += ((g[n / l][m / l][r] - g[n / l][m / l][l - 1]) % mod + mod) % mod) %= mod;
    }
    cout << ans << '\n';
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(NULL);cout.tie(NULL);
    init();
    int T;
    cin >> T;
    while (T -- ) solve();
    return 0;
}

__EOF__

本文作者：Yip.Chip
本文链接：https://www.cnblogs.com/YipChipqwq/p/-/shu-lun-fen-kuai2.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！
posted @   YipChip  阅读(35)  评论(0)  编辑  收藏  举报
相关博文：
· 数论分块（一）
· 类欧几里得算法
· 数论函数从入门到进门
· 数论\to 模范
· 积性函数学习笔记
阅读排行：
· 单线程的Redis速度为什么快？
· 展开说说关于C#中ORM框架的用法！
· Pantheons：用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效
点击右上角即可分享
微信分享提示