[spojDIVCNT1]Counting Divisors
定义
约定1:以下分数都是最简,且令$\frac{1}{0}$有意义,其大于其余分数,并称平行于$y$轴的直线斜率为$-\frac{1}{0}$
分数加:对于分数$a=\frac{a_{1}}{a_{0}}$、$b=\frac{b_{1}}{b_{0}}$,定义$a\oplus b=\frac{a_{1}+b_{1}}{a_{0}+b_{0}}$
Farey neighbor:对于两个非负最简分数$a=\frac{a_{1}}{a_{0}}$和$b=\frac{b_{1}}{b_{0}}$,其为Farey neighbor当且仅当$a_{0}b_{1}-a_{1}b_{0}=1$
约定2:二叉树根深度为0
Stern-Brocot Tree:这是一棵满二叉树,其每一个节点的元素是一个分数
其中根的分数为$\frac{1}{1}$,第$h$层($h\ge 1$)的分数通过如下方式构造:
对于前$h-1$层,其中序遍历依次为$a_{1},a_{2},...,a_{2^{h}-1}$,再令$a_{0}=\frac{0}{1}$,$a_{2^{h}}=\frac{1}{0}$,则第$h$层的元素从左到右依次为$a_{0}\oplus a_{1},a_{1}\oplus a_{2},...,a_{2^{h}-1}\oplus a_{2^{h}}$
约定3:$x_{A}$和$y_{A}$分别表示点$A$的$x$和$y$坐标,$A$为整点当且仅当$x_{A},y_{A}\in Z$
约定4:"图形内"不包含边界,严格在某一条线上指不包含该线端点的部分
题解
题目即求$\sum_{i=1}^{N}\lfloor\frac{N}{i}\rfloor$,这个用数论分块做可以做到$o(T\sqrt{N})$,但还是无法通过
考虑这个式子的几何意义:对于$xy=N$在第一象限内的曲线$L$,其与$x$轴和$y$轴所构成的图形内以及严格在$L$上的整点数
定义$S(P,a,b)$:令$A$和$B$分别表示$L$上导数为$-a$和$-b$的两点,$t_{A}$为$L$过点$A$的切线,则$S(P,a,b)$为$t_{A}$、$t_{B}$与$L$所构成的图形内部以及严格在$L$的$AB$上的整点数,并令$T$为$t_{A}$与$t_{B}$的交点
(虽然与$P$无关,但没有关系)
设置$S(P,a,b)$的定义域为:
1.$a$和$b$为Farey neighbor
2.$P$是整点,过点$P$作两条斜率为$-a$和$-b$的直线,记作$p_{a}$和$p_{b}$,称$p_{a}$对应的直线为$t_{A}$,$p_{b}$对应的直线为$t_{B}$,则满足——
(1)$p_{a},p_{b}$在其对应的直线下方(可以重合)
(2)这四条直线所构成的无限大的图形内没有整点
(3)$t_{A},t_{B}$分别满足:与对应其的直线重合或其上没有整点
初始即求$S(P(0,0),\frac{0}{1},\frac{1}{0})$,符合定义域
$S(P,a,b)$的计算
考虑求$S(P,a,b)$,设切点别为$A$和$B$,$a=\frac{a_{1}}{a_{0}}$,$b=\frac{b_{1}}{b_{0}}$,令$c=a\oplus b$,选取点$C$使得$x_{B}<x_{C}<x_{A}$且$f'(x_{C})=-c$,由于$f'(x)$(在$x>0$)是一个单调递增且连续的函数,因此显然存在
设$t_{C}$交$t_{A}$和$t_{B}$分别于$P_{1}$和$P_{2}$,那么所求的问题也就分为了三个部分$S(P',a,c)$、$S(P',c,b)$和在$\Delta TP_{1}P_{2}$内部以及严格在$P_{1}P_{2}$上的整点数,递归即可
但现在还有三个问题:
1.第三部分的计算
2.$P'$的选择以及递归结束条件
3.时间复杂度
问题1
首先,根据$P$的第2和3个条件,可以看作与$P$相交,具体来说再令$P_{1}$和$P_{2}$分别为$t_{C}$与$p_{a}$和$p_{b}$的交点,那么所求的也就是$\Delta PP_{1}P_{2}$以及严格在$P_{1}P_{2}$上的整点数
对于上面的问题,考虑构建一个坐标系:
1.$P$为原点,$PP_{1}$和$PP_{2}$分别为$x$轴和$y$轴(同时方向也按照这个,即$P_{1}$和$P_{2}$都在$x$和$y$轴正半轴上)
2.令$\vec{i}=(a_{0},-a_{1})$,$\vec{j}=(-b_{0},b_{1})$,则单位长度即为两向量模长
不难发现,$x$轴方向上单位长度的向量即为$\vec{i}$,类似地$y$轴上是$j$,那么对于$(x,y)$,其所对应的新坐标系中的点$(u,v)$即满足$u\vec{i}+v\vec{j}+P=(x,y)$
由于是恰好是两个方程两个未知数,且显然不矛盾或等价,因此两坐标系中的点一一对应
接下来,我们还要证明两坐标系中的整点一一对应,换言之若$(u,v)$是整点则$(x,y)$是整点,若$(x,y)$是整点则$(u,v)$也是整点
设$P(x_{0},y_{0})$,由于$P$是整点,$x_{0},y_{0}$也是整数,那么前者显然成立
后者考虑暴力解出方程组,即
$$
\begin{cases}ua_{0}-vb_{0}+x_{0}=x\\-ua_{1}+vb_{1}+y_{0}=y\end{cases}\iff\begin{cases}u=\frac{b_{1}(x-x_{0})+b_{0}(y-y_{0})}{a_{0}b_{1}-a_{1}b_{0}}\\v=\frac{a_{1}(x-x_{0})+a_{0}(y-y_{0})}{a_{0}b_{1}-a_{1}b_{0}}\end{cases}
$$
由于$a$和$b$是Farey neighbor,根据定义即$a_{0}b_{1}-a_{1}b_{0}=1$,因此$u$和$v$也是整数
由于整点一一对应,我们将$\Delta PP_{1}P_{2}$放到这个新的坐标系中去,考虑$P_{1}P_{2}$斜率:
$$
-ua_{1}+vb_{1}+y_{0}=-\frac{a_{1}+b_{1}}{a_{0}+b_{0}}(ua_{0}-vb_{0}+x_{0})+b
$$
简单化简,即
$$
v=\frac{(a_{0}+b_{0})a_{1}-(a_{1}+b_{1})a_{0}}{(a_{0}+b_{0})b_{1}-(a_{1}+b_{1})b_{0}}u+b'
$$
展开后可以发现,这个斜率恰好为-1,同时$P_{1}$与$x$轴有交点且在$P_{1}P_{2}$上,因此在新坐标系中有$P_{1}(b',0)$,类似地$P_{2}(0,b')$,那么不难求出整点数为$\frac{\lfloor b'\rfloor(\lfloor b'\rfloor-1)}{2}$
具体式子:$C(\sqrt{\frac{N}{c}},\sqrt{Nc})$,$b=c\cdot x_{C}+y_{C}=2y_{C}$,$b'=(a_{0}+b_{0})(b-y_{0})-(a_{1}+b_{1})x_{0}$
问题2
选择$P'(\lfloor b'\rfloor,0)$和$P'(0,\lfloor b'\rfloor)$即可,正确性显然
边界条件需要考虑曲线$L$在上面所构造的新坐标系中的一些性质,首先就是对应的方程:
$$
(ua_{0}-vb_{0}+x_{0})(-ua_{1}+vb_{1}+y_{0})=N
$$
构造$w_{0}$和$w_{1}$,使得
$$
\begin{cases}x_{0}=a_{0}w_{0}-b_{0}w_{1}\\y_{0}=b_{1}w_{1}-a_{1}w_{0}\end{cases}\iff\begin{cases}w_{0}=\frac{b_{1}x_{0}+b_{0}y_{0}}{a_{0}b_{1}-a_{1}b_{0}}\\w_{1}=\frac{a_{1}x_{0}+a_{0}y_{0}}{a_{0}b_{1}-a_{1}b_{0}}\end{cases}
$$
由于$a_{0}b_{1}-a_{1}b_{0}=1$,即
$$
\begin{cases}w_{0}=b_{1}x_{0}+b_{0}y_{0}\\w_{1}=a_{1}x_{0}+a_{0}y_{0}\end{cases}
$$
接下来,代入即可得到
$$
(a_{0}(u+w_{0})-b_{0}(v+w_{1}))(b_{1}(v+w_{1})-a_{1}(u+w_{0}))=N
$$
化简后,即
$$
(a_{0}b_{1}+a_{1}b_{0})(u+w_{0})(v+w_{1})-a_{0}a_{1}(u+w_{0})^{2}-b_{0}b_{1}(v+w_{1})^{2}=N
$$
这个形式已经比较简洁了,直接暴力使用求根公式,即
$$
v=\frac{(a_{0}b_{1}+a_{1}b_{0})(u+w_{0})\pm\sqrt{(a_{0}b_{1}+a_{1}b_{0})^{2}(u+w_{0})^{2}-4b_{0}b_{1}(a_{0}a_{1}(u+w_{0})^{2}-N)}}{2b_{0}b_{1}}-w_{1}
$$
对后面的根号中提出$(u+w_{0})^{2}$的系数,恰好是$(a_{0}b_{1}-a_{1}b_{0})^{2}=1$,即
$$
v=\frac{(a_{0}b_{1}+a_{1}b_{0})(u+w_{0})\pm \sqrt{(u+w_{0})^{2}-4b_{0}b_{1}N}}{2b_{0}b_{1}}-w_{1}
$$
对于正负号,联系图像不难发现另一个解是在$L$上$B$的左侧,其的$v$较大,因此应该取负,即
$$
V(u)=\frac{(a_{0}b_{1}+a_{1}b_{0})(u+w_{0})-\sqrt{(u+w_{0})^{2}-4b_{0}b_{1}N}}{2b_{0}b_{1}}-w_{1}
$$
对其求导,即
$$
V'(u)=\frac{a_{0}b_{1}+a_{1}b_{0}-\frac{u+w_{0}}{\sqrt{(u+w_{0})^{2}-4b_{0}b_{1}N}}}{2b_{0}b_{1}}
$$
考虑后式,即
$$
\frac{u+w_{0}}{\sqrt{(u+w_{0})^{2}-4b_{0}b_{1}N}}=\sqrt{1+\frac{4b_{0}b_{1}N}{(u+w_{0})^{2}-4b_{0}b_{1}N}}
$$
因此$V(u)$的导数单调递增,换言之其下凸
进一步的,考虑$A$所对应的点必然是$v$最小的点,根据下凸那么$A$所对应的$u$之前的点都是单调递减的,这也就是我们所关心的范围
这个范围中,由于$x$轴和$y$轴是不算的,因此若$u=1$时$V'(u)>0$或$V(u)<1$,结果即为0
在这个基础上,我们再设置一个边界条件:
若$a=\frac{0}{1}$或$b=\frac{1}{0}$(以$a=\frac{0}{1}$为例),当$b_{0}>N^{\frac{1}{3}}$,暴力枚举$y\in [1,\lfloor y_{B}\rfloor]$进行统计
问题3
下面我们进行复杂度分析,同样按照上面的边界条件分为两类
1.当$a=\frac{0}{1}$或$b=\frac{1}{0}$(以$a=\frac{0}{1}$为例),$b_{0}$是不断增加1的,且由于$b_{0}>N^{\frac{1}{3}}$时即结束,那么状态只有$o(N^{\frac{1}{3}})$个,同时枚举$y$的时间复杂度为$y_{B}=\frac{N}{x_{B}}=\sqrt{bN}=\sqrt{\frac{N}{b_{0}}}$,计算复杂度也是$o(N^{\frac{1}{3}})$
2.当$a\ne \frac{0}{1}$且$b\ne \frac{1}{0}$,将所有状态,再分为两类:
(1)直接返回的状态,假设有$N_{0}$个
(2)递归下去的状态,假设有$N_{2}$个
考虑这棵搜索树,不难证明$N_{0}=N_{2}+1$,即我们的复杂度是$o(N_{0}+N_{2})=o(N_{2})$
对于$N_{2}$,由于每一个$(a,b)$最多被搜索一次,也就是有多少对$(a,b)$是第二类
考虑$(a,b)$是第二类状态的条件,必要条件为新坐标系中$(u,v)=(0,1)或(1,0)$都在新坐标系中$V(u)$与$x$轴、$y$轴内部
这两个点所对应的点的$x$坐标即为$x_{0}-b_{0}$和$x_{0}+a_{0}$,其恰好在$PP_{1}$或$PP_{2}$上,该图形的$x$坐标范围为$[x_{P_{2}},x_{P_{1}}]$,因此即有
$$
x_{P_{2}}\le x_{0}-b_{0}<x_{0}+a_{0}\le x_{P_{1}}
$$
根据$P$的第二个性质,有$x_{P_{1}}\le x_{A}+1$以及$x_{B}-1\le x_{P_{2}}$,否则必然存在整点,再根据$x_{A}=\sqrt{\frac{N}{a}}$以及$x_{B}=\sqrt{\frac{N}{b}}$,代入即
$$
\sqrt{\frac{N}{b}}-1\le x_{0}-b_{0}<x_{0}+a_{0}\le \sqrt{\frac{N}{a}}+1
$$
由于不确定$x_{0}$,不能直接判断,但可以推出$\sqrt{\frac{N}{a}}-\sqrt{\frac{N}{b}}+2\ge a_{0}+b_{0}$,状态数即为
$$
\sum_{a_{0}b_{1}-a_{1}b_{0}=1}[\sqrt{\frac{N}{a}}-\sqrt{\frac{N}{b}}+2\ge a_{0}+b_{0}]
\\=\sum_{a_{0}b_{1}-a_{1}b_{0}=1}[\frac{\sqrt{a_{0}b_{1}}-\sqrt{a_{1}b_{0}}}{\sqrt{a_{1}b_{1}}}\ge \frac{a_{0}+b_{0}-2}{\sqrt{N}}]
$$
注意到$a_{0}b_{1}=a_{1}b_{0}+1$,有$\sqrt{w+1}-\sqrt{w}=\frac{1}{\sqrt{w}+\sqrt{w+1}}$,即
$$
\sum_{a_{0}b_{1}-a_{1}b_{0}=1}[\frac{\sqrt{a_{0}b_{1}}-\sqrt{a_{1}b_{0}}}{\sqrt{a_{1}b_{1}}}\ge \frac{a_{0}+b_{0}-2}{\sqrt{N}}]\\\le \sum_{a_{0}b_{1}-a_{1}b_{0}=1}[\frac{1}{\sqrt{a_{1}^{2}b_{0}b_{1}}}\ge \frac{a_{0}+b_{0}-2}{\sqrt{N}}]\\=\sum_{a_{0}b_{1}-a_{1}b_{0}=1}[a_{1}^{2}b_{0}b_{1}(a_{0}+b_{0}-2)^{2}\le N]
$$
由于$a_{0},b_{0}\ge 1$,因此$(a_{0}+b_{0}-2)^{2}\ge 2a_{0}b_{0}$,即
$$
\sum_{a_{0}b_{1}-a_{1}b_{0}=1}[a_{1}^{2}b_{0}b_{1}(a_{0}+b_{0}-2)^{2}\le N]\\\le \sum_{a_{0}b_{1}-a_{1}b_{0}=1}[2a_{0}a_{1}^{2}b_{0}^{2}\frac{1+a_{1}b_{0}}{a_{0}}\le N]\\\le \sum_{a_{0}b_{1}-a_{1}b_{0}=1}[2a_{1}^{3}b_{0}^{3}\le N]\\\le \sum_{i=1}^{\lfloor N^{\frac{1}{3}}\rfloor}\sigma(i)\cdot \sigma(i+1)
$$
又注意到有$xy\le \frac{x^{2}+y^{2}}{2}$,即
$$
\sum_{i=1}^{\lfloor N^{\frac{1}{3}}\rfloor}\sigma(i)\cdot \sigma(i+1)\le \sum_{i=1}^{\lfloor N^{\frac{1}{3}}\rfloor}\sigma(i)^{2}
$$
后面的这个东西有结论,$\sum_{i=1}^{n}\sigma(i)^{2}$是$o(n\log^{3}n)$级别的,因此总复杂度即$o(N^{\frac{1}{3}}\log^{3}N)$
似乎还可以更精确地证明其是$o(N^{\frac{1}{3}}\log N)$的
代码实现
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define ll long long 4 #define lll __int128 5 #define eps 1e-7 6 int t,a[105]; 7 ll n; 8 lll ans; 9 bool check(lll k){ 10 return k>(lll)1<<63; 11 } 12 lll sqr(lll k){ 13 return k*k; 14 } 15 void write(lll ans){ 16 a[0]=0; 17 while (ans){ 18 a[++a[0]]=ans%10; 19 ans/=10; 20 } 21 if (!a[0])printf("0"); 22 for(int i=a[0];i;i--)printf("%d",a[i]); 23 printf("\n"); 24 } 25 void dfs(lll x0,lll y0,int a1,int a0,int b1,int b0){ 26 if ((!a1)&&(a0==1)){ 27 if ((lll)b0*b0*b0>n){ 28 ll b=floor(2*sqrt(n)*sqrt(b0)+eps); 29 for(int i=1;(lll)i*i*b0<=n;i++)ans+=n/i-(b-1LL*i*b0); 30 return; 31 } 32 } 33 else{ 34 ll s=(ll)a0*b1+(ll)a1*b0; 35 lll w0=b1*x0+b0*y0+1,w1=a1*x0+a0*y0+1; 36 if (check(w0))return; 37 lll D=sqr(w0)-(lll)4*b0*b1*n; 38 if ((D<0)||(D>sqr(w0)/sqr(s))||(s*w0<(lll)2*b0*b1*w1))return; 39 if ((!check(s*w0-(lll)2*b0*b1*w1))&&(sqr(s*w0-(lll)2*b0*b1*w1)<D))return; 40 } 41 int c1=a1+b1,c0=a0+b0; 42 double c=1.0*c1/c0,b=2*sqrt(n)*sqrt(c); 43 lll bb=floor((a0+b0)*b+eps)-(a0+b0)*y0-(a1+b1)*x0; 44 ans+=bb*(bb-1)/2; 45 dfs(bb*a0+x0,-bb*a1+y0,a1,a0,c1,c0); 46 if ((b1==1)&&(b0==0))ans=ans*2-bb*(bb-1)/2; 47 else dfs(-bb*b0+x0,bb*b1+y0,c1,c0,b1,b0); 48 } 49 int main(){ 50 scanf("%d",&t); 51 while (t--){ 52 scanf("%lld",&n); 53 ans=0; 54 dfs(0,0,0,1,1,0); 55 write(ans); 56 } 57 }
这份代码具有较大的浮点误差,只能保证在$N\le 10^{14}$的正确性,且递归常数较大
优化
考虑上面这个过程中所有结束状态的$p_{a}$和$p_{b}$,若将其以此连在一起,其在曲线$L$的下方(可以有重合的点),且之间没有任何整点
同时,如果能够确定这条折线,由于所有点都是整点,根据皮克定理,不难确定这个折线下方的整点数,也就是所求的答案
事实上,这也就是上面这个算法的本质,我们不妨考虑直接构造这条直线
另外,由于$L$是下凸的,在$L$上方两点连线一定在$L$的外部,而内部没有这个性质,因此考虑构造在$L$外部的折线
具体实现可以参考代码,代码之后会再说明一下
(代码大概就是抄了一下论文中所给的代码)
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define ll long long 4 #define lll __int128 5 struct frac{ 6 int x,y; 7 frac operator +(const frac &a)const{ 8 return frac{x+a.x,y+a.y}; 9 } 10 }; 11 int t,a[105]; 12 ll n; 13 lll x,y,ans; 14 stack<frac>st; 15 bool check(lll x,lll y){ 16 return x*y>n; 17 } 18 void write(lll ans){ 19 a[0]=0; 20 while (ans){ 21 a[++a[0]]=ans%10; 22 ans/=10; 23 } 24 if (!a[0])printf("0"); 25 for(int i=a[0];i;i--)printf("%d",a[i]); 26 printf("\n"); 27 } 28 int main(){ 29 scanf("%d",&t); 30 while (t--){ 31 scanf("%lld",&n); 32 ll K=floor(sqrt(n)); 33 ans=0; 34 x=n/K,y=n/x+1; 35 st.push(frac{1,0}); 36 st.push(frac{1,1}); 37 while (!st.empty()){ 38 frac k=st.top(); 39 st.pop(); 40 while (check(x+k.x,y-k.y)){ 41 ans+=x*k.y+(k.y+1)*(k.x-1)/2; 42 x+=k.x,y-=k.y; 43 } 44 if (y*y*y<=n)break; 45 frac o; 46 while (1){ 47 o=k,k=st.top(); 48 if (check(x+k.x,y-k.y))break; 49 st.pop(); 50 } 51 while (1){ 52 frac s=k+o; 53 if (check(x+s.x,y-s.y))st.push(k=s); 54 else{ 55 if ((x+s.x)*(x+s.x)*k.y>=(lll)n*k.x)break; 56 o=s; 57 } 58 } 59 } 60 for(int i=1;i<y;i++)ans+=n/i; 61 ans=ans*2-K*K; 62 write(ans); 63 } 64 }
这个过程其实与前面的过程是一致的,实际上可以看作是一个将dfs的栈手动实现的过程
从dfs的角度考虑,栈中相邻两个元素即为之前的一次dfs中的两个分数
之后相当于不断递归,其中有一个比较特殊的地方就是关于
if ((x+s.x)*(x+s.x)*k.y>=(lll)n*k.x)break;
这里判定的其实是$-f'(x+s.x)\ge \frac{k.y}{k.x}$,这个判定是类似于前面判定$U(1)\ge 1$的,也就是说没有整点就结束
关于没有整点,其实就等价于没有斜率(绝对值)比我更大的直线,否则那一个点就在其下方且是整点,因此要求$-f'(x+s.x)\ge \frac{k.y}{k.x}$
另外还有一些细节问题:
1.从$K=\lfloor\sqrt{N}\rfloor$,因为根据对称性有$\sum_{i=1}^{N}\lfloor\frac{N}{i}\rfloor=2\sum_{i=1}^{K}\lfloor\frac{N}{i}\rfloor-K^{2}$
2.做第一个优化可以使得我们只需要判定$y\le N^{\frac{1}{3}}$的时候(不需要判定$x\le N^{\frac{1}{3}}$),根据前面这里要暴力处理
这样的时间复杂度与之前的复杂度相同,即为$o(N^{\frac{1}{3}}\log^{3}N)$,且没有浮点误差、常数更小