[JZOJ5134][SDOI省队集训2017]三元组

题目大意

求$\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c[(i,j)=1][(i,k)=1][(j,k)=1]$

推式子

首先假设a<=b<=c。
第一步转化为
$\sum_{i=1}^a\sum_{j=1,(j,i)=1}^b\sum_{k=1,(k,i)=1}^c[(j,k)=1]$
注意到$[n=1]=\sum_{d|n}\mu(d)$
所以转化为
$\sum_{i=1}^a\sum_{j=1,(j,i)=1}^b\sum_{k=1,(k,i)=1}^c\sum_{d|j,d|k}\mu(d)$
交换一下枚举顺序
$\sum_{i=1}^a\sum_{d=1,(d,i)=1}^b\mu(d)*\sum_{j=1,(j,i)=1,d|j}^b\sum_{k=1,(k,i)=1,d|k}^c1$
后面可以直接计算。
$\sum_{i=1}^a\sum_{d=1,(d,i)=1}^b\mu(d)*(\sum_{j=1,(j,i)=1}^{\frac{b}{d}}1)*(\sum_{k=1,(k,i)=1}^{\frac{c}{d}}1)$
我们可以直接枚举i，然后考虑计算。
这个形式其实和一道叫循环之美的题差不多哦！
首先观察后两维可以分块。
于是我们需要计算一个和k互质的<=n的i的莫比乌斯函数和或者个数。
计算这个即可。
容易想到洲阁筛。
设S(i,r)表示<=r的与k的前i个质因子互质的莫比乌斯函数和。
$S(i,r)=S(i-1,r)-\mu(pi)*S(i,r/pi)$
注意因为莫比乌斯函数要不含平方因子，所以后面那个第一维仍然是i，即要求它只有一个pi。
设T(i,r)表示<=r的与与k的前i个质因子互质的数的个数。
$T(i,r)=T(i-1,r)-T(i-1,r/pi)$
这样递归加记忆化有点慢。
我们考虑处理g[i]表示将i每个质因子的次数改成1的数。
一个优化是一开始枚举i的时候相同g的可以利用之前的答案。
于是我们只需要计算莫比乌斯函数不为0的i。
然后设f[i]表示i的最大质因数。
我们把S的意义改为<=r的与i互质的莫比乌斯函数和，T同理。
这样的好处是可以全局记忆化。
注意到状态只有n^1.5，可以开数组存。
为什么n^1.5？因为只有分块边界会进去，而不断除以一个质数的过程可以看做常数级别。
递归会比较慢，可以改为递推的形式提前计算。
详见代码。

#include<cstdio>
#include<algorithm>
#include<cmath>
#include<map>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const int maxn=50000+10,maxnn=700+10,mo=1000000007;
int pri[maxn],mu[maxn],sum[maxn],f[maxn],g[maxn],num[maxn],p[maxn];
int S[maxnn][maxn],T[maxnn][maxn],id[maxn],di[maxn];
bool bz[maxn],pd[maxn],dp[maxn];
int i,j,k,l,r,t,n,m,tot,pot,top,cnt,ans,a,b,c;
void work(int k){
    int i,j,t,l,r,sum=0;
    k=g[k];
    if (num[k]){
        (ans+=num[k])%=mo;
        //if (ans>=mo) ans-=mo;
        return;
    }
    i=1;t=0;
    while (i<=b){
        j=min(b/(b/i),c/(c/i));
        l=S[id[j]][k];
        r=(ll)T[di[b/i]][k]*T[di[c/i]][k]%mo;
        (sum+=(ll)(l-t)*r%mo)%=mo;
        //if (sum>=mo) sum-=mo;
        t=l;
        i=j+1;
    }
    (ans+=sum)%=mo;
    //if (ans>=mo) ans-=mo;
    num[k]=sum;
}
int main(){
    freopen("triple.in","r",stdin);freopen("triple.out","w",stdout);
    mu[1]=g[1]=1;
    fo(i,2,maxn-10){
        if (!bz[i]){
            pri[++top]=i;
            f[i]=g[i]=i;
            mu[i]=-1;
        }
        fo(j,1,top){
            if ((ll)i*pri[j]>maxn-10) break;
            bz[i*pri[j]]=1;
            f[i*pri[j]]=f[i];
            if (i%pri[j]==0){
                g[i*pri[j]]=g[i];
                break;
            }
            g[i*pri[j]]=g[i]*pri[j];
            mu[i*pri[j]]=-mu[i];
        }
    }
    fo(i,1,maxn-10) sum[i]=sum[i-1]+mu[i];
    scanf("%d%d%d",&a,&b,&c);
    if (a>b) swap(a,b);
    if (b>c) swap(b,c);
    if (a>b) swap(a,b);
    i=1;
    while (i<=b){
        j=min(b/(b/i),c/(c/i));
        pd[j]=1;
        dp[b/i]=1;dp[c/i]=1;
        i=j+1;
    }
    fo(i,1,c)
        if (i==g[i]) p[++cnt]=i;
    fd(i,c,1){
        if (pd[i]){
            fo(j,1,top){
                if (i<pri[j]) break;
                pd[i/pri[j]]=1;
            }
        }
        if (dp[i]){
            fo(j,1,top){
                if (i<pri[j]) break;
                pd[i/pri[j]]=1;
            }
        }
    }
    fo(i,1,c)
        if (pd[i]){
            id[i]=++tot;
            if (i==1){
                fo(j,1,cnt) S[tot][p[j]]=1;
                continue;
            }
            S[tot][1]=sum[i];
            fo(j,2,cnt){
                k=p[j];
                S[tot][k]=(S[tot][k/f[k]]-mu[f[k]]*S[id[i/f[k]]][k])%mo;
            }
        }
    fo(i,1,c)
        if (dp[i]){
            di[i]=++pot;
            if (i==1){
                fo(j,1,cnt) T[pot][p[j]]=1;
                continue;
            }
            T[pot][1]=i;
            fo(j,2,cnt){
                k=p[j];
                T[pot][k]=(T[pot][k/f[k]]-T[di[i/f[k]]][k/f[k]])%mo;
            }
        }
    fo(i,1,a) work(i);
    (ans+=mo)%=mo;
    printf("%d\n",ans);
}