[spojDIVCNT1]Counting Divisors

定义

约定1：以下分数都是最简，且令$frac{1}{0}$有意义，其大于其余分数，并称平行于$y$轴的直线斜率为$-frac{1}{0}$

分数加：对于分数$a=frac{a_{1}}{a_{0}}$、$b=frac{b_{1}}{b_{0}}$，定义$aoplus 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$层（$hge 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}lfloorfrac{N}{i} floor$，这个用数论分块做可以做到$o(Tsqrt{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=aoplus 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)$即满足$uvec{i}+vvec{j}+P=(x,y)$

由于是恰好是两个方程两个未知数，且显然不矛盾或等价，因此两坐标系中的点一一对应

接下来，我们还要证明两坐标系中的整点一一对应，换言之若$(u,v)$是整点则$(x,y)$是整点，若$(x,y)$是整点则$(u,v)$也是整点

设$P(x_{0},y_{0})$，由于$P$是整点，$x_{0},y_{0}$也是整数，那么前者显然成立

后者考虑暴力解出方程组，即
$$
egin{cases}ua_{0}-vb_{0}+x_{0}=x\-ua_{1}+vb_{1}+y_{0}=yend{cases}iffegin{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' floor(lfloor b' floor-1)}{2}$

具体式子：$C(sqrt{frac{N}{c}},sqrt{Nc})$，$b=ccdot x_{C}+y_{C}=2y_{C}$，$b'=(a_{0}+b_{0})(b-y_{0})-(a_{1}+b_{1})x_{0}$

问题2

选择$P'(lfloor b' floor,0)$和$P'(0,lfloor b' floor)$即可，正确性显然

边界条件需要考虑曲线$L$在上面所构造的新坐标系中的一些性质，首先就是对应的方程：
$$
(ua_{0}-vb_{0}+x_{0})(-ua_{1}+vb_{1}+y_{0})=N
$$
构造$w_{0}$和$w_{1}$，使得
$$
egin{cases}x_{0}=a_{0}w_{0}-b_{0}w_{1}\y_{0}=b_{1}w_{1}-a_{1}w_{0}end{cases}iffegin{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$，即
$$
egin{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})pmsqrt{(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}}$，暴力枚举$yin [1,lfloor y_{B} floor]$进行统计

问题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 e frac{0}{1}$且$b e 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}-1le x_{P_{2}}$，否则必然存在整点，再根据$x_{A}=sqrt{frac{N}{a}}$以及$x_{B}=sqrt{frac{N}{b}}$，代入即
$$
sqrt{frac{N}{b}}-1le x_{0}-b_{0}<x_{0}+a_{0}le sqrt{frac{N}{a}}+1
$$
由于不确定$x_{0}$，不能直接判断，但可以推出$sqrt{frac{N}{a}}-sqrt{frac{N}{b}}+2ge a_{0}+b_{0}$，状态数即为

$$
sum_{a_{0}b_{1}-a_{1}b_{0}=1}[sqrt{frac{N}{a}}-sqrt{frac{N}{b}}+2ge 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}} floor}sigma(i)cdot sigma(i+1)
$$
又注意到有$xyle frac{x^{2}+y^{2}}{2}$，即
$$
sum_{i=1}^{lfloor N^{frac{1}{3}} floor}sigma(i)cdot sigma(i+1)le sum_{i=1}^{lfloor N^{frac{1}{3}} floor}sigma(i)^{2}
$$
后面的这个东西有结论，$sum_{i=1}^{n}sigma(i)^{2}$是$o(nlog^{3}n)$级别的，因此总复杂度即$o(N^{frac{1}{3}}log^{3}N)$

似乎还可以更精确地证明其是$o(N^{frac{1}{3}}log N)$的

代码实现

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define lll __int128
#define eps 1e-7
int t,a[105];
ll n;
lll ans;
bool check(lll k){
    return k>(lll)1<<63;
}
lll sqr(lll k){
    return k*k;
}
void write(lll ans){
    a[0]=0;
    while (ans){
        a[++a[0]]=ans%10;
        ans/=10;
    }
    if (!a[0])printf("0");
    for(int i=a[0];i;i--)printf("%d",a[i]);
    printf("
");
}
void dfs(lll x0,lll y0,int a1,int a0,int b1,int b0){
    if ((!a1)&&(a0==1)){
        if ((lll)b0*b0*b0>n){
            ll b=floor(2*sqrt(n)*sqrt(b0)+eps);
            for(int i=1;(lll)i*i*b0<=n;i++)ans+=n/i-(b-1LL*i*b0);
            return;
        }
    }
    else{
        ll s=(ll)a0*b1+(ll)a1*b0;
        lll w0=b1*x0+b0*y0+1,w1=a1*x0+a0*y0+1;
        if (check(w0))return;
        lll D=sqr(w0)-(lll)4*b0*b1*n;
        if ((D<0)||(D>sqr(w0)/sqr(s))||(s*w0<(lll)2*b0*b1*w1))return;
        if ((!check(s*w0-(lll)2*b0*b1*w1))&&(sqr(s*w0-(lll)2*b0*b1*w1)<D))return;
    }
    int c1=a1+b1,c0=a0+b0;
    double c=1.0*c1/c0,b=2*sqrt(n)*sqrt(c);
    lll bb=floor((a0+b0)*b+eps)-(a0+b0)*y0-(a1+b1)*x0;
    ans+=bb*(bb-1)/2;
    dfs(bb*a0+x0,-bb*a1+y0,a1,a0,c1,c0);
    if ((b1==1)&&(b0==0))ans=ans*2-bb*(bb-1)/2;
    else dfs(-bb*b0+x0,bb*b1+y0,c1,c0,b1,b0);
}
int main(){
    scanf("%d",&t);
    while (t--){
        scanf("%lld",&n);
        ans=0;
        dfs(0,0,0,1,1,0);
        write(ans);
    }
}

这份代码具有较大的浮点误差，只能保证在$Nle 10^{14}$的正确性，且递归常数较大

优化

考虑上面这个过程中所有结束状态的$p_{a}$和$p_{b}$，若将其以此连在一起，其在曲线$L$的下方（可以有重合的点），且之间没有任何整点

同时，如果能够确定这条折线，由于所有点都是整点，根据皮克定理，不难确定这个折线下方的整点数，也就是所求的答案

事实上，这也就是上面这个算法的本质，我们不妨考虑直接构造这条直线

另外，由于$L$是下凸的，在$L$上方两点连线一定在$L$的外部，而内部没有这个性质，因此考虑构造在$L$外部的折线

具体实现可以参考代码，代码之后会再说明一下

（代码大概就是抄了一下论文中所给的代码）

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define lll __int128
struct frac{
    int x,y;
    frac operator +(const frac &a)const{
        return frac{x+a.x,y+a.y};
    }
};
int t,a[105];
ll n;
lll x,y,ans;
stack<frac>st;
bool check(lll x,lll y){
    return x*y>n;
}
void write(lll ans){
    a[0]=0;
    while (ans){
        a[++a[0]]=ans%10;
        ans/=10;
    }
    if (!a[0])printf("0");
    for(int i=a[0];i;i--)printf("%d",a[i]);
    printf("
");
}
int main(){
    scanf("%d",&t);
    while (t--){
        scanf("%lld",&n);
        ll K=floor(sqrt(n));
        ans=0; 
        x=n/K,y=n/x+1;
        st.push(frac{1,0});
        st.push(frac{1,1});
        while (!st.empty()){
            frac k=st.top();
            st.pop();
            while (check(x+k.x,y-k.y)){
                ans+=x*k.y+(k.y+1)*(k.x-1)/2;
                x+=k.x,y-=k.y;
            }
            if (y*y*y<=n)break;
            frac o;
            while (1){
                o=k,k=st.top();
                if (check(x+k.x,y-k.y))break;
                st.pop();
            }
            while (1){
                frac s=k+o;
                if (check(x+s.x,y-s.y))st.push(k=s);
                else{
                    if ((x+s.x)*(x+s.x)*k.y>=(lll)n*k.x)break;
                    o=s;
                }
            }
        }
        for(int i=1;i<y;i++)ans+=n/i;
        ans=ans*2-K*K;
        write(ans);
    }
}

这个过程其实与前面的过程是一致的，实际上可以看作是一个将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=lfloorsqrt{N} floor$，因为根据对称性有$sum_{i=1}^{N}lfloorfrac{N}{i} floor=2sum_{i=1}^{K}lfloorfrac{N}{i} floor-K^{2}$

2.做第一个优化可以使得我们只需要判定$yle N^{frac{1}{3}}$的时候（不需要判定$xle N^{frac{1}{3}}$），根据前面这里要暴力处理

这样的时间复杂度与之前的复杂度相同，即为$o(N^{frac{1}{3}}log^{3}N)$，且没有浮点误差、常数更小