【SPOJ-DIVCNT1】Counting Divisors（Stern-Brocot tree）（凸壳拟合曲线）（pick定理）（数论）

传送门

---

看到出题人是min_25，先orz为敬。

题解：

首先答案可能会爆long long，需要用int128，但是windows不能用，这个就看个人造化了。

我看的题解是这篇：(https://yhx-12243.github.io/OI-transit/records/spojDIVCNT1.html)，由于代码本身很短，所以看完题解后自己写的时候，我基本上已经把代码背下来了。。。相似度可能有点高，没办法了。。。

首先我们知道：

$$\sum _{i=1}^n{\sigma }_0(i)=\sum _{i=1}^n\lfloor \frac{n}{i}\rfloor$$

于是就可开心地整除分块TLE了。

直接做显然是不科学的，考虑第一步转化，考虑$n$以内的每个数的每对因子拆分都有至少一个在$\sqrt{n}$以内：

$$\begin{array}{rlrl}\sum _{i=1}^n{\sigma }_0(i)& =& & \sum _{i=1}^n\sum _{j=1}^n[i\cdot j\le n]\\ & =& & \sum _{i=1}^{\sqrt{n}}\lfloor \frac{n}{i}\rfloor +\sum _{j=1}^{\sqrt{n}}\lfloor \frac{n}{j}\rfloor -\sum _{i=1}^{\sqrt{n}}\sum _{j=1}^{\sqrt{n}}1\\ & =& & 2\sum _{i=1}^{\sqrt{n}}\lfloor \frac{n}{i}\rfloor -\lfloor \sqrt{n}{\rfloor }^2\end{array}$$

然而到这一步仍然没有什么屁用，如果你会证明整除分块复杂度是严格$\Theta (\sqrt{n})$而不是上界$O(\sqrt{n})$的话你就知道，$\sqrt{n}$以下的每一个数整除$n$得到的结果都不同，所以整除分块算上面的东西复杂度还是$\Theta (\sqrt{n})$的，过不了。

现在需要优化求这个东西：
$$\sum _{i=1}^{\sqrt{n}}\lfloor \frac{n}{i}\rfloor$$

仔细考虑一下发现是$y=\frac{n}{x}，x=0，y=0，y=\sqrt{n}$这几条线围成封闭图形内部整点个数，算上边界，不算$x$轴。

（盗一张图）：

那么我们要求的就是紫色区域（以下称为R区域）的整点个数，我们发现在$y=\sqrt{n}$以下的图形中，R图形的补形是一个凸壳，考虑去把图形切割成一条条横条拟合这个凸壳，那么我们就需要对这些横条求整点个数。

然后我们有一种相当神奇的做法，确定一个初始点，用Stern-Brocot tree来寻找斜率对点进行平移拟合凸壳。（不知道这是什么的请自行维基）

在拟合过程中，我们必须始终保证点不会跑到紫色区域的边界上，所以初始点最好设置为$(\lfloor \frac{n}{\lfloor \sqrt{n}\rfloor }\rfloor ,\lfloor \sqrt{n}\rfloor +1)$

接下来就是神仙操作，找斜率。

很显然地，斜率始终为负，我们考虑维护斜率绝对值，在初始点的时候需要平移的斜率不会小于-1，并且斜率单调不减。我们用一个单调栈来维护斜率，栈中的斜率全部用SBT上的结点权值来表示，并且斜率绝对值从栈底向栈顶单调增。

初始栈中有两个元素：$(1,0),(1,1)$，每次从栈顶开始取向量，如果加上当前向量加上它不会越界，将它关于$x$轴的对称向量与当前点相加，由于是SBT上的元素，一定是既约分数，我们可以直接用pick定理算出这一次扫过的横条中的整点个数。

接下来需要对栈进行调整，直接弹栈顶，直到上一个栈顶元素会让当前点跑进R区域且当前栈顶元素不会，将这两个向量记作$L,R$。

然后我们需要找到一个更加精确的斜率进行接下来的拟合，直接在SBT上向下走。

设$M=L+R$，显然$M$是$L,R$中某一个在SBT上的儿子。

如果当前的儿子结点不会导致越界就直接$R=M$，并把$M$推入栈中，继续迭代。

如果当前$M$会导致越界，考虑两种情况，假设当前点横坐标为$x_0$，我们看$f(x)=\frac{n}{x}$在$x=x_0+M_x$处的斜率和$R$的斜率

1. $|f^{\prime }(x_0+M_{x}x)|\le slope(R)$，由 $f^{\prime }(x)=-\frac{n}{x^2}$ 可得该条件等价于 $n\cdot R_x\le (x_0+M_x{)}^2{R}_y$，则之后再怎么迭代都已经凉凉了，直接break。
2. 否则，我们令$L=M$继续迭代。

复杂度？有篇论文证明了有效斜率只有$O(n^{\frac{1}{3}})$种，所以总复杂度不超过$O(n^{\frac{1}{3}}\log n)$

不过在18年论文集里面朱老大给出了一个简洁得多的$O(n^{\frac{1}{3}}{\log }^{3}n)$的证明，也可以参考一下。

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

#ifdef zxyoi
#define lll long long
#else 
#define lll __int128
#endif

using std::cerr;

inline void print(lll x){
	static char buf[31],len;
	do{buf[++len]=x%10,x/=10;}while(x);
	while(len)putchar(buf[len--]^48);
}

struct Point{
	ll x,y;
	Point(){}
	Point(ll _x,ll _y):x(_x),y(_y){}
	friend Point operator+(cs Point &a,cs Point &b){return Point(a.x+b.x,a.y+b.y);}
};

ll n;

inline bool in(ll x,ll y){return x*y>n;}
inline bool judge(ll x,cs Point &p){return (lll)n*p.x<=(lll)x*x*p.y;}

Point st[100007];int tp;

inline lll solve(){
	ll sq3=cbrt(n),sq=sqrt(n),x=n/sq,y=sq+1;
	lll res=0;Point L,R,M;
	st[1]=Point(1,0),st[tp=2]=Point(1,1);
	while(true){
		for(L=st[tp--];in(x+L.x,y-L.y);x+=L.x,y-=L.y)
		res+=x*L.y+(L.y+1)*(L.x-1)/2;
		if(y<=sq3)break;
		for(R=st[tp];!in(x+R.x,y-R.y);R=st[--tp])L=R;
		for(M=L+R;true;M=L+R)
		if(in(x+M.x,y-M.y))st[++tp]=R=M;
		else {
			if(judge(x+M.x,R))break;
			L=M;
		}
	}
	for(int re i=1;i<y;++i)res+=n/i;
	return res*2-sq*sq;
}

signed main(){
#ifdef zxyoi
	freopen("DIVCNT1.in","r",stdin);
#endif
	int T;scanf("%d",&T);
	while(T--){scanf("%lld",&n);print(solve()),putchar('\n');}
	return 0;
}