poj 3461 Oulipo 后缀数组

缘起

【1】中用后缀树切的时候，MLE掉了. 因为MLE消耗内存甚巨. 所以这次采用后缀数组的姿势来切.

分析

题目见【1】好了.

先来说说使用后缀数组的切法. 首先我们使用da或者dc3得到了原串的sa之后我们能干什么? 首先sa是什么?

sa[1,…,n]得到的就是n个后缀（因为sa[0]才是空串, 所以sa[1,…,n]并没有包含空串）的字典序升序排序. 等等! 升序!! 是不是有点二分的感觉出来了? 没错, 我们使用二分先查找前缀包含模式串的后缀的最大索引, 再使用二分查找前缀包含模式串的后缀的最小索引, 则这两个索引的差值就是要求的个数. 复杂度是 O(n)+O(模式串长度*logn)

注意，查找部分的复杂度降低到了logn, 比kmp的On复杂度要低. 但是这里构建后缀树就要花费On时间. 所以总体复杂度依旧是O(n)的. 但是在精确匹配领域, 后缀树、后缀数组都和kmp、bm、karp-rabin算法无法匹敌的。因=究其根本是因为前者需要对文本串建立点儿什么, 而后者仅仅需要对模式串建立点儿什么. 而实际问题中，文本串的长度远大于模式串的长度. 所以前者的性能远不及后者. 后者复杂度是O(n+m), 而且常数几乎就是1.

举个例子

aabaaaab, 它的后缀数组是

1
2
3
4
5
6
7
8
9
aaaab   	即sa[1]=3
aaab		即sa[2]=4
aab			即sa[3]=5
aabaaaab	即sa[4]=0
ab			即sa[5]=6
abaaaab		即sa[6]=1
b			即sa[7]=7
baaaab		即sa[8]=2
之所以是sa[1,...,8]是因为【2】中的后缀数组算法最后得到的sa的索引范围是1~n,但是值在[0,n-1]范围内

假设模式串是aaa, 则sa[1,…,n]中包含aaa的后缀最小为sa[1], 最大为sa[2], 即以模式串为前缀的后缀的在sa中的最大索引是2, 在sa中的最小索引是1, 所以模式串出现的次数是2-1+1=2

于是不难写下如下代码（使用da）

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
#include "stdafx.h"
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
#define LOCAL

const int maxs=1000005, maxp = 10005;
char s[maxs],p[maxp];
int k,n,sa[maxs],r[maxs],t[maxs];

bool c0(int i, int j)
{
	if (r[i]!=r[j])
	{
		return r[i]<r[j];
	}
	int ri = i+k<=n?r[i+k]:-1, rj = j+k<=n?r[j+k]:-1;
	return ri<rj;
}

void da()
{
	k = 1, n = 0;
	while (s[n])
	{
		sa[n] =n;
		r[n] = s[n]-'A'+1;
		n++;
	}
	r[n]=0, sa[n]=n;
	while(k<=n)
	{
		sort(sa, sa+n+1, c0);
		t[sa[0]] = 0;
		for (int i =1;i<=n;i++)
		{
			t[sa[i]] = t[sa[i-1]]+c0(sa[i-1], sa[i]);
		}
		memcpy(r,t,sizeof(int)*(n+1));
		k<<=1;
	}
}

int ck(int x) // S[sa[x],...] 这个后缀是否包含模式串p? 返回0表示包含, 返回1表示不包含并且严格比模式串大, 返回-1表示不包含并且严格比模式串小
{
	int i = sa[x], j = 0;
	for (; p[j];i++,j++)
	{
		if (s[i]!=p[j])
		{
			return s[i]>p[j]?1:-1;
		}
	}
	return 0;
}

int kk(bool flag) // 二分求出最早和最晚包含模式串的后缀在sa中的索引, flag=true表示求出包含模式串t在sa中的最大索引,false表示求出包含模式串t在sa中的的最小索引
{
	int lo = 1, hi = n, ans = 0, mid;
	while(lo<=hi)
	{
		mid = lo+hi>>1;
		int check = ck(mid);
		if (!check)
		{
			ans = mid;
			flag?lo = mid+1:hi = mid-1;
		}
		else
		{
			if (check==1) // 不包含并且严格比模式串大
			{
				hi = mid-1;
			}
			else // 不包含并且严格比模式串小
			{
				lo = mid+1;
			}
		}
	}
	return ans;
}

int main() {
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
//	freopen("d:\\my.out", "w", stdout);
#endif
	int kase;
	scanf("%d", &kase);
	while(kase--)
	{
		scanf("%s%s",p,s);
		da();
		int lo = kk(false), hi = kk(true);
		printf("%d\n", hi?hi-lo+1:0);
	}
	return 0;
}

ac情况

遗憾的被T掉了. 其实不难知道 da的复杂度是 O(nlogn)的，尽管常数较小，但是logn（以2为底的对数）达到了惊人的20，即da建立后缀数组的复杂度至少是两千万这个级别的（我还没乘以常数）.

但是和ac代码对拍, 不会wa，仅仅是被t掉了.

基本可以断言——da只能打到1e5级别的数据. 1e6就有些费劲了.

试试dc3, 注意, 因为在【2】中dc3和da的输出结构都是完全一样的. 所以可以简单的将上面被T掉的代码换成dc3的.

下面的代码依旧被T掉的了, 但是对拍的话, 还是不会wa的.

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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
//#include "stdafx.h"
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
//#define LOCAL
#define F(x) ((x)/3+((x)%3==1?0:tb))
#define G(x) ((x)<tb?(x)*3+1:((x)-tb)*3+2)
const int maxs=1000005, maxp = 10005;
char s[maxs],pp[maxp];
int n, sa[maxs*3], r[maxs*3], wa[maxs*3], wb[maxs*3], wv[maxs*3], ws[maxs*3];

void sort(int *r, int *a, int *b, int n, int m)
{
	for (int i = 0;i<n; i++)
	{
		wv[i] = r[a[i]];
	}
	for (int i = 0;i<m;i++)
	{
		ws[i] = 0;
	} 
	for (int i =0; i<n; i++)
	{
		ws[wv[i]]++;
	}
	for (int i = 1;i<m;i++)
	{
		ws[i]+=ws[i-1];
	}
	for (int i = n-1; ~i; i--)
	{
		b[--ws[wv[i]]]=a[i];
	} 
}

int c0(int *r, int a, int b)
{
	return r[a]==r[b]&&r[a+1]==r[b+1]&&r[a+2]==r[b+2];
}

bool c12(int k, int *r, int a, int b)
{
	if (k==2) 
	{
		return r[a]<r[b]||r[a]==r[b]&&c12(1,r,a+1,b+1);
	}
	else 
	{
		return r[a]<r[b]||r[a]==r[b] && wv[a+1]<wv[b+1];
	}
}

void dc3(int *r, int *sa, int n, int m)
{
	int *rn = r+n, *san = sa+n, ta = 0, tb = (n+1)/3, tbc=0, p; 
	r[n] = r[n+1] =0;
	for (int i=0; i<n; i++) 
	{
		if (i%3!=0)
		{
			wa[tbc++] = i;
		}
	} 
	sort(r+2, wa, wb, tbc, m);
	sort(r+1, wb,wa, tbc, m);
	sort(r, wa,wb,tbc, m);
	rn[F(wb[0])]=0; p=1; 
	for (int i = 1; i<tbc; i++)
	{
		rn[F(wb[i])] = c0(r, wb[i-1], wb[i])?p-1:p++;
	} 
	if (p<tbc)
	{
		dc3(rn, san, tbc, p); 
	} 
	else 
	{
		for (int i =0 ; i<tbc; i++)
		{
			san[rn[i]] = i;
		}
	}
	for (int i = 0; i<tbc; i++)
	{
		if (san[i]<tb)
		{
			wb[ta++] = san[i]*3;
		}
	} 
	if (n%3==1) 
	{
		wb[ta++] = n-1;
	}
	sort(r,wb,wa,ta,m);
	for (int i =0;i<tbc; i++)
	{
		wv[wb[i]=G(san[i])] = i; 
	} 
	int i=0,j=0;
	for (p=0; i<ta && j<tbc; p++) 
	{
		sa[p] = c12(wb[j]%3, r, wa[i], wb[j])?wa[i++]:wb[j++]; 
	}
	while (i<ta)
	{
		sa[p++]=wa[i++];
	}
	while(j<tbc)
	{
		sa[p++]=wb[j++];
	}
}


int ck(int x) 
{
	int i = sa[x], j = 0;
	for (; pp[j];i++,j++)
	{
		if (s[i]!=pp[j])
		{
			return s[i]>pp[j]?1:-1;
		}
	}
	return 0;
}

int kk(bool flag)
{
	int lo = 1, hi = n, ans = 0, mid;
	while(lo<=hi)
	{
		mid = lo+hi>>1;
		int check = ck(mid);
		if (!check)
		{
			ans = mid;
			flag?lo = mid+1:hi = mid-1;
		}
		else
		{
			if (check==1) 
			{
				hi = mid-1;
			}
			else
			{
				lo = mid+1;
			}
		}
	}
	return ans;
}

int main() {
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
	//	freopen("d:\\my.out", "w", stdout);
#endif
	int kase;
	scanf("%d", &kase);
	while(kase--)
	{
		scanf("%s%s",pp,s);
		n = 0;
		while(s[n])
		{
			r[n] = s[n]-'A'+1;
			n++;
		}
		r[n]=0;
		dc3(r,sa,n+1,27);
		int lo = kk(false), hi = kk(true);
		printf("%d\n", hi?hi-lo+1:0);
	}
	return 0;
}

上面的代码没有注释, 详细注释参见 【2】.

但是本地跑的时候明显能感觉得到dc3比da在1e6级别的表现上出色很多.

参考

【1】https://yfsyfs.github.io/2019/10/06/poj-3461-Oulipo-%E5%90%8E%E7%BC%80%E6%A0%91-KMP/

【2】https://yfsyfs.github.io/2019/10/07/%E5%90%8E%E7%BC%80%E6%95%B0%E7%BB%84%E5%AD%A6%E4%B9%A0/