Large Power Mod - Extreme Efficiency
TL;DR

Use strictness annotations at the right point to get rid of thunks:

largepowmod :: Integer -&gt; Integer -&gt; Integer -&gt; Integer
largepowmod num pow modbase = (operate num opList) `mod` modbase
  where
    opList = drop 1 $ binarizeE pow
    operate !k    []   = k
    operate !k (0:ops) = operate (k^2         `rem` modbase) ops
    operate !k (1:ops) = operate ((k^2 * num) `rem` modbase) ops
--          ^


Note that this is equivalent to using foldl'.

Proper benchmarks

First of all, let's get reliable benchmarks. We can use Criterion for this:

import Criterion.Main
import Text.Printf (printf)

main :: IO()
main = defaultMain [
  bgroup "largepowmod" $ 
    flip map ([1,3..5] ++ [9999]) $ \i -&gt;
      bench (printf "10^1000 + %d" i) $   
        let n = 10^1000 + i :: Integer
        in whnf (largepowmod 2 ((n - 1) `div` 2)) n
  ]


On my machine, I get about 16.4ms per call with -O2:

benchmarking largepowmod/10^1000 + 1
time                 16.14 ms   (16.07 ms .. 16.22 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.18 ms   (16.13 ms .. 16.23 ms)
std dev              124.6 μs   (99.83 μs .. 164.8 μs)

benchmarking largepowmod/10^1000 + 3
time                 16.44 ms   (16.26 ms .. 16.63 ms)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 16.35 ms   (16.29 ms .. 16.42 ms)
std dev              168.3 μs   (135.2 μs .. 223.2 μs)

benchmarking largepowmod/10^1000 + 5
time                 16.24 ms   (16.19 ms .. 16.30 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.35 ms   (16.30 ms .. 16.41 ms)
std dev              138.4 μs   (89.59 μs .. 182.9 μs)

benchmarking largepowmod/10^1000 + 9999
time                 16.45 ms   (16.34 ms .. 16.58 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.37 ms   (16.30 ms .. 16.43 ms)
std dev              165.2 μs   (129.3 μs .. 235.4 μs)


But let's not stop here. Let's have a look at the profile:

$ ghc --make Test.hs -O2 -prof -auto-all &amp;&amp; ./Test +RTS -s -p &amp;&amp; head -n20 Test.prof

        Tue Sep 27 13:32 2016 Time and Allocation Profiling Report  (Final)

           Test +RTS -s -p -RTS --output benchmark.html

        total time  =       19.65 secs   (19652 ticks @ 1000 us, 1 processor)
        total alloc = 8,760,057,056 bytes  (excludes profiling overheads)

COST CENTRE         MODULE                %time %alloc

largepowmod.operate Main                   91.9   83.5
getGCStats          Criterion.Measurement   5.1   10.6
getOverhead         Criterion.Monad         2.2    1.5
applyGCStats        Criterion.Measurement   0.1    1.5


                                                                                                        individual     inherited
COST CENTRE                                  MODULE                                   no.     entries  %time %alloc   %time %alloc

MAIN                                         MAIN                                     273           0    0.0    0.0   100.0  100.0
 main                                        Main                                     550           0    0.0    0.0   100.0   99.9
  main.\                                     Main                                     745           0    0.0    0.0    91.9   83.5
   largepowmod                               Main                                     746        1208    0.0    0.0    91.9   83.5
    largepowmod.operate                      Main                                     752     4011768   91.9   83.5    91.9   83.5
    largepowmod.opList                       Main                                     751        1208    0.0    0.0     0.0    0.0


So most of the time is used in operate. It's allocating most of the data. However, that's a little bit of a red herring, since binarizeE is inlined at this point, and poor operate is merely allocating the list that will be generated lazily.

Getting stricter

However, we're not strict enough. In operate we create a large chunk of intermediate data:

operate k    []   = k
operate k (0:ops) = operate (k^2         `rem` modbase) ops
operate k (1:ops) = operate ((k^2 * num) `rem` modbase) ops
--                           ^^^^^^^^^^^^^^^^^^^^^^^^^


Haskell is lazy, therefore (k^2 * num) `rem` modbase isn't evaluated until it's actually needed. Therefore we allocate \$\mathcal O(\log_2 p)\$ additional terms. This isn't necessary. We can provide a small fix to make sure that k has been fully evaluated:

largepowmod' :: Integer -&gt; Integer -&gt; Integer -&gt; Integer
largepowmod' num pow modbase = (operate num opList) `mod` modbase
  where
    opList = drop 1 $ binarizeE pow
    operate !k    []   = k
    operate !k (0:ops) = operate (k^2         `rem` modbase) ops
    operate !k (1:ops) = operate ((k^2 * num) `rem` modbase) ops


This needs -XBangPatterns. How does this fare? We change our main:

main :: IO()
main = defaultMain [
  bgroup "largepowmod" $ 
    flip map ([1,3..5] ++ [9999]) $ \i -&gt;
      bench (printf "10^1000 + %d" i) $   
        let n = 10^1000 + i :: Integer
        in nf (largepowmod 2 ((n - 1) `div` 2)) n
 , bgroup "largepowmod'" $ 
    flip map ([1,3..5] ++ [9999]) $ \i -&gt;
      bench (printf "10^1000 + %d" i) $   
        let n = 10^1000 + i :: Integer
        in nf (largepowmod' 2 ((n - 1) `div` 2)) n
  ]


And run it again (without -prof):

benchmarking largepowmod/10^1000 + 1
time                 16.14 ms   (16.07 ms .. 16.22 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.18 ms   (16.13 ms .. 16.23 ms)
std dev              124.6 μs   (99.83 μs .. 164.8 μs)

benchmarking largepowmod/10^1000 + 3
time                 16.44 ms   (16.26 ms .. 16.63 ms)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 16.35 ms   (16.29 ms .. 16.42 ms)
std dev              168.3 μs   (135.2 μs .. 223.2 μs)

benchmarking largepowmod/10^1000 + 5
time                 16.24 ms   (16.19 ms .. 16.30 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.35 ms   (16.30 ms .. 16.41 ms)
std dev              138.4 μs   (89.59 μs .. 182.9 μs)

benchmarking largepowmod/10^1000 + 9999
time                 16.45 ms   (16.34 ms .. 16.58 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 16.37 ms   (16.30 ms .. 16.43 ms)
std dev              165.2 μs   (129.3 μs .. 235.4 μs)

benchmarking largepowmod'/10^1000 + 1
time                 14.54 ms   (14.49 ms .. 14.60 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 14.51 ms   (14.46 ms .. 14.54 ms)
std dev              101.0 μs   (64.28 μs .. 177.1 μs)

benchmarking largepowmod'/10^1000 + 3
time                 14.53 ms   (14.49 ms .. 14.58 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 14.53 ms   (14.50 ms .. 14.56 ms)
std dev              78.32 μs   (62.67 μs .. 95.66 μs)

benchmarking largepowmod'/10^1000 + 5
time                 14.42 ms   (14.36 ms .. 14.48 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 14.56 ms   (14.51 ms .. 14.65 ms)
std dev              180.4 μs   (92.90 μs .. 298.2 μs)

benchmarking largepowmod'/10^1000 + 9999
time                 14.60 ms   (14.49 ms .. 14.85 ms)
                     0.998 R²   (0.993 R² .. 1.000 R²)
mean                 14.61 ms   (14.53 ms .. 14.87 ms)
std dev              314.4 μs   (69.97 μs .. 618.3 μs)


We went from 16.14ms to 14.77ms (0.89 of original time), or if we extrapolate this to your 110s to 98s, which is closer to Mathematica's time.