posit 2022.1.0.0 → 2022.2.0.0
raw patch · 9 files changed
+285/−145 lines, 9 filesdep ~posit
Dependency ranges changed: posit
Files
- ChangeLog.md +14/−0
- README.md +9/−6
- posit.cabal +28/−5
- src/Posit.hs +75/−58
- src/Posit/Internal/PositC.hs +18/−2
- stack.yaml +3/−2
- test/Test/Algorithms.hs +36/−18
- test/TestFusedRewrite.hs +48/−0
- test/TestPosit.hs +54/−54
ChangeLog.md view
@@ -1,5 +1,19 @@ # Changelog for Posit Numbers +# posit-2022.2.0.0++ * Updated divide (/) to have one less round operation, previous default was `a * recip b`, divide two rational numbers then round `viaRational2 (/) a b`, it sort of fuses the old default, reducing error+ * New flag for Fused Operations applied automatically via rewrite rules `--flag posit:do-rewrite`+ * Test for Fused Operations `stack test posit:test-posit-fusedRewrite`, with no fusing+ * Test for Fused Operations `stack test posit:test-posit-fusedRewrite --flag posit:do-rewrite`, to execute with the rewrite rules for fusing+ * Added several `fms` or `fma` functions to `asin`, `acos`, `asinh`, `acosh`. No apparent change in the elementary functions.+ * User defined rewrite rules seem to work before but break in `ghc-9.6.3`. I really hope it gets fixed. See: (How Cool is Rulz!)[https://discourse.haskell.org/t/how-cool-is-rulz/7738/9]+ * Added new `approxEq` function to the `AltFloating` class+ * Implemented `approxEq` by converting to the next lower `ES` and then using the normal `(==)` function, using a new `Prev es` type family+ * Added new `Prev es` type family to the `PositF es` Constraint Synonym+ * Changed naming in `AltFloating` to free up the name space and not use greek letters, `eps` is now `machEps`, `psi` is now `goldenRatio`+ * Removed `gamma`, `sinc`, `expm1` from `AltFloating` making way for updates to the Posit Standard in the future for special functions+ # posit-2022.1.0.0 * Optimized the `Show` instance to properly handle the tapered precision
README.md view
@@ -1,4 +1,4 @@-# posit 2022.1.0.0+# posit 2022.2.0.0 The [Posit Standard 2022](https://posithub.org/docs/posit_standard-2.pdf), and [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf), @@ -51,15 +51,18 @@ ``` class AltFloating p where- eps :: p -- Machine Epsilon near 1.0- phi :: p- gamma :: p -> p- sinc :: p -> p- expm1 :: p -> p+ machEps :: p -- Machine Epsilon near 1.0+ approxEq :: p -> p -> Bool+ goldenRatio :: p hypot2 :: p -> p -> p hypot3 :: p -> p -> p -> p hypot4 :: p -> p -> p -> p -> p ```++The 'FusedOps' class has been updated with rewrite rules, enabled by the+`--flag posit:do-rewrite` flag. When the flag is enabled, the `ghc`+compiler will match the pattern and then automatically rewrite with the+fused operation. ``` class Num a => FusedOps a where
posit.cabal view
@@ -1,7 +1,7 @@ cabal-version: 1.12 name: posit-version: 2022.1.0.0+version: 2022.2.0.0 category: Numeric, Math description: The Posit Number format attempting to conform to the Posit Standard Versions 3.2 and 2022. Where Real numbers are approximated by `Maybe Rational` and sampled in a similar way to the projective real line. homepage: https://github.com/waivio/posit#readme@@ -12,12 +12,15 @@ license: BSD3 license-file: LICENSE build-type: Simple-tested-with: GHC == 8.10.4,- GHC == 8.10.7,+tested-with: GHC == 8.10.7, GHC == 9.0.2, GHC == 9.2.8,- GHC == 9.4.6,- GHC == 9.6.2+ GHC == 9.4.8,+ GHC == 9.6.3+ -- fails rewrite rules with ghc-9.6.3 bummer+ -- Things before ghc-8.10.7 fail, don't care! ++ synopsis: Posit Numbers extra-source-files: README.md@@ -38,6 +41,10 @@ manual: True default: False +flag do-rewrite+ description: Build with Rewrite Rules for Fused Operations.+ manual: True+ default: False library exposed-modules:@@ -53,6 +60,10 @@ -- Compiler options ghc-options: -Wall -O2 + if flag(do-rewrite)+ cpp-options: -DO_REWRITE+ ghc-options: -fforce-recomp -ddump-rule-firings+ if flag(do-liquid) ghc-options: -fplugin=LiquidHaskell -fplugin-opt=LiquidHaskell:--fast -fplugin-opt=LiquidHaskell:--no-termination -fplugin-opt=LiquidHaskell:--max-case-expand=4 -fplugin-opt=LiquidHaskell:--short-names @@ -100,6 +111,7 @@ hs-source-dirs: test ghc-options: -O2 -threaded -rtsopts -with-rtsopts=-N+ cpp-options: -DO_REWRITE build-depends: base >= 4.7 && < 5, posit >= 2022.0.1,@@ -125,6 +137,17 @@ type: exitcode-stdio-1.0 hs-source-dirs: test main-is: TestReadShow.hs+ ghc-options: -Wall -O2+ build-depends:+ base >= 4.7 && < 5,+ posit+ default-language: Haskell2010++-- Test the Rewrite Rules for FusedOps+benchmark test-posit-fusedRewrite+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: TestFusedRewrite.hs ghc-options: -Wall -O2 build-depends: base >= 4.7 && < 5,
src/Posit.hs view
@@ -140,7 +140,7 @@ -- Perhaps on the chopping block if we are moving to ElementaryFunctions -- Imports for implementing the Transcendental Functions import GHC.Natural (Natural) -- Import the Natural Numbers ℕ (u+2115) for some of the Transcendental Functions-import Data.Ratio ((%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions+import Data.Ratio () -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions, no more (%) now. -- for NFData instance import Control.DeepSeq (NFData, rnf)@@ -230,19 +230,43 @@ -- I'm num trying to get this definition: instance PositC es => Num (Posit es) where -- Addition- (+) = viaRational2 (+)+ (+) = positADD -- Multiplication- (*) = viaRational2 (*)+ (*) = positMULT -- 'abs', Absolute Value, it's like a magnitude of sorts, abs of a posit is the same as abs of the integer representation- abs = viaIntegral abs+ abs = positABS -- 'signum' it is a kind of an representation of directionality, the sign of a number for instance- signum = viaRational signum+ signum = positSIGNUM -- 'fromInteger' rounds the integer into the closest posit number fromInteger int = R $ fromInteger int -- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation- negate = viaIntegral negate+ negate = positNEG+ -- '(-)', minus, but explisit with an implementaiton that will fuse+ (-) = positSUB -- +--+-- To be able to have rewrite rules to function the instance needs to be some un-inlined function+{-# NOINLINE [1] positADD #-}+positADD :: forall es. PositC es => Posit es -> Posit es -> Posit es+positADD = viaRational2 (+)+{-# NOINLINE [1] positMULT #-}+positMULT :: forall es. PositC es => Posit es -> Posit es -> Posit es+positMULT = viaRational2 (*)+{-# NOINLINE [1] positNEG #-}+positNEG :: forall es. PositC es => Posit es -> Posit es+positNEG = viaIntegral negate+{-# NOINLINE [1] positABS #-}+positABS :: forall es. PositC es => Posit es -> Posit es+positABS = viaIntegral abs+{-# NOINLINE [1] positSIGNUM #-}+positSIGNUM :: forall es. PositC es => Posit es -> Posit es+positSIGNUM = viaRational signum+{-# NOINLINE [1] positSUB #-}+positSUB :: forall es. PositC es => Posit es -> Posit es -> Posit es+positSUB = \ p1 p2 -> positADD p1 (positNEG p2)++ -- deriving via Integral Class, for the Integral representation of the posit viaIntegral :: PositC es => (IntN es -> IntN es) -> Posit es -> Posit es viaIntegral f (Posit int) = Posit $ f int@@ -291,6 +315,8 @@ mid = (e2 - e1) / 2 predicate | e2 >= e1 = (<= e3 + mid) | otherwise = (>= e3 + mid)+ fromEnum _ = error "Please do not use 'fromEnum', it is size limited to Int, and can be machine dependant. Please advocate for the function to be size polymorphic of a FixedWidthInteger."+ toEnum _ = error "Please do not use 'toEnum', it is size limited to Int, and can be machine dependant. Please advocate for the function to be size polymorphic of a FixedWidthInteger." -- @@ -300,11 +326,19 @@ -- How the Frac do I get this definition: instance PositC es => Fractional (Posit es) where fromRational = R- + + _ / 0 = NaR+ a / b = viaRational2 (/) a b+ recip 0 = NaR- recip p = viaRational recip p+ recip p = positRECIP p -- +{-# NOINLINE [1] positRECIP #-}+positRECIP :: forall es. PositC es => Posit es -> Posit es+positRECIP = viaRational recip++ -- Rational Instances; Num & Ord Instanced => Real -- -- I for real want this definition:@@ -404,6 +438,24 @@ fsm :: a -> a -> a -> a +#ifdef O_REWRITE+{-# RULES+"posit/fdot4" forall a0 a1 a2 a3 b0 b1 b2 b3. positADD (positADD (positADD (positMULT a0 b0) (positMULT a1 b1)) (positMULT a2 b2)) (positMULT a3 b3) = fdot4 a0 a1 a2 a3 b0 b1 b2 b3 -- (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)+"posit/fsum4" forall a b c d. positADD a (positADD b (positADD c d)) = fsum4 a b c d -- a + b + c + d+"posit/fdot3" forall a1 a2 a3 b1 b2 b3. positADD (positADD (positMULT a1 b1) (positMULT a2 b2)) (positMULT a3 b3) = fdot3 a1 a2 a3 b1 b2 b3 -- (a1 * b1) + (a2 * b2) + (a3 * b3)+"posit/fsum3" forall a b c. positADD a (positADD b c) = fsum3 a b c -- a + b + c+"posit/fsmSub" forall a b c. positSUB a (positMULT b c) = fsm a b c -- a - (b * c)+"posit/fsm" forall a b c. positADD a (positNEG (positMULT b c)) = fsm a b c -- a - (b * c)+"posit/fsmSwaped" forall a b c. positADD (positNEG (positMULT b c)) a = fsm a b c -- negate (b * c) + a+"posit/fma" forall a b c. positADD (positMULT a b) c = fma a b c -- (a * b) + c+"posit/fmaSwaped" forall a b c. positADD c (positMULT a b) = fma a b c -- c + (a * b)+"posit/fam" forall a b c. positMULT (positADD a b) c = fam a b c -- (a + b) * c+"posit/famSwaped" forall a b c. positMULT c (positADD a b) = fam a b c -- c * (a + b)+"posit/fmmsSub" forall a b c d. positSUB (positMULT a b) (positMULT c d) = fmms a b c d -- (a * b) - (c * d)+"posit/fmms" forall a b c d. positADD (positMULT a b) (positNEG (positMULT c d)) = fmms a b c d -- (a * b) - (c * d)+"posit/fmmsSwapped" forall a b c d. positADD (positNEG (positMULT c d)) (positMULT a b) = fmms a b c d -- negate (c * d) + (a * b)+ #-}+#endif -- Rational Instance instance FusedOps Rational where@@ -768,7 +820,7 @@ | x == -1 = -approx_pi/2 | otherwise = approx_atan w where- w = x / approx_sqrt (1 - x^2)+ w = x / approx_sqrt (fsm 1 x x) -- (1 - x^2) approx_acos :: PositC es => Posit es -> Posit es@@ -780,7 +832,7 @@ | x > 0 = approx_atan invw | otherwise = error "Prove it covers for Rational Numbers." where- invw = approx_sqrt (1 - x^2) / x+ invw = approx_sqrt (fsm 1 x x) / x -- (1 - x^2) approx_atan :: PositC es => Posit es -> Posit es@@ -805,14 +857,14 @@ approx_asinh :: PositC es => Posit es -> Posit es approx_asinh NaR = NaR-approx_asinh x = approx_log $ x + approx_sqrt (x^2 + 1)+approx_asinh x = approx_log $ x + approx_sqrt (fma x x 1) -- (x^2 + 1) approx_acosh :: PositC es => Posit es -> Posit es approx_acosh NaR = NaR approx_acosh x | x < 1 = NaR- | otherwise = approx_log $ x + approx_sqrt (x^2 - 1)+ | otherwise = approx_log $ x + approx_sqrt (fma x x (-1)) -- (x^2 - 1) approx_atanh :: forall es. PositC es => Posit es -> Posit es@@ -980,24 +1032,21 @@ -- ========================================================= class AltFloating p where- eps :: p- phi :: p- gamma :: p -> p- sinc :: p -> p- expm1 :: p -> p+ machEps :: p+ approxEq :: p -> p -> Bool+ goldenRatio :: p hypot2 :: p -> p -> p hypot3 :: p -> p -> p -> p hypot4 :: p -> p -> p -> p -> p ---instance PositF es => AltFloating (Posit es) where- phi = 1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374847540880753868917521266338 -- approx_phi 1.6- eps = succ 1.0 - 1.0- gamma = approx_gamma- sinc = approx_sinc- expm1 x =- let b = approx_atanh $ x / 2- in (2 * b) / (1 - b)+instance forall es. PositF es => AltFloating (Posit es) where+ machEps = succ 1.0 - 1.0+ approxEq a b =+ let a' = convert a :: Posit (Prev es)+ b' = convert b :: Posit (Prev es)+ in a' == b'+ goldenRatio = 1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374847540880753868917521266338 -- approx_phi 1.6 hypot2 a b = let a' :: Posit (Next es) = convert a b' :: Posit (Next es) = convert b in convert (approx_sqrt $ a'^2 + b'^2) :: Posit es@@ -1013,46 +1062,14 @@ ----approx_gamma :: forall es. PositC es => Posit es -> Posit es-approx_gamma z = approx_sqrt(2 * approx_pi) * (z `approx_pow` (z - 0.5)) * approx_exp (negate z) * (1 + series)- where- series :: Posit es- series = sum $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n | n <- [1..len]] -- zipWith (\x y -> ) a001163 a001164- lenA = length a001163- lenB = length a001164- len = if lenA == lenB- then lenA- else error "Seiries Numerator and Denominator do not have the same length."-------- Looks like 1 ULP for 0.7813-approx_sinc :: PositC es => Posit es -> Posit es-approx_sinc NaR = NaR-approx_sinc 0 = 1 -- Why the hell not!-approx_sinc theta = approx_sin theta / theta------ -- =====================================================================--- Useful Constants+-- Useful Constants used internally for Elementary Functions -- ===================================================================== -- -- Use the constant, for performance lnOf2 :: PositC es => Posit es lnOf2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202--------a001163 :: [Integer] -- Numerator-a001163 = [1, 1, -139, -571, 163879, 5246819, -534703531, -4483131259, 432261921612371, 6232523202521089, -25834629665134204969, -1579029138854919086429, 746590869962651602203151, 1511513601028097903631961, -8849272268392873147705987190261, -142801712490607530608130701097701]-a001164 :: [Integer] -- Denominator-a001164 = [12, 288, 51840, 2488320, 209018880, 75246796800, 902961561600, 86684309913600, 514904800886784000, 86504006548979712000, 13494625021640835072000, 9716130015581401251840000, 116593560186976815022080000, 2798245444487443560529920000, 299692087104605205332754432000000, 57540880724084199423888850944000000] -- twoMsqrt3 :: PositC es => Posit es
src/Posit/Internal/PositC.hs view
@@ -42,7 +42,8 @@ IntN, FixedWidthInteger(), Max,- Next+ Next,+ Prev ) where import Prelude hiding (exponent,significand)@@ -156,6 +157,21 @@ Next IV_2022 = V_2022 Next V_2022 = V_2022 +type family Prev (es :: ES)+ where+ Prev Z_3_2 = Z_3_2+ Prev I_3_2 = Z_3_2+ Prev II_3_2 = I_3_2+ Prev III_3_2 = II_3_2+ Prev IV_3_2 = III_3_2+ Prev V_3_2 = IV_3_2+ Prev Z_2022 = Z_2022+ Prev I_2022 = Z_2022+ Prev II_2022 = I_2022+ Prev III_2022 = II_2022+ Prev IV_2022 = III_2022+ Prev V_2022 = IV_2022+ -- | The 'FixedWidthInteger' is a Constraint Synonym that contains all -- of the constraints provided by the 'IntN' Type Family. It is a super -- class for the Posit Class.@@ -442,7 +458,7 @@ nBytes = 2^5 -- | `PositF` Constraint Synonym for things that need both the word size and the next higher word size-type PositF es = (PositC es, PositC (Next es))+type PositF es = (PositC es, PositC (Next es), PositC (Prev es)) -- =====================================================================
stack.yaml view
@@ -1,8 +1,9 @@ # This file is attempting to maintain the working Liquid Haskell versions # that coorispond to a specific GHC or Stackage version -# resolver: nightly-2023-07-06 # ghc-9.6.2-resolver: lts-21.9 # ghc-9.4.6+# resolver: nightly-2023-10-16 # ghc-9.6.3+resolver: lts-21.22 # ghc-9.4.8+# resolver: lts-21.19 # ghc-9.4.7 # resolver: lts-20.26 # ghc-9.2.8 # resolver: lts-19.33 # ghc-9.0.2 # resolver: lts-18.28 # ghc-8.10.7
test/Test/Algorithms.hs view
@@ -278,7 +278,7 @@ -- -- Gauss–Legendre algorithm, Seems only accurate to 2-3 ULP, but really slow-funPi1 :: forall es. (PositC es, PositC (Next es)) => Posit es+funPi1 :: forall es. (PositF es) => Posit es funPi1 = go 0 3 1 (recip.sqrt $ 2) (recip 4) 1 where go :: Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es@@ -296,15 +296,15 @@ -- Borwein's algorithm, with quintic convergence, -- gets to 7 ULP in 4 iterations, but really slow due to expensive function evaluations -- quite unstable and will not converge if sqrt is not accurate, which means log must be accurate-funPi2 :: forall es. (PositC es, PositC (Next es)) => Posit es-funPi2 = recip $ go 0 0 0 0.5 (5 / phi^3)+funPi2 :: forall es. (PositF es) => Posit es+funPi2 = recip $ go 0 0 0 0.5 (5 / goldenRatio^3) where- go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es+ go :: forall es. (PositF es) => Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es go !prevA !prevS !n !a !s | prevA == a = a | prevS == s = a- | abs (prevA - a) <= 3*eps = a -- P256 or Posit128, will not reach a fixed point where `prevA == a` it sort of oscelates until divergence occurs, if we test for less than 3*eps it can stop early- | abs (prevS - s) <= 3*eps = a+ | abs (prevA - a) <= 3*machEps = a -- P256 or Posit128, will not reach a fixed point where `prevA == a` it sort of oscelates until divergence occurs, if we test for less than 3*eps it can stop early+ | abs (prevS - s) <= 3*machEps = a | otherwise = let x = 5 / s - 1 y = (x - 1)^2 + 7@@ -343,10 +343,10 @@ -- Borwin's Quadradic Alogrithm 1985-funPi5 :: forall es. (PositC es, PositC (Next es)) => Posit es+funPi5 :: forall es. (PositF es) => Posit es funPi5 = recip $ go 0 0 1 (6 - 4 * sqrt 2) (sqrt 2 - 1) where- go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es+ go :: forall es. (PositF es) => Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es go !prevA !prevY !n a y | prevA == a = a | prevY == y = a@@ -362,7 +362,7 @@ -- ULP: -97 -- Borwin's Cubic Algirthm-funPi6 :: forall es. (PositC es, PositC (Next es)) => Posit es+funPi6 :: forall es. (PositF es) => Posit es funPi6 = recip $ go 0 0 1 (1/3) ((sqrt 3 - 1) / 2) where go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es@@ -435,7 +435,7 @@ -- where the functions are equal, they are not equal. It a class of -- functions with the charicteristic of being a band pass filter in the -- frequency space.--- Shannon wavelet+{- Shannon wavelet funPsiSha1 :: Posit256 -> Posit256 funPsiSha1 NaR = NaR funPsiSha1 t = 2 * sinc (2 * t) - sinc t@@ -445,7 +445,7 @@ funPsiSha2 :: Posit256 -> Posit256 funPsiSha2 NaR = NaR funPsiSha2 t = sinc (t/2) * cos (3*pi*t/2)---+-} -- Shannon wavelet, same as funPsiSha1 but with a factor of pi, with the -- Law: funPsiSha1.(pi*) === funPsiSha3@@ -507,7 +507,7 @@ -funGammaSeriesFused :: forall es. (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaSeriesFused :: forall es. (PositF es) => Posit es -> Posit es funGammaSeriesFused z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series) where series :: Posit es@@ -560,7 +560,7 @@ -- ---funGammaRamanujan :: (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaRamanujan :: (PositF es) => Posit es -> Posit es funGammaRamanujan z = sqrt pi * (x / exp 1)**x * (8*x^3 + 4*x^2 + x + (1/30))**(1/6) where x = z - 1@@ -568,24 +568,24 @@ ---funGammaCalc :: (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaCalc :: (PositF es) => Posit es -> Posit es funGammaCalc z = sqrt (2*pi / z) * ((z / exp 1) * sqrt (z * sinh (recip z) + recip (810 * z^6)))**z -funGammaNemes :: (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaNemes :: (PositF es) => Posit es -> Posit es funGammaNemes z = sqrt (2*pi / z) * (recip (exp 1) * (z + recip (12 * z - recip (10 * z))))**z -funGammaYang :: (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaYang :: (PositF es) => Posit es -> Posit es funGammaYang z = sqrt (2 * pi * x) * (x / exp 1)**x * (x * sinh (recip x))**(x/2) * exp (fromRational (7 % 324) * recip (x^3 * (35 * x^2 + 33))) where x = z - 1 -funGammaChen :: (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaChen :: (PositF es) => Posit es -> Posit es funGammaChen z = sqrt (2 * pi * x) * (x / exp 1)**x * (1 + recip (12*x^3 + (24/7)*x - 0.5))**(x^2 + fromRational (53 % 210)) where x = z - 1 -funGammaXminus1 :: (PositC es, PositC (Next es)) => Posit es -> Posit es+funGammaXminus1 :: (PositF es) => Posit es -> Posit es funGammaXminus1 x = go (x - 1) where go z = sqrt (2 * pi) * exp z ** (negate z) * z ** (z + 0.5)@@ -621,7 +621,25 @@ -- +approx_gamma :: forall es. PositF es => Posit es -> Posit es+approx_gamma z = sqrt(2 * pi) * (z ** (z - 0.5)) * exp (negate z) * (1 + series)+ where+ series :: Posit es+ series = sum $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n | n <- [1..len]] -- zipWith (\x y -> ) a001163 a001164+ lenA = length a001163+ lenB = length a001164+ len = if lenA == lenB+ then lenA+ else error "Seiries Numerator and Denominator do not have the same length."+-- ++-- Looks like 1 ULP for 0.7813+approx_sinc :: PositF es => Posit es -> Posit es+approx_sinc NaR = NaR+approx_sinc 0 = 1 -- Why the hell not!+approx_sinc theta = sin theta / theta+--
+ test/TestFusedRewrite.hs view
@@ -0,0 +1,48 @@++--------------------------------------------------------------------------------------------+-- | Posit Numbers+-- Copyright : (C) 2022-2023 Nathan Waivio+-- License : BSD3+-- Maintainer : Nathan Waivio <nathan.waivio@gmail.com>+-- Stability : Stable+-- Portability : Portable+--+-- Test Suite for a Library implementing standard Posit Numbers+-- to test the Fused Operations Rewrite Rules.+-- +---------------------------------------------------------------------------------------------+++import Posit++main :: IO () +main = do+ print $ "== Let's first take a look at fsm, when examined you shall find that (-) doesn't inline early enough to triger fsm =="+ print $ "posit/fsm pi - (exp(1) * goldenRatio) :: Posit256 : " ++ show (pi - ((exp 1 :: Posit256) * goldenRatio))+ print $ "posit/fsmIn pi + negate (exp(1) * goldenRatio) :: Posit256 : " ++ show (pi + (negate (exp 1 :: Posit256) * goldenRatio))+ print $ "posit/fsmInSwaped negate (exp(1) * goldenRatio) + pi :: Posit256 : " ++ show (negate ((exp 1 :: Posit256) * goldenRatio) + pi)+ print $ "posit/fsmDef `fsm pi exp(1) goldenRatio` :: Posit256 : " ++ show (fsm pi (exp 1 :: Posit256) goldenRatio)+ print $ "== Let's now consider the goodness of fma, it all fuses =="+ print $ "posit/fma (pi * exp(1)) + goldenRatio :: Posit256 : " ++ show ((pi * (exp 1 :: Posit256)) + goldenRatio)+ print $ "posit/fmaSwaped goldenRatio + (pi * exp(1)) :: Posit256 : " ++ show (goldenRatio + (pi * (exp 1 :: Posit256)))+ print $ "posit/fmaDef `fma pi exp(1) goldenRatio` :: Posit256 : " ++ show (fma pi (exp 1 :: Posit256) goldenRatio)+ print $ "== Let's now consider the goodness of fam, it all fuses =="+ print $ "posit/fam (pi + (exp 1)) * goldenRatio :: Posit256 : " ++ show ((pi + (exp 1 :: Posit256)) * goldenRatio)+ print $ "posit/famSwaped goldenRatio * (pi + (exp 1)) :: Posit256 : " ++ show (goldenRatio * (pi + (exp 1 :: Posit256)))+ print $ "posit/famDef `fam pi exp(1) goldenRatio` :: Posit256 : " ++ show (fam pi (exp 1 :: Posit256) goldenRatio)+ print $ "== Let's enjoy `fmms` it now fuses! =="+ print $ "posit/fmms (pi * exp(1)) - (goldenRatio * (sqrt 2)) :: Posit256 : " ++ show ((pi * (exp 1 :: Posit256)) - (goldenRatio * (sqrt 2)))+ print $ "posit/fmmsIn (pi * exp(1)) + negate (goldenRatio * (sqrt 2)) :: Posit256 : " ++ show ((pi * (exp 1 :: Posit256)) + negate (goldenRatio * (sqrt 2)))+ print $ "posit/fmmsInSwap (negate (goldenRatio * (sqrt 2)) + (pi * exp(1))) :: Posit256 : " ++ show (negate (goldenRatio * (sqrt 2)) + (pi * (exp 1 :: Posit256)))+ print $ "posit/fmmsDef `fmms pi exp(1) goldenRatio (sqrt 2)` :: Posit256 : " ++ show (fmms pi (exp 1 :: Posit256) goldenRatio (sqrt 2))+ print $ "== Let's see about fused dot products =="+ print $ "posit/fdot3 (pi * exp(1)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)) :: Posit256 : " ++ show ((pi * (exp 1 :: Posit256)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)))+ print $ "posit/fdot3Def `fdot3 pi goldenRatio (sqrt goldenRatio) (exp 1) (sqrt 2) (sqrt 5)` :: Posit256 : " ++ show (fdot3 pi goldenRatio (sqrt goldenRatio) (exp 1) (sqrt 2) (sqrt 5) :: Posit256)+ print $ "posit/fdot3 (pi * exp(1)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)) :: P256 : " ++ show ((pi * (exp 1 :: P256)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)))+ print $ "posit/fdot3Def `fdot3 pi goldenRatio (sqrt goldenRatio) (exp 1) (sqrt 2) (sqrt 5)` :: P256 : " ++ show (fdot3 pi goldenRatio (sqrt goldenRatio) (exp 1) (sqrt 2) (sqrt 5) :: P256)+ print $ "== Let's see about fused dot products =="+ print $ "posit/fdot4 (pi * exp(1)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)) + ((recip pi) * (recip goldenRatio)) :: Posit256 : " ++ show ((pi * (exp 1 :: Posit256)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)) + ((recip pi) * (recip goldenRatio)))+ print $ "posit/fdot4Def `fdot4 pi goldenRatio (sqrt goldenRatio) (recip pi) (exp 1) (sqrt 2) (sqrt 5) (recip goldenRatio)` :: Posit256 : " ++ show (fdot4 pi goldenRatio (sqrt goldenRatio) (recip pi) (exp 1) (sqrt 2) (sqrt 5) (recip goldenRatio) :: Posit256)+ print $ "posit/fdot4 (pi * exp(1)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)) + ((recip pi) * (recip goldenRatio)) :: P256 : " ++ show ((pi * (exp 1 :: P256)) + (goldenRatio * (sqrt 2)) + ((sqrt goldenRatio) * (sqrt 5)) + ((recip pi) * (recip goldenRatio)))+ print $ "posit/fdot4Def `fdot4 pi goldenRatio (sqrt goldenRatio) (recip pi) (exp 1) (sqrt 2) (sqrt 5) (recip goldenRatio)` :: P256 : " ++ show (fdot4 pi goldenRatio (sqrt goldenRatio) (recip pi) (exp 1) (sqrt 2) (sqrt 5) (recip goldenRatio) :: P256)+
test/TestPosit.hs view
@@ -63,62 +63,62 @@ print $ "exp(1)**(pi*sqrt 163) :: Posit8 " ++ show (exp(1 :: Posit8) ** (pi * sqrt 163)) -- print $ "exp(1)**(pi*sqrt 163) :: P8 " ++ show (exp(1 :: P8) ** (pi * sqrt 163)) -- -- | 'EPS'- print $ "Machine epsilon Posit8 ~1.0: " ++ show (eps :: Posit8) -- succ (Posit int) = Posit (succ int)- print $ "Machine epsilon Posit16 ~1.0: " ++ show (eps :: Posit16) -- - print $ "Machine epsilon Posit32 ~1.0: " ++ show (eps :: Posit32) -- - print $ "Machine epsilon Posit64 ~1.0: " ++ show (eps :: Posit64) -- - print $ "Machine epsilon Posit128 ~1.0: " ++ show (eps :: Posit128) -- - print $ "Machine epsilon Posit256 ~1.0: " ++ show (eps :: Posit256) -- - print $ "Machine epsilon P8 ~1.0: " ++ show (eps :: P8) -- succ (Posit int) = Posit (succ int)- print $ "Machine epsilon P16 ~1.0: " ++ show (eps :: P16) -- - print $ "Machine epsilon P32 ~1.0: " ++ show (eps :: P32) -- - print $ "Machine epsilon P64 ~1.0: " ++ show (eps :: P64) -- - print $ "Machine epsilon P128 ~1.0: " ++ show (eps :: P128) -- - print $ "Machine epsilon P256 ~1.0: " ++ show (eps :: P256) -- + print $ "Machine epsilon Posit8 ~1.0: " ++ show (machEps :: Posit8) -- succ (Posit int) = Posit (succ int)+ print $ "Machine epsilon Posit16 ~1.0: " ++ show (machEps :: Posit16) -- + print $ "Machine epsilon Posit32 ~1.0: " ++ show (machEps :: Posit32) -- + print $ "Machine epsilon Posit64 ~1.0: " ++ show (machEps :: Posit64) -- + print $ "Machine epsilon Posit128 ~1.0: " ++ show (machEps :: Posit128) -- + print $ "Machine epsilon Posit256 ~1.0: " ++ show (machEps :: Posit256) -- + print $ "Machine epsilon P8 ~1.0: " ++ show (machEps :: P8) -- succ (Posit int) = Posit (succ int)+ print $ "Machine epsilon P16 ~1.0: " ++ show (machEps :: P16) -- + print $ "Machine epsilon P32 ~1.0: " ++ show (machEps :: P32) -- + print $ "Machine epsilon P64 ~1.0: " ++ show (machEps :: P64) -- + print $ "Machine epsilon P128 ~1.0: " ++ show (machEps :: P128) -- + print $ "Machine epsilon P256 ~1.0: " ++ show (machEps :: P256) -- print $ "Does (1 - 1) == 0 ?: " ++ show ((1 - 1) == (0 :: Posit256)) -- [(1 - 1) == zero | zero = 0 :: Posit es, es <- Z .. V] print $ "truth :: Posit256: " ++ show (1.2720196495140689642524224617374914917156080418400962486166403825392975755360680118303842149884602585385141476367280265057103381 :: Posit256)- print $ "sqrt phi :: Posit256: " ++ show (sqrt phi :: Posit256)- print $ "sqrt phi :: Posit128: " ++ show (sqrt phi :: Posit128)- print $ "sqrt phi :: Posit64: " ++ show (sqrt phi :: Posit64)- print $ "sqrt phi :: Posit32: " ++ show (sqrt phi :: Posit32)- print $ "sqrt phi :: Posit16: " ++ show (sqrt phi :: Posit16)- print $ "sqrt phi :: Posit8: " ++ show (sqrt phi :: Posit8)+ print $ "sqrt goldenRatio :: Posit256: " ++ show (sqrt goldenRatio :: Posit256)+ print $ "sqrt goldenRatio :: Posit128: " ++ show (sqrt goldenRatio :: Posit128)+ print $ "sqrt goldenRatio :: Posit64: " ++ show (sqrt goldenRatio :: Posit64)+ print $ "sqrt goldenRatio :: Posit32: " ++ show (sqrt goldenRatio :: Posit32)+ print $ "sqrt goldenRatio :: Posit16: " ++ show (sqrt goldenRatio :: Posit16)+ print $ "sqrt goldenRatio :: Posit8: " ++ show (sqrt goldenRatio :: Posit8) print $ "truth :: P256: " ++ show (1.2720196495140689642524224617374914917156080418400962486166403825392975755360680118303842149884602585385141476367280265057103381 :: P256)- print $ "sqrt phi :: P256: " ++ show (sqrt phi :: P256)- print $ "sqrt phi :: P128: " ++ show (sqrt phi :: P128)- print $ "sqrt phi :: P64: " ++ show (sqrt phi :: P64)- print $ "sqrt phi :: P32: " ++ show (sqrt phi :: P32)- print $ "sqrt phi :: P16: " ++ show (sqrt phi :: P16)- print $ "sqrt phi :: P8: " ++ show (sqrt phi :: P8)+ print $ "sqrt goldenRatio :: P256: " ++ show (sqrt goldenRatio :: P256)+ print $ "sqrt goldenRatio :: P128: " ++ show (sqrt goldenRatio :: P128)+ print $ "sqrt goldenRatio :: P64: " ++ show (sqrt goldenRatio :: P64)+ print $ "sqrt goldenRatio :: P32: " ++ show (sqrt goldenRatio :: P32)+ print $ "sqrt goldenRatio :: P16: " ++ show (sqrt goldenRatio :: P16)+ print $ "sqrt goldenRatio :: P8: " ++ show (sqrt goldenRatio :: P8) {- let truthPosit256 = 0.8956731517052878608869612167009786079379812529831641161347143256 :: Posit256 -- 0.89566032673209158354178209470474131001971567786620187475744721557 :: Posit256 -- 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: Posit256 let truthP256 = 0.8956731517052878608869612167009786079379812529831641161347143256 :: P256 -- 0.89566032673209158354178209470474131001971567786620187475744721557 :: P256 -- 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: P256- eval "Standard: gamma(phi) :: Posit256 " (gamma (phi)) truthPosit256- eval "Standard: gamma(phi) :: P256 " (gamma (phi)) truthP256- eval "Fused Gamma: gamma(phi) :: Posit256 " (funGammaSeriesFused (phi)) truthPosit256- eval "Fused Gamma: gamma(phi) :: P256 " (funGammaSeriesFused (phi)) truthP256- eval "Ramanujan Gamma: gamma(phi) :: Posit256 " (funGammaRamanujan (phi)) truthPosit256- eval "Ramanujan Gamma: gamma(phi) :: P256 " (funGammaRamanujan (phi)) truthP256- eval "Calc Gamma: gamma(phi) :: Posit256 " (funGammaCalc (phi)) truthPosit256- eval "Calc Gamma: gamma(phi) :: P256 " (funGammaCalc (phi)) truthP256- eval "Nemes Gamma: gamma(phi) :: Posit256 " (funGammaNemes (phi)) truthPosit256- eval "Nemes Gamma: gamma(phi) :: P256 " (funGammaNemes (phi)) truthP256- eval "Yang Gamma: gamma(phi) :: Posit256 " (funGammaYang (phi)) truthPosit256- eval "Yang Gamma: gamma(phi) :: P256 " (funGammaYang (phi)) truthP256- eval "Chen Gamma: gamma(phi) :: Posit256 " (funGammaChen (phi)) truthPosit256- eval "Chen Gamma: gamma(phi) :: P256 " (funGammaChen (phi)) truthP256- eval "Gamma (x - 1): gamma(phi) :: Posit256 " (funGammaXminus1 (phi)) truthPosit256- eval "Gamma (x - 1): gamma(phi) :: P256 " (funGammaXminus1 (phi)) truthP256- eval "Calcuation of gamma(phi) using lngamma :: Posit256" (funGammaViaLngamma (phi)) truthPosit256- eval "Calcuation of gamma(phi) using lngamma :: P256" (funGammaViaLngamma (phi)) truthP256- eval "Wolfram alpha: gamma(phi) :: Posit256 " truthPosit256 truthPosit256- eval "Wolfram alpha: gamma(phi) :: P256 " truthP256 truthP256+ eval "Standard: gamma(goldenRatio) :: Posit256 " (gamma (goldenRatio)) truthPosit256+ eval "Standard: gamma(goldenRatio) :: P256 " (gamma (goldenRatio)) truthP256+ eval "Fused Gamma: gamma(goldenRatio) :: Posit256 " (funGammaSeriesFused (goldenRatio)) truthPosit256+ eval "Fused Gamma: gamma(goldenRatio) :: P256 " (funGammaSeriesFused (goldenRatio)) truthP256+ eval "Ramanujan Gamma: gamma(goldenRatio) :: Posit256 " (funGammaRamanujan (goldenRatio)) truthPosit256+ eval "Ramanujan Gamma: gamma(goldenRatio) :: P256 " (funGammaRamanujan (goldenRatio)) truthP256+ eval "Calc Gamma: gamma(goldenRatio) :: Posit256 " (funGammaCalc (goldenRatio)) truthPosit256+ eval "Calc Gamma: gamma(goldenRatio) :: P256 " (funGammaCalc (goldenRatio)) truthP256+ eval "Nemes Gamma: gamma(goldenRatio) :: Posit256 " (funGammaNemes (goldenRatio)) truthPosit256+ eval "Nemes Gamma: gamma(goldenRatio) :: P256 " (funGammaNemes (goldenRatio)) truthP256+ eval "Yang Gamma: gamma(goldenRatio) :: Posit256 " (funGammaYang (goldenRatio)) truthPosit256+ eval "Yang Gamma: gamma(goldenRatio) :: P256 " (funGammaYang (goldenRatio)) truthP256+ eval "Chen Gamma: gamma(goldenRatio) :: Posit256 " (funGammaChen (goldenRatio)) truthPosit256+ eval "Chen Gamma: gamma(goldenRatio) :: P256 " (funGammaChen (goldenRatio)) truthP256+ eval "Gamma (x - 1): gamma(goldenRatio) :: Posit256 " (funGammaXminus1 (goldenRatio)) truthPosit256+ eval "Gamma (x - 1): gamma(goldenRatio) :: P256 " (funGammaXminus1 (goldenRatio)) truthP256+ eval "Calcuation of gamma(goldenRatio) using lngamma :: Posit256" (funGammaViaLngamma (goldenRatio)) truthPosit256+ eval "Calcuation of gamma(goldenRatio) using lngamma :: P256" (funGammaViaLngamma (goldenRatio)) truthP256+ eval "Wolfram alpha: gamma(goldenRatio) :: Posit256 " truthPosit256 truthPosit256+ eval "Wolfram alpha: gamma(goldenRatio) :: P256 " truthP256 truthP256 -} let truth = 5.0431656433600286513118821892854247103235901754138463603020001967777869609108929428415187821843384653305404495551887666992776792 :: Posit256- eval "Standard: exp(phi):" (exp (phi)) truth- eval "Taylor: exp(phi):" (funExp2 funExpTaylor (phi / log 2)) truth- eval "Tuma: exp(phi):" (funExp2 funExpTuma (phi / log 2)) truth- eval "Wolfram Alpha: exp(phi):" truth truth+ eval "Standard: exp(goldenRatio):" (exp (goldenRatio)) truth+ eval "Taylor: exp(goldenRatio):" (funExp2 funExpTaylor (goldenRatio / log 2)) truth+ eval "Tuma: exp(goldenRatio):" (funExp2 funExpTuma (goldenRatio / log 2)) truth+ eval "Wolfram Alpha: exp(goldenRatio):" truth truth let truth = 2.6881171418161354484126255515800135873611118773741922415191608615280287034909564914158871097219845710811670879190576068697e43 :: Posit256 eval "Standard: exp(100):" (exp (100)) truth eval "Taylor: exp(100):" (funExp2 funExpTaylor (100 / log 2)) truth@@ -140,10 +140,10 @@ eval "Tuma: exp(-1000):" (funExp2 funExpTuma (-1000 / log 2)) truth eval "Wolfram Alpha: exp(-1000):" truth truth let truth = 0.4812118250596034474977589134243684231351843343856605196610181688401638676082217744120094291227234749972318399582936564112725683 :: Posit256- eval "Standard: log(phi):" (log (phi)) truth- eval "Taylor: log(phi):" (funLogDomainReduction funLogTaylor (phi)) truth- eval "Tuma: log(phi):" (funLogDomainReduction funLogTuma (phi)) truth- eval "Wolfram Alpha: log(phi):" truth truth+ eval "Standard: log(goldenRatio):" (log (goldenRatio)) truth+ eval "Taylor: log(goldenRatio):" (funLogDomainReduction funLogTaylor (goldenRatio)) truth+ eval "Tuma: log(goldenRatio):" (funLogDomainReduction funLogTuma (goldenRatio)) truth+ eval "Wolfram Alpha: log(goldenRatio):" truth truth let truth = -4.6051701859880913680359829093687284152022029772575459520666558019351452193547049604719944101791965966839355680845724972668190 :: Posit256 eval "Standard: log(1/100):" (log (1/100)) truth eval "Taylor: log(1/100):" (funLogDomainReduction funLogTaylor (1/100)) truth@@ -155,8 +155,8 @@ eval "Tuma: log(1/1000):" (funLogDomainReduction funLogTuma (1/1000)) truth eval "Wolfram Alpha: log(1/1000):" truth truth let truth = 4.5347571611551792889915884948567915637887680293971326427244942079650289300980475282698882636812383679690567084677326507550787791 :: Posit256- eval "Standard: phi**pi:" ((phi) ** pi) truth- eval "Wolfram Alpha: phi**pi:" truth truth+ eval "Standard: goldenRatio**pi:" ((goldenRatio) ** pi) truth+ eval "Wolfram Alpha: goldenRatio**pi:" truth truth let tPiPosit256 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: Posit256 let tPiP256 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: P256 let tPiPosit128 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: Posit128