1

So, I tried implementing my own ModExp implementation (I know I should use the existing one, but it's a task for school) and I'm getting a stack overflow when using large numbers (>= 2^1024)

modExp :: Integer -> Integer -> Integer -> Integer
modExp x y m
    | y == 0 = 1
    | y `mod` 2 == 0 = modExp ((x ^ 2) `mod` m) (y `div` 2) m
    | otherwise = x * modExp x (y - 1) m `mod` m

What is the problem? I heard that you have to use tail recursion to avoid stack overflows, but here I have two situations where the recursion must be called. How can I fix that?

1
  • 2
    The common trick to get rid of nontail recursion is to shove forcibly whatever you need to do with the result of the function to the function itself. In your case you need to multiply the result by x, so add an extra argument and pass your factor to it. Multiply inside, but before calling recursively. Commented Apr 9, 2017 at 16:57

1 Answer 1

4

Your function is singly recursive, so if we ensure it also is tail recursive we can maintain its memory size. The trick is not to produce more outstanding operations. The guards make the function strict in y, so we know that isn't some arbitrary tree of calculations waiting to happen. But x varies with call depth and isn't required until we produce the return value, and the return value has three cases: constant, tail recursive, and recursive with two operations (* and mod) outside the call. Those operations will go on some sort of call stack.

Solving this is probably a two step operation: First, make the function wholly tail recursive. Second, make it strict, such that it evaluates eagerly rather than building more and more complex expressions.

Tail recursive variant:

modExp :: Integer -> Integer -> Integer -> Integer
modExp x y m = modExp' x y m 1
    where
        modExp' x y m acc
            | y == 0 = acc
            | y `mod` 2 == 0 = modExp' ((x ^ 2) `mod` m) (y `div` 2) m acc
            | otherwise = modExp' x (y - 1) m (x * acc `mod` m)

It may be useful in some cases to enforce strictness in x and acc as well. That can be done with seq:

modExp :: Integer -> Integer -> Integer -> Integer
modExp x y m = modExp' x y m 1
    where
        modExp' x y m acc
            | y == 0 = acc
            | y `mod` 2 == 0 = let x' = (x ^ 2) `mod` m
                                   y' = y `div` 2
                               in x' `seq` modExp' x' y' m acc
            | otherwise = let acc' = x * acc `mod` m
                          in acc' `seq` modExp' x (y - 1) m acc'

A slightly simpler variant could apply seq to x or acc themselves, which in principle might leave one layer of calculation unfinished, but no more. Either way enforces constant space.

Unfortunately I didn't reproduce your stack overflow using ghc 7.10.3, so I'm not quite certain which lengths will prove necessary.

Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.