285484026

bzoj 3309 DZY Loves Math 莫比乌斯函数反演 整数分块
问题为:然后就可以用线性筛递推以后用整数分块计算答案。#include <bits/stdc++.h>...
扫描右侧二维码阅读全文
26
2019/04

bzoj 3309 DZY Loves Math 莫比乌斯函数反演 整数分块

问题为:$$ans=\sum^{n}_{i=1}\sum^{n}_{j=1}f(gcd(i,j))\ \ \ 1\leq n\leq10^7$$ 对上式作以下推导 \begin{gather*} &=\sum^{n}_{d=1}\sum^{n}_{i=1}\sum^{n}_{j=1}f(d)[gcd(i,j)==d]\\ &=\sum^{n}_{d=1}f(d)\sum^{n}_{i=1}\sum^{n}_{j=1}[gcd(i,j)==d]\\ \end{gather*} 对最后一节进行莫比乌斯反演得 \begin{align*} 设F(n,d)&=\sum^{n}_{i=1}\sum^{n}_{j=1}[d|gcd(i,j)]\\ F(n,d)&=\sum^{\lfloor\frac{n}{d}\rfloor}_{i=1}\sum^{\lfloor\frac{n}{d}\rfloor}_{j=1}[1|gcd(i,j)]\\ F(n,d)&=\lfloor\frac{n}{d}\rfloor\lfloor\frac{n}{d}\rfloor\\ \end{align*} \begin{align*} \sum^{n}_{i=1}\sum^{n}_{j=1}[gcd(i,j)==d]&=\sum^{\lfloor\frac{n}{d}\rfloor}_{x=1}\mu(x)F(n,dx)\\ \sum^{n}_{i=1}\sum^{n}_{j=1}[gcd(i,j)==d]&=\sum^{\lfloor\frac{n}{d}\rfloor}_{x=1}\mu(x)\lfloor\frac{n}{xd}\rfloor\lfloor\frac{n}{xd}\rfloor\\ \end{align*} \begin{align*} ans&=\sum^{n}_{d=1}f(d)\sum^{\lfloor\frac{n}{d}\rfloor}_{x=1}\mu(x)\lfloor\frac{n}{xd}\rfloor\lfloor\frac{n}{xd}\rfloor\\ 设T&=dx并带入上式的中得\\ ans&=\sum^{n}_{T=1}\lfloor\frac{n}{xd}\rfloor\lfloor\frac{n}{xd}\rfloor\sum_{d|T}\mu(\frac{T}{d})f(d)\\ 设g(T)&=\sum_{d|T}\mu(\frac{T}{d})f(d)\\ ans&=\sum^{n}_{T=1}\lfloor\frac{n}{xd}\rfloor\lfloor\frac{n}{xd}\rfloor g(T) \end{align*}

$$设x=p_1^{k_1}……p_n^{k_n}$$

$$ g(x)=\left\{ \begin{aligned} & (-1)^{k+1} &k_1=k_2=……=k_n \\ & 0 &其它 \\ \end{aligned} \right. $$

然后就可以用线性筛递推以后用整数分块计算答案。

#include <bits/stdc++.h>

const int MAXN = 1e7;

using namespace std;

int T,n,m,cnt;
int p[MAXN+5],g[MAXN+5],a[MAXN+5],b[MAXN+5];
bool vis[MAXN+5];

inline void init()
{
    for(int i = 2;i <= MAXN;i++){
        if(!vis[i]){
            p[++cnt] = i;
            a[i] = 1;b[i] = i;
            g[i] = 1;
        }
        for(int j = 1;j <= cnt && p[j]*i <= MAXN;j++){
            vis[p[j]*i] = true;
            if(i%p[j] == 0){
                 a[p[j]*i] = a[i] + 1;
                 b[p[j]*i] = b[i] * p[j];
                 int d = i / b[i];
                 if(d == 1) g[p[j]*i] = 1;
                 else {
                     if(a[p[j]*i] == a[d]) g[p[j]*i] = -g[d];
                     else g[p[j]*i] = 0;
                 }
                 break;
            }
            else {
                a[p[j]*i] = 1;b[p[j]*i] = p[j];
                if(a[i] == 1) g[p[j]*i] = -g[i];
                else g[p[j]*i] = 0;
            }
        }
    }
    for(int i = 2;i < MAXN;i++) g[i] += g[i-1];
}

int main()
{
    init();
    scanf("%d",&T);
    while(T--){
        scanf("%d %d",&n,&m);
        int R;
        long long ans = 0;
        for(int i = 1;i <= min(n,m);i = R + 1)
        {
            R = min(n/(n/i),m/(m/i));
            ans += 1LL*(n/i)*(m/i)*(g[R]-g[i-1]);
        }
        printf("%lld\n",ans);
    }
    return 0;
} 
Last modification:April 26th, 2019 at 11:39 pm
If you think my article is useful to you, please feel free to appreciate

Leave a Comment