Solution for Range Minimum Query & Lowest Common Ancestor

req-lca.cs

/*
    Solves both Range Minimum Query & Lowest Common Ancestor
    in O(n) space & O(1) query time following Farach-Colton and Bender's algorithm.
    Reference: https://www.ics.uci.edu/~eppstein/261/BenFar-LCA-00.pdf
*/
using System.Collections.Generic;
using System.Linq;

public class RMQLin
{
    IList<int>[] cartesianTree;
    IList<int> eulerTour;
    IList<int> firstVisit; // representatives; holds indices in eulerTour
    IList<int> depths; // depths of the nodes
    IList<int> logTable; // precomputed logarithms
    int[,] sparseTable;
    int blockSize;
    int blockCount;
    int[][,] blocks;
    IList<int> blockMasks;

    public RMQLin(int[] items)
    {
        int root = this.ConstructCartesianTree(items);
        this.PrecomputeLca(root, items.Length);
    }

    public int RMQ(int i, int j)
    {
        int result = this.LCA(i, j);
        return result;
    }

    int ConstructCartesianTree(IList<int> A)
    {
        var parentChildMap = new int[A.Count];

        for (int i = 0; i < parentChildMap.Length; i++)
            parentChildMap[i] = -1;

        // Stores the indices of the rightmost nodes.
        var stack = new Stack<int>();

        for (int i = 0; i < A.Count; i++)
        {
            int last = -1;

            // Move through the right spine towards the root
            // until finding an element smaller than the current one
            while (stack.Any() && A[stack.Peek()] >= A[i])
                last = stack.Pop();

            if (stack.Any())
                parentChildMap[i] = stack.Peek();

            if (last >= 0)
                parentChildMap[last] = i;

            stack.Push(i);
        }

        int root = -1;

        this.cartesianTree = new List<int>[A.Count];
        for (int i = 0; i < parentChildMap.Length; i++)
        {
            if (parentChildMap[i] == -1)
            {
                root = i;
                continue;
            }

            if (this.cartesianTree[parentChildMap[i]] == null)
                this.cartesianTree[parentChildMap[i]] = new List<int>();

            this.cartesianTree[parentChildMap[i]].Add(i);
        }

        return root;
    }

    void DepthFirstSearch(int current, int nodeDepth)
    {
        this.firstVisit[current] = this.eulerTour.Count;
        this.eulerTour.Add(current);
        this.depths[current] = nodeDepth;

        // Reached a leaf node
        if (cartesianTree[current] == null)
            return;

        foreach (int child in this.cartesianTree[current])
        {
            this.DepthFirstSearch(child, nodeDepth + 1);
            this.eulerTour.Add(current);
        }
    }

    int MinByDepth(int i, int j)
    {
        return this.depths[this.eulerTour[i]] <= this.depths[this.eulerTour[j]] ? i : j;
    }

    void PrecomputeLca(int root, int itemCount)
    {
        // Get Euler Tour Indices of first occurrances
        this.firstVisit = new int[itemCount];
        this.depths = new int[itemCount];
        this.eulerTour = new List<int>(itemCount * 2);
        this.DepthFirstSearch(root, 0);

        // Precompute all log values
        int eulerTourLen = this.eulerTour.Count;
        this.logTable = new List<int>();
        this.logTable.Add(-1);

        for (int i = 1; i <= eulerTourLen; i++)
            this.logTable.Add(logTable[i / 2] + 1);

        // Determine block size
        this.blockSize = this.logTable[eulerTourLen] / 2;
        if (this.blockSize < 1)
            this.blockSize = 1;

        this.blockCount = (eulerTourLen + this.blockSize - 1) / this.blockSize;

        // Precompute minimum of each block and build sparse table
        this.sparseTable = new int[this.blockCount, this.logTable[this.blockCount] + 1];

        // Inital Step
        for (int i = 0, positionInBlock = 0, blockIndex = 0; i < eulerTourLen; i++, positionInBlock++)
        {
            if (positionInBlock == this.blockSize)
            {
                positionInBlock = 0;
                blockIndex++;
            }

            // The first element of the block is the initial minimum
            if (positionInBlock == 0)
                this.sparseTable[blockIndex, 0] = i;

            // A new minimum is found. Overwrite the existing.
            if (MinByDepth(i, this.sparseTable[blockIndex, 0]) == i)
                this.sparseTable[blockIndex, 0] = i;
        }

        // Fill in the dynamic programming table
        for (int j = 1; j <= this.logTable[this.blockCount]; j++)
        {
            for (int i = 0; i + (1 << j) <= this.blockCount; i++)
            {
                int leftIndex = this.sparseTable[i, j - 1];
                int rightIndex = this.sparseTable[i + (1 << (j - 1)), j - 1];
                this.sparseTable[i, j] = MinByDepth(leftIndex, rightIndex);
            }
        }

        // Precompute mask for each block
        this.blockMasks = new int[this.blockCount];

        for (int i = 0, posInBlock = 0, blockIndex = 0; i < eulerTourLen; i++, posInBlock++)
        {
            if (posInBlock == this.blockSize)
            {
                posInBlock = 0;
                blockIndex++;
            }

            if (posInBlock > 0)
            {
                // If depth increases the bit is 1, otherwise 0 so skip the step
                if (MinByDepth(i - 1, i) == i - 1)
                {
                    this.blockMasks[blockIndex] += (1 << (posInBlock - 1));
                }
            }
        }

        // Precompute RMQ for each unique block
        int possibilities = 1 << (this.blockSize - 1);
        this.blocks = new int[possibilities][,];

        for (int blockIndex = 0; blockIndex < this.blockCount; blockIndex++)
        {
            int mask = this.blockMasks[blockIndex];

            if (this.blocks[mask] == null)
            {
                this.blocks[mask] = new int[this.blockSize, this.blockSize];
                int currentBlockStartsAt = blockIndex * this.blockSize;

                for (int left = 0; left < this.blockSize; left++)
                {
                    this.blocks[mask][left, left] = left;

                    for (int right = left + 1; right < this.blockSize; right++)
                    {
                        int minIndexSoFar = currentBlockStartsAt + this.blocks[mask][left, right - 1];
                        int rightIndex = currentBlockStartsAt + right;

                        this.blocks[mask][left, right] = MinByDepth(minIndexSoFar, rightIndex) - currentBlockStartsAt;
                    }
                }
            }
        }
    }

    int RMQInBlock(int blockIndex, int left, int right)
    {
        int mask = this.blockMasks[blockIndex];
        int blockStartsAt = blockIndex * this.blockSize;

        return blockStartsAt + this.blocks[mask][left, right];
    }

    int LCA(int first, int second)
    {
        int left = this.firstVisit[first];
        int right = this.firstVisit[second];

        if (left > right)
        {
            right = left;
            left = this.firstVisit[second];
        }

        int leftBlockIndex = left / blockSize;
        int rightBlockIndex = right / blockSize;

        if (leftBlockIndex == rightBlockIndex)
        {
            int rmqInBlock = RMQInBlock(leftBlockIndex, left % blockSize, right % blockSize);
            return this.eulerTour[rmqInBlock];
        }

        int ans1 = RMQInBlock(leftBlockIndex, left % this.blockSize, this.blockSize - 1);
        int ans2 = RMQInBlock(rightBlockIndex, 0, right % this.blockSize);
        int ans = MinByDepth(ans1, ans2);

        if (leftBlockIndex + 1 < rightBlockIndex)
        {
            int blockSeqLenPower = this.logTable[rightBlockIndex - leftBlockIndex - 1];
            int ans3 = this.sparseTable[leftBlockIndex + 1, blockSeqLenPower];
            int ans4 = this.sparseTable[rightBlockIndex - (1 << blockSeqLenPower), blockSeqLenPower];
            ans = MinByDepth(ans, MinByDepth(ans3, ans4));
        }

        return this.eulerTour[ans];
    }
}