HDU3939 Sticks and Right Triangle 毕达哥拉斯三元组+容斥原理

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3939

题目大意：给出一个n，找出不超过n的毕达哥拉斯三元组的数目。

分析：这题和POJ1305差不多，只是数据范围大了一些（n为10^12），显然我们不能再去暴力枚举了。

大致看一下m和n的范围可知，m<sqrt(l) ;n<sqrt(l-m*m)=lim 且 n<m;考虑到m和n一奇一偶：

（1）m为偶数时：

    ①如果m<=lim，那么n的取值就是[1,m)中与m互素的数的个数，即m的欧拉函数值；

    ②如果m>lim，我们可以对m进行素因子分解，然后通过容斥原理来计算[1,lim)内与m互素的数的个数，这点可以参考POJ1091 跳蚤的解法.

（2）m为奇数时（这种情况下n是与m互素的偶数）：

    ①如果m<=lim，那么n的选择就是[1,m/2)内与m互素的数的个数（这样我们乘以2以后得到的就一定是偶数了）；

    ②如果m>lim，那么n的选择就是[1,lim/2)内与m互素的数的个数了（原理同上）.

实现代码如下：

#include <iostream>#include <cstdio>#include <cstring>#include <cmath>using namespace std;typedef long long LL;#define maxn 1000010int phi[maxn]; //打欧拉表int prime[maxn],cnt; //打素数表，并纪录10^6内的素数的个数bool flag[maxn]; //flag[i]为0表示i为素数int p[100],num; //纪录一个数的素因子，并纪录因子数numLL ans,l; //ans为不超过l的毕达哥拉斯三元组的数目void init(){   //打素数表    cnt=0;    memset(flag,0,sizeof(flag));    for(int i=2;i<maxn;i++)      if(!flag[i])      {          prime[cnt++]=i;          for(int j=i+i;j<maxn;j+=i)              flag[j]=1;      }    //打欧拉表    for(int i=0;i<maxn;i++) phi[i]=i;    for(int i=2;i<maxn;i+=2) phi[i]>>=1;    for(int i=3;i<maxn;i+=2)      if(phi[i]==i)        for(int j=i;j<maxn;j+=i)          phi[j]=phi[j]-phi[j]/i;}void divide(int x){//分解质因子    num=0;    for(int i=0;i<cnt&&prime[i]*prime[i]<=x;i++)      if(x%prime[i]==0)      {          p[num++]=prime[i];          while(x%prime[i]==0) x/=prime[i];      }    if(x>1) p[num++]=x;}void dfs(int id,int mul,int tot,int x){//搜索[1,x]与m互素的数的个数    if(id==num)    {        if(tot&1) ans=ans-x/mul;        else ans=ans+x/mul;        return ;    }    dfs(id+1,mul*p[id],tot+1,x);    dfs(id+1,mul,tot,x);}int main(){    init();    int t;    cin>>t;    while(t--)    {        scanf("%lld",&l);        int temp=sqrt(l+0.0);        ans=0;        for(int i=1;i<=temp;i++) //枚举m        {            int lim=sqrt(l-(LL)i*i+0.0);  //m确定后n的最大值            if(i&1) //m为奇数            {                divide(i);                if(i<=lim) dfs(0,1,0,i>>1);                else dfs(0,1,0,lim>>1);            }            else //m为偶数            {                if(i<=lim) ans+=phi[i];                else                {                    divide(i);                    dfs(0,1,0,lim);                }            }        }        printf("%lld\n",ans);    }    return 0;}