「雅礼集训 2017 Day11」DIV

一道基本的杜教筛。

题目大意

给定 n∈Zn\in\mathbb{Z}nZz∈Z[ı]z\in \mathbb{Z}[\imath]zZ[ı],求

∑i=1n∑z∣i,Re(z)>0Re(z)\sum_{i=1}^n\sum_{z\mid i,\mathrm{Re}(z)>0}\mathrm{Re}(z)i=1nzi,Re(z)>0Re(z)

答案模 100453580910045358091004535809 输出。

n≤1011n\le10^{11}n1011

题目链接

LYOJ #219

题解

首先考虑这个复数肯定要稍微转化一下。设 n=(a+bı)(c+dı)n=(a+b\imath)(c+d\imath)n=(a+bı)(c+dı),分别考虑实部和虚部,得到:

acbd=nad+bc=0

而后者其实就是

ab=−cd\frac ab=-\frac cdba=dc

我们设这个比例为 pq\displaystyle \frac pqqp,其中 (p,q)=1(p,q)=1(p,q)=1,则一定有如下形式:

n=x(p+qı)⋅y(p−qı) n=x(p+q\imath)\cdot y(p-q\imath)n=x(p+qı)y(pqı)

容易发现这个式子前后两项是对称的,因此我们先只考虑前面一项,即 Im(z)>0\mathrm{Im}(z)>0Im(z)>0zzz 对答案的贡献。这时上式等价于

n=xy(p2+q2)n=xy\left(p^2+q^2\right)n=xy(p2+q2)

因此一对整数 p,qp,qp,q 对答案有贡献当且仅当 p2+q2∣np^2+q^2\mid np2+q2n,且贡献大小为

p⋅∑x∣np2+q2x=p⋅σ(np2+q2)p\cdot\sum_{x\mid\frac n{p^2+q^2}}x=p\cdot\sigma\left(\frac n{p^2+q^2}\right)pxp2+q2nx=pσ(p2+q2n)

其中 σ(n)\sigma(n)σ(n) 代表 nnn 的约数之和。枚举所有的 p,qp,qp,q,对答案的总贡献就是

∑(p,q)=1p∑p2+q2∣iσ(ip2+q2)\sum_{(p,q)=1}p\sum_{p^2+q^2\mid i}\sigma\left(\frac i{p^2+q^2}\right)(p,q)=1pp2+q2iσ(p2+q2i)

考虑到求和的项中出现了 np2+q2\displaystyle\frac{n}{p^2+q^2}p2+q2n,我们枚举 k=p2+q2k=p^2+q^2k=p2+q2,有总贡献为

k(p,q)=1p2+q2=kpkiσ(ik)

可以发现最后一层求和的结果实际上是一个 σ\sigmaσ 的前缀和,因此记 D(n)=∑i=1nσ(i)D(n)=\sum_{i=1}^n\sigma(i)D(n)=i=1nσ(i),表达式简化为

k(p,q)=1p2+q2=kpD(nk)

由于 nk\displaystyle\frac nkkn 的取值只有 O(n1/2)O\left(n^{1/2}\right)O(n1/2) 种,考虑枚举之,这时每次考虑的 kkk 形成一个连续区间,如果能快速计算 (p,q)=1p2+q2=kp 的前缀和与 D(n)D(n)D(n),我们就可以快速求出答案了。

f(n)=(p,q)=1p2+q2=npF(n)=(p,q)=1p2+q2np

F(n)F(n)F(n)f(n)f(n)f(n) 的前缀和。考虑去除 pppqqq 互质的限制,设

G(n)=∑p2+q2≤npG(n)=\sum_{p^2+q^2\le n}pG(n)=p2+q2np

则若枚举 d=(p,q)d=(p,q)d=(p,q),则可以发现 F(n)F(n)F(n)G(n)G(n)G(n) 间存在如下关系:

G(n)=∑d≥1F(⌊nd2⌋)dG(n)=\sum_{d\ge1}F\left(\left\lfloor\frac n{d^2}\right\rfloor\right)dG(n)=d1F(d2n)d

d=1d=1d=1 的情况单独考虑,可以得到 F(n)F(n)F(n) 的递推式:

F(n)=G(n)−∑d≥2F(⌊nd2⌋)dF(n)=G(n)-\sum_{d\ge2}F\left(\left\lfloor\frac n{d^2}\right\rfloor\right)dF(n)=G(n)d2F(d2n)d

G(n)G(n)G(n) 可以依下式在 O(n1/2)O\left(n^{1/2}\right)O(n1/2) 的时间复杂度内计算:

G(n)=∑p=1⌊n1/2⌋p⋅⌊n−p2⌋G(n)=\sum_{p=1}^{\left\lfloor n^{1/2}\right\rfloor}p\cdot\left\lfloor \sqrt{n-p^2}\right\rfloorG(n)=p=1n1/2pnp2

则计算 F(n)F(n)F(n) 的时间复杂度为 O(n1/2)O\left(n^{1/2}\right)O(n1/2)

而关于 D(n)D(n)D(n) 有经典结论

D(n)=∑i=1nσ(i)=∑d=1nd⌊nd⌋D(n)=\sum_{i=1}^n\sigma(i)=\sum_{d=1}^nd\left\lfloor\frac n{d}\right\rfloorD(n)=i=1nσ(i)=d=1nddn

如果枚举 ⌊nd⌋\displaystyle\left\lfloor\frac n{d}\right\rfloordn,则 D(n)D(n)D(n) 也可以在 O(n1/2)O\left(n^{1/2}\right)O(n1/2) 的时间内求出。

若直接计算所有的函数值,则总时间复杂度为 O(n3/4)O\left(n^{3/4}\right)O(n3/4)。考虑利用杜教筛的思想,预处理前 n2/3n^{2/3}n2/3 个点的函数值,则 D(n)D(n)D(n) 可以利用 Euler 筛在 O(n2/3)O\left(n^{2/3}\right)O(n2/3) 的时间内预处理,而 F(n)F(n)F(n) 则可以通过直接枚举 pppqqqO(n2/3logn)O\left(n^{2/3}\log n\right)O(n2/3logn) 的时间内简单预处理出。具体实现可以参考代码。

这样,我们就计算出了满足 Im(z)>0\mathrm{Im}(z)>0Im(z)>0 的所有 zzz 的贡献。Im(z)<0\mathrm{Im}(z)<0Im(z)<0 时的情况是对称的,贡献与 Im(z)>0\mathrm{Im}(z)>0Im(z)>0 时相等。而当 Im(z)=0\mathrm{Im}(z)=0Im(z)=0 时,zzz 为实数,这一部分的贡献为 ∑i=1nσ(i)=D(n)\sum_{i=1}^n\sigma(i)=D(n)i=1nσ(i)=D(n)。将这几部分的贡献相加,我们就在 O(n2/3logn)O\left(n^{2/3}\log n\right)O(n2/3logn) 的时间复杂度内计算出了答案。

UPD:发现这题 Project Euler 上有原题 Project Euler 153 Investigating Gaussian Integers,不过是 O(n)O(n)O(n) 可以跑的……大家快去水一水呀

代码

#include <bits/stdc++.h>
#define debug(x) std::cerr << #x << " : " << (x) << std::endl

typedef long long LL;

const int MAXR = 100031, MAXT = 5000000, P = 1004535809, I2 = 502267905;

inline LL mul(LL a, LL b) {
    return a %= P, b %= P, a * b % P;
}

LL n = 0, rt = 0, c = 0;

bool ip[MAXT];
int p[MAXT];
LL ds[MAXT], fs[MAXT], d[2 * MAXR], f[2 * MAXR];
int ind[2][MAXR];

void sieve() {
    ds[1] = 1;
    for (int i = 2, cp = 0; i < rt; i++) {
        if (!ip[i]) p[cp++] = i, ds[i] = i + 1;
        for (int j = 0; j < cp && i * p[j] < rt; j++) {
            ip[i * p[j]] = 1;
            if (i % p[j]) ds[i * p[j]] = ds[i] * (1 + p[j]);
            else {
                ds[i * p[j]] = (p[j] + 1) * ds[i] - p[j] * ds[i / p[j]];
                break;
            }
        }
    }
    for (LL i = 1; i * i < rt; i++) {
        LL t = i * i;
        for (LL j = 1; j * j + t < rt; j++) if (std::__gcd(i, j) == 1) fs[j * j + t] += i;
    }
    for (int i = 1; i < rt; i++) ds[i] = (ds[i] + ds[i - 1]) % P, fs[i] = (fs[i] + fs[i - 1]) % P;
}

int &id(LL x) {
    return x < MAXR ? ind[0][x] : ind[1][n / x];
}

void calc() {
    for (LL i = 1; i <= n; i++) {
        LL z = n / (n / i);
        if (z < rt) d[c] = ds[z], f[c] = fs[z];
        else {
            for (LL a = 1; a * a <= z; a++) f[c] = (f[c] + a * (LL)floor(sqrt(z - a * a))) % P;
            for (LL a = 2; a * a <= z; a++) f[c] = (f[c] - mul(f[id(z / (a * a))], a) + P) % P;
            for (LL j = 1; j <= z; j++) {
                LL t = z / (z / j);
                d[c] = (d[c] + mul(mul(mul(j + t, t - j + 1), I2), z / t)) % P;
                j = t;
            }
        }
        id(i = z) = c++;
    }
    id(0) = c;
}

int main() {
    scanf("%lld", &n);
    rt = exp(log(n) / 3 * 2) + 1;
    sieve(), calc();
    LL ans = 0;
    for (LL i = 1; i <= n; i++) {
        LL t = n / (n / i);
        ans = (ans + mul(f[id(t)] - f[id(i - 1)] + P, d[id(n / t)])) % P;
        i = t;
    }
    printf("%lld\n", (ans * 2 + d[id(n)]) % P);
    return 0;
}

5 thoughts on “「雅礼集训 2017 Day11」DIV

  1. 真的是...基本吗... 我可能是太弱了

  2. 钊爷,如何jjj是一个⌊n/j⌋\lfloor n/j \rfloorn/j的分割点, 那么⌊j/a⌋\lfloor j/a \rfloorj/a也一定是分割点咩?

    • 是的。 假设 n,i,j∈Zn,i,j\in\mathbb{Z}n,i,jZ,那么有 ⌊⌊ni⌋/j⌋=⌊nij⌋\left\lfloor\left\lfloor\displaystyle\frac n i\right\rfloor/j\right\rfloor=\left\lfloor\displaystyle\frac n {ij}\right\rfloorin/j=ijn

发表评论

电子邮件地址不会被公开。 必填项已用*标注