传送门
看到出题人是min_25,先orz为敬。
题解:
首先答案可能会爆long long,需要用int128,但是windows不能用,这个就看个人造化了。
我看的题解是这篇:(https://yhx-12243.github.io/OI-transit/records/spojDIVCNT1.html),由于代码本身很短,所以看完题解后自己写的时候,我基本上已经把代码背下来了。。。相似度可能有点高,没办法了。。。
首先我们知道:
∑ i = 1 n σ 0 ( i ) = ∑ i = 1 n ⌊ n i ⌋ sum_{i=1}^nsigma_0(i)=sum_{i=1}^nlfloorfrac{n}{i}rfloor i=1∑nσ0(i)=i=1∑n⌊in⌋
于是就可开心地整除分块TLE了。
直接做显然是不科学的,考虑第一步转化,考虑 n n n以内的每个数的每对因子拆分都有至少一个在 n sqrt n n以内:
∑ i = 1 n σ 0 ( i ) = ∑ i = 1 n ∑ j = 1 n [ i ⋅ j ≤ n ] = ∑ i = 1 n ⌊ n i ⌋ + ∑ j = 1 n ⌊ n j ⌋ − ∑ i = 1 n ∑ j = 1 n 1 = 2 ∑ i = 1 n ⌊ n i ⌋ − ⌊ n ⌋ 2 begin{aligned} sum_{i=1}^nsigma_0(i)&=&&sum_{i=1}^{n}sum_{j=1}^n[icdot jleq n]\ &=&&sum_{i=1}^{sqrt n}lfloorfrac{n}{i}rfloor+sum_{j=1}^{sqrt n}lfloorfrac{n}{j}rfloor-sum_{i=1}^{sqrt n}sum_{j=1}^{sqrt n}1\ &=&&2sum_{i=1}^{sqrt n}lfloorfrac{n}{i}rfloor-lfloorsqrt nrfloor^2 end{aligned} i=1∑nσ0(i)===i=1∑nj=1∑n[i⋅j≤n]i=1∑n⌊in⌋+j=1∑n⌊jn⌋−i=1∑nj=1∑n12i=1∑n⌊in⌋−⌊n⌋2
然而到这一步仍然没有什么屁用,如果你会证明整除分块复杂度是严格 Θ ( n ) Theta(sqrt n) Θ(n)而不是上界 O ( n ) O(sqrt n) O(n)的话你就知道, n sqrt n n以下的每一个数整除 n n n得到的结果都不同,所以整除分块算上面的东西复杂度还是 Θ ( n ) Theta(sqrt n) Θ(n)的,过不了。
现在需要优化求这个东西:
∑
i
=
1
n
⌊
n
i
⌋
sum_{i=1}^{sqrt n}lfloorfrac{n}{i}rfloor
i=1∑n⌊in⌋
仔细考虑一下发现是 y = n x , x = 0 , y = 0 , y = n y=frac{n}{x},x=0,y=0,y=sqrt n y=xn,x=0,y=0,y=n这几条线围成封闭图形内部整点个数,算上边界,不算 x x x轴。
(盗一张图):
那么我们要求的就是紫色区域(以下称为R区域)的整点个数,我们发现在 y = n y=sqrt n y=n以下的图形中,R图形的补形是一个凸壳,考虑去把图形切割成一条条横条拟合这个凸壳,那么我们就需要对这些横条求整点个数。
然后我们有一种相当神奇的做法,确定一个初始点,用Stern-Brocot tree来寻找斜率对点进行平移拟合凸壳。(不知道这是什么的请自行维基)
在拟合过程中,我们必须始终保证点不会跑到紫色区域的边界上,所以初始点最好设置为 ( ⌊ n ⌊ n ⌋ ⌋ , ⌊ n ⌋ + 1 ) (lfloorfrac{n}{lfloorsqrt nrfloor}rfloor,lfloorsqrt nrfloor+1) (⌊⌊n⌋n⌋,⌊n⌋+1)
接下来就是神仙操作,找斜率。
很显然地,斜率始终为负,我们考虑维护斜率绝对值,在初始点的时候需要平移的斜率不会小于-1,并且斜率单调不减。我们用一个单调栈来维护斜率,栈中的斜率全部用SBT上的结点权值来表示,并且斜率绝对值从栈底向栈顶单调增。
初始栈中有两个元素: ( 1 , 0 ) , ( 1 , 1 ) (1,0),(1,1) (1,0),(1,1),每次从栈顶开始取向量,如果加上当前向量加上它不会越界,将它关于 x x x轴的对称向量与当前点相加,由于是SBT上的元素,一定是既约分数,我们可以直接用pick定理算出这一次扫过的横条中的整点个数。
接下来需要对栈进行调整,直接弹栈顶,直到上一个栈顶元素会让当前点跑进R区域且当前栈顶元素不会,将这两个向量记作 L , R L,R L,R。
然后我们需要找到一个更加精确的斜率进行接下来的拟合,直接在SBT上向下走。
设 M = L + R M=L+R M=L+R,显然 M M M是 L , R L,R L,R中某一个在SBT上的儿子。
如果当前的儿子结点不会导致越界就直接 R = M R=M R=M,并把 M M M推入栈中,继续迭代。
如果当前 M M M会导致越界,考虑两种情况,假设当前点横坐标为 x 0 x_0 x0,我们看 f ( x ) = n x f(x)=frac{n}{x} f(x)=xn在 x = x 0 + M x x=x_0+M_x x=x0+Mx处的斜率和 R R R的斜率
- ∣ f ′ ( x 0 + M x x ) ∣ ≤ s l o p e ( R ) |f'(x_0+M_xx)|leq slope(R) ∣f′(x0+Mxx)∣≤slope(R),由 f ′ ( x ) = − n x 2 f'(x) = - dfrac n {x^2} f′(x)=−x2n 可得该条件等价于 n ⋅ R x ≤ ( x 0 + M x ) 2 R y ncdot R_x leq (x_0+M_x)^2 R_y n⋅Rx≤(x0+Mx)2Ry,则之后再怎么迭代都已经凉凉了,直接break。
- 否则,我们令 L = M L=M L=M继续迭代。
复杂度?有篇论文证明了有效斜率只有 O ( n 1 3 ) O(n^{frac{1}{3}}) O(n31)种,所以总复杂度不超过 O ( n 1 3 log n ) O(n^{frac{1}{3}}log n) O(n31logn)
不过在18年论文集里面朱老大给出了一个简洁得多的 O ( n 1 3 log 3 n ) O(n^{frac{1}{3}}log^3 n) O(n31log3n)的证明,也可以参考一下。
代码:
#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;
}
最后
以上就是幽默美女最近收集整理的关于【SPOJ-DIVCNT1】Counting Divisors(Stern-Brocot tree)(凸壳拟合曲线)(pick定理)(数论)传送门题解:的全部内容,更多相关【SPOJ-DIVCNT1】Counting内容请搜索靠谱客的其他文章。
发表评论 取消回复