51nod 1220 约数之和

定义除数函数： $$\sigma_{k}(n)=\sum_{a|n}a^k$$

令$d=\sigma_1$

题目要求计算：

$$\sum_{i=1}^n\sum_{j=1}^nd(ij)$$

之前做BZOJ的时候做过一题。

http://blog.csdn.net/zlh_hhhh/article/details/77849859

BZOJ 4176 Lucas的数论

对于 $$\sigma_0(nm)=\sum_{a|nm}1=\sum_{a|n}\sum_{b|m}[a\perp b]$$

$\sigma_1$也有类似结论：

$$\sigma_1(nm)=d(nm)=\sum_{a|nm}a\ \ =\sum_{a|n}\sum_{b|m}\frac{bn}{a}[a\perp b]$$

证明类似于Lucas的数论：

$$d(nm)=\sum_{a|nm}a\\n=\prod_{k=1}^rP_k^{x_k}\\m=\prod_{k=1}^{r}P_{k}^{y_k}$$

$$d(nm)=\sum_{a|nm}a\\=\sum_{k=1}^r\sum_{c\big|\frac{nm}{P_k^{x_k+y_k}}}\sum_{a|P_k^{x_k+y_k}}ac\\=\sum_{a_1\big|P_1^{x_1+y_1}}\sum_{a_2\big|P_2^{x_2+y_2}}...\sum_{a_r\big|P_r^{x_r+y_r}}\Big(\prod_{i=1}^ra_i\Big)$$

对于某个素数的$P_k$ $$\sum_{a\big | P_k^{x_k+y_k}}a=\sum_{a|P_{k}^{x_k}}\sum_{b|P_{k}^{y_k}}[a\perp b]\frac{P_{k}^{x_k}b}{a}$$

所以：

$$d(nm)=\sum_{a|n}\sum_{b|m}\frac{nb}{a}[a \perp b]$$

$$answer=\sum_{i=1}^n\sum_{j=1}^n\sum_{a|i}\sum_{b|j}\frac{ib}{a}[a \perp b]$$

记： $$f(k)=\sum_{i=1}^n\sum_{j=1}^n\sum_{a|i}\sum_{b|j}\frac{ib}{a}[gcd(a ,b)=k]\\F(k)=\sum_{k|a}f(a)$$

$$F(k)=\sum_{i\leq n,k|i}^n\sum_{j\leq n,k|j}^n\sum_{a|\frac{i}{k}}\sum_{b|\frac{j}{k}}\frac{ib}{a}\\=k\sum_{i\leq n,k|i}^n\sum_{j\leq n,k|j}^n\sum_{a|\frac{i}{k}}\sum_{b|\frac{j}{k}}\frac{\frac{i}{k}b}{a}\\=k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{a|i}\sum_{b|j}\frac{ib}{a}\\=k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}\Big(\sum_{a|i}\frac{i}{a}\Big)\Big(\sum_{b|j}b\Big)\\=k\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}d(i)d(j)=kS_d\big(\big\lfloor\frac{n}{k}\big\rfloor\big)^2$$

对于某一算术函数$g$ , $S_g$为：

$$S_g(n)=\sum_{i=1}^ng(i)$$

所以：

$$answer=f(1)=\sum_{k=1}^n\mu(k)kS_d\big(\big\lfloor\frac{n}{k}\big\rfloor\big)^2$$

令：$h(n)=\mu(n)n$

我们知道：$\varphi^-=h*I$

注：不太明白的，可以查看这个：http://blog.csdn.net/zlh_hhhh/article/details/77857141

所以：$h*I*\varphi=h*N=e$

对于$S_d$.我么即可可以用杜教筛来算。也可以直接$O(\sqrt{n}$分段算。

杜教筛的话：

$$d(n)=\sum_{a|n}a$$

可以看出：$d=N*I$

所以：$d*\mu=N*I*\mu=N$

$O(\sqrt{n})$算的话：

$$\sum_{i=1}^nd(i)=\sum_{i}^n\sum_{a|i}a=\sum_{a=1}^n\sum_{a|i,i\leq n}a=\sum_{a=1}^na\big\lfloor\frac{n}{a}\big\rfloor\\=\sum_{L=1}^mL\Big(S_N\big(\big\lfloor\frac{n}{L}\big\rfloor\big)-S_N\big(\big\lfloor\frac{n}{L+1}\big\rfloor\big)\Big)-\sum_{i=1}^{\lfloor\frac{n}{m+1}\rfloor}i\big\lfloor\frac{n}{i}\big\rfloor$$

代码：

#include <algorithm>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define MAXN 1000006
using namespace std;
typedef long long LL;
const int P=1e9+7;

int mu[MAXN];
int phi_[MAXN];
int Sp[MAXN];
int Sd[MAXN];
int M[MAXN];//梅藤斯函数
int _Sp[MAXN];
int _Sd[MAXN];
int _M[MAXN];

void init()
{
    mu[1]=-1;
    for(int i=1;i<MAXN;i++)
    {
        mu[i]=-mu[i];
        for(int j=i<<1;j<MAXN;j+=i) mu[j]+=mu[i];
    }
    for(int i=1;i<MAXN;i++)
    {
        M[i]=M[i-1]+mu[i];
        phi_[i]=i*mu[i];
        Sp[i]=phi_[i]+Sp[i-1];
        for(int j=i;j<MAXN;j+=i) Sd[j]+=i;
        Sd[i]+=Sd[i-1];
        if(Sd[i]>=P)Sd[i]-=P;
        if(Sp[i]<0)Sp[i]+=P;
        if(Sp[i]>=P)Sp[i]-=P;
        if(M[i]>=P)M[i]-=P;
        if(M[i]<0)M[i]+=P;
    }
}

int Sn(int n)
{
    if(n&1)
        return (LL)((n+1)>>1)*n%P;
    else
        return (LL)(n>>1)*(n+1)%P;
}

void clat_d(int n,int d)
{
    int ans=0,m=sqrt(n)+1;
    for(int L=1;L<m;L++)
    {
        int f1=n/L,f2=n/(L+1);
        if(f1<MAXN)
            f1=M[f1];
        else
            f1=_M[L*d];
        if(f2<MAXN)
            f2=M[f2];
        else
            f2=_M[(L+1)*d];
        ans+=(LL)Sd[L]*((f1-f2+P)%P)%P;
        if(ans>=P)ans-=P;
    }
    for(int i=n/m;i>1;i--)
    {
        int u=n/i;
        if(u<MAXN)
            u=Sd[u];
        else
            u=_Sd[d*i];
        ans+=mu[i]*u;
        if(ans<0)ans+=P;
        if(ans>=P)ans-=P;
    }
    _Sd[d]=(Sn(n)-ans+P)%P;
}

void clat_m(int n,int d)
{
    int ans=0,m=sqrt(n)+1;
    for(int L=1;L<m;L++)
    {
        int u=n/L-n/(L+1);
        ans+=(LL)M[L]*u%P;
        if(ans>=P)ans-=P;
    }
    for(int i=n/m;i>1;i--)
    {
        int u=n/i;
        if(u<MAXN)
            ans+=M[u];
        else
            ans+=_M[i*d];
        if(ans>=P)ans-=P;
    }
    _M[d]=(1-ans+P)%P;
}

void clat_p(int n,int d)
{
    int ans=0,m=sqrt(n)+1;
    for(int L=1;L<m;L++)
    {
        ans+=(LL)Sp[L]*((Sn(n/L)-Sn(n/(L+1))+P)%P)%P;
        if(ans>=P)ans-=P;
    }
    for(int i=n/m;i>1;i--)
    {
        int u=n/i;
        if(u<MAXN)
            u=Sp[u];
        else
            u=_Sp[i*d];
        ans+=(LL)i*u%P;
        if(ans>=P)ans-=P;
    }
    _Sp[d]=(1-ans+P)%P;
}

void slove_Sd(int n)
{
    if(n<MAXN)return ;
    for(int i=n/MAXN;i;i--) clat_d(n/i,i);
}

void slove_M(int n)//最先计算
{
    if(n<MAXN)return ;
    for(int i=n/MAXN;i;i--) clat_m(n/i,i);
}

void slove_Sp(int n)
{
    if(n<MAXN)return ;
    for(int i=n/MAXN;i;i--) clat_p(n/i,i);
}

int main ()
{
    init();
    int n;
    scanf("%d",&n);
    slove_M(n);
    slove_Sd(n);
    slove_Sp(n);
    int m=sqrt(n)+1,ans=0;
    for(int L=1;L<m;L++)
    {
        int f1=n/L,f2=n/(L+1),u=(LL)Sd[L]*Sd[L]%P;
        if(f1<MAXN)
            f1=Sp[f1];
        else
            f1=_Sp[L];
        if(f2<MAXN)
            f2=Sp[f2];
        else
            f2=_Sp[L+1];
        ans+=(LL)u*(f1-f2+P)%P;
        if(ans>=P)ans-=P;
        if(ans<0)ans+=P;
    }
    for(int i=n/m;i;i--)
    {
        int u=n/i;
        if(u<MAXN)
            u=Sd[u];
        else
            u=_Sd[i];
        u=(LL)u*u%P;
        ans+=((LL)phi_[i]*u%P+P)%P;
        if(ans>=P)ans-=P;
    }
    printf("%d\n",ans);
}