1. Ordered Statistics Tree Part 1: product of finger trees

    One thing about finger trees is that we often want to argument a few different data. For example, I might want to know the number of elements in a finger tree, while still able to do something else.

    Certainly this is not hard easy, since we can just create a pair of two measurable types, add a sum type of integers. Just let the finger tree to measure both of them.

    In order to do that, we should write a function that takes the product of measures.

    productMeasure:: (Measured v a, Measured u b) => (a -> v)->(b -> u)->((a,b)->(v,u)).

    Try to do this by Control.Arrow. It should be a really short answer.

     
  2. The $m$th element of the union of $k$ non decreasing sequence.

    Problem

    Let $f_i:\mathbb{N}\to X$ for all $1\leq i\leq k$, where $X$ is totally ordered, and if $a\leq b$ then $f_i(a)\leq f_i(b)$.

    Write a function mthElementNonDecreasingSequences :: Ord a=>[(Integer->a)]->Integer->a that returns the $m$th element of the list that is created from combining all these sequence and sort them.

    Solution

    mthElementNonDecreasingSequences :: Ord a=>[Integer->a]->Integer->a 
    mthElementNonDecreasingSequences xs m = solve (zip xs (repeat (0,m))) (m+1)
      where solve :: Ord a=> [(Integer->a,(Integer,Integer))]->Integer->a
            solve [x] m = func x (left x+m-1)
            solve xs m  = solve bounds t
             where (t,op,op2)
                    | s<=m      = (m - diff v,minimum,removeDown)
                    | otherwise = (m, maximum, removeUp)
                   s = sum $ map diff xs
                   z = zip (map takeMid xs) [0..]
                   index = snd $ op z
                   bounds = filter (\x-> left x <= right x) $ take index xs ++ [ op2 v ] ++ drop (index+1) xs
                   v = xs!!index
            mid (_,(a,b))   = (a+b) `div` 2
            left (_,(a,_))  = a
            right (_,(_,b)) = b
            func (f,_)      = f
            takeMid a = func a (mid a)
            removeDown a = (func a, (mid a + 1,right a))
            removeUp   a = (func a, (left a, mid a - 1))
            diff a = mid a - left a + 1
    

    This solution takes $O(k^2 \log m)$ time.

     
  3. Alias method

    School is going to start, I don’t have as much time in updating these…

    Problem

    Implements the alias method.

    Solution

    import Data.Ratio
    import Data.Array
    import Data.List (partition)
    import System.Random
    
    type DiscreteDistribution a = Array Int (Rational, a, a)
    
    discreteDistribution :: [(a,Integer)] -> DiscreteDistribution a
    discreteDistribution xs = listArray (0,length xs-1) (uncurry buildTable $ partition ((<h).snd) xs')
      where xs' = map (\(x,y)->(x,y%1)) xs
            s = foldl (\y (_,x)-> x + y) (0%1) xs'
            n = fromIntegral (length xs)%1
            h = s / n
            buildTable [] ys = map (\(b,y)->(1%1,b,b)) ys
            buildTable ((a,x):xs) ((b,y):ys) = (x / h, a, b):sol
              where v = y-(h-x) 
                    sol
                     | v >= h = buildTable xs ((b,v):ys)
                     | v < h  = buildTable ((b,v):xs) ys
    
    randomElement :: (RandomGen g) => DiscreteDistribution a -> g -> (a, g)
    randomElement a g = (if u <= numerator v then x else y,g'')
      where (r,g')  = randomR (bounds a) g
            (v,x,y) = a!r
            (u,g'') = randomR (1, denominator v) g'
    
     
  4. Periodic Sequence

    A sequence is periodic if there exist a $p$, such that $a_i=a_{i+p}$ for all $i\geq 0$.

    Problem

    Implement the data PeriodicSequence a using Data.Sequence. Such that it is an Applicative, and it models function composition. eval f n :: PeriodicSequence a ->Int->a that evaluate the sequence at position n.

    Solution

    module PeriodicSequence(PeriodicSequence, eval, periodicSequence) where 
    import Prelude hiding (zipWith, foldr)
    import qualified Prelude as P
    import Data.Sequence hiding (take, length)
    import Data.Foldable
    import Data.Monoid
    import Control.Applicative
    data PeriodicSequence a = PS Int (Seq a) deriving (Eq, Show, Ord)
    
    instance Functor PeriodicSequence where
        fmap f (PS n t) = PS n (fmap f t)
    
    instance (Monoid a) => Monoid (PeriodicSequence a) where
        mempty  = PS 1 (singleton mempty)
        mappend x y = fmap (<>) x <*> y
    
    instance Applicative PeriodicSequence where
        pure f = PS 1 (singleton f)
        (<*>) (PS n xs) (PS m ys) = PS l (zipWith ($) xs' ys')
            where xs' = fromList $ take l $ cycle (toList xs)
                  ys' = fromList $ take l $ cycle (toList ys)
                  l = lcm n m
    
    eval :: PeriodicSequence a -> Int -> a
    eval (PS n xs) m = index xs (m `rem` n)
    
    periodicSequence :: [a]->PeriodicSequence a
    periodicSequence xs = PS (length xs) (fromList xs)
    
     
  5. Lattice Points on a Plane

    A point $(x,y)$ on the plane is a lattice point if both $x$ and $y$ are integers.

    Problem

    Implement latticePoints2D :: [(Integer,Integer)], a infinite list of lattice points, such that every lattice point is in the list.

    Solution

    import Control.Arrow
    import Control.Monad
    import Data.List
    import Data.Tuple
    latticePoints :: [(Integer,Integer)]                     --        90                  180                    270
    latticePoints = (0,0):concat (transpose $ map (`map` v) [id, swap . first negate, join (***) negate, swap . second negate])
      where v = takeLine 0 quadrant
            takeLine n xss = concat before ++ takeLine (n+1) (after ++ rest)
              where (prefix,rest) = splitAt n xss
                    (before, after) = unzip $ map (splitAt 1) prefix 
            quadrant = [[(x,y)|x<-[1..]]|y<-[0..]]
    
     
  6. First true in a monotonic predicate over the naturals

    Let False < True. A predicate $p$ is monotonic if for all $x\leq y$, $p(x)\leq p(y)$.

    Problem

    Write a function firstTrueMonotonicPredicate :: Integral a =>(a->Bool)->a that output the first position where a monotonic predicate over the natural numbers become true. The algorithm should evaluate the function $O(\log n)$ time, where $n$ is the output. The function does not terminate if all the values are False.

    Solution

    monotonicPredicate :: Integral a => (a -> Bool) -> a
    monotonicPredicate f = binSearch 0 ub
      where ub = 1 + expSearch 1
            expSearch i
             | f i       = i
             | otherwise = expSearch (2*i)
            -- search from [i,j)
            binSearch i j
             | j - i < 3 = if f i then i else i+1
             | f mid     = binSearch i (mid+1)
             | otherwise = binSearch (mid+1) j
             where mid = (i+j) `div` 2
    
     
  7. The minimum weight path of length $k$ in a DAG with Monge property

    Problem

    Consider we have a complete directed acyclic graph of $n$ vertices, such that the topological sort order is $v_1,\ldots,v_n$. There is a weight function $w(v_i,v_j)$, such that it has the Monge property, i.e. $w(v_i,v_j) + w(v_{i+1},v_{j + 1}) \leq w(v_i,v_{j+1}) + w(v_{i+1},v_{j})$ holds for all $1 < i + 1 < j \leq n$.

    The weight of a path is the sum of the weight of all edges in the path. We want to find a path of length $k$ from $v_1$ to $v_n$, such that the weight is minimized.

    The function would have the signature minWeightKLinkedPath :: (Integral a, Num b)=>a->a->(a->a->b)->[a]. The input is $n$, $k$, $w$ and the output is a minimum weight path of length $k$.

    Solve the problem in $O(kn)$ time.

    Hint

    There is an algorithm to solve it in $O(k(n+m))$ time on stackoverflow, which can be modified to run in $O(kn^2)$ time.

    Apply the technique in the talk Revisiting The Monge Property. The exercise on totally monotone matrices already implemented the subroutine!

    Extensions

    We can actually implement the most general form of the problem in Revisiting The Monge Property, and also reduce the memory usage from $O(kn)$ to $O(n)$. Looks like no one has done it with Haskell yet.

    Applications

    Consider the linear partition problem that got featured on hacker news. It can be reduced to this problem and produce a $O(kn)$ time solution. In the application $k$ is just $cn$ for some constant $c$, thus in theory there is a even better (in theory) solution that takes $n2^{O(\sqrt{\log k \log \log n})}$ time discovered by Schieber.

    Solution

    minWeightkEdgePath :: (Ord a, Num a) => Int -> Int -> (Int -> Int -> a) -> [Int]
    minWeightkEdgePath k n w = backtrack (k-1) n
      where 
        backtrack (-1) _ = []
        backtrack x y = y':backtrack (x-1) y'
          where y' = snd (d!(x,y))
        d = array ((0,0),(k-1,n)) [ ((i,j), f i j) | i<-[0..k-1],j<-[0..n]]
        f 0 i = (w 0 i,0)
        -- uncomment the following line to use the slower algorithm
        -- f k i = minimum [(m k j i, j)|j<-[0..n]]
        f k i = (m k j i,j)
          where j = p!(k,i)
        p = array ((1,0),(k,n)) [((z,i),t) |z<-[1..k],(i,t)<-zip [0..n] (columnMinima (m z) (n+1) (n+1))]
        m k j i = fst (d!(k-1,j)) + w j i
    
     
  8. Test if a integer is a perfect power

    Problem

    A integer $n$ is a perfect power if $n=m^k$ for some integer $m$ and $k>1$.

    Write a program to check if a number is a perfect power.

    Solution

    import Math.NumberTheory.Primes.Heap
    
    -- Check if n = m^k. We do a normal search on k.
    lg :: Integral a => a->a
    lg n = monotonicPredicate (\k-> 2^k>n)
    
    isPerfectPower :: Integral a => a -> Bool
    isPerfectPower n = any (isExponent n) $ take (fromIntegral $ lg n) primes
    
    isExponent :: Integral a => a -> a -> Bool
    isExponent n m = cand ^ m == n
      where cand = monotonicPredicate (\k -> k^m >= n) 
    
    monotonicPredicate :: Integral a => (a -> Bool) -> a
    monotonicPredicate f = binSearch 0 ub
      where ub = 1 + expSearch 1
            expSearch i
             | f i       = i
             | otherwise = expSearch (2*i)
            -- search from [i,j)
            binSearch i j
             | j - i < 3 = if f i then i else i+1
             | f mid     = binSearch i (mid+1)
             | otherwise = binSearch (mid+1) j
             where mid = (i+j) `div` 2
    
     
  9. Sum over a consecutive subsequence

    Problem

    Solve the Interval Product problem by create a set of functions that solve the following more abstract problem using Data.FingerTree.

    Let’s consider a sequence of elements in a monoid of $n$ elements, we are interested in the ability to modify an element in position i, insert/delete at position i, and calculate the sum of all elements between position i and j. All these operations should take $O(\log n)$ amortized time under the assumption that all monoid operation takes $O(1)$ time.

    Motivation: This exercise shows that finger trees can replace binary indexed (Fenwick) trees. The abstract problem can be used to solve problems like Potentiometers and Frequent Values.

    Solution

    MonoidSequence.hs

    {-# LANGUAGE FlexibleInstances #-}
    {-# LANGUAGE MultiParamTypeClasses #-}
    {-# LANGUAGE FlexibleContexts #-}
    module MonoidSequence(MonoidSequence, splitPosition, update, insert, delete, infixSum, fromList) where
    import Data.FingerTree hiding (fromList)
    import qualified Data.FingerTree as F
    import Data.Monoid
    type MonoidSequence a = FingerTree (Sum Int, a) (S a)
    
    newtype S a = S a
    instance Monoid a => Measured (Sum Int, a) (S a) where
      measure (S x) = (Sum 1, x)
    
    splitPosition :: Measured (Sum Int, a) (S a) => Int -> MonoidSequence a -> (MonoidSequence a, MonoidSequence a)
    splitPosition i = split (\(Sum x,_)-> x > i)
    
    update :: Measured (Sum Int, a) (S a) => Int -> a -> MonoidSequence a -> MonoidSequence a
    update i entry tree = insert i entry (delete i tree)
    insert :: Measured (Sum Int, a) (S a) => Int -> a -> MonoidSequence a -> MonoidSequence a
    insert i entry tree = x >< singleton (S entry) >< y
      where (x,y) = splitPosition i tree
    delete :: Measured (Sum Int, a) (S a) => Int -> MonoidSequence a -> MonoidSequence a
    delete i tree = x >< y
      where (x,z) = splitPosition i tree
            (_,y) = splitPosition 1 z
    
    -- sum from position i to position j not inclusive
    infixSum :: Measured (Sum Int, a) (S a) => Int -> Int -> MonoidSequence a -> a
    infixSum i j tree = snd $ measure y
      where (_,z) = splitPosition i tree
            (y,_) = splitPosition (j-i) z
    
    fromList :: Monoid a => [a]-> MonoidSequence a
    fromList xs = F.fromList $ map S xs
    

    IntervalProduct.hs

    -- Interval Product 
    -- ICPC Latin American Regional 2012 Problem I
    -- ACM-ICPC Live Archive Problem 6139
    
    import MonoidSequence
    import Data.Monoid
    import Data.List.Split
    
    data Sign = N | Z | P deriving (Eq)
    
    instance Show Sign where
      show N = "-"
      show P = "+"
      show Z = "0"
    
    instance Monoid Sign where
      mempty = P
      mappend Z _ = Z
      mappend _ Z = Z
      mappend P x = x
      mappend x P = x
      mappend N N = P
    
    toSign :: Int -> Sign
    toSign x
     | x == 0 = Z
     | x > 0  = P
     | x < 0  = N
    
    main :: IO ()
    main = interact process
    
    process :: String -> String
    process input = sol l
      where 
            l = words input
            sol (n':k':xs) = solve (fromList values) operations ++ sol (drop (n+3*k) xs)
              where 
                (n,k) = (read n', read k')
                values = map (toSign.read) (take n xs)
                operations = map (\[a,b,c]-> (a, read b, read c)) $ chunksOf 3 $ take (3*k) (drop n xs)
            sol _ = ""
    
    solve :: MonoidSequence Sign -> [(String, Int, Int)] -> String
    solve values ((a,b,c):operations)
     | a == "P" = show (infixSum (b-1) c values) ++ solve values operations
     | a == "C" = solve (update (b-1) (toSign c) values) operations
    solve _ [] = "\n"
    
     
  10. Column Minima in a totally monotone matrix

    Let type Matrix a = Int->Int->a to denote a matrix. $A=\{a_{i,j}\}$ for $1\leq i\leq n, 1\leq j\leq m$ is a totally monotone matrix is if for all $i < i’$ and $j < j’$, we have $a_{i,j} > a_{i’,j} \Rightarrow a_{i,j’} > a_{i’,j’}$. We want an algorithm that uses $O(n+m)$ evaluation of the entries of the matrix to find out the row position of each column’s minima.

    The SMAWK algorithm solves this problem. Here is an video which features the algorithm around 37 minutes into the video. David Eppstein have a python implementation of the algorithm.

    Implement columnMinima :: Ord a => Matrix a->Int->Int->[Int], which takes a totally monotone matrix with it’s width and height, it find the list of column minima’s row position.

    Example

    SMAWK.hs

    module SMAWK (columnMinima) where 
    -- matrix is totally monotone.
    matrix i j = (a!!i)!!j
      where a = [[8,5,5,13,53,68],
                  [m,5,1,5,37,50],
                  [m,m,5,1,17,26],
                  [m,m,m,13,13,18],
                  [m,m,m,m,53,50],
                  [m,m,m,m,m,68]]
            m = 10000
    
    ghci> columnMinima matrix 6 6
    [0,0,1,2,3,3]
    

    Solution

    My solution is just a Haskell translation of Eppstein’s code. It requires evenOddSplit and horizontalPath.

    module SMAWK (columnMinima) where import qualified Data.IntMap as M import Data.Array

    type Coord = M.IntMap Int
    type Matrix a = Int->Int->a
    
    columnMinima :: Ord a => Matrix a->Int->Int->[Int]
    columnMinima m w h = [(M.!) pos x|x<-[0..w-1]]
      where pos = smawk [0..h-1] [0..w-1] m
    
    -- Rows, Cols, Matrix, Maximum
    smawk :: Ord a => [Int]->[Int]->Matrix a ->Coord
    smawk rows cols m
     | null cols  = M.empty
     | otherwise  = interpolate row2
     where row2 = reverse $ reduce rows []
           minima = smawk row2 (snd $ evenOddSplit cols) m
           cols'  = listArray (0,n-1) cols
           n = length cols
           -- the reduce part
           reduce :: [Int]->[Int]->[Int]
           reduce xs ys -- ys is stack
            | null xs = ys
            | null ys = cont
            | m y ci > m x ci = reduce xs (tail ys)
            | t /= n = cont
            | otherwise = reduce (tail xs) ys
            where x = head xs
                  y = head ys
                  ci = cols'!(t-1)
                  cont = reduce (tail xs) (x:ys)
                  t = length ys
           interpolate :: [Int]->Coord
           interpolate rows = M.fromList (map f broken) `M.union` minima
             where
               broken = zip (fst $ evenOddSplit cols)
                             (horizontalPath rows ([head rows]++[M.findWithDefault (last rows) (cols'!(c+1)) minima | c<-[0,2..n-2]]++[last rows]))
               f (a,b) = (a, snd $ minimum $ map g b)
                 where g i = (m i a, i)