[BZOJ 1041] 圆上的整点

本博客中部分文章已被迁移至新博客，但本文由于内容较平凡并未迁移且不在迁移计划中。

---

Description
　　求一个给定的圆($x^2+y^2=r^2$)，在圆周上有多少个点的坐标是整数
Input
　　只有一个正整数$n$,$n\le 2000000000$
Output
　　整点个数
Sample Input
4
Sample Output
4

初见这道题时候果断不会。无奈上网找题解。
我在网上找到的做法是这样的：http://hzwer.com/1457.html
超级麻烦的推导过程让没有耐心的我直接放弃这道题。

后来胡乱翻数论书，偶然看见一个神奇的定理，它大概长这样：

设$n\ge 1$，则不定方程$x^2+y^2=n$的解数为：
$$N(p)=4\sum _{d|p}h(d)$$其中函数$h(d)$定义为：
KaTeX parse error: No such environment: eqnarray at position 8: \begin{̲e̲q̲n̲a̲r̲r̲a̲y̲}̲h(d)= \begin{ca…

至于这个定理的证明，中心思想就是把不定方程 $x^2+y^2=n$ 的非负本原解和二次同余方程 $s^2\equiv -1(\mathrm{mod}n)$ 的解一一对应起来。证明细节从略。

再仔细看一看我们可以发现下面几件事：

• $若d=2^q\cdot p，其中q\ge 0，p为奇数，那么N(d)=N(p)$
• $h(d)=1当且仅当4|d-1，h(d)=-1当且仅当4|d-3$
• $N(d)表示d除4余1的因子与除4余3的因子个数之差的四倍$

也就是说只需要分别求出 $d$ 除4余1与除4余3的因子个数
首先明确暴力搜索因子肯定会T，因此很容易想到把 $r$ 标准分解，然后把所有素因子按照mod4分类，在注意到在mod4下有 $3\times 3=1$ 即可算出所需要的结果。

下一步的想法就是：
如果$r^2$的标准分解形式为$r^2=\Pi {p_i}^{2{\alpha }_i}\ast \Pi {q_i}^{2{\beta }_i}$，其中$4|p_i-1，4|q_i-3$，并记$A=\Pi {p_i}^{2{\alpha }_i},B=\Pi {q_i}^{2{\beta }_i}$，那么$r^2$的因子一定可表示为$B$的某一因子与$A$的某一因子的乘积，又由于$A$的所有因子均满足mod4=1，由乘法原理，$r^2$的除4余1因子个数等于$B$的除4余1因子个数乘上$A$的因子个数

下面给出dfs的实现方法

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
typedef long long LL;
 
int prime[2][2340],cnt[2];
//prime[0]记录mod4=1的素数，用cnt[0]记录个数,prime[1]和cnt[1]同理
int po[2340];
//po[i]记录prime[1][i]在r中的幂次
int tmp,ans1=1,ans;
bool vis[44725];
LL r;
 
void dfs(int u,int mod)  //当前讨论prime[1][u],当前余数为mod
{
    if(u>cnt[1])
    {
        if(mod==1)  ans++;
        else  ans--;
        return;
    }
    dfs(u+1,mod);
    for(int i=1;i<=po[u];i++)
    {
        if(i&1)  dfs(u+1,mod^2);
        else  dfs(u+1,mod);
    }
}
void init()
{
    int mx=sqrt(r);
    for(int i=3;i<=mx;i+=2)  if(!vis[i])
    {
        if((i>>1)&1)  prime[1][++cnt[1]]=i;
        else  prime[0][++cnt[0]]=i;
        for(int j=i*2;j<=mx;j+=i)  vis[j]=1;
    }
    for(int i=1;i<=cnt[0];i++,tmp=0)//mod4=1的素数的幂次后期不需要，因此不用存
    {
        while(!(r%prime[0][i]))  tmp++,r/=prime[0][i];
        ans1*=1+(tmp<<1);
    }
    for(int i=1;i<=cnt[1];i++)
    {
        while(!(r%prime[1][i]))  po[i]++,r/=prime[1][i];
        po[i]<<=1;
    }
    //此时r有三种情况1,mod4=1的素数，mod4=3的素数
    if((r>>1)&1)  po[++cnt[1]]=2;
    else  if(r>1)  ans1*=3;
}
void solve()
{
    dfs(1,1);
    printf("%d",ans*ans1<<2);
}
 
int main()
{
    scanf("%lld",&r);
    while(!(r&1))  r>>=1;
    init();  solve();
    return 0;
}

代码意义很明确了吧？那就不做过多说明了。
评测结果：运行时间8ms，比上边给链接的那个(188ms)快好多。

本来到这就结束了。但我翻了翻评测记录，发现有若干4ms的程序以及几个0ms的。我当时就不服了。
我们发现其实po数组大部分数是0，于是乎就有了一个小小的优化

void solve()
{
	sort(po+1,po+cnt[1]+1,greater<int>());
	while(!po[cnt[1]])  cnt[1]--;
    dfs(1,1);
    printf("%d",ans*ans1<<2);
}

运行时间从8ms提高到4ms
后边我就真的没辙了。然后经由一位学长提醒，我发现这题能dp，而且是入门dp

简而言之，用$dp[0][i]$表示使用$q_1,q_2,\cdot \cdot \cdot ,q_i$能产生的$B$的mod4=1的因数个数，$dp[1][i]$表示使用$q_1,q_2,\cdot \cdot \cdot ,q_i$能产生的$B$的mod4=3的因数个数,于是有状态转移方程：

dp[1][i]=dp[0][i-1]*(po[i]+1)/2+dp[1][i-1]*(po[i]/2+1);
dp[0][i]=dp[0][i-1]*(po[i]/2+1)+dp[1][i-1]*(po[i]+1)/2;

这里放一个给出这种思路的学长的blog的链接

运行时间：0ms

Update On 2017.10.9

这一定是这道题的最优解法
——才怪。

事实上在今年的中国东南数学联赛之前我是这么认为的。
不过后来，在东南联赛的考场上，我意外地发现，东南联赛高一组的D1T3本质上就是这个定理。
既然它能出成数竞题（貌似还是个IMOSL），那么就一定是有比DP更优的解法。
我们再来分析下这个DP的式子。这次，我们用数列语言表述它：

记数列 { a n = d p [ 0 ] [ n ] } , { b n = d p [ 1 ] [ n ] } , { α n = p o [ n ] } , \{a_n=dp[0][n]\},\{b_n=dp[1][n]\},\{\alpha_n=po[n]\}, {an​=dp[0][n]},{bn​=dp[1][n]},{αn​=po[n]},那么有
$$a_i=a_{i-1}\cdot \left(\left[\frac{{\alpha }_i}{2}\right]+1\right)+b_{i-1}\cdot \left[\frac{{\alpha }_i+1}{2}\right]$$
$$b_i=a_{i-1}\cdot \left[\frac{{\alpha }_i+1}{2}\right]+b_{i-1}\cdot \left(\left[\frac{{\alpha }_i}{2}\right]+1\right)$$
上述两式相减，得到
$$a_i-b_i=(a_{i-1}-b_{i-1})\cdot \left(\left[\frac{{\alpha }_i}{2}\right]+1-\left[\frac{{\alpha }_i+1}{2}\right]\right)$$
由题目限制，${\alpha }_i$始终为偶数，因此上式等号右边第二项的值恒为1，于是就有$a_n-b_n=1$

再结合前面的全部讨论，我们不难得出结论，所求答案就是$A$的因数个数的4倍。于是我们就可以在分解素因子的同时解决这个问题。
时间复杂度$O(\sqrt{n})$

附：神仙LCA的做法。这里用到了代数数论中的二次代数数的相关知识，虽然用到了不同的道具，但所得的结果和我这个初等方法一致。