首页 > itarticle > b444. 期望試驗 快速冪次

b444. 期望試驗 快速冪次

admin 11月 11, 2020 0

contents

Problem

背景

曾經某 M 被期望值坑，就只是在計算 $x^y mod z$ 時偷偷替換成 $x^{y-1} times x mod z$，結果得到 Time Limit Exceeded。

根據分析 $y = 16$ 時，用二進制表示為 $(10000)_{2}$，若變成 $y = 15$，就會變成 $(01111)_{2}$，通常快速求冪的乘法次數與二進制的 1 個數成正比，所以速度就慢非常多。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
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; 
}
UINT64 mpow(UINT64 x, UINT64 y, UINT64 mod) {
	UINT64 ret = 1;
	while (y) {
		if (y&1)
			ret = mul(ret, x, mod);
		y >>= 1, x = mul(x, x, mod);
	}
	return ret;
}

問題描述

讓我們來一場 $x^y mod z$ 的期望值試驗吧，基礎目標是減少乘法次數。

其中 $1 le x, z le 10^{18}$, $0 le y le 2^{2^{20}}$，在這場試驗中展現你的優化吧。

Sample Input

1
2
3
4
 
5 01 7
5 10 7
5 0101 7
3 0110 7

Sample Output

1
2
3
4
 
5
4
3
1

Solution

詳細可以參考《資訊安全 - 近代加密 快速冪次計算》那一篇。

Algorithm Table Size #squaring Average #Multiplication
Right-To-Left 1: $x^{2^i}$ $n$ $n/2$
Left-To-Right 1: $x$ $n$ $n/2$
Left-To-Right(2-bits) 3: $x$, $x^2$, $x^3$ $n$ $3n/8$
Left-To-Right(sliding) 2: $x$, $x^3$ $n$ $n/3$

減少乘法次數，但以上期望乘法次數是跟 1 的個數有關，雖然最好是從 $n/2$ 降到 $n/3$，並不表示速度會真的快上 $1.5$ 倍左右，畢竟還有所謂的基礎乘法次數需求，根據實驗下來大約能快個 10% 到 20% 之間，加上 -Ofast 編譯此時的差異又會再少一點，看起來實作方法影響很嚴重。

例如在不加編譯優化參數下

1
2
3
4
 
if a[i] == 0 && a[i+1] == 0
else if a[i] == 0 && a[i+1] == 1
else if a[i] == 1 && a[i+1] == 0
else

上述做法會比下述來得快上許多

1
2
3
4
5
6
 
if a[i] == 0
	if a[i+1] == 0
	else
else
	if a[i+1] == 0
	else

最後，產出一個 cheat 版本，使用 L-to-R-2bits 的概念下去擴充，使用 loop unrolling 進行加速，由於會發生不被整除的問題，小測資就靠 L-to-R-sliding 的方案去解決。

在 zerojudge 主機上平台上，隨機測資下的運作情況如下：

Algorithm Time
Right-To-Left 5.4s
Left-To-Right(2-bits) 4.9s
Left-To-Right(sliding) 4.8s
Left-To-Right-sliding-cheat 4.2s

R-to-L

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
 
#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; 
}
UINT64 mpowR2L(UINT64 x, char y[], UINT64 z) {
	int n;
	for (n = 1; y[n]; n <<= 1);
	UINT64 ret = 1;
	for (int i = n-1; i >= 0; i--) {
		if (y[i] == '1')
			ret = mul(ret, x, z);
		x = mul(x, x, z);
	}
	return ret;
}
char y[(1<<20) + 5];
int main() {
	long long x, z;
	while (scanf("%lld %s %lld", &x, y, &z) == 3) {
		printf("%llun", mpowR2L(x, y, z));
	}
	return 0;
}

L-to-R-2bits

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;
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; 
}
UINT64 mpowL2R2(UINT64 x, char y[], UINT64 z) {
	UINT64 x2 = mul(x, x, z);
	UINT64 x3 = mul(x2, x, z);
	UINT64 ret = 1;
	for (int i = 0; y[i]; i += 2) {
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		if (y[i] == '1' && y[i+1] == '1') {
			ret = mul(ret, x3, z);
		} else if (y[i] == '1' && y[i+1] == '0') {
			ret = mul(ret, x2, z); 
		} else if (y[i+1] == '1') {
			ret = mul(ret, x, z);
		}
	}
	return ret;
}
char y[(1<<20) + 5];
int main() {
	long long x, z;
	while (scanf("%lld %s %lld", &x, y, &z) == 3) {
		printf("%llun", mpowL2R2(x, y, z));
	}
	return 0;
}

L-to-R-sliding

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
 
#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; 
}
UINT64 mpowL2RS(UINT64 x, char y[], UINT64 z) {
	UINT64 x3 = mul(mul(x, x, z), x, z);
	UINT64 ret = 1;
	for (int i = 0; y[i]; ) {
		ret = mul(ret, ret, z);
		if (y[i] == '1' && y[i+1] == '1') {
			ret = mul(ret, ret, z);
			ret = mul(ret, x3, z);
			i += 2;
		} else if (y[i] == '1') {
			ret = mul(ret, x, z);
			i ++;
		} else if (y[i+1] == '0') {
			ret = mul(ret, ret, z);
			i += 2; 
		} else {
			i++;
		}
	}
	return ret;
}
char y[(1<<20) + 5];
int main() {
	long long x, z;
	while (scanf("%lld %s %lld", &x, y, &z) == 3) {
		printf("%llun", mpowL2RS(x, y, z));
	}
	return 0;
}

L-to-R-sliding-cheat

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
 
#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; 
}
UINT64 mpowL2RS(UINT64 x, char y[], UINT64 z) {
	UINT64 x3 = mul(mul(x, x, z), x, z);
	UINT64 ret = 1;
	for (int i = 0; y[i]; ) {
		ret = mul(ret, ret, z);
		if (y[i] == '1' && y[i+1] == '1') {
			ret = mul(ret, ret, z);
			ret = mul(ret, x3, z);
			i += 2;
		} else if (y[i] == '1') {
			ret = mul(ret, x, z);
			i ++;
		} else if (y[i+1] == '0') {
			ret = mul(ret, ret, z);
			i += 2; 
		} else {
			i++;
		}
	}
	return ret;
}
#define PREPROC 8
UINT64 mpowCHEAT(UINT64 x, char y[], UINT64 z) {
	int n;
	for (n = 1; y[n]; n <<= 1);
	if (n < 1<<PREPROC)
		return mpowL2RS(x, y, z);
	UINT64 X[1<<PREPROC] = {1};
	for (int i = 1; i < (1<<PREPROC); i++)
		X[i] = mul(X[i-1], x, z);
	UINT64 ret = 1;
	for (int i = 0, v; y[i]; i += PREPROC) {
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		ret = mul(ret, ret, z);
		v = (y[i]-'0')<<7|(y[i+1]-'0')<<6|(y[i+2]-'0')<<5|(y[i+3]-'0')<<4|(y[i+4]-'0')<<3|(y[i+5]-'0')<<2|(y[i+6]-'0')<<1|(y[i+7]-'0'); 
		ret = mul(ret, X[v], z);
	}
	return ret;
}
char y[(1<<20) + 5];
int main() {
	long long x, z;
	while (scanf("%lld %s %lld", &x, y, &z) == 3) {
		printf("%llun", mpowCHEAT(x, y, z));
	}
	return 0;
}
微海报分享

近期文章

近期评论

标签

热门

文章归档

分类目录

功能