[BestCoder Round #68][HDU5608]function

题目大意

已知函数$f(x)$满足：

$$n^2-3n+2=\sum_{d|n}f(d)$$
给定 $n$，请计算$\sum_{i=1}^nf(i)$对 $10^9+7$的结果。
一个测试点有 $T$组数据。

$T\le500,n\le10^9$，只有五组数据的 $n>10^6$

---

题目分析

设$G(x)=x^2-3x+2$，由题意得$f$函数和$\mathbb1$函数的狄利克雷卷积是$G$函数。
考虑使用杜教筛套路，设$S(n)=\sum_{i=1}^nf(i)$：

$$S(n)=\sum_{i=1}^nG(i)-\sum_{i=2}^nS\left(\lfloor\frac n i\rfloor\right)$$
其中 $G(i)$的前缀和使用数列求和公式直接得到，后面的分块。
直接计算的时间复杂度是 $\mathrm O(n^\frac34)$的，会T掉。
考虑算出小于等于 $X$的部分，对后面的使用杜教筛。怎么计算前面的部分呢？
对$G(n)=\sum_{d|n}f(d)$使用莫比乌斯反演，得到 $f(n)=\sum_{d|n}\mu(\frac nd)G(d)$。这个显然可以像埃拉托色尼筛法一样，枚举 $\mu(i)$然后计算其对所有 $i$的倍数的$f$的贡献。时间复杂度是 $\mathrm O(X\log X)$的。
$X$取$n^\frac23$的话时间复杂度是 $\mathrm O(n^\frac23\log n^{\frac23}+Tn^\frac23)$的。当然你也可以玄学调整 $X$的大小，单组询问的话将$\log$视作 $\frac14$次方，求导得到 $X$最优是$n^\frac47$。

---

代码实现

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>

using namespace std;

const int P=1000000007;
const int L=1000000;
const int MAXN=31622;

int pri[L+50],mu[L+50],f[L+50],sum[L+50];
bool mark[L+50];
bool vis[MAXN];
int S[MAXN];
int itwo,isix,T,N,s,l;

int quick_power(int x,int y)
{
    int ret=1;
    for (;y;y>>=1,x=1ll*x*x%P) if (y&1) ret=1ll*ret*x%P;
    return ret;
}

void pre()
{
    mark[1]=1,mu[1]=1,f[1]=1;
    for (int i=1;i<=L;++i)
    {
        if (!mark[i]) f[pri[++pri[0]]=i]=i,mu[i]=-1;
        for (int j=1,k;j<=pri[0];++j)
        {
            if (1ll*pri[j]*i>L) break;
            mark[k=pri[j]*i]=1,f[k]=pri[j];
            mu[k]=f[k]==f[i]?0:-mu[i];
            if (!(i%pri[j])) break;
        }
    }
    for (int i=1;i<=L;++i)
    {
        for (int j=1;i*j<=L;++j)
            (((sum[i*j]+=mu[j]*((1ll*i*i-3*i+2)%P))%=P)+=P)%=P;
        (sum[i]+=sum[i-1])%=P;
    }
}

int id(int x){return N/x;}

int sieve(int n)
{
    if (n<=L) return sum[n];
    int idn=id(n);
    if (vis[idn]) return S[idn];
    int res=((1ll*n*(n+1)%P*((n*2+1)%P)%P*isix%P-1ll*n*(n+1)%P*3%P*itwo%P)%P+n*2%P)%P;
    for (int st=2,en,x;st<=n;st=en)
    {
        x=n/st,en=n/x+1;
        res=(res-1ll*sieve(x)*(en-st)%P+P)%P;
    }
    return vis[idn]=1,S[idn]=res;
}

int main()
{
    itwo=quick_power(2,P-2),isix=quick_power(6,P-2),pre();
    freopen("function.in","r",stdin),freopen("function.out","w",stdout);
    for (scanf("%d",&T);T--;scanf("%d",&N),memset(vis,0,sizeof vis),s=trunc(sqrt(N)),printf("%d\n",sieve(N)));
    fclose(stdin),fclose(stdout);
    return 0;
}