字尾陣列——羅穗騫倍增演算法程式碼詳解
首先解釋一下用到的幾個陣列。
陣列sa:構造完成前表示關鍵字陣列,下標表示名次,值表示關鍵字的首字元位置,值相同的時候名次根據在原串中相對位置的先後決定;構造完成後表示字尾陣列,下標表示名次,值表示字尾的首字元位置。
陣列x:表示rank陣列,下標表示關鍵字的位置,值表示關鍵字大小(rank),相同的值有相同的rank。初始化為字串r的每個字元大小(此時x並不代表rank,只借助其值比較相對大小)。在每次迭代後,根據sa重新賦值,並代表rank。
陣列y:排序後的第二關鍵字陣列,下標表示名次,值代表第二關鍵字的首字元位置。
for(i = 0; i < m; i++) ws[i] = 0; for(i = 0; i < n; i++) ws[x[i] = r[i]]++; for(i = 1; i < m; i++) ws[i] += ws[i-1]; for(i = n-1; i >= 0; i--) sa[--ws[x[i]]] = i;
這段程式碼是計數排序,是對每個字元排序用來初始化陣列sa和x。這段程式碼不難懂,原理我就不講了,重點在於這個計數排序有個特點:
先根據x的值排序,x值相等時根據出現先後次序排序。
for(j = 1,p = 1; p < n ; j <<= 1,m = p)
進入迴圈,開始構造字尾陣列,j表示關鍵字長度。
這裡是對第二關鍵字排序。由於第一關鍵字和第二關鍵字都是長度為j的關鍵字,所以這裡完全可以使用之前已經算好的第一關鍵字,但是第二關鍵字和第一關鍵字有些不同。看圖。for(p = 0, i = n - j; i < n; i++) y[p++]=i; for(i = 0; i < n; i++) if(sa[i] >= j) y[p++] = sa[i] - j;
上面的藍線指向的是一個關鍵字。以長度為4時為例,下標超過4(圖中rank為2)時的關鍵字將不存在第二關鍵字。所以從n-j到n的關鍵字的第二關鍵字都應該為0,排序的時候應該放在前面。下面這段程式碼就表達了這樣的意思。
for(p = 0, i = n - j; i < n; i++)
y[p++]=i;
其餘部分可以利用第一關鍵字的排序進行。但是第二關鍵字的下標和第一關鍵字的下標是不一樣的,第一關鍵字的值對應下標就是該關鍵字的下標,第二關鍵字的值的下標減去j才是對應關鍵字的下標。舉個例子,上面的第一個組合關鍵字是(4,2)下標為0。其中4是第一關鍵字下標0;2是第二關鍵字下標是4,需要減去長度4才能對應到組合關鍵字的下標。
for(i = 0; i < n; i++)
if(sa[i] >= j)
y[p++] = sa[i] - j;
for(i = 0; i < n; i++)
wv[i] = x[y[i]];
for(i = 0; i < m; i++)
ws[i] = 0;
for(i = 0; i < n; i++)
ws[wv[i]]++;
for(i = 1; i < m; i++)
ws[i] += ws[i-1];
for(i = n-1; i >= 0; i--)
sa[--ws[wv[i]]] = y[i];
這段程式碼還是計數排序,其中這句wv[i] = x[y[i]];是最讓人費解的。我認為這也是這段程式碼最巧妙所在。
首先思考一下接下來我們應該幹什麼?很明顯是根據合併第一第二關鍵字,然後排序。
先根據第一關鍵字排序,第一關鍵字相等時根據第二關鍵字大小排序。
但是這段程式碼看上去並沒有這麼做,它只是做了一個計數排序而已。
還記得這個計數排序的特點:先根據x的值排序,x值相等時根據出現先後次序排序。
x裡面存了上次關鍵字的排序,在本次排序中即是第一關鍵字的排序,x的值排序==第一關鍵字排序,這裡的計數排序做的是對的。那麼第二關鍵字呢?
前面對第二關鍵字進行了排序,在這裡wv[i] = x[y[i]];根據第二關鍵字的順序重新改變了第一關鍵字的順序,也就是說在本次計數排序中,出現先後次序排序==第二關鍵字大小排序。
換句話說,我們先單獨對第二關鍵字排序,根據這個順序改變第一關鍵字的順序,由於在計數排序時首先按照第一關鍵字的值排序,而第一關鍵字的值沒有改變所以首先還是根據第一關鍵字排序,改變的是第一關鍵字相同的時候,出現在前面的第二關鍵字排在前面。
做到這裡就完成了第一第二關鍵字的合併,得到了合併以後的關鍵字順序,它可以用於下次迭代。
每次新的迭代要用到rank陣列x,由於有了關鍵字排序陣列sa,要得到rank陣列也很容易。但是對於相同的值,rank應該相同,所以要判斷一下合併以後的關鍵字是否相同。這裡都很好懂,我就不細講了。
for(t = x,x = y,y = t,p = 1,x[sa[0]] = 0,i = 1; i < n;i++)
x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
程式碼完整版:
int wa[maxn],wb[maxn],wv[maxn],ws[maxn];
int cmp(int *r , int a, int b, int l)
{
return r[a] == r[b] && r[a+l] == r[b+l];
}
void da (int *r , int *sa , int n, int m)
{
int i, j, p, *x = wa, *y = wb , *t;
for(i = 0; i < m; i++)
ws[i] = 0;
for(i = 0; i < n; i++)
ws[x[i] = r[i]]++;
for(i = 1; i < m; i++)
ws[i] += ws[i-1];
for(i = n-1; i >= 0; i--)
sa[--ws[x[i]]] = i;
for(j = 1,p = 1; p < n ; j <<= 1,m = p)
{
for(p = 0, i = n - j; i < n; i++)
y[p++]=i;
for(i = 0; i < n; i++)
if(sa[i] >= j)
y[p++] = sa[i] - j;
for(i = 0; i < n; i++)
wv[i] = x[y[i]];
for(i = 0; i < m; i++)
ws[i] = 0;
for(i = 0; i < n; i++)
ws[wv[i]]++;
for(i = 1; i < m; i++)
ws[i] += ws[i-1];
for(i = n-1; i >= 0; i--)
sa[--ws[wv[i]]] = y[i];
for(t = x,x = y,y = t,p = 1,x[sa[0]] = 0,i = 1; i < n;i++)
x[sa[i]]=cmp(y,sa[i-1],sa[i],j)?p-1:p++;
}
}