洛谷3768 简单的数学题 莫比乌斯反演+杜教筛

Description

---

由于出题人懒得写背景了，题目还是简单一点好。
输入一个整数n和一个整数p，你需要求出$\sum_{i=1}^n\sum_{j=1}^n i\times j\times gcd(i,j)\pmod p$

对于20%的数据， $n \leq 1000$ 。
对于30%的数据， $n \leq 5000$ 。
对于60%的数据， $n \leq 10^6$ ，时限1s。
对于另外20%的数据， $n \leq 10^9$ ，时限3s。
对于最后20%的数据， $n \leq 10^{10}$ ，时限6s。
对于100%的数据，$5 \times 10^8 \leq p \leq 1.1 \times 10^9$ 且p为质数。

Solution

---

感觉对杜教筛的理解又深一点点

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{n}i\times j\times gcd(i,j)=\sum_{d=1}^{n}d^3\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}i\times j\times [gcd(i,j)==1]$$
第一波操作照例是枚举gcd，然后再套μ得到 $$ans=\sum_{d=1}^{n}d^3\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}x^2\mu(x)S^2(\lfloor\frac{n}{dx}\rfloor)$$
其中$S(n)=\sum_{i=1}^{n}i$
令$T=dx$枚举$T$得到 $$ans=\sum_{T=1}^{n}S^2(\lfloor\frac{n}{T}\rfloor)T^2\sum_{d|T}\mu(d)\frac{T}{d}$$
有一个结论是$\sum_{d|T}\mu(d)\frac{T}{d}=\varphi(T)$，知道了这个就可以继续往下推惹
$$ans=\sum_{T=1}^{n}S^2(\lfloor\frac{n}{T}\rfloor)T^2\varphi(T)$$
这里已经可以分块求，但是还需解决后面的前缀和问题，考虑杜教筛
设$f(i)=i^2\varphi(i)$，$g(i)=i^2$，$h(n)=\sum_{i=1}^{n}f(i)$
容易得到 $$\sum_{i=1}^{n}\sum_{d|i}g(d)\times f(\frac{i}{d})=\sum_{i=1}^{n}\sum_{d|i}d^2\times ({\frac{i}{d}})^2\times \varphi(\frac{i}{d})=\sum_{i=1}^{n}i^2\sum_{d|i}\varphi(d)=\sum_{i=1}^{n}i^3$$ 同时
$$\sum_{i=1}^{n}\sum_{d|i}g(d)\times f(\frac{i}{d})=\sum_{i=1}^{n}\sum_{d|i}d^2\times ({\frac{i}{d}})^2\times \varphi(\frac{i}{d})=\sum_{d=1}^{n}d^2\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}i^2\times \varphi(i)=\sum_{d=1}^{n}d^2\times h(\lfloor\frac{n}{d}\rfloor)$$
当$i=1$时拉出来得到 $$h(n)=\sum_{i=1}^{n}i^3-\sum_{d=2}^{n}d^2h(\lfloor\frac{n}{d}\rfloor)$$
那么就可以杜教筛直接上了

需要指出的是，$\sum_{i=1}^{n}i^3=({\sum_{i=1}^{n}i})^2$，我怕不是小学奥数题都不会了
这题的模数不一定是质数因此不能用费马小定理求逆元，考虑到n不是特别大可以直接除

Code

---

#include <stdio.h>
#include <map>

typedef long long LL;
const int N=10000000;

std:: map <LL,LL> map;

int prime[N+5],MOD;
bool not_prime[N+5];
LL f[N+5];

void pre_work() {
    f[1]=1;
    for (int i=2;i<=N;i++) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            f[i]=i-1;
        }
        for (int j=1;i*prime[j]<=N&&j<=prime[0];j++) {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) {
                f[i*prime[j]]=f[i]*prime[j];
                break;
            }
            f[i*prime[j]]=f[i]*f[prime[j]];
        }
    }
    for (int i=1;i<=N;i++) f[i]=(f[i-1]+f[i]*i%MOD*i%MOD)%MOD;
}

LL S(LL n) {
    n=n%MOD;
    return (LL)n*(n+1)/2%MOD;
}

LL S2(LL n)
{
    n=n%MOD;
    if ((LL)n*(n+1)%3==0) return (LL)n*(n+1)/6%MOD*(n*2+1)%MOD;
    else return (LL)(n*2+1)/3*((LL)n*(n+1)/2%MOD)%MOD;
}

LL get_sum(LL n) {
    if (n<=N) return f[n];
    if (map[n]) return map[n];
    LL ret=S(n); ret=ret*ret%MOD;
    for (LL i=2,j;i<=n;i=j+1) {
        j=n/(n/i);
        ret=(ret-(S2(j)-S2(i-1)+MOD)%MOD*get_sum(n/i)%MOD+MOD)%MOD;
    }
    return map[n]=ret;
}

void solve(LL n) {
    LL ans=0;
    for (LL i=1,j;i<=n;i=j+1) {
        j=n/(n/i);
        LL sum1=S(n/i);
        ans=(ans+(sum1*sum1)%MOD*(get_sum(j)-get_sum(i-1)+MOD)%MOD)%MOD;
    }
    printf("%lld\n", (ans+MOD)%MOD);
}

int main(void) {
    LL n; scanf("%d%lld",&MOD,&n);
    pre_work();
    solve(n);
    return 0;
}