Suffix Trees are very useful in numerous string processing and computational biology problems.

The task is to create a function which implements Ukkonen’s algorithm to create a useful Suffix Tree as described:

Part 1
Part 2
Part 3
Part 4
Part 5
Part 6

Using Arithmetic-geometric mean/Calculate Pi generate the first 1000, 10000, and 100000 decimal places of pi. Using your implementation with an alphabet of 0 through 9 (plus $ say to make the tree explicit) find the longest repeated string in each list. Time your results and demonstrate that your implementation is linear (i.e. that 10000 takes approx. 10 times as long as 1000). You may vary the size of the lists of decimal places of pi to give reasonable answers.

-- demo/rosetta/Ukkonens_Suffix_Tree.exw
with javascript_semantics
integer maxChar = 'z'
 
sequence children = {},
         suffixLinks = {},
         starts = {},
         endIndices = {},
         suffixIndices = {},
         leaves = {}

function new_leaf(integer v=0)
    leaves = append(leaves,v)
    return length(leaves)
end function

string text
integer splitEndIdx,
        rootEndIdx,
        leafEndIdx = new_leaf(),
        root = NULL,
        lastNewNode,
        activeNode,
        activeEdge = -1,
        activeLength = 0,
        remainingSuffixCount = 0,
        size = -1
 
function newNode(integer start, finishIdx, bool bKids=false)
    children = append(children,iff(bKids?repeat(NULL,maxChar):0))
    suffixLinks = append(suffixLinks,root)
    starts = append(starts,start)
    endIndices = append(endIndices,finishIdx)
    suffixIndices = append(suffixIndices,-1)
    return length(children)
end function
 
function edgeLength(integer n)
    return iff(n==root?0:leaves[endIndices[n]] - starts[n] + 1)
end function
 
function walkDown(integer currNode)
    integer l = edgeLength(currNode)
    if activeLength >= l then
        activeEdge += l
        activeLength -= l
        activeNode = currNode
        return true
    end if
    return false
end function
 
procedure extendSuffixTree(integer pos)
    leaves[leafEndIdx] = pos
    remainingSuffixCount += 1
    lastNewNode = NULL
    while remainingSuffixCount > 0 do
        if activeLength == 0 then
            activeEdge = pos
        end if
        integer ta = text[activeEdge]
        bool ca0 = children[activeNode]=0
        integer next = iff(ca0?NULL:children[activeNode][ta])
        if next==null then
            if ca0 then
                children[activeNode] = repeat(NULL,maxChar)
            end if
            children[activeNode][ta] = newNode(pos, leafEndIdx)
            if lastNewNode!=NULL then
                suffixLinks[lastNewNode] = activeNode
                lastNewNode = NULL
            end if
        else
            if walkDown(next) then
                continue
            end if
            integer tp = text[pos]
            if text[starts[next]+activeLength] == tp then
                if lastNewNode!=NULL and activeNode!=root then
                    suffixLinks[lastNewNode] = activeNode
                    lastNewNode = NULL
                end if
                activeLength += 1
                exit
            end if
            integer temp = starts[next] + activeLength - 1
            splitEndIdx = new_leaf(temp)
            integer splitnode = newNode(starts[next], splitEndIdx, true)
            ta = text[activeEdge]
            children[activeNode][ta] = splitnode
            children[splitnode][tp] = newNode(pos, leafEndIdx)
            starts[next] += activeLength
            children[splitnode][text[starts[next]]] = next
            if lastNewNode!=NULL then
                suffixLinks[lastNewNode] = splitnode
            end if
            lastNewNode = splitnode
        end if
        remainingSuffixCount -= 1
        if activeNode==root and activeLength>0 then
            activeLength -= 1
            activeEdge = pos - remainingSuffixCount + 1
        elsif activeNode!=root then
            activeNode = suffixLinks[activeNode]
        end if
    end while
end procedure
 
procedure setSuffixIndexByDFS(integer n, labelHeight)
    if n!=NULL then
        if children[n]=0 then
            suffixIndices[n] = size - labelHeight
        else
            bool leaf = true
            for i=1 to maxChar do
                integer nci = children[n][i]
                if nci!=NULL then
                    leaf = false
                    setSuffixIndexByDFS(nci, labelHeight+edgeLength(nci))
                end if
            end for
            if leaf then ?9/0 end if -- (sanity check)
        end if
    end if
end procedure
 
procedure buildSuffixTree()
    size = length(text)
    rootEndIdx = new_leaf(-1)
    root = newNode(-1, rootEndIdx)
    activeNode = root
    for i=1 to size do
        extendSuffixTree(i)
    end for
    integer labelHeight = 0
    setSuffixIndexByDFS(root, labelHeight)
end procedure
 
procedure doTraversal(integer n, labelHeight, maxHeightIdx, substringStartIndex)
    if n!=NULL then
        integer nsi = suffixIndices[n], newHeight
        if nsi == -1 then
            for i=1 to maxChar do
                integer nci = children[n][i]
                if nci!=NULL then
                    newHeight = labelHeight+edgeLength(nci)
                    doTraversal(nci, newHeight, maxHeightIdx, substringStartIndex)
                end if
            end for
        elsif nsi > -1 then
            newHeight = labelHeight-edgeLength(n)
            if leaves[maxHeightIdx]<newHeight then
                leaves[maxHeightIdx] = newHeight
                leaves[substringStartIndex] = nsi
            end if
        end if
    end if
end procedure
 
function getLongestRepeatedSubstring()
    integer maxHeightIdx = new_leaf(),
            substringStartIndex = new_leaf()
    doTraversal(root, 0, maxHeightIdx, substringStartIndex)
    integer maxHeight = leaves[maxHeightIdx],
            start = leaves[substringStartIndex]
    string t = iff(maxHeight=0?"No repeated substring"
                              :text[start+1..start+maxHeight])
    return t
end function
 
constant tests = {"CAAAABAAAABD$",
                  "GEEKSFORGEEKS$",
                  "AAAAAAAAAA$",
                  "ABCDEFG$",
                  "ABABABA$",
                  "ATCGATCGA$",
                  "banana$",
                  "abcpqrabpqpq$",
                  "pqrpqpqabab$"}
printf(1,"Longest Repeated Substring in:\n")
for i=1 to length(tests) do
    text = tests[i]
    buildSuffixTree()
    printf(1,"  %s is: %s\n", {text,getLongestRepeatedSubstring()})
end for
printf(1,"\n")
 
include mpfr.e
string piStr
if platform()=JS then
    mpfr pi = mpfr_init(0,-1001) -- (set precision to 1,000 dp, plus the "3.")
    mpfr_const_pi(pi)
    piStr = mpfr_get_fixed(pi,1000) -- (all we can really manage under pwa/p2js)
else
    -- gmp crashes when I try 100,000 dp, so just get from file
    piStr = get_text(`E:\downloads\misc\arm\pi-10million.txt`)
end if
piStr = piStr[3..$] -- discard leading "3."
maxChar = '9'
for i=3 to iff(platform()=JS?3:6) do
    atom t0 = time()
    integer n = power(10,i)
    text = piStr[1..n] & "$"
    buildSuffixTree()
    string r = getLongestRepeatedSubstring(),
           e = elapsed(time()-t0)
    printf(1,"  first %,d d.p. of Pi is: %s (%s)\n", {n,r,e})
end for
