Competitive programming in Haskell: folding folds

Now that the semester is over—and I will be on sabbatical in the fall!—you can expect a lot more competitive programming in Haskell posts. In a previous post, I challenged you to solve Origami. j0sejuan took me up on the challenge, as did Aaron Allen and Ryan Yates; if you still want to try it, go do it before reading on!

In the problem, we start with a square sheet of paper and are given a series of folds to perform in sequence; each fold is specified as a line, and we fold whatever is on one side of the line across onto the other side. Given some query points, we have to compute how thick the resulting origami design is at each point.

Lines

The first order of business is some computational geometry relating to lines in 2D (this code can all be found in Geom.hs. Here I am following Victor Lecomte’s excellent Handbook of geometry for competitive programmers, which I think I’ve mentioned before. I’ll try to give a bit of explanation, but if you want full explanations and proofs you should consult that document.

The equation of a line \(ax + by = c\) can be thought of as the set of all points \((x,y)\) whose dot product with the vector \((a,b)\) is a constant \(c\). This will in fact be a line perpendicular to the vector \((a,b)\), where \(c\) determines the distance of the line from the origin. Alternatively, we can think of the vector \((b,-a)\), which is perpendicular to \((a,b)\) and thus parallel to the line; the line now consists of all points \((x,y)\) whose 2D cross product with \((b,-a)\) is the constant \(c\) (since \((b,-a) \times (x,y) = by - (-a)x = ax + by\); note that the order matters). Either representation would work, but I will follow Lecomte in choosing the second: we represent a line by a vector giving its direction, and a scalar offset.

data L2 s = L2 { getDirection :: !(V2 s), getOffset :: !s }
type L2D = L2 Double

There are a few ways to construct a line: from an equation \(ax + by = c\), or from two points which lie on the line. The first is easy, given the above discussion. For the second, given points \(p\) and \(q\), we can easily construct the direction of the line as \(v = q - p\). Then to get the constant \(c\), we simply use the fact that \(c\) is the cross product of the direction vector with any point on the line, say, \(p\) (of course \(q\) would also work).

lineFromEquation :: Num s => s -> s -> s -> L2 s
lineFromEquation a b c = L2 (V2 b (-a)) c

lineFromPoints :: Num s => P2 s -> P2 s -> L2 s
lineFromPoints p q = L2 v (v `cross` p)
  where
    v = q ^-^ p

Now we can write some functions to decide where a given point lies with respect to a line. First, the side function computes \(ax + by - c = (b,-a) \times (x,y) - c\) for any point \(p = (x,y)\).

side :: Num s => L2 s -> P2 s -> s
side (L2 v c) p = cross v p - c

Of course, for points that lie on the line, this quantity will be zero. We can also classify points \(p\) as lying to the left or right of the line (looking in the direction of \(v\)) depending on whether side l p is positive or negative, respectively.

onLine :: (Num s, Eq s) => L2 s -> P2 s -> Bool
onLine l p = side l p == 0

leftOf :: (Num s, Ord s) => L2 s -> P2 s -> Bool
leftOf l p = side l p > 0

rightOf :: (Num s, Ord s) => L2 s -> P2 s -> Bool
rightOf l p = side l p < 0

The last piece we will need to solve the problem is a way to reflect a point across a line. toProjection l p computes the vector perpendicular to \(l\) which points from \(p\) to \(l\), and reflectAcross works by adding toProjection l p to p twice. I won’t derive the definition of toProjection, but the basic idea is to start with a vector perpendicular to the direction of the line (perp v) and scale it by a factor related to side l p. (Intuitively, it makes sense that \(ax + by - c\) tells us something about the distance from \((x,y)\) to the line; the farther away \((x,y)\) is from the line, the farther \(ax + by\) is from \(c\).) See Lecomte for the full details.

toProjection :: Fractional s => L2 s -> P2 s -> V2 s
toProjection l@(L2 v _) p = (-side l p / normSq v) *^ perp v

project :: Fractional s => L2 s -> P2 s -> P2 s
project l p = p ^+^ toProjection l p

reflectAcross :: Fractional s => L2 s -> P2 s -> P2 s
reflectAcross l p = p ^+^ (2 *^ toProjection l p)

Folding origami

Finally we can solve the problem! First, some imports and input parsing.

{-# LANGUAGE RecordWildCards #-}

import           Control.Arrow
import qualified Data.ByteString.Lazy.Char8 as C

import           Geom
import           ScannerBS

main = C.interact $
  runScanner tc >>> solve >>> map (show >>> C.pack) >>> C.unlines

data TC = TC { steps :: [L2D], holes :: [P2D] }

tc = TC <$> numberOf (lineFromPoints <$> p2 double <*> p2 double) <*> numberOf (p2 double)

solve :: TC -> [Int]
solve TC{..} = map countLayers holes
  where

For countLayers, the idea is to work backwards from a given query point to find all its preimages, that is, the points that will eventually map to that point under the folds. Then we can just count how many of those points lie (strictly) inside the original square of paper.

    inSquare (V2 x y) = 0 < x && x < 1000 && 0 < y && y < 1000

For a given point and fold, there are two possibilities, depending on which side of the fold line the point falls on. If the point falls on the fold or to the right of it, then it has no preimages (we always fold from right to left, so after the fold, there will be no paper on the right side of the line, and the problem specifies that points exactly on a folded edge do not count). Hence we can just discard such a point. On the other hand, if the point lies on the left side of the line, then the point has two preimages: the point itself, and its reflection across the fold line.

    preimage :: L2D -> P2D -> [P2D]
    preimage l p
      | leftOf l p = [p, reflectAcross l p]
      | otherwise  = []

So we keep a set of points, starting with the singleton query point, and for each fold (in order from last to first) we find the preimage of every point in the set under the fold. We actually use lists of points instead of sets, because (1) we won’t ever get any collisions (actually, the more I think about this, the less sure I am!) and (2) it lets us use the actual list monad instead of making some ad-hoc Set monad operations.

    countLayers :: P2D -> Int
    countLayers q = length . filter inSquare $ foldr (\l -> (>>= preimage l)) [q] steps

It is very satisfying to use a fold to process a list of folds!

Next time: Please, Go First

For next time, I invite you to solve Please, Go First.