后缀数组

2. #include <iostream>
3. #include <string>
4. #include <algorithm>
5. using namespace std;
7. #define Min(a,b) (a)<(b)?(a):(b)
8. const int N = 2000;
9. const int D = 10;
10. const int M = 260;
11. int n;
12. char s[N] = "abcabdeabf";
13. int cnt[N], mem[4][N], *rank, *nrank, *sa, *nsa, lcp[20][N];
14. // lcp[i][j]: longest commen prefix ( suffix(sa[k+1]), suffix(sa[k]) ) j <= k < j+2^i
16. void suffix()
17. {
18.  int i, j, k;
19.  rank = mem[0];
20.  nrank = mem[1];
21.  sa = mem[2];
22.  nsa = mem[3];
24.  for(i = 0; i < n; i++) cnt[s[i]]++;
25.  for(i = 1; i < M; i++) cnt[i] += cnt[i-1];
26.  for(i = n-1; i >= 0; i--) sa[--cnt[s[i]]] = i;
27.  for(rank[0]=0, i=1; i < n; i++)
28.  {
29.   rank[sa[i]] = rank[sa[i-1]];
30.   if(s[sa[i]]!=s[sa[i-1]]) rank[sa[i]]++;
31.  }
33.  for(k = 1; k<n && rank[sa[n-1]] < n-1; k*=2)
34.  {
35.   for(i = 0; i < n; i++) cnt[rank[sa[i]]] = i+1;
36.   for(i = n-1; i >= 0; i--) if(sa[i]-k>=0)
37.    nsa[--cnt[rank[sa[i]-k]]] = sa[i]-k;
39.   // max(sa[i]-k)=n-k-1 , therefore i = n-k;
40.   for(i = n-k; i < n; i++)
41.    nsa[--cnt[rank[i]]] = i;
43.   for(nrank[nsa[0]], i=1; i < n; i++)
44.   {
45.    nrank[nsa[i]] = nrank[nsa[i-1]];
46.    if(rank[nsa[i]] != rank[nsa[i-1]]
47.     || rank[nsa[i]+k] != rank[nsa[i-1]+k])
48.    nrank[nsa[i]]++;
49.   }
50.   swap(rank, nrank);
51.   swap(sa, nsa);
52.  }
53. }
55. void print()
56. {
57.  for(int i=0;i<n;i++) printf("%d %s/n", sa[i], s + sa[i]);puts("");
58.  puts(s);for(int i=0;i<n;i++) printf("%d ", rank[i]);puts("");
59. }
61. void get_lcp_rmq()
62. {
63.  int i, j, k;
64.  for(i=0,k=0; i<n; i++)
65.  {
66.   if(rank[i]==n-1) lcp[0][rank[i]]=k=0;
67.   else
68.   {
69.    if(k>0)k--;
70.    j = sa[rank[i]+1];
71.    for(;s[i+k]==s[j+k];k++) ;
72.    lcp[0][rank[i]]=k;
73.   }
74.  }
76.  for(i=0,k=1; k<n; k*=2, i++)
77.  {
78.   for(j = 0; j+k < n; j++)
79.    lcp[i+1][j] = Min(lcp[i][j], lcp[i][j+k]);
80.  }
81. }
83. int rmq(int a, int b)
84. {
85.  int i, j, k;
87.  if(a==b) return n-1-a;
88.  a = rank[a];
89.  b = rank[b];
90.  if(a>b) swap(a,b);
92.  int t=b-a;
93.  for(i=0,k=1; 2*k<t; i++,k*=2) ;
94.  return Min(lcp[i][a], lcp[i][b-k]); // not  b-k+1
95.  // lcp[0][i]: LCP(i,i+1)
96. }
100. int main()
101. {
102.  int i, j, k;
104.  n = strlen(s);
105.  s[n++]=0;
106.  memset(cnt, 0, sizeof(cnt));
107.  memset(mem, 0, sizeof(mem));
109.  suffix();
110.  print();
112.  get_lcp_rmq();
113.  printf("%d/n", rmq(2,3));
115.     return 0;
116. }