解題區/解題區 - POJ

POJ 1811 - Prime Test

contents

  1. 1. Problem
  2. 2. Sample Input
  3. 3. Sample Output
  4. 4. Solution

Problem

給一個大整數 $N \le 2^{54}$，判斷是否為質數，若不是則進行因數分解，並且輸出最小的質因數。

備註：大整數分解不可使用一般的 $O(\sqrt{N})$ 算法完成。

Sample Input

1
2
3
2
5
10

Sample Output

1
2
Prime
2

Solution

以下的寫法還不是最快的，而且沒有加上最小質因數的剪枝。若在 pollard-rho 算法中，使用 Brent’s algorithm 取代 Floyd’s algorithm 進行環偵測，同時減少模數使用會變得更快 (常數上)。

討論區提供卡邁克爾數 (Carmichael Number)，讓費馬素性檢驗 (Fermat primality test) 判斷失誤，若用盧卡斯-萊默檢驗法 (Miller–Rabin primality test) 則不受影響。這讓我想到寫過鄒賽爾數 (Zeisel Number) 生成，由數個質數構成大數 $N$ 能測試更多的效能瓶頸吧。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#include <stdio.h> 
#include <vector>
#include <map>
#include <algorithm>
using namespace std;
#define MAXL (50000>>5)+1
#define GET(x) (mark[x>>5]>>(x&31)&1)
#define SET(x) (mark[x>>5] |= 1<<(x&31))
int mark[MAXL];
int P[50000], Pt = 0;
void sieve() {
    register int i, j, k;
    SET(1);
    int n = 46340;
    for (i = 2; i <= n; i++) {
        if (!GET(i)) {
            for (k = n/i, j = i*k; k >= i; k--, j -= i)
                SET(j);
            P[Pt++] = i;
        }
    }
}
long long mul(unsigned long long a, unsigned long long b, unsigned long long mod) { 
    long long ret = 0; 
    for (a %= mod, b %= mod; b != 0; b >>= 1, a <<= 1, a = a >= mod ? a - mod : a) {
        if (b&1) {
            ret += a;
            if (ret >= mod)	ret -= mod;
        }
    } 
    return ret; 
}
void exgcd(long long x, long long y, long long &g, long long &a, long long &b) {
    if (y == 0)
        g = x, a = 1, b = 0;
    else
        exgcd(y, x%y, g, b, a), b -= (x/y) * a;
}
long long llgcd(long long x, long long y) {
    if (x < 0)    x = -x;
    if (y < 0)    y = -y;
    if (!x || !y)    return x + y;
    long long t;
    while (x%y)
        t = x, x = y, y = t%y;
    return y;
}
long long inverse(long long x, long long p) {
    long long g, b, r;
    exgcd(x, p, g, r, b);
    if (g < 0)	r = -r;
    return (r%p + p)%p;
}
long long mpow(long long x, long long y, long long mod) { // mod < 2^32 
    long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
long long mpow2(long long x, long long y, long long mod) {
    long long ret = 1;
    while (y) {
        if (y&1) 
            ret = mul(ret, x, mod);
        y >>= 1, x = mul(x, x, mod);
    }
    return ret % mod;
}
int isPrime(long long p) { // implements by miller-babin
    if (p < 2)	return 0;
    if (p == 2)	return 1;
    if (!(p&1))	return 0;
    long long q = p-1, a, t;
    int k = 0, b = 0;
    while (!(q&1))	q >>= 1, k++;
    for (int it = 0; it < 2; it++) {
        a = rand()%(p-4) + 2;
        t = mpow2(a, q, p);
        b = (t == 1) || (t == p-1);
        for (int i = 1; i < k && !b; i++) {
            t = mul(t, t, p);
            if (t == p-1)
                b = 1;
        }
        if (b == 0)
            return 0;
    }
    return 1;
}
long long pollard_rho(long long n, long long c) {
    register long long x = 2, y = 2, d;
    do { 
        x = mul(x, x, n)+c;
        if (x >= n)	x -= n;
        y = mul(y, y, n)+c;
        if (y >= n)	y -= n;
        y = mul(y, y, n)+c;
        if (y >= n)	y -= n;
        d = llgcd(x-y, n);
        if(d > 1)    return d;
    } while(true);
    return n;
}
void factorize(int n, vector<long long> &f) {
    for (int i = 0; i < Pt && P[i]*P[i] <= n; i++) {
    	if (n%P[i] == 0) {
    		while (n%P[i] == 0)
    			f.push_back(P[i]), n /= P[i];
    	}
    }
    if (n != 1)	f.push_back(n);
}
void llfactorize(long long n, vector<long long> &f, int i = 2) {
    if (n == 1)	
        return ;	
    if (n < 1e+9) {
        factorize(n, f);
        return ;
    }
    if (isPrime(n)) {
        f.push_back(n);
        return ;
    }
    long long d = n;
    for (; d == n; i++)
        d = pollard_rho(n, i);
    llfactorize(d, f, i);
    llfactorize(n/d, f, i);
}
int main() {
    sieve();
    int testcase;
    scanf("%d", &testcase);
    while (testcase--) {
        long long n;
        scanf("%lld", &n);
        
        vector<long long> f;
            
        llfactorize(n, f);
        sort(f.begin(), f.end());
        
        if (f.size() == 1)
            puts("Prime");
        else
            printf("%lld\n", f[0]);
    }
    return 0;
}