POJ1811 Prime Test

最新推荐文章于 2024-09-12 13:19:42 发布

tju_virus

最新推荐文章于 2024-09-12 13:19:42 发布

阅读量379

点赞数

分类专栏：数学文章标签：算法 c

本文链接：https://blog.csdn.net/tju_virus/article/details/7859678

版权

数学专栏收录该内容

9 篇文章 0 订阅

订阅专栏

经典的数学题

题意：判断一个数是不是素数，如果是输出Prime，否则输出这个数的质因数中最小的一个

题解：用Miller_rabbin算法判断这个数是不是素数，如果不是就用Pollard_rho算法进行因式分解

以下为代码，但是在poj上不知道为什么如果我用G++交会报runtime error而用C++交就能过，求大牛指教

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <map>
#include <iostream>
using namespace std;
#define MAXN 10
#define C 16381
typedef __int64  I64;
map<I64,int> fac;

//a*b可能会越界，用这个函数来实现乘法
I64 multi(I64 a, I64 b, I64 n){
    I64 tmp = a % n, s = 0;
    while(b){
        if(b & 1) s = (s + tmp) % n;
        tmp = (tmp + tmp) % n;
        b >>= 1;
    }
    return s;
}

//快速幂取模a的b次方
I64 Pow(I64 a, I64 b, I64 n){
    I64 tmp = a % n, s = 1;
    while(b){
        if(b & 1) s = multi(s, tmp, n);
        tmp = multi(tmp, tmp, n);
        b >>= 1;
    }
    return s;
}

//验证a的n-1次方
int witness(I64 a, I64 n){
    I64 u = n - 1, t = 0, i, x, y;
    while(!(u & 1))    u >>= 1, t++;
    x = Pow(a, u, n);
    for(i = 1; i <= t; i ++){
        y = multi(x, x, n);
        if(y == 1 && x != 1 && x != n -1) return 1;
        x = y;
    }
    if(x != 1) return 1;
    return 0;
}

//Miller_Rabbin法判断n是否为素数
int test(I64 n){
    I64 a;
    int i;
    if(n == 2) return 1;
    if(n < 2 || !(n & 1)) return 0;
    srand((I64)time(0));
    for(i = 0; i < MAXN; i ++){
        a = ((I64) rand()) % (n - 2) + 2;
        if(witness(a, n)) return 0;
    }
    return 1;
}

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

//pollard_rho算法核心
I64 pollard_rho(I64 n, I64 c){
    I64 x, y, d, i = 1, k = 2;
    srand((I64)time(0));
    x = ((I64) rand()) % (n - 1) + 1;
    y = x;
    while(1){
        i ++;
        x = (multi(x, x, n) + c) % n;
        d = gcd(y - x + n, n);
        if(d != 1 && d != n) return d;
        if(y == x) return n;
        if(i == k) y = x, k <<= 1;
    }
}

//把n分解质因数，c为循环调用随机算法的次数
void find(I64 n, I64 c){
    I64 r;
    if(n <= 1) return;
    if(test(n)){
        fac[n] ++;
        return;
    }
    r = pollard_rho(n, c--);
    find(n / r, c);
    find(r, c);
}

//调用find(n, C);
//一般C用120就差不多了 

int main()
{
	int t;
	I64 n;
	scanf("%d",&t);
	while(t--)
	{
		scanf("%I64d",&n);
		fac.clear();
		if(test(n))
			printf("Prime\n");
		else
		{
			find(n,120);
			cout << fac.begin()->first << endl;
		}
	}
	return 0;
}