# Fastest modular exponentiation

The challenge is to write code in any open source, available on linux, language of your choice to perform modular exponentiation. The input will be two randomly chosen 2048 bit positive integers numbers x and y and a 2048 bit prime z. Here is a sample snippet of python to compare to.

def pow_mod(x, y, z):
n = 1
while y:
if y & 1:
n = n * x % z
y >>= 1
x = x * x % z
return n


The code should accept a file with each of the three numbers on separate lines and output the result to standard out. You may not use any libraries for modular exponentiation of course.

The code that runs the fastest on average over 100 runs on my computer wins.

For those who want really accurate timing, the best way may be for the timing code to actually be in the code you provide and for your code to simply repeat 100 times (without cheating :) ). This avoids any problems with start up overheads.

• A good solution will probably use some addition-chain or addition-subtraction chain wizardry. For example (I don't know if this particular approach is good, but...) cr.yp.to/ecdh/diffchain-20060219.pdf – boothby Oct 4 '13 at 16:53
• @Ray Ubuntu 64 bit 13.04 on amd fx8350 with 8GB RAM. – user9206 Oct 6 '13 at 7:56

## C++/asm

Runs in about 7ms on my laptop. Less than a factor of 2 away from GMP, which is heavily optimized for this task.

Uses Montgomery reduction to compute the modular reductions quickly. Uses a python script to generate x86 assembly to do the 2048 bit ops (multiply, add, etc.) in great swathes of straightline code.

Though it's simple, it should be pretty fast if built with GHC:

{-# OPTIONS_GHC -O2 #-}
{-# LANGUAGE BangPatterns #-}

module Main where

import Data.Bits

seqI :: Integer -> Integer
seqI x = seq x x

seqB :: Bool -> Bool
seqB x = seq x x

getInteger :: IO Integer
getInteger = fmap (seqI . read) getLine

modExp' :: Integer -> Integer -> Integer -> Integer -> Integer
modExp' r b p m = let nextR = seqI $mod (seqI$ r*b) m
nextB = seqI $mod (seqI$ b*b) m
nextP = seqI $shiftR p 1 incl = seqB$ testBit p 0
cont  = seqB $nextP /= 0 in if incl then if cont then modExp' nextR nextB nextP m else nextR else if cont then modExp' r nextB nextP m else r modExp :: Integer -> Integer -> Integer -> Integer {-# INLINE modExp #-} modExp !b !p !m = modExp' 1 b p m main = do b <- getInteger p <- getInteger m <- getInteger print$ modExp b p m


### Prime modulus optimization

If the modulus is guaranteed to be a prime, but the exponent might be greater than or equal to the modulus, the last line should be changed to:

          print \$ modExp b (mod p (m - 1)) m


### Handling multiple inputs

If you'd like it to handle multiple inputs per run, add import System.IO under the other import and the following two lines just under the print line:

          isDone <- isEOF
if isDone then return () else main


### Running

Once it's compiled, it can be run as modexp < file.

• I compiled it with ghc -O2. It takes about 0.02 seconds but maybe it makes more sense to have the timing in the code itself. – user9206 Oct 4 '13 at 18:22
• It might have issues with startup time. I've added a section on how to get it to run, taking 3 lines at a time, all the way until the file ends. You can use cat *.whatever | modexp if you have the test inputs separated into many files. – Olathe Oct 4 '13 at 19:13