[BZOJ 1041] 圆上的整点

由于博主水平有限，文中如有错误请指正，谢谢
本思路为博主原创（至少我没在网上找到和我思路一致的），转载请注明出处

—————————-分割线———————————-

先摆原题

Description
　　求一个给定的圆(x2+y2=r2)，在圆周上有多少个点的坐标是整数
Input
　　只有一个正整数n,n≤2000000000
Output
　　整点个数
Sample Input
4
Sample Output
4

初见这道题时候果断不会。无奈上网找题解。
我在网上找到的做法是这样的：http://hzwer.com/1457.html
感觉我此生所学都用进去了qwq。加上严重的阅读障碍症，我果断放弃这道题。

后来胡乱翻数论书，偶然看见一个神奇的定理，也就是我们今天的主角：

设n≥1，则不定方程x2+y2=n的解数为：

N(p)=4∑d|ph(d)

其中函数h(d)定义为：

h(d)=⎧⎩⎨⎪⎪1,0,(−1)(d−1)/2,d=12|delse

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

• 若d=2q⋅p，其中q≥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∗3=1即可算出所需要的结果。

下一步的想法就是：
如果r2的标准分解形式为r2=Πpi2αi∗Πqi2βi，其中4|pi−1，4|qi−3，并记A=Πpi2αi,B=Πqi2βi，那么r2的因子一定可表示为B的某一因子与A的某一因子的乘积，又由于A的所有因子均满足mod4=1，由乘法原理，r2的除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>());    //greater需要包含头文件functional,当然你也可以写一个cmp函数    while(!po[cnt[1]])  cnt[1]--;    dfs(1,1);    printf("%d",ans*ans1<<2);}

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

简而言之，用dp[0][i]表示使用q1,q2,⋅⋅⋅,qi能产生的B的mod4=1的因数个数，dp[1][i]表示使用q1,q2,⋅⋅⋅,qi能产生的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;

不过这个dp的思路不是我想到的，就不贴代码了。
甩一个给出这种思路的dalao的blog的链接：
http://blog.csdn.net/qq_34564984/article/details/73274937

运行时间：0ms

Update On 2017.10.9

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

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

记数列{an=dp[0][n]},{bn=dp[1][n]},{αn=po[n]},那么有

ai=ai−1∗([αi2]+1)+bi−1∗[αi+12]

bi=ai−1∗[αi+12]+bi−1∗([αi2]+1)

上述两式相减，得到

ai−bi=(ai−1−bi−1)∗([αi2]+1−[αi+12])

由题目限制，αi始终为偶数，因此上式等号右边第二项的值恒为1，于是就有an−bn=1

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

附：LCAjulao的做法。这里用到了代数数论中的二次代数数的相关知识，虽然用到了不同的道具，但所得的结果和我这个初等方法一致。
http://blog.csdn.net/liu_cheng_ao/article/details/69296100