给定正整数b，求最大的整数a，满足a*(a+b) 为完全平方数

http://www.cnblogs.com/snake-hand/archive/2013/06/05/3120010.html

题意：给定正整数b，求最大的整数a，满足a*(a+b) 为完全平方数，1 <= b <= 10^9

 

解题思路：

我们设a，b的最大公约数为g，则a与a+b的最大公约数也为g，因为最大公约数有性质：gcd(a,b)=gcd(a,a+b)

这样我们就可以进一步化简有a*(a+b)=g^2*a1*(a1+b1)，其中a1与b1就一定互素了，因为已经约去最大公约数g了，所以得到a1与a1+b1也互素。

由于要为完全平方数，所以可以设x^2=a1，y^2=a1+b1   进而推出：b1=y^2-x^2=(y-x)*(y+x)

现在我们令n=y+x，m=y-x  很明显n>m，由于a1与a1+b1互素，所以有gcd(x^2,y^2)=1，进而有gcd(x,y)=1

 

到了这里本题就可以这样做了：

先求出b的所有约数，这样g就等于b/b的约数，然后又分别求出b的约数的约数，假设为arr2[j]，由于n>m，我们只需要枚举到sqrt(arr1[i])即可，然后解

出x=(n-m)/2，y=(n+m)/2，但是前提是要保证都能整除，然后再判断gcd(x,y)=1,两层for循环，记录最大值即可。至于怎样求一个数因子，前面的文章已

经给出，就是素因子分解加上dfs，由于因子一般不会很多，所以一般没有问题。建议不要看我的代码，很乱的。

#include <iostream>
#include <string.h>
#include <algorithm>
#include <stdio.h>

using namespace std;

typedef long long LL;
const int N=1000010;
const int M=1050;

bool prime
;
LL p
;
LL pr1[M];
LL kk1[M];

LL pr2[M];
LL kk2[M];
LL k=0;
LL c1,r1;
LL arr1[M];
LL c2,r2;
LL arr2[M];

void isprime()
{
LL i,j;
memset(prime,true,sizeof(prime));
for(i=2;i<N;i++)
{
if(prime[i])
{
p[k++]=i;
for(j=i+i;j<N;j+=i)
{
prime[j]=false;
}
}
}
}

void CalFactor1(LL n)
{
LL t=n,i,a;c1=0;
for(i=0;p[i]*p[i]<=n;i++)
{
a=0;
if(n%p[i]==0)
{
pr1[c1]=p[i];
while(n%p[i]==0)
{
a++;
n/=p[i];
}
kk1[c1]=a;
c1++;
}
}
if(n>1)
{
pr1[c1]=n;
kk1[c1]=1;
c1++;
}
}

void dfs1(LL dep, LL product)
{
if ( dep == c1 )
{
arr1[r1++]=product;
return;
}
for ( LL i = 0; i <= kk1[dep]; ++i )
{
dfs1(dep + 1, product);
product *= pr1[dep];
}
}

void CalFactor2(LL n)
{
LL t=n,i,a;c2=0;
for(i=0;p[i]*p[i]<=n;i++)
{
a=0;
if(n%p[i]==0)
{
pr2[c2]=p[i];
while(n%p[i]==0)
{
a++;
n/=p[i];
}
kk2[c2]=a;
c2++;
}
}
if(n>1)
{
pr2[c2]=n;
kk2[c2]=1;
c2++;
}
}

void dfs2(LL dep, LL product)
{
if ( dep == c2 )
{
arr2[r2++]=product;
return;
}
for ( LL i = 0; i <= kk2[dep]; ++i )
{
dfs2(dep + 1, product);
product *= pr2[dep];
}
}

LL gcd(LL a,LL b)
{
return b? gcd(b,a%b):a;
}

int main()
{
LL m,n,i,j,g,b,val,t;
isprime();
cin>>t;
while(t--)
{
cin>>b;
r1=r2=0;
CalFactor1(b);
dfs1(0,1);
sort(arr1,arr1+r1);
LL max=0;
for(i=0;i<r1;i++)
{
r2=0;
g=b/arr1[i];
memset(arr2,0,sizeof(arr2));
CalFactor2(arr1[i]);
dfs2(0,1);
for(j=0;j<r2;j++)
{
if(arr2[j]*arr2[j]<arr1[i])
{
n=arr1[i]/arr2[j];
m=arr2[j];
if((n-m)%2==0&&(n+m)%2==0&&gcd((n+m)/2,(n-m)/2)==1)
{
val=g*((n-m)/2)*((n-m)/2);
if(val>max) max=val;
}
}
}
}
cout<<max<<endl;
}
return 0;
}