问题描述:
给定一个数,问它是否能够表示成两个整数的平方和。
解题思路:
费马平方和定理的表述是:奇素数能表示为两个平方数之和的充分必要条件是该素数被4除余1。
于是一个素数p = 4n+1或者p = 2则必定能表示成两个平方数之和。
同样可以推导出如果两个数均能表示成两个平方数之和,则他们的乘积必定能表示成两个平方数之和:
p = a^2 + b^2;
q = c^2 + d^2;
p*q = (a^2 + b^2)(c^2 + d^2) = (ac)^2 + (ad)^2 + (bc)^2 + (bd)^2 = (ac + bd)^2 + (ac - bd)^2
只要将所求数分解素因子,然后对于每个素数看是否满足条件,一旦有一个素因子不满足费马平方和定理则直接输出"NO",全部满足则输出"YES"
代码如下:
#include <iostream>
#include <cstdlib>
#include <cmath>
#define gcc 10007
typedef __int64 Int;
#define MAX ((Int)1<<63)-1
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
using namespace std;
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
Int p[10] =
{2, 3, 5, 7, 11, 13, 17, 19, 23, 29};
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
Int Gcd(Int a, Int b)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
Int m = 1;
while(m)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
m = a % b;
a = b;
b = m;
}
return a;
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
//计算a*b%n
inline Int Produc_Mod(Int a, Int b, Int mod)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
Int sum = 0;
while(b)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
if(b & 1) sum = (sum + a) % mod;
a = (a + a) % mod;
b /= 2;
}
return sum;
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
//计算a^b%n
inline Int Power(Int a, Int b, Int mod)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
Int sum = 1;
while(b)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
if(b & 1) sum = Produc_Mod(sum, a, mod);
a = Produc_Mod(a, a, mod);
b /= 2;
}
return sum;
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
//Rabin_Miller判素
bool Rabin_Miller(Int n)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
int i, j, k = 0;
Int u, m, buf;
//将n-1分解为m*2^k
if(n == 2)
return true;
if(n < 2 || !(n & 1))
return false;
m = n-1;
while(!(m & 1))
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
k++;
m /= 2;
}
for(i = 0; i < 9; i++)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
if(p[i] >= n)
return true;
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
u = Power(p[i], m, n);
if(u == 1)
continue;
for(j = 0; j < k; j++)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
buf = Produc_Mod(u, u, n);
//看是否有非平凡因子存在
if(buf == 1 && u != 1 && u != n-1)
return false;
u = buf;
}
//如果p[i]^(n-1) % n != 1 那么 n为合数
if(u-1)
return false;
}
return true;
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
Int Pollard_rho(Int n)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
while(1)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
int i = 1;
Int x = rand() % (n-1) + 1;
Int y = x;
Int k = 2;
Int d;
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
do
{
i++;
d = Gcd(n + y - x, n);
if(d > 1 && d < n)
return d;
if(i == k)
y = x, k *= 2;
x = (Produc_Mod(x, x, n) + n - gcc) % n;
}while(y != x);
}
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
Int prime[10000];
int top;
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)
void Prime_Divisor(Int key)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
if( Rabin_Miller(key) )
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
prime[ top ++ ] = key;
}else
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
Int buf = Pollard_rho(key);
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if( Rabin_Miller(buf) )
{
prime[ top ++ ] = buf;
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
}else
{
Prime_Divisor(buf);
}
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
if( Rabin_Miller(key / buf) )
{
prime[ top ++ ] = key / buf;
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
}else
{
Prime_Divisor(key / buf);
}
}
}
int main()
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedBlockStart.gif)
![](http://www.cppblog.com/Images/OutliningIndicators/ContractedBlock.gif)
{
int t, i;
Int n;
scanf("%d", &t);
while(t--)
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
scanf("%I64d", &n);
if(n < 0)
printf("NO\n");
else if(n == 0 || n == 1)
printf("YES\n");
else
![](http://www.cppblog.com/Images/OutliningIndicators/ExpandedSubBlockStart.gif)
{
top = 0;
Prime_Divisor(n);
![](http://www.cppblog.com/Images/OutliningIndicators/InBlock.gif)
for(i = 0; i < top; i++)
if( prime[i] != 2 && (prime[i] - 1) % 4 )
break;
if(i < top)
printf("NO\n");
else
printf("YES\n");
}
}
}
![](http://www.cppblog.com/Images/OutliningIndicators/None.gif)