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 +

a761. 罐頭的記憶體

Problem

模擬作業系統的 dynamic memory allocation,

  • 宣告連續一段記憶體配置給某個程序,並且標記辨識碼。
  • 若宣告失敗,則輸出碰撞區域中,辨識碼最小的數字。
  • 詢問記憶體存取時,是否合法,若合法輸出程序的辨識碼。

Sample Input

1
2
3
4
5
6
7
8
7
load 29
map 1 25 39
map 2 23 24
map 3 22 25
map 4 25 40
store 33
store 22

Sample Output

1
2
3
4
5
6
7
fail to load from 29
region 1 created
region 2 created
fail to create region 3, overlap with region 1
fail to create region 4, overlap with region 1
store region to 1
fail to store to 22

Solution

溫馨提示:此題的測資隨機,單純離散化套用二分查找,編號最小使用線性搜索即可通過。by 蔡神 asas

儘管直觀作法可以通過,來討論看看別種作法吧。

  • 線段樹 + 延遲標記 + 離線處理
    離散結果最多有 $O(N + Q)$ 的離散點,單筆操作複雜度為 $O(\log Q)$,空間複雜度 $O(Q)$,空間消耗量相較於其他算法最多 (附加延遲標記的使用)。
  • 分塊算法 + 離線處理
    利用分塊 Unrolled Linked List 的快取優勢來降複雜度,空間複雜度仍然是 $O(Q)$,單筆操作複雜度為 $O(\sqrt{Q})$。空間複雜度的常數比線段樹少一半以上,別忘詢問離線是共同的空間成本。
  • 分塊在線
    直接宣告 $O(2^{32})$ 進行配置,因此分塊大小為 $O(\sqrt{2^{32}})$。將配置連續區段用區間 $[l, r, pid]$ 的方式儲存在塊之中,用一個平衡樹維護每一個塊的狀態。由於要輸出編號最小的,詢問複雜度至少為 $O(65536)$,為了加速詢問可以利用剪枝 (塊的最小值仍大於當前的答案) 來完成。空間複雜度跟實際情況有關,這一做法在當前數據中是最快,記憶體最小的一種辦法,因為在線很多詢問點是沒必要實際存在的。

由於這一題沒有釋放記憶體操作,從均攤上可以清楚,分塊離線的修改離散點的次數最多為 $O(Q)$,也就是每一個點均被塗色僅僅一次,剩下就是在詢問上加速。

線段樹

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
#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1<<21, INF = 0x3f3f3f3f;
class SegmentTree {
public:
struct Node {
int val, del;
void init(int a = 0, int b = 0) {
val = a, del = b;
}
} nodes[MAXN];
void pushDown(int k, int l, int r) {
if (nodes[k].del != INF) {
nodes[k<<1].val = min(nodes[k<<1].val, nodes[k].del);
nodes[k<<1].del = min(nodes[k<<1].del, nodes[k].del);
nodes[k<<1|1].val = min(nodes[k<<1|1].val, nodes[k].del);
nodes[k<<1|1].del = min(nodes[k<<1|1].del, nodes[k].del);
nodes[k].del = INF;
}
}
void pushUp(int k, int l, int r) {
nodes[k].val = min(nodes[k<<1].val, nodes[k<<1|1].val);
}
void build(int k, int l, int r) {
nodes[k].init(INF, INF);
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);
}
void assignUpdate(int k, int l, int r, int val) {
nodes[k].del = min(nodes[k].del, val);
nodes[k].val = min(nodes[k].val, val);
}
void assign(int k, int l, int r, int x, int y, int val) {
if (x <= l && r <= y) {
assignUpdate(k, l, r, val);
return;
}
pushDown(k, l, r);
int mid = (l + r)/2;
if (x <= mid) assign(k<<1, l, mid, x, y, val);
if (y > mid) assign(k<<1|1, mid+1, r, x, y, val);
pushUp(k, l, r);
}
int query(int k, int l, int r, int x) {
if (x <= l && r <= x || nodes[k].val == INF)
return nodes[k].val;
pushDown(k, l, r);
int mid = (l + r)/2, ret;
if (x <= mid) ret = query(k<<1, l, mid, x);
else ret = query(k<<1|1, mid+1, r, x);
pushUp(k, l, r);
return ret;
}
int query(int k, int l, int r, int x, int y) {
if (x <= l && r <= y || nodes[k].val == INF)
return nodes[k].val;
pushDown(k, l, r);
int mid = (l + r)/2, ret = INF;
if (x <= mid)
ret = min(ret, query(k<<1, l, mid, x, y));
if (y > mid)
ret = min(ret, query(k<<1|1, mid+1, r, x, y));
pushUp(k, l, r);
return ret;
}
} tree;
const int MAXQ = 524288;
char cmd[MAXQ][8];
int args[MAXQ][3];
int main() {
int N;
while (scanf("%d", &N) == 1) {
map<int, int> R;
for (int i = 0; i < N; i++) {
scanf("%s", &cmd[i]);
int a, b, c;
if (cmd[i][0] == 'l')
scanf("%d", &a), R[a] = 0;
else if (cmd[i][0] == 'm')
scanf("%d %d %d", &a, &b, &c), R[b] = R[c] = 0;
else
scanf("%d", &a), R[a] = 0;
args[i][0] = a, args[i][1] = b, args[i][2] = c;
}
int size = 0;
for (auto &x : R)
x.second = ++size;
tree.build(1, 1, size);
for (int i = 0; i < N; i++) {
int a, b, c, t;
a = args[i][0], b = args[i][1], c = args[i][2];
if (cmd[i][0] == 'l') {
t = tree.query(1, 1, size, R[a]);
if (t == INF)
printf("fail to load from %d\n", a);
else
printf("load from region %d\n", t);
} else if (cmd[i][0] == 'm') {
t = tree.query(1, 1, size, R[b], R[c]);
if (t != INF)
printf("fail to create region %d, overlap with region %d\n", a, t);
else
printf("region %d created\n", a), tree.assign(1, 1, size, R[b], R[c], a);
} else {
t = tree.query(1, 1, size, R[a]);
if (t == INF)
printf("fail to store to %d\n", a);
else
printf("store to region %d\n", t);
}
}
}
return 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
#include <bits/stdc++.h>
using namespace std;
const int INF = 0x3f3f3f3f, MAXQ = 524288;
class Unrolled {
public:
int PILE, mask, shift;
vector< vector<int> > nodes;
vector<int> dirty;
void alloc(int size) {
nodes.clear(), dirty.clear();
for (PILE = 1, shift = 0; PILE * PILE < size; PILE <<= 1, shift++);
mask = PILE - 1;
for (int l = 0, r; l < size; l = r+1) {
r = min(l+PILE-1, size-1);
nodes.push_back(vector<int>(r-l+1, INF));
dirty.push_back(INF);
}
}
int operator[](const int x) const {
return nodes[x>>shift][x&mask];
}
int empty(int l, int r) {
int bl = l>>shift, br = r>>shift;
int ret = INF;
if (bl == br) {
for (int i = l; i <= r; i++)
ret = min(ret, (*this)[i]);
return ret;
}
for (int i = bl+1; i < br; i++)
ret = min(ret, dirty[i]);
for (int i = (bl+1) * PILE-1; i >= l; i--)
ret = min(ret, (*this)[i]);
for (int i = (br) * PILE; i <= r; i++)
ret = min(ret, (*this)[i]);
return ret;
}
void fill(int l, int r, int val) {
int bl = l>>shift, br = r>>shift;
for (int i = bl; i <= br; i++)
dirty[i] = min(dirty[i], val);
for (int i = l, bi; i <= r; i++)
nodes[i>>shift][i&mask] = val;
}
} mem;
char cmd[MAXQ][8];
int args[MAXQ][3];
int main() {
int N;
while (scanf("%d", &N) == 1) {
map<int, int> R;
for (int i = 0; i < N; i++) {
scanf("%s", &cmd[i]);
int a, b, c;
if (cmd[i][0] == 'l')
scanf("%d", &a), R[a] = 0;
else if (cmd[i][0] == 'm')
scanf("%d %d %d", &a, &b, &c), R[b] = R[c] = 0;
else
scanf("%d", &a), R[a] = 0;
args[i][0] = a, args[i][1] = b, args[i][2] = c;
}
int size = 0;
for (auto &x : R)
x.second = size++;
mem.alloc(size);
for (int i = 0; i < N; i++) {
int a, b, c, t;
a = args[i][0], b = args[i][1], c = args[i][2];
if (cmd[i][0] == 'l') {
t = mem[R[a]];
if (t == INF)
printf("fail to load from %d\n", a);
else
printf("load from region %d\n", t);
} else if (cmd[i][0] == 'm') {
t = mem.empty(R[b], R[c]);
if (t != INF)
printf("fail to create region %d, overlap with region %d\n", a, t);
else
printf("region %d created\n", a), mem.fill(R[b], R[c], a);
} else {
t = mem[R[a]];
if (t == INF)
printf("fail to store to %d\n", a);
else
printf("store to region %d\n", t);
}
}
}
return 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
#include <bits/stdc++.h>
using namespace std;
const int INF = 0x3f3f3f3f;
class Unrolled {
public:
struct Interval {
long long l, r;
int pid;
Interval(long long a = 0, long long b = 0, int c = 0):
l(a), r(b), pid(c) {}
bool operator<(const Interval &x) const {
return l < x.l || (l == x.l && r > x.r);
}
bool include(int x) const {
return l <= x && x <= r;
}
};
long long PILE, mask;
int shift;
vector< set<Interval> > nodes;
vector<int> dirty;
void alloc(long long size) {
nodes.clear(), dirty.clear();
for (PILE = 1, shift = 0; PILE * PILE < size; PILE <<= 1, shift++);
mask = PILE - 1;
for (long long l = 0, r; l < size; l = r+1) {
r = min(l+PILE-1, size-1);
nodes.push_back(set<Interval>());
dirty.push_back(INF);
}
}
int operator[](const long long x) const {
const set<Interval> &s = nodes[x>>shift];
auto it = s.upper_bound(Interval(x, x));
it = it == s.begin() ? it : --it;
return it == s.end() || !(it->include(x)) ? INF : it->pid;
}
int empty(long long l, long long r) {
int bl = l>>shift, br = r>>shift;
int ret = INF;
if (bl == br) {
set<Interval> &s = nodes[bl];
auto st = s.upper_bound(Interval(l, l));
st = st == s.begin() ? st : --st;
for (auto it = st; it != s.end(); it++) {
if (max(l, it->l) <= min(r, it->r))
ret = min(ret, it->pid);
if (it->l > r)
break;
}
return ret;
}
for (int i = bl+1; i < br; i++)
ret = min(ret, dirty[i]);
for (auto it = nodes[bl].rbegin(); it != nodes[bl].rend() && ret > dirty[bl]; it++) {
if (max(l, it->l) <= min(r, it->r)) {
ret = min(ret, it->pid);
} else {
break;
}
}
for (auto it = nodes[br].begin(); it != nodes[br].end() && ret > dirty[br]; it++) {
if (max(l, it->l) <= min(r, it->r)) {
ret = min(ret, it->pid);
} else {
break;
}
}
return ret;
}
void fill(long long l, long long r, int val) {
int bl = l>>shift, br = r>>shift;
for (int i = bl; i <= br; i++)
dirty[i] = min(dirty[i], val);
if (bl == br) {
nodes[bl].insert(Interval(l, r, val));
return ;
}
for (int i = bl+1; i < br; i++)
nodes[i].insert(Interval(i*PILE, (i+1)*PILE-1, val));
nodes[bl].insert(Interval(l, (bl+1) * PILE-1, val));
nodes[br].insert(Interval(br*PILE, r, val));
}
} mem;
int main() {
int N;
char cmd[8];
while (scanf("%d", &N) == 1) {
mem.alloc(1LL<<31);
for (int i = 0; i < N; i++) {
scanf("%s", &cmd);
int a, b, c, t;
if (cmd[0] == 'l') {
scanf("%d", &a);
t = mem[a];
if (t == INF)
printf("fail to load from %d\n", a);
else
printf("load from region %d\n", t);
} else if (cmd[0] == 'm') {
scanf("%d %d %d", &a, &b, &c);
t = mem.empty(b, c);
if (t != INF)
printf("fail to create region %d, overlap with region %d\n", a, t);
else
printf("region %d created\n", a), mem.fill(b, c, a);
} else {
scanf("%d", &a);
t = mem[a];
if (t == INF)
printf("fail to store to %d\n", a);
else
printf("store to region %d\n", t);
}
}
}
return 0;
}
Read More +

b436. 圖片的梯度再進化

Problem

可以參考 中央大學影像處理實驗室 課程 講義第八章 邊擷取

除了索貝爾運算 Sobel operator,還有如 Prewitt operator、Kirsch operator … 等,來進行邊擷取。

Sample Input

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

Sample Output

1
2
3
2 2
191 191 191 128 128 128
128 128 128 64 64 64

Solution

比較麻煩是在處理邊緣擴充,寫一個函數去特判延伸的 pixel 位置。

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
#include <bits/stdc++.h>
using namespace std;
class IMAGE {
public:
struct Pixel {
double r, g, b, a;
Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
r(x), g(y), b(z), a(w) {}
void read() {
scanf("%lf %lf %lf", &r, &g, &b);
}
Pixel operator-(const Pixel &x) const {
return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
}
Pixel operator+(const Pixel &x) const {
return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
}
Pixel operator*(const double x) const {
return Pixel(r*x, g*x, b*x, a*x);
}
Pixel operator/(const double x) const {
return Pixel(r/x, g/x, b/x, a/x);
}
void print() {
printf("%d %d %d", (int)round(r), (int)round(g), (int)round(b));
}
};
int W, H;
static const int MAXN = 256;
Pixel data[MAXN][MAXN], tmp[MAXN][MAXN];
void read() {
scanf("%d %d", &W, &H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].read();
}
Pixel pabs(Pixel x) {
return Pixel(fabs(x.r), fabs(x.g), fabs(x.b));
}
Pixel getPixel(int x, int y) {
if (x >= 0 && y >= 0 && x < H && y < W)
return data[x][y];
if (y < 0) return data[min(max(x, 0), H-1)][0];
if (y >= W) return data[min(max(x, 0), H-1)][W-1];
if (x < 0) return data[0][min(max(y, 0), W-1)];
if (x >= H) return data[H-1][min(max(y, 0), W-1)];
return Pixel(0, 0, 0);
}
void border() {
const int dx[] = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
const int dy[] = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
const int xw[] = {-1, 0, 1, -2, 0, 2, -1, 0, 1};
const int yw[] = {-1, -2, -1, 0, 0, 0, 1, 2, 1};
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
Pixel Dx(0, 0, 0), Dy(0, 0, 0);
for (int k = 0; k < 9; k++) {
Dx = Dx + getPixel(i+dx[k], j+dy[k]) * xw[k];
Dy = Dy + getPixel(i+dx[k], j+dy[k]) * yw[k];
}
Dx = pabs(Dx), Dy = pabs(Dy);
tmp[i][j] = (Dx + Dy)/8.0;
}
}
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j] = tmp[i][j];
}
void print() {
printf("%d %d\n", W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
}
} test;
int main() {
test.read();
test.border();
test.print();
return 0;
}
Read More +

b434. 圖片的梯度

Problem

可以參考 中央大學影像處理實驗室 課程 講義第八章 邊擷取

邊緣擷取可以靠相鄰像素的差異來判斷是否是邊緣,若顏色差異很大則表示是邊。然而圖片像素是不連續的離散點,因此沒有辦法利用微積分學到的微分來得精準,而是要採用差分 (減法),至於遮罩要取多大是個問題,一般而言還會搭配降噪的矩陣來完成。

如 Prewitt operator、Sobel operator、Kirsch operator … 等,矩陣都有搭配降噪處理。之所以用矩陣表示,是為了讓 GPU 快速地計算數值。

Sample Input

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

Sample Output

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

Solution

按照公式模擬即可。

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
#include <bits/stdc++.h>
using namespace std;
class IMAGE {
public:
struct Pixel {
double r, g, b, a;
Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
r(x), g(y), b(z), a(w) {}
void read() {
scanf("%lf %lf %lf", &r, &g, &b);
}
Pixel operator-(const Pixel &x) const {
return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
}
Pixel operator+(const Pixel &x) const {
return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
}
Pixel operator*(const double x) const {
return Pixel(r*x, g*x, b*x, a*x);
}
Pixel operator/(const double x) const {
return Pixel(r/x, g/x, b/x, a/x);
}
void print() {
printf("%d %d %d", (int)round(r), (int)round(g), (int)round(b));
}
};
int W, H;
static const int MAXN = 256;
Pixel data[MAXN][MAXN];
void read() {
scanf("%d %d", &W, &H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].read();
}
Pixel pabs(Pixel x) {
return Pixel(fabs(x.r), fabs(x.g), fabs(x.b));
}
void border() {
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
Pixel dx(0, 0, 0), dy(0, 0, 0);
if (i+1 < H)
dx = data[i+1][j] - data[i][j];
if (j+1 < W)
dy = data[i][j+1] - data[i][j];
dx = pabs(dx), dy = pabs(dy);
data[i][j] = (dx + dy)/2.0;
}
}
}
void print() {
printf("%d %d\n", W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
}
} test;
int main() {
test.read();
test.border();
test.print();
return 0;
}
Read More +

b427. 漸層色彩

Problem

給予兩個顏色的像素,在一張空白影像中使用漸層效果。

  • 水平由左到右的漸層
  • 從中心點擴散的漸層

Sample Input

1
3 3 0 0 0 0 255 255 255 255 255

Sample Output

1
2
3
4
3 3
0 0 0 255 128 128 128 255 255 255 255 255
0 0 0 255 128 128 128 255 255 255 255 255
0 0 0 255 128 128 128 255 255 255 255 255

Solution

線性內插,特別小心中心點 (x, y) 會是一個浮點數情況。

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
#include <bits/stdc++.h>
using namespace std;
class IMAGE {
public:
struct Pixel {
double r, g, b, a;
Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
r(x), g(y), b(z), a(w) {}
void read() {
scanf("%lf %lf %lf %lf", &r, &g, &b, &a);
}
Pixel operator-(const Pixel &x) const {
return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
}
Pixel operator+(const Pixel &x) const {
return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
}
Pixel operator*(const double x) const {
return Pixel(r*x, g*x, b*x, a*x);
}
Pixel operator/(const double x) const {
return Pixel(r/x, g/x, b/x, a/x);
}
void print() {
printf("%d %d %d %d", (int)round(r), (int)round(g), (int)round(b), (int)round(a));
}
};
int W, H;
static const int MAXN = 300;
Pixel data[MAXN][MAXN];
void read() {
scanf("%d %d", &W, &H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].read();
}
void alloc() {
int TYPE;
Pixel st, ed;
scanf("%d %d %d", &W, &H, &TYPE);
st.read(), ed.read();
gradient(st, ed, TYPE);
}
void gradient(Pixel st, Pixel ed, int TYPE) {
if (TYPE != 0 && TYPE != 1)
return;
if (TYPE == 0) {
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j] = W-1 ? st + (ed - st) * j / (W-1) : st;
} else {
double cx = (H-1)/2.0, cy = (W-1)/2.0;
double cr = hypot(cx, cy);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j] = cr ? st + (ed - st) * hypot(i-cx, j-cy) / cr : st;
}
}
void print() {
printf("%d %d\n", W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
}
} test;
int main() {
test.alloc();
test.print();
return 0;
}
Read More +

b426. 宇宙光明體

Problem

模擬兩張圖片的合成。

將前景圖片疊在背景圖片之上,按照權重比例來融合重疊部分得像素。

Sample Input

1
2
3
4
5
0 0 0.5
1 1
255 255 255
1 1
0 0 0

Sample Output

1
2
1 1
128 128 128

Solution

原來是宇宙光明體,還以為是宇宙大覺者呢。

按照公式純模擬即可。

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
#include <bits/stdc++.h>
using namespace std;
const double eps = 0;
class IMAGE {
public:
struct Pixel {
long double r, g, b, a;
Pixel(long double x = 0, long double y = 0, long double z = 0, long double w = 0):
r(x), g(y), b(z), a(w) {}
void read() {
int p1, p2, p3;
scanf("%d %d %d", &p1, &p2, &p3);
r = p1, g = p2, b = p3;
}
Pixel operator-(const Pixel &x) const {
return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
}
Pixel operator+(const Pixel &x) const {
return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
}
Pixel operator*(const long double x) const {
return Pixel(r*x, g*x, b*x, a*x);
}
Pixel operator/(const long double x) const {
return Pixel(r/x, g/x, b/x, a/x);
}
void print() {
double p1 = r, p2 = g, p3 = b, p4 = a;
printf("%.0lf %.0lf %.0lf", p1 + eps, p2 + eps, p3 + eps);
}
};
int W, H;
static const int MAXN = 256;
Pixel data[MAXN][MAXN];
void read() {
scanf("%d %d", &W, &H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].read();
}
void add(IMAGE &other, int X, int Y, double OPACITY) {
for (int i = 0; i < other.H && i+X < H; i++) {
for (int j = 0; j < other.W && j+Y < W; j++) {
if (i+X >= 0 && j+Y >= 0)
data[i+X][j+Y] = other.data[i][j] * OPACITY + data[i+X][j+Y] * (1 - OPACITY);
}
}
}
void print() {
printf("%d %d\n", W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
}
} back, fore;
int main() {
double OPACITY;
int X, Y;
scanf("%d %d %lf", &X, &Y, &OPACITY);
fore.read();
back.read();
back.add(fore, Y, X, OPACITY);
back.print();
return 0;
}
Read More +

b425. 抽菸動作 請勿模仿

Problem

替影像打上馬賽克效果,馬賽克效果是把區域內部的每一個像素改變,改變的方法為從舊的圖像中,搜索以其為中心的 11 x 11 方形內部所有像素的平均。

這一題要馬賽克區域為圓形,接著找到方形內部的像素平均即可。

Sample Input

1
2
3
4
0 1 2
1 2
1 2 3
4 5 6

Sample Output

1
2
3
1 2
3 4 5
3 4 5

Solution

馬賽克的方案也許還有其他,從這裡可以知道平均造就出新的像素顏色。那麼可以列出方程式,所謂的解馬賽克估計也是去找線性聯立方程組,利用高斯消去法就能求解,可惜的是會有多組解的情況吧。

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
#include <bits/stdc++.h>
using namespace std;
class IMAGE {
public:
struct Pixel {
double r, g, b, a;
Pixel(double x = 0, double y = 0, double z = 0, double w = 0):
r(x), g(y), b(z), a(w) {}
void read() {
scanf("%lf %lf %lf", &r, &g, &b);
}
Pixel operator-(const Pixel &x) const {
return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
}
Pixel operator+(const Pixel &x) const {
return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
}
Pixel operator*(const double x) const {
return Pixel(r*x, g*x, b*x, a*x);
}
Pixel operator/(const double x) const {
return Pixel(r/x, g/x, b/x, a/x);
}
void print() {
printf("%d %d %d", (int)round(r), (int)round(g), (int)round(b));
}
};
int W, H;
static const int MAXN = 300;
Pixel data[MAXN][MAXN], tmp[MAXN][MAXN];
void read() {
scanf("%d %d", &W, &H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].read();
}
void blur(long long X, long long Y, long long R) {
if (R < 0) return ;
for (int i = 0; i < H; i++) {
for (int j = 0; j < W; j++) {
if ((i-X)*(i-X) + (j-Y)*(j-Y) <= R*R) {
Pixel B(0, 0, 0);
int cnt = 0;
for (int a = -5; a <= 5; a++) {
for (int b = -5; b <= 5; b++) {
if (i+a < H && j+b < W && i+a >= 0 && j+b >= 0) {
B = B + data[i+a][j+b];
cnt++;
}
}
}
tmp[i][j] = B / cnt;
} else {
tmp[i][j] = data[i][j];
}
}
}
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j] = tmp[i][j];
}
void print() {
printf("%d %d\n", W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
}
} test;
int main() {
long long X, Y, R;
scanf("%lld %lld %lld", &X, &Y, &R);
test.read();
test.blur(Y, X, R);
test.print();
return 0;
}
Read More +

b424. 圖片縮放

Problem

此題有精準度誤差,未能提交成功。

進行圖片縮放,有兩種方案來補足缺失或者是簡化的像素計算。不管縮放如何,在新的坐標系的點可以對應到舊坐標系中,並且會有四個像素點 (正方形) 包含這個對應點,有可能會有所缺失,但至少包含一個點。

  • 最近鄰點法 (nearest-neighbor approach)
    找正方形中哪個點最靠近新的對應點,就令其像素顏色等於最靠近的點。
  • 雙線性內插 (bi-linear interpolation)
    顏色由相鄰四個點決定,使用雙性內插來得到新的像素顏色。

特別注意,雖然邊長根據縮放四捨五入,但在計算時先不考慮最終結果的長寬,否則縮放倍率也會跟著四捨五入。

可以參考 中央大學影像處理實驗室 課程 講義第二章

Sample Input

1
2
3
2 2 0
2 1
1 2 3 255 4 5 6 255

Sample Output

1
2
3
4 2
1 2 3 255 3 4 5 255 4 5 6 255 4 5 6 255
1 2 3 255 3 4 5 255 4 5 6 255 4 5 6 255

Solution

除了陷阱的邊長四捨五入外,精準度根據雙性內插影響。

雙線性內插可以由下列兩種方案,或者還有其他的

  • 兩次 lerp 插值,先水平後垂直或者先垂直後水平
  • 特別的 $f(x, y) = (1-a)(1-b) g(j, k) + a(1-b) g(j+1, k) + (1-a) b g(j, k+1) + a b g(j+1, k+1)$

特別發現到,第二種方案沒有用到除法運算,但乘法運算的優先順序會影響到誤差。然而我的代碼是混用造成難以評估的誤差。關於第二種方案的坐標系,可以參考講義中的圖示。

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
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-18;
class IMAGE {
public:
struct Pixel {
long double r, g, b, a;
Pixel(long double x = 0, long double y = 0, long double z = 0, long double w = 0):
r(x), g(y), b(z), a(w) {}
void read() {
int p1, p2, p3, p4;
scanf("%d %d %d %d", &p1, &p2, &p3, &p4);
r = p1, g = p2, b = p3, a = p4;
}
Pixel operator-(const Pixel &x) const {
return Pixel(r-x.r, g-x.g, b-x.b, a-x.a);
}
Pixel operator+(const Pixel &x) const {
return Pixel(r+x.r, g+x.g, b+x.b, a+x.a);
}
Pixel operator*(const long double x) const {
return Pixel(r*x, g*x, b*x, a*x);
}
Pixel operator/(const long double x) const {
return Pixel(r/x, g/x, b/x, a/x);
}
void print() {
double p1 = r, p2 = g, p3 = b, p4 = a;
printf("%.0lf %.0lf %.0lf %.0lf", p1 + eps, p2 + eps, p3 + eps, p4 + eps);
}
};
int W, H;
static const int MAXN = 512;
Pixel data[MAXN][MAXN], tmp[MAXN][MAXN];
void read() {
scanf("%d %d", &W, &H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].read();
}
int isValid(int x, int y) {
return x >= 0 && y >= 0 && x < H && y < W;
}
void resize(long double SX, long double SY, int trans = 0) {
if (trans != 0 && trans != 1)
return ;
int rW = (int) round(SX * W);
int rH = (int) round(SY * H);
if (rW > 256 || rH > 256 || rW <= 0 || rH <= 0)
return ;
for (int i = 0; i < rH; i++) {
for (int j = 0; j < rW; j++) {
long double x = i / SY;
long double y = j / SX;
int lx, rx, ly, ry;
lx = floor(x), rx = ceil(x);
ly = floor(y), ry = ceil(y);
Pixel v[] = {data[lx][ly], data[lx][ry], data[rx][ly], data[rx][ry]};
int px[] = {lx, lx, rx, rx};
int py[] = {ly, ry, ly, ry};
if (trans == 0) {
Pixel t1, t2;
if (isValid(px[0], py[0]) && isValid(px[1], py[1]) && isValid(px[2], py[2]) && isValid(px[3], py[3])) {
long double a = x - lx, b = y - ly;
tmp[i][j] = v[0] * (1 - a) * (1 - b) + v[2] * a * (1 - b)
+ v[1] * (1 - a) * b + v[3] * a * b;
} else if (isValid(px[0], py[0]) && isValid(px[2], py[2])) {
if (rx != lx)
tmp[i][j] = v[0] + (v[2] - v[0]) * ((x - lx) / (rx - lx));
else
tmp[i][j] = v[0];
} else if (isValid(px[0], py[0]) && isValid(px[1], py[1])) {
if (ry != ly)
tmp[i][j] = v[0] + (v[1] - v[0]) * ((y - ly) / (ry - ly));
else
tmp[i][j] = v[0];
} else if (isValid(px[0], py[0])) {
tmp[i][j] = v[0];
} else {
puts("WTF");
}
} else {
int c = -1;
double mndist;
for (int k = 0; k < 4; k++) {
if (!isValid(px[k], py[k]))
continue;
double d = (x-px[k])*(x-px[k])+(y-py[k])*(y-py[k]);
if (c == -1 || mndist > d)
c = k, mndist = d;
}
tmp[i][j] = data[px[c]][py[c]];
}
}
}
W = rW, H = rH;
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j] = tmp[i][j];
}
void print() {
printf("%d %d\n", W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++)
data[i][j].print(), printf("%c", j == W-1 ? '\n' : ' ');
}
} test;
int main() {
double X, Y;
int TYPE;
scanf("%lf %lf %d", &X, &Y, &TYPE);
test.read();
test.resize(X, Y, TYPE);
test.print();
return 0;
}
Read More +

b423. 魔術橡皮擦

Problem

提供影像去背、相鄰去色兩種方案。

  • 去背對整張影像找像素差小於某一個範圍內,就將其 alpha 值設成 0。
  • 相鄰去色對鄰近四個方向上下左右進行搜索,將像素差低於某個閥值的都加入到集合中進行 alpha 值設成 0。

Sample Input

1
2
3
4
0 0 0 0
3 2
1 2 3 255 4 5 6 255 7 8 9 255
10 11 12 255 13 14 15 255 16 17 18 255

Sample Output

1
2
3
3 2
1 2 3 0 4 5 6 255 7 8 9 255
10 11 12 255 13 14 15 255 16 17 18 255

Solution

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
#include <bits/stdc++.h>
using namespace std;
int X, Y, TYPE, DIFF;
namespace Image {
const int MAXN = 256;
int n, m, rgba[MAXN][MAXN][4], used[MAXN][MAXN] = {};
int RGB_diff(int x, int y, int a, int b) {
return (rgba[x][y][0]-rgba[a][b][0])*(rgba[x][y][0]-rgba[a][b][0])+
(rgba[x][y][1]-rgba[a][b][1])*(rgba[x][y][1]-rgba[a][b][1])+
(rgba[x][y][2]-rgba[a][b][2])*(rgba[x][y][2]-rgba[a][b][2]);
}
void gTransparent(int x, int y, long long diff) {
if (x < 0 || y < 0 || x >= m || y >= n || diff < 0)
return ;
diff = diff * diff;
for (int i = 0; i < m; i++)
for (int j = 0; j < n; j++)
if (RGB_diff(x, y, i, j) <= diff)
rgba[i][j][3] = 0;
}
void lTransparent(int x, int y, int ox, int oy, long long diff) {
if (ox < 0 || oy < 0 || ox >= m || oy >= n || diff < 0)
return ;
static const int dx[] = {0, 0, 1, -1};
static const int dy[] = {1, -1, 0, 0};
if (x < 0 || y < 0 || x >= m || y >= n)
return ;
if (RGB_diff(x, y, ox, oy) > diff * diff || used[x][y])
return ;
rgba[x][y][3] = 0, used[x][y] = 1;
for (int i = 0; i < 4; i++)
lTransparent(x+dx[i], y+dy[i], ox, oy, diff);
}
void print() {
printf("%d %d\n", n, m);
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++)
printf("%d %d %d %d%c", rgba[i][j][0], rgba[i][j][1], rgba[i][j][2], rgba[i][j][3], j == n-1 ? '\n' : ' ');
}
}
}
int main() {
scanf("%d %d %d %d", &X, &Y, &TYPE, &DIFF);
scanf("%d %d", &Image::n, &Image::m);
for (int i = 0; i < Image::m; i++)
for (int j = 0; j < Image::n; j++)
for (int k = 0; k < 4; k++)
scanf("%d", &Image::rgba[i][j][k]);
swap(X, Y);
if (TYPE == 0)
Image::gTransparent(X, Y, DIFF);
if (TYPE == 1)
Image::lTransparent(X, Y, X, Y, DIFF);
Image::print();
return 0;
}
Read More +

b422. Colorful Life and Monochromatic Life

Problem

將一張 rgb 表示的彩色影像變成灰階影像。

套用公式 round((r + g + b)/3) 得到新的像素色彩。

Sample Input

1
2
3
3 2
1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18

Sample Output

1
2
3
3 2
2 2 2 5 5 5 8 8 8
11 11 11 14 14 14 17 17 17

Solution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <bits/stdc++.h>
using namespace std;
int main() {
int n, m;
while (scanf("%d %d", &n, &m) == 2) {
printf("%d %d\n", n, m);
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
int r, g, b, w;
scanf("%d %d %d", &r, &g, &b);
w = (int) round((r + g + b)/3.0);
printf("%d %d %d%c", w, w, w, j == n-1 ? '\n' : ' ');
}
}
}
return 0;
}
Read More +