From 2c9b0976646af1638292cfb19e8c56751e9fd27b Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Sat, 17 Mar 2018 20:56:41 +0000 Subject: [PATCH 1/8] Improve comments in Math.NumberTheory.Moduli.Jacobi --- Math/NumberTheory/Moduli/Jacobi.hs | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index b1f873f03..4a70f3024 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -1,12 +1,14 @@ -- | -- Module: Math.NumberTheory.Moduli.Jacobi --- Copyright: (c) 2011 Daniel Fischer, 2017 Andrew Lelechenko +-- Copyright: (c) 2011 Daniel Fischer, 2017-2018 Andrew Lelechenko -- Licence: MIT -- Maintainer: Andrew Lelechenko -- Stability: Provisional -- Portability: Non-portable (GHC extensions) -- --- Jacobi symbol. +-- +-- is a generalization of the Legendre symbol, useful for primality +-- testing and integer factorization. -- {-# LANGUAGE CPP #-} @@ -27,7 +29,8 @@ import Data.Semigroup import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils --- | Type for result of 'jacobi'. +-- | Represents three possible values of +-- . data JacobiSymbol = MinusOne | Zero | One deriving (Eq, Ord, Show) @@ -47,11 +50,17 @@ negJS = \case Zero -> Zero One -> MinusOne --- | Jacobi symbol of two numbers. --- The \"denominator\" must be odd and positive, this condition is checked. +-- | of two arguments. +-- The lower argument (\"denominator\") must be odd and positive, +-- this condition is checked. -- --- If both numbers have a common prime factor, the result --- is @0@, otherwise it is ±1. +-- If arguments have a common factor, the result +-- is 'Zero', otherwise it is 'MinusOne' or 'One'. +-- +-- > > jacobi 1001 9911 +-- > Zero -- arguments have a common factor 11 +-- > > jacobi 1001 9907 +-- > MinusOne {-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol, Int -> Int -> JacobiSymbol, Word -> Word -> JacobiSymbol @@ -63,9 +72,8 @@ jacobi a b | b == 1 = One | otherwise = jacobi' a b -- b odd, > 1 --- Invariant: b > 1 and odd --- | Jacobi symbol of two numbers without validity check of --- the \"denominator\". +-- | Similar to 'jacobi', but the condition on the lower argument +-- (\"denominator\") is __not__ checked. {-# SPECIALISE jacobi' :: Integer -> Integer -> JacobiSymbol, Int -> Int -> JacobiSymbol, Word -> Word -> JacobiSymbol From dedb5456c8ec1627b945db3313993206a3e4d01e Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Wed, 21 Mar 2018 22:53:56 +0000 Subject: [PATCH 2/8] Benchmark jacobi and jacobi' --- Math/NumberTheory/Moduli/Jacobi.hs | 3 +++ arithmoi.cabal | 1 + benchmark/Bench.hs | 2 ++ benchmark/Math/NumberTheory/JacobiBench.hs | 26 ++++++++++++++++++++++ 4 files changed, 32 insertions(+) create mode 100644 benchmark/Math/NumberTheory/JacobiBench.hs diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index 4a70f3024..85722a782 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -25,6 +25,7 @@ import Data.Bits #if __GLASGOW_HASKELL__ < 803 import Data.Semigroup #endif +import Numeric.Natural import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils @@ -62,6 +63,7 @@ negJS = \case -- > > jacobi 1001 9907 -- > MinusOne {-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol, + Natural -> Natural -> JacobiSymbol, Int -> Int -> JacobiSymbol, Word -> Word -> JacobiSymbol #-} @@ -75,6 +77,7 @@ jacobi a b -- | Similar to 'jacobi', but the condition on the lower argument -- (\"denominator\") is __not__ checked. {-# SPECIALISE jacobi' :: Integer -> Integer -> JacobiSymbol, + Natural -> Natural -> JacobiSymbol, Int -> Int -> JacobiSymbol, Word -> Word -> JacobiSymbol #-} diff --git a/arithmoi.cabal b/arithmoi.cabal index e59156cc9..5d5c06467 100644 --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -180,6 +180,7 @@ benchmark criterion other-modules: Math.NumberTheory.ArithmeticFunctionsBench Math.NumberTheory.GCDBench + Math.NumberTheory.JacobiBench Math.NumberTheory.MertensBench Math.NumberTheory.PowersBench Math.NumberTheory.PrimesBench diff --git a/benchmark/Bench.hs b/benchmark/Bench.hs index d9db574bc..f7a428f70 100644 --- a/benchmark/Bench.hs +++ b/benchmark/Bench.hs @@ -4,6 +4,7 @@ import Gauge.Main import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions import Math.NumberTheory.GCDBench as GCD +import Math.NumberTheory.JacobiBench as Jacobi import Math.NumberTheory.MertensBench as Mertens import Math.NumberTheory.PowersBench as Powers import Math.NumberTheory.PrimesBench as Primes @@ -13,6 +14,7 @@ import Math.NumberTheory.SieveBlockBench as SieveBlock main = defaultMain [ ArithmeticFunctions.benchSuite , GCD.benchSuite + , Jacobi.benchSuite , Mertens.benchSuite , Powers.benchSuite , Primes.benchSuite diff --git a/benchmark/Math/NumberTheory/JacobiBench.hs b/benchmark/Math/NumberTheory/JacobiBench.hs new file mode 100644 index 000000000..563ac913a --- /dev/null +++ b/benchmark/Math/NumberTheory/JacobiBench.hs @@ -0,0 +1,26 @@ +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.JacobiBench + ( benchSuite + ) where + +import Data.Bits +import Gauge.Main +import Numeric.Natural + +import Math.NumberTheory.Moduli.Jacobi + +doBench :: (Integral a, Bits a) => (a -> a -> JacobiSymbol) -> a -> a +doBench func lim = sum [ x + y | y <- [3, 5 .. lim], x <- [0..y], func x y == One ] + +benchSuite :: Benchmark +benchSuite = bgroup "Jacobi" + [ bench "jacobi/Int" $ nf (doBench jacobi :: Int -> Int) 2000 + , bench "jacobi'/Int" $ nf (doBench jacobi' :: Int -> Int) 2000 + , bench "jacobi/Word" $ nf (doBench jacobi :: Word -> Word) 2000 + , bench "jacobi'/Word" $ nf (doBench jacobi' :: Word -> Word) 2000 + , bench "jacobi/Integer" $ nf (doBench jacobi :: Integer -> Integer) 2000 + , bench "jacobi'/Integer" $ nf (doBench jacobi' :: Integer -> Integer) 2000 + , bench "jacobi/Natural" $ nf (doBench jacobi :: Natural -> Natural) 2000 + , bench "jacobi'/Natural" $ nf (doBench jacobi' :: Natural -> Natural) 2000 + ] From 0bc40842ad708dd958a77dc74211206275534de0 Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Wed, 21 Mar 2018 23:09:38 +0000 Subject: [PATCH 3/8] Remove useless SPECIALISE pragmas --- Math/NumberTheory/Moduli/Jacobi.hs | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index 85722a782..3070069cf 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -106,10 +106,6 @@ jacobi' a b _ -> jacOL One b a -- numerator positive and smaller than denominator -{-# SPECIALISE jacPS :: JacobiSymbol -> Integer -> Integer -> JacobiSymbol, - JacobiSymbol -> Int -> Int -> JacobiSymbol, - JacobiSymbol -> Word -> Word -> JacobiSymbol - #-} jacPS :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol jacPS j a b | evenI a = case shiftToOddCount a of @@ -120,10 +116,6 @@ jacPS j a b | otherwise = jacOL (if rem4 a .&. rem4 b == 3 then (negJS j) else j) b a -- numerator odd, positive and larger than denominator -{-# SPECIALISE jacOL :: JacobiSymbol -> Integer -> Integer -> JacobiSymbol, - JacobiSymbol -> Int -> Int -> JacobiSymbol, - JacobiSymbol -> Word -> Word -> JacobiSymbol - #-} jacOL :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol jacOL j a b | b == 1 = j @@ -135,24 +127,12 @@ jacOL j a b -- For large Integers, going via Int is much faster than bit-fiddling -- on the Integer, so we do that. -{-# SPECIALISE evenI :: Integer -> Bool, - Int -> Bool, - Word -> Bool - #-} evenI :: Integral a => a -> Bool evenI n = fromIntegral n .&. 1 == (0 :: Int) -{-# SPECIALISE rem4 :: Integer -> Int, - Int -> Int, - Word -> Int - #-} rem4 :: Integral a => a -> Int rem4 n = fromIntegral n .&. 3 -{-# SPECIALISE rem8 :: Integer -> Int, - Int -> Int, - Word -> Int - #-} rem8 :: Integral a => a -> Int rem8 n = fromIntegral n .&. 7 From 98ee5fe16d5d392ddd075a48589ad40595d4f8b5 Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Wed, 21 Mar 2018 23:11:25 +0000 Subject: [PATCH 4/8] Rework jacobi utilities --- Math/NumberTheory/Moduli/Jacobi.hs | 88 +++++++++++++----------------- 1 file changed, 38 insertions(+), 50 deletions(-) diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index 3070069cf..07326da26 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -20,14 +20,12 @@ module Math.NumberTheory.Moduli.Jacobi , jacobi' ) where -import Data.Array.Unboxed import Data.Bits #if __GLASGOW_HASKELL__ < 803 import Data.Semigroup #endif import Numeric.Natural -import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils -- | Represents three possible values of @@ -68,11 +66,11 @@ negJS = \case Word -> Word -> JacobiSymbol #-} jacobi :: (Integral a, Bits a) => a -> a -> JacobiSymbol +jacobi _ 1 = One jacobi a b - | b < 0 = error "Math.NumberTheory.Moduli.jacobi: negative denominator" - | evenI b = error "Math.NumberTheory.Moduli.jacobi: even denominator" - | b == 1 = One - | otherwise = jacobi' a b -- b odd, > 1 + | b < 0 = error "Math.NumberTheory.Moduli.jacobi: negative denominator" + | evenI b = error "Math.NumberTheory.Moduli.jacobi: even denominator" + | otherwise = jacobi' a b -- b odd, > 1 -- | Similar to 'jacobi', but the condition on the lower argument -- (\"denominator\") is __not__ checked. @@ -82,59 +80,49 @@ jacobi a b Word -> Word -> JacobiSymbol #-} jacobi' :: (Integral a, Bits a) => a -> a -> JacobiSymbol +jacobi' 0 _ = Zero +jacobi' 1 _ = One jacobi' a b - | a == 0 = Zero - | a == 1 = One - | a < 0 = let n | rem4 b == 1 = One - | otherwise = MinusOne - -- Blech, minBound may pose problems - (z,o) = shiftToOddCount (abs $ toInteger a) - s | evenI z || unsafeAt jac2 (rem8 b) == 1 = n - | otherwise = negJS n - in s <> jacobi' (fromInteger o) b - | a >= b = case a `rem` b of - 0 -> Zero - r -> jacPS One r b - | evenI a = case shiftToOddCount a of - (z,o) -> let r | rem4 o .&. rem4 b == 1 = One - | otherwise = MinusOne - s | evenI z || unsafeAt jac2 (rem8 b) == 1 = r - | otherwise = negJS r - in jacOL s b o - | otherwise = case rem4 a .&. rem4 b of - 3 -> jacOL MinusOne b a - _ -> jacOL One b a + | a < 0 = let n = if rem4is3 b then MinusOne else One + -- Blech, minBound may pose problems + (z, o) = shiftToOddCount (abs $ toInteger a) + s = if evenI z || rem8is1or7 b then n else negJS n + in s <> jacobi' (fromInteger o) b + | a >= b = case a `rem` b of + 0 -> Zero + r -> jacPS One r b + | evenI a = case shiftToOddCount a of + (z, o) -> let r = if rem4is3 o && rem4is3 b then MinusOne else One + s = if evenI z || rem8is1or7 b then r else negJS r + in jacOL s b o + | otherwise = jacOL (if rem4is3 a && rem4is3 b then MinusOne else One) b a -- numerator positive and smaller than denominator jacPS :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol -jacPS j a b - | evenI a = case shiftToOddCount a of - (z,o) | evenI z || unsafeAt jac2 (rem8 b) == 1 -> - jacOL (if rem4 o .&. rem4 b == 3 then (negJS j) else j) b o - | otherwise -> - jacOL (if rem4 o .&. rem4 b == 3 then j else (negJS j)) b o - | otherwise = jacOL (if rem4 a .&. rem4 b == 3 then (negJS j) else j) b a +jacPS acc a b + | evenI a = case shiftToOddCount a of + (z, o) + | evenI z || rem8is1or7 b -> jacOL (if rem4is3 o && rem4is3 b then negJS acc else acc) b o + | otherwise -> jacOL (if rem4is3 o && rem4is3 b then acc else negJS acc) b o + | otherwise = jacOL (if rem4is3 a && rem4is3 b then negJS acc else acc) b a -- numerator odd, positive and larger than denominator jacOL :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol -jacOL j a b - | b == 1 = j - | otherwise = case a `rem` b of - 0 -> Zero - r -> jacPS j r b +jacOL acc _ 1 = acc +jacOL acc a b = case a `rem` b of + 0 -> Zero + r -> jacPS acc r b -- Utilities --- For large Integers, going via Int is much faster than bit-fiddling --- on the Integer, so we do that. -evenI :: Integral a => a -> Bool -evenI n = fromIntegral n .&. 1 == (0 :: Int) +-- Sadly, GHC do not optimise `Prelude.even` to a bit test automatically. +evenI :: Bits a => a -> Bool +evenI n = not (n `testBit` 0) -rem4 :: Integral a => a -> Int -rem4 n = fromIntegral n .&. 3 +-- For an odd input @n@ test whether n `rem` 4 == 1 +rem4is3 :: Bits a => a -> Bool +rem4is3 n = n `testBit` 1 -rem8 :: Integral a => a -> Int -rem8 n = fromIntegral n .&. 7 - -jac2 :: UArray Int Int -jac2 = array (0,7) [(0,0),(1,1),(2,0),(3,-1),(4,0),(5,-1),(6,0),(7,1)] +-- For an odd input @n@ test whether (n `rem` 8) `elem` [1, 7] +rem8is1or7 :: Bits a => a -> Bool +rem8is1or7 n = n `testBit` 1 == n `testBit` 2 From 9bf60c33a8d89d35e54373b7f57225945a3f80d1 Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Thu, 22 Mar 2018 21:56:29 +0000 Subject: [PATCH 5/8] Make accumulator strict --- Math/NumberTheory/Moduli/Jacobi.hs | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index 07326da26..5a625b78a 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -11,8 +11,9 @@ -- testing and integer factorization. -- -{-# LANGUAGE CPP #-} -{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE LambdaCase #-} module Math.NumberTheory.Moduli.Jacobi ( JacobiSymbol(..) @@ -99,7 +100,7 @@ jacobi' a b -- numerator positive and smaller than denominator jacPS :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol -jacPS acc a b +jacPS !acc a b | evenI a = case shiftToOddCount a of (z, o) | evenI z || rem8is1or7 b -> jacOL (if rem4is3 o && rem4is3 b then negJS acc else acc) b o @@ -108,8 +109,8 @@ jacPS acc a b -- numerator odd, positive and larger than denominator jacOL :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol -jacOL acc _ 1 = acc -jacOL acc a b = case a `rem` b of +jacOL !acc _ 1 = acc +jacOL !acc a b = case a `rem` b of 0 -> Zero r -> jacPS acc r b From f8bc5afd98c5d5e92b000b091a7d617f730d926e Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Thu, 22 Mar 2018 22:02:19 +0000 Subject: [PATCH 6/8] Do not expose jacobi' with unchecked input --- Math/NumberTheory/Moduli/Jacobi.hs | 8 -------- Math/NumberTheory/Moduli/Sqrt.hs | 4 ++-- Math/NumberTheory/Primes/Testing/Probabilistic.hs | 2 +- benchmark/Math/NumberTheory/JacobiBench.hs | 4 ---- 4 files changed, 3 insertions(+), 15 deletions(-) diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index 5a625b78a..63c0530ed 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -18,7 +18,6 @@ module Math.NumberTheory.Moduli.Jacobi ( JacobiSymbol(..) , jacobi - , jacobi' ) where import Data.Bits @@ -73,13 +72,6 @@ jacobi a b | evenI b = error "Math.NumberTheory.Moduli.jacobi: even denominator" | otherwise = jacobi' a b -- b odd, > 1 --- | Similar to 'jacobi', but the condition on the lower argument --- (\"denominator\") is __not__ checked. -{-# SPECIALISE jacobi' :: Integer -> Integer -> JacobiSymbol, - Natural -> Natural -> JacobiSymbol, - Int -> Int -> JacobiSymbol, - Word -> Word -> JacobiSymbol - #-} jacobi' :: (Integral a, Bits a) => a -> a -> JacobiSymbol jacobi' 0 _ = Zero jacobi' 1 _ = One diff --git a/Math/NumberTheory/Moduli/Sqrt.hs b/Math/NumberTheory/Moduli/Sqrt.hs index c8dd5b924..66d85be06 100644 --- a/Math/NumberTheory/Moduli/Sqrt.hs +++ b/Math/NumberTheory/Moduli/Sqrt.hs @@ -43,7 +43,7 @@ import Math.NumberTheory.Utils (shiftToOddCount, splitOff) -- the result is @Nothing@. sqrtModP :: Integer -> Integer -> Maybe Integer sqrtModP n 2 = Just (n `mod` 2) -sqrtModP n prime = case jacobi' n prime of +sqrtModP n prime = case jacobi n prime of MinusOne -> Nothing Zero -> Just 0 One -> Just (sqrtModP' (n `mod` prime) prime) @@ -215,7 +215,7 @@ findNonSquare n where primelist = [3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67] ++ sieveFrom (68 + n `rem` 4) -- prevent sharing - search (p:ps) = case jacobi' p n of + search (p:ps) = case jacobi p n of MinusOne -> p _ -> search ps search _ = error "Should never have happened, prime list exhausted." diff --git a/Math/NumberTheory/Primes/Testing/Probabilistic.hs b/Math/NumberTheory/Primes/Testing/Probabilistic.hs index 69c931a67..bee231c95 100644 --- a/Math/NumberTheory/Primes/Testing/Probabilistic.hs +++ b/Math/NumberTheory/Primes/Testing/Probabilistic.hs @@ -156,7 +156,7 @@ lucasTest n square = isPossibleSquare2 n && r*r == n r = integerSquareRoot n d = find True 5 - find !pos cd = case jacobi' (n `rem` cd) cd of + find !pos cd = case jacobi (n `rem` cd) cd of MinusOne -> if pos then cd else (-cd) Zero -> if cd == n then 1 else 0 One -> find (not pos) (cd+2) diff --git a/benchmark/Math/NumberTheory/JacobiBench.hs b/benchmark/Math/NumberTheory/JacobiBench.hs index 563ac913a..186be668b 100644 --- a/benchmark/Math/NumberTheory/JacobiBench.hs +++ b/benchmark/Math/NumberTheory/JacobiBench.hs @@ -16,11 +16,7 @@ doBench func lim = sum [ x + y | y <- [3, 5 .. lim], x <- [0..y], func x y == On benchSuite :: Benchmark benchSuite = bgroup "Jacobi" [ bench "jacobi/Int" $ nf (doBench jacobi :: Int -> Int) 2000 - , bench "jacobi'/Int" $ nf (doBench jacobi' :: Int -> Int) 2000 , bench "jacobi/Word" $ nf (doBench jacobi :: Word -> Word) 2000 - , bench "jacobi'/Word" $ nf (doBench jacobi' :: Word -> Word) 2000 , bench "jacobi/Integer" $ nf (doBench jacobi :: Integer -> Integer) 2000 - , bench "jacobi'/Integer" $ nf (doBench jacobi' :: Integer -> Integer) 2000 , bench "jacobi/Natural" $ nf (doBench jacobi :: Natural -> Natural) 2000 - , bench "jacobi'/Natural" $ nf (doBench jacobi' :: Natural -> Natural) 2000 ] From 8f7600a17b10fcaa3f94cac0d855af5866c0837e Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Thu, 22 Mar 2018 22:12:47 +0000 Subject: [PATCH 7/8] Specialize shiftToOddCount for Natural --- Math/NumberTheory/Utils.hs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs index cc3f63f7c..d4b831365 100644 --- a/Math/NumberTheory/Utils.hs +++ b/Math/NumberTheory/Utils.hs @@ -29,6 +29,7 @@ import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals +import GHC.Natural import Data.Bits @@ -42,6 +43,7 @@ uncheckedShiftR (W# w#) (I# i#) = W# (uncheckedShiftRL# w# i#) "shiftToOddCount/Int" shiftToOddCount = shiftOCInt "shiftToOddCount/Word" shiftToOddCount = shiftOCWord "shiftToOddCount/Integer" shiftToOddCount = shiftOCInteger +"shiftToOddCount/Natural" shiftToOddCount = shiftOCNatural #-} {-# INLINE [1] shiftToOddCount #-} shiftToOddCount :: Integral a => a -> (Int, a) @@ -75,6 +77,18 @@ shiftOCInteger n@(Jn# bn#) = case bigNatZeroCount bn# of 0# -> (0, n) z# -> (I# z#, n `shiftRInteger` z#) +-- | Specialised version for @'Natural'@. +-- Precondition: argument nonzero (not checked). +shiftOCNatural :: Natural -> (Int, Natural) +shiftOCNatural n@(NatS# i#) = + case shiftToOddCount# i# of + (# z#, w# #) + | isTrue# (z# ==# 0#) -> (0, n) + | otherwise -> (I# z#, NatS# w#) +shiftOCNatural n@(NatJ# bn#) = case bigNatZeroCount bn# of + 0# -> (0, n) + z# -> (I# z#, NatJ# (bn# `shiftRBigNat` z#)) + -- | Count trailing zeros in a @'BigNat'@. -- Precondition: argument nonzero (not checked, Integer invariant). bigNatZeroCount :: BigNat -> Int# From 5f9b43d8b99bc99f1d7e95abd5f94f26edb106db Mon Sep 17 00:00:00 2001 From: Bodigrim Date: Wed, 11 Apr 2018 23:39:26 +0100 Subject: [PATCH 8/8] More tests, including minBound edge case --- Math/NumberTheory/Moduli/Jacobi.hs | 5 +- .../Math/NumberTheory/Moduli/JacobiTests.hs | 70 +++++++++++++++++-- 2 files changed, 68 insertions(+), 7 deletions(-) diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs index 63c0530ed..befdb08db 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -77,10 +77,9 @@ jacobi' 0 _ = Zero jacobi' 1 _ = One jacobi' a b | a < 0 = let n = if rem4is3 b then MinusOne else One - -- Blech, minBound may pose problems - (z, o) = shiftToOddCount (abs $ toInteger a) + (z, o) = shiftToOddCount (negate a) s = if evenI z || rem8is1or7 b then n else negJS n - in s <> jacobi' (fromInteger o) b + in s <> jacobi' o b | a >= b = case a `rem` b of 0 -> Zero r -> jacPS One r b diff --git a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs index 04e3348b6..f321c52cd 100644 --- a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs @@ -23,6 +23,7 @@ import Data.Bits #if __GLASGOW_HASKELL__ < 803 import Data.Semigroup #endif +import Numeric.Natural import Math.NumberTheory.Moduli hiding (invertMod) import Math.NumberTheory.TestUtils @@ -40,29 +41,90 @@ jacobiProperty3 (AnySign a) (MyCompose (Positive (Odd n))) = case jacobi a n of Zero -> a `gcd` n /= 1 One -> a `gcd` n == 1 +doesProductOverflow :: Integral a => a -> a -> Bool +doesProductOverflow x y = abs (toInteger (x * y)) < abs (toInteger x * toInteger y) + -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 4 jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> (MyCompose Positive Odd) a -> Bool -jacobiProperty4 (AnySign a) (AnySign b) (MyCompose (Positive (Odd n))) = jacobi (a * b) n == jacobi a n <> jacobi b n +jacobiProperty4 (AnySign a) (AnySign b) (MyCompose (Positive (Odd n))) = + doesProductOverflow a b || + jacobi (a * b) n == jacobi a n <> jacobi b n + +jacobiProperty4_Int :: AnySign Int -> AnySign Int -> (MyCompose Positive Odd) Int -> Bool +jacobiProperty4_Int = jacobiProperty4 + +jacobiProperty4_Word :: AnySign Word -> AnySign Word -> (MyCompose Positive Odd) Word -> Bool +jacobiProperty4_Word = jacobiProperty4 jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> (MyCompose Positive Odd) Integer -> Bool jacobiProperty4_Integer = jacobiProperty4 +jacobiProperty4_Natural :: AnySign Natural -> AnySign Natural -> (MyCompose Positive Odd) Natural -> Bool +jacobiProperty4_Natural = jacobiProperty4 + -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 5 jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> (MyCompose Positive Odd) a -> (MyCompose Positive Odd) a -> Bool -jacobiProperty5 (AnySign a) (MyCompose (Positive (Odd m))) (MyCompose (Positive (Odd n))) = jacobi a (m * n) == jacobi a m <> jacobi a n +jacobiProperty5 (AnySign a) (MyCompose (Positive (Odd m))) (MyCompose (Positive (Odd n))) = + doesProductOverflow m n || + jacobi a (m * n) == jacobi a m <> jacobi a n + +jacobiProperty5_Int :: AnySign Int -> (MyCompose Positive Odd) Int -> (MyCompose Positive Odd) Int -> Bool +jacobiProperty5_Int = jacobiProperty5 + +jacobiProperty5_Word :: AnySign Word -> (MyCompose Positive Odd) Word -> (MyCompose Positive Odd) Word -> Bool +jacobiProperty5_Word = jacobiProperty5 jacobiProperty5_Integer :: AnySign Integer -> (MyCompose Positive Odd) Integer -> (MyCompose Positive Odd) Integer -> Bool jacobiProperty5_Integer = jacobiProperty5 +jacobiProperty5_Natural :: AnySign Natural -> (MyCompose Positive Odd) Natural -> (MyCompose Positive Odd) Natural -> Bool +jacobiProperty5_Natural = jacobiProperty5 + -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 6 jacobiProperty6 :: (Integral a, Bits a) => (MyCompose Positive Odd) a -> (MyCompose Positive Odd) a -> Bool jacobiProperty6 (MyCompose (Positive (Odd m))) (MyCompose (Positive (Odd n))) = gcd m n /= 1 || jacobi m n <> jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then One else MinusOne) +-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 7 +jacobiProperty7 :: (Integral a, Bits a) => (MyCompose Positive Odd) a -> Bool +jacobiProperty7 (MyCompose (Positive (Odd n))) = + jacobi (-1) n == if n `mod` 4 == 1 then One else MinusOne + +jacobiProperty7_Int :: (MyCompose Positive Odd) Int -> Bool +jacobiProperty7_Int = jacobiProperty7 + +jacobiProperty7_Integer :: (MyCompose Positive Odd) Integer -> Bool +jacobiProperty7_Integer = jacobiProperty7 + +-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 8 +jacobiProperty8 :: (Integral a, Bits a) => (MyCompose Positive Odd) a -> Bool +jacobiProperty8 (MyCompose (Positive (Odd n))) = + even n || + jacobi 2 n == if n `mod` 8 == 1 || n `mod` 8 == 7 then One else MinusOne + +jacobiProperty9 :: (Integral a, Bits a, Bounded a) => (MyCompose Positive Odd) a -> Bool +jacobiProperty9 (MyCompose (Positive (Odd n))) = + jacobi m n == jacobi (toInteger m) (toInteger n) + where + m = minBound + +jacobiProperty9_Int :: (MyCompose Positive Odd) Int -> Bool +jacobiProperty9_Int = jacobiProperty9 + testSuite :: TestTree testSuite = testGroup "Jacobi" [ testSameIntegralProperty "same modulo n" jacobiProperty2 , testSameIntegralProperty "consistent with gcd" jacobiProperty3 - , testSmallAndQuick "multiplicative 1" jacobiProperty4_Integer - , testSmallAndQuick "multiplicative 2" jacobiProperty5_Integer + , testSmallAndQuick "multiplicative 1 Int" jacobiProperty4_Int + , testSmallAndQuick "multiplicative 1 Word" jacobiProperty4_Word + , testSmallAndQuick "multiplicative 1 Integer" jacobiProperty4_Integer + , testSmallAndQuick "multiplicative 1 Natural" jacobiProperty4_Natural + , testSmallAndQuick "multiplicative 2 Int" jacobiProperty5_Int + , testSmallAndQuick "multiplicative 2 Word" jacobiProperty5_Word + , testSmallAndQuick "multiplicative 2 Integer" jacobiProperty5_Integer + , testSmallAndQuick "multiplicative 2 Natural" jacobiProperty5_Natural , testSameIntegralProperty "law of quadratic reciprocity" jacobiProperty6 + , testSmallAndQuick "-1 Int" jacobiProperty7_Int + , testSmallAndQuick "-1 Integer" jacobiProperty7_Integer + , testIntegralProperty "2" jacobiProperty8 + , testSmallAndQuick "minBound Int" jacobiProperty9_Int ]