SP33039 AFS3 - Amazing Factor Sequence

· · 题解

本文主要源于 Min_25 的博客

我们考虑一般的情况

S_k(n)=\sum_{i=1}^n\sigma_k(i)=\sum_{i=1}^ni^k\lfloor\frac{n}{i}\rfloor

介绍一种 \Omega(k^2n^{\frac{1}{3}}\log n) 的做法我们发现当 k = 0 的时候这个复杂度就是 0 !

所以在讨论这个一般情况的时候请默认 k\geq 1

首先让我们回忆 DIVCNT1 是怎么解决的 详细地可以去看一下 这个博客 这里我将只提及主要思想

对于 DIVCNT1 等价于双曲线下整点计数

我们考虑用一个凸包去拟合这个区域 得到凸包之后可以通过 Pick 定理 得到内部的顶点数量

现在我们的问题转化成了如何选择这个凸包 然后我们需要一个性质:

本质不同的斜率个数是 \Omega(n^{\frac{1}{3}}) 的

证明有 论文

考虑怎么得到这些斜率 这时候有一个奇妙的数据结构出现了 Stern-Brocot Tree 这个数据结构可以不重不漏地构造既约真分数 接下来我们可以通过在 SBT 上二分得到斜率 这个过程用一个单调栈维护即可 所以我们最终复杂度是 \Omega(n^{\frac{1}{3}}\log n)

下面我们默认读者已经完全明白了 DIVCNT1 是怎么做的 因为 SBT 部分是几乎相同的

回到本题 我们照样去构造曲线用凸包拟合

我们考虑这样记录 S_k(n) 上的点 (a,b,T_0(a,b),T_1(a,b),\cdots, T_k(a,b))

其中 T_k(a,b)=\sum_{i=0}^{b-1}i\lfloor\frac{ai}{b}\rfloor^k

那么当我们用凸包去拟合 S_k(n) 的时候对于 一个梯形区域 我们可以直接通过 T_0(a,b),T_1(a,b),\cdots,T_k(a,b) 和 P_0(b),P_1(b),\cdots,P_k(b) 在 \mathcal O(k) 内得到 (Pick 定理) 其中 P_k(b)=\sum_{i=0}^{b-1}i^k 而我们计算 P_k(b) 需要 \mathcal O(k^2)

然后对于 SBT 部分 直接维护:

(a,b,T_0(a,b),T_1(a,b),\cdots,T_k(a,b))+(c,d,T_0(a,b),T_1(c,d),\cdots,T_k(c,d))\\\to\\(a + c,b + d,T_0(a+c,b+d),T_1(a + c,b + d),\cdots,T_k(a + c,b + d))

每次在 SBT 上二分时直接暴力卷积即可 这部分 \mathcal O(k^2)

由于这题细节过多 给出代码:

#include <bits/stdc++.h>
#define int long long
#define ll __int128
using namespace std;

template <typename T>

void read(T &a)
{
    T x = 0,f = 1;
    char ch = getchar();
    while (ch < '0' || ch > '9')
    {
        if (ch == '-') f = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9')
    {
        x = (x << 1) + (x << 3) + (ch ^ 48);
        ch = getchar();
    }
    a = x * f;
}

template <typename T>

void write(T x)
{
    if (x < 0) putchar('-'),x = -x;
    if (x < 10) return (void) putchar(x + '0');
    write(x / 10);
    putchar(x % 10 + '0');
}

template <typename T>

void writeln(T x)
{
    write(x);
    putchar('\n');
}

template <typename T>

void writes(T x)
{
    write(x);
    putchar(' ');
}

template <typename T,typename... Args> 

void read(T &t,Args &... args)
{
    read(t);
    read(args...);
}

template <typename T,typename... Args>

void writes(T t,Args ... args)
{
    writes(t);
    writes(args...);
}

template <typename T,typename... Args>

void writeln(T t,Args ... args)
{
    writes(t);
    writes(args...);
    putchar('\n');
}

int n;
struct frac
{
    ll p,q,x,y,z;
};

signed main()
{
    int t;
    read(t);
    while (t--)
    {
        read(n);
        std::vector<frac> sk;
        sk.push_back((frac) {1,0,0,0,0});
        sk.push_back((frac) {1,1,0,0,0});
        ll tmp = sqrt(n + 0.5);
        ll ans = 0,p = n / tmp,q = tmp + 1;
        while ("QAQ")
        {
            frac qaq = sk.back();
            while ((p + qaq.p) * (q - qaq.q) > n)
            {
                ans += qaq.q * (2 * q - qaq.q - 1) / 2 * p + (p + 1) * p / 2 * qaq.q;
                ans += (qaq.p + p + q - qaq.q) * qaq.x + (qaq.x + qaq.y) / 2 - qaq.z + (q - qaq.q) * (qaq.p - 1) + (2 * p + qaq.p) * (qaq.p - 1) / 2;
                p += qaq.p;
                q -= qaq.q;
            }
            if (q <= (int) cbrt(n + 0.5)) break;
            frac qwq = sk.back();
            while (!sk.empty() && (p + qwq.p) * (q - qwq.q) <= n) qaq = sk.back(),sk.pop_back(),qwq = sk.back();
            while ("QwQ")
            {
                frac ovo = (frac) {qwq.p + qaq.p,
                    qwq.q + qaq.q,
                    qwq.x + qaq.x + qwq.q * qaq.p,
                    qwq.y + qaq.y + 2 * qwq.q * qaq.x + qwq.q * qwq.q * qaq.p,
                    qwq.z + qaq.z + qwq.p * qaq.x + (qaq.p - 1) * qaq.p / 2 * qwq.q + qwq.p * qwq.q * qaq.p};
                if ((p + ovo.p) * (q - ovo.q) > n) sk.push_back(ovo),qwq = ovo;
                else
                {
                    if ((p + ovo.p) * (p + ovo.p) * qwq.q >= n * qwq.p) break;
                    qaq = ovo;
                }
            }
        }
        for (int i = 1; i < q; i++)
            ans += (ll) n / i * i + (ll) (n / i + 1) * (n / i) / 2;
        writeln(ans - (ll) tmp * tmp * (tmp + 1) / 2 - (ll) n * (n + 1) / 2);
    }
    return 0;
}