51nod 1220 约数之和 杜教筛

题意

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)的结果。
2 <= N <= 10^9

分析

设\(\sigma(i)\)表示i的约数和。
这题要用到一个结论，就是\[\sigma(ij)=\sum\limits_{p|i}\sum\limits_{q|j}[(p,q)=1]\frac{pj}{q}\]。
证明的话，可以考虑每个素因子分开考虑。设\(r\)为素数，\(i=r^a,j=r^b\)。当\(p=r^0\)时\(\frac{pj}{q}\)可以取到\(r^0..r^b\)；反之当\(q=r^0\)时\(\frac{pj}{q}\)可以取到\(r^{b+1}..r^{a+b}\)。
那么有\[ans=\sum_{i=1}^n\sum_{j=1}^n\sum_{p|i}\sum_{q|j}[(p,q)=1]\frac{pj}{q}\]
反演一下可以得到\[ans=\sum_{d=1}^n\mu(d)\sum_{i=1}^n\sum_{j=1}^n\sum_{p|i}\sum_{q|j}[d|(p,q)]\frac{pj}{q}\]
\[=\sum_{d=1}^n\mu(d)\sum_{d|p}\sum_{d|q}\frac{p}{q}\sum_{p|i}\sum_{q|j}j\]
\[=\sum_{d=1}^n\mu(d)\sum_{d|p}\sum_{d|q}\frac{p}{q}\lfloor\frac{n}{p}\rfloor q\frac{\lfloor\frac{n}{q}\rfloor(\lfloor\frac{n}{q}\rfloor+1)}{2}\]
\[=\sum_{d=1}^n\mu(d)\sum_{d|p}p\lfloor\frac{n}{p}\rfloor\sum_{d|q}\frac{\lfloor\frac{n}{q}\rfloor(\lfloor\frac{n}{q}\rfloor+1)}{2}\]
\[=\sum_{d=1}^n\mu(d)d(\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}p\lfloor\frac{n}{dp}\rfloor)(\sum_{q=1}^{\lfloor\frac{n}{q}\rfloor}\frac{\lfloor\frac{n}{dq}\rfloor(\lfloor\frac{n}{dq}\rfloor+1)}{2})\]
这已经是一个分块的式子了，考虑把下取整里面的东西换一下。
\[\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\]
\[\sum_{i=1}^n\frac{\lfloor\frac{n}{i}\rfloor(\lfloor\frac{n}{i}\rfloor+1)}{2}\]
也就是说我们要考虑如何快速求这两个东西。
其中第一个式子的下取整可以看做n以内有多少个数是i的倍数，那么可以得到\[\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor=\sum_{i=1}^n\sigma(i)\]
对第二个式子进行变形：
\[\sum_{i=1}^n\frac{\lfloor\frac{n}{i}\rfloor(\lfloor\frac{n}{i}\rfloor+1)}{2}=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j=\sum_{i=1}^ni\lfloor\frac{n}{i}\rfloor\]
也就是说上面两个式子的值居然是相等的！
那么就有\[ans=\sum_{d=1}^n\mu(d)d(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sigma(i))^2\]
前面那部分的前缀和可以用杜教筛来求，后面那部分就小的预处理大的分块求即可。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;

typedef long long LL;

const int N=2000005;
const int MOD=1000000007;

int n,prime[N],yue[N],low[N],s[N],tot;
bool not_prime[N];
map<int,int> w;

void get_prime(int n)
{
    yue[1]=s[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i]) prime[++tot]=i,yue[i]=i+1,s[i]=-i,low[i]=i;
        for (int j=1;j<=tot&&i*prime[j]<=n;j++)
        {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                low[i*prime[j]]=low[i]*prime[j];
                yue[i*prime[j]]=(LL)yue[i]*prime[j]%MOD+yue[i/low[i]];
                yue[i*prime[j]]-=yue[i*prime[j]]>=MOD?MOD:0;
                break;
            }
            low[i*prime[j]]=prime[j];
            yue[i*prime[j]]=(LL)yue[i]*yue[prime[j]]%MOD;
            s[i*prime[j]]=s[i]*s[prime[j]];
        }
    }
    for (int i=1;i<=n;i++) s[i]+=s[i]<0?MOD:0,s[i]+=s[i-1],s[i]-=s[i]>=MOD?MOD:0;
    for (int i=1;i<=n;i++) yue[i]+=yue[i-1],yue[i]-=yue[i]>=MOD?MOD:0;
}

int get_s(int n)
{
    if (n<=2000000) return s[n];
    if (w[n]) return w[n];
    int ans=1;
    for (int i=2,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=MOD-(LL)((LL)last*(last+1)/2-(LL)i*(i-1)/2)%MOD*get_s(n/i)%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return w[n]=ans;
}

int get_yue(int n)
{
    if (n<=2000000) return yue[n];
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ans+=(LL)((LL)last*(last+1)/2-(LL)i*(i-1)/2)%MOD*(n/i)%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return ans;
}

int solve(int n)
{
    int ans=0;
    for (int i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);int w=get_yue(n/i);
        ans+=(LL)(get_s(last)+MOD-get_s(i-1))*w%MOD*w%MOD;
        ans-=ans>=MOD?MOD:0;
    }
    return ans;
}

int main()
{
    get_prime(2000000);
    scanf("%d",&n);
    printf("%d",solve(n));
    return 0;
}

转载于:https://www.cnblogs.com/beginend/p/8240908.html