bodigrim / mod Goto Github PK

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.

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.

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.

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

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?

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.

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.

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.

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