题目描述
定义 (d(n)) 为 (n) 的正因数的个数,比如 (d(2) = 2, d(6) = 4)。
令 $ S_1(n) = sum_{i=1}^n d(i) $
给定 (n),求 (S_1(n))。
输入格式
第一行包含一个正整数 (T) ((T leq 10^5)),表示数据组数。
接下来的 (T) 行,每行包含一个正整数 (n) ((n < 2^{63}))。
输出格式
对于每个 (n),输出一行一个整数,表示 (S_1(n)) 的值。
题解
显然 $ S_1(n) = sum_{i=1}^n left lfloor frac ni
ight
floor (
整除分块可以做到)O left( sqrt n
ight)(,这是要死的复杂度。
但实际上,这玩意是有个优化套路的:
显然有) sum_{i=1}^n left lfloor frac ni
ight
floor = 2 sum_{i=1}^{lfloor sqrt n
floor} left lfloor frac ni
ight
floor - left lfloor sqrt n
ight
floor^2 $
那么只需求 (sumlimits_{i=1}^{lfloor sqrt n
floor} left lfloor dfrac ni
ight
floor) 的值
它显然等于 (f(x) = dfrac n x)、直线 (y = sqrt n)、(x) 轴和 (y) 轴所围成的区域 (R) 中整点 (格点) 的个数 (含边界,(x) 轴除外)。我们把它变成了一个区间数点问题。
但实际上这个函数是个凸函数,因此区域 (R) 的补集 (R') 为一个凸集,我们很显然的想到用凸包去拟合,然后Pick定理随便搞一搞。
怎么找凸包呢?我们用SBT大力求斜率,然后弄个单调栈去搞。
记 (B = left lfloor sqrt n
ight
floor),我们从点 (P_0 left( left lfloor dfrac nB
ight
floor, B+1
ight)) 开始,一步一步寻找凸包上所有的点。
(步骤一)我们在栈中加入两个分数 (为斜率的绝对值) (dfrac 01) 和 (dfrac 11),代表向量 ((1, 0)) 和 ((1, 1)),由于 (f' left( sqrt n
ight) = -1),因此这部分所有斜率都在 (-1 sim 0) 之间。 此外,这个栈需要满足自顶向上是单调递增的。
(步骤二)接下来,我们取出栈顶向量,将 (P_0) 持续与这个向量 (关于 (x) 轴的对称,下略) 相加,直到这个点 (P_k) 不在区域 (R') 中 (即在区域 (R) 中)。由于这些分数是在 Stern-Brocot 树中产生的,因此一定是既约分数 (即分子与分母互素)。
因此每加一步,我们可以计算出这个横条的面积:设向量为 ((u, v)),上一个点 ((P_{k-1})) 的横坐标为 (x),则面积为 (x v + dfrac 12 (v + 1) (u - 1))。
(步骤三)接着我们要对栈进行调整。由于函数的斜率的绝对值单调递减,因此栈中的分数也需要单调递减。
故我们需要把值过大的分数弹出栈外,直到栈顶和它下面的元素与 (P_k) 相加后,前者在 (R') 外,后者在 (R') 内。把这两个向量记作 (l) 和 (r)。
(步骤四)然后我们在SBT上二分(不是说大力求嘛).
记 (l = dfrac {y_l} {x_l}, r = dfrac {y_r} {x_r}) ((l > r)),则 (m = dfrac {y_m} {x_m} = dfrac {y_l + y_r} {x_l + x_r})。令 (M = P_k + m) (即 (P_k) 按照向量 (m) 平移后的点),如果 (M) 在 (R') 中,则将 (r) 压入栈后令 (r gets m),继续二分;否则,分以下两种情况讨论(二轮判断):
一:如果 (left| f' left( x_k + x_m
ight)
ight| leq r) (其中 (x_k) 为 (P_k) 的横坐标)——由 (f'(x) = - dfrac n {x^2}) 可得该条件等价于 (n x_r leq (x_k + x_m)^2 y_r) ——则容易得到如果再迭代下去的话,所有的 (P_k + m) 都不会落在 (R') 内,也就能说已经二分完毕了,因此只需保留当前的栈重新回到步骤 2 进入下一轮迭代。
二:如果 (left| f' left( x_k + x_m
ight)
ight| > r),则接下来的向量还有可能落入 (R') 中,因此令 (l gets m) 后继续二分。
这个做法的复杂度上界不超过 (O left( n^{1/3} log n
ight)) ,因为斜率只有不超过(O left( n^{1/3}
ight))个。(至于怎么证明文文也不会)。
后期函数的斜率非常小,都是 (dfrac 1k) 的形式,会退化为 (O(y)),因此可以在 (y leq sqrt[3]n) 的部分中直接暴力计算
Code
#include <bits/stdc++.h>
const int N=1000005;
typedef long long ll;
struct qwq {
ll x, y;
qwq (ll x0 = 0, ll y0 = 0) : x(x0), y(y0) {}
inline qwq operator + (const qwq &B) const {return qwq(x + B.x, y + B.y);}
};
ll n;
int top;
qwq stk[N];
inline void push(qwq x) { stk[++top]=x;}
inline void putint(__int128 x) {
static char buf[36];
if (!x) {putchar(48); return;} int i = 0;
for (; x; buf[++i] = x % 10 | 48, x /= 10);
for (; i; --i) putchar(buf[i]);
}
inline bool check(ll x, ll y) {return n < x * y;}
inline bool judge(ll x, qwq v) {return (__int128)n * v.x <= (__int128)x * x * v.y;}
__int128 S1() {
int i, crn = cbrt(n);
ll srn = sqrt(n), x = n / srn, y = srn + 1;
__int128 ret = 0;
qwq L, R, M;
push(qwq(1, 0)); push(qwq(1, 1));
while(true) {
for (L = stk[top--]; check(x + L.x, y - L.y); x += L.x, y -= L.y)
ret += x * L.y + (L.y + 1) * (L.x - 1) / 2;
if (y <= crn) break;
for (R = stk[top]; !check(x + R.x, y - R.y); R = stk[--top]) L = R;
for (; M = L + R, 1; )
if (check(x + M.x, y - M.y)) push(R = M);
else {
if (judge(x + M.x, R)) break;
L = M;
}
}
for (i = 1; i < y; ++i) ret += n / i;
return ret*2-srn*srn;
}
int main() {
int T;
for (scanf("%d", &T); T; --T) {scanf("%lld", &n); putint(S1());puts("");}
return 0;
}