码迷,mamicode.com
首页 > 其他好文 > 详细

BZOJ2693 jzptab 莫比乌斯反演

时间:2017-02-27 01:06:50      阅读:202      评论:0      收藏:0      [点我收藏+]

标签:click   tor   names   分享   src   clu   zoj   sed   log   

题意:给定N,M,求$\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {lcm(i,j)} }$,多组询问

题解:

设$F(x,y) = \sum\limits_{i = 1}^x {\sum\limits_{j = 1,\gcd (i,j) = 1}^y {ij} }$  $S(x,y) = \frac{{x(x + 1)}}{2}\frac{{y(y + 1)}}{2}$

我们枚举最大公因数d,则有\[\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {lcm(i,j)} }  = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {\frac{{ij}}{{\gcd (i,j)}}} }  = \frac{{{d^2}F(\left\lfloor {\frac{N}{d}} \right\rfloor ,\left\lfloor {\frac{M}{d}} \right\rfloor )}}{d} = dF(\left\lfloor {\frac{N}{d}} \right\rfloor ,\left\lfloor {\frac{M}{d}} \right\rfloor )\]

反演得到\[F(x,y) = \sum\limits_{i = 1}^{\min (x,y)} {{i^2}\mu (i)S(\left\lfloor {\frac{x}{i}} \right\rfloor ,\left\lfloor {\frac{y}{i}} \right\rfloor )} \]

继续优化\[Ans = \sum\limits_{d = 1}^{\min (N,M)} {[d\sum\limits_{i = 1}^{\min (N,M)} {{i^2}\mu (i)S(\left\lfloor {\frac{x}{{di}}} \right\rfloor ,\left\lfloor {\frac{y}{{di}}} \right\rfloor )} ]}  = \sum\limits_{D = 1}^{\min (N,M)} S (\left\lfloor {\frac{N}{D}} \right\rfloor ,\left\lfloor {\frac{M}{D}} \right\rfloor )\sum\limits_{i|D} {\frac{D}{i}{i^2}\mu (i)} \]

后面那个和式可以通过线性筛预处理出来,后面分块就好

技术分享
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>
using namespace std;
#define ll long long

const int P=20101009;
const int MAXN=10000000+2;
int T;
ll N,M,mu[MAXN],prime[MAXN],f[MAXN];
bool flag[MAXN]; 

void Pre(ll N){
    f[1]=1;
    for(int i=2,j=0;i<=N;i++){
        if(!flag[i]) prime[++j]=i,f[i]=(i-(ll)i*i)%P;
        for(ll k=1;k<=j && i*prime[k]<=N;k++){
            flag[i*prime[k]]=1;
            if(i%prime[k]==0){
                f[i*prime[k]]=(f[i]*prime[k])%P;
                break;
            }
            f[i*prime[k]]=(f[prime[k]]*f[i])%P;
        }
    }
    for(int i=1;i<=N;i++) f[i]=(f[i]+f[i-1])%P;
}

ll Calc(ll x,ll y){
    ll ret1,ret2;
    ret1=((x*(x+1))>>1)%P;
    ret2=((y*(y+1))>>1)%P;
    return (ret1*ret2)%P;
}

int main(){
    Pre(MAXN);

    cin >> T;
    while(T--){
        ll ans=0;
        scanf("%lld %lld",&N,&M);
        if(N<M) swap(N,M);

        for(ll i=1,j;i<=M;i=j+1){
            j=min(N/(N/i),M/(M/i));
            ans=(ans+((f[j]-f[i-1])*Calc(N/j,M/j))%P)%P;
        }

        printf("%lld\n",(ans+P)%P);
    }

    return 0;
}
View Code

 

BZOJ2693 jzptab 莫比乌斯反演

标签:click   tor   names   分享   src   clu   zoj   sed   log   

原文地址:http://www.cnblogs.com/WDZRMPCBIT/p/6464159.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!