Programming Praxis – Modular Arithmetic

Yesterday’s Programming Praxis problem is about making a convenient way to do modular arithmetic. While the most elegant way to do this in Haskell is probably via a Num instance, my attempts at achieving this have failed because either the typechecker cannot infer everything anymore, or because you cannot use the same syntax for modular square roots, congruency and the rest of the operations. Maybe a better Haskell hacker than me can come up with a working solution.

In the meantime, we’re going to use a solution that’s consistent and plays nice with the type checker, even if it’s a tad less elegant.

First our imports:

import Data.Bits
import Data.List

We’ll be needing these two functions for divisions.

euclid :: Integral a => a -> a -> a
euclid x y = f 1 0 x 0 1 y where
f a b g u v w | w == 0    = mod a y
| otherwise = f u v w (a-q*u) (b-q*v) (g-q*w)
where q = div g w

inv :: Integral a => a -> a -> a
inv x m | gcd x m == 1 = euclid x m
| otherwise    = error "divisor must be coprime to modulus"

All the basic operations are pretty simple:

(<=>) :: Integral a => a -> a -> a -> Bool
a <=> b = \m -> mod a m == mod b m

(<+>), (<->), (<*>), (</>) :: Integral a => a -> a -> a -> a
a <+> b = \m -> mod (a + mod b m) m
a <-> b = a <+> (-b)
a <*> b = \m -> mod (a * mod b m) m
a </> b = \m -> (a <*> inv b m) m

For exponentiation we get to recycle the expm function, which was used in previous exercises.

expm :: Integer -> Integer -> Integer -> Integer
expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 .
filter (flip testBit 0 . snd) .
zip (iterate (flip mod m . (^ 2)) b) \$
takeWhile (> 0) \$ iterate (`shiftR` 1) e

And a little syntactic shortcut for expm:

(<^>) :: Integer -> Integer -> Integer -> Integer
(<^>) = expm

This modular square root implementation is undoubtedly a lot slower than the provided one. However, implementing that one would be boring (just translating Scheme to Haskell, with not much opportunity for size reduction) and would double or triple the size of the code. Therefore we use a more naive algorithm, that works well enough when the modulus is small. Copying the provided algorithm to Haskell is left as an exercise for the reader.

sqrtm :: Integral a => a -> a -> [a]
sqrtm x m = case [n | n <- [1..m], (n <*> n) m == mod x m] of
(_:_:_:_) -> []
s         -> s

The modulo function is just syntactic sugar for reverse function application.

modulo :: a -> (a -> b) -> b
modulo = flip (\$)

This results in test code that looks as follows:

main :: IO ()
main = do mapM_ (modulo 12) [
print . (17 <=> 5),
print . (8 <+> 9),
print . (4 <-> 9),
print . (3 <*> 7),
print . (9 </> 7)]
mapM_ (modulo 13) [
print . (6 <^> 2),
print . (7 <^> 2),
print . sqrtm 10]

Not terrible, but I believe a better solution exists. If anyone has one, please let me know.

5 Responses to “Programming Praxis – Modular Arithmetic”

1. Lii Says:

I been looking for a good solution to this problem, but there seems to be surprisingly few around.

I dont immediately see how your solution would look in a more complicated example. Do you have to write out the same modulo number more than once?

For example in this hideous function to add two points on an elliptic curve:

type Point = (Integer, Integer)
inf = (-1, -1)

eadd :: Point -> Point -> Eps -> Point
eadd p1@(x1, y1) p2@(x2, y2) eps@Eps{pVal = mo, bVal = b}
| p1 == inf = p2
| p2 == inf = p1
| (x1 == x2 && y1 /= y2) || (y1 == 0 && y2 == 0) = inf
| True =
let
m = if x1 == x2 then
mdiv (madd (mmult 3 (mexp x1 2 mo) mo) b mo) (mmult 2 y1 mo) mo
else
mdiv (msub y2 y1 mo) (msub x2 x1 mo) mo
x = msub (msub (mexp m 2 mo) x1 mo) x2 mo
y = msub (mmult m (msub x1 x mo) mo) y1 mo
in
(x, y)

• Lii Says:

It had some indentation when I pasted it in. Stupid WordPress.

type Point = (Integer, Integer)
inf = (-1, -1)

eadd :: Point -> Point -> Eps -> Point
eadd p1@(x1, y1) p2@(x2, y2) eps@Eps{pVal = mo, bVal = b}
| p1 == inf = p2
| p2 == inf = p1
| (x1 == x2 && y1 /= y2) || (y1 == 0 && y2 == 0) = inf
| True =
let
m = if x1 == x2 then
mdiv (madd (mmult 3 (mexp x1 2 mo) mo) b mo) (mmult 2 y1 mo) mo
else
mdiv (msub y2 y1 mo) (msub x2 x1 mo) mo
x = msub (msub (mexp m 2 mo) x1 mo) x2 mo
y = msub (mmult m (msub x1 x mo) mo) y1 mo
in
(x, y)

• Lii Says:

2. Lii Says:

After more investigation it seems like it is not possible to compose your mod operations to a compound expression:

n = modulo 13 (5 6 7)

The reason is that you functions take integers as arguments, but return functions, so the result of one operation cant be the argument to the next.

Hm, what would be a good solution to that…?

3. Lii Says:

Here is one approach where the operation take and return functions. The disadvantage is that you cant use ints directly in the expressions, you have to wrap them in the ugly i functions.

tfunc = ((i 5) *% (i 6) +% (i 7)) 13

(=%) :: TModNum -> TModNum -> (Integer -> Bool)
a =% b = \m -> (a m) == (b m)

type TModNum = (Integer -> Integer)

(+%), (-%), (*%), (/%) :: TModNum -> TModNum -> TModNum

a +% b = \m -> mod ((a m) + (b m)) m
a -% b = \m -> mod ((a m) - (b m)) m
a *% b = \m -> mod ((a m) * (b m)) m
a /% b = \m -> mod ((a m) * (inv (b m) m)) m

i :: Integer -> (Integer -> Integer)
i n = \m -> mod n m