packages feed

integer-logarithms (empty) → 1

raw patch · 13 files changed

+1072/−0 lines, 13 filesdep +QuickCheckdep +arraydep +basesetup-changed

Dependencies added: QuickCheck, array, base, ghc-prim, integer-gmp, integer-logarithms, nats, smallcheck, tasty, tasty-hunit, tasty-quickcheck, tasty-smallcheck

Files

+ LICENSE view
@@ -0,0 +1,16 @@+Copyright (c) 2011 Daniel Fischer, 2017 Oleg Grenrus++Permission is hereby granted, free of charge, to any person obtaining a copy of this software and+ associated documentation files (the "Software"), to deal in the Software without restriction,+ including without limitation the rights to use, copy, modify, merge, publish, distribute,+ sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is+ furnished to do so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all copies or+substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT+LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+ IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN+ CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ Setup.hs view
@@ -0,0 +1,5 @@+module Main where++import Distribution.Simple++main = defaultMain
+ integer-logarithms.cabal view
@@ -0,0 +1,101 @@+name:               integer-logarithms+version:            1+cabal-version:      >= 1.10+author:             Daniel Fischer+copyright:          (c) 2011 Daniel Fischer+license:            MIT+license-file:       LICENSE+maintainer:         Oleg Grenrus <oleg.grenrus@iki.fi>+build-type:         Simple+stability:          Provisional+homepage:           https://github.com/phadej/integer-logarithms+bug-reports:        https://github.com/phadej/integer-logarithms/issues++synopsis:           Integer logarithms.+description:+  "Math.NumberTheory.Logarithms" and "Math.NumberTheory.Powers.Integer"+  from the arithmoi package.+  .+  Also provides "GHC.Integer.Logarithms.Compat" and+  "Math.NumberTheory.Power.Natural" modules, as well as some+  additional functions in migrated modules.++category:           Math, Algorithms, Number Theory++tested-with         :+  GHC==7.0.4,+  GHC==7.2.2,+  GHC==7.4.2,+  GHC==7.6.3,+  GHC==7.8.4,+  GHC==7.10.3,+  GHC==8.0.1++extra-source-files  : readme.md++flag check-bounds+  description:  Replace unsafe array operations with safe ones+  default:      False+  manual:       True++library+  default-language: Haskell2010+  hs-source-dirs: src+  build-depends:+    base >= 4.3 && < 4.10,+    array >= 0.3 && < 0.6,+    ghc-prim < 0.6,+    integer-gmp < 1.1+  if impl(ghc >= 7.10)+    cpp-options: -DBase48+  else+    build-depends: nats >= 1.1 && <1.2++  exposed-modules:+    Math.NumberTheory.Logarithms+    Math.NumberTheory.Powers.Integer+    Math.NumberTheory.Powers.Natural+    GHC.Integer.Logarithms.Compat+  other-extensions:+    BangPatterns+    CPP+    MagicHash++  ghc-options: -O2 -Wall+  if flag(check-bounds)+    cpp-options: -DCheckBounds++source-repository head+  type:     git+  location: https://github.com/phadej/integer-logarithms++test-suite spec+  type:                 exitcode-stdio-1.0+  hs-source-dirs:       test-suite+  ghc-options:          -Wall+  main-is:              Test.hs+  default-language:     Haskell2010+  other-extensions:+    StandaloneDeriving+    FlexibleContexts+    FlexibleInstances+    GeneralizedNewtypeDeriving+    MultiParamTypeClasses+  build-depends:+    base,+    integer-logarithms,+    tasty >= 0.10 && < 0.12,+    tasty-smallcheck >= 0.8 && < 0.9,+    tasty-quickcheck >= 0.8 && < 0.9,+    tasty-hunit >= 0.9 && < 0.10,+    QuickCheck >= 2.9 && < 2.10,+    smallcheck >= 1.1 && < 1.2+  if !impl(ghc >= 7.10)+    build-depends: nats >= 1.1 && <1.2++  other-modules:+    Math.NumberTheory.LogarithmsTests+    Math.NumberTheory.Powers.IntegerTests+    Math.NumberTheory.Powers.NaturalTests+  other-modules:+    Math.NumberTheory.TestUtils
+ readme.md view
@@ -0,0 +1,3 @@+# integer-logarithms++`Math.NumberTheory.Logarithms` splitted out of [`arithmoi`](http://hackage.haskell.org/package/arithmoi)
+ src/GHC/Integer/Logarithms/Compat.hs view
@@ -0,0 +1,152 @@+-- |+-- Module:      GHC.Integer.Logarithms.Compat+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Low level stuff for integer logarithms.+{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-}+module GHC.Integer.Logarithms.Compat+    ( -- * Functions+      integerLogBase#+    , integerLog2#+    , wordLog2#+    ) where++#if __GLASGOW_HASKELL__ >= 702++-- Stuff is already there+import GHC.Integer.Logarithms++#else++-- We have to define it here+#include "MachDeps.h"++import GHC.Base+import GHC.Integer.GMP.Internals++#if (WORD_SIZE_IN_BITS != 32) && (WORD_SIZE_IN_BITS != 64)+#error Only word sizes 32 and 64 are supported.+#endif+++#if WORD_SIZE_IN_BITS == 32++#define WSHIFT 5+#define MMASK 31++#else++#define WSHIFT 6+#define MMASK 63++#endif++-- Reference implementation only, the algorithm in M.NT.Logarithms is better.++-- | Calculate the integer logarithm for an arbitrary base.+--   The base must be greater than 1, the second argument, the number+--   whose logarithm is sought; should be positive, otherwise the+--   result is meaningless.+--+-- > base ^ integerLogBase# base m <= m < base ^ (integerLogBase# base m + 1)+--+-- for @base > 1@ and @m > 0@.+integerLogBase# :: Integer -> Integer -> Int#+integerLogBase# b m = case step b of+                        (# _, e #) -> e+  where+    step pw =+      if m < pw+        then (# m, 0# #)+        else case step (pw * pw) of+               (# q, e #) ->+                 if q < pw+                   then (# q, 2# *# e #)+                   else (# q `quot` pw, 2# *# e +# 1# #)++-- | Calculate the integer base 2 logarithm of an 'Integer'.+--   The calculation is much more efficient than for the general case.+--+--   The argument must be strictly positive, that condition is /not/ checked.+integerLog2# :: Integer -> Int#+integerLog2# (S# i) = wordLog2# (int2Word# i)+integerLog2# (J# s ba) = check (s -# 1#)+  where+    check i = case indexWordArray# ba i of+                0## -> check (i -# 1#)+                w   -> wordLog2# w +# (uncheckedIShiftL# i WSHIFT#)++-- | This function calculates the integer base 2 logarithm of a 'Word#'.+--   @'wordLog2#' 0## = -1#@.+{-# INLINE wordLog2# #-}+wordLog2# :: Word# -> Int#+wordLog2# w =+  case leadingZeros of+   BA lz ->+    let zeros u = indexInt8Array# lz (word2Int# u) in+#if WORD_SIZE_IN_BITS == 64+    case uncheckedShiftRL# w 56# of+     a ->+      if a `neWord#` 0##+       then 64# -# zeros a+       else+        case uncheckedShiftRL# w 48# of+         b ->+          if b `neWord#` 0##+           then 56# -# zeros b+           else+            case uncheckedShiftRL# w 40# of+             c ->+              if c `neWord#` 0##+               then 48# -# zeros c+               else+                case uncheckedShiftRL# w 32# of+                 d ->+                  if d `neWord#` 0##+                   then 40# -# zeros d+                   else+#endif+                    case uncheckedShiftRL# w 24# of+                     e ->+                      if e `neWord#` 0##+                       then 32# -# zeros e+                       else+                        case uncheckedShiftRL# w 16# of+                         f ->+                          if f `neWord#` 0##+                           then 24# -# zeros f+                           else+                            case uncheckedShiftRL# w 8# of+                             g ->+                              if g `neWord#` 0##+                               then 16# -# zeros g+                               else 8# -# zeros w++-- Lookup table+data BA = BA ByteArray#++leadingZeros :: BA+leadingZeros =+    let mkArr s =+          case newByteArray# 256# s of+            (# s1, mba #) ->+              case writeInt8Array# mba 0# 9# s1 of+                s2 ->+                  let fillA lim val idx st =+                        if idx ==# 256#+                          then st+                          else if idx <# lim+                                then case writeInt8Array# mba idx val st of+                                        nx -> fillA lim val (idx +# 1#) nx+                                else fillA (2# *# lim) (val -# 1#) idx st+                  in case fillA 2# 8# 1# s2 of+                      s3 -> case unsafeFreezeByteArray# mba s3 of+                              (# _, ba #) -> ba+    in case mkArr realWorld# of+        b -> BA b++#endif
+ src/Math/NumberTheory/Logarithms.hs view
@@ -0,0 +1,330 @@+-- |+-- Module:      Math.NumberTheory.Logarithms+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Integer Logarithms. For efficiency, the internal representation of 'Integer's+-- from integer-gmp is used.+--+{-# LANGUAGE CPP #-}+{-# LANGUAGE MagicHash #-}+module Math.NumberTheory.Logarithms+    ( -- * Integer logarithms with input checks+      integerLogBase+    , integerLog2+    , integerLog10++    , naturalLogBase+    , naturalLog2+    , naturalLog10++    , intLog2+    , wordLog2++      -- * Integer logarithms without input checks+      --+      -- | These functions are total, however, don't rely on the values with out-of-domain arguments.+    , integerLogBase'+    , integerLog2'+    , integerLog10'++    , intLog2'+    , wordLog2'+    ) where++import GHC.Exts++import Data.Bits+import Data.Array.Unboxed+import Numeric.Natural++import GHC.Integer.Logarithms.Compat+#if Base48+import GHC.Integer.GMP.Internals (Integer (..))+import GHC.Natural+#endif++#if CheckBounds+import Data.Array.IArray (IArray, (!))+#else+import Data.Array.Base (unsafeAt)+#endif++import Math.NumberTheory.Powers.Integer+import Math.NumberTheory.Powers.Natural++-- | Calculate the integer logarithm for an arbitrary base.+--   The base must be greater than 1, the second argument, the number+--   whose logarithm is sought, must be positive, otherwise an error is thrown.+--   If @base == 2@, the specialised version is called, which is more+--   efficient than the general algorithm.+--+--   Satisfies:+--+-- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1)+--+-- for @base > 1@ and @m > 0@.+integerLogBase :: Integer -> Integer -> Int+integerLogBase b n+  | n < 1       = error "Math.NumberTheory.Logarithms.integerLogBase: argument must be positive."+  | n < b       = 0+  | b == 2      = integerLog2' n+  | b < 2       = error "Math.NumberTheory.Logarithms.integerLogBase: base must be greater than one."+  | otherwise   = integerLogBase' b n++-- | Calculate the integer logarithm of an 'Integer' to base 2.+--   The argument must be positive, otherwise an error is thrown.+integerLog2 :: Integer -> Int+integerLog2 n+  | n < 1       = error "Math.NumberTheory.Logarithms.integerLog2: argument must be positive"+  | otherwise   = I# (integerLog2# n)++-- | Cacluate the integer logarithm for an arbitrary base.+--   The base must be greater than 1, the second argument, the number+--   whose logarithm is sought, must be positive, otherwise an error is thrown.+--   If @base == 2@, the specialised version is called, which is more+--   efficient than the general algorithm.+--+--   Satisfies:+--+-- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1)+--+-- for @base > 1@ and @m > 0@.+naturalLogBase :: Natural -> Natural -> Int+naturalLogBase b n+  | n < 1       = error "Math.NumberTheory.Logarithms.naturalLogBase: argument must be positive."+  | n < b       = 0+  | b == 2      = naturalLog2' n+  | b < 2       = error "Math.NumberTheory.Logarithms.naturalLogBase: base must be greater than one."+  | otherwise   = naturalLogBase' b n++-- | Calculate the natural logarithm of an 'Natural' to base 2.+--   The argument must be non-zero, otherwise an error is thrown.+naturalLog2 :: Natural -> Int+naturalLog2 n+  | n < 1       = error "Math.NumberTheory.Logarithms.naturalLog2: argument must be non-zero"+  | otherwise   = I# (naturalLog2# n)++-- | Calculate the integer logarithm of an 'Int' to base 2.+--   The argument must be positive, otherwise an error is thrown.+intLog2 :: Int -> Int+intLog2 (I# i#)+  | isTrue# (i# <# 1#)  = error "Math.NumberTheory.Logarithms.intLog2: argument must be positive"+  | otherwise           = I# (wordLog2# (int2Word# i#))++-- | Calculate the integer logarithm of a 'Word' to base 2.+--   The argument must be positive, otherwise an error is thrown.+wordLog2 :: Word -> Int+wordLog2 (W# w#)+  | isTrue# (w# `eqWord#` 0##)  = error "Math.NumberTheory.Logarithms.wordLog2: argument must not be 0."+  | otherwise                   = I# (wordLog2# w#)++-- | Same as 'integerLog2', but without checks, saves a little time when+--   called often for known good input.+integerLog2' :: Integer -> Int+integerLog2' n = I# (integerLog2# n)++-- | Same as 'naturalLog2', but without checks, saves a little time when+--   called often for known good input.+naturalLog2' :: Natural -> Int+naturalLog2' n = I# (naturalLog2# n)++-- | Same as 'intLog2', but without checks, saves a little time when+--   called often for known good input.+intLog2' :: Int -> Int+intLog2' (I# i#) = I# (wordLog2# (int2Word# i#))++-- | Same as 'wordLog2', but without checks, saves a little time when+--   called often for known good input.+wordLog2' :: Word -> Int+wordLog2' (W# w#) = I# (wordLog2# w#)++-- | Calculate the integer logarithm of an 'Integer' to base 10.+--   The argument must be positive, otherwise an error is thrown.+integerLog10 :: Integer -> Int+integerLog10 n+  | n < 1     = error "Math.NumberTheory.Logarithms.integerLog10: argument must be positive"+  | otherwise = integerLog10' n++-- | Calculate the integer logarithm of an 'Integer' to base 10.+--   The argument must be not zero, otherwise an error is thrown.+naturalLog10 :: Natural -> Int+naturalLog10 n+  | n < 1     = error "Math.NumberTheory.Logarithms.naturalaLog10: argument must be non-zero"+  | otherwise = naturalLog10' n++-- | Same as 'integerLog10', but without a check for a positive+--   argument. Saves a little time when called often for known good+--   input.+integerLog10' :: Integer -> Int+integerLog10' n+  | n < 10      = 0+  | n < 100     = 1+  | otherwise   = ex + integerLog10' (n `quot` integerPower 10 ex)+    where+      ln = I# (integerLog2# n)+      -- u/v is a good approximation of log 2/log 10+      u  = 1936274+      v  = 6432163+      -- so ex is a good approximation to integerLogBase 10 n+      ex = fromInteger ((u * fromIntegral ln) `quot` v)++-- | Same as 'naturalLog10', but without a check for a positive+--   argument. Saves a little time when called often for known good+--   input.+naturalLog10' :: Natural -> Int+naturalLog10' n+  | n < 10      = 0+  | n < 100     = 1+  | otherwise   = ex + naturalLog10' (n `quot` naturalPower 10 ex)+    where+      ln = I# (naturalLog2# n)+      -- u/v is a good approximation of log 2/log 10+      u  = 1936274+      v  = 6432163+      -- so ex is a good approximation to naturalLogBase 10 n+      ex = fromInteger ((u * fromIntegral ln) `quot` v)++-- | Same as 'integerLogBase', but without checks, saves a little time when+--   called often for known good input.+integerLogBase' :: Integer -> Integer -> Int+integerLogBase' b n+  | n < b       = 0+  | ln-lb < lb  = 1     -- overflow safe version of ln < 2*lb, implies n < b*b+  | b < 33      = let bi = fromInteger b+                      ix = 2*bi-4+                      -- u/v is a good approximation of log 2/log b+                      u  = logArr `unsafeAt` ix+                      v  = logArr `unsafeAt` (ix+1)+                      -- hence ex is a rather good approximation of integerLogBase b n+                      -- most of the time, it will already be exact+                      ex = fromInteger ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v)+                  in case u of+                      1 -> ln `quot` v      -- a power of 2, easy+                      _ -> ex + integerLogBase' b (n `quot` integerPower b ex)+  | otherwise   = let -- shift b so that 16 <= bi < 32+                      bi = fromInteger (b `shiftR` (lb-4))+                      -- we choose an approximation of log 2 / log (bi+1) to+                      -- be sure we underestimate+                      ix = 2*bi-2+                      -- u/w is a reasonably good approximation to log 2/log b+                      -- it is too small, but not by much, so the recursive call+                      -- should most of the time be caught by one of the first+                      -- two guards unless n is huge, but then it'd still be+                      -- a call with a much smaller second argument.+                      u  = fromIntegral $ logArr `unsafeAt` ix+                      v  = fromIntegral $ logArr `unsafeAt` (ix+1)+                      w  = v + u*fromIntegral (lb-4)+                      ex = fromInteger ((u * fromIntegral ln) `quot` w)+                  in ex + integerLogBase' b (n `quot` integerPower b ex)+    where+      lb = integerLog2' b+      ln = integerLog2' n++-- | Same as 'naturalLogBase', but without checks, saves a little time when+--   called often for known good input.+naturalLogBase' :: Natural -> Natural -> Int+naturalLogBase' b n+    | n < b       = 0+  | ln-lb < lb  = 1     -- overflow safe version of ln < 2*lb, implies n < b*b+  | b < 33      = let bi = fromIntegral b+                      ix = 2*bi-4+                      -- u/v is a good approximation of log 2/log b+                      u  = logArr `unsafeAt` ix+                      v  = logArr `unsafeAt` (ix+1)+                      -- hence ex is a rather good approximation of integerLogBase b n+                      -- most of the time, it will already be exact+                      ex = fromNatural ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v)+                  in case u of+                      1 -> ln `quot` v      -- a power of 2, easy+                      _ -> ex + naturalLogBase' b (n `quot` naturalPower b ex)+  | otherwise   = let -- shift b so that 16 <= bi < 32+                      bi = fromNatural (b `shiftR` (lb-4))+                      -- we choose an approximation of log 2 / log (bi+1) to+                      -- be sure we underestimate+                      ix = 2*bi-2+                      -- u/w is a reasonably good approximation to log 2/log b+                      -- it is too small, but not by much, so the recursive call+                      -- should most of the time be caught by one of the first+                      -- two guards unless n is huge, but then it'd still be+                      -- a call with a much smaller second argument.+                      u  = fromIntegral $ logArr `unsafeAt` ix+                      v  = fromIntegral $ logArr `unsafeAt` (ix+1)+                      w  = v + u*fromIntegral (lb-4)+                      ex = fromNatural ((u * fromIntegral ln) `quot` w)+                  in ex + naturalLogBase' b (n `quot` naturalPower b ex)+    where+      lb = naturalLog2' b+      ln = naturalLog2' n++-- Lookup table for logarithms of 2 <= k <= 32+-- In each row "x , y", x/y is a good rational approximation of log 2  / log k.+-- For the powers of 2, it is exact, otherwise x/y < log 2/log k, since we don't+-- want to overestimate integerLogBase b n = floor $ (log 2/log b)*logBase 2 n.+logArr :: UArray Int Int+logArr = listArray (0, 61)+          [ 1 , 1,+            190537 , 301994,+            1 , 2,+            1936274 , 4495889,+            190537 , 492531,+            91313 , 256348,+            1 , 3,+            190537 , 603988,+            1936274 , 6432163,+            1686227 , 5833387,+            190537 , 683068,+            5458 , 20197,+            91313 , 347661,+            416263 , 1626294,+            1 , 4,+            32631 , 133378,+            190537 , 794525,+            163451 , 694328,+            1936274 , 8368437,+            1454590 , 6389021,+            1686227 , 7519614,+            785355 , 3552602,+            190537 , 873605,+            968137 , 4495889,+            5458 , 25655,+            190537 , 905982,+            91313 , 438974,+            390321 , 1896172,+            416263 , 2042557,+            709397 , 3514492,+            1 , 5+          ]++-------------------------------------------------------------------------------+-- Unsafe+-------------------------------------------------------------------------------++#if CheckBounds+unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e+unsafeAt = (!)`+#endif++-------------------------------------------------------------------------------+-- Natural helpers+-------------------------------------------------------------------------------++fromNatural :: Num a => Natural -> a+fromNatural = fromIntegral++naturalLog2# :: Natural -> Int#+#if Base48+naturalLog2# (NatS# b) = wordLog2# b+naturalLog2# (NatJ# n) = integerLog2# (Jp# n)+#else+naturalLog2# n = integerLog2# (toInteger n)+#endif++#if __GLASGOW_HASKELL__ < 707+-- The times they are a-changing. The types of primops too :(+isTrue# :: Bool -> Bool+isTrue# = id+#endif
+ src/Math/NumberTheory/Powers/Integer.hs view
@@ -0,0 +1,62 @@+-- |+-- Module:      Math.NumberTheory.Powers.Integer+-- Copyright:   (c) 2011-2014 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Potentially faster power function for 'Integer' base and 'Int'+-- or 'Word' exponent.+--+{-# LANGUAGE CPP          #-}+{-# LANGUAGE MagicHash    #-}+{-# LANGUAGE BangPatterns #-}+module Math.NumberTheory.Powers.Integer+    ( integerPower+    , integerWordPower+    ) where++import GHC.Exts+import GHC.Integer.Logarithms.Compat (wordLog2#)++-- | Power of an 'Integer' by the left-to-right repeated squaring algorithm.+--   This needs two multiplications in each step while the right-to-left+--   algorithm needs only one multiplication for 0-bits, but here the+--   two factors always have approximately the same size, which on average+--   gains a bit when the result is large.+--+--   For small results, it is unlikely to be any faster than '(^)', quite+--   possibly slower (though the difference shouldn't be large), and for+--   exponents with few bits set, the same holds. But for exponents with+--   many bits set, the speedup can be significant.+--+--   /Warning:/ No check for the negativity of the exponent is performed,+--   a negative exponent is interpreted as a large positive exponent.+integerPower :: Integer -> Int -> Integer+integerPower b (I# e#) = power b (int2Word# e#)++-- | Same as 'integerPower', but for exponents of type 'Word'.+integerWordPower :: Integer -> Word -> Integer+integerWordPower b (W# w#) = power b w#++power :: Integer -> Word# -> Integer+power b w#+  | isTrue# (w# `eqWord#` 0##) = 1+  | isTrue# (w# `eqWord#` 1##) = b+  | otherwise           = go (wordLog2# w# -# 1#) b (b*b)+    where+      go 0# l h = if isTrue# ((w# `and#` 1##) `eqWord#` 0##) then l*l else (l*h)+      go i# l h+        | w# `hasBit#` i#   = go (i# -# 1#) (l*h) (h*h)+        | otherwise         = go (i# -# 1#) (l*l) (l*h)++-- | A raw version of testBit for 'Word#'.+hasBit# :: Word# -> Int# -> Bool+hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##)++#if __GLASGOW_HASKELL__ < 707+-- The times they are a-changing. The types of primops too :(+isTrue# :: Bool -> Bool+isTrue# = id+#endif
+ src/Math/NumberTheory/Powers/Natural.hs view
@@ -0,0 +1,63 @@+-- |+-- Module:      Math.NumberTheory.Powers.Natural+-- Copyright:   (c) 2011-2014 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Potentially faster power function for 'Natural' base and 'Int'+-- or 'Word' exponent.+--+{-# LANGUAGE CPP          #-}+{-# LANGUAGE MagicHash    #-}+{-# LANGUAGE BangPatterns #-}+module Math.NumberTheory.Powers.Natural+    ( naturalPower+    , naturalWordPower+    ) where++import GHC.Exts+import Numeric.Natural+import GHC.Integer.Logarithms.Compat (wordLog2#)++-- | Power of an 'Natural' by the left-to-right repeated squaring algorithm.+--   This needs two multiplications in each step while the right-to-left+--   algorithm needs only one multiplication for 0-bits, but here the+--   two factors always have approximately the same size, which on average+--   gains a bit when the result is large.+--+--   For small results, it is unlikely to be any faster than '(^)', quite+--   possibly slower (though the difference shouldn't be large), and for+--   exponents with few bits set, the same holds. But for exponents with+--   many bits set, the speedup can be significant.+--+--   /Warning:/ No check for the negativity of the exponent is performed,+--   a negative exponent is interpreted as a large positive exponent.+naturalPower :: Natural -> Int -> Natural+naturalPower b (I# e#) = power b (int2Word# e#)++-- | Same as 'naturalPower', but for exponents of type 'Word'.+naturalWordPower :: Natural -> Word -> Natural+naturalWordPower b (W# w#) = power b w#++power :: Natural -> Word# -> Natural+power b w#+  | isTrue# (w# `eqWord#` 0##) = 1+  | isTrue# (w# `eqWord#` 1##) = b+  | otherwise           = go (wordLog2# w# -# 1#) b (b*b)+    where+      go 0# l h = if isTrue# ((w# `and#` 1##) `eqWord#` 0##) then l*l else (l*h)+      go i# l h+        | w# `hasBit#` i#   = go (i# -# 1#) (l*h) (h*h)+        | otherwise         = go (i# -# 1#) (l*l) (l*h)++-- | A raw version of testBit for 'Word#'.+hasBit# :: Word# -> Int# -> Bool+hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##)++#if __GLASGOW_HASKELL__ < 707+-- The times they are a-changing. The types of primops too :(+isTrue# :: Bool -> Bool+isTrue# = id+#endif
+ test-suite/Math/NumberTheory/LogarithmsTests.hs view
@@ -0,0 +1,146 @@+-- |+-- Module:      Math.NumberTheory.LogarithmsTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Logarithms+--++{-# LANGUAGE CPP       #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.LogarithmsTests+  ( testSuite+  ) where++import Test.Tasty++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif+import Numeric.Natural++import Math.NumberTheory.Logarithms+import Math.NumberTheory.TestUtils++-- | Check that 'integerLogBase' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@.+integerLogBaseProperty :: Positive Integer -> Positive Integer -> Bool+integerLogBaseProperty (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n+  where+    l = toInteger $ integerLogBase b n++-- | Check that 'integerLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+integerLog2Property :: Positive Integer -> Bool+integerLog2Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n+  where+    l = toInteger $ integerLog2 n++integerLog2HugeProperty :: Huge (Positive Integer) -> Bool+integerLog2HugeProperty (Huge (Positive n)) = 2 ^ l <= n && 2 ^ (l + 1) > n+  where+    l = toInteger $ integerLog2 n++-- | Check that 'integerLog10' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@.+integerLog10Property :: Positive Integer -> Bool+integerLog10Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n+  where+    l = toInteger $ integerLog10 n++-- | Check that 'naturalLogBase' returns the largest natural @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@.+naturalLogBaseProperty :: Positive Natural -> Positive Natural -> Bool+naturalLogBaseProperty (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n+  where+    l = fromIntegral $ naturalLogBase b n++-- | Check that 'naturalLog2' returns the largest natural @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+naturalLog2Property :: Positive Natural -> Bool+naturalLog2Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n+  where+    l = fromIntegral $ naturalLog2 n++naturalLog2HugeProperty :: Huge (Positive Natural) -> Bool+naturalLog2HugeProperty (Huge (Positive n)) = 2 ^ l <= n && 2 ^ (l + 1) > n+  where+    l = fromIntegral $ naturalLog2 n++-- | Check that 'naturalLog10' returns the largest natural @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@.+naturalLog10Property :: Positive Natural -> Bool+naturalLog10Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n+  where+    l = fromIntegral $ naturalLog10 n++-- | Check that 'intLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+intLog2Property :: Positive Int -> Bool+intLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)+  where+    l = intLog2 n++-- | Check that 'wordLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+wordLog2Property :: Positive Word -> Bool+wordLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)+  where+    l = wordLog2 n++-- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@.+integerLogBase'Property :: Positive Integer -> Positive Integer -> Bool+integerLogBase'Property (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n+  where+    l = toInteger $ integerLogBase' b n++-- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@ for @b@ > 32 and @n@ >= @b@ ^ 2.+integerLogBase'Property2 :: Positive Integer -> Positive Integer -> Bool+integerLogBase'Property2 (Positive b') (Positive n') = b ^ l <= n && b ^ (l + 1) > n+  where+    b = b' + 32+    n = n' + b ^ 2 - 1+    l = toInteger $ integerLogBase' b n++-- | Check that 'integerLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+integerLog2'Property :: Positive Integer -> Bool+integerLog2'Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n+  where+    l = toInteger $ integerLog2' n++-- | Check that 'integerLog10'' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@.+integerLog10'Property :: Positive Integer -> Bool+integerLog10'Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n+  where+    l = toInteger $ integerLog10' n++-- | Check that 'intLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+intLog2'Property :: Positive Int -> Bool+intLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)+  where+    l = intLog2' n++-- | Check that 'wordLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@.+wordLog2'Property :: Positive Word -> Bool+wordLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2)+  where+    l = wordLog2' n++testSuite :: TestTree+testSuite = testGroup "Logarithms"+  [ testSmallAndQuick "integerLogBase"  integerLogBaseProperty+  , testSmallAndQuick "integerLog2"     integerLog2Property+  , testSmallAndQuick "integerLog2Huge" integerLog2HugeProperty+  , testSmallAndQuick "integerLog10"    integerLog10Property+  , testSmallAndQuick "naturalLogBase"  naturalLogBaseProperty+  , testSmallAndQuick "naturalLog2"     naturalLog2Property+  , testSmallAndQuick "naturalLog2Huge" naturalLog2HugeProperty+  , testSmallAndQuick "naturalLog10"    naturalLog10Property+  , testSmallAndQuick "intLog2"         intLog2Property+  , testSmallAndQuick "wordLog2"        wordLog2Property++  , testSmallAndQuick "integerLogBase'" integerLogBase'Property+  , testSmallAndQuick "integerLogBase' with base > 32 and n >= base ^ 2"+      integerLogBase'Property2+  , testSmallAndQuick "integerLog2'"    integerLog2'Property+  , testSmallAndQuick "integerLog10'"   integerLog10'Property+  , testSmallAndQuick "intLog2'"        intLog2'Property+  , testSmallAndQuick "wordLog2'"       wordLog2'Property+  ]
+ test-suite/Math/NumberTheory/Powers/IntegerTests.hs view
@@ -0,0 +1,41 @@+-- |+-- Module:      Math.NumberTheory.Powers.IntegerTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.Integer+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.IntegerTests+  ( testSuite+  ) where++import Test.Tasty++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.Powers.Integer+import Math.NumberTheory.TestUtils++-- | Check that 'integerPower' == '^'.+integerPowerProperty :: Integer -> Power Int -> Bool+integerPowerProperty a (Power b) = integerPower a b == a ^ b++-- | Check that 'integerWordPower' == '^'.+integerWordPowerProperty :: Integer -> Power Word -> Bool+integerWordPowerProperty a (Power b) = integerWordPower a b == a ^ b++testSuite :: TestTree+testSuite = testGroup "Integer"+  [ testSmallAndQuick "integerPower"     integerPowerProperty+  , testSmallAndQuick "integerWordPower" integerWordPowerProperty+  ]
+ test-suite/Math/NumberTheory/Powers/NaturalTests.hs view
@@ -0,0 +1,42 @@+-- |+-- Module:      Math.NumberTheory.Powers.NaturalTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.Natural+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.NaturalTests+  ( testSuite+  ) where++import Test.Tasty++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.Powers.Natural+import Numeric.Natural+import Math.NumberTheory.TestUtils++-- | Check that 'naturalPower' == '^'.+naturalPowerProperty :: Natural -> Power Int -> Bool+naturalPowerProperty a (Power b) = naturalPower a b == a ^ b++-- | Check that 'naturalWordPower' == '^'.+naturalWordPowerProperty :: Natural -> Power Word -> Bool+naturalWordPowerProperty a (Power b) = naturalWordPower a b == a ^ b++testSuite :: TestTree+testSuite = testGroup "Natural"+  [ testSmallAndQuick "naturalPower"     naturalPowerProperty+  , testSmallAndQuick "naturalWordPower" naturalWordPowerProperty+  ]
+ test-suite/Math/NumberTheory/TestUtils.hs view
@@ -0,0 +1,96 @@+-- |+-- Module:      Math.NumberTheory.TestUtils+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--++{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE StandaloneDeriving         #-}++{-# OPTIONS_GHC -fno-warn-orphans #-}++module Math.NumberTheory.TestUtils+  ( module Test.SmallCheck.Series+  , Power (..)+  , Huge (..)+  , testSmallAndQuick+  ) where++import Test.SmallCheck.Series (cons2)+import Test.Tasty+import Test.Tasty.SmallCheck as SC+import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative)+import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series, generate)++import Control.Applicative+import Data.Word+import Numeric.Natural++testSmallAndQuick+  :: SC.Testable IO a+  => QC.Testable a+  => String -> a -> TestTree+testSmallAndQuick name f = testGroup name+  [ SC.testProperty "smallcheck" f+  , QC.testProperty "quickcheck" f+  ]++-------------------------------------------------------------------------------+-- Serial monadic actions++instance Monad m => Serial m Word where+  series =+    generate (\d -> if d >= 0 then pure 0 else empty) <|> nats+    where+      nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty++instance Monad m => Serial m Natural where+  series =+    generate (\d -> if d >= 0 then pure 0 else empty) <|> nats+    where+      nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty++-------------------------------------------------------------------------------+-- Power++newtype Power a = Power { getPower :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real)++instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where+  series = Power <$> series `suchThatSerial` (> 0)++instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where+  arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0)+  shrink (Power x) = Power <$> filter (> 0) (shrink x)++suchThatSerial :: Series m a -> (a -> Bool) -> Series m a+suchThatSerial s p = s >>= \x -> if p x then pure x else empty++-------------------------------------------------------------------------------+-- Huge++newtype Huge a = Huge { getHuge :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real)++instance (Num a, Arbitrary a) => Arbitrary (Huge a) where+  arbitrary = do+    Positive l <- arbitrary+    ds <- vector (l :: Int)+    return $ Huge $ foldl1 (\acc n -> acc * 2^(63 :: Int) + n) ds++-- | maps 'Huge' constructor over series+instance Serial m a => Serial m (Huge a) where+  series = fmap Huge series++-------------------------------------------------------------------------------+-- Positive from smallcheck++instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where+  arbitrary = Positive <$> (arbitrary `suchThat` (> 0))+  shrink (Positive x) = Positive <$> filter (> 0) (shrink x)
+ test-suite/Test.hs view
@@ -0,0 +1,15 @@+import Test.Tasty++import qualified Math.NumberTheory.LogarithmsTests as Logarithms+import qualified Math.NumberTheory.Powers.IntegerTests as PowerInteger+import qualified Math.NumberTheory.Powers.NaturalTests as PowerNatural++main :: IO ()+main = defaultMain tests++tests :: TestTree+tests = testGroup "All"+    [ Logarithms.testSuite+    , PowerInteger.testSuite+    , PowerNatural.testSuite+    ]