Warshall‘s algorithm 演算法的實現及優化（修改版）

用 Warshall’s 演算法計算傳遞閉包

• 離散數學定義：

t® = R u R^2 u R^3 u… 其中R^(n+1) = R^n 複合 R

矩陣表示：

M（R） = M + M^2 + M^3 +…+M^n(其中加為邏輯加）

所以我們只要按照這個公式每次更新M，最後的Mn就是傳遞閉包。

• Warshall演算法：

(1)置新矩陣A＝M；
(2)i＝1；
(3)對所有j如果A[j，i]＝1，則對k＝1，2，…，n，A[j，k]＝A[j，k]∨A[i，k]；
(4)i加1；（i是行，j是列）
(5)如果i≤n，則轉到步驟3)，否則停止。

• 時間複雜度為：O(nnn)
• 用R的無窮閉包時間複雜度為O(nnn*(n - 1))

對於每個相通的j - > i,我們可以從這個相通關係出發，看看能不能通過這條相通的j - > i,更新一下j - >k。對所有的可通關係都更新一遍M，最後的結果就是傳遞閉包。

實現及優化：

下面展示一些 經典程式碼片。

void computeAPSP(const int n) {
    /* calculate shortest paths from every vertex to every vertex */
    for (int k = 0; k < n; k++) 
    {
        for (int i = 0; i < n; i++) 
        {
            for (int j = 0; j < n; j++) 
            {
                a[i][j] = min( a[i][j], a[i][k] + a[k][j] );
            }
        }
    }
}

利用矩陣的對稱性優化：

void computeAPSP(const int n) 
{
    for (int k = 0; k < n; k++) 
    {
        for (int i = 0; i < n; i++) 
        {
            if (k != i) 
            {

                const int a_ki = (k < i) ? a[i][k] : a[k][i];

                for (int j = 0; j < min(k, i); j++)
                    a[i][j] = min( a[i][j], a_ki + a[k][j] );

                for (int j = k + 1; j < i; j++)
                    a[i][j] = min( a[i][j], a_ki + a[j][k] );

            }
        }
    }
}

只使用矩陣的下三角部分進行優化：

void computeAPSP(const int n) 
{
    for (int k = 0; k < n; k++) 
    {
        for (int i = 0; i < n; i++) 
        {
            if (k != i) 
            {

                const int a_ki = (k < i) ? a[i][k] : a[k][i];

                for (int j = 0; j < min(k, i); j++)
                    a[i][j] = min( a[i][j], a_ki + a[k][j] );

                for (int j = k + 1; j < i; j++)
                    a[i][j] = min( a[i][j], a_ki + a[j][k] );

            }
        }
    }
}

跳過不存在的路徑的優化：

void computeAPSP(const int n) 
{
    for (int k = 0; k < n; k++) 
    {
        for (int i = 0; i < n; i++) 
        {
            if (k != i) 
            {

                const int a_ki = (k < i) ? a[i][k] : a[k][i];

                // skip if no path
                if (a_ki == POSITIVE_INFINITY)    continue;

                for (int j = 0; j < min(k, i); j++)
                    a[i][j] = min( a[i][j], a_ki + a[k][j] );

                for (int j = k + 1; j < i; j++)
                    a[i][j] = min( a[i][j], a_ki + a[j][k] );

            }
        }
    }
}

避免大量呼叫數學函式進行優化：

void computeAPSP(const int n) 
{
    for (int k = 0; k < n; k++) 
    {
        for (int i = 0; i < n; i++) 
        {
            if (k != i) 
            {

                const int a_ki = (k < i) ? a[i][k] : a[k][i];

                // skip if no path
                if (a_ki == POSITIVE_INFINITY)    continue;

                for (int j = 0; j < min(k, i); j++) 
                {
                    const int s_kj = a_ki + a[k][j];
                    if( s_kj < a[i][j] )    a[i][j] = s_kj;
                }

                for (int j = k + 1; j < i; j++) 
                {
                    const int s_jk = a_ki + a[j][k];
                    if( s_jk < a[i][j] )    a[i][j] = s_jk;
                }

            }
        }
    }
}