HDU 3939-Sticks and Right Triangle(毕达哥拉斯三元组+欧拉函数)

题目地址：HDU 3939
题意：给出勾股方程X^2+Y^2=Z^2，满足X,Y,Z不超过L，问存在多少个解满足X,Y,Z两两互素。
思路：其实就是求有多少个毕达哥拉斯本原三元组，因为数据很大，所以不能像上一道题一样，直接暴力枚举。那下面我们来分析一下：
由X^2+Y^2+Z^2我们可以知道，方程的解为X=m^2-n^2，Y=2*m*n，Z=m^2+n^2，其中m>n，m与n一奇一偶，且gcd(m,n)=1。由式子我们可以知道Z最大，这样只需要Z=m^2+n^2<=L即可。所以我们可以确定一下m和n的大体范围，m<=sqrt(L) ;n<=sqrt(l-m*m)=Nmax。
1）当m为偶数,此时n为奇数：
1.如果m<=Nmax，那么n的取值就是[1,m)中与m互素的数的个数(m的欧拉函数值)；
2.如果m>Nmax，对m进行素因子分解，然后通过容斥原理来计算[1,Nmax)内与m互素的个数。
2）当m为奇数时，此时n为偶数：
1.如果m<=Nmax，那么n的选择就是[1,m/2)内与m互素的数的个数（这样我们乘以2以后得到的就一定是偶数了）；
2.如果m>Nmax，那么n的选择就是[1,Nmax/2)内与m互素的个数。

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <iostream>
#include <sstream>
#include <algorithm>
#include <set>
#include <queue>
#include <stack>
#include <map>
#include <bitset>
#pragma comment(linker, "/STACK:102400000,102400000")
using namespace std;
typedef __int64 LL;
const int inf=0x3f3f3f3f;
const double pi= acos(-1.0);
const double esp=1e-6;
using namespace std;
const int Maxn=1e6+10;
bitset<Maxn>pri;
LL prime[Maxn];
LL sprime[Maxn];
LL phi[Maxn];
LL k,cnt;
void is_prime()
{
    pri.set();
    for(LL i=2; i<Maxn; i++) {
        if(pri[i]) {
            prime[k++]=i;
            for(LL j=i+i; j<Maxn; j+=i)
                pri[j]=0;
        }
    }
}

void Divide(LL n)
{
    cnt=0;
    for(LL i=0;i<k&&prime[i]*prime[i]<=n; i++) {
        if(n%prime[i]==0) {
            sprime[cnt++]=prime[i];
            while(n%prime[i]==0)
                n/=prime[i];
        }
    }
    if(n>1)
        sprime[cnt++]=n;
}
LL Ex(LL n)
{
    LL ans=0;
    LL tmp,flag;
    LL i,j;
    for(i=1; i<(LL)(1<<cnt); i++) {
        tmp=1;
        flag=0;
        for(j=0; j<cnt; j++)
            if(i&((LL)(1<<j))) {
                flag++;
                tmp*=sprime[j];
            }
        if(flag&1)
            ans+=n/tmp;
        else
            ans-=n/tmp;
    }
    return ans;
}
void Euler()
{
    LL i,j;
    memset(phi,0,sizeof(phi));
    phi[1]=1;
    for(i=2;i<=Maxn;i++)
    {
        if(!phi[i])
        {
            for(j=i;j<=Maxn;j+=i)
            {
                if(!phi[j])
                    phi[j]=j;
                phi[j]=phi[j]/i*(i-1);
            }
        }
    }
}
int main()
{
    int T;
    LL L,res;
    is_prime();
    Euler();
    scanf("%d",&T);
    while(T--)
    {
        res=0;
        scanf("%lld",&L);
        LL t=(LL)sqrt(1.0*L);
        for(LL m=1;m<=t;m++)
        {
            LL Nmax=(LL)sqrt(L-m*m);
            if(m&1)
            {
                Divide(m);
                if(m<=Nmax){
                    res+=m/2-Ex(m/2);
                }
                else{
                    res+=Nmax/2-Ex(Nmax/2);
                }
            }
            else
            {
                if(m<=Nmax) res+=phi[m];
                else
                {
                    Divide(m);
                    res+=Nmax-Ex(Nmax);
                }
            }
        }
        printf("%lld\n",res);
    }
    return 0;
}