【学术篇】51nod 1238 最小公倍数之和
题目大意
求
n∑i=1n∑j=1lcm(i,j)n∑i=1n∑j=1lcm(i,j)
一看就是辣鸡反演一类的题目, 那就化式子呗..
n∑i=1n∑j=1lcm(i,j)=n∑i=1n∑j=1ijgcd(i,j)=n∑i=1n∑j=1n∑k=1ijk[gcd(i,j)=k]=n∑k=1⌊nk⌋∑i=1⌊nk⌋∑j=1ijk[gcd(i,j)=1]=n∑k=1k⌊nk⌋∑i=1i⌊nk⌋∑j=1j[gcd(i,j)=1]=n∑k=1k⌊nk⌋∑i=1i⌊nk⌋∑j=1j⌊nk⌋∑d=1[d|i][d|j]μ(d)=n∑k=1k⌊nk⌋∑d=1μ(d)⌊nk⌋∑i=1[d|i]i⌊nk⌋∑j=1[d|j]j设i=dq,j=dq,原式=n∑k=1k⌊nk⌋∑d=1μ(d)⌊nk⌋∑p=1dp⌊nk⌋∑q=1dq=n∑k=1k⌊nk⌋∑d=1μ(d)∗d2⌊nk⌋∑p=1p⌊nk⌋∑q=1q∵n∑i=1i=n(n+1)2,∴原式=n∑k=1k⌊nk⌋∑d=1μ(d)∗d2(⌊nk⌋(⌊nk⌋+1)2)2设t=kd,原式=n∑t=1(⌊nt⌋(⌊nt⌋+1)2)2∑d|tμ(d)∗td设g(x)=∑d|xμ(d)∗xd,则原式=n∑t=1(⌊nt⌋(⌊nt⌋+1)2)2g(t)
然后前面这一堆我们可以分块O(1)搞出来, 所以问题的关键就变成了求g(x)的前缀和.
因为题目中n≤1010, 所以这里要用O(n23)的杜教筛.
但是杜教筛的框架我们是知道的,
设要求前缀和的函数为f(x),其前缀和为sf(x)设有好求前缀和的积性函数g(x)和h(x),使得(f∗g)(x)=h(x),设h(x)的前缀和为sh(x),则有sf(x)=sh(x)−∑nd=2sf(⌊nd⌋)g(d)g(1)
所以我们的任务又变成了找到合适的g(x)和h(x).
在某篇歪果仁的blog上, 他说:
With a little luck, we can notice that Id(l) = g * Id2(l),
- 注: id(x)=x(以下称为n)id2(x)=x2(以下称为n2)
But 这运气也太好了点吧? 我怎么就找不到这种函数呢?!
不过我们还是要考虑一下的... (这波化式子的时候注意乘号⋅和卷积号∗)
g(x)=∑d|xμ(d)xd=∑d|xd2μ(d)xd=(n2⋅μ∗n)(x)(g∗n2)(x)=((n2⋅μ∗n)∗n2)(x)⋯①积性函数有一个性质,如果两边的乘积都含有某一项的时候,可以把这一项提出来.所以我们可以把x2都提出来,(证明可以用定义式推咯 )①=((n2⋅μ∗n2)∗n)(x)=(n2⋅(μ∗1)∗n)(x)=(n2⋅ϵ∗n)(x)=n(x)
所以我们就证明出来了(g∗n2)(x)=n(x)
这样就可以直接杜教筛了= =
然后就是注意细节了OvO
比如非常坑人的1010∗(1e9+7)会爆long long!!
所以要开unsigned...
就这样吧= =
代码:
#include <map>
#include <cstring>
#include <iostream>
using namespace std;
typedef unsigned long long LL;
#define ri register int
#define rll register unsigned long long
const int p=1e9+7;
int qpow(int a,int b,int s=1){
for(;b;b>>=1,a=1LL*a*a%p)
if(b&1) s=1LL*s*a%p;
return s;
}
const int inv2=qpow(2,p-2),inv4=qpow(4,p-2),inv6=qpow(6,p-2);
const int N=4800000,_N=N+10;
const int SQ=233333;
LL n;
map<LL,int> mp;
map<LL,int>::iterator it;
int g[_N],pr[N/2],tot;
bool notp[_N];
void shai(){
g[1]=1; ri i,j; rll k;
for(i=2;i<=N;++i){
if(!notp[i])
pr[++tot]=i,g[i]=(-1LL*i*(i-1)%p+p)%p;
for(j=1;j<=tot&&(k=1LL*i*pr[j])<=N;++j){
notp[k]=1;
if(i%pr[j]==0){g[k]=1LL*g[i]*pr[j]%p;break;}
g[k]=1LL*g[i]*g[pr[j]]%p;
}
}
for(i=2;i<=N;++i) g[i]=(g[i]+g[i-1])%p;
}
inline int calcsqr(LL x){x%=p;return x*(x+1)%p*(2*x+1)%p*inv6%p;}
inline int calcsum(LL x){x%=p;return x*x%p*(x+1)%p*(x+1)%p*inv4%p;}
int calc(LL x){
if(x<=N) return g[x];
it=mp.find(x);
if(it!=mp.end()) return it->second;
int ans=1LL*x%p*(x+1)%p*inv2%p;
LL last=0;
for(rll i=2;i<=x;i=last+1){
last=x/(x/i);
ans=(ans-1LL*calc(x/i)*(calcsqr(last)-calcsqr(i-1)+p)%p+p)%p;
}
return mp[x]=(ans%p+p)%p;
}
int main(){
shai();cin>>n;
LL last=0; int ans=0;
for(rll i=1;i<=n;i=last+1){
last=n/(n/i);
ans=(ans+1LL*calcsum(n/i)*(calc(last)-calc(i-1)+p)%p)%p;
}
cout<<ans;
}
这破题我连推带写加总结写了整整一天!!!
窝还是太弱了..
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
· 没有源码,如何修改代码逻辑?
· 一个奇形怪状的面试题:Bean中的CHM要不要加volatile?
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· PowerShell开发游戏 · 打蜜蜂
· 在鹅厂做java开发是什么体验
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战