arithmoi 0.2.0.5 → 0.2.0.6
raw patch · 3 files changed
+57/−42 lines, 3 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- Changes +3/−0
- Math/NumberTheory/Moduli.hs +53/−41
- arithmoi.cabal +1/−1
Changes view
@@ -1,3 +1,6 @@+0.2.0.6:+ Performance tweaks for powerModInteger (~10%) and+ invertMod (~25%). 0.2.0.5: Fix bug in psieveFrom 0.2.0.4:
Math/NumberTheory/Moduli.hs view
@@ -28,7 +28,6 @@ import Data.Array.Unboxed import Data.Array.Base (unsafeAt) -import Math.NumberTheory.GCD (extendedGCD) import Math.NumberTheory.Utils (shiftToOddCount) -- | Invert a number relative to a modulus.@@ -43,13 +42,26 @@ -- If @gcd number modulus > 1@, the result is @Nothing@. invertMod :: Integer -> Integer -> Maybe Integer invertMod k 0 = if k == 1 || k == (-1) then Just k else Nothing-invertMod k m = case extendedGCD k' m' of- (1, u, _) -> Just (if u < 0 then m' + u else u)- _ -> Nothing+invertMod k m = wrap $ go False 1 0 m' k' where m' = abs m- k' | k >= m' || k < 0 = k `mod` m'- | otherwise = k+ k' | r < 0 = r+m'+ | otherwise = r+ where+ r = k `rem` m'+ wrap x = case (x*k') `rem` m' of+ 1 -> Just x+ _ -> Nothing+ -- Calculate modular inverse of k' modulo m' by continued fraction expansion+ -- of m'/k', say [a_0,a_1,...,a_s]. Let the convergents be p_j/q_j.+ -- Starting from j = -2, the arguments of go are+ -- (p_j/q_j) > m'/k', p_{j+1}, p_j, and n, d with n/d = [a_{j+2},...,a_s].+ -- Since m'/k' = p_s/q_s, and p_j*q_{j+1} - p_{j+1}*q_j = (-1)^(j+1), we have+ -- p_{s-1}*k' - q_{s-1}*m' = (-1)^s * gcd m' k', so if the inverse exists,+ -- it is either p_{s-1} or -p_{s-1}, depending on whether s is even or odd.+ go !b _ po _ 0 = if b then po else (m'-po)+ go b !pn po n d = case n `quotRem` d of+ (q,r) -> go (not b) (q*pn+po) pn d r -- | Jacobi symbol of two numbers. -- The \"denominator\" must be odd and positive, this condition is checked.@@ -191,7 +203,7 @@ -- | Specialised worker without input checks. Makes the same assumptions -- as the general version 'powerMod''. powerModInteger' :: Integer -> Integer -> Integer -> Integer-powerModInteger' base expo md = go e1 w1 1 base+powerModInteger' base expo md = go w1 1 base e1 where w1 = fromInteger expo e1 = expo `shiftR` 64@@ -202,55 +214,55 @@ -- thus it is faster to split each Word64 into the constituent 32-bit -- Words and process those separately. -- The code becomes ugly, unfortunately.- go :: Integer -> Word64 -> Integer -> Integer -> Integer- go 0 !w !a !s = end w a s- go e w a s = inner1 0 a s+ go :: Word64 -> Integer -> Integer -> Integer -> Integer+ go !w !a !s 0 = end a s w+ go w a s e = inner1 a s 0 where wl :: Word !wl = fromIntegral w wh :: Word !wh = fromIntegral (w `shiftR` 32)- inner1 32 !au !sq = inner2 0 au sq- inner1 i au sq- | testBit wl i = inner1 (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)- | otherwise = inner1 (i+1) au ((sq*sq) `rem` md)- inner2 32 !au !sq = go (e `shiftR` 64) (fromInteger e) au sq- inner2 i au sq- | testBit wh i = inner2 (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)- | otherwise = inner2 (i+1) au ((sq*sq) `rem` md)- end w !a !s- | wh == 0 = fin wl a s- | otherwise = innerE 0 a s+ inner1 !au !sq 32 = inner2 au sq 0+ inner1 au sq i+ | testBit wl i = inner1 ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)+ | otherwise = inner1 au ((sq*sq) `rem` md) (i+1)+ inner2 !au !sq 32 = go (fromInteger e) au sq (e `shiftR` 64)+ inner2 au sq i+ | testBit wh i = inner2 ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)+ | otherwise = inner2 au ((sq*sq) `rem` md) (i+1)+ end !a !s w+ | wh == 0 = fin a s wl+ | otherwise = innerE a s 0 where wl :: Word !wl = fromIntegral w wh :: Word !wh = fromIntegral (w `shiftR` 32)- innerE 32 !au !sq = fin wh au sq- innerE i au sq- | testBit wl i = innerE (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)- | otherwise = innerE (i+1) au ((sq*sq) `rem` md)- fin :: Word -> Integer -> Integer -> Integer- fin 1 !a !s = (a*s) `rem` md- fin w a s- | testBit w 0 = fin (w `shiftR` 1) ((a*s) `rem` md) ((s*s) `rem` md)- | otherwise = fin (w `shiftR` 1) a ((s*s) `rem` md)+ innerE !au !sq 32 = fin au sq wh+ innerE au sq i+ | testBit wl i = innerE ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)+ | otherwise = innerE au ((sq*sq) `rem` md) (i+1)+ fin :: Integer -> Integer -> Word -> Integer+ fin !a !s 1 = (a*s) `rem` md+ fin a s w+ | testBit w 0 = fin ((a*s) `rem` md) ((s*s) `rem` md) (w `shiftR` 1)+ | otherwise = fin a ((s*s) `rem` md) (w `shiftR` 1) #else -- WORD_SIZE_IN_BITS == 64, otherwise things wouldn't compile anyway -- Shorter code since we need not split each 64-bit word.- go :: Integer -> Word -> Integer -> Integer -> Integer- go 0 !w !a !s = end w a s- go e w a s = inner 0 a s+ go :: Word -> Integer -> Integer -> Integer -> Integer+ go !w !a !s 0 = end a s w+ go w a s e = inner a s 0 where- inner 64 !au !sq = go (e `shiftR` 64) (fromInteger e) au sq- inner i au sq- | testBit w i = inner (i+1) ((au*sq) `rem` md) ((sq*sq) `rem` md)- | otherwise = inner (i+1) au ((sq*sq) `rem` md)- end 1 !a !s = (a*s) `rem` md- end w a s- | testBit w 0 = end (w `shiftR` 1) ((a*s) `rem` md) ((s*s) `rem` md)- | otherwise = end (w `shiftR` 1) a ((s*s) `rem` md)+ inner !au !sq 64 = go (fromInteger e) au sq (e `shiftR` 64)+ inner au sq i+ | testBit w i = inner ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)+ | otherwise = inner au ((sq*sq) `rem` md) (i+1)+ end !a !s 1 = (a*s) `rem` md+ end a s w+ | testBit w 0 = end ((a*s) `rem` md) ((s*s) `rem` md) (w `shiftR` 1)+ | otherwise = end a ((s*s) `rem` md) (w `shiftR` 1) #endif
arithmoi.cabal view
@@ -1,5 +1,5 @@ name : arithmoi-version : 0.2.0.5+version : 0.2.0.6 cabal-version : >= 1.6 author : Daniel Fischer copyright : (c) 2011 Daniel Fischer