解題區/解題區 - Zerojudge

b442. 快取實驗 矩陣乘法

Problem

背景

記得 b439: 快取置換機制 提到的快取置換機制嗎？現在來一場實驗吧！

題目描述

相信不少人都已經實作所謂的矩陣乘法，計算兩個方陣大小為 $N \times N$ 的矩陣 $A, B$。為了方便起見，提供一個偽隨機數的生成，減少在輸入處理浪費的時間，同時也減少上傳測資的辛苦。

根據種子 $c = S1$ 生成矩陣 $A$，種子 $c = S2$ 生成矩陣 $B$，計算矩陣相乘 $A \times B$，為了方便起見，使用 hash 函數進行簽章，最後輸出一個值。由於會牽涉到 overflow 問題，此題作為快取實驗就不考慮這個，overflow 問題大家都會相同。

請利用快取優勢修改代碼如下：

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
#include <bits/stdc++.h>
using namespace std;
class Matrix {
public:
    vector< vector<int> > data;
    int row, col;
    Matrix(int n = 1, int m = 1) {
        data = vector< vector<int> >(n, vector<int>(m, 0));
        row = n, col = m;
    }
    void rand_gen(int c) {
        int x = 2, n = row * col;
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                x = ((long long) x * x + c + i + j)%n;
                data[i][j] = x;
            }
        }
    }
    Matrix multiply(Matrix x) {
        Matrix ret(row, x.col);
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < x.col; j++) {
                int sum = 0;	// overflow
                for (int k = 0; k < col; k++)
                    sum += data[i][k] * x.data[k][j];
                ret.data[i][j] = sum;
            }
        }
        return ret;
    }
    void print() {
        for (int i = 0; i < row; i++) {
            printf("[");
            for (int j = 0; j < col; j++) {
                printf(" %d", data[i][j]);
            }
            printf(" ]\n");
        }
    }
    unsigned long signature() {
        unsigned long h = 0;
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                h = hash(h + data[i][j]);
            }
        }
        return h;
    }
private:
    unsigned long hash(unsigned long x) {
        return (x * 2654435761LU);
    }
};
int main() {
    int N, S1, S2;
    while (scanf("%d %d %d", &N, &S1, &S2) == 3) {
        Matrix A(N, N), B(N, N), C;
        A.rand_gen(S1);
        B.rand_gen(S2);
        C = A.multiply(B);		
//		A.print();
//		B.print();
//		C.print();
        printf("%lu\n", C.signature());
    }
    return 0;
}

Sample Input

1
2
2 2 5
2 2 5

Sample Output

1
2
573770929
573770929

Solution

在越大的矩陣中，快取記憶體置換是很慢的，即便用 L1, L2, main memory 分層快取，速度差異很明顯，盡可能在線性 $O(n)$ 運算上都不造成記憶體置換 (搬移)，雖然複雜度都是 $O(n^3)$，相信沒人會去實作 $O(n^{2.807})$ 的 Strassen 算法，又或者是現今更快的 Coppersmith-Winograd $O(n^{2.3727})$。

作業系統題目第二彈，考驗快取應用。在 ZJ 主機上速度可以快二十倍，在本機上只快了十倍。

  1. 蔡神 asas 修改輸入，而我選擇了轉置，速度落後。
  2. 柳州 liouzhou_101 同步簽章，而我選擇分開函數，速度落後。
  3. 廖氏如神 pcsh710742 下刀計組，等等，我看見了什麼。

現在已經快了四十倍，就只是因為編譯器的優化參數變成手動。詳細實作和效能分析，需要的知識不只是作業系統，同時涵蓋計算機組織的 CPU stall 問題，參考 Optimizing Large Matrix-Vector Multiplications

先來個一般解，在計算矩陣 $A \times B$ 前，先將 $B$ 轉置，接著修改計算方式 

1
2
for (int k = 0; k < col; k++) 
    sum += data[i][k] * x.data[j][k];

在這一個迴圈中，陣列採用 row-major 的方式儲存，所以第二維度的區域基本上都在快取中，miss 的機會大幅度降低。

接下來要提案的做法是結合上述三位的做法，代碼快，但不能當 Matrix 模板使用，僅僅做為一個效能體現。共用記憶體，採用指針的方式減少複製，一樣進行轉置，但這會破壞原有的 $B$ 矩陣配置。同步進行簽章，捨棄掉 $C$ 矩陣的儲存，使用 LOOP UNROLLING，由於 Zerojudge 的優化沒有開到 -O2、-O3、-Ofast，關於這一部分的優化由使用者自己來。

LOOP UNROLLING 的優化在於 branch 指令，由於 pipeline 會讓效能提升，但遇到 branch 時必須捨棄掉偷偷載入的指令、算到一半的結果，效能會下降，因此使用 LOOP UNROLLING 的好處在於減少 branch 次數。

關於 data prefetch 可以用 __builtin_prefetch() 來完成，根據廖氏如神所言，這個概念可以預先載入防止 stall (pipeline hazard) 的拖延，速度並不會提升太多，有可能是硬體已經完成這一部分，甚至用別的架構去克服。

轉置解

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
#include <bits/stdc++.h>
using namespace std;
class Matrix {
public:
    vector< vector<int> > data;
    int row, col;
    Matrix(int n = 1, int m = 1) {
        data = vector< vector<int> >(n, vector<int>(m, 0));
        row = n, col = m;
    }
    void rand_gen(int c) {
        int x = 2, n = row * col;
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                x = ((long long) x * x + c + i + j)%n;
                data[i][j] = x;
            }
        }
    }
    Matrix multiply(Matrix &x) {
        Matrix ret(row, x.col);
        x.transpose();
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < x.col; j++) {
                int sum = 0;	// overflow
                int *a = &data[i][0], *b = &x.data[j][0];
                for (int k = 0; k < col; k++)
                    sum += *a * *b, a++, b++;
                ret.data[i][j] = sum;
            }
        }
        x.transpose();
        return ret;
    }
    void transpose() {
        for (int i = 0; i < row; i++)
        	for (int j = i+1; j < col; j++)
            	swap(data[i][j], data[j][i]);
    }
    void print() {
        for (int i = 0; i < row; i++) {
            printf("[");
            for (int j = 0; j < col; j++)
                printf(" %d", data[i][j]);
            printf(" ]\n");
        }
    }
    unsigned long signature() {
        unsigned long h = 0;
        for (int i = 0; i < row; i++)
            for (int j = 0; j < col; j++)
                h = hash(h + data[i][j]);
        return h;
    }
private:
    inline unsigned long hash(unsigned long x) {
        return (x * 2654435761LU);
    }
};
int main() {
    int N, S1, S2;
    while (scanf("%d %d %d", &N, &S1, &S2) == 3) {
        Matrix A(N, N), B(N, N), C;
        A.rand_gen(S1);
        B.rand_gen(S2);
        C = A.multiply(B);
        printf("%lu\n", C.signature());
    }
    return 0;
}

LOOP UNROLL 解

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
#include <bits/stdc++.h>
using namespace std;
#define LOOP_UNROLL 8
class Matrix {
public:
    vector< vector<int> > data;
    int row, col;
    Matrix(int n = 1, int m = 1) {
        row = n, col = m;
        data = vector< vector<int> >(n, vector<int>(m, 0));
    }
    void rand_gen(int c) {
        int x = 2, n = row * col;
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < col; j++) {
                x = ((long long) x * x + c + i + j)%n;
                data[i][j] = x;
            }
        }
    }
    unsigned long multiply(Matrix &x) {
        x.transpose();
        unsigned long h = 0;
        for (int i = 0; i < row; i++) {
            for (int j = 0; j < x.col; j++) {
                register int sum = 0;
                int *a = &data[i][0], *b = &x.data[j][0], k;
                for (k = 0; k+LOOP_UNROLL < col; k += LOOP_UNROLL) {
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                    sum += *a * *b, a++, b++;
                }
                for (; k < col; k++)
                    sum += *a * *b, a++, b++;
                h = hash(h + sum);
            }
        }
        return h;
    }
    void transpose() {
        for (int i = 0; i < row; i++)
        	for (int j = i+1; j < col; j++)
            	swap(data[i][j], data[j][i]);
    }
    void print() {
        for (int i = 0; i < row; i++) {
            printf("[");
            for (int j = 0; j < col; j++)
                printf(" %d", data[i][j]);
            printf(" ]\n");
        }
    }
private:
    unsigned long hash(unsigned long x) {
        return (x * 2654435761LU);
    }
} A(1000, 1000), B(1000, 1000);
int main() {
    int N, S1, S2;
    while (scanf("%d %d %d", &N, &S1, &S2) == 3) {
        A.row = A.col = B.row = B.col = N;
        A.rand_gen(S1);
        B.rand_gen(S2);
        printf("%lu\n", A.multiply(B));
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b439. 快取置換機制

Problem

背景

編寫程式不僅僅是在數學分析上達到高效率，儘管數學分析的複雜度不是最好，理解電腦運作模式也能有效地讓程式變快。例如 資料結構 Unrolled linked list (常翻譯成 塊狀鏈表) 便是利用此一優勢，讓速度顯著地提升，之所以能追上不少常數大的平衡樹操作運用的技巧就是快取效能改善。

講一個更簡單的例子，宣告整數陣列的兩個方案：

方案一

1
const int LARGE_SIZE = 1<<30;  int A[LARGE_SIZE], B[LARGE_SIZE];

方案二

1
const int LARGE_SIZE = 1<<30; struct DATA { int A, B; } E[LARGE_SIZE];

演算法的複雜度倘若一樣，若在 A, B 相當高頻率的替換，則快取操作必須不斷地將內容置換，若 A, B 在算法中是獨立運算，則方案一的寫法會來得更好，反之取用方案二會比較好。最常見的運作可以在軟體模擬矩陣乘法中見到，預先將矩陣轉置，利用快取優勢速度可以快上 8 倍以上。

問題描述

給予一個記憶體空間大小為 $M$，使用 Least Recently Used (LRU) 策略進行置換，LRU 的策略為將最久沒有被使用過的空間替換掉，也就是需要從硬碟讀取 $disk[i]$ 到 $memory$ 時，發現記憶體都已經用光，則把不常用的 $mem[j]$ 寫入 $disk[k]$，再將 $disk[i]$ 內容寫入 $mem[j]$。

下圖是一個 $M = 4$ 簡單的範例：

1
2
3
4
5
6
7
8
9
10
       +---+     +---+     +---+     +---+     +---+     +---+     +---+     +---+
mem[0] | 1 |     | 1 |     | 1 |     | 1 |     | 1 |     | 1 |     | 1 |     | 1 |
       +---+     +---+     +---+     +---+     +---+     +---+     +---+     +---+
mem[1] |   |     | 2 |     | 2 |     | 2 |     | 2 |     | 2 |     | 2 |     | 2 |
       +---+ +-> +---+ +-> +---+ +-> +---+ +-> +---+ +-> +---+ +-> +---+ +-> +---+
mem[2] |   |     |   |     | 3 |     | 3 |     | 3 |     | 3 |     | 5 |     | 5 |
       +---+     +---+     +---+     +---+     +---+     +---+     +---+     +---+
mem[3] |   |     |   |     |   |     | 4 |     | 4 |     | 4 |     | 4 |     | 3 |
       +---+     +---+     +---+     +---+     +---+     +---+     +---+     +---+ 
         1         2         3         4         1         2         5         3

依序使用 1, 2, 3, 4, 1, 2, 5, 3 的配置情況，特別是在 5 使用的時候，會將記憶體中上一次最晚使用的 3 替換掉。

現在寫程式去模擬置換情況。

Sample Input

1
2
3
4
5
6
7
8
9
10
11
12
4
0 1
1 2 514
1 3 101
1 4 50216
0 1
0 2
1 5 6
0 3
0 5
0 2

Sample Output

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

Solution

軟體仿作，使用一個 hash 和一個 list 進行維護，而不是使用 priority queue 來維護，用 list 取代之。當進行取址時，順道把對應 list 指針移動，所有步驟期望複雜度都是在 $O(1)$ 完成，若使用 priority queue 會掉到 $O(\log n)$ 去進行更新。

換出作業系統題目，來一個最常見到的快取 LRU 算法。這個問題在 Leetcode OJ 上面也有，公司面試有機會遇到。

但實作時被自己坑，比較柳柳州給予的測試，速度居然慢個二十多倍。原因是這樣子的

1
2
map[key] = value;
return map.find(key);

替換成以下寫法，速度就起飛了。

1
return map.insert({key, value}).first;

也就是說插入的時候順便把指針拿到，避免第二次搜索。

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
#include <bits/stdc++.h>
using namespace std;
typedef pair<int, list<int>::iterator> PIT_TYPE;
typedef unordered_map<int, pair<int, list<int>::iterator> >::iterator HIT_TYPE;
class LRUCache {
public:	
    int memIdx, size;
    list<int> time;
    vector<int> mem;
    unordered_map<int, PIT_TYPE> addr;
    LRUCache(int capacity) {
    	size = capacity;
        mem.resize(capacity, 0);
        addr.clear(), time.clear();
        memIdx = 0;
    }
    
    pair<int, int> get(int key) {
    	it = addr.find(key);
        if (it == addr.end())
            it = replace(key), mem[it->second.first] = 0;
        it->second.second = recent(it->second.second);
        return make_pair(it->second.first, mem[it->second.first]);
    }
    
    HIT_TYPE set(int key, int value) {
        it = addr.find(key);
        if (it == addr.end())
        	it = replace(key);
        it->second.second = recent(it->second.second);
        mem[it->second.first] = value;
        return it;
    }
private:
    HIT_TYPE it;
    HIT_TYPE replace(int key) {
        int mpos = -1, trash;
        list<int>::iterator lpos;
        if (addr.size() == size) {
            trash = time.front(), time.pop_front();
            it = addr.find(trash);
            mpos = it->second.first;
            addr.erase(it);
        } else {
            mpos = memIdx++;
        }
        lpos = time.insert(time.end(), key);
        return addr.insert({key, make_pair(mpos, lpos)}).first;
    }
    list<int>::iterator recent(list<int>::iterator p) {
        int key = *p;
        time.erase(p);
        return time.insert(time.end(), key);
    }
};
int main() {
    int M, cmd, x, y;
    pair<int, int> t;
    scanf("%d", &M);
    
    LRUCache LL(M);
    while (scanf("%d", &cmd) == 1) {
        if (cmd == 0) {
            scanf("%d", &x);
            pair<int, int> t = LL.get(x);
            printf("%d %d\n", t.first, t.second);
        } else {
            scanf("%d %d", &x, &y);
            LL.set(x, y);
        }
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b435. 尋找原根

Problem

給定一個模數 $m$，找到其 primitive root $g$。

primitive root $g$ 滿足

$$g^{\phi(m)} \equiv 1 \mod m \\ g^{\gamma} \not\equiv 1 \mod m \text{ , for } 1 \le \gamma < \phi(m)$$

意即在模 $m$ 下，循環節 $ord_m(g) = \phi(m)$，若 $m$ 是質數，則 $\phi(m) = m-1$，也就是說循環節長度是 $m-1$，更可怕的是組成一個 $[1, m-1]$ 的數字排列。

現在假定 $m$ 是質數，請求出一個最小的 primitive root $g$。

Sample Input

1
2
3
4
2
3
5
7

Sample Output

1
2
3
4
1
2
2
3

Solution

primitive root 在密碼學中，用在 Diffie-Hellman 的選定中的穩定度，倘若基底是一個 primitive root 將會提升破解的難度，因為對應情況大幅增加 (值域比較廣，因為不循環)。

由於歐拉定理，可以得知 $a^{\phi(p)} \equiv 1 \mod p$，那麼可以堆導得到若 $a^x \equiv 1 \mod p$，則滿足 $x | \phi(p)$，否則與歐拉矛盾，藉此可以作為一個篩選的檢查手法，對於所有 $phi(p)$ 的因數都進行檢查，最後我們可以得到 一般檢測 的作法。

為了加速驗證，由於 $x | \phi(p)$，若最小的 $x$ 存在 $a^x \equiv 1 \mod p$，那麼 $kx$ 也會符合等式，因此只需要檢查 $x = (p-1)/\text{prime-factor}$ 的情況，如此一來速度就快上非常多。最後得到 加速版本 ，此做法由 liouzhou_101 提供想法。

一般檢測

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 <bits/stdc++.h>
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[5500], Pt = 0;
void sieve() {
    register int i, j, k;
    SET(1);
    int n = 45825;
    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;
        }
    }
}
vector< pair<int, int> > factor(int n) {
    vector< pair<int, int> > R;
    
    for(int i = 0, j; i < Pt && P[i] * P[i] <= n; i++) {
        if(n%P[i] == 0) {
            for(j = 0; n%P[i] == 0; n /= P[i], j++);
            R.push_back(make_pair(P[i], j));
        }
    }
    if(n != 1)	R.push_back(make_pair(n, 1));
    return R;
}
void gen_factor(int idx, int x, vector< pair<int, int> > &f, vector<int> &g) {
    if (idx == f.size()) {
        g.push_back(x);
        return ;
    }
    for (int i = 0, a = 1; i <= f[idx].second; i++, a *= f[idx].first) 
        gen_factor(idx+1, x*a, f, g);
}
long long mpow(long long x, long long y, long long mod) {
    long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
int primitive_root(int p) {
    if (p == 2)	return 1;
    vector< pair<int, int> > f = factor(p-1);
    vector<int> g;
    gen_factor(0, 1, f, g);
    g.erase(g.begin()), g.pop_back();	// remove 1, p-1
    random_shuffle(g.begin(), g.end());
    for (int i = 2; i <= p; i++) {
        int ok = 1;
        for (auto &x: g) {
            if (mpow(i, x, p) == 1) {
                ok = 0;
                break;
            }
        }
        if (ok)	return i;
    }
    return -1;
}
int main() {
    sieve();
    int p;
    while (scanf("%d", &p) == 1) 
        printf("%d\n", primitive_root(p));
    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
#include <bits/stdc++.h>
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[5500], Pt = 0;
void sieve() {
    register int i, j, k;
    SET(1);
    int n = 45825;
    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;
        }
    }
}
vector< pair<int, int> > factor(int n) {
    vector< pair<int, int> > R;
    
    for(int i = 0, j; i < Pt && P[i] * P[i] <= n; i++) {
        if(n%P[i] == 0) {
            for(j = 0; n%P[i] == 0; n /= P[i], j++);
            R.push_back(make_pair(P[i], j));
        }
    }
    if(n != 1)	R.push_back(make_pair(n, 1));
    return R;
}
long long mpow(long long x, long long y, long long mod) {
    long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
int primitive_root(int p) {
    if (p == 2)	return 1;
    vector< pair<int, int> > f = factor(p-1);
    for (int i = 2; i <= p; i++) {
        int ok = 1;
        for (auto &x: f) {
            if (mpow(i, (p-1)/x.first, p) == 1) {
                ok = 0;
                break;
            }
        }
        if (ok)	return i;
    }
    return -1;
}
int main() {
    sieve();
    int p;
    while (scanf("%d", &p) == 1) 
        printf("%d\n", primitive_root(p));
    return 0;
}
Read More +
解題區/解題區 - UVa

UVa 11476 - Factorizing Larget Integers

Problem

給你一個正整數 $N$ ($1 < N \le10^{18}$)，請你把 $N$ 質因數分解。

注意：大整數分解

Sample Input

1
2
3
4
3
60
36
10007

Sample Output

1
2
3
60 = 2^2 * 3 * 5
36 = 2^2 * 3^2
10007 = 10007

Solution

2015/07/11 第二版，加速三倍，加快模乘法運算，減少模數利用減法。

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
155
#include <bits/stdc++.h> 
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 || !(p&1))	return 0;
    if (p == 2)				return 1;
    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) {
    long long x = 2, y = 2, i = 1, k = 2, d;
    while (true) {
        x = (mul(x, x, n) + c);
        if (x >= n)	x -= n;
        d = llgcd(x - y, n);
        if (d > 1) return d;
        if (++i == k) y = x, k <<= 1;
    }
    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) {
    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 (int i = 2; d == n; i++)
        d = pollard_rho(n, i);
    llfactorize(d, f);
    llfactorize(n/d, f);
}
int main() {
    sieve();
    int testcase;
    scanf("%d", &testcase);
    while (testcase--) {
        long long n;
        scanf("%lld", &n);
        
        vector<long long> f;
        map<long long, int> r;
            
        llfactorize(n, f);
        for (auto &x : f)
            r[x]++;
            
        printf("%lld =", n);
        for (auto it = r.begin(); it != r.end(); it++) {
            if (it != r.begin())
                printf(" *");
            printf(" %lld", it->first);
            if (it->second > 1)
                printf("^%d", it->second);
        }
        puts("");
        
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b432. 趣味加分

Problem

背景

廖氏如神 (pcsh710742) 在作業中遇到一題，用手算會算到天昏地暗，而問題是這樣子的：

$a^{13337} \equiv n \mod 2^{64}$

其中 $n= 2015 + 2 \times (\text{The last 2 digit of your student ID number})$，請找到 $a < 2^{64}$ 的其中一組解。

問題描述

模數 $2^{64}$ 看起來不大，但對於 Ghz 為 CPU 運算速度單位的電腦而言還是要跑非常久的。因此將問題簡化：

$a^{23333} \equiv n \mod 2^{20}$

現在給予一個 $n$，請求出一組 $a$，測資中保證答案唯一。

Sample Input

1
2
3
4
5
268275
888817
89215
63495
976477

Sample Output

1
2
3
4
5
387
817
639
487
909

Solution

除了偶數無解外，奇數都至少有一個解。而這一題的題目數據恰好奇數都只有一解，那麼就不必處理多組解或者是字典順序最小的，只要專心找到符合的解即可。

  • 暴力建表 $O(2^{20})$ 建完，之後直接查找。
  • 篩選正解，依序窮舉從最低位到最高位 (二進制下) 為 0 還是 1，由於次方會不斷地推移，高位結果不影響低位的對應。窮舉時保留低位符合的解，並且不斷地篩選掉不可能的解方案，複雜度 $O(20 k)$，$k$ 是難以估計的數字。
  • 快速假解，類似篩選正解的做法，但只保留其中一組解進行，複雜度 $O(20)$。這個解法是有點毛病的，但目前找不到反例。

快速假解

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
#include <bits/stdc++.h> 
using namespace std;
long long mul(long long a, long long b, 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; 
}
unsigned long long mpow(unsigned long long x, unsigned long long y, unsigned long long mod) {
    unsigned long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
// find a^23333 = n \mod 2^32
int main() {
    const long long M = 1LL<<20;
    const long long E = 23333;
    long long n;
    while (scanf("%lld", &n) == 1) {
        long long a = 0;
        for (int i = 0; i < 32; i++) {
            long long t = mpow(a, E, M), mask = (1LL<<(i+1))-1;
            if ((t&mask) == (n&mask)) {
                
            } else {
                a |= 1LL<<i;
            }			
        }
        printf("%lld\n", a);
    }
    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
#include <bits/stdc++.h> 
using namespace std;
unsigned long long mpow(unsigned long long x, unsigned long long y, unsigned long long mod) {
    unsigned long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
// find a^23333 = n \mod 2^32
int main() {
    const long long M = 1LL<<20;
    const long long E = 23333;
    long long n;
    while (scanf("%lld", &n) == 1) {
        vector<long long> A(1, 0);
        for (int i = 0; i < 20; i++) {
            vector<long long> nA;
            for (auto &a : A) {
                long long t, mask;
                t = mpow(a, E, M), mask = (1LL<<(i+1))-1;
                if ((t&mask) == (n&mask)) {
                    nA.push_back(a);
                } 
                t = mpow(a|(1LL<<i), E, M), mask = (1LL<<(i+1))-1;
                if ((t&mask) == (n&mask)) {
                    nA.push_back(a|(1LL<<i));
                }			
            }
            A = nA;
        }
        assert(A.size() > 0 && mpow(A[0], E, M) == n);
        printf("%lld\n", A[0]);
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b431. 中國餘數

Problem

背景

中國餘數定理 (Chinese Remainder Theorem，簡稱 CRT) 經常是工程學裡面常用的一種轉換域，很多人不知道當初在大學離散數學中學這個做什麼，但是在不少計算的設計都會運用到 CRT。由於電腦 CPU 架構中的運算單位是 32-bits 或者 64-bits (也許在未來會更長)，但值域高達 128-bits 或者 512-bits 以上模擬運算成了麻煩之處。

回顧中國餘數定理 CRT

$$(S): \left\{\begin{matrix} x \equiv a_1 \mod m_1 \\ x \equiv a_2 \mod m_2 \\ \vdots \\ x \equiv a_n \mod m_n \end{matrix}\right.$$
  • $m_1, m_2, \cdots , m_n$ 任兩數互質，意即 $\forall i \neq j, gcd(m_i, m_j) = 1$
  • 對於任意整數 $a_1, a_2, \cdots , a_n$ 方程組 $(S)$ 均有解，意即找得到一個 $x$ 的參數解。

構造法解 CRT

  1. 設$M = m_1 \times m_2 \times \cdots \times m_n = \prod_{i=1}^{n} m_i$
  2. 設$M_i = M / m_i$
  3. 設$t_i = M_i^{-1} \mod m_i<span>$，意即$</span><!-- Has MathJax -->t_i M_i \equiv 1 \mod m_i$
  4. 方程組 $(S)$ 的通解形式為： $x = a_1 t_1 M_1 + a_2 t_2 M_2 + \cdots + a_n t_n M_n + kM = kM + \sum_{i = 1}^{n} a_i t_i M_i, k \in \mathbb{Z}$
  5. 若限定 $0 \le x < M$，則 $x$ 只有一解。

很多人會納悶通解為什麼長那樣，原因很簡單，要滿足方程組每一條式子，勢必對於$a_i t_i M_i$ 要滿足$x \equiv a_i \mod m_i$ ，因此 $a_i t_i M_i \equiv a_i (t_i M_i) \mod m_i \equiv a_i \mod m_i$ 成立，但是$\forall i \neq j$，滿足$a_i t_i M_i \equiv a_i (t_i M_i) \mod m_j \equiv 0 \mod m_j$

問題描述

來個簡單運用，來計算簡單的 RSA 加解密，特化其中的數學運算。

$M \equiv C^d \mod n$ $n = p \times q$，其中 $p, q$ 是兩個不同的質數，已知 $C, d, p, q$，請求出 $M$。

## Sample Input ##
1
2
88 7 17 11
11 23 17 11

Sample Output

1
2
11
88

Solution

RSA 可以預先處理

  • $c_p = q \times (q^{-1} \mod p)$
  • $c_q = p \times (p^{-1} \mod q)$

還原的算法則是 $M = Mp \times c_p + Mq \times c_q \mod N$

由於拆分後的 bit length 少一半，乘法速度快 4 倍，快速冪次方快 2 倍 (次方的 bit length 少一半)，但是要算 2 次，最後共計快 4 倍。CPU 的乘法想必不會用快速傅立葉 FFT 來達到乘法速度為 $O(n \log n)$。

特別小心「 bit length 少一半 」必須在 $gcd(C, p) = 1$ 時才成立，互質機率機率非常高，仍然有不成立的時候，這情況下速度不是加快 400%。請參照一般解。

例如 C = 27522, d = 17132, p = 2, q = 17293，若使用歐拉定理計算 $Mp = C^{d \mod (p-1)} \mod p = 1$ 事實上 $Mp = C^d \mod p = 0$。有一個特性尚未被利用，對於模質數 $p$ 而言，所有數 $x$ 在模 $p$ 下的循環長度 $L | (p-1)$，最後可以套用 $Mp = C^{d \mod (p-1) + (p-1)} \mod p$ 來完成。請參照循環解，如此一來就不必先判斷 $gcd(C, p) = 1$。

番外

但是利用 CRT 計算容易受到硬體上攻擊，因為會造成 $p, q$ 在分解過程中獨立出現，當初利用 $N$ 很難被分解的特性來達到資訊安全，但是卻因為加速把 $p, q$ 存放道不同時刻的暫存器中。

其中一種攻擊，計算得到 $M = Mp \times q \times (q^{-1} \mod p) + Mq \times p \times (p^{-1} \mod q) \mod N$ 當擾亂後面的式子 (提供不穩定的計算結果)。得到 $M' = Mp \times q \times (q^{-1} \mod p) + (Mq + \Delta) \times p \times (p^{-1} \mod q) \mod N$

接著 $(M' - M) = (\Delta' \times p) \mod N$，若要求 $p$ 的方法為 $gcd(M' - M, N) = gcd(\Delta' \times p, N) = p$，輾轉相除的概念呼之欲出，原來 $p$ 會被這麼夾出，當得到兩個 $p, q$，RSA 算法就會被破解。

一般解

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
#include <bits/stdc++.h> 
using namespace std;
long long mul(long long a, long long b, 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) {
    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) {
    long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
int main() {
    long long C, d, p, q;
    long long N, M, Cp, Cq, Mp, Mq;
    while (scanf("%lld %lld %lld %lld", &C, &d, &p, &q) == 4) {
        N = p * q;
        Mp = mpow(C%p, llgcd(C, p) == 1 ? d%(p-1) : d, p);
        Mq = mpow(C%q, llgcd(C, q) == 1 ? d%(q-1) : d, q);
        Cp = q*inverse(q, p)%N;
        Cq = p*inverse(p, q)%N;
        M = (mul(Mp, Cp, N) + mul(Mq, Cq, N))%N;
        printf("%lld\n", M);
    }
    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
#include <bits/stdc++.h> 
using namespace std;
typedef unsigned long long UINT64;
UINT64 mul(UINT64 a, UINT64 b, UINT64 mod) { 
    UINT64 ret = 0; 
    for (a = a >= mod ? a%mod : a, b = b >= mod ? b%mod : b; 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 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) {
    long long ret = 1;
    while (y) {
        if (y&1) 
            ret = (ret * x)%mod;
        y >>= 1, x = (x * x)%mod;
    }
    return ret % mod;
}
int main() {
    long long C, d, p, q;
    long long N, M, Cp, Cq, Mp, Mq;
    while (scanf("%lld %lld %lld %lld", &C, &d, &p, &q) == 4) {
        N = p * q;
        Mp = mpow(C%p, d%(p-1) + (p-1), p);
        Mq = mpow(C%q, d%(q-1) + (q-1), q);
        Cp = mul(q, inverse(q, p), N);
        Cq = mul(p, inverse(p, q), N);
        M = (mul(Mp, Cp, N) + mul(Mq, Cq, N))%N;
        printf("%lld\n", M);
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b430. 簡單乘法

Problem

背景

在早期密碼世界中， 各種運算都先講求速度，不管是在硬體、軟體利用各種數學定義來設計加密算法就為了加快幾倍的速度，但在近代加密中，加速方法會造成硬體實作攻擊，速度和安全，你選擇哪一個呢。

題目描述

$$a b \equiv x \mod n$$

已知 $a, b, n$，求出 $x$。

Sample Input

1
2
3
4
3 5 7
2 4 3
2 0 2
5 1 4

Sample Output

1
2
3
4
1
2
0
1

Solution

利用加法代替乘法避免溢位，加法取模換成減法加快速度。

由於這一題範圍在在 $10^{18}$，用 long long 型態沒有問題，若在 $2^{63}$ 請替換成 unsigned long long。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#include <bits/stdc++.h> 
using namespace std;
long long mul(long long a, long long b, 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; 
}
int main() {
    long long a, b, n;
    while (scanf("%lld %lld %lld", &a, &b, &n) == 3)
        printf("%lld\n", mul(a, b, n));
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b429. 離散對數

Problem

背景

高中時上過對數，了解 $a^x = b$，則 $x = \log_{a}{b}$。這個問題很簡單，但是 log() 又是怎麼運作，當時是用查表法，不久在大學就會學到泰勒級數，藉由電腦運算，計算越多次就能更逼近真正的結果。

離散對數的形式如下：

$$a^x \equiv b \mod n$$

已知 $a, b, n$，通常會設定 $0 \le a, b < n$。這問題的難處在於要怎麼解出 $x$，沒有 log() 可以這麼迅速得知。

為什麼需要離散對數？不少的近代加密算法的安全強度由這個問題的難度決定，例如 RSA 加密、Diffie-Hellman 金鑰交換 … 等，實際運用需要套用許多數論原理。然而，加密機制要保證解得回來，通常會保證 $gcd(a, n) = 1$，讓乘法反元素 (逆元) 存在。

問題描述

$$a^x \equiv b \mod n$$

已知 $a, b, n$，解出最小的 $x$，若不存在解則輸出 NOT FOUND。

Sample Input

1
2
3
4
2 1 5
2 2 5
3 5 17
4 2 17

Sample Output

1
2
3
4
0
1
5
NOT FOUND

Solution

解決問題 $y = g^x \mod p$，當已知 $y, g, p$，要解出 $x$ 的難度大為提升，不像國高中學的指數計算，可以藉由 log() 運算來完成，離散對數可以的複雜度相當高，當 $p$ 是一個相當大的整數時，通常會取用到 256 bits 以上，複雜度則會在 $O(2^{100})$ 以上。

實際上有一個有趣的算法 body-step, giant-step algorithm，中文翻譯為 小步大步算法 ，在 ACM-ICPC 競賽中也可以找到一些題目來玩玩，算法的時間複雜度為 $O(\sqrt{p})$，空間複雜度也是 $O(\sqrt{p})$。相信除了這個外，還有更好的算法可以完成。

小步大步算法其實很類似塊狀表的算法，分塊處理，每一塊的大小為 $\sqrt{p}$，為了找尋答案計算其中一塊的所有資訊，每一塊就是一小步，接著就是利用數學運算，拉動數線，把這一塊往前推動 (或者反過來把目標搜尋結果相對往塊的地方推動)。因此需要走 $\sqrt{p}$ 大步完成。

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
#include <bits/stdc++.h> 
using namespace std;
// Baby-step Giant-step Algorithm
// a x + by = g
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 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 BSGS(long long P, long long B, long long N) {
    unordered_map<long long, int> R;
    long long sq = (long long) sqrt(P);
    long long t = 1, f;
    for (int i = 0; i < sq; i++) {		
        if (t == N)	
            return i; 		
        if (!R.count(t))
            R[t] = i;
        t = (t * B) % P;
    }
    
    f = inverse(t, P);
    for (int i = 0; i <= sq+1; i++) {
        if (R.count(N))
            return i * sq + R[N];
        N = (N * f) % P;
    }
    return -1;
}
int main() {
    long long P, B, N; // find B^L = N mod P
    while (scanf("%lld %lld %lld", &B, &N, &P) == 3) {
        long long L = BSGS(P, B, N);
        if (L == -1)
            puts("NOT FOUND");
        else
            printf("%lld\n", L);
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b433. 中間相遇法

Problem

背景

加密的每一道過程若存在反運算，則存在解密的程序。絕對安全的加密沒有反運算，那解密也是沒有辦法做到的。通常加解密中一定會用到互斥或 (XOR) 運算，從表格來看，單獨看任何一行一列都恰好分配一半的 0、一半的 1。對於密碼來說，明文跟密文的關係，不管知道的明文還是密文的部分，猜中的機率只有 $1/2$，只要位元數一多，猜中的機率近乎 $(\frac{1}{2})^n \cong 0$。

「只用 XOR key 是不行的，再做一次不就回來了嗎？」

「那用循環位移和加法的 overflow 破壞線性結構！」

「但記得要能解密回來哦。」

於是某 M 提供一個簡單的加密運算，明文 $M$，加密金鑰 $key$

  • 密文 $C = ((M \text{<<<} key) + key) \text{ XOR } key$
  • 明文 $M = ((C \text{ XOR } key) + (\sim key) + 1) \text{>>>} key$

其中每一個運算都在 16-bit CPU 上運作，$\text{<<<}$ 表示循環左移 (circular shift)，$\sim$ 表示 1 補數。

用抽象化表示加解密流程

  • $C = E_{key}(M)$
  • $M = D_{key}(C)$

問題描述

現在某 M 正用他自己設計的加解密協定跟未來的自己溝通，但發現到這種加解密，如果對方知道某 M 一定會傳送哪一個明文，接著只要匹配密文，窮舉 $O(2^{16})$ 來找到加密金鑰就會破功。

在電影《模仿遊戲》中，德國納粹 Enigma 密碼機，訊息中的結尾一定會附加「希特勒萬歲！(Heil Hitler!)」匹配密文後，窮舉金鑰就能破解。而某 M 常常傳送「萌萌哒！(Moe Moe Ta!)」只需要 $O(2^{16})$ 就能被破解。

於是某 M 強化他的加密協定，希望破解者至少要在 $O(2^{32})$ 的時間內找到，拖延破解時間。

$C = E_{key_2}(E_{key_1}(M))$ $M = D_{key_1}(D_{key_2}(C))$

現在知道某 M 傳送的其中一組 $(C, M)$，請告訴某 M 有多少組 $(key_1, key_2)$ 的可能性。萬萬沒想到，中間相遇法能在 $O(2^{16} \times 16)$ 解決這個問題。

Sample Input

1
2
40380 23333
30767 13657

Sample Output

1
2
114688
45568

Solution

套用中間相遇法$E_{key_1}(M) = D_{key_2}(C)$ 分別對等式兩邊進行計算，查看相碰情況。

複雜度為窮舉所有可能的 $key_1$ 和 $key_2$，若用平衡樹進行碰撞檢查 $O(2^{16} \times 16)$，利用 HASH 可以降到 $O(2^{16})$，由於 $O(2^{16})$ 是記憶體可以宣告的，則直接使用數組來完成。

Hash

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
#include <bits/stdc++.h> 
using namespace std;
typedef unsigned short int UINT16;
class SimpleIdea {
public:
    UINT16 key;
    SimpleIdea(UINT16 k = 0):
        key(k) {}
    UINT16 encrypt(UINT16 m) {
        return (rotate_left(m, key&15) + key)^key;
    }
    UINT16 decrypt(UINT16 m) {
        return rotate_left((m^key)+(~key)+1, 16 - (key&15));
    }
private:
    UINT16 rotate_left(UINT16 x, UINT16 n) {
    	return  (x << n) | (x >> (16-n));
    }
} test;
int main() {
    int C, M;
    while (scanf("%d %d", &C, &M) == 2) {
        unordered_map<int, int> H1, H2;
        for (int i = 0; i < (1<<16); i++) {
            int key1 = i, key2 = i;
            test.key = key1;
            H1[test.encrypt(M)]++;
            test.key = key2;
            H2[test.decrypt(C)]++;
        }
        
        int ret = 0;
        for (auto &x : H1) {
            if (H2.count(x.first))
                ret += H2[x.first] * x.second;
        }
        printf("%d\n", ret);
    }
    return 0;
}

Array

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
#include <bits/stdc++.h> 
using namespace std;
typedef unsigned short int UINT16;
class SimpleIdea {
public:
    UINT16 key;
    SimpleIdea(UINT16 k = 0):
        key(k) {}
    UINT16 encrypt(UINT16 m) {
        return (rotate_left(m, key&15) + key)^key;
    }
    UINT16 decrypt(UINT16 m) {
        return rotate_left((m^key)+(~key)+1, 16 - (key&15));
    }
private:
    UINT16 rotate_left(UINT16 x, UINT16 n) {
    	return  (x << n) | (x >> (16-n));
    }
} test;
int main() {
    int C, M;
    while (scanf("%d %d", &C, &M) == 2) {
        int H[1<<16] = {}, ret = 0;
        for (int i = 0; i < (1<<16); i++) {
            test.key = i;
            H[test.encrypt(M)]++;
        }
        for (int i = 0; i < (1<<16); i++) {
            test.key = i;
            ret += H[test.decrypt(C)];
        }
        printf("%d\n", ret);
    }
    return 0;
}
Read More +
解題區/解題區 - Zerojudge

b428. 凱薩加密

Problem

背景

曾幾何時，基礎題庫已經成了不基礎的題庫。小小新手們寫個題目，不少拿了 TLE、CE 求助無門，就再也不想打開 Zerojudge。高中生哪有寫這麼困難的題目，高中生都不像高中生。在某 M 那個年代寫的題目非常簡單，沒有特別變化處理，更別說多麼高檔的資料結構，暴力算法 (naive algorithm) 就能輕鬆切題。

「年代變了呢，現在的高中生要寫出比大學生的某 M 更困難的題目」

重溫解題的那份初心吧！

題目描述

在西元前就存在的一種加密－凱薩加密為目前最早發現的替換加密 (substitution cipher)。其原理很簡單，將一段明文往替換成往後數的第 $k$ 個英文字母。

若用數學式表示凱薩加密和解密，如下：

加密 $C = E_K(P) = (P + k) \mod 26$
解密 $P = D_K(P) = (C - k) \mod 26$

例如 $k = 3$ 時，發生的情況如下：

明文字母表：ABCDEFGHIJKLMNOPQRSTUVWXYZ
密文字母表：DEFGHIJKLMNOPQRSTUVWXYZABC

從數學的觀點來看，每一個字母就是一個數字。

A = 0, B = 1, C = 2, …,X = 23, Y = 24, Z = 25

Sample Input

1
2
3
4
5
6
ABCDEFGHIJKLMNOPQRSTUVWXYZ 
DEFGHIJKLMNOPQRSTUVWXYZABC
DLQXABXEEQMEUQYLZPEK
YGLSVWSZZLHZPLTGUKZF
Z
Z

Sample Output

1
2
3
3
21
0

Solution

只需要看第一個字符的變換方式即可。

1
2
3
4
5
6
7
8
9
#include <bits/stdc++.h>
using namespace std;
int main() {
    char s1[1024], s2[1024];
    while (scanf("%s %s", s1, s2) == 2)
        printf("%d\n", (((s2[0] - 'A') - (s1[0] - 'A'))%26 + 26)%26);
    return 0;
}
Read More +