285484026

bzoj2693 jzptab 莫比乌斯反演+整除分块
\begin{gather*} &ans=\sum^{n}_{i=1}\sum^{m}_{j=1}\frac{ij...
扫描右侧二维码阅读全文
23
2019/05

bzoj2693 jzptab 莫比乌斯反演+整除分块

$$ans=\sum^{n}_{i=1}\sum^{m}_{j=1}lcm(i,j)$$

\begin{gather*} &ans=\sum^{n}_{i=1}\sum^{m}_{j=1}\frac{ij}{gcd(i,j)}\\ &ans=\sum^{n}_{d=1}\sum^{m}_{i=1}\sum^{n}_{j=1}\frac{ij}{d}[gcd(i,j)==d]\\ &ans=\sum^{n}_{d=1}\frac{1}{d}\sum^{m}_{i=1}i\sum^{n}_{j=1}j[gcd(i,j)==d]\\ \end{gather*} \begin{gather*} 设f(n,d)=\sum^{n}_{i=1}i\sum^{m}_{j=1}j[gcd(i,j)==d]\\ 设F(n,d)=\sum^{n}_{i=1}i\sum^{m}_{j=1}j[d|gcd(i,j)]\\ \end{gather*} \begin{align*} F(n,d)&=d^2\frac{(\lfloor\frac{n}{d}\rfloor+1)*\lfloor\frac{n}{d}\rfloor*(\lfloor\frac{m}{d}\rfloor+1)*\lfloor\frac{m}{d}\rfloor}{4} \\ 由莫比乌&斯反演得f(n,d)=\sum^{\frac{n}{d}}_{x=1}\mu(x)*F(n,dx)\\ f(n,d)&=\sum^{\frac{n}{d}}_{x=1}\mu(x)*d^2*x ^2\frac{(\lfloor\frac{n}{d}\rfloor+1)*\lfloor\frac{n}{d}\rfloor*(\lfloor\frac{m}{d}\rfloor+1)*\lfloor\frac{m}{d}\rfloor}{4}\\ \end{align*} \begin{gather*} ans=\sum^{n}_{d=1}\frac{1}{d}\sum^{\frac{n}{d}}_{x=1}\mu(x)*d^2*x ^2\frac{(\lfloor\frac{n}{dx}\rfloor+1)*\lfloor\frac{n}{dx}\rfloor*(\lfloor\frac{m}{dx}\rfloor+1)*\lfloor\frac{m}{dx}\rfloor}{4}\\ ans=\sum^{n}_{T=1}\frac{(\lfloor\frac{n}{T}\rfloor+1)*\lfloor\frac{n}{T}\rfloor*(\lfloor\frac{m}{T}\rfloor+1)*\lfloor\frac{m}{T}\rfloor}{4}\sum_{d|T}d*\mu(\frac{T}{d})*(\frac{T}{d})^2\\ ans=\sum^{n}_{T=1}\frac{(\lfloor\frac{n}{T}\rfloor+1)*\lfloor\frac{n}{T}\rfloor*(\lfloor\frac{m}{T}\rfloor+1)*\lfloor\frac{m}{T}\rfloor}{4}\sum_{d|T}(\frac{T}{d})*\mu(d)*d^2\\ ans=\sum^{n}_{T=1}\frac{(\lfloor\frac{n}{T}\rfloor+1)*\lfloor\frac{n}{T}\rfloor*(\lfloor\frac{m}{T}\rfloor+1)*\lfloor\frac{m}{T}\rfloor}{4}\sum_{d|T}(\frac{T}{d})*\mu(d)*d^2\\ 设g(T)=\sum_{d|T}(\frac{T}{d})*\mu(d)*d^2\\ g(T)=T\sum_{d|T}\mu(d)*d\\ ans=\sum^{n}_{T=1}\frac{(\lfloor\frac{n}{T}\rfloor+1)*\lfloor\frac{n}{T}\rfloor*(\lfloor\frac{m}{T}\rfloor+1)*\lfloor\frac{m}{T}\rfloor}{4}g(T)\\ \end{gather*}

对于该表达式可以用线性筛预处理出一个前缀和,然后使用整除分块处理每个询问。

代码:

#include <bits/stdc++.h>

const int MAXN = 1e7,MOD = 100000009;

using namespace std;

int T,n,m,cnt,inv;
long long f[MAXN+5],p[MAXN+5];
bool vis[MAXN+5];

inline long long pow_mod(long long a,long long b){
    long long ret = 1;
    while(b){
        if(b&1) ret = ret * a % MOD;
        a = a * a % MOD;
        b >>= 1;
    }
    ret %= MOD;
    return ret;
}

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

int main()
{
    inv = pow_mod(4,99328719);
    init();
    //printf("%d",f[MAXN]);
    scanf("%d",&T);
    while(T--){
        scanf("%d %d",&n,&m);

        long long ret,ans = 0,R;
        for(long long i = 1;i <= min(n,m);i = R + 1)
        {
            R = min(n/(n/i),m/(m/i));
            ret = f[R] - f[i-1];
            ans += 1LL*(n/i)*(n/i+1)%MOD*(m/i)%MOD*(m/i+1)%MOD*ret%MOD*inv%MOD;
            ans %= MOD;
        }
        ans %= MOD;
        ans += MOD;
        ans %= MOD;
        printf("%lld\n",ans);        
    }

    return 0;
}
Last modification:May 23rd, 2019 at 08:26 am
If you think my article is useful to you, please feel free to appreciate

Leave a Comment