随笔 - 97, 文章 - 22, 评论 - 81, 引用 - 0
数据加载中……

Pku 3209 From Pythagoras to …(数论)

问题描述:
给定一个数,问它是否能够表示成两个整数的平方和。
解题思路:
费马平方和定理的表述是:奇素数能表示为两个平方数之和的充分必要条件是该素数被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

using namespace std;

Int p[
10= {2357111317192329};

Int Gcd(Int a, Int b)
{
    Int m 
= 1;
    
while(m)
    
{
        m 
= a % b;
        a 
= b;
        b 
= m;
    }

    
return a;
}


//计算a*b%n
inline Int Produc_Mod(Int a, Int b, Int mod)
{
    Int sum 
= 0;
    
while(b)
    
{
        
if(b & 1) sum = (sum + a) % mod;
        a 
= (a + a) % mod;
        b 
/= 2;
    }

    
return sum;
}


//计算a^b%n
inline Int Power(Int a, Int b, Int mod)
{
    Int sum 
= 1;
    
while(b)
    
{
        
if(b & 1) sum = Produc_Mod(sum, a, mod);
        a 
= Produc_Mod(a, a, mod);
        b 
/= 2;
    }

    
return sum;
}


//Rabin_Miller判素
bool Rabin_Miller(Int n)
{
    
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))
    
{
        k
++;
        m 
/= 2;
    }

    
for(i = 0; i < 9; i++)
    
{
        
if(p[i] >= n)
            
return true;

        u 
= Power(p[i], m, n);
        
if(u == 1)
            
continue;
        
for(j = 0; j < k; j++)
        
{
            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;
}


Int Pollard_rho(Int n)
{
    
while(1)
    
{
        
int i = 1;
        Int x 
= rand() % (n-1+ 1;
        Int y 
= x;
        Int k 
= 2;
        Int d;
        
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);
    }

}


Int prime[
10000];
int top;

void Prime_Divisor(Int key)
{
    
if( Rabin_Miller(key) )
    
{
        prime[ top 
++ ] = key;
    }
else
    
{
        Int buf 
= Pollard_rho(key);

        
if( Rabin_Miller(buf) ){
            prime[ top 
++ ] = buf;
        }
else{
            Prime_Divisor(buf);
        }



        
if( Rabin_Miller(key / buf) ){
            prime[ top 
++ ] = key / buf;
        }
else{
            Prime_Divisor(key 
/ buf);
        }

    }

}

int main()
{
    
int t, i;
    Int n;
    scanf(
"%d"&t);
    
while(t--)
    
{
        scanf(
"%I64d"&n);
        
if(n < 0)
            printf(
"NO\n");
        
else if(n == 0 || n == 1)
            printf(
"YES\n");
        
else 
        
{
            top 
= 0;
            Prime_Divisor(n);

            
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");
        }

    }

}

posted on 2009-02-10 21:04 英雄哪里出来 阅读(331) 评论(0)  编辑 收藏 引用 所属分类: ACM


只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理