Luogu5176 公约数 莫比乌斯反演、线性筛

传送门


好像是我们联考时候的题目?

一个结论:\(\gcd(ij,ik,jk) \times \gcd(i,j,k) = \gcd(i,j) \times \gcd(i,k) \times \gcd(j,k)\),证明:由于\(\gcd\)是积性函数,所以我们分成每个质因子考虑。假设对于质数\(p\)\(i,j,k\)的素数唯一分解中\(p\)的质数为\(a \geq b \geq c\),那么如果只考虑\(p\)这一个质因子,\(\gcd(ij,ik,jk) = p^{b+c} , \gcd(i,j,k) = p^c , \gcd(i,j) = p^b , \gcd(i,k) = p^c , \gcd(j,k) = p^c\),那么两边的乘积都是\(p^{b+2c}\)

那么稍微化简一下,原式就等于\(\sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{k=1}^p \gcd(i,j)^2 + \gcd(i,k)^2 + \gcd(j,k)^2 = pf(n,m)+mf(n,p)+nf(m,p)\),其中\(f(n,m) = \sum\limits_{i=1}^n \sum\limits_{j=1}^m gcd(i,j)^2\)

接下来求\(f(n,m)\)(默认\(n \leq m\)):

\(\begin{align*} f(n,m) & = \sum\limits_{i=1}^n \sum\limits_{j=1}^m \gcd(i,j)^2 \\ & = \sum\limits_{p=1}^n p^2 \sum\limits_{i=1}^{\frac{n}{p}} \sum\limits_{j=1}^{\frac{m}{p}} \sum\limits_{d | i , d | j} \mu(d) = \sum\limits_{p=1}^n p^2 \sum\limits_{d = 1}^{\frac{n}{p}} \mu(d) \frac{n}{dp} \frac{m}{dp} = \sum\limits_{T = 1}^n \frac{n}{T} \frac{m}{T} \sum\limits_{p | T} p^2 \mu(\frac{T}{p}) \end{align*}\)

可以发现\(\frac{n}{T} \frac{m}{T}\)可以数论分块,那么我们需要求\(g(T) = \sum\limits_{p | T} p^2 \mu(\frac{T}{p})\)的前缀和。

如果数据范围很大\(id^2 * \mu\)可以杜教筛,而\(n,m,p \leq 2 \times 10^7\)显然是可以线性筛的。考虑当前的数\(T\)乘上一个质数\(p\)之后\(g(Tp)\)相比\(g(T)\)的变化。

如果\(p \not\mid T\),那么\(g(T)\)\(g(Tp)\)中的贡献是\(-1\),而把\(p\)加入到\(g(T)\)中,\(p^2g(T)\)\(g(Tp)\)的贡献是\(1\),所以\(g(Tp) = (p^2-1)g(T)\)

如果\(p \mid T\),那么\(g(T)\)\(g(Tp)\)中的贡献就是\(0\),即\(g(Tp) = p^2g(T)\)

那么可以线性筛出所有的\(g(T)\)。复杂度\(O(n + T\sqrt{n})\)

#include<bits/stdc++.h>
//this code is written by Itst
using namespace std;

const int _ = 2e7 + 3 , MOD = 1e9 + 7;
int prm[_] , val[_] , cnt;
bool nprm[_];

void init(){
    val[1] = 1;
    for(int i = 2 ; i <= 2e7 ; ++i){
        if(!nprm[i]){
            prm[++cnt] = i;
            val[i] = (1ll * i * i - 1) % MOD;
        }
        for(int j = 1 ; prm[j] * i <= 2e7 ; ++j){
            nprm[prm[j] * i] = 1;
            if(i % prm[j] == 0){
                val[i * prm[j]] = 1ll * val[i] * prm[j] % MOD * prm[j] % MOD;
                break;
            }
            val[i * prm[j]] = (1ll * val[i] * prm[j] % MOD * prm[j] % MOD - val[i] + MOD) % MOD;
        }
    }
    for(int i = 2 ; i <= 2e7 ; ++i)
        val[i] = (val[i] + val[i - 1]) % MOD;
}

long long calc(int x , int y){
    int sum = 0;
    for(int i = 1 , pi ; i <= x && i <= y ; i = pi + 1){
        pi = min(x / (x / i) , y / (y / i));
        sum = (sum + 1ll * (val[pi] - val[i - 1] + MOD) * (x / i) % MOD * (y / i)) % MOD;
    }
    return sum;
}

int main(){
#ifndef ONLINE_JUDGE
    freopen("in","r",stdin);
    //freopen("out","w",stdout);
#endif
    init(); int T;
    for(scanf("%d" , &T) ; T ; --T){
        int N , M , P;
        scanf("%d %d %d" , &N , &M , &P);
        printf("%lld\n" , (calc(N , M) * P + calc(M , P) * N + calc(N , P) * M) % MOD);
    }
    return 0;
}

转载于:https://www.cnblogs.com/Itst/p/10770677.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值