51nod1238 最小公倍数之和 V3 莫比乌斯函数 杜教筛
题意:求∑ni=1∑nj=1lcm(i,j)∑ni=1∑nj=1lcm(i,j).
题解:虽然网上很多题解说用mu卡不过去,,,不过试了一下貌似时间还挺充足的,。。也许有时间用phi试试?
因为是用的莫比乌斯函数求的,所以推导比大部分题解多。。。而且我写式子一般都比较详细,所以可能看上去很多式子,实际上是因为每一步都写了,几乎没有跳过的。所以应该都可以看懂的。
末尾的e函数是指的e[1]=1,e[x]=0(x!=1)这样一个函数
n∑i=1n∑j=1lcm(i,j)
n∑i=1n∑i=1ijgcd(i,j)
枚举gcd
n∑d=1⌊nd⌋∑i=1⌊nd⌋∑j=1[gcd(i,j)==d]ijd
因为(ijd2d=ijd),所以:
n∑d=1⌊nd⌋∑i=1⌊nd⌋∑j=1[gcd(i,j)==d]ijd
n∑d=1d⌊nd⌋∑i=1⌊nd⌋∑j=1ij[gcd(i,j)==1]
n∑d=1d⌊nd⌋∑i=1⌊nd⌋∑j=1ij∑k|gcd(i,j)μ(k)
枚举k,再枚举k的倍数。
n∑d=1d⌊nd⌋∑k=1μ(k)⌊ndk⌋∑i=1ik⌊ndk⌋∑j=1jk
设S(n)=∑ni=1i
n∑d=1d⌊nd⌋∑k=1μ(k)k2S(ndk)
枚举T=dk
n∑T=1S(nT)2∑k|Tμ(k)k2Tk
n∑T=1S(nT)2∑k|Tμ(k)kT
n∑T=1S(nk)2⋅T∑k|Tμ(k)k
设f(T)=T∑k|Tμ(k)k,卷上id2,因为S(nk)可以数论分块,所以我们只需要快速求出区间[l,r]内的f之和即可,显然求出f的前缀和即可解决问题
(f∗id2)(n)=∑i|nf(i)n2i2=∑i|ni∑k|iμ(k)kn2i2
∑i|n∑k|iμ(k)kn2i=n∑i|n∑k|iμ(k)kni
设h(i)=∑k|iμ(k)k,则原式:
n∑i|nh(i)ni=n(h∗id)(n)
(f∗id2)(n)=n(h∗id)(n)
h(n)=∑k|nμ(k)k=(μ⋅id)∗1
f∗id2=n[(μ⋅id)∗1∗id]=n[(μ⋅id)∗id∗1]
其中(μ⋅id)∗id=∑i|nμ(i)ini=n∑i|nμ(i)=e
所以
n[(μ⋅id)∗id∗1]=n[e∗1]=n
带入杜教筛的式子:
g(1)S(n)=n∑i=1(f∗g)(i)−n∑i=2g(i)S(ni)
=n∑i=1i−n∑i=2i2S(ni)
然后直接上杜教就可以了.
其实还有一个问题。。。一开始预处理的前缀和怎么求?
要知道前缀和,首先要求出f.
因为f(T)=T∑k|Tμ(k)k,所以如果我们可以快速求出∑k|Tμ(k)k,然后就只需要再O(n)的乘上T就可以了.
我们先预处理出μ(k),然后对于每一个k,枚举它的倍数,统计贡献。那么复杂度为 n1+n2+...+nn=nlogn(此处的n为原题面的23次方,即要预处理的f个数)
#include<bits/stdc++.h>
using namespace std;
#define R register int
#define LL long long
#define RL register LL
#define AC 3000
#define ac 5000000
#define p 1000000007LL
//#define h(x) ((x <= block) ? sum[x] : S[n / x])
LL n, ans, block;
LL mu[ac], S[AC], sum[ac], inv[AC];
int pri[ac], tot;
bool z[ac], vis[AC];
inline LL h(LL x)
{
return ((x <= block) ? sum[x] : S[n / x]);
}
inline void up(LL & a, LL b)
{
a += b;
if(a >= p) a -= p;
if(a <= -p) a += p;
}
LL count(LL l, LL r){
return (r - l + 1) % p * ((r + l) % p) % p * inv[2] % p;
}
void pre()
{
scanf("%lld", &n), block = pow(n, 0.66666);
mu[1] = 1;
for(R i = 2; i <= block; i ++)
{
if(!z[i]) pri[++ tot] = i, mu[i] = -1;
for(R j = 1; j <= tot; j ++)
{
int now = pri[j];
if(i * now > block) break;
z[i * now] = true;
if(!(i % now)) break;
mu[i * now] = - mu[i];
}
}
inv[1] = 1;
for(R i = 2; i <= 10; i ++) inv[i] = (p - p / i) * inv[p % i] % p;
for(R i = 1; i <= block; i ++)//枚举mu(i)
for(R j = 1; j; j ++)//枚举i的倍数
{
if(j * i > block) break;
up(sum[i * j], mu[i] * i % p);
}
for(R i = 1; i <= block; i ++) sum[i] = sum[i] * i % p;
for(R i = 1; i <= block; i ++) up(sum[i], sum[i - 1]);//算出f数组后还要统计前缀和
}
LL get(LL x)
{
x %= p;
return x * (x + 1) % p * (2 * x + 1) % p * inv[6] % p;
}
void cal(LL x)
{
if(x <= block || vis[n / x]) return ;
LL rnt = count(1, x);
for(RL i = 2, lim, now; i <= x; i = lim + 1)
{
lim = x / (x / i), now = x / i, cal(now);
up(rnt, - ((get(lim) - get(i - 1)) % p * h(now) % p));
}
S[n / x] = rnt, vis[n / x] = true;
}
void work()
{
for(RL i = 1, lim, now, x; i <= n; i = lim + 1)
{
lim = n / (n / i), now = n / i, x = count(1, now);
up(ans, (x % p * x % p * ((h(lim) - h(i - 1)) % p) % p));
}
printf("%lld\n", (ans + p) % p);
}
int main()
{
//freopen("in.in", "r", stdin);
pre();
cal(n);
work();
// fclose(stdin);
return 0;
}
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· AI与.NET技术实操系列:基于图像分类模型对图像进行分类
· go语言实现终端里的倒计时
· 如何编写易于单元测试的代码
· 10年+ .NET Coder 心语,封装的思维:从隐藏、稳定开始理解其本质意义
· .NET Core 中如何实现缓存的预热?
· 25岁的心里话
· 闲置电脑爆改个人服务器(超详细) #公网映射 #Vmware虚拟网络编辑器
· 基于 Docker 搭建 FRP 内网穿透开源项目(很简单哒)
· 零经验选手,Compose 一天开发一款小游戏!
· 一起来玩mcp_server_sqlite,让AI帮你做增删改查!!