diff --git a/Numeric/FastMath.hs b/Numeric/FastMath.hs --- a/Numeric/FastMath.hs +++ b/Numeric/FastMath.hs @@ -1,25 +1,18 @@ --- | Compile-time optimisations for 'Float' and 'Double' that break IEEE-754 --- compatibility. --- --- Namely, this otherwise empty module contains RULES that rewrite @x-x@, --- @x*0@ and @0*x@ to @0@, which is incorrect (according to IEEE-754) when --- @x@ is @NaN@. +-- | This module loads all rewrite rules. Unless you know that some rules +-- will be unsafe for your application, this is the module you should load. +-- Importing this module is roughly equivalent to gcc's @-ffast-math@ +-- compilation flag. -- --- At the time of writing, @base-4.3.1.0:GHC/Base.lhs@ erroneously includes --- these rules for 'Float's, but not for 'Double's. This has been reported --- as GHC bug #5178: . - -module Numeric.FastMath () where - -import GHC.Exts - -{-# RULES -"minusFloat x x" forall x#. minusFloat# x# x# = 0.0# -"timesFloat x 0.0" forall x#. timesFloat# x# 0.0# = 0.0# -"timesFloat 0.0 x" forall x#. timesFloat# 0.0# x# = 0.0# +-- The best way to figure out what optimizations these modules do is by +-- looking at the source code. RULES pragmas are surprisingly readable. -"minusDouble x x" forall x#. (-##) x# x# = 0.0## -"timesDouble 0.0 x" forall x#. (*##) 0.0## x# = 0.0## -"timesDouble x 0.0" forall x#. (*##) x# 0.0## = 0.0## - #-} +module Numeric.FastMath + ( module Numeric.FastMath.Approximation + , module Numeric.FastMath.NaN + , module Numeric.FastMath.SignedZeros + ) + where +import Numeric.FastMath.Approximation () +import Numeric.FastMath.NaN () +import Numeric.FastMath.SignedZeros () diff --git a/Numeric/FastMath/Approximation.hs b/Numeric/FastMath/Approximation.hs new file mode 100644 --- /dev/null +++ b/Numeric/FastMath/Approximation.hs @@ -0,0 +1,256 @@ +-- | This module contains rewrite rules that may change the lowest order bits +-- of a computation. They take advantage of: +-- +-- * distributivity +-- +-- * repeated addition/multiplication +-- +-- * exponentiation rules +-- +-- All of these RULES should be safe in the presence of `NaN` and `Infinity` +-- +-- Importing this module is similar to compiling with gcc's +-- @-funsafe-math-operations@. +-- +module Numeric.FastMath.Approximation + where + +import GHC.Exts +import Prelude + +--------------------------------------- +-- distributivity +-- +-- NOTE: these rules are sufficient to capture the property +-- +-- > x*y1+x*y2+x*y3 == x*(y1+y2+y3) +-- +-- because they will be applied recursively during the optimization passes + +{-# RULES + +"double *,+ distribute A" forall x y1 y2. (x *## y1) +## (x *## y2) + = x *## (y1 +## y2) + +"double *,+ distribute B" forall x y1 y2. (y1 *## x) +## (x *## y2) + = x *## (y1 +## y2) + +"double *,+ distribute C" forall x y1 y2. (y1 *## x) +## (y2 *## x) + = x *## (y1 +## y2) + +"double *,+ distribute D" forall x y1 y2. (x *## y1) +## (y2 *## x) + = x *## (y1 +## y2) + + + +"double *,- distribute A" forall x y1 y2. (x *## y1) -## (x *## y2) + = x *## (y1 -## y2) + +"double *,- distribute B" forall x y1 y2. (y1 *## x) -## (x *## y2) + = x *## (y1 -## y2) + +"double *,- distribute C" forall x y1 y2. (y1 *## x) -## (y2 *## x) + = x *## (y1 -## y2) + +"double *,- distribute D" forall x y1 y2. (x *## y1) -## (y2 *## x) + = x *## (y1 -## y2) + + + +"double /,+ distribute" forall x y1 y2. (y1 *## x) +## (y2 *## x) + = (y1 +## y2) /## x + +"double /,- distribute" forall x y1 y2. (y1 /## x) -## (y2 /## x) + = (y1 -## y2) /## x + + #-} + + + +{-# RULES + +"float *,+ distribute A" forall x y1 y2. (x `timesFloat#` y1) `plusFloat#` (x `timesFloat#` y2) + = x `timesFloat#` (y1 `plusFloat#` y2) + +"float *,+ distribute B" forall x y1 y2. (y1 `timesFloat#` x) `plusFloat#` (x `timesFloat#` y2) + = x `timesFloat#` (y1 `plusFloat#` y2) + +"float *,+ distribute C" forall x y1 y2. (y1 `timesFloat#` x) `plusFloat#` (y2 `timesFloat#` x) + = x `timesFloat#` (y1 `plusFloat#` y2) + +"float *,+ distribute D" forall x y1 y2. (x `timesFloat#` y1) `plusFloat#` (y2 `timesFloat#` x) + = x `timesFloat#` (y1 `plusFloat#` y2) + + + +"float *,- distribute A" forall x y1 y2. (x `timesFloat#` y1) `minusFloat#` (x `timesFloat#` y2) + = x `timesFloat#` (y1 `minusFloat#` y2) + +"float *,- distribute B" forall x y1 y2. (y1 `timesFloat#` x) `minusFloat#` (x `timesFloat#` y2) + = x `timesFloat#` (y1 `minusFloat#` y2) + +"float *,- distribute C" forall x y1 y2. (y1 `timesFloat#` x) `minusFloat#` (y2 `timesFloat#` x) + = x `timesFloat#` (y1 `minusFloat#` y2) + +"float *,- distribute D" forall x y1 y2. (x `timesFloat#` y1) `minusFloat#` (y2 `timesFloat#` x) + = x `timesFloat#` (y1 `minusFloat#` y2) + + + +"float /,+ distribute" forall x y1 y2. (y1 `timesFloat#` x) `plusFloat#` (y2 `timesFloat#` x) + = (y1 `plusFloat#` y2) `divideFloat#` x + +"float /,- distribute" forall x y1 y2. (y1 `divideFloat#` x) `minusFloat#` (y2 `divideFloat#` x) + = (y1 `minusFloat#` y2) `divideFloat#` x + + #-} + +--------------------------------------- +-- fancy distributing +-- +-- NOTE: I'm not yet sure if all of these are a great idea to have on by +-- default due to stability issues... + +{-# RULES + +"double **,* distribute" forall x y1 y2. (y1 **## x) *## (y2 **## x) = (y1 *## y2) *## x + +"double **,log distribute" forall x y. logDouble# (x **## y) = y *## (logDouble# x) + + #-} + +--------------------------------------- +-- Repeated addition +-- +-- NOTE: It is important that these rules should fire after the distributivity +-- rules. This ensures that +-- +-- > x*x+x*y +-- +-- gets simplified to +-- +-- > x*(x+y) +-- +-- rather than +-- +-- > x+x+x*y +-- +{-# RULES + +"double mulToAdd 2" [0] forall x . x *## 2.0## = x +## x +"double mulToAdd 3" [0] forall x . x *## 3.0## = x +## x +## x +"double mulToAdd 4" [0] forall x . x *## 4.0## = x +## x +## x +## x + + #-} + +{-# RULES + +"float mulToAdd 2" [0] forall x . timesFloat# x 2.0# = plusFloat# x x +"float mulToAdd 3" [0] forall x . timesFloat# x 3.0# = plusFloat# x (plusFloat# x x) +"float mulToAdd 4" [0] forall x . timesFloat# x 4.0# = plusFloat# x (plusFloat# x (plusFloat# x x)) + + #-} + +--------------------------------------- +-- left associate / commute + +-- NOTE: phase controls are needed to prevent infinite loops when interacting +-- with the repeated multiplication rules. +-- +-- We should slightly prefer commuting rather than associating because it doesn't +-- change the floating point results + +{-# RULES + +"double commute left *" [~2] forall x1 x2 x3. (*##) x1 ((*##) x2 x3) = (*##) ((*##) x2 x3) x1 +"double associate left *" [~1] forall x1 x2 x3. (*##) x1 ((*##) x2 x3) = (*##) ((*##) x1 x2) x3 + +"double commute left +" [~2] forall x1 x2 x3. (+##) x1 ((+##) x2 x3) = (+##) ((+##) x2 x3) x1 +"double associate left +" [~1] forall x1 x2 x3. (+##) x1 ((+##) x2 x3) = (+##) ((+##) x1 x2) x3 + + #-} + +{-# RULES + +"float commute left *" [~2] forall x1 x2 x3. timesFloat# x1 (timesFloat# x2 x3) = timesFloat# (timesFloat# x2 x3) x1 +"float associate left *" [~1] forall x1 x2 x3. timesFloat# x1 (timesFloat# x2 x3) = timesFloat# (timesFloat# x1 x2) x3 + +"float commute left +" [~2] forall x1 x2 x3. plusFloat# x1 (plusFloat# x2 x3) = plusFloat# (plusFloat# x2 x3) x1 +"float associate left +" [~1] forall x1 x2 x3. plusFloat# x1 (plusFloat# x2 x3) = plusFloat# (plusFloat# x1 x2) x3 + + #-} + +--------------------------------------- +-- Repeated multiplication + +-- FIXME: I can't get thise rules to work for more than 4 repeats without +-- causing an infinite loop in the simplifier + +{-# RULES + +"double repmul 4" [1] forall x . ((x *## x) *## x) *## x + = let xx = (x *## x) in (xx *## xx) + + #-} + +-- "double repmul 5" forall x . x *## x *## x *## x *## x +-- = let xx = x *## x in xx *## xx *## x +-- +-- "double repmul 6" forall x . x *## x *## x *## x *## x *## x +-- = let xx = x *## x in xx *## xx *## xx +-- +-- "double repmul 7" forall x . x *## x *## x *## x *## x *## x *## x +-- = let xx = x *## x in xx *## xx *## xx *## x +-- +-- "double repmul 8" forall x . x *## x *## x *## x *## x *## x *## x *## x +-- = let xxx = (let xx = x *## x in xx *## xx) in xxx *## xxx + +{-# RULES + +"double repmul 4" forall x . timesFloat# x (timesFloat# x (timesFloat# x x)) + = let xx = timesFloat# x x in timesFloat# xx xx + + #-} + +-- "double repmul 5" forall x . timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x x))) +-- = let xx = timesFloat# x x in timesFloat# x (timesFloat# xx xx) +-- +-- "double repmul 6" forall x . timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x x)))) +-- = let xx = timesFloat# x x in timesFloat# xx (timesFloat# xx xx) +-- +-- "double repmul 7" forall x . timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x x))))) +-- = let xx = timesFloat# x x in timesFloat# x (timesFloat# xx (timesFloat# xx xx)) +-- +-- "double repmul 8" forall x . timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x (timesFloat# x x)))))) +-- = let xxx = (let xx = timesFloat# x x in timesFloat# xx xx) in timesFloat# xxx xxx + + +--------------------------------------- +-- Exponentiation + +{-# RULES +"double **0" forall x . x **## 0.0## = 1.0## +"double **1" forall x . x **## 1.0## = x +"double **2" forall x . x **## 2.0## = x *## x +"double **3" forall x . x **## 3.0## = x *## x *## x +"double **4" forall x . x **## 4.0## = let xx = x *## x in xx *## xx +"double **8" forall x . x **## 8.0## = let xxx = (let xx = x *## x in xx *## xx) in xxx *## xxx + +"double **(1/2)" forall x## . x## **## 0.500## = sqrtDouble# x## +"double **(1/4)" forall x## . x## **## 0.250## = sqrtDouble# (sqrtDouble# x##) +"double **(1/8)" forall x## . x## **## 0.125## = sqrtDouble# (sqrtDouble# (sqrtDouble# x##)) + #-} + +{-# RULES +"float **0" forall x# . powerFloat# x# 0.0# = 1.0# +"float **1" forall x# . powerFloat# x# 1.0# = x# +"float **2" forall x# . powerFloat# x# 2.0# = timesFloat# x# x# +"float **3" forall x# . powerFloat# x# 3.0# = timesFloat# (timesFloat# x# x#) x# +"float **4" forall x# . powerFloat# x# 4.0# = let xx# = (timesFloat# x# x#) in timesFloat# xx# xx# +"float **8" forall x# . powerFloat# x# 8.0# = let xxx# = (let xx# = (timesFloat# x# x#) in timesFloat# xx# xx#) in timesFloat# xxx# xxx# + +"float **(1/2)" forall x# . powerFloat# x# 0.500# = sqrtFloat# x# +"float **(1/4)" forall x# . powerFloat# x# 0.250# = sqrtFloat# (sqrtFloat# x#) +"float **(1/8)" forall x# . powerFloat# x# 0.125# = sqrtFloat# (sqrtFloat# (sqrtFloat# x#)) + #-} + diff --git a/Numeric/FastMath/Infinitesimal.hs b/Numeric/FastMath/Infinitesimal.hs deleted file mode 100644 --- a/Numeric/FastMath/Infinitesimal.hs +++ /dev/null @@ -1,11 +0,0 @@ --- | Also rewrite @0/x@ to @+0@, which should really be @-0@ for negative @x@. - -module Numeric.FastMath.Infinitesimal () where - -import GHC.Exts - -{-# RULES -"divideFloat 0.0 x" forall x#. divideFloat# 0.0# x# = 0.0# -"divideDouble 0.0 x" forall x#. (/##) 0.0## x# = 0.0## - #-} - diff --git a/Numeric/FastMath/NaN.hs b/Numeric/FastMath/NaN.hs new file mode 100644 --- /dev/null +++ b/Numeric/FastMath/NaN.hs @@ -0,0 +1,35 @@ +-- | This module contains rules that break the way NaN is handled for "Float" +-- and "Double" types. Still, these rules should be safe in the vast majority of +-- applications. +-- +-- Importing this module is similar to compiling with gcc's @-fno-signaling-nans@ +-- and @-ffinite-math-only@. +-- +module Numeric.FastMath.NaN + where + +import GHC.Exts + +{-# RULES +"minusDouble x x" forall x. (-##) x x = 0.0## + +"timesDouble 0 x" forall x. (*##) 0.0## x = 0.0## +"timesDouble x 0" forall x. (*##) x 0.0## = 0.0## + +"divideDouble x 1" forall x. (/##) x 1.0## = x +"divideDouble x -1" forall x. (/##) x (-1.0##) = negateDouble# x +"divideDouble 0 x" forall x. (/##) 0.0## x = 0.0## + #-} + +{-# RULES +"minusFloat x x" forall x. minusFloat# x x = 0.0# + +"timesFloat x 0" forall x. timesFloat# x 0.0# = 0.0# +"timesFloat 0 x" forall x. timesFloat# 0.0# x = 0.0# + +"divideFloat x 1" forall x. divideFloat# x 1.0# = x +"divideFloat x -1" forall x. divideFloat# x (-1.0#) = negateFloat# x +"divideFloat 0 x" forall x. divideFloat# 0.0# x = 0.0# + + #-} + diff --git a/Numeric/FastMath/SignedZeros.hs b/Numeric/FastMath/SignedZeros.hs new file mode 100644 --- /dev/null +++ b/Numeric/FastMath/SignedZeros.hs @@ -0,0 +1,24 @@ +-- | IEEE 754 math makes a distrinction between -0.0 and +0.0. This module +-- contains RULES that ignore this distinction. +-- +-- Importing this module is similar to compiling with gcc's +-- @-fno-signed-zeros@. + +module Numeric.FastMath.SignedZeros () where + +import GHC.Exts + +{-# RULES + +"minusDouble 0 x" forall x. (-##) 0.0## x = negateDouble# x +"divideDouble 0 x" forall x. (/##) 0.0## x = 0.0## + #-} + +{-# RULES + +"minusFloat 0 x" forall x. minusFloat# 0.0# x = negateFloat# x +"divideFloat 0 x" forall x. divideFloat# 0.0# x = 0.0# + + #-} + + diff --git a/fast-math.cabal b/fast-math.cabal --- a/fast-math.cabal +++ b/fast-math.cabal @@ -1,34 +1,42 @@ name: fast-math -version: 0.1 +version: 1.0 synopsis: Non IEEE-754 compliant compile-time floating-point optimisations description: - The "Numeric.FastMath" module brings into scope RULES for 'Float's and - 'Double's that rewrite @x-x@, @0*x@ and @x*0@ to @0@. This disagrees - with IEEE-754 when @x@ is @NaN@, but is acceptable for most - applications. - . - Importing "Numeric.FastMath.Infinitesimal" also rewrites @0/x@ to @0@. + The "Numeric.FastMath" module brings into scope many unsafe @RULES@ for + 'Float's and 'Double's that can greatly improve run time performance. + It is roughly equivalent to gcc\'s @-ffast-math@ compiler flag. + Optimisation (at least @-O1@) must be enabled for any @RULES@ to take effect. . - Optimisation (at least @-O1@) must be enabled for any RULES to take effect. + These rules are unsafe because they don't strictly adhere to the + IEEE-754 regulations and may subtly change the results of your numeric computations. + See the on github for more details. + license: BSD3 license-file: LICENSE -author: Liyang HU +author: Liyang HU and Mike Izbicki maintainer: fast-math@liyang.hu copyright: © 2011, Liyang HU category: Math, Numeric build-type: Simple -cabal-version: >= 1.2.3 +cabal-version: >= 1.10 +source-repository head + type: git + location: http://github.com/liyang/fast-math + library + default-language: Haskell2010 build-depends: base >= 3 && < 5 exposed-modules: Numeric.FastMath - Numeric.FastMath.Infinitesimal - extensions: + Numeric.FastMath.Approximation + Numeric.FastMath.NaN + Numeric.FastMath.SignedZeros + default-extensions: NoImplicitPrelude MagicHash - ghc-options: -Wall -fno-warn-orphans - --- vim: et ts=4 sw=4: + ghc-options: + -Wall + -fno-warn-orphans