[51Nod 1238]最小公倍数之和 V3 杜教筛

传送门·Teleport!

题目大意

给定N,求
$$ \sum_{i=1}^n \sum_{j=1}^n \frac{ij}{gcd(i,j)} $$
N<=10^10,答案对10^9+7取模。

解题报告

好 毒瘤 题啊!

首先设[latex] Ans = \sum_{i=1}^n \sum_{j=1}^n \frac{ij}{gcd(i,j)} [/latex]
然后设[latex] F(n) = \sum_{i=1}^n \sum_{j=1}^{i} \frac{ij}{gcd(i,j)} [/latex]
那么有
$$ Ans = 2F(n)-\frac{(n+1)n}{2} $$
我们知道小于i并且与i互质的数之和的公式为
$$ \frac{1}{2}\varphi(i) \times i $$
因为如果j和i互质,那么i-j也一定和i互质,那么这些和i互质的数字是成对出现的,并且每一对的和为i。
因为互质的数字一共有phi(i)个,对数就是1/2phi(i)个,就出来了。
特别的,如果i=1需要特殊考虑。
然后化简[latex] F(n) [/latex]
$$ F(n) = \sum_{i=1}^n \sum_{j=1}^{i} \frac{ij}{gcd(i,j)} = \sum_{i=1}^n \frac{1}{2} (\sum_{d|i} \varphi(\frac{i}{d}) \times \frac{i}{d} + 1) $$
那么有:
$$ Ans = 2F(n) – \frac{(n+1)n}{2} = \
\sum_{i=1}^n i \times \sum_{d|i} {\varphi(\frac{i}{d})\times \frac{i}{d}} = \sum_{i=1}^n i \times \sum_{d|i} {\varphi(d)\times d} $$

然后介绍杜教筛的套路:
定理:
如果有数论积性函数[latex]f(x)[/latex]
那么对于任意数论函数[latex]g(x)[/latex]有:
$$ \sum_{i=1}^{n} \sum_{d|i} f(d)g(\frac{i}{d}) = \sum_{i=1}^n S(\lfloor \frac{n}{i} \rfloor) g(i) $$
其中,[latex]S(n)=\sum_{i=1}^n f(i)[/latex] (说白了就是数论函数f(x)的前缀和。)
如果我们能快速计算等式左边的内容的话,用上面式子移向之后就可以推出[latex]S(n)[/latex]的递推式了。
那么关键在于对数论函数[latex]g(x)[/latex]的选择。

回到题目,发现如果令[latex]f(x)=\varphi(d)\times d[/latex],[latex]g(x)=1[/latex],那么
$$ Ans = \sum_{i=1}^n i \times \sum_{d|i} {\varphi(d)\times d} = \sum_{i=1}^n i \times \sum_{d|i} f(d)g(\frac{i}{d}) $$
继续推导(这一步不是很懂。。):
$$ = \sum_{i=1}^n i \sum_{j=1}^{\lfloor \frac{n}{i} \rfloor} i \times f(j)g(j) $$
然后我们发现n/i只有根号n个取值,第一个sigma就可以分块处理,接下来再用杜教筛处理后面部分:
$$ \sum_{i=1}^n f(i)g(i) = \sum_{i=1}^n \varphi(i)\times i^2 $$
我们在用杜教筛的那个套路,取[latex]g(x)=x^2[/latex]来消掉[latex]i^2[/latex]
然后:
$$ \sum_{i=1}^n \sum_{d|i} \varphi(d) \times d^2 \times \frac{i^2}{d^2} = \sum_{i=1}^n S(\lfloor \frac{n}{i} \rfloor) \times i^2 $$
化简移项:
$$ \sum_{i=1}^n i^3 – \sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor) \times i^2 = S(n)$$
然后立方和是可以公式O(1)算的,后面递推可以一样分块处理,我们就可以记忆化递归地计算S(n)了,这里发现需要计算的S(n)和前面的分块需要的S(n)有很多重复,复杂度就得到保障(但是我并不会分析复杂度。。。)。

然后算就可以了。。。。

最后说一句,这题是个多倍经验题,肛完之后可以去做51Nod上的最大公因数之和,平均公倍数之和,AC之后杜教筛也能算是入门了吧。。
最进玩JeremyGuo出的一堆杜教筛都要被玩坏了。。。

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <climits>
#include <cstdlib>
#include <map>
typedef long long LL;
using namespace std;
const int MAXL = 5005000;
const int MOD = 1000000007;
LL phi[MAXL], sum[MAXL];
bool prime[MAXL];
int prec[MAXL/6], ptot, inilen;

void init (int len) {
    inilen = len;
    phi[1] = 1;
    for(int i = 2; i<=len; i++) {
        if(!prime[i]) phi[i] = i-1, prec[++ptot] = i;
        for(int j = 1; j<=ptot && prec[j]*i <= len; j++) {
            prime[i*prec[j]] = true;
            if(i%prec[j] == 0) {phi[i*prec[j]] = phi[i]*prec[j]; break;}
            phi[i*prec[j]] = phi[i]*(prec[j]-1);
        }
    }
    for(int i = 1; i<=len; i++)
        sum[i] = (sum[i-1]+phi[i]*i%MOD*i%MOD)%MOD;
}
const LL inv2 = 500000004, inv6 = 166666668;
LL Sum1 (LL x) {
    return x%MOD*((x+1)%MOD)%MOD*inv2%MOD;
}
LL Sum2 (LL x) {
    return x%MOD*((x+1)%MOD)%MOD*((2*x+1)%MOD)%MOD*inv6%MOD;
}
LL Sum3 (LL x) {
    return Sum1(x)*Sum1(x)%MOD;
}
map<LL, LL> rec;
LL getSumFunc (LL x) {
    if(x <= inilen) return sum[x];
    if(rec.find(x) != rec.end()) return rec[x];
    LL ret = Sum3(x);
    for(LL i = 2; i<=x; i++) {
        LL nxt = x/(x/i);
        LL intl = Sum2(i-1), intr = Sum2(nxt);
        LL sig = (intr-intl+MOD)%MOD;
        ret = ((ret-sig*getSumFunc(x/i))%MOD+MOD)%MOD;
        i = nxt;
    }
    return rec[x] = ret;
}

LL solve (LL x) {
    LL ret = 0;
    for(LL i = 1; i<=x; i++) {
        LL nxt = x/(x/i);
        LL sig = (Sum1(nxt)-Sum1(i-1)+MOD)%MOD;
        ret = (ret+sig*getSumFunc(x/i)%MOD)%MOD;
        i = nxt;
    }
    return ret;
}

int main () {
    LL x;
    init(MAXL-100);
    scanf("%lld", &x);
    printf("%lld\n", solve(x));
}