数论分块学习笔记

准备开始复习莫比乌斯反演,杜教筛这一部分,先复习一下数论分块

0.随便说说

数论分块可以快速计算一些含有除法向下取整的和式,即形如 ∑ i = 1 n f ( i ) g ( ⌊ n i ⌋ ) \sum_{i=1}^{n}f(i)g(\lfloor\frac{n}{i}\rfloor) i=1nf(i)g(⌊in⌋)的式子。
利用的原理是 ⌊ n i ⌋ \lfloor\frac{n}{i}\rfloor in的不同的值不超过 2 n 2\sqrt{n} 2n 个。
当我们可以在 O ( 1 ) O(1) O(1)的时间快速处理出 ∑ i = l r f ( i ) \sum_{i=l}^{r}f(i) i=lrf(i)或提前预处理出 f ( x ) f(x) f(x)的前缀和时,上述式子可在 O ( n ) O(\sqrt{n}) O(n )的时间计算出来。

1.代码实现

怎么找每个块是个问题,有个结论:
设块的左端点为 ⌊ n l ⌋ \lfloor\frac{n}{l}\rfloor ln,右端点为 ⌊ n r ⌋ \lfloor\frac{n}{r}\rfloor rn,则 r = ⌊ n ⌊ n l ⌋ ⌋ r=\lfloor\frac{n}{\lfloor\frac{n}{l}\rfloor}\rfloor r=lnn

证明也挺好证的 设 k = ⌊ n i ⌋ , k=\lfloor\frac{n}{i}\rfloor, k=in, k ≤ n i , k\le \frac{n}{i}, kin,
因此 ⌊ n k ⌋ ≥ ⌊ n n i ⌋ = i \lfloor\frac{n}{k}\rfloor\ge\lfloor\frac{n}{\frac{n}{i}}\rfloor=i kninn=i,即 i ≤ ⌊ n k ⌋ = ⌊ n ⌊ n i ⌋ ⌋ i\le\lfloor\frac{n}{k}\rfloor=\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor ikn=inn因此右端点 r = i m a x = ⌊ n ⌊ n i ⌋ ⌋ r=i_{max}=\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor r=imax=inn

因此每个块为 i = l i=l i=l i = ⌊ n ⌊ n i ⌋ ⌋ i=\lfloor\frac{n}{\lfloor\frac{n}{i}\rfloor}\rfloor i=inn

在这里插入图片描述

2.例题

先顺手把 O I    W i k i OI\,\,Wiki OIWiki上的三个例题做了吧

UVA11526 H(n)

题面
用洛谷的题面了,这题就是求 ∑ i = 1 n ⌊ n i ⌋ \sum_{i=1}^{n}\lfloor\frac{n}{i}\rfloor i=1nin,就是板子题,相当于 f ( x ) = 1 , g ( n / i ) = n / i f(x)=1,g(n/i)=n/i f(x)=1,g(n/i)=n/i。注意如果 n = 2147483647 n=2147483647 n=2147483647,最后一次 n x t + 1 nxt+1 nxt+1会爆 i n t int int U V A UVA UVA神奇 o j oj oj会报 R E RE RE,所以就都开 l o n g    l o n g long\,\,long longlong就行,时间复杂度 O ( T n ) O(T\sqrt{n}) O(Tn )

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int t,n;
int main(){
    cin>>t;
    while(t--){
        cin>>n;
        ll nxt=0,ans=0;
        for(ll i=1;i<=n;i=nxt+1){
            nxt=n/(n/i);
            ans+=(nxt-i+1)*(n/i);
        }
        cout<<ans<<endl;
    }
}

P2261 [CQOI2007] 余数求和

题面
Solution
O I OI OI时期的博客有这道题,题解挂链接了,随手一推就是这样,减号左边是 n k , nk, nk,右侧用数论分块做

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll n,k,ans,nxt,sum;
inline ll s(ll n){
    return n*(n+1)/2;
}
int main(){
    cin>>n>>k;
    ans=n*k;
    for(int i=1;i<=min(n,k);i=nxt+1){
        nxt=min(k/(k/i),n);
        sum+=(s(nxt)-s(i-1))*(k/i);
    }
    cout<<ans-sum<<endl;
}
P3455 [POI2007] ZAP-Queries

题面

最后用个二维数论分块就可
闲得没事 这题想多打几个空格

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 5e4 + 100;
int cnt, pri[maxn], mu[maxn], s_mu[maxn];
bool vis[maxn];
inline void read(int &x) {
    int s = 0, w = 1; char ch = getchar();
    while (ch < '0' || ch > '9') {if(ch == '-') w = -1; ch = getchar(); }
    while (ch >= '0' && ch <= '9') {s = (s << 3) + (s << 1) + (ch & 15); ch = getchar(); }
    x = s * w;
}
void setup() {
    mu[1] = 1;
    for (int i = 2; i <= maxn - 100; i++) {
        if(!vis[i]) pri[++cnt] = i, mu[i] = -1;
        for (int j = 1; j <=cnt && i * pri[j] <= maxn - 100; j++) {
            vis[i * pri[j]] = true;
            if(i % pri[j]) mu[i * pri[j]] = -mu[i];
            else {
                mu[i * pri[j]] = 0;
                break;
            }
        }
    }
    for (int i = 1; i <= maxn - 100; i++) s_mu[i] = s_mu[i - 1] + mu[i];
}
int T, a, b, d, nxt;
int main() {
    setup();
    read(T);
    while(T--) {
        read(a), read(b), read(d);
        a /= d, b /= d;
        if(a > b) a ^= b, b ^= a, a ^= b;
        ll ans = 0;
        for (int i = 1; i <= a; i = nxt + 1) {
            nxt = min(a / (a / i), b / (b / i));
            ans += 1ll * (s_mu[nxt] - s_mu[i - 1]) * (a / i) * (b / i);
        }
        printf("%lld\n", ans);
    }
}

放两道套路题,这一类题都是以 ∑ i = 1 n ∑ j = 1 n f ( g c d ( i , j ) ) \sum_{i=1}^{n}\sum_{j=1}^{n}f(gcd(i,j)) i=1nj=1nf(gcd(i,j))形式的,我们要将 g c d ( i , j ) gcd(i,j) gcd(i,j)提出来,作如下变化:

然后求出 f ( x ) , ϕ ( x ) f(x),\phi(x) f(x),ϕ(x)的前缀和并应用数论分块解决。
后面两题的数据范围不一样,第一题是 n ≤ 1 e 7 n\le 1e7 n1e7,第二题是 n ≤ 1 e 9 n\le 1e9 n1e9,前者可以应用朴素的筛法求出欧拉函数前缀和,后者要用到杜教筛进行求解。

bzoj4804 欧拉心算

题面
请添加图片描述
预处理出欧拉函数前缀和,应用数论分块即可。
写错了,最后是 s u m ϕ ( n ) sum\phi(n) sumϕ(n)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll maxn = 1e7 + 10;
int cnt, pri[maxn], euc[maxn];
ll s[maxn];
bool vis[maxn];
void setup(){
    euc[1] = 1;
    for (int i = 2; i <= maxn - 10; i++) {
        if (!vis[i]) pri[++cnt] = i, euc[i] = i - 1;
        for (int j = 1; j <= cnt && i * pri[j] <= maxn - 10; j++) {
            vis[i * pri[j]] = true;
            if (i % pri[j]) euc[i * pri[j]] = euc[i] * (pri[j] - 1);
            else {
                euc[i * pri[j]] = pri[j] * euc[i];
                break;
            }
        }
    }
    for (int i = 1; i <= maxn - 10; i++) s[i] = s[i - 1] + euc[i];
}
int T, n;
int main() {
    setup();
    scanf("%d", &T);
    while (T--) {
        scanf("%d", &n);
        int nxt = 0; ll ans = 0;
        for (int i = 1; i <= n; i = nxt + 1) {
            nxt = n / (n / i);
            ans += (s[nxt] - s[i - 1]) * s[n / i];
        }
        printf("%lld\n", 2ll * ans - s[n]);
    }
}

HDU7325 GCD Magic

题面
考场上考到得,打了 150 150 150行, M L E MLE MLE了一次, W A WA WA了一次,最后过了
M L E MLE MLE是因为对 s b H D U O J sbHDUOJ sbHDUOJ不信任,预处理的 1 e 7 1e7 1e7的欧拉函数前缀和,后来改成了 2 e 6 2e6 2e6 W A WA WA了,看了一会儿发现是因为杜教筛欧拉函数前缀和没模 m o d mod mod导致后面计算答案时候爆 l o n g    l o n g long\,\,long longlong了,改过来交一发对了,太不容易了,这要再错真不知道要 d e b u g debug debug到哪年去。
思路如下,就不打公式了,太多了,手写了。
请添加图片描述
上代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 998244353;
const int maxn = 2000010;
inline void read(int &x)
{
    int s = 0, w = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9')
    {
        if (ch == '-')
            w = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9')
    {
        s = (s << 3) + (s << 1) + (ch & 15);
        ch = getchar();
    }
    x = s * w;
}
ll ans, s[maxn], frac[20], fracinv[20];
int t, n, k, lst;
ll pri[maxn], euc[maxn], cur, mu[maxn], sum_mu[maxn];
bool vis[maxn];
map<ll, ll> mp_mu;

ll S_mu(ll x)
{
    if (x < maxn)
        return sum_mu[x];
    if (mp_mu[x])
        return mp_mu[x];
    ll ret = (ll)1;
    for (ll i = 2, j; i <= x; i = j + 1)
    {
        j = x / (x / i);
        ret -= S_mu(x / i) * (j - i + 1);
    }
    return mp_mu[x] = ret % mod;
}

ll S_phi(ll x)
{
    ll ret = (ll)0;
    ll j;
    for (ll i = 1; i <= x; i = j + 1)
    {
        j = x / (x / i);
        ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
    }
    return ((ret - 1) / 2 + 1) % mod;
}
inline ll qpow(ll x, ll y)
{
    ll re = 1LL;
    while (y)
    {
        if (y & 1)
            (re *= x) %= mod;
        (x *= x) %= mod;
        y >>= 1;
    }
    return re;
}
inline ll inv(ll x)
{
    return qpow(x, mod - 2);
}
void setup()
{
    frac[0] = fracinv[0] = 1LL;
    for (int i = 1; i <= 15; i++)
        frac[i] = 1ll * frac[i - 1] * i % mod;
    fracinv[15] = inv(frac[15]);
    for (int i = 14; i; i--)
        fracinv[i] = 1ll * fracinv[i + 1] * (i + 1) % mod;
    euc[1] = 1;
    mu[1] = 1;
    for (int i = 2; i < maxn; i++)
    {
        if (!vis[i])
        {
            pri[++cur] = i;
            mu[i] = -1;
            euc[i] = i - 1;
        }
        for (int j = 1; j <= cur && i * pri[j] < maxn; j++)
        {
            vis[i * pri[j]] = true;
            if (i % pri[j])
                mu[i * pri[j]] = -mu[i], euc[i * pri[j]] = euc[i] * euc[pri[j]];
            else
            {
                mu[i * pri[j]] = 0;
                euc[i * pri[j]] = euc[i] * pri[j];
                break;
            }
        }
    }

    for (int i = 1; i < maxn; i++)
        s[i] = (1ll * s[i - 1] + euc[i]) % mod;
    for (int i = 1; i < maxn; i++)
        sum_mu[i] = (1ll * sum_mu[i - 1] + mu[i]) % mod;
}
inline ll C(ll n, ll m)
{
    if (n < m)
        return 0;
    return 1ll * frac[n] * fracinv[m] % mod * fracinv[n - m] % mod;
}
inline ll seq(ll n, ll k)
{
    if (k == 0)
        return n;
    ll a1 = (1LL << k) % mod, q = (1LL << k) % mod;
    return 1LL * a1 * ((qpow(q, n) - 1 + mod) % mod) % mod * inv((q - 1 + mod) % mod) % mod;
}
inline ll sg(ll n, ll k)
{
    if (n == 0)
        return 0LL;
    ll ans = 0, flag = 1;
    for (int i = 0; i <= k; i++)
    {
        ans = (1ll * ans + 1ll * ((flag + mod) % mod) * C(k, i) % mod * seq(n, k - i) % mod) % mod;
        flag *= (-1ll);
    }
    return ans % mod;
}
int main()
{
    setup();
    read(t);
    while (t--)
    {
        read(n), read(k);
        ans = 0ll;
        lst = 0;
        for (int i = 1; i <= n; i = lst + 1)
        {
            lst = n / (n / i);
            ll snow = 0;
            if (n / i >= maxn - 10)
                snow = S_phi(n / i);
            else
                snow = s[n / i];
            ans = (1ll * ans + (1ll * sg(lst, k) - sg(i - 1, k) + mod) % mod * snow % mod) % mod;
        }
        printf("%lld\n", (2LL * ans % mod - sg(n, k) + mod) % mod);
    }
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值