最近被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); }