“21 天好习惯”第一期－12

P1829 [国家集训队]Crash的数字表格 / JZPTAB

求值：

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

易得原式如下：

${\sum }_{i=1}^n{\sum }_{j=1}^m\frac{i\ast j}{gcd(i,j)}$

枚举最大公因数：

${\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{d|i,d|j,gcd(\frac{i}{d},\frac{i}{d})=1}\frac{i\ast j}{d}$

非常经典的$gcd$式子的化法：

${\sum }_{d=1}d\ast {\sum }_{i=1}^{\frac{n}{d}}{\sum }_{j=1}^{\frac{m}{d}}[gcd(i,j)=1]i\ast j$

式子的后半段出现了互质数对之积之和，考率先单独拿出来,记

$sum(n,m)={\sum }_{i=1}^n{\sum }_{j=1}^m[gcd(i,j)=1]i\ast j$

有：

$sum(n,m)={\sum }_{i=1}^n{\sum }_{j=1}^m\epsilon (gcd(i,j))\ast i\ast j$

$={\sum }_{i=1}^n{\sum }_{j=1}^m{\sum }_{d|gcd(i,j)}\mu (d)\ast i\ast j$

$={\sum }_{d=1}\mu (d)\ast d^2{\sum }_{i=1}^{\frac{n}{d}}{\sum }_{j=1}^{\frac{m}{d}}i\ast j$

观察上式，前半段可以预处理前缀和；后半段又是一个范围内数对之和，记

$g(n,m)={\sum }_{i=1}^n{\sum }_{j=1}^{m}i\ast j=\frac{n(n+1)}{2}\ast \frac{m\ast (m+1)}{2}$

可以$o(1)$求解。

至此：

$sum(n,m)={\sum }_{d=1}\mu (d)\ast d^2\ast g(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )$

我们可以数论分块求解这部分。

回到定义$sum(n,m)$的地方，则原式为：

${\sum }_{d=1}d\ast sum(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )$

这又是一个可以数论分块求解的式子。

Code:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=1e7+7,mod=20101009;
int mu[maxn];
ll f[maxn];
int g(ll x,ll y) {
	return (x*(x+1)/2%mod)*(y*(y+1)/2%mod)%mod;
}
int sum(int n,int m) {
	int res(0);
	for(int d=1,x,y,r;d<=n;++d) {
		x=n/d,y=m/d;
    	r=min(n/x,m/y);
    	res=(res+(f[r]-f[d-1])*g(x,y))%mod;
		d=r;
	}
	return res;
}
bool vis[maxn];
vector<int>prime;
inline void init(int n) {
    f[1]=mu[1]=1;
    for(int i=2;i<=n;++i) {
        if(!vis[i]) {
            prime.emplace_back(i); mu[i]=-1;
        }
        for(auto &j:prime) {
            if(i*j>n) break;
            vis[i*j]=1;
            if(i%j==0) break;
            mu[i*j]=-mu[i];
        }
        f[i]=(f[i-1]+(ll)i*i*mu[i])%mod;
    }
}
signed main() {
	cin.sync_with_stdio(false), cin.tie(nullptr),cout.tie(nullptr);
	int n,m;
    cin>>n>>m;
    if(n>m) swap(n,m);//便于处理
	init(n);
    int ans(0);
    for(int d=1,x,y,r;d<=n;++d) {
    	x=n/d,y=m/d;
    	r=min(n/x,m/y);
    	ans=(ans+(ll)(r+d)*(r-d+1)/2*sum(x,y))%mod;
    	d=r;
	}
	cout<<(ans+mod)%mod<<'\n';
	return 0;
}

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-v4QRI83N-1636734163952)(https://uploadfiles.nowcoder.com/images/20210905/137609403_1630833968576/8748513304B8F3238731469EDE09F8A1 “图片标题”)]

这题线性筛比普通筛快很多。