fast-mult (empty) → 0.1.0.0

raw patch · 6 files changed

+366/−0 lines, 6 filesdep +basedep +ghc-primdep +integer-gmpsetup-changed

Dependencies added: base, ghc-prim, integer-gmp, strict-base

Files

LICENSE +30/−0
README.md +1/−0
Setup.hs +2/−0
fast-mult.cabal +26/−0
src/Data/FastMult.hs +3/−0
src/Data/FastMult/Internal.hs +304/−0

+ LICENSE view

@@ -0,0 +1,30 @@+Copyright Clinton Mead (c) 2017++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Clinton Mead nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

+ README.md view

@@ -0,0 +1,1 @@+# fast-mult

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

+ fast-mult.cabal view

@@ -0,0 +1,26 @@+name:                fast-mult+version:             0.1.0.0+synopsis:            Numeric type with asymptotically faster multiplications.+description:         This numeric type internally reorders multiplications to achieve+                     asymptotically faster multiplication of large numbers of small integers in particular.+                     See the module docs for more detail.+homepage:            https://github.com/clintonmead/fast-mult#readme+license:             BSD3+license-file:        LICENSE+author:              Clinton Mead+maintainer:          clintonmead@gmail.com+copyright:           Copyright: (c) 2017 Clinton Mead+category:            Web+build-type:          Simple+extra-source-files:  README.md+cabal-version:       >=1.10++library+  hs-source-dirs:      src+  exposed-modules:     Data.FastMult, Data.FastMult.Internal+  build-depends:       base >= 4.9 && < 5, integer-gmp, ghc-prim, strict-base+  default-language:    Haskell2010++source-repository head+  type:     git+  location: https://github.com/clintonmead/fast-mult

+ src/Data/FastMult.hs view

@@ -0,0 +1,3 @@+module Data.FastMult (FastMult, FastMultSeq, simplify) where++import Data.FastMult.Internal

+ src/Data/FastMult/Internal.hs view

@@ -0,0 +1,304 @@+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE UnboxedTuples #-}+{-# LANGUAGE BangPatterns #-}++module Data.FastMult.Internal (FastMult(FastMult), FastMultSeq, simplify) where++#include "MachDeps.h"++import Prelude hiding (Integer)++import GHC.Integer.GMP.Internals (BigNat, Integer(S#, Jp#, Jn#), sizeofBigNat#, timesBigNat, bigNatToWord, wordToBigNat, wordToBigNat2)+import Data.Bits (FiniteBits, countLeadingZeros, finiteBitSize, xor, complement)+import GHC.Base (timesWord2#, Word(W#), int2Word#, Int(I#), eqWord#, (>=#), negateInt#)+import Data.Ord (comparing)+import Data.Strict.List ( List((:!)) )+import qualified Data.Strict.List as Strict+import GHC.TypeLits (Nat, KnownNat, natVal)+import GHC.Conc.Sync (par)+import Data.Proxy (Proxy)+import Data.Word (Word)+import Data.Ratio ((%))+import Data.Foldable (foldl')++-- We use 'Word' here not 'Bool' because it unpacks. I'm not sure if this is an optimisation.+newtype Sign = Sign Word deriving Show++pattern Pos = Sign 0+pattern Neg = Sign (-1)++multSigns :: Sign -> Sign -> Sign+multSigns (Sign x) (Sign y) = Sign (x `xor` y)++negateSign :: Sign -> Sign+negateSign (Sign x) = Sign (complement x)+++data BigNatWithScale (n :: Nat) where+  BigNatWithScale :: KnownNat n => {-# UNPACK #-} !Word -> BigNat -> BigNatWithScale n++-- Internal debug only:+instance Show (BigNatWithScale n) where+  show (BigNatWithScale scale bigNat) = "(BigNatWithScale scale = " ++ show scale ++ ", num = " ++ show (Jp# bigNat) ++ ")"++getBigNat :: BigNatWithScale (n :: Nat) -> BigNat+getBigNat (BigNatWithScale _ x) = x++{-|+  'FastMult' is a Numeric type that can be used in any place a 'Num a' is required.+  It represents a standard integer using three components, which multiplied together represent the stored number:++  1. The number's sign+  2. An unsigned machine word.+  3. A (possibly empty) list of 'BigNat's, which are the internal type for 'Integer's which are too large to fit in a machine word.++  Each 'BigNat' in the list has a scale. It's scale is the log base 2 of the number of words to store the machine word, minus 1.++  Note that we never store BigNats with length of only one machine word in this list, we instead convert them to an ordinary+  unsigned machine word and multiply them by item 2 in the list above. Only then if the result overflows we place them in this+  'BigNat' list.++  This is a few examples of "MachineWords -> Scale"++  2 -> 0+  3 -> 1+  4 -> 1+  5 -> 2+  6..8 -> 2+  9..16 -> 3+  17..32 -> 4++  etc.++  Note this "scale" has the very nice property that multipling 'BigNat's of scale @x@ always results in a 'BigNat' of scale @x+1@.++  The list of 'BigNat's only ever contains one 'BigNat' of each "scale". As the size of 'BigNat's increases exponentially with scale,+  this list should always be relatively small. The 'BigNat' list is always sorted as well, smallest to largest.++  When we multiply two 'FastMult's, we merge the BigNat lists. This is basically a simple merge of sorted list,+  but with one significant change. Note that we said that the 'BigNat' list cannot contain two 'BigNat's of the same scale.+  So if find that a 'BigNat' in the left hand list of the multiplication is the same scale as a 'BigNat' in right hand list,+  we multiply these two 'BigNat's to create a 'BigNat' one "scale" larger. We then continue the merge, including this new BigNat.++  As a result, we only ever multiply numbers of the same "scale", that is, no more than double the length of one another.++  Why do we do this? Well, an ordinary product, say @product [1..1000000]@, towards the end of the list involves multiplications+  of a very large number by a machine word. These take @O(n)@ time. So the whole product takes @O(n^2)@ time.++  If we instead did the following:++  @+    product x y = product x mid * product mid y+      mid = (x + y) `div` 2++    (suitible base case here)+  @++  We find that this runs a lot faster. The reason is that with this approach we're minimising products involving very large numbers,+  and importantly, multiplying two @n@ length numbers doesn't take @O(n^2)@ but more like @O(n*log(n))@ time.+  For this reason it's better to do a few multiplication of large numbers by large numbers,+  instead of lots of multiplications of large numbers by small numbers.++  But to do this I've had to redefine product. What if you don't want to change the algorithm, but just want to use one that's+  already been written, perhaps inefficiently. Well this is where 'FastMult' is useful. Instead of making the algorithm smarter,+  'FastMult' just makes numbers smarter. The numbers themselves reorder the multiplications so you don't have too.++  As well as having the advantage of speeding up existing algorithms, 'FastMult' dynamically behaves differently based on+  what numbers it's actually multiplying and always maintains the invariant that multiplications will not be performed between+  numbers greater than twice the size each other.++  At this point I haven't mentioned the meaning of the `FastMult` type parameter @n@'. 'FastMult' can also add paralellism to+  your multiplication algorithms. However, sparking new GHC threads has a cost, so we only want to do it for large multiplications.+  Multiplications of @scale > n@ will spark a new thread, so @n = 0@ will spark new threads for any multiplication+  involving at least 3 machine words. This is probably too small, you can experiment with different numbers.+  Note that @n@ represents the scale, not size, so for example setting @n=4@ will only spark threads for multiplications involving+  at least 33 machine words.++  How well parallelism works (or if it works at all) hasn't been tested yet however.++  We include an ordinary machine word in the type as an optimisation for single machine word numbers.+  This is because multiplying 'BigNat's involves calling GMP using a C call, which is a large overhead for small multiplications.++  To use 'FastMult', all you have to do is import it's type, not it's implementation.+  If you're not interested in parallelism, just import 'FastMultSeq'.++  For example, just compare in GHCi:++  @+  product [1..100000]+  @++  and:++  @+  product [1::FastMultSeq..100000]+  @++  and you should find the latter completes much faster.++  Converting to and from 'Integer's can be done with the+  'toInteger' and 'fromInteger' class methods from 'Integral' and 'Num' respectively.+-}+data FastMult (n :: Nat) where+  FastMult :: KnownNat n => {-# UNPACK #-} !Sign -> {-# UNPACK #-} !Word -> !(Strict.List (BigNatWithScale n)) -> FastMult n++{-|+  A type synonym for a fully sequential 'FastMult'. The parameter is supposed to be 'WORD_MAX', but I couldn't find that+  defined, anyway what's important is that anything of scale smaller than @0xFFFFFFFF@ will be sequential, which is everything.+-}+type FastMultSeq = FastMult 0xFFFFFFFF++data BigNatMultResult (n :: Nat) where+  ScaleLT :: BigNatMultResult n+  ScaleEQ :: (KnownNat n) => BigNatWithScale n -> BigNatMultResult n+  ScaleGT :: BigNatMultResult n++singletonStrictList :: a -> Strict.List a+singletonStrictList x = x :! Strict.Nil++instance KnownNat n => Eq (FastMult n) where+  x == y = toInteger x == toInteger y++instance KnownNat n => Ord (FastMult n) where+  x `compare` y = toInteger x `compare` toInteger y++instance KnownNat n => Enum (FastMult n) where+  toEnum = fromIntegral+  fromEnum = fromIntegral++instance KnownNat n => Num (FastMult n) where+  fromInteger = \case+    (S# prim_i) -> case (prim_i >=# 0#) of+      1# -> FastMult Pos (W# (int2Word# prim_i)) Strict.Nil+      0# -> FastMult Neg (W# (int2Word# (negateInt# prim_i))) Strict.Nil+    (Jp# x) -> fromBigNat Pos x+    (Jn# x) -> fromBigNat Neg x+    where+      fromBigNat :: Sign -> BigNat -> FastMult n+      fromBigNat sign x = case (W# (int2Word# (sizeofBigNat# x)) - 1) of+        0 -> FastMult sign (W# (bigNatToWord x)) Strict.Nil+        size -> FastMult sign 1 (singletonStrictList (BigNatWithScale (logBase2Int size) x))+      logBase2Int :: Word -> Word+      logBase2Int x = WORD_SIZE_IN_BITS - 1 - (fromIntegral (countLeadingZeros x))+  (FastMult s1 w1 l1) * (FastMult s2 w2 l2) =+    let+      multBigNatWithScale :: forall n. BigNatWithScale n -> BigNatWithScale n -> BigNatMultResult n+      multBigNatWithScale (BigNatWithScale s1 n1) (BigNatWithScale s2 n2) =+          case (s1 `compare` s2) of+            EQ -> result `seqOrPar` (ScaleEQ (BigNatWithScale (s1 + 1) result)) where+              result = n1 `timesBigNat` n2+              seqOrPar = if s1 <= maxSeq then seq else par+            LT -> ScaleLT+            GT -> ScaleGT+          where+            maxSeq = fromIntegral (natVal (undefined :: Proxy n))+      sr = multSigns s1 s2+      (# wu_prim, wl_prim #) =+        let+          !(W# w1_prim) = w1+          !(W# w2_prim) = w2+        in+          timesWord2# w1_prim w2_prim+      merge :: Strict.List (BigNatWithScale n) -> Strict.List (BigNatWithScale n) -> Strict.List (BigNatWithScale n)+      merge xl Strict.Nil = xl+      merge Strict.Nil yl = yl+      merge xl@(x:!xs) yl@(y:!ys) = case multBigNatWithScale x y of+        ScaleEQ result -> mergeWithCarry result xs ys+        ScaleLT -> x :! merge xs yl+        ScaleGT -> y :! merge xl ys++      mergeWithCarry :: BigNatWithScale n -> Strict.List (BigNatWithScale n) -> Strict.List (BigNatWithScale n) -> Strict.List (BigNatWithScale n)+      mergeWithCarry carry xl Strict.Nil = mergeOneCarry carry xl+      mergeWithCarry carry Strict.Nil yl = mergeOneCarry carry yl+      mergeWithCarry carry xl@(x:!xs) yl@(y:!ys) = case multBigNatWithScale x y of+        ScaleEQ result -> mergeWithCarry carry xs ys+        ScaleLT -> contCarry x xs yl+        ScaleGT -> contCarry y ys xl+        where+          contCarry x xs yl =+            case multBigNatWithScale carry x of+              ScaleEQ result -> mergeWithCarry result xs yl+              ScaleLT -> carry :! x :! merge xs yl+              ScaleGT -> error $+                "Carry should never be larger than first element. This should never happen. Report as bug.\n" +++                "Details:\n" +++                "carry =\n" +++                show carry ++ "\n" +++                "xl =\n" +++                show xl ++ "\n" +++                "yl =\n" +++                show yl ++ "\n"++      mergeOneCarry carry Strict.Nil = singletonStrictList carry+      mergeOneCarry carry xl@(x:!xs) = case multBigNatWithScale carry x of+        ScaleLT -> carry :! xl+        ScaleEQ result -> mergeOneCarry result xs+        ScaleGT -> error $+          "Carry should never be larger than first element. This should never happen. Report as bug.\n" +++          "Details:\n" +++          "carry =\n" +++          show carry ++ "\n" +++          "xl =\n" +++          show xl ++ "\n"+    in+      case eqWord# wu_prim (int2Word# 0#) of+        0# -> FastMult sr (W# wl_prim) (merge l1 l2)+        _  -> FastMult sr 1 (mergeWithCarry (BigNatWithScale 0 (wordToBigNat2 wu_prim wl_prim)) l1 l2)+  (+) = binaryViaInteger (+)+  (-) = binaryViaInteger (-)+  abs (FastMult _ word l) = FastMult Pos word l+  signum (FastMult Pos _ _) = FastMult Pos 1 Strict.Nil+  signum (FastMult Neg _ _) = FastMult Neg 1 Strict.Nil+  negate (FastMult sign word l) = FastMult (negateSign sign) word l++binaryViaInteger f x y = fromInteger (toInteger x `f` toInteger y)+unaryViaInteger f = fromInteger . f . toInteger++instance KnownNat n => Real (FastMult n) where+  toRational x = (toInteger x) % 1++instance KnownNat n => Integral (FastMult n) where+  toInteger (FastMult sign (W# word_prim) l) = case sign of+    Pos -> Jp# result+    Neg -> Jn# result+    where+      result = foldl' (\x y -> x `timesBigNat` getBigNat y) (wordToBigNat word_prim) l+  x `quotRem` y = let (x_r, y_r) = (toInteger x `quotRem` toInteger y) in (fromInteger x_r, fromInteger y_r)++instance KnownNat n => Show (FastMult n) where+  show = show . toInteger++instance KnownNat n => Read (FastMult n) where+  readsPrec p s = map (\(x,y) -> (fromInteger x,y)) (readsPrec p s)++{-|+  'simplify' returns a 'FastMult' the same as it's argument but "simplified".++  To explain this, consider the following for @x :: FastMult@:++  @+  f x = (show x, x + 1)+  @++  It will multiply out @x@ twice, once for the addition, and once for 'show'. Note that the list of 'BigInt's in @x@ is generally+  a small number, as only one 'BigInt' is stored for each scale, and the sizes of scales increase exponentially, but there+  may be some multiplications required nevertheless. A better way to write this is as follows:++  @+  f x = let y = simplify x in (show y, y + 1)+  @++  This will ensure that @x@ is multiplied out only once.++  Unfortunately using 'simplify' stops your algorithms from being generic,+  so it might be better to define simplify as 'id' with a rewrite rule. I'll think about this.+-}+simplify :: KnownNat n => FastMult n -> FastMult n+simplify = fromInteger . toInteger