hanoi.hs

{-# LANGUAGE ScopedTypeVariables #-}

-- | Solve the r-peg Towers of Hanoi puzzle, using the Frame-Stewart algorithm.

module Hanoi where

import Prelude
  ( Bool(..), Int, Integer, Integral, Eq(..), Ord(..), Ordering(..)
  , error, fst, length, map, otherwise, snd, (+), (++), (-), (*), ($)
  , minimum, head, last
  )

import Control.Arrow ((&&&))
import Data.Array (Array,Ix,(!),listArray,range)
import Data.Bits (Bits(..),FiniteBits(..))
import Data.Function (on)
import Data.List (groupBy)

type Discs    = Int                      -- ^ Number of discs.
type Pegs     = Int                      -- ^ Number of pegs.
type K        = Int                      -- ^ Frame-Stewart parameter @K@.
type Memo res = Array (Discs, Pegs) res  -- ^ Memoised results, indexed by discs, pegs.
type Move peg = (peg, peg)               -- ^ Move a disc @(from,to)@.
type Moves    = Integer                  -- ^ Number of moves.

-- | Solve the @r@-peg Towers of Hanoi puzzle, using the Frame-Stewart algorithm.
--
hanoi :: Memo K     -- ^ Table of 'K'-parameters, at least @(n,r)@ big.
      -> Discs      -- ^ Number of discs, @n >= 1@.
      -> [peg]      -- ^ Peg lables, length @r >= 3@: (@start:finish:spare:spares@).
      -> [Move peg] -- ^ Solution, as list of moves: @[(from,to)..]@.
--
hanoi bestk n pegs = steps (n, length pegs) pegs where

  steps (1,_) (a:b:_) = [(a,b)]
  steps (n,r) (a:b:c:t)
      =  steps   (k,r)   (a:c:b:t)
      ++ steps (n-k,r-1) (a:b:t)
      ++ steps   (k,r)   (c:b:a:t)
    where
      k = bestk ! (n,r)

type Minimise = (K -> Moves) -> (K,K) -> (K,K)

-- | Generate memo tables for the @r@-peg Towers of Hanoi puzzle,
--   using the Frame-Stewart algorithm.
--
memos :: Minimise             -- ^ Use this minimise function
      -> Discs                -- ^ Maximum number of discs, @n >= 2@.
      -> Pegs                 -- ^ Maximum number of pegs, @r >= 3@.
      -> (Memo K, Memo Moves) -- ^ Optimal 'K'-parameters, moves.
--
memos minimise n r = (bestk, mover) where

  -- An optimal 'K'-parameter for the given discs and pegs, memoised.
  bestk :: Memo K
  bestk = memo ((2,3),(n,r)) f where
    f (n,3) = n - 1
    f (n,r) = fst $ minimise (movek (n,r)) (1,n-1)

  -- Number of moves to solve the @n@-disc, @r@-peg Towers of Hanoi,
  -- using the given 'K'-parameter at the top level of the recursion,
  -- and the minimal moves for lower levels of recursion.
  movek :: (Discs, Pegs) -> K -> Moves
  movek (n,r) k = 2 * moveb (k,r) + moveb (n-k,r-1)

  -- Minimal number of moves, given discs and pegs, base case.
  moveb :: (Discs, Pegs) -> Moves
  moveb (1,_) = 1
  moveb i = mover ! i

  -- Minimal number of moves, given discs and pegs, recursive case,
  -- memoised.
  mover :: Memo Moves
  mover = memo ((2,3),(n,r)) f where
    f i = movek i (bestk ! i)

-- | Memoise a function as a (bounded) array.
--
memo :: Ix i => (i,i) -> (i -> e) -> Array i e
memo bounds f = listArray bounds $ map f $ range bounds

-- | Find @k@ that minimises convex @f k@ on @[a..b]@, by linear search.
--
slowmin :: forall i r. (Integral i, Ord r)
        => (i -> r)            -- ^ Convex function @f@.
        -> (i,i)               -- ^ Interval @(a,b)@ (inclusive) to search.
        -> (i,i)               -- ^ Sub-interval @(h,k)@ which minimises @f@.
--
slowmin f (a,b) = snd $ minimum
  $ map (\g -> (fst $ head g, (snd $ head g, snd $ last g)))
  $ groupBy ((==) `on` fst) [ (f k, k) | k <- [a..b] ]

-- | Find @k@ that minimises convex @f k@ on @[a..b]@, by binary search.
--
minimise :: forall i r. (ILog2 i, Ord r)
         => (i -> r)            -- ^ Convex function @f@.
         -> (i,i)               -- ^ Interval @(a,b)@ (inclusive) to search.
         -> (i,i)               -- ^ Sub-interval @(h,k)@ which minimises @f@.
--
minimise f (a,b) = case compare u v of
  LT -> minl (a,b) u
  GT -> minr (a,b) v
  EQ -> flat (a,b) u
  where
    u = f a
    v = f b

    -- Begin the search for the minimum value.
    -- Completely handles functions that are increasing or decreasing only.
    -- Delegates concave and flat functions to @refine@.
    --
    init :: (i -> i -> i -> r -> (i,i))
         -> ((i,i) -> (i,i))
         -> Bool
         -> Bool
         -> (i,i)
         -> r
         -> (i,i)
    init rec done l r (a,b) v
      | g <= 1 = done (a,b)
      | otherwise = case compare u v of
          LT -> refine a b c c h u False False
          GT -> rec a b c v
          EQ -> refine a b c c h u l r
      where
        g = b - a
        h = bisect g
        c = a + h
        u = f c

    -- Specialised versions of init.
    minl = init (\a b c -> minl (a,c)) (fst &&& fst) True False
    minr = init (\a b c -> minr (c,b)) (snd &&& snd) False True
    flat = init (error "non-convex function") (fst &&& snd) True True

    -- Given at least three points, on each iteration, either:
    --   - narrow the search range @(a,b), if we find a point less than
    --     previously seen, or
    --   - narrow the grain size (@g@) by half, which means that we've
    --     sampled twice as many points, and found none less than
    --     what we've previously seen.
    --
    refine :: i      -- Absolute x-value of left-hand end (@a@).
           -> i      -- Absolute x-value of right-hand end (@b@).
           -> i      -- Absolute x-value of middle (@c@, with @0 < b - c <= g@).
           -> i      -- Absolute x-value of one grain from the left (@d@, with @d == a + g@)
           -> i      -- Granularity tested so far (@g@), a power of two.
           -> r      -- Minimum y-value seen so far (@v@).
           -> Bool   -- Is the left-hand end flat (@l@)?
           -> Bool   -- Is the right-hand end flat (@r@)?
           -> (i,i)  -- Range over which @f@ is minimised.
    refine a b c d g v l r
      | g == 1 = (if l then a else d, if r then b else c)
      | otherwise = case compare u v of
          GT -> scan e d e False
          LT -> refine a d e e h u False False
          EQ -> scan a e e l
      where
        h = unsafeShiftR g 1
        e = a + h
        u = f e

        -- Scan intervening points at half the current grain size.
        --
        scan :: i      -- Absolute x-value of left-hand end (@a@).
             -> i      -- Absolute x-value of one half-grain from the left (@e@, with @e == a + h@).
             -> i      -- Absolute x-value scanned so far (@s@).
             -> Bool   -- Is the left-hand end flat?
             -> (i,i)  -- Range over which @f@ is minimised.
        scan a e s l
          | t >= b    = refine a b c e h v l r
          | t > c     = case compare u v of
              LT -> refine c b t t h u False False
              GT -> refine a t (t-h) e h v l False
              EQ -> refine a b t e h v l r
          | u < v     = refine (t-h) (t+h) t t h u False False
          | otherwise = scan a e t l
          where
            t = s + g
            u = f t

-- | Integer logarithms.
--
class (Bits i, Integral i) => ILog2 i where

  -- | Given @i > 0@, find the smallest n such that @i <= 2^(n+1)@.
  --
  iLog2 :: i -> Int

instance ILog2 Int where iLog2 = finiteILog2

-- | Given @i > 0@, find the smallest @n@ such that @i <= 2^(n+1)@.
--
finiteILog2 :: forall i. (FiniteBits i, Integral i) => i -> Int
finiteILog2 i
  | i <= 1 = -1
  | otherwise = search 0 (finiteBitSize i)
  where
    -- invariant: @2^offset < i <= 2^(offset+width)@
    search :: Int -> Int -> Int
    search offset 1 = offset
    search offset width
      | test < i  = search middle half
      | otherwise = search offset half
      where
        half = unsafeShiftR width 1
        middle = offset + half
        test = bit middle

-- | Given @i > 1@, find the largest power of 2 less than @i@.
--
bisect :: ILog2 i => i -> i
bisect i = bit (iLog2 i)