Code

LCA和RSA

Tarjan离线求解LCA和RSA加密算法。

Tarjan离线求解LCA

LCA（Least Common Ancestors），最近公共祖先，即在有根树中寻找到两个结点u，v距离和最短的结点，亦即离根最远的结点。

  • 算法基于DFS和并查集，时间复杂度由DFS的V和查询次数M决定，O（V+M）。在下面给出的代码样例中采用链式结构存储树，采用邻接表结构双向存储查询端点，也可以通过两个边集数组实现。

  • 算法流程：

    • 初始化：设置par数组（并查集），par[i] = i。离线所有查找请求，存储于无向图邻接表中。
    • 从根结点开始tarjan（DFS），对于每次新探求到的顶点u，对照每个查询请求（u，v），判断这条询问请求中的另一个结点v是否被搜索过。如果v尚未被搜索则暂时跳过该条请求（在我的样例中将这条请求留到了（v，u）处理），如果v已经被搜索过，则LCA（u，v）= Find（v），这里的Find是并查集操作。
    • 搜索并处理u包含的所有子树（每棵子树搜索完成意味着所有子树内的LCA询问都已经被解决）。在v的所有子树被搜索完成后，回溯时合并祖先，将所有自v处搜索出的顶点的par值设为u。
    • tarjan所有结点。
    • 核心部分的伪代码
1
2
3
4
5
6
7
8
9
tarjan(u)
     vis[u] <- true
     for each (u, v) in questions
         if vis[v]
             ans(u, v) <- Find(v)
     for each (u, v) in tree
         if !vis[v]
             tarjan(v)
             par[v] <- u

  • 对LCA（u，v）= Find（v）给出证明：因为v在遍历到u（当前结点）前被遍历过，而又由某个结点拓展到u，所以以最近公共祖先为根并且包含v的子树已经被搜索过，若有一个从当前结点u到结点v的查询，这时当前这棵子树的祖先就是最近公共祖先（注意当前的含义，当对v的搜索不断回溯时，Find（v）的值会在不断变化）。

  • 样例

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
#include <stdio.h>
#include <memory.h>
#include <stdlib.h>

typedef struct Node {
	int num;
	struct Node * left, * right;
}node;

typedef struct granode{
	int v;
	int ans;
	struct granode *next;
}granode; 

int V, E, Q;
int * vis, * par;
node * tree, * root = NULL;
granode ** list;

int Find( int u ){
	return par[u] == u ? u : Find( par[u] );
}

void insert( int u, int v ){
	granode * temp = ( granode * ) malloc ( sizeof( granode ) );
	temp->v = v;
	temp->next = list[u];
	list[u] = temp;
}

void init(){
	scanf( "%d %d %d", &V, &E, &Q );
	int i, father, child, u, v;
	tree = ( node * ) malloc ( sizeof( node ) * ( V + 1 ) );
	par = ( int * ) malloc ( sizeof( int ) * ( V + 1 ) );
	vis = ( int * ) malloc ( sizeof( int ) * ( V + 1 ) );
	list = ( granode ** ) malloc ( sizeof( granode* ) * ( V + 1 ) );
	memset( vis, 0, sizeof( int ) * ( V + 1 ) );
	for ( i = 1 ; i <= V ; i++ ){
		tree[i].num = i;
		tree[i].left = tree[i].right = NULL;
		list[i] = NULL;
		par[i] = i;
	}
	for ( i = 0 ; i < E ; i++ ){
		scanf( "%d %d", &father, &child );
		vis[child] = 1;
		if ( tree[father].left == NULL )
			tree[father].left = tree + child;
		else 
			tree[father].right = tree + child;
	}
	for ( i = 1 ; i <= V ; i++ )
		if ( !vis[i] )
			root = tree + i;
	memset( vis, 0, sizeof( int ) * ( V + 1 ) );
	for ( i = 1 ; i <= Q ; i++ ){
		scanf( "%d %d", &u, &v );
		insert( u, v );
		insert( v, u );
	}
}

void tarjan( node * n ){
	vis[ n->num ] = 1;
	int u = n->num;
	granode * temp = list[u];
	while ( temp != NULL ){
		if ( vis[ temp->v ] )
			temp->ans = Find( temp->v );
		temp = temp->next;
	}
	if ( n->left != NULL )
		if ( !vis[ n->left->num ] ){
			tarjan( n->left );
			par[ n->left->num ] = u;
		}
	if ( n->right != NULL )
		if ( !vis[ n->right->num ] ){
			tarjan( n->right );
			par[ n->right->num ] = u;
		}
}

void print(){
	int i;
	for ( i = 1 ; i <= V ; i++ ){
		granode * temp = list[i];
		while ( temp != NULL ) {
			if ( temp->ans )
				printf( "%d %d %d\n", i, temp->v, temp->ans );
			temp = temp->next;
		}
	}
}

int main(){
	init();
	printf("Root : %d\n",root->num );
	tarjan( root );
	print();
	return 0;
}
  • 运行结果
    LCA的Tarjan离线算法演示

RSA

选取两个不同的大素数p、q，并计算N = p×q，选取小素数d，并求出e = ~d，使 d × e ≡ 1 mod (p-1)(q-1)。对于任意明文 A < N有密文B =A^d mod N，解密时A = B^e mod N。因此d和~d形成了非对称密钥关系，使用公钥d加密，私钥e解密，即使第三方获取了密文B、公钥d和模数N，因为目前的RSA大多以大于1024位的大数为基础，因此对N进行素数因子分解近似不可能（不考虑量子算法），因此无法获取p，q，也就无法推算e。

大数运算

目前多为32和64bit编译器，对于1024位以上的大数运算通常将其作为数组处理。如果选取十进制，2^1024会产生一个长度300多位的数组，对于四则运算的任何一种操作都需要在两个数百位的数组上进行，效率非常低，并且将其转换成二进制非常麻烦。通常的改进方法是将大数用n进制数组表示，这里的n为2的幂，对32bit系统可以取n = 2 ^ 32（0x100000000），将1024位大数用该进制表示即为32位，每一位的取值是2^32进制的0～0xffffffff。这刚好是32位系统中long的长度。此时对该数组进行的任意运算所需的循环规模不超过32次，并且因为进制是2的次幂，与二进制的转换非常容易。

  • 存储方式：例如0xffffffff ffffffff，转换成十进制是18446744073709551615，格式类似于十进制中的99，十六进制中的0xFF。而0x00000001 00000000 00000000，转换成十进制是18446744073709551615 + 1 = 18446744073709551616，格式类似十进制的100，十六进制的0x0100。在数组中，i从0到n分别按由低位到高位的顺序存储。对于一个数组A表示的大数，其值通过如下方式计算：for i from 0 to n { ans <- ans + A[i] × base ^ i }。显然各个数字间的运算中间结果可能超出32位长度，因此需要64位的数据（long long，__int 64）类型存储（部分32位编译系统不支持64位，可使用2^16进制，并将32位数据类型用来存储中间结果；或者将2^32进制改成4×2^28进制，也就是0x40000000进制，这样两个数字运算最大结果为0x3fffffff × 0x3fffffff，不超过0xffffffffffffff，可以用double类型存储，但不能将中间结果简单的拆分成高位和低位）。

假设存在大数A和B，A.length和B.length表示各自长度，base = 0x100000000，0 <= A[i],B[i] < base。

  • 加法：处理进位情况和长度变化即可 
1
2
3
4
5
6
7
8
__int64  temp;
int pre = 0, len = max( A.length, B.length );
for ( i = 0 ; i <= len ; i++ ){
	temp = A[i] + B[i] + pre;
	res[i] = temp % 0x100000000;
	pre = temp / 0x100000000;
}
res.length = pre == 0 ? len : len + 1;
  • 减法：对中间结果是负数的需要借位 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int pre = 0;
int len = max( A.length, B.length );
for ( i = 0 ; i <= len ; i++ ){
	res[i] = A[i] - B[i] - pre;
	if ( res[i] >= 0 )
		pre = 0;
	else {
		res[i] += 0x100000000;	// 借位
		pre = 1;
	}
}
res.length = len;
while ( res[res.length] == 0 )
	res.length--;		//去除大数前端多余的0
  • 乘法：竖式运算，假设B.length < A.length，则res[i] = Sum{ A[i-j] × B[j], j from 0 to B.length }。res.length初始可以设置为A.length + B.length - 1，考虑到进位，最后可以根据情况增加一位。 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
__int temp;
int pre = 0;
res.length = A.length + B.length - 1;
for ( i = 0 ; i <= res.length ; i++ ){
	temp = pre;
	for ( j = 0 ; j <= B.length ; j++ ){
		if ( i >= j && i <= j + A.length ){
			temp += A[i-j] * B[j];		// 竖乘
			res[i] = temp % 0x100000000;
			pre = temp / 0x100000000;
		}
	}
}
if ( pre )
	res[++res.length] = pre;	// 最高位进位
  • 除法：通过蒙哥马利算法可以通过移位操作代替除法求得模乘运算的结果

米勒·拉宾算法

  • 算法流程
    • 根据RSA中的N计算出一个奇数M，使得N = 1 + M × 2^r 。
    • 选择一个小于N的随机数C，对于任意的i < r，如果C^( ( M × 2^i ) % N) == N - 1或者C^M % N == 1，则N通过随机数C的测试。
    • 重复上一步五次，如果全部通过则判定N为素数

若N通过一次测试，则N为合数的概率为25%，当N通过t次测试，N为合数的概率为0.25^t。t为5时N为素数的概率已经大于99.9%。实际操作中通常用300个小素数对N进行测试，如果N通过小素数测试再使用米勒·拉宾测试，提高测试速度。

欧几里德方程

根据辗转相除法，当gcd(a,b) = 1，有sa + tb = 1。将求解a和b的gcd过程写出，逆向即可求得~b，使得b×~b ≡ 1（mod a）。

  • 辗转相除法的逆运算
    求解gcd（101,4620），有如下：
    4620 = 45 × 101 + 75
    101 = 1 × 75 + 26
      75 = 2 × 26 + 23
      26 = 1 × 23 + 3
       3 = 1 × 2 + 1
       2 = 2 × 1
    对上述操作求逆操作，为：
    1 = 3 - 1 × 2
      = 3 - 1 ×（23 - 7 × 3）= -1 × 23 + 8 × 3
      = -1 × 23 + 8 ×（26 - 1 × 23）= 8 × 26 - 9 × 23
      = 8 × 26 - 9 ×（75 - 2 × 26）= -9 × 75 + 26 × 26
      = -9 × 75 + 26 ×（101 - 1 × 75）= 26 × 101 - 35 × 75
      = 26 × 101 - 35 ×（4620 - 45 × 101）= -35 × 4620 + 1601 × 101
  • 下列代码传入的参数r最后为求得的~b。
1
2
3
4
5
6
7
8
9
10
void get( int a, int b, int *l, int *r){	/* 逆向欧几里得求 1 = st + qa */ 
	if ( a % b == 1 ){
		*l = 1;
		*r = -a / b;
	}
	else {
		get( b, a % b, l, r );
		(*r) = (*r) * ( - a / b ) + *l; 				/*  r为求得的modular  */ 
	}	
}
  • 非递归的代码如下( sa + tb = 1 )，返回的t即为~b。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int get( int a, int b ){
	int t = 0, s = 1, temp;
	while ( b ! = 0 ){
		temp = s;
		s = t - s * b / a;
		t = temp;
		temp = b;
		b = a % b;
		a = temp;
		if ( t < 0 )
			t += a;
	}
	return t;
}

Montgomery幂模算法

Montgomery（蒙哥马利）算法的优点在于减少了取模的次数，并且用移位取代除法。举例：进制base = 13，对于32 = 2 × 13 + 7。要求32×13^k mod 7的值，可以对32进行如下操作：32 = 32 + 7 × q，这样不影响取模的结果。当32 + 7 × q是13的整数倍的时候就可以用移位操作除以13，直到求出结果。

  • 伪代码
1
2
3
4
5
6
7
8
9
func() a^b mod c
	res <- 1
	pow <- b
	while a > 0
		pow <- ( pow * pow ) mod c
		if a & 1 = 1		// 以二进制为例，对于n进制有同样的增值运算
			res <- ( res * pow ) mod c
		a >> 1			// 对于任意进制同样移位
	return res

原创作品，允许转载，转载时无需告知，但请务必以超链接形式标明文章原始出处(https://forec.github.io/2015/09/10/Algorithms-Must1/) 、作者信息（Forec）和本声明。

分享到