码迷,mamicode.com
首页 > 其他好文 > 详细

[spojDIVCNT1]Counting Divisors

时间:2021-02-19 13:39:21      阅读:0      评论:0      收藏:0      [点我收藏+]

标签:col   直接   scan   树根   枚举   lld   时间   本质   使用   

定义

约定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)$的

代码实现

技术图片
 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 }
View Code

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

优化

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

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

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

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

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

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

技术图片
 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 }
View Code

这个过程其实与前面的过程是一致的,实际上可以看作是一个将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)$,且没有浮点误差、常数更小

[spojDIVCNT1]Counting Divisors

标签:col   直接   scan   树根   枚举   lld   时间   本质   使用   

原文地址:https://www.cnblogs.com/PYWBKTDA/p/14411352.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!