【GDSOI 2016】第一题 互补约数

Description

简化题意：求

∑i=1n∑j|igcd(j,ij)

Input

输入包含一行，一个正整数 n。

Output

输出只有一行， F(n)。

Sample Input

Sample Input1：
10

Sample Input2：
1000

Sample Output

Sample Output1：
32

Sample Output2：
12776

Data Constraint

n<=1011

Solution

两个sigma外加gcd，直接考虑反演

∑i=1n∑j|igcd(j,ij)

这里是枚举i和i的约数，我们转换为枚举i和i的倍数

∑i=1n∑j=1⌊ni⌋gcd(i,j)

套路反演

设

f[d]=∑i=1n∑j=1⌊ni⌋gcd(i,j)=d

设

g[d]=∑i=1⌊nd⌋f[id]

反演得

f[d]=∑i=1⌊nd⌋g[id]∗μ(i)

ans=∑d=1nd∑i=1⌊nd⌋g[id]∗μ(i)

---

设 T=id，我们发现gcd(i,j)的值不可能大于n√，所以

ans=∑T=1n√g[T]∗∑i|TTi∗μ(i)

将∑n√T=1∑i|TTi∗μ(i)预处理，设为a(i)
为了方便预处理，设i∗j=T变形为

a(i)=∑i=1n√∑i=1⌊n√d⌋j∗μ(i)

我们知道，如果gcd(x,y)为T的倍数，那么x,y也是T的倍数
那么x=x′∗T,y=y′∗T
因为x∗y<=n
所以x′
枚举x′∗T∗y′∗T<=n
所以

g[T]=∑x=1⌊nT2⌋⌊nT2x⌋

---

于是，你会惊奇的发现，答案变成了

ans=∑T=1n√g[T]∗a(t)

很简单对吧
但还是过不了
求g[t]的时候用分块优化一下，而t本身不能分块（为什么？）

时间复杂度O(n34)

Code

#include<cstdio>#include<algorithm>#include<cstring>#include<cmath>#define ll long long#define fo(i,a,b) for(int i=a;i<=b;i++)#define N 401000using namespace std;ll n,mu[N],a[N],bz[N],s[N];int main(){    freopen("gcd.in","r",stdin);freopen("gcd.out","w",stdout);    scanf("%lld",&n);    ll m=sqrt(n);    mu[1]=1;    fo(i,2,m)    {        if(!bz[i]) s[++s[0]]=i,mu[i]=-1;        fo(j,1,s[0])        {            int k=i*s[j];            if(k>m) break;            bz[k]=1;            if(i%s[j]==0) {mu[k]=0;break;}            mu[k]=-mu[i];        }    }    fo(i,1,m) fo(j,1,m/i) a[i*j]+=j*mu[i];    ll ans=0;    for(ll t=1;t<=m;t++)    {        ll jy=0,q=n/(t*t);        for(ll x=1,y=1;x<=q;y=x+1,x=min((q/(q/y)),q))         {            jy+=(x-y+1)*(q/x);            if(x>=q) break;        }        ans+=(jy*a[t]);    }    printf("%lld",ans);}