POJ 1811 - Prime Test

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;
}
Read More +

POJ 1151 - Atlantis

Problem

給一堆平行兩軸的矩形,請問面積的聯集為何。

Sample Input

1
2
3
4
2
10 10 20 20
15 15 25 25.5
0

Sample Output

1
2
Test case #1
Total explored area: 180.00

Solution

利用掃描線算法,從 x 軸小掃描到大,維度 y 軸上的覆蓋情況。

若單純用離散方法在 2D 網格上填充,算法複雜度是 $O(n^2)$,最慘情況還會再增加到 $O(n^4)$,平均上來說離散方法是可行,但是記憶體是 $O(n^2)$

由於 ITSA 的那題 n 會高達 30000,於是拿這一題來練手,掃描線算法複雜度 $O(n \log n)$,利用線段樹操作時,維護某個區間的完全覆蓋數,當完全覆蓋數等於 0 時,則計算子樹的覆蓋長度總和,雖然沒辦法完全覆蓋,則部分覆蓋可以從子樹計算得到。

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
#include <stdio.h>
#include <string.h>
#include <map>
#include <vector>
#include <math.h>
#include <algorithm>
using namespace std;
const int MAXN = 262144;
class SegmentTree {
public:
struct Node {
int cover; // #cover
double L, R; // value in real line, cover [L, R]
double clen; // cover length
void init(int a = 0, double b = 0, double c = 0, double d = 0) {
cover = a, L = b, R = c, clen = d;
}
} nodes[MAXN];
double Y[MAXN];
void pushDown(int k, int l, int r) {
}
void pushUp(int k, int l, int r) {
if (nodes[k].cover > 0) {
nodes[k].clen = nodes[k].R - nodes[k].L;
} else {
if (l == r)
nodes[k].clen = 0;
else
nodes[k].clen = nodes[k<<1].clen + nodes[k<<1|1].clen;
}
}
void build(int k, int l, int r) {
nodes[k].init(0, Y[l], Y[r+1], 0);
if (l == r)
return ;
int mid = (l + r)/2;
build(k<<1, l, mid);
build(k<<1|1, mid+1, r);
pushUp(k, l, r);
}
// operator, add
void addUpdate(int k, int l, int r, int val) {
nodes[k].cover += val;
if (nodes[k].cover > 0) {
nodes[k].clen = nodes[k].R - nodes[k].L;
} else {
if (l == r)
nodes[k].clen = 0;
else
nodes[k].clen = nodes[k<<1].clen + nodes[k<<1|1].clen;
}
}
void add(int k, int l, int r, int x, int y, int val) {
if (x <= l && r <= y) {
addUpdate(k, l, r, val);
return;
}
pushDown(k, l, r);
int mid = (l + r)/2;
if (x <= mid)
add(k<<1, l, mid, x, y, val);
if (y > mid)
add(k<<1|1, mid+1, r, x, y, val);
pushUp(k, l, r);
}
// query
double r_sum;
void qinit() {
r_sum = 0;
}
void query(int k, int l, int r, int x, int y) {
if (x <= l && r <= y) {
r_sum += nodes[k].clen;
return ;
}
pushDown(k, l, r);
int mid = (l + r)/2;
if (x <= mid)
query(k<<1, l, mid, x, y);
if (y > mid)
query(k<<1|1, mid+1, r, x, y);
pushUp(k, l, r);
}
} tree;
struct Event {
double x, yl, yr;
int val;
Event(double a = 0, double b = 0, double c = 0, int d = 0):
x(a), yl(b), yr(c), val(d) {}
bool operator<(const Event &a) const {
return x < a.x;
}
};
double Lx[32767], Ly[32767], Rx[32767], Ry[32767], K[32767];
int N;
double areaCoverAll() {
vector<Event> events;
vector<double> VY;
map<double, int> RY;
for (int i = 0; i < N; i++) {
double x1, x2, y1, y2;
x1 = Lx[i], x2 = Rx[i];
y1 = Ly[i], y2 = Ry[i];
events.push_back(Event(x1, y1, y2, 1));
events.push_back(Event(x2, y1, y2, -1));
VY.push_back(y1), VY.push_back(y2);
}
sort(events.begin(), events.end());
sort(VY.begin(), VY.end());
VY.resize(unique(VY.begin(), VY.end()) - VY.begin());
for (int i = 0; i < VY.size(); i++) {
tree.Y[i] = VY[i];
RY[VY[i]] = i;
}
tree.build(1, 0, VY.size()-2);
double area = 0, prevX = 0;
for (int i = 0, j; i < events.size(); ) {
if (i > 0) {
tree.qinit();
tree.query(1, 0, VY.size()-2, 0, VY.size()-2);
area += (events[i].x - prevX) * tree.r_sum;
}
j = i;
while (j < events.size() && events[j].x <= events[i].x) {
tree.add(1, 0, VY.size()-2, RY[events[j].yl], RY[events[j].yr] - 1, events[j].val);
j++;
}
prevX = events[i].x;
i = j;
}
return area;
}
int main() {
int testcase, cases = 0;
while (scanf("%d", &N) == 1 && N) {
for (int i = 0; i < N; i++)
scanf("%lf %lf %lf %lf", &Lx[i], &Ly[i], &Rx[i], &Ry[i]);
printf("Test case #%d\n", ++cases);
printf("Total explored area: %.2f\n", areaCoverAll());
puts("");
}
return 0;
}
/*
2
10 10 20 20
15 15 25 25.5
*/
Read More +

POJ 1166 - The Clocks

Description

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
|-------| |-------| |-------|
| | | | | | |
|---O | |---O | | O |
| | | | | |
|-------| |-------| |-------|
A B C
|-------| |-------| |-------|
| | | | | |
| O | | O | | O |
| | | | | | | | |
|-------| |-------| |-------|
D E F
|-------| |-------| |-------|
| | | | | |
| O | | O---| | O |
| | | | | | | |
|-------| |-------| |-------|
G H I

(Figure 1)

There are nine clocks in a 3*3 array (figure 1). The goal is to return all the dials to 12 o’clock with as few moves as possible. There are nine different allowed ways to turn the dials on the clocks. Each such way is called a move. Select for each move a number 1 to 9. That number will turn the dials 90’ (degrees) clockwise on those clocks which are affected according to figure 2 below.

Move Affected clocks
1 ABDE
2 ABC
3 BCEF
4 ADG
5 BDEFH
6 CFI
7 DEGH
8 GHI
9 EFHI

(Figure 2)

Input

Your program is to read from standard input. Nine numbers give the start positions of the dials. 0=12 o’clock, 1=3 o’clock, 2=6 o’clock, 3=9 o’clock.

Output

Your program is to write to standard output. Output a shortest sorted sequence of moves (numbers), which returns all the dials to 12 o’clock. You are convinced that the answer is unique.

Sample Input

1
2
3
3 3 0
2 2 2
2 1 2

Sample Output

1
4 5 8 9

Solution

原來之前寫過這一題,請參閱 d730: 升旗典礼 ——加强版

用 BFS 從 0 0 0 0 0 0 0 0 0 開始倒推回去,在 mark 的部份,用 4 進制 bitmask。

這題在 POJ 似乎只有一組測資,在 Zerojudge 則有非常多筆測資。

而事實上網路上還有另外一個高斯消去法,先預求反矩陣,之後直接帶入計算即可。沒有仔細去研究,不過聽說那解法限定一定有合法解,如何判定不合法解又是另外一個問題。總之,用了 BFS 過了就沒仔細去探討其他解法。

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
#include <stdio.h>
#include <queue>
using namespace std;
char mark[1<<18] = {}, ans[1<<18], best[1<<18];
char dirp[4] = {-3, 1, 1, 1};
int next[1<<18] = {};
int pow4[10] = {1};
int move[10][7] = { {},
{1, 2, 4, 5, 0},
{1, 2, 3, 0},
{2, 3, 5, 6, 0},
{1, 4, 7, 0},
{2, 4, 5, 6, 8, 0},
{3, 6, 9, 0},
{4, 5, 7, 8, 0},
{7, 8, 9, 0},
{5, 6, 8, 9, 0}
};
void trans(int N, char s[]) {
int i = 0;
while(N) {
s[i++] = N%4, N /= 4;
}
}
void build() {
int a, b, c, qt = 0, t, p;
queue<int> Q;
best[0] = 0, mark[0] = 1;
Q.push(0);
while(!Q.empty()) {
int state = Q.front();
Q.pop();
char buf[10] = {};
trans(state, buf);
for(int i = 1; i < 10; i++) {
int nstate = state;
for(int j = 0; move[i][j]; j++)
nstate -= pow4[move[i][j] - 1] * dirp[buf[move[i][j] - 1]];
if(mark[nstate] == 0) {
mark[nstate] = 1;
next[nstate] = state;
best[nstate] = best[state] + 1;
ans[nstate] = i + '0';
Q.push(nstate);
}
}
}
}
int print_f = 0;
void Print(int state) {
if(state) {
Print(next[state]);
if(print_f)
putchar(' ');
print_f = 1;
putchar(ans[state]);
}
}
int main() {
for(int i = 1; i < 10; i++)
pow4[i] = pow4[i-1] * 4;
build();
while(true) {
int state = 0, x;
for(int i = 0; i < 9; i++) {
if(scanf("%d", &x) != 1)
return 0;
state |= x * pow4[i];
}
print_f = 0;
Print(state);
puts("");
}
return 0;
}
Read More +