51 NOD 1238 最小公倍数之和 V3

原题链接

 

最近被51NOD的数论题各种刷……（NOI快到了我在干什么啊！

然后发现这题在网上找不到题解……那么既然A了就来骗一波访问量吧……

（然而并不怎么会用什么公式编辑器，写得丑也凑合着看吧……

$$
ANS=sum_{i=1}^{n}sum_{j=1}^{n} frac{i*j}{gcd(i,j)}
$$
$$
=sum_{d=1}^{n} d*sum_{i=1}^{leftlfloorfrac{n}{d} ight floor} sum_{j=1}^{leftlfloorfrac{n}{d} ight floor} i*j*[gcd(i,j)==1]
$$
$$
=sum_{d=1}^{n} d*sum_{d'=1}^{leftlfloorfrac{n}{d} ight floor} f(leftlfloorfrac{n}{d*d'} ight floor)*d'^2*μ(d')
$$
$$
=sum_{d'=1}^{n} d'^2*μ(d')*sum_{d=1}^{leftlfloorfrac{n}{d'} ight floor} d*f(leftlfloorfrac{n}{d*d'} ight floor)
$$

这里$f(x)=(frac{x*(x+1)}{2})^2$，即1到x中数字之间两两乘积之和。

 

可以看出不同的$leftlfloorfrac{n}{d'} ight floor$只有根号n个，对于某一个$leftlfloorfrac{n}{d'} ight floor$，不同的$leftlfloorfrac{n}{d*d'} ight floor$也只有根号n个，那么把它们压在一起处理就好了。

但是要做到这一点需要快速求出$d'^2*μ(d')$的前缀和，这时就要用上杜教筛了。

记$g(x)=x^2*μ(x)$，$S(x)=sum_{i=1}^{x} g(i)$

那么

$  sum_{i=1}^{n} sum_{x|i} g(x)*(frac{i}{x})^2$

$=sum_{i=1}^{n}  i*sum_{x|i} μ(x)$

$=sum_{i=1}^{n} [i==1] $

$= 1$

同时

$  sum_{i=1}^{n} sum_{x|i} g(x)*(frac{i}{x})^2$

$=sum_{i=1}^{n} sum_{x=1}^{[frac{n}{i}]} g(x)*i^2$

$=sum_{i=1}^{n} i^2*S[frac{n}{i}]$

可以得到$S(n)=1-sum_{i=2}^{n} i^2*S[frac{n}{i}]$，预处理+哈希维护一波就行了

但是不知道是我常数太大，还是方法没别人优越，卡了很久常数才卡过去。

#include<ctime>
#include<cstdio>
#include<algorithm>
#define ll long long
#define N 1000001
#define K 500000004
using namespace std;
ll read_p,read_ca,read_f;
inline ll read(){
    read_p=0;read_ca=getchar();read_f=1;
    while(read_ca<'0'||read_ca>'9') {if (read_ca=='-') read_f=-1;read_ca=getchar();}
    while(read_ca>='0'&&read_ca<='9') read_p=read_p*10+read_ca-48,read_ca=getchar();
    return read_p*read_f;
}
const int HA=1e6+7;
const int MOD=1e9+7;
int MMH=0,p[N],num=0,mu[N],_ff[N],l[HA],Num,f[N],s[N],U=0;
ll n,P;
bool bo[N];
struct na{int z,ne;ll y;}b[HA];
inline void M(int &x){while (x>=MOD) x-=MOD;while(x<0)x+=MOD;}
inline ll MM(ll x){return x<MOD&&x>-MOD?x:x%MOD;}
inline void in(int x,ll y,int z){b[++Num].y=y;b[Num].z=z;b[Num].ne=l[x];l[x]=Num;}
inline int _f(ll x){
    if (x>MOD) x%=MOD;
    if (x<N) return _ff[x];
    return MM(1LL*x*(x+1)>>1);
}
inline ll sqr(ll x){return x*x%MOD;}
inline int F(ll x){
    int MMH=0;ll P;
    for (register ll i=1,j;i<=x;i=j+1) j=x/(P=x/i),M(MMH+=MM(sqr(_f(P))*(_f(j)-_f(i-1))));
    return MMH;
}
inline int Q(ll x){if (x>MOD) x%=MOD;return MM(MM(x*(x+1))*(x+x+1))*166666668%MOD;}
inline int mmh(ll x){
    if (x<N) return s[x];
    register ll i,j;ll P;int o;
    for (i=l[o=x%HA];i;i=b[i].ne) if (b[i].y==x) return b[i].z;
    int MMH=0;
    for (i=2;i<=x;i=j+1) j=x/(P=x/i),M(MMH=MM(1LL*(Q(j)-Q(i-1))*mmh(P)+MMH));
    M(MMH=1-MMH);in(o,x,MMH);return MMH;
}
int main(){
    register int i,j;mu[1]=1;
    for (i=2;i<N;i++){
        if (!bo[i]) p[++num]=i,mu[i]=-1;
        for (j=1;j<=num&&i*p[j]<N;j++){
            bo[i*p[j]]=1;
            if (i%p[j]) mu[i*p[j]]=-mu[i];else {mu[i*p[j]]=0;break;}
        }
    }
    for (i=1;i<N;i++) M(f[i]=1LL*i*i*mu[i]%MOD+MOD),M(s[i]=s[i-1]+f[i]),_ff[i]=MM(1LL*i*(i+1)>>1);
    n=read();
    for (ll i=1,j;i<=n;i=j+1) j=n/(P=n/i),M(MMH+=MM(1LL*(mmh(j)-mmh(i-1))*F(P)));
    printf("%d
",MMH);
}