Competitive programming in Haskell: permutations

In my previous post I challenged you to solve The Power of Substitution. This problem presents us with a substitution cipher, and asks how many times we would have to iterate the encoding process in order to turn a given message into a given encryption.

A non-solution

Several commenters mentioned that they tried some sort of straightforward brute force approach: just iterate the encoding process and count how many iterations are needed to reach the specified encryption. This certainly works for the provided sample inputs. However, I hinted that this approach likely wouldn’t fit within Kattis’s time limit for the other, secret, test cases.

How did I know this solution would exceed the time limit? It is not just an issue of using efficient data structures! Even if someone told me they coded this straightforward solution in C++, I would have the same reaction. The problem specifies that the answer should be at most 10^9—and, as we will see, it is not hard to come up with test cases where the answer is indeed that large—and the simple fact is that counting to 10^9 and applying the encryption at each step (which requires iterating through a list of length 200) is definitely going to take more than 1 second. A useful rule of thumb that I learned from Competitive Programming 3 is 10^8 operations per second. (Of course your processor can carry out much more than 10^8 instructions per second, but in practice this rule of thumb seems to work remarkably well for predicting run times up to the right order of magnitude.)

Cycle decomposition of permutations

The given encoding, of course, is a permutation (the problem specifies only that it will be one-to-one, but any one-to-one endofunction on a finite set must in fact be a permtuation). Let’s call it p. If we start with an arbitrary m and repeatedly apply p—that is, m, p(m), p^2(m), p^3(m), \dots—what happens? Of course, because of the pigeonhole principle, the sequence must eventually repeat. But actually, something stronger is true: because p is a permutation, the first repeated value must be m itself. For suppose that p^i(m) = p^j(m) was the first repeated value in the sequence. But then since p is one-to-one, it must be the case that p^{i-1}(m) = p^{j-1}(m) as well, which means p^{i-2}(m) = p^{j-2}(m), and so on up to m = p^0(m) = p^{j-i}(m).

So in fact if we start at an arbitrary m and iterate p, we will find a cycle that includes m (including the possibility of a trivial length-1 cycle if p(m) = m). If there are other elements not included in m’s cycle, we can pick any one of them and repeat the process to find another cycle (which can’t possibly overlap at all with m’s cycle—do you see why?). In general, any permutation can be decomposed in this way into a collection of disjoint cycles.

CRT to the rescue

This idea of cycle decomposition is the key to unlocking the problem. Think about what happens to a particular letter m_i in the message, which we eventually want to become c_i. This will happen after applying the permutation some small number of times j, such that p^j(m_i) = c_i. (In general, of course, it would be quite possible that m_i and c_i are not in the same cycle at all, and so m_i will never turn into c_i no matter how many times we apply p; but the problem statement guarantees that this will not be the case.)

The problem, of course, is that all the other letters may not be encrypted properly after only j encryption steps, in which case we need to keep going until all the cycles line up. Suppose m_i and c_i are in a cycle of length n. That means applying the encryption n times to c_i will result in c_i again (and not before). Thus, we will get c_i after j steps and then every n steps thereafter (j, j+n, j+2n, \dots); in other words, the number of encryption steps must be equivalent to j \pmod n.

Every position in the message yields a similar modular equivalence, giving us a system of up to 200 simultaneous modular equivalences which we can solve using the generalized Chinese Remainder Theorem. Incidentally, this is why the solution can be so large—if we have cycles of sizes n_1, n_2, \dots, n_q, then it could take up to \mathrm{lcm}(n_1, n_2, \dots, n_q) iterations for them to all line up. In the special case that all the cycle sizes are relatively prime, this is just their product. So for example we could have cycles of sizes 2, 3, 5, 7, 11, 13, 17, 19, 23, which add up to exactly 100; the product of these is the primorial 23\# = 223\,092\,870.

My solution

Here is my solution. We’re going to use an unboxed array to represent the permutation, and the implementation of GCRT from a previous post.

{-# LANGUAGE RecordWildCards #-}

import           Control.Arrow
import           Data.Array.Unboxed

import           NumberTheory
import           Scanner

Here’s main, along with a data type to represent a single test case and a Scanner for reading one in. I like using record syntax to help me remember which field is which, combined with the RecordWildCards extension to unpack the data structure and get the field names as local variables.

main = interact $
  runScanner (numberOf tc) >>> map (solve >>> show) >>> unlines

data TC = TC { message :: [Int], crypt :: [Int], subst :: [Int] }
  deriving Show

tc = do
  l <- int
  TC <$> l `times` int <*> l `times` int <*> 100 `times` int

We’re going to represent a permutation as an unboxed array, which gives us nice constant-time lookup. I often use Data.Array or Data.Array.Unboxed to represent read-only information (which comes up a lot more than you might think!), giving all the benefits of fast array access with none of the complications of mutability.

type Perm = UArray Int Int

The distance between i and j is simply the number of times we have to apply p to get from i to j. Of course this definition would hang if i and j are not part of the same cycle, but we know they will be. Given dist, we can also find the length of a cycle containing i as one more than the distance from p(i) to i. (We can’t just ask for the distance from i to itself since that would return 0.)

dist :: Perm -> Int -> Int -> Int
dist p i j = length $ takeWhile (/= j) (iterate (p!) i)

cycleLen :: Perm -> Int -> Int
cycleLen p i = succ $ dist p (p!i) i

Finally, we can put these pieces together: create an array for the permutation, zip together the message and desired encryption, generating a modular equivalence for each, and solve the resulting system using gcrt.

solve :: TC -> Int
solve (TC{..}) = pick . gcrt $ zipWith modEqn message crypt
    p :: UArray Int Int
    p = listArray (1,100) subst

    modEqn m c = (fromIntegral (dist p m c), fromIntegral (cycleLen p m))
    pick (Just (z,k)) = fromIntegral (z `mod` k)

Incidentally, this code inspired me to create a Util.hs in my comprog-hs repository containing (for now) fi as an alias for fromIntegral, and both to apply a function to both elements of a tuple (sadly lens is not available in the Kattis environment). Then we can just write modEqn m c = both fi (dist p m c, cycleLen p m).

Solving bigger instances

The above solution works because the alphabet is quite small (only 100). However, it’s actually quite wasteful. For example, suppose that the given message consists of 200 copies of the number 1; then we will recompute the length of 1’s cycle 200 times. It’s easy to imagine a variant of this problem where both the message length and the alphabet size could be much larger. Then my solution above would be too slow. For example, suppose the permutation consists of one giant cycle of length 10^5, and the message also has length 10^5. We would traverse the entire cycle for every single character in the message, for a total of about 10^{10} operations—much too slow. This post has gotten long enough, but in another post I will show an alternative solution which I believe would work quickly enough even for such large inputs (assuming that the input was restricted such that the answer was still of a reasonable size!). The idea is to precompute the cycle decomposition of the permutation (in time proportional to the size of the alphabet), storing the information in such a way that for each pair of letters in the message and desired encryption, we can find the distance between them and the length of their cycle in constant time.

Next up: geometry

Next, I’d like to spend a few posts on the topic of geometry. I find that geometry problems work particularly well in Haskell (I don’t think I’ve solved a single geometry problem in Java). Let’s kick things off with a problem on the easier side:

Have fun!

About Brent

Assistant Professor of Computer Science at Hendrix College. Functional programmer, mathematician, teacher, pianist, follower of Jesus.
This entry was posted in competitive programming, haskell and tagged , , , , , , , . Bookmark the permalink.

5 Responses to Competitive programming in Haskell: permutations

  1. Justin says:

    Really appreciate this write up. As usual, I’m amazed to see how clean a solution can be.

  2. shaurya gupta says:

    Been waiting for this!
    Also realized I did an unnecessary step in the last problem where I tried to reach to the first point inside the loop before applying the CRT logic.

    Here is my solution to today’s problem –

    The main part was parsing (I used your scanner, I will start working on a CP utils module soon) and using printf to control the precision while printing the real numbers.

    type Coord = (Double, Double)
    type Angle = Double
    parseInput :: String -> [[(Angle, Double)]]
    parseInput = runScanner $ numberOf $ numberOf $ (,) . (\d -> d/180*pi)  double  double
    solve :: [(Angle, Double)] -> String
    solve xs = let (x, y) = snd . foldl move (pi/2, (0, 0)) $ xs
                 in printf "%.4f %.4f" x y
    move :: (Angle, Coord) -> (Angle, Double) -> (Angle, Coord)
    move (theta, (x, y)) (dTheta, dist) = let theta' = dTheta + theta
                                              x' = x + dist * cos theta'
                                              y' = y + dist * sin theta'
                                            in (theta', (x', y'))
  3. Globules says:

    Here is the interesting part of my solution. I use complex numbers to do the vector additions.

    import Data.Complex (Complex(..), mkPolar)
    import Data.List (foldl', scanl')
    -- Treat each move as a complex number in polar coordinates.  We begin with a 90
    -- degree rotation because our starting position is along the Y axis.  The final
    -- position is given by the sum of the numbers.
    solve :: [(Double, Double)] -> [Double]
    solve mvs = let (x :+ y) = sum' $ map fromPolar $ normAngles ((90, 0) : mvs)
                    sum' = foldl' (+) 0
                in [x, y]
    -- Normalize the angles so that they're all relative to the X axis, rather than
    -- relative to the previous angle.
    normAngles :: [(Double, Double)] -> [(Double, Double)]
    normAngles = scanl' (\(a, _) (a', l') -> (a + a', l')) (0, 0)
    -- A complex number from a polar representation of (degrees, length).
    fromPolar :: (Double, Double) -> Complex Double
    fromPolar (deg, len) = mkPolar len (pi/180 * deg)
  4. Pingback: Competitive programming in Haskell: cycle decomposition with mutable arrays | blog :: Brent -> [String]

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.