来自FallDream的博客,未经允许,请勿转载,谢谢。
------------------------------------------------------------------------------------------
题意:求$$sum_{i=1}^{n}sum_{j=1}^{n}lcm(i,j) \ nleqslant 10^{10}$$
题解:题目即求$$sum_{i=1}^{n}sum_{j=1}^{n}frac{i*j}{gcd(i,j)}$$
$$=sum_{d=1}^{n}d*sum_{i=1}^{lfloor n/d floor}sum_{j=1}^{lfloor n/d floor}i*j*[gcd(i,j)=1]$$
已知$$sum_{i=1}^{n}i*[gcd(i,n)=1]=frac{n*varphi(n)}{2}$$
所以所求即为$$sum_{d=1}^{n}d*sum_{i=1}^{lfloor n/d floor}i*i*varphi(i)$$
$lfloorfrac{n}{d} floor$只有$sqrt(n)$种取值,那么我们考虑快速求出$g(i)=i^{2}*varphi(i)$的前缀和$S(i)$。
$$sum_{n|d}varphi(d)=n$$
$$sum_{i=1}^{n}sum_{d|n}varphi(d)=frac{n(n+1)}{2}$$
$$sum_{i=1}^{n}sum_{d|n}varphi(d)*i^{2}=frac{n^{2}*(n+1)^{2}}{4}$$
$$sum_{i=1}^{n}sum_{d|n}g(d)*(frac{i}{d})^{2}=frac{n^{2}*(n+1)^{2}}{4}$$
$$sum_{i=1}^{n}sum_{d=1}^{lfloor n/i floor}g(d)*i^{2}=frac{n^{2}*(n+1)^{2}}{4}$$
$$S(i)=frac{n^{2}*(n+1)^{2}}{4}-sum_{i=2}^{n}i^{2}*S(lfloor n/i floor)$$
这个可以在$Oleft(n^{frac{2}{3}} ight)$时间内做完。此题得解。
-------------------------------------------------------------------
我好菜啊,推了好久.....
-----
#include<iostream> #include<cstdio> #include<cstring> #include<cmath> #include<map> #define MAXN 5000000 #define ll long long #define mod 1000000007 #define ditoly 6666666 #define inv2 500000004 #define inv4 250000002 #define inv6 166666668 using namespace std; inline ll read() { ll x=0,f=1;char ch=getchar(); while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=getchar();} while(ch>='0'&&ch<='9'){x=x*10+ch-'0'; ch=getchar();} return x*f; } struct mymap{ ll x,ans;int next; }e[8000000]; ll phi[MAXN+5],n,ans=0; int s[MAXN/5],num=0,head[ditoly+5]; bool b[MAXN+5]; inline ll getcube(ll x){x%=mod;return x*(x+1)%mod*(x<<1|1)%mod*inv6%mod;} inline void ins(ll x,ll sum) { int j=x%ditoly; e[++num]=(mymap){x,sum,head[j]};head[j]=num; } inline ll getsq(ll x){x%=mod;x=x*(x+1)%mod;return x*x%mod*inv4%mod;} ll calc(ll x) { if(x<=MAXN)return phi[x]; for(int i=head[x%ditoly];i;i=e[i].next) if(e[i].x==x)return e[i].ans; ll last,sum=getsq(x); for(ll i=2;i<=x;i=last+1) { last=x/(x/i); sum-=(getcube(last)-getcube(i-1)+mod)%mod*calc(x/i)%mod; while(sum<0)sum+=mod; } ins(x,sum); return sum; } int main() { n=read();phi[1]=1; for(int i=2;i<=MAXN;i++) { if(!b[i]) phi[i]=i-1,s[++num]=i; for(int j=1;s[j]*i<=MAXN;j++) { b[s[j]*i]=1; if(i%s[j]==0){phi[s[j]*i]=phi[i]*s[j];break;} phi[s[j]*i]=phi[i]*(s[j]-1); } }num=0; for(int i=2;i<=MAXN;i++) phi[i]=(phi[i-1]+1LL*i*i%mod*phi[i]%mod)%mod; for(ll i=1,last;i<=n;i=last+1) { last=n/(n/i);ll x=(n/i)%mod; ans+=x*(x+1)%mod*inv2%mod*((calc(last)-calc(i-1)+mod)%mod)%mod; while(ans>=mod)ans-=mod; } cout<<(ans+mod)%mod; return 0; }