1. 程式人生 > >POJ2778 DNA Sequence 題解(AC自動機+矩陣快速冪)

POJ2778 DNA Sequence 題解(AC自動機+矩陣快速冪)

(題目描述略)

本題對於時間的要求比較嚴格,這意味著樸素的搜尋演算法是很難用優化手段通過的。為了解決這個問題,我們需要用一種不同於列舉的演算法。

幾乎任何一個問題都有其對應的圖論模型,這個問題也不例外。我們可以將此問題轉化為在其對應的圖上求可接受的長度確定的路徑條數。由於本題只有 AGCT 四個選擇,所以其規模是我們可以接受的。

為了將這個問題轉化為圖論,我們首先要將其等價的圖建立起來。根據定義,一個圖由點集和邊集組成。對於本題,我們若已知序列 DNA[1 .. n],欲推匯出 DNA[1 .. n+1],我們可以將 DNA[1 .. n] 和 DNA[1 .. n+1] 作為點,即狀態,將 DNA[n+1] 作為邊,即狀態的轉移,記權值為 1 。當一個狀態序列中包含禁阻子串時,我們將這個狀態稱為禁阻狀態,即不合法的狀態,而所有從其他狀態轉移到此狀態(此點的入邊)和從此狀態轉移到其他狀態(此點的出邊)都稱為禁阻轉移,即不合法的轉移。我們將所有禁阻狀態和禁阻轉移從我們所構建的圖中移去,則剩餘圖的所有狀態必全為可接受的狀態,即有效狀態。我們再在此圖中統計給定長度的路徑條數,即為所求。

但是,以上做法不僅沒有解決時間的問題,還增加了空間的問題。為了解決空間不足的問題,我們可以去除不需要的狀態。設存在一個由有效狀態 DNA[1 .. n] 轉移到另一個有效狀態 DNA[1 .. n+1],且 n 大於最長的禁阻子串的長度,則在有效狀態 DNA[1 .. n+1] 中 DNA[1] 為多餘狀態,因為其既不參與之後的轉移,又不可能與之後的字元構成禁阻子串。推廣開來,我們可以預處理出最長的禁阻子串的長度 length,則狀態序列的最長長度即為 length - 1,當由有效狀態 DNA[n-length+2 .. n] 轉移到狀態 DNA[n-length+3 .. n+1] 時,我們只需判斷狀態序列 DNA[n-length+3 .. n+1] 是否包含禁阻子串,若否則其必為有效狀態。

儘管有上文的空間優化,其空間複雜度仍為指數級別。為了進一步優化空間,我們需要儘可能地壓縮狀態。由於禁阻子串的個數很小,我們可以發現在狀態序列的某幾位上某幾個字元是幾乎等價的。我們將這些狀態合併,便可以使狀態總數維持在可接受的範圍之內。但是,這一步說起來容易,做起來卻是比較複雜的。為了更高效地完成這一操作,我們將使用一種線性時間複雜度和線性空間複雜度的資料結構,即本題的正解 AC自動機

AC自動機是一種可以接受給定字串集合內所有字串的確定有窮狀態自動機,可以線上性時間內構造。其具體實現方法在這裡不做過多介紹。根據禁阻序列所構造的AC自動機其實已經完成了狀態的壓縮,我們只需刪去所有禁阻狀態和禁阻轉移便可實現空間的優化。關於以上結論的解釋如下,因為AC自動機中,我們每插入一個節點,便對其父節點進行了一次狀態的劃分,即將父節點所對應的剩餘有效狀態劃分成通向新禁阻狀態的轉移和通向新有效狀態的轉移,而在當前狀態下所有通向有效狀態的轉移是等價的。即使一個節點的狀態轉移被劃分多次,剩餘有效狀態轉移仍然等價。例如對於禁阻序列 {“AT”,”AC”,”AG”,”AA”},當插入 “AT” 和 “AC” 完成後,狀態 “A” 被劃分成 ‘T’ 轉移和 ‘C’ 轉移,而此時狀態 “A” 的 ‘G’ 轉移和 ‘A’ 轉移是等價的。基於以上原理,一個和問題等價的圖論模型就建立起來了。

接下來我們要解決的是時間的問題。圖建立好之後,我們要使用一種高效的手段計算以源點為路徑起點的長度為 n 的不同路徑的總數。在這裡,我們仍然用狀態轉移的思想進行分析。我們記路徑長度為 i,起始點 j,終止點 k 的不同路徑總數為 f[i, j, k], 則 f[i, j, k] = sum(f[i-1, j, l] * f[1, l, k]) (1 ≤ l ≤ m,m 為圖中點的總數)(很像 Floyd)。翻譯成中文,意思即為將長度為 i 的路徑劃分為長度為 i - 1 的路徑和長度為 1 的路徑, 列舉劃分點求和。當然,你也可以將其劃分為長度 i - 2 的路徑和長度 2 的路徑,或長度為 i - 3 的路徑和長度 3 的路徑,等等。 然而為了更好地探究其性質,劃分成長度 i - 1 和長度 1 的路徑方便一些。在列出上文的動態轉移方程之後,我們驚奇地發現這與一個數據結構的運算非常相似,這個資料結構就是 矩陣,而這種運算就是 矩陣乘法

關於矩陣及矩陣乘法,我在這裡也不再做過多的介紹,只說明一點,那就是矩陣乘法的運演算法則:

ApqBqr=CprC(i,j)=k=1qA(i,k)B(k,j)同時,矩陣乘法不滿足交換律,滿足結合律,這使得矩陣乘法可以用 矩陣快速冪 進行運算。我們將圖用鄰接矩陣的形式儲存,記為 A,即 A[i, j] 表示從 i 到 j 的路徑條數,初始路徑長度為 1,將 A 看做矩陣,mi=1An(1,i)即為所求。這樣,就使得時間複雜度也降低到了我們可以接受的範圍之內。

最後還有幾個小問題。其一是結構體賦值的問題。對於矩陣乘法的運算,由於在矩陣快速冪中不止出現一次,最好寫子程式或過載乘法運算子。這樣不可避免地出現用結構體定義矩陣的情況,而結構體賦值便成了一個糾結的問題。其實,在 C/CPP 中,結構體賦值是允許的,其原理應該是逐位賦值。於是,過載賦值運算子便成了不必要的了。

其二是整型溢位的問題。雖然題目描述中提到最後結果對 100,000 取模,但事實上矩陣乘法運算過程中仍有整型溢位的可能。為解決這個問題,long long 超長整型需要在運算過程中使用。

其三是矩陣 1 的問題。對於矩陣乘法,我們需要定義一個初始矩陣 C,使得任意矩陣 A * C = A,即類似於數 1 的效果。這樣的矩陣是存在的,其為:

C(i,j)={10i = ji ≠ j

程式碼如下:

#include"stdio.h"
#include"string.h"
bool end[105];
char s[15];
int f[105],fail[105],m,next[105][4],q[105],top=1;
struct Matrix
{
    long long data[105][105];
}mat;
int reflect(char c)
{
    switch(c)
    {
        case 'A':return 0;
        case 'C':return 1;
        case 'G':return 2;
        case 'T':return 3;
    }
}
void Insert(char str[])
{
    int now=0;
    for(int i=0;str[i];now=next[now][reflect(str[i++])])
        if(next[now][reflect(str[i])]==-1)
            next[now][reflect(str[i])]=top++;
    end[now]=false;
}
void Build()
{
    top=1;
    for(int i=0;i<4;i++)
        if(next[0][i]>-1)
            fail[next[0][i]]=0,q[top++]=next[0][i];
        else
            next[0][i]=0;
    for(int i=1;i<top;i++)
        for(int j=0;j<4;j++)
            if(next[q[i]][j]>-1)
                fail[next[q[i]][j]]=next[fail[q[i]]][j],q[top++]=next[q[i]][j];
            else
                next[q[i]][j]=next[fail[q[i]]][j];
}
Matrix Multiply(Matrix A,Matrix B)
{
    Matrix C;
    for(int i=0;i<m;i++)
        for(int j=0;j<m;j++)
        {
            C.data[i][j]=0;
            for(int k=0;k<m;k++)
                C.data[i][j]=(A.data[i][k]*B.data[k][j]%100000+C.data[i][j])%100000;
        }
    return C;
}
Matrix pow(Matrix A,int n)
{
    Matrix C;
    for(int i=0;i<m;i++)
        for(int j=0;j<m;j++)
            C.data[i][j]=(i==j)?1:0;
    for(Matrix B=A;n>0;n>>=1,B=Multiply(B,B))
        if(n&1)
            C=Multiply(B,C);
    return C;
}
int add(int n)
{
    int sum=0;
    mat=pow(mat,n);
    for(int i=0;i<m;i++)
        sum=(mat.data[0][i]+sum)%100000;
    return sum;
}
int main()
{
    int n;
    scanf("%d %d\n",&m,&n);
    memset(end,true,sizeof(end));
    memset(next,-1,sizeof(next));
    while(m--)
        gets(s),Insert(s);
    Build();
    m=0;
    memset(f,-1,sizeof(f));
    for(int i=0;i<top;i++)
    {
        for(int temp=i;temp>0&&end[i];temp=fail[temp])
            end[i]&=end[temp];
        if(end[i])
            f[i]=m++;
    }
    memset(mat.data,0,sizeof(mat.data));
    for(int i=0;i<top;i++)
        if(end[i])
            for(int j=0;j<4;j++)
                if(f[next[i][j]]>-1)
                    mat.data[f[i]][f[next[i][j]]]++;
    printf("%d",add(n));
    return 0;
}

對於樣例 {“AT”,”AC”,”AG”,”AA”},所建矩陣

(3010)矩陣快速冪後結果(27090)

本題非常經典,建議編寫。由於本題樣例過弱,不建議用樣例檢驗程式的正確性。