题解：SP26073 DIVCNT1 - Counting Divisors

本来是 5.15 写完的，但是重新搞了一下。

修改代码风格，之前的太石了。

---

\[\gdef \so { \implies } \gdef \P { \mathbb{P} } \gdef \id { \mathrm{id} } \gdef \lcm { \mathrm{lcm} } \gdef \or { \mathbin{ \hbox{ \rm or } } } \gdef \xor { \mathbin{ \hbox{ \rm xor } } } \gdef \and { \mathbin{ \hbox{ \rm and } } } \gdef \inv #1{ { #1 ^ {-1} } } \gdef \abs #1{ \left \lvert #1 \right \rvert } \gdef \paren #1{ \left \lparen #1 \right \rparen } \gdef \brack #1{ \left \lbrack #1 \right \rbrack } \gdef \brace #1{ \left \lbrace #1 \right \rbrace } \gdef \ceil #1{ \left \lceil #1 \right \rceil } \gdef \floor #1{ \left \lfloor #1 \right \rfloor } \gdef \Vert #1{ \left \lVert #1 \right \rVert } \gdef \vector #1#2{ \begin{bmatrix} #1 \\ #2 \end{bmatrix} } % \gdef \FirstStirling #1#2{ \begin{bmatrix} #1 \\ #2 \end{bmatrix} } % \gdef \SecondStirling #1#2{ \begin{Bmatrix} #1 \\ #2 \end{Bmatrix} } \gdef \zc #1#2{\lfloor\frac{#1}{#2}\rfloor} \gdef \dzc #1#2{\lfloor\dfrac{#1}{#2}\rfloor} \gdef \eps {\varepsilon} \gdef \ph {\varphi} \gdef \ba {\begin{aligned}} \gdef \ea {\end{aligned}} \gdef \lcm {\mathrm{lcm}} \gdef \d {\cdot} \gdef \etc {\cdots} \gdef \Etc {\cdots\cdots} \gdef \O {\mathcal{O}} \gdef \w {\omega} \gdef \sl {\displaystyle} \gdef \inv {^{-1}} \gdef \med{\operatorname{mediant}} \gdef \mediant#1#2#3#4{\operatorname{mediant}\left(\dfrac{#1}{#2},\dfrac{#3}{#4}\right)} \gdef \fkh #1{\left[#1\right]} \gdef \zkh #1{\left[#1\right]} \gdef \xkh #1{\left(#1\right)} \gdef \doublearrow {\Rightarrow} \gdef \sothat {\Longrightarrow} \]

题目等价于 \(y=\dfrac nx\) 函数下方有多少个在一象限的点。

前置知识：Trick for H(n)。了解了这个解法，我们就能做这题了。

\[\ba\sum_{i=1}^n\sigma_0(i)&=\sum_{i=1}^n\sum_{ab=i}1\\&=\sum_{ab\le n}1\\&=\sum_{i=1}^n\floor{\dfrac ni}\\&=-\floor{\sqrt n}^2+2\sum_{i=1}^{\floor{\sqrt n}}\floor{\frac ni}\ea \]

于是计算 \(\sl \sum_{i=1}^{\floor{\sqrt n}}\floor{\frac ni}\)。

这一部分考虑化曲为直。将双曲线变成一些线段。

抽象理解。用一坨线段拟合这个双曲线。

如图，橙色是用来拟合的线段。

容易发现斜率递减。

首先我们知道对于满足 \(\gcd(u,v)=1\) 的向量 \((u,-v)\) 以及点 \((x,y)\)，可以用皮克定理计算。
:::info[计算过程]{open}

左边那一部分显然是 \(xv\)。

右边不太好计算：

考虑皮克定理：

\[\ba S&=I+\dfrac B2-1\\\sothat I&=S-\dfrac B2+1\ea \]

来计算 \(B\)：

\[\ba B&=B_{AH}+B_{BH}+B_{AB}-3\\&=(v+1)+(u+1)+2-3\\&=u+v+1\ea \]

所以：

\[\ba I&=\dfrac {uv-(u+v+1)+2}2\\&=\dfrac{uv-u-v+1}2\ea \]

\[\ba\Delta ans&=(x-1)v+I+B\red{-2}\\&=(x-1)v+\dfrac{uv-u-v+1}2+u+v+1-2\\&=(x-1)v+\dfrac{uv+u+v+1}2-1\ea \]

注意标红的 \(\red{-2}\)：因为 \(A,B\) 点都是在要求的区域外的。具体原因参考接下来的过程。

:::

剩下的考虑使用 Stern-Brocot Tree 维护任意斜率向量。

二分。首先弄一个斜率单调栈储存向量。

取出栈顶 \(L\)。

首先一直尝试走 \(L\)，直到走进双曲线下方或双曲线上（记为区域 \(X\)）

然后一直弹栈。\(R\) 为栈顶。直到出现当前位置加上向量 \(R\) 不在 \(X\) 内，但是加上 \(L\) 就进入了 \(X\) 内。

判断死循环：

如果 \(R\) 的斜率 \(\dfrac{\Delta y}{\Delta x}\) 小于等于 \(f(x)=n\d x^{-1}\) 在 \(x+x_{mid}\) 处的斜率，则无解。如图所示，二分永远不会结束。

其斜率为 \(f'(x)=-n\d x^{-2}\)。判断条件为 \(\dfrac{\Delta y}{\Delta x}le-n\d x^{-2}\) 即 \(\Delta y\d x^2\le-n\Delta x\)。

#include<bits/stdc++.h>
using namespace std;
#define int long long
struct point {
    int x, y;
    point(int a, int b){x = a, y = b;}
};
struct vec {
    int dx, dy;
    //dy<0， 存储 |dy|
    vec(){dx = dy = 0;}
    vec(int da, int db){dx = da, dy = db;}
    inline friend vec operator + (vec A, vec B){
        return {A.dx + B.dx, A.dy + B.dy};
    }
};
point operator + (point a, vec V){return {a.x + V.dx, a.y - V.dy};}
// point operator+(vec V,point a){return {a.x+V.dx,a.y+V.dy};}
#define out(P) ((__int128)((P).x)*(P).y>n)
vec st[1000001];
int Top = 0;
inline void push(vec x){st[++Top] = x;}
inline void pop(){Top--;}
inline void clear(){Top = 0;}
vec top(){return st[Top];}
#ifdef __linux__
#define getchar getchar_unlocked
#define putchar putchar_unlocked
#else
#define getchar _getchar_nolock
#define putchar _putchar_nolock
#endif
inline int read(){
    int 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();}
    return x * f;
}
inline void write(__int128 x){
    if(x < 0){putchar('-'); x =-x;}
    if(x > 9)write(x / 10);
    putchar(x % 10 + '0');
}
inline bool check(int x, vec A, int n){return(__int128)x * x * (-A.dy) <=-(__int128)n * A.dx;}
__int128 solve(int n){
    clear();
    push({1, 0});
    push({1, 1});
    int A = cbrt(n), B = sqrt(n);
    __int128 ans = 0;
    point now(n / B, B + 1);
    while(1){
        vec L = top();
        pop();
        while(out(now + L)){
            ans += (__int128)(now.x - 1) * L.dy + ((__int128)L.dx * L.dy + L.dx + L.dy + 1) / 2 - 1;
            now = now + L;
        }
        if(now.y <= A)break;
        vec R = top();
        while(!out(now + R)){
            L = R;
            pop();
            R = top();
        }
        while(1){
            vec mid = L + R;
            if(out(now + mid)){
                R = mid;
                push(mid);
            }
            else {
                if(check(now.x + mid.dx, R, n)){//R 是最缓的一个
                    break;
                    //没救了
                }
                else L = mid;
            }
        }
    }
    for (int i = 1; i < now.y; i++)ans += n / i;
    return ans * 2 - (__int128)B * B;
}
#define endl '\n'
signed main(){
    ios::sync_with_stdio(0); cin.tie(0);
    int T = read();
    while(T--){
        write(solve(read()));
        putchar('\n');
    }
    return 0;
}