P2522 HAOI2011 Problem b [莫比乌斯反演,数论分块]

最新推荐文章于 2022-02-26 16:04:00 发布

phython96

最新推荐文章于 2022-02-26 16:04:00 发布

阅读量280

点赞数

CC 4.0 BY-SA版权

分类专栏： ACM-ICPC训练题解数论与组合数学系列

本文链接：https://blog.youkuaiyun.com/weixin_37517391/article/details/83757798

ACM-ICPC训练题解同时被 2 个专栏收录

241 篇文章

订阅专栏

数论与组合数学系列

21 篇文章

订阅专栏

P2522 HAOI2011

题意

对于给出的n个询问，每次求有多少个数对 $(x, y)$ ，满足 $a \leq x \leq b$ ， $c \leq y \leq d$ ，且 $g c d (x, y) = k$ ， $g c d (x, y)$ 函数为 $x$ 和 $y$ 的最大公约数.

题解

即求式子 $∑x=ab∑y=cd[gcd(x,y)=k]\sum_{x=a}^b\sum_{y=c}^d[gcd(x,y)=k]$ .

记 $f(n,m)=∑x=1n∑y=1m[gcd(x,y)=k]f(n,m)=\sum_{x=1}^n\sum_{y=1}^m[gcd(x,y)=k]$

根据二维前缀和公式,可以将式子转换成:

$∑x=ab∑y=cd[gcd(x,y)=k]=f(b,d)+f(a−1,c−1)−f(a−1,d)−f(b,c−1)\sum_{x=a}^b\sum_{y=c}^d[gcd(x,y)=k]=f(b,d)+f(a-1,c-1)-f(a-1,d)-f(b,c-1)$

因此我们只要能得到 $f (n, m)$ 的计算方法即可.

求 $f (n, m)$ 的套路非常明显:莫比比乌斯反演

由于 $∑k∣d∑x=1n∑y=1m[gcd(x,y)=d]=⌊nk⌋⌊mk⌋\sum_{k|d}\sum_{x=1}^n\sum_{y=1}^m[gcd(x,y)=d] = \lfloor \frac{n}{k} \rfloor \lfloor \frac{m}{k} \rfloor$ .

反演得到 $f(n,m)=∑k∣dμ(dk)⌊nd⌋⌊md⌋f(n,m)=\sum_{k|d}\mu(\frac{d}{k})\lfloor \frac{n}{d} \rfloor \lfloor \frac{m}{d} \rfloor$

令 $t = d / k$ .

$f(n,m)=∑t=1n/dμ(t)⌊nkt⌋⌊mkt⌋f(n,m)=\sum_{t=1}^{n/d}\mu(t)\lfloor \frac{n}{kt} \rfloor \lfloor \frac{m}{kt} \rfloor$

对上式子进行分块计算,可以将时间复杂度从 $O (n)$ 降至 $O(n)O(\sqrt{n})$ .

总结

对于形如 $f(x)=∑x∣dμ(dx)g(⌊nd⌋)f(x)=\sum_{x|d}\mu(\frac{d}{x})g(\lfloor \frac{n}{d} \rfloor)$ 这样的式子,我们都可以用 $t = d / x$ 代换后数论分块进行加速.

$f(x)=∑t=1n/xμ(t)g(⌊nxt⌋)f(x)=\sum_{t=1}^{n/x}\mu(t)g(\lfloor \frac{n}{xt} \rfloor)$

时间复杂度从 $O (n)$ 降至 $O(n)O(\sqrt{n})$ .

代码

#include <iostream>
#include <algorithm>
#include <cstring>
#define pr(x) std::cout << #x << ':' << x << std::endl
#define rep(i,a,b) for(int i = a;i <= b;++i)

const int N = 50000;

int prime[N+10],zhi[N+10],mu[N+10],pcnt;

void sieve() {
    zhi[1] = mu[1] = 1;
    for(int i = 2;i <= N;++i) {
        if(!zhi[i]) {
            mu[i] = -1;
            prime[pcnt++] = i;
        }
        for(int j = 0;j < pcnt && prime[j]*i <= N;++j) {
            zhi[i*prime[j]] = 1;
            if(i % prime[j] == 0) {
                mu[i*prime[j]] = 0;
                break;
            }
            else 
                mu[i*prime[j]] = -mu[i];
        }
    }
    for(int i = 1;i <= N;++i) mu[i] += mu[i-1];
}
int a,b,c,d,k,T;
int calc(int n,int m) {
    int ans = 0;
    int lim = std::min(n/k,m/k);
    for(int i = 1,nx1,nx2,nxt;i <= lim;i=nxt+1) {
        nx1 = n/(n/i);
        nx2 = m/(m/i);
        nxt = nx1>nx2?nx2:nx1;
        ans += (mu[nxt]-mu[i-1])*(n/i/k)*(m/i/k);
    }
    return ans;
}
int main() {
    std::ios::sync_with_stdio(false);
    sieve();
    std::cin >> T;
    while(T--) {
        std::cin >> a >> b >> c >> d >> k;
        int ans = calc(b,d)+calc(a-1,c-1)-calc(a-1,d)-calc(b,c-1);
        std::cout << ans << std::endl;
    }
    return 0;
}