BZOJ 2693: jzptab(莫比乌斯反演)

最新推荐文章于 2019-06-25 21:12:00 发布

weixin_30548917

最新推荐文章于 2019-06-25 21:12:00 发布

阅读量122

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/sdfzsyq/p/10279408.html

本文介绍了一种解决数学问题的算法，即求解从1到n和从1到m的所有数对的最小公倍数之和。通过将问题转化为求最大公约数的乘积和，再利用数论分块和积性函数性质，最终简化为等差数列求和问题。代码使用C++实现，包含预处理和主要计算部分。

传送门

解题思路

\[ ans=\sum\limits_{i=1}^n\sum\limits_{j=1}^mlcm(i,j) \]
\[ ans=\sum\limits_{i=1}^n\sum\limits_{j=1}^m\frac{i*j}{gcd(i,j)} \]
\[ ans=\sum\limits_{d=1}^nd\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{m}{d}}i*j*[gcd(i,j)=1] \]
\[ ans=\sum\limits_{d=1}^nd\sum\limits_{d'=1}^{\frac{n}{d}}\mu(d')*d'^2\sum\limits_{i=1}^{\frac{n}{dd'}}\sum\limits_{j=1}^{\frac{n}{dd'}}i*j \]
发现后面是一个等差数列求和，记为\(sum\)
\[ ans=\sum\limits_{T=1}^nsum(n/T,m/T)\sum\limits_{d=1}^n\mu(d)*d^2*\frac{T}{d} \]
前面的可以数论分块，主要是求后面的。注意\(d\)和\(\frac{T}{d}\)不能直接约分，因为后者为向下取整。发现后面的东西其实是一个积性函数，可以晒出来的。设后面的为\(f\)，首先\(f(p)=p-p*p\)，然后考虑\(f(i*p)\)，如果\(i\)中有\(p\)这个因子，那么\(p\)对\(i\)没有新的约数的贡献，因为新的约数中一定含有\(p^2\)这个因子，而对应的\(\mu\)为\(0\)，\(p\)对\(i\)的贡献为\(i\)中\(\frac{i}{d}\)的\(i\)扩大，所以这种情况下\(f(i*p)=f(i)*p\)。否则\(f(i*p)=f(i)*f(p)\)。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>

using namespace std;
typedef long long LL;
const int N=1e7+5;
const int MOD=1e8+9;

inline int rd(){
    int x=0,f=1;char ch=getchar();
    while(!isdigit(ch)) f=ch=='-'?0:1,ch=getchar();
    while(isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
    return f?x:-x;
}

inline int min(int x,int y){return x<y?x:y;}

int T,n,m,prime[N/10],cnt,f[N];
bool vis[N];

inline void prework(){
    f[1]=1;int lim=1e7;
    for(int i=2;i<=lim;i++){
        if(!vis[i]) {prime[++cnt]=i; f[i]=(i-(LL)i*i%MOD+MOD)%MOD;}
        for(int j=1;j<=cnt && (LL)prime[j]*i<=lim;j++){
            vis[i*prime[j]]=1;
            if(!(i%prime[j])) {
                f[i*prime[j]]=(LL)f[i]*prime[j]%MOD;
                break;
            }
            f[i*prime[j]]=(LL)f[i]*f[prime[j]]%MOD;
        }
    }
    for(int i=1;i<=lim;i++) f[i]=(f[i-1]+f[i])%MOD;
}

inline int calc(int x,int y){
    return (LL)((LL)(x+1)*x/2%MOD)*((LL)(y+1)*y/2%MOD)%MOD;
}

inline void work(){
    int ans=0;
    n=rd(),m=rd();if(n>m) swap(n,m);
    for(int l=1,r;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=(ans+(LL)calc(n/l,m/l)*(f[r]-f[l-1]+MOD)%MOD)%MOD;
    }
    printf("%d\n",ans);
}

int main(){
    prework();T=rd();
    while(T--) work();  
    return 0;
}

转载于:https://www.cnblogs.com/sdfzsyq/p/10279408.html