51nod 1220 约数之和

Posted 2020-11-30 tongseli

约数之和

Problem

d(k)表示k的所有约数的和。d(6) = 1 + 2 + 3 + 6 = 12。

定义S(N) = ∑1<=i<=N ∑1<=j<=N d(i*j)。

例如：S(3) = d(1) + d(2) + d(3) + d(2) + d(4) + d(6) + d(3) + d(6) + d(9) = 59，S(1000) = 563576517282。

给出正整数N，求S(N)，由于结果可能会很大，输出Mod 1000000007(10^9 + 7)的结果。

输入

输入一个数N(2 <= N <= 10^9)。

输出

输出S(N) Mod 1000000007(10^9 + 7)的结果。

输入样例

1000

输出样例

576513341

Guess

首先得会把(d(ij)) 写成(gcd) 的形式

然后搞一波莫比乌斯变换

再用和式化一化

上个杜教筛，打个分块

好像就能解决了！

Solve

就是推公式呗。

这里这个(d(ij)) 与之前的不太一样，这里是约数之和的函数。

必须要知道这个式子才可以往后推：

[d(ij)=sumlimits_{x|i}sumlimits_{y|j}[gcd(x,y)==1]xcdotfrac{j}{y}]

（意会一下这个式子，大概和之前的除数函数相类似，枚举质因子分析的感觉）
（我不懂怎么来的）

这样的话

[ egin{align} S(n) &= sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}d(ij) &= sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}sumlimits_{x|i}sumlimits_{y|j}[gcd(x,y)==1]xcdotfrac{j}{y} & 来一波莫比乌斯变换吧！ &=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}sumlimits_{x|i}sumlimits_{y|j}sumlimits_{d|gcd(x,y)}mu (d)cdot xcdotfrac{j}{y} &=sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}sumlimits_{x|i}sumlimits_{y|j}sumlimits_{d=1}^nmu (d)[d|gcd(x,y)]cdot xcdotfrac{j}{y} &我们知道d|gcd(x,y)=d|x 且d|y &=sumlimits_{d=1}^nmu (d)sumlimits_{i=1}^{n}sumlimits_{j=1}^{n}sumlimits_{x|i}sumlimits_{y|j}[d|x][d|y]cdot xcdotfrac{j}{y}&对于枚举1到n的数i的因子x的因子d&其实就是枚举1到n的因子d的倍数x的倍数i&=sumlimits_{d=1}^nmu (d)sumlimits_{d|x}sumlimits_{d|y}frac{x}{y}sumlimits_{x|i}sumlimits_{y|j}j&这里的sumlimits_{y|j}j=y+2y+3y+cdots+lfloorfrac{n}{y} floor y&=sumlimits_{d=1}^nmu (d)sumlimits_{d|x}sumlimits_{d|y}frac{x}{y}lfloorfrac{n}{x} floorcdot y frac{lfloorfrac{n}{y} floor(lfloorfrac{n}{y} floor+1)}{2} &=sumlimits_{d=1}^nmu (d)sumlimits_{d|x}xlfloorfrac{n}{x} floorsumlimits_{d|y}frac{lfloorfrac{n}{y} floor(lfloorfrac{n}{y} floor+1)}{2} &=sumlimits_{d=1}^nmu (d)sumlimits_{x=1}^{lfloorfrac{n}{d} floor}dxlfloorfrac{n}{dx} floorsumlimits_{y=1}^{lfloorfrac{n}{y} floor}frac{lfloorfrac{n}{dy} floor(lfloorfrac{n}{dy} floor+1)}{2} &=sumlimits_{d=1}^nmu (d)dsumlimits_{x=1}^{lfloorfrac{n}{d} floor}xlfloorfrac{n}{dx} floorsumlimits_{y=1}^{lfloorfrac{n}{y} floor}frac{lfloorfrac{n}{dy} floor(lfloorfrac{n}{dy} floor+1)}{2} \end{align} ]
而对于中间的和号可以看成是(f(lfloorfrac{n}{d} floor))

[f(n)=sumlimits_{x=1}^{n}xlfloorfrac{n}{x} floor]

对于后面的和号看成是(g(lfloorfrac{n}{d} floor))

[ egin{align} g(n)&=sumlimits_{y=1}^{n}frac{lfloorfrac{n}{y} floor(lfloorfrac{n}{y} floor+1)}{2}&等差数列前n项和&=sumlimits_{y=1}^{n}sumlimits_{i=1}^{lfloorfrac{n}{y} floor}i&=sumlimits_{y=1}^{n}ylfloorfrac{n}{y} floor end{align} ]
所以这两个式子是一样的。

[S(n)=sumlimits_{d=1}^nmu (d)cdot dcdot (sumlimits_{x=1}^{lfloorfrac{n}{d} floor}xlfloorfrac{n}{dx} floor)^2 ]

这样已经可以写了，但是我用的是这个式子

[S(n)=sumlimits_{d=1}^nmu (d)cdot dcdot (sumlimits_{x=1}^{lfloorfrac{n}{d} floor}sumlimits_{i=1}^{lfloorfrac{n}{dx} floor}i)^2 ]

1. 后面的就是分块

2. 前面的用杜教筛（其实就是构造(h=f*g)）

3. 外面再套一个分块

补充杜教筛的推导

要求(A(n)=sumlimits_{i=1}^nicdot mu(i))

我们可以设(f(i)=icdot mu(i))

然后开始构造
[ egin{align} h&=f*g&=sumlimits_{d|n}dcdot mu(d) cdot g(n/d)&可以取g(n/d)=n/d&=nsumlimits_{d|n}mu(d)&然后我们知道mu *u=I&=n[n==1]&=[n==1]=I end{align} ]
套一下杜教筛的公式
[A(n)=1-sumlimits_{i=2}^ndcdot A(lfloorfrac{n}{d} floor)]

Code

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;

const int mod = 1e9 + 7;
const int maxn = 6e6;
const int N = 6e6 + 10;
const int inv2 = (mod+1)/2;

int prime[N];
int vis[N];
int mu[N];
ll sum[N];


void sieve()
{
    int p=0;
    mu[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!vis[i]){
            prime[++p]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=p&&i*prime[j]<=maxn;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)break;
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=maxn;i++){
        sum[i]=(sum[i-1]+1ll*mu[i]*i)%mod;
    }
}

unordered_map<int, ll> mp;

ll S(int l, int r){
    return (1ll*(1ll*(l+r)%mod*(r-l+1)%mod)%mod*inv2)%mod;
}

ll djs(int n){
    if(n<=maxn)return sum[n];
    if(mp[n])return mp[n];
    ll ans=1+mod;
    for(int l=2,r;l>=0&&l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans-1ll*S(l,r)%mod*djs(n/l)%mod+mod)%mod;
    }
    return mp[n]=ans;
}

ll F(int x){
    ll ans=0;
    for(int l=1,r;l<=x;l=r+1){
        r=x/(x/l);
        ans+=1ll*(r-l+1)*S(1,x/l)%mod;
    }
    return ans%mod;
}

int main()
{
    int n;
    scanf("%d",&n);
    sieve();
    ll ans=0;
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans%mod+((djs(r)-djs(l-1))%mod)%mod*(F(n/l)%mod*F(n/l)%mod)%mod)%mod;
        ans=(ans+mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}