bodigrim / mod Goto Github PK

Modular arithmetic, promoting moduli to the type level

Home Page: http://hackage.haskell.org/package/mod

License: MIT License

Haskell 99.52% C 0.48%

modular-inversions moduli exponentiation modular-arithmetic modular-exponentiation inversion finite-numbers number-theory

mod's Introduction

mod

Modular arithmetic, promoting moduli to the type level, with an emphasis on performance. Originally a part of the arithmoi package.

> :set -XDataKinds
> 4 + 5 :: Mod 7
2
> 4 - 5 :: Mod 7
6
> 4 * 5 :: Mod 7
6
> 4 / 5 :: Mod 7
5
> 4 ^ 5 :: Mod 7
2

Competitors

There are other Haskell packages, employing the very same idea of moduli on the type level, namely modular, modular-arithmetic and finite-field. One can also use finite-typelits, which covers some elementary modular arithmetic as well. Unfortunately, all of them fall behind in terms of performance. Here is a brief comparison:

Discipline	`mod`	`modular`	`modular-arithmetic`	`finite-typelits`	`finite-field`
Addition	Fast	Slow	Slow	Slow	Slow
Small `(*)`	Fast	Slow	Slow	Slow	Slow
Inversion	Fast	N/A	Slow	N/A	Slow
Power	Fast	Slow	Slow	Slow	Slow
Overflows	Safe	Safe	Unsafe	Safe	Safe

Addition. All competing implementations of the modular addition involve divisions, while mod completely avoids this costly operation. This makes a difference even for small numbers; e. g., sum [1..10^7] becomes 5x faster. For larger integers the speed up is even more significant, because the computational complexity of division is not linear.
Small (*). When a modulus fits in a machine word (which is quite a common case on 64-bit architectures), mod implements the modular multiplication as a couple of CPU instructions and neither allocates intermediate arbitrary-precision values, nor calls libgmp at all. For computations like product [1..10^7] this gives a 3x boost to performance in comparison to other libraries.
Inversion. This package relies on libgmp for modular inversions. Even for small arguments it is about 5x faster than the native implementation of modular inversion in modular-arithmetic.
Power. This package relies on libgmp for modular exponentiation. Even for small arguments it is about 2x faster than competitors.
Overflows. At first glance modular-arithmetic is more flexible than mod, because it allows to specify the underlying representation of a modular residue, e. g., Mod Integer 100, Mod Int 100, Mod Word8 100. We argue that this is a dangerous freedom, vulnerable to overflows. For instance, 20 ^ 2 :: Mod Word8 100 returns 44 instead of the expected 0. Even less expected is that 50 :: Mod Word8 300 appears to be 6 (remember that type-level numbers are always Natural).

What is the difference between `mod` and `finite-typelits`?

mod is specifically designed to represent modular residues for mathematical applications (wrapping-around finite numbers) and provides modular inversion and exponentiation.

The main focus of finite-typelits is on non-wrapping-around finite numbers, like indices of arrays in vector-sized. It features a Num instance only for the sake of overloading numeric literals. There is no lawful way to define Num except modular arithmetic, but from finite-typelits' viewpoint this is a by-product.

Citius, altius, fortius!

If you are looking for an ultimate performance and your moduli fit into Word, try Data.Mod.Word, which is a drop-in replacement of Data.Mod, offering better performance and much less allocations.

Benchmarks

Here are some relative benchmarks (less is better), which can be reproduced by running cabal bench.

Discipline	`Data.Mod.Word`	`Data.Mod`	`modular`	`modular-arithmetic`	`finite-typelits`	`finite-field`
Sum	0.44x	1x	16.6x	8.9x	14.7x	14.2x
Product	0.95x	1x	7.8x	4.5x	7.0x	7.0x
Inversion	0.54x	1x	N/A	3.2x	N/A	1.8x
Power	0.29x	1x	2.0x	1.2x	1.4x	1.5x

What's next?

This package was cut out of arithmoi to provide modular arithmetic with a light dependency footprint. This goal certainly limits the scope of the API to the bare minimum. If you need more advanced tools (the Chinese remainder theorem, cyclic groups, modular equations, etc.) please refer to the Math.NumberTheory.Moduli module.

mod's People

Contributors

Stargazers

Watchers

Forkers

incertia konsumlamm

mod's Issues

Add Read and Real instances?

Benchmark with -O2

Compare against finite-field

Add flag to build with integer-simple

It would be nice to have a flag to enable building with integer-simple, even if it was slower.

Reason for asking: HLS builds some static binaries on an alpine distribution that forces integer-simple in GHC.

Write Data.Mod.Word for modulos < 2^64

Implement modular inversion using extended binary gcd (e. g., https://link.springer.com/content/pdf/10.1007%2F3-540-36400-5_6.pdf)

Understanding the `GcdDomain`, `Euclidean`, `Field` instances

GcdDomain

divide returns a Maybe (Mod m), so why does it always return Just (x / y), even when y is not invertible? I also don't really understand why gcd/lcm are defined to always be 1, is this just an arbitrary choice (assuming that every element is invertible)?

Euclidean

Again, why does rem x y always return 0, when not every y is invertible? Is there even a definition of remainder that makes sense for integers modulo m?

Field

ℤ/mℤ is a field if and only if m is prime, so I think that should at least be documented. Adding an IsPrime m constraint would probably be too expensive.

Wrong results with GHC-9.4.4

With GHC-9.4.4 and optimization enabled (-O1 or -O2) I have been observing erroneous results for newtypes over Mod q for $q \geq 2^{64}$.

The behavior is somewhat non-deterministic and doesn't show up all the time. I have been able to reproduce it consistently (both on MacOS and Linux) with the following program:

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE TypeApplications #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Main (main) where

import Data.Mod (Mod)
import GHC.TypeNats (KnownNat)
import Numeric.Natural (Natural)
import Test.QuickCheck (Arbitrary, arbitrary)
import Test.QuickCheck.Instances ()

newtype Fq = Fq (Mod 18446744073709551616 {- 2^64 -}) deriving (Num)

instance KnownNat n => Arbitrary (Mod n) where
    arbitrary = fromIntegral <$> arbitrary @Natural

instance Arbitrary Fq where
    arbitrary = Fq <$> arbitrary

main :: IO ()
main = let Fq x = -1 in print x

(cf. https://gist.github.com/larskuhtz/c216e2fd9e6aac580ea9640d554dddf8)

This program is expected to print 18446744073709551615, but it prints 0 when compiled with optimizations enabled. I tested this with the latest version of the mod package.

The redundant Arbitrary instances are required in this example for the bug to manifest. However, I've also observed the issue in larger production code bases that didn't include Arbitrary instances for Mod and Fq.

Update benchmarks in README

They were generated for GHC 8.10, and results in GHC 9.0 are quite different. Need to wait for GHC 9.4 and regenerate.

Behaviour of `Mod 0`

Mathematically, ℤ/0ℤ is isomorphic to ℤ. I understand that implementing this behaviour here would require a lot of special casing and is probably not desirable.

However, even the way Mod 0 is treated in this package is inconsistent. It is mostly treated as having no values, for example trying to construct a Mod 0 with fromInteger always fails with an exception. However, e.g. zero is unconditionally defined as Mod 0, so that zero :: Mod 0 is a valid (non-bottom) value.

I'm not sure what the best way to address this is. One possibility is to add a bunch of 1 <= m constraints (which is what the modular-arithmetic package does for example). Another possibility would be to simply document that the Mod 0 type shouldn't be used.

Speed up fromInteger

Along these lines:

fromIntegerMod :: Natural -> Integer -> Natural
fromIntegerMod (NatS# m#) (S# x#) =
  if isTrue# (x# >= 0#)
    then NatS# (int2Word x# `remWord#` m#)
    else NatS# ...

instance KnownNat m => Num (Mod m) where
  fromInteger x = mx
    where
      mx = Mod $ fromIntegerMod (natVal mx) x

`Show` instance

I'm confused as to why the show instance always uses parentheses, e.g.:

 >>> -1 :: Mod 10
(9 `modulo` 10)

I'd expect that to only happen when the operator precedence requires it (using showsPrec). Furthermore, why is modulo used here? There is no function of that name in this package or in base. I think mod would be better here, as it mirrors the base function and the math notation. I also would like something like "9 mod 10" (without the backticks), but the convention is to produce (more or less) valid Haskell expressions.

Instances of unboxed vectors

Does the ^% rewrite rule fire?

A simple example:

{-# language ScopedTypeVariables #-}
{-# language RankNTypes #-}

import Data.Mod
import Data.Proxy
import GHC.Natural
import GHC.TypeNats hiding ( Mod(..) )

e,n,msg :: Natural
e = 35954996171344944356088924593542576361517455332243067
n = 53932494257017416534133387405057892857322884506118267
msg = 40511415361412953977729088020603939101623107

main = print $ rsa (e,n) msg

rsa :: (Natural, Natural) -> Natural -> Natural
rsa (e,n) m = case someNatVal n of
  SomeNat (_ :: Proxy p) -> unMod ((fromIntegral m ^ e) :: Mod p)

No matter what optimization flags I use, the exponentiation is always done by GHC.Real.^ (as I can see by inspecting the STG). I expect that the rewrite rules for (^%) should have fired, leading to the exponentiation being done by powModNatural. Are they firing?