Skip to main content
2 of 3
fixed inefficiency of suffix array sorting algorithm
Dmitry
  • 4.6k
  • 1
  • 19
  • 32

The only approach I can suggest is to use an LCP array.

Since my first implementation has extremely slow sorting of suffic array, I decided to base my new solution on the Sais implementation of the induced sorting based suffix array construction algorithm. It has a complexity of \$O(n)\$.

Here is my new implementation.


###Simplified Sais algorithm implementation Creates sorted int[] suffix array.
It has unreadable code, but it is very fast.

internal interface IBaseArray
{
    int this[int i] { get; set; }
}


internal sealed class IntArray : IBaseArray
{
    private readonly int[] m_array;
    private readonly int m_pos;

    public IntArray(int[] array, int pos)
    {
        m_pos = pos;
        m_array = array;
    }

    public int this[int i]
    {
        get { return m_array[i + m_pos]; }
        set { m_array[i + m_pos] = value; }
    }
}

/// <summary>
/// An implementation of the induced sorting based suffix array construction algorithm.
/// </summary>
public static class Sais
{
    private const int MINBUCKETSIZE = 256;

    private static void getCounts(IBaseArray t, IBaseArray c, int n, int k)
    {
        int i;
        for (i = 0; i < k; ++i)
        {
            c[i] = 0;
        }
        for (i = 0; i < n; ++i)
        {
            c[t[i]] = c[t[i]] + 1;
        }
    }

    private static void getBuckets(IBaseArray c, IBaseArray b, int k, bool end)
    {
        int i, sum = 0;
        if (end)
        {
            for (i = 0; i < k; ++i)
            {
                sum += c[i];
                b[i] = sum;
            }
        }
        else
        {
            for (i = 0; i < k; ++i)
            {
                sum += c[i];
                b[i] = sum - c[i];
            }
        }
    }

    /* sort all type LMS suffixes */

    private static void LMSsort(IBaseArray t, int[] sa, IBaseArray c, IBaseArray B, int n, int k)
    {
        int i;
        int c0, c1;
        /* compute SAl */
        if (c == B)
        {
            getCounts(t, c, n, k);
        }
        getBuckets(c, B, k, false); /* find starts of buckets */
        int j = n - 1;
        int b = B[c1 = t[j]];
        --j;
        sa[b++] = (t[j] < c1) ? ~j : j;
        for (i = 0; i < n; ++i)
        {
            if (0 < (j = sa[i]))
            {
                if ((c0 = t[j]) != c1)
                {
                    B[c1] = b;
                    b = B[c1 = c0];
                }
                --j;
                sa[b++] = (t[j] < c1) ? ~j : j;
                sa[i] = 0;
            }
            else if (j < 0)
            {
                sa[i] = ~j;
            }
        }
        /* compute SAs */
        if (c == B)
        {
            getCounts(t, c, n, k);
        }
        getBuckets(c, B, k, true); /* find ends of buckets */
        for (i = n - 1, b = B[c1 = 0]; 0 <= i; --i)
        {
            if (0 < (j = sa[i]))
            {
                if ((c0 = t[j]) != c1)
                {
                    B[c1] = b;
                    b = B[c1 = c0];
                }
                --j;
                sa[--b] = (t[j] > c1) ? ~(j + 1) : j;
                sa[i] = 0;
            }
        }
    }

    private static int LMSpostproc(IBaseArray t, int[] sa, int n, int m)
    {
        int i, j, p, q;
        int qlen, name;
        int c1;

        /* compact all the sorted substrings into the first m items of SA
            2*m must be not larger than n (proveable) */
        for (i = 0; (p = sa[i]) < 0; ++i)
        {
            sa[i] = ~p;
        }
        if (i < m)
        {
            for (j = i, ++i;; ++i)
            {
                if ((p = sa[i]) < 0)
                {
                    sa[j++] = ~p;
                    sa[i] = 0;
                    if (j == m)
                    {
                        break;
                    }
                }
            }
        }

        /* store the length of all substrings */
        i = n - 1;
        j = n - 1;
        int c0 = t[n - 1];
        do
        {
            c1 = c0;
        } while ((0 <= --i) && ((c0 = t[i]) >= c1));
        for (; 0 <= i;)
        {
            do
            {
                c1 = c0;
            } while ((0 <= --i) && ((c0 = t[i]) <= c1));
            if (0 <= i)
            {
                sa[m + ((i + 1) >> 1)] = j - i;
                j = i + 1;
                do
                {
                    c1 = c0;
                } while ((0 <= --i) && ((c0 = t[i]) >= c1));
            }
        }

        /* find the lexicographic names of all substrings */
        for (i = 0, name = 0, q = n, qlen = 0; i < m; ++i)
        {
            p = sa[i];
            int plen = sa[m + (p >> 1)];
            bool diff = true;
            if ((plen == qlen) && ((q + plen) < n))
            {
                for (j = 0; (j < plen) && (t[p + j] == t[q + j]); ++j)
                {
                }
                if (j == plen)
                {
                    diff = false;
                }
            }
            if (diff)
            {
                ++name;
                q = p;
                qlen = plen;
            }
            sa[m + (p >> 1)] = name;
        }

        return name;
    }

    /* compute SA and BWT */

    private static void induceSA(IBaseArray t, int[] sa, IBaseArray c, IBaseArray B, int n, int k)
    {
        int b, i, j;
        int c0, c1;
        /* compute SAl */
        if (c == B)
        {
            getCounts(t, c, n, k);
        }
        getBuckets(c, B, k, false); /* find starts of buckets */
        j = n - 1;
        b = B[c1 = t[j]];
        sa[b++] = ((0 < j) && (t[j - 1] < c1)) ? ~j : j;
        for (i = 0; i < n; ++i)
        {
            j = sa[i];
            sa[i] = ~j;
            if (0 < j)
            {
                if ((c0 = t[--j]) != c1)
                {
                    B[c1] = b;
                    b = B[c1 = c0];
                }
                sa[b++] = ((0 < j) && (t[j - 1] < c1)) ? ~j : j;
            }
        }
        /* compute SAs */
        if (c == B)
        {
            getCounts(t, c, n, k);
        }
        getBuckets(c, B, k, true); /* find ends of buckets */
        for (i = n - 1, b = B[c1 = 0]; 0 <= i; --i)
        {
            if (0 < (j = sa[i]))
            {
                if ((c0 = t[--j]) != c1)
                {
                    B[c1] = b;
                    b = B[c1 = c0];
                }
                sa[--b] = ((j == 0) || (t[j - 1] > c1)) ? ~j : j;
            }
            else
            {
                sa[i] = ~j;
            }
        }
    }

    private static int computeBWT(IBaseArray t, int[] sa, IBaseArray c, IBaseArray B, int n, int k)
    {
        int b, i, j, pidx = -1;
        int c0, c1;
        /* compute SAl */
        if (c == B)
        {
            getCounts(t, c, n, k);
        }
        getBuckets(c, B, k, false); /* find starts of buckets */
        j = n - 1;
        b = B[c1 = t[j]];
        sa[b++] = ((0 < j) && (t[j - 1] < c1)) ? ~j : j;
        for (i = 0; i < n; ++i)
        {
            if (0 < (j = sa[i]))
            {
                sa[i] = ~(c0 = t[--j]);
                if (c0 != c1)
                {
                    B[c1] = b;
                    b = B[c1 = c0];
                }
                sa[b++] = ((0 < j) && (t[j - 1] < c1)) ? ~j : j;
            }
            else if (j != 0)
            {
                sa[i] = ~j;
            }
        }
        /* compute SAs */
        if (c == B)
        {
            getCounts(t, c, n, k);
        }
        getBuckets(c, B, k, true); /* find ends of buckets */
        for (i = n - 1, b = B[c1 = 0]; 0 <= i; --i)
        {
            if (0 < (j = sa[i]))
            {
                sa[i] = (c0 = t[--j]);
                if (c0 != c1)
                {
                    B[c1] = b;
                    b = B[c1 = c0];
                }
                sa[--b] = ((0 < j) && (t[j - 1] > c1)) ? ~t[j - 1] : j;
            }
            else if (j != 0)
            {
                sa[i] = ~j;
            }
            else
            {
                pidx = i;
            }
        }
        return pidx;
    }

    /* find the suffix array SA of T[0..n-1] in {0..k-1}^n
       use a working space (excluding T and SA) of at most 2n+O(1) for a constant alphabet */

    private static int sais_main(IBaseArray t, int[] sa, int fs, int n, int k, bool isbwt)
    {
        IBaseArray c, B;
        int i;
        int name, pidx = 0;
        int c1;
        uint flags;

        if (k <= MINBUCKETSIZE)
        {
            c = new IntArray(new int[k], 0);
            if (k <= fs)
            {
                B = new IntArray(sa, n + fs - k);
                flags = 1;
            }
            else
            {
                B = new IntArray(new int[k], 0);
                flags = 3;
            }
        }
        else if (k <= fs)
        {
            c = new IntArray(sa, n + fs - k);
            if (k <= (fs - k))
            {
                B = new IntArray(sa, n + fs - k * 2);
                flags = 0;
            }
            else if (k <= (MINBUCKETSIZE * 4))
            {
                B = new IntArray(new int[k], 0);
                flags = 2;
            }
            else
            {
                B = c;
                flags = 8;
            }
        }
        else
        {
            c = B = new IntArray(new int[k], 0);
            flags = 4 | 8;
        }

        /* stage 1: reduce the problem by at least 1/2
           sort all the LMS-substrings */
        getCounts(t, c, n, k);
        getBuckets(c, B, k, true); /* find ends of buckets */
        for (i = 0; i < n; ++i)
        {
            sa[i] = 0;
        }
        int b = -1;
        i = n - 1;
        int j = n;
        int m = 0;
        int c0 = t[n - 1];
        do
        {
            c1 = c0;
        } while ((0 <= --i) && ((c0 = t[i]) >= c1));
        for (; 0 <= i;)
        {
            do
            {
                c1 = c0;
            } while ((0 <= --i) && ((c0 = t[i]) <= c1));
            if (0 <= i)
            {
                if (0 <= b)
                {
                    sa[b] = j;
                }
                b = --B[c1];
                j = i;
                ++m;
                do
                {
                    c1 = c0;
                } while ((0 <= --i) && ((c0 = t[i]) >= c1));
            }
        }
        if (1 < m)
        {
            LMSsort(t, sa, c, B, n, k);
            name = LMSpostproc(t, sa, n, m);
        }
        else if (m == 1)
        {
            sa[b] = j + 1;
            name = 1;
        }
        else
        {
            name = 0;
        }

        /* stage 2: solve the reduced problem
           recurse if names are not yet unique */
        if (name < m)
        {
            if ((flags & 4) != 0)
            {
                c = null;
                B = null;
            }
            if ((flags & 2) != 0)
            {
                B = null;
            }
            int newfs = (n + fs) - (m * 2);
            if ((flags & (1 | 4 | 8)) == 0)
            {
                if ((k + name) <= newfs)
                {
                    newfs -= k;
                }
                else
                {
                    flags |= 8;
                }
            }
            for (i = m + (n >> 1) - 1, j = m * 2 + newfs - 1; i >= m; --i)
            {
                if (sa[i] != 0)
                {
                    sa[j--] = sa[i] - 1;
                }
            }
            IBaseArray ra = new IntArray(sa, m + newfs);
            sais_main(ra, sa, newfs, m, name, false);

            i = n - 1;
            j = m * 2 - 1;
            c0 = t[n - 1];
            do
            {
                c1 = c0;
            } while ((0 <= --i) && ((c0 = t[i]) >= c1));
            for (; 0 <= i;)
            {
                do
                {
                    c1 = c0;
                } while ((0 <= --i) && ((c0 = t[i]) <= c1));
                if (0 <= i)
                {
                    sa[j--] = i + 1;
                    do
                    {
                        c1 = c0;
                    } while ((0 <= --i) && ((c0 = t[i]) >= c1));
                }
            }

            for (i = 0; i < m; ++i)
            {
                sa[i] = sa[m + sa[i]];
            }
            if ((flags & 4) != 0)
            {
                c = B = new IntArray(new int[k], 0);
            }
            if ((flags & 2) != 0)
            {
                B = new IntArray(new int[k], 0);
            }
        }

        /* stage 3: induce the result for the original problem */
        if ((flags & 8) != 0)
        {
            getCounts(t, c, n, k);
        }
        /* put all left-most S characters into their buckets */
        if (1 < m)
        {
            getBuckets(c, B, k, true); /* find ends of buckets */
            i = m - 1;
            j = n;
            int p = sa[m - 1];
            c1 = t[p];
            do
            {
                int q = B[c0 = c1];
                while (j > q)
                {
                    sa[--j] = 0;
                }
                do
                {
                    sa[--j] = p;
                    if (--i < 0)
                    {
                        break;
                    }
                    p = sa[i];
                } while ((c1 = t[p]) == c0);
            } while (0 <= i);

            while (0 < j)
            {
                sa[--j] = 0;
            }
        }
        if (isbwt == false)
        {
            induceSA(t, sa, c, B, n, k);
        }
        else
        {
            pidx = computeBWT(t, sa, c, B, n, k);
        }
        return pidx;
    }


    /// <summary>
    /// Constructs the suffix array of a given sequence in linear time.
    /// </summary>
    /// <param name="t">input sequence</param>
    /// <param name="k">alphabet size</param>
    /// <returns>output suffix array</returns>
    public static int[] sufsort(int[] t, int k)
    {
        if (t == null)
            throw new ArgumentNullException("t");

        if (k <= 0)
            throw new ArgumentOutOfRangeException("k");

        // Length of the given string
        int n = t.Length;

        // Output suffix array
        int[] sa = new int[n];

        return n <= 1 || sais_main(new IntArray(t, 0), sa, 0, n, k, false) == 0 ? sa : null;
    }
}

###LongestCommonPhraseInfo class A class that will hold a result

public sealed class LongestCommonPhraseInfo
{
    public readonly string[] CommonPhraseWords;
    public readonly int WordIndex1;
    public readonly int WordIndex2;

    public LongestCommonPhraseInfo(string[] commonPhraseWords, int wordIndex1, int wordIndex2)
    {
        CommonPhraseWords = commonPhraseWords;
        WordIndex1 = wordIndex1;
        WordIndex2 = wordIndex2;
    }

    public string Phrase
    {
        get { return String.Join(" ", CommonPhraseWords); }
    }
}

###The main class Calls the Sais algorithm class to create a suffix array and creates the corresponding LCP array. It contains a public method GetLongestCommonPhrase and a tiny utility method Between<T>:

/// <summary>
/// Checks if <paramref name="x"/> value is between values 
/// <paramref name="a"/> and <paramref name="b"/> (or 
/// <paramref name="b"/> and <paramref name="a"/>).
/// </summary>
/// <typeparam name="T">Type of arguments.</typeparam>
private static bool Between<T>(T x, T a, T b)
    where T : IComparable<T>
{
    return a.CompareTo(b) <= 0
               ? x.CompareTo(a) >= 0 && x.CompareTo(b) <= 0
               : x.CompareTo(a) <= 0 && x.CompareTo(b) >= 0;
}

public static LongestCommonPhraseInfo GetLongestCommonPhrase(string text1, string text2)
{
    const string Sentinel1 = "\x00";
    const string Sentinel2 = "\x01";

    // Split texts.
    string[] textWords1 = Regex.Split(text1, @"\s+|-+", RegexOptions.Compiled | RegexOptions.Singleline);
    string[] textWords2 = Regex.Split(text2, @"\s+|-+", RegexOptions.Compiled | RegexOptions.Singleline);

    // Combine texts into single array.
    string[] textWords = new string[textWords1.Length + textWords2.Length + 2];
    Array.Copy(textWords1, textWords, textWords1.Length);
    textWords[textWords1.Length] = Sentinel1;
    Array.Copy(textWords2, 0, textWords, textWords1.Length + 1, textWords2.Length);
    textWords[textWords.Length - 1] = Sentinel2;

    // Get distinct words of text.
    string[] alphabet = textWords.Distinct().ToArray();

    // Create temp dictionary.
    // Key: alphabet element.
    // Value: index in alphabet.
    Dictionary<string, int> wordsIndex = new Dictionary<string, int>(alphabet.Length);
    for (int i = 0; i < alphabet.Length; i++)
        wordsIndex[alphabet[i]] = i;

    // Convert each word of the text to its index.
    int[] indexedText = Array.ConvertAll(textWords, w => wordsIndex[w]);

    // Call the Sais algorithm to create int[] suffix array.
    int[] sa = Sais.sufsort(indexedText, alphabet.Length);
    if (sa == null)
        return null;

    // If succeededб create LCP array.
    int[] lcps = new int[sa.Length];
    int prev = sa[0];
    int maxLcp = -1;
    int maxLcpIndex = -1;

    for (int i = 1; i < lcps.Length; i++)
    {
        int cur = sa[i];

        if (Between(textWords1.Length, prev, cur))
        {
            // If `prev` and `cur` suffixes belong to different text parts.

            // Calculate the LCP.
            int lcp = 0;
            while (indexedText[prev + lcp] == indexedText[cur + lcp])
            {
                lcp++;
            }
            // Find the maximum LCP.
            if (lcp > maxLcp)
            {
                maxLcp = lcp;
                maxLcpIndex = i;
            }
        }

        prev = cur;
    }

    int[] tmp = new int[maxLcp];
    Array.Copy(indexedText, sa[maxLcpIndex], tmp, 0, maxLcp);

    return new LongestCommonPhraseInfo(Array.ConvertAll(tmp, i => alphabet[i]),
        Math.Min(sa[maxLcpIndex], sa[maxLcpIndex - 1]),
        Math.Max(sa[maxLcpIndex], sa[maxLcpIndex - 1]));
}

Usage sample:

var book1 = File.ReadAllText("50503-0.txt").ToLower();
var book2 = File.ReadAllText("50511-0.txt").ToLower();

    // Read the books.
    string book1 = File.ReadAllText(@"D:\50503-0.txt").ToLower();
    string book2 = File.ReadAllText(@"D:\50511-0.txt").ToLower();

    // Make the books longer.
    string text1 = String.Concat(book1, " word1 ", book1, " word2 ", book1, " word3 ", book1);
    string text2 = String.Concat(book2, " word4 ", book2, " word5 ", book2, " word6 ", book2);

    // Go!
    Stopwatch sw = new Stopwatch();
    sw.Start();

    // Find the longest phrase.
    var longestPhrase = GetLongestCommonPhrase(text1, text2);

    sw.Stop();

    if (longestPhrase != null)
    {
        Console.WriteLine("Common Phrase Words Count: {0}", longestPhrase.CommonPhraseWords.Length);
        Console.WriteLine("Common Phrase: \"{0}...\"", longestPhrase.Phrase.Substring(0, 64));
        Console.WriteLine("Time elapsed (ms): {0}", sw.ElapsedMilliseconds);
    }
    else
    {
        Console.WriteLine("Failed!");
    }

Output for enlarged texts (on my PC):

Common Phrase Words Count: 2999
Common Phrase: "updated editions will replace the previous one the old editions ..."
Time elapsed (ms): 4537

For non-enlarged texts:

Time elapsed (ms): 1241
Dmitry
  • 4.6k
  • 1
  • 19
  • 32