diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,16 @@ # Changelog for Posit Numbers +# posit-3.2.0.1 + + * Refactored `IntN` Type Family to be a closed type family instead of an associated type family + * Refactored `IntN` constraints to use `ConstraintKinds` and made that to be a Super Class of `PositC` to improve the encapsulation the Constraints of the internal implementation + * Refactored `PositC` to make use of `ConstrainedClassMethods` vastly reducing code duplication + * Eliminated the `FlexableContexts` Language Extension from Posit.hs Interface, since the `InN` constraints no longer bleed into that file + * Added test of Heegner numbers (almost integers) + * Added test of various Gamma Function approximations + * Improved function names in the Orphan Instance for `Storable` ( `Word128` ) + * Improved documentation + ## posit-3.2.0.0 * Posit Standard 3.2 [Posit Standard] (https://posithub.org/docs/posit_standard.pdf) diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,10 +1,11 @@ -# posit 3.2.0.0 +# posit 3.2.0.1 The [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf), -where Real numbers are approximated by Maybe Rational. The Posit type -is mapped to a 2's complement integer type; smoothly and with tapering -precision, in a similar way to the projective real line. The 'posit' -library implements the following standard classes: +where Real numbers are approximated by Maybe Rational. The Posit +Numbers are a drop in replacement for `Float` or `Double` mapped to a +2's complement integer type; smoothly and with tapering precision, in a +similar way to the projective real line. The 'posit' library implements +the following standard classes: * Show * Eq @@ -24,7 +25,7 @@ * Floating -- Mathematical functions such as logarithm, exponential, trigonometric, and hyperbolic functions. Warning! May induce trance. The Posits are indexed by the type (es :: ES) where exponent size and -word size are related. In `posit-3.2.0.0` es is instantiated as Z, I, +word size are related. In `posit-3.2.0.1` es is instantiated as Z, I, II, III, IV, V. The word size (in bits) of the value is `= 8 * 2^es`, that is `2^es` bytes. The Types: 'Posit8', 'Posit16', 'Posit32', 'Posit64', 'Posit128', and 'Posit256' are implemented and include a diff --git a/posit.cabal b/posit.cabal --- a/posit.cabal +++ b/posit.cabal @@ -1,7 +1,7 @@ cabal-version: 1.12 name: posit -version: 3.2.0.0 +version: 3.2.0.1 description: The Posit Number format. Please see the README on GitHub at homepage: https://github.com/waivio/posit#readme bug-reports: https://github.com/waivio/posit/issues @@ -55,7 +55,7 @@ ghc-options: -Wall -O2 if flag(do-liquid) - ghc-options: -fplugin=LiquidHaskell -fplugin-opt=LiquidHaskell:--fast -fplugin-opt=LiquidHaskell:--max-case-expand=4 -fplugin-opt=LiquidHaskell:--no-termination -fplugin-opt=LiquidHaskell:--short-names + 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 if flag(do-no-storable) cpp-options: -DO_NO_STORABLE @@ -64,7 +64,7 @@ cpp-options: -DO_NO_ORPHANS if flag(do-liquid) - cpp-options: -DO_LIQUID -DO_NO_STORABLE -DO_NO_READ -DO_NO_SHOW + cpp-options: -DO_LIQUID -DO_NO_STORABLE if flag(do-test) cpp-options: -DO_TEST diff --git a/src/Posit.hs b/src/Posit.hs --- a/src/Posit.hs +++ b/src/Posit.hs @@ -21,7 +21,6 @@ {-# LANGUAGE BangPatterns #-} -- Added Strictness for some fixed point algorithms {-# LANGUAGE PatternSynonyms #-} -- for a nice NaR interface {-# LANGUAGE FlexibleInstances #-} -- To make instances for each specific type [Posit8 .. Posit256] -{-# LANGUAGE FlexibleContexts #-} -- Allow non-type variables in the constraints {-# LANGUAGE TypeApplications #-} -- To apply types: @Type, it seems to select the specific class instance, when GHC is not able to reason about things, commenting this out shows an interesting interface {-# LANGUAGE MultiParamTypeClasses #-} -- To convert between Posit Types {-# LANGUAGE ScopedTypeVariables #-} -- To reduce some code duplication @@ -31,6 +30,7 @@ {-# OPTIONS_GHC -Wno-type-defaults #-} -- Turn off noise {-# OPTIONS_GHC -Wno-unused-top-binds #-} -- Turn off noise + -- ---- -- Posit numbers implementing: -- @@ -54,7 +54,8 @@ -- ---- module Posit -(-- * Main Exported Types +(Posit(), + -- * Main Exported Types Posit8, -- |An 8-bit Posit number with 'es' ~ 'Z' Posit16, -- |An 16-bit Posit number with 'es' ~ 'I' Posit32, -- |An 32-bit Posit number with 'es' ~ 'II' @@ -72,8 +73,10 @@ -- * Posits are Convertable between different Posit representations Convertible(..), +#ifndef O_NO_SHOW -- * Additional functions to show the Posit in different formats AltShow(..), +#endif -- * Additional Special Functions AltFloating(..), @@ -86,6 +89,7 @@ viaRational4, viaRational6, viaRational8, + #ifdef O_TEST -- * Alternative algorithms for test purposes funExp, @@ -93,6 +97,13 @@ funExpTaylor, funLogTaylor, funExpTuma, + funGammaSeriesFused, + funGammaRamanujan, + funGammaCalc, + funGammaNemes, + funGammaYang, + funGammaChen, + funGammaXminus1, funLogTuma, funLogDomainReduction, funPi1, @@ -103,14 +114,14 @@ funPsiSha2, funPsiSha3 #endif + ) where import Prelude hiding (rem) -- Imports for Show and Read Instances -import Data.Scientific (Scientific - ,scientificP +import Data.Scientific (scientificP ,fromRationalRepetendUnlimited ,formatScientific ,FPFormat(Generic)) -- Used to print/show and read the rational value @@ -130,7 +141,7 @@ -- Imports for Storable Instance import Foreign.Storable (Storable, sizeOf, alignment, peek, poke) -- Used for Storable Instances of Posit -import Foreign.Ptr (Ptr, plusPtr, castPtr) -- Used for dealing with Pointers for the Posit Storable Instance +import Foreign.Ptr (Ptr, castPtr) -- Used for dealing with Pointers for the Posit Storable Instance -- would like to: @@ -138,7 +149,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 (Rational, (%)) -- 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 import Debug.Trace (trace) -- temporary for debug purposes @@ -148,7 +159,7 @@ -- ===================================================================== -- The machine implementation of the Posit encoding/decoding -import Posit.Internal.PositC (ES(..), PositC(..)) -- The main internal implementation details +import Posit.Internal.PositC -- The main internal implementation details -- |Base GADT rapper type, that uses the Exponent Size kind to index the various implementations @@ -156,7 +167,7 @@ Posit :: PositC es => !(IntN es) -> Posit es -- |Not a Real Number, the Posit is like a Maybe type, it's either a real number or not -pattern NaR :: (PositC es) => Posit es +pattern NaR :: PositC es => Posit es pattern NaR <- (Posit (decode -> Nothing)) where NaR = Posit unReal -- @@ -182,7 +193,7 @@ #ifndef O_NO_SHOW -- Show -- -instance forall es. (PositC es) => Show (Posit es) where +instance PositC es => Show (Posit es) where show NaR = "NaR" show (R r) = formatScientific Generic (Just $ decimalPrec @es) (fst.fromRationalRepetendUnlimited $ r) -- @@ -193,7 +204,7 @@ -- Two Posit Numbers are Equal if their Finite Precision Integer representation is Equal -- -- All things equal I would rather write it like this: -instance forall es. (Eq (IntN es)) => Eq (Posit es) where +instance PositC es => Eq (Posit es) where (Posit int1) == (Posit int2) = int1 == int2 -- @@ -202,7 +213,7 @@ -- Two Posit Numbers are ordered by their Finite Precision Integer representation -- -- Ordinarily I would only like one instance to cover them all -instance forall es. (Ord (IntN es), PositC es) => Ord (Posit es) where +instance PositC es => Ord (Posit es) where compare (Posit int1) (Posit int2) = compare int1 int2 -- @@ -211,7 +222,7 @@ -- Num -- -- I'm num trying to get this definition: -instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Num (Posit es) where +instance PositC es => Num (Posit es) where -- Addition (+) = viaRational2 (+) -- Multiplication @@ -227,7 +238,7 @@ -- -- deriving via Integral Class, for the Integral representation of the posit -viaIntegral :: forall es. PositC es => (IntN es -> IntN es) -> Posit es -> Posit es +viaIntegral :: PositC es => (IntN es -> IntN es) -> Posit es -> Posit es viaIntegral f (Posit int) = Posit $ f int -- @@ -235,7 +246,7 @@ -- Enum-ish, A Posit has a Successor and Predecessor so its an ordinal number, as per Posit standard next, prior -- The Posit Standard requires 2's complement integer overflow to be ignored -instance forall es. (Num (IntN es), Enum (IntN es), Ord (IntN es), PositC es) => Enum (Posit es) where +instance PositC es => Enum (Posit es) where -- succ (Posit int) = Posit (int + 1) succ = viaIntegral (+1) -- succ = viaIntegral succ -- Non-compliant, runtime error pred NaR, and worse it is Int64 for types of greater precision, probably because of Preludes gross abomination of toEnum/fromEnum @@ -281,7 +292,7 @@ -- Fractional Instances; (Num => Fractional) -- -- How the Frac do I get this definition: -instance forall es. (Num (IntN es), Ord (IntN es), Eq (IntN es), PositC es) => Fractional (Posit es) where +instance PositC es => Fractional (Posit es) where fromRational = R recip 0 = NaR @@ -291,7 +302,7 @@ -- Rational Instances; Num & Ord Instanced => Real -- -- I for real want this definition: -instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Real (Posit es) where +instance PositC es => Real (Posit es) where toRational NaR = error "Your input is Not a Real or Rational (NaR) number, please try again!" toRational (R r) = r -- @@ -299,25 +310,25 @@ -- Implementing instances via Rational Data Type's instance, -- The function checks for NaR, to protect against the runtime error 'toRational' would generate if called with a NaR value -- Unary::Arity NaR guarded pass through with wrapping and unwrapping use of a Rational function -viaRational :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational) -> Posit es -> Posit es +viaRational :: PositC es => (Rational -> Rational) -> Posit es -> Posit es viaRational _ NaR = NaR viaRational f (R r) = fromRational $ f r -- Binary NaR guarded pass through with wrapping and unwrapping use of a Rational function -viaRational2 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es +viaRational2 :: PositC es => (Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es viaRational2 _ NaR _ = NaR viaRational2 _ _ NaR = NaR viaRational2 f (R r1) (R r2) = R $ r1 `f` r2 -- Ternary NaR guarded pass through with wrapping and unwrapping use of a Rational function -viaRational3 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es +viaRational3 :: PositC es => (Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es viaRational3 _ NaR _ _ = NaR viaRational3 _ _ NaR _ = NaR viaRational3 _ _ _ NaR = NaR viaRational3 f (R r1) (R r2) (R r3) = R $ f r1 r2 r3 -- Quaternary NaR guarded pass through with wrapping and unwrapping use of a Rational function -viaRational4 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es +viaRational4 :: PositC es => (Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es viaRational4 _ NaR _ _ _ = NaR viaRational4 _ _ NaR _ _ = NaR viaRational4 _ _ _ NaR _ = NaR @@ -325,7 +336,7 @@ viaRational4 f (R r0) (R r1) (R r2) (R r3) = R $ f r0 r1 r2 r3 -- Senary NaR guarded pass through with wrapping and unwrapping use of a Rational function -viaRational6 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es +viaRational6 :: PositC es => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es viaRational6 _ NaR _ _ _ _ _ = NaR viaRational6 _ _ NaR _ _ _ _ = NaR viaRational6 _ _ _ NaR _ _ _ = NaR @@ -335,7 +346,7 @@ viaRational6 f (R a1) (R a2) (R a3) (R b1) (R b2) (R b3) = R $ f a1 a2 a3 b1 b2 b3 -- Octonary NaR guarded pass through with wrapping and unwrapping use of a Rational function -viaRational8 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es +viaRational8 :: PositC es => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es viaRational8 _ NaR _ _ _ _ _ _ _ = NaR viaRational8 _ _ NaR _ _ _ _ _ _ = NaR viaRational8 _ _ _ NaR _ _ _ _ _ = NaR @@ -351,7 +362,7 @@ -- Bounded, bounded to what?!? To the ℝ! NaR is out of bounds!!! -- -- I'm bound to want this definition: -instance forall es. PositC es => Bounded (Posit es) where +instance PositC es => Bounded (Posit es) where -- 'minBound' the most negative number represented minBound = Posit mostNegVal -- 'maxBound' the most positive number represented @@ -411,7 +422,7 @@ -- -- -instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => FusedOps (Posit es) where +instance PositC es => FusedOps (Posit es) where -- Fused Subtract Multiply fsm = viaRational3 fsm -- Fuse Multiply Add @@ -460,7 +471,7 @@ class Convertible a b where convert :: a -> b -instance forall es1 es2. (PositC es1, PositC es2, Ord (IntN es1), Ord (IntN es2), Num (IntN es1), Num (IntN es2)) => Convertible (Posit es1) (Posit es2) where +instance (PositC es1, PositC es2) => Convertible (Posit es1) (Posit es2) where convert NaR = NaR convert (R r) = R r -- @@ -484,7 +495,7 @@ -- -- -instance forall es. (Show (IntN es), Ord (IntN es), Num (IntN es), PositC es) => AltShow (Posit es) where +instance PositC es => AltShow (Posit es) where displayBinary (Posit int) = displayBin int displayIntegral (Posit int) = show int @@ -494,7 +505,7 @@ displayDecimal = viaShowable (fst.fromRationalRepetendUnlimited) -- -viaShowable :: forall es a. (Show a, Ord (IntN es), Num (IntN es), PositC es) => (Rational -> a) -> Posit es -> String +viaShowable :: (Show a, PositC es) => (Rational -> a) -> Posit es -> String viaShowable _ NaR = "NaR" viaShowable f (R r) = show $ f r #endif @@ -505,7 +516,7 @@ -- ===================================================================== -- -instance forall es. (PositC es) => Read (Posit es) where +instance PositC es => Read (Posit es) where readPrec = parens $ do x <- lexP @@ -528,7 +539,7 @@ -- #ifndef O_NO_STORABLE -- -instance forall es. (Storable (IntN es), PositC es) => Storable (Posit es) where +instance PositC es => Storable (Posit es) where sizeOf _ = fromIntegral $ nBytes @es alignment _ = fromIntegral $ nBytes @es peek ptr = do @@ -545,12 +556,12 @@ -- ===================================================================== -- -instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => RealFrac (Posit es) where +instance PositC es => RealFrac (Posit es) where -- properFraction :: Integral b => a -> (b, a) properFraction = viaRationalErrTrunkation "NaR value is not a RealFrac" properFraction -- -viaRationalErrTrunkation :: forall es a. (Num (IntN es), (Ord (IntN es)), PositC es, Integral a) => String -> (Rational -> (a, Rational)) -> Posit es -> (a, Posit es) +viaRationalErrTrunkation :: PositC es => String -> (Rational -> (a, Rational)) -> Posit es -> (a, Posit es) viaRationalErrTrunkation err _ NaR = error err viaRationalErrTrunkation _ f (R r) = let (int, r') = f r @@ -560,7 +571,7 @@ -- === Real Float === -- ===================================================================== -- -instance forall es. (Eq (IntN es), Ord (IntN es), Num (IntN es), Floating (Posit es), PositC es) => RealFloat (Posit es) where +instance (Floating (Posit es), PositC es) => RealFloat (Posit es) where isIEEE _ = False isDenormalized _ = False isNegativeZero _ = False @@ -773,6 +784,8 @@ | otherwise = funPhi (Posit x') where (Posit x') = (px^2 + 2*px) / (px^2 + 1) + -- LiquidHaskell is telling me this is unsafe if px is imaginary + -- lucky for us Posit256 is not imaginary -- calculate atan(1/2^n) @@ -1131,6 +1144,7 @@ in go next ((a' + b')^2 / (4 * t')) a' b' t' p' -- +#ifndef O_NO_SHOW -- 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 @@ -1148,6 +1162,7 @@ s' = 25 / ((z + x/z + 1)^2 * s) in go a (n+1) (trace (show a') a') s' -- +#endif -- Bailey–Borwein–Plouffe (BBP) formula, to 1-2 ULP, and blazing fast, converges in 60 iterations @@ -1253,7 +1268,9 @@ -- Interestingly enough, wikipedia defines two alternative solutions -- for the Shannon Wavelet, eventhough there are infinite solutions --- where the functions are equal, they are not equal +-- 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 funPsiSha1 :: Posit256 -> Posit256 funPsiSha1 NaR = NaR @@ -1266,18 +1283,21 @@ funPsiSha2 t = funSinc (t/2) * cos (3*pi*t/2) -- --- Shannon wavelet +-- Shannon wavelet, same as funPsiSha1 but with a factor of pi, with the +-- Law: funPsiSha1.(pi*) === funPsiSha3 +-- or : funPsiSha1 === funpsiSha3.(/pi) +-- Posit256 seems to hold to a few ULP funPsiSha3 :: Posit256 -> Posit256 funPsiSha3 NaR = NaR funPsiSha3 0 = 1 -- Why the hell not! funPsiSha3 t = let pit = pi * t - invpit = recip $ pit + invpit = recip pit in invpit * (sin (2 * pit) - sin pit) -- --- funPsiSha1.(pi*) === funPsiSha3 + -- Using the CORDIC domain reduction and some approximation function funLogDomainReduction :: (Posit256 -> Posit256) -> Posit256 -> Posit256 funLogDomainReduction _ NaR = NaR @@ -1346,6 +1366,17 @@ len = if lenA == lenB then lenA else error "Seiries Numerator and Denominator do not have the same length." + +funGammaSeriesFused :: Posit256 -> Posit256 +funGammaSeriesFused z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series) + where + series :: Posit256 + series = fsumL $ 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." -- funGammaCalc :: Posit256 -> Posit256 @@ -1365,4 +1396,7 @@ where x = z - 1 - +funGammaXminus1 :: Posit256 -> Posit256 +funGammaXminus1 x = go (x - 1) + where + go z = sqrt (2 * pi) * exp z ** (negate z) * z ** (z + 0.5) diff --git a/src/Posit/Internal/PositC.hs b/src/Posit/Internal/PositC.hs --- a/src/Posit/Internal/PositC.hs +++ b/src/Posit/Internal/PositC.hs @@ -20,6 +20,8 @@ {-# LANGUAGE AllowAmbiguousTypes #-} -- The Haskell/GHC Type checker seems to have trouble things in the PositC class {-# LANGUAGE ScopedTypeVariables #-} -- To reduce some code duplication {-# LANGUAGE FlexibleContexts #-} -- To reduce some code duplication by claiming the type family provides some constraints, that GHC can't do without fully evaluating the type family +{-# LANGUAGE ConstrainedClassMethods #-} -- Allows constraints on class methods so default implementations of methods with Type Families can be implemented +{-# LANGUAGE ConstraintKinds #-} -- Simplify all of the constraints into a combinded constraint for the super class constraint {-# LANGUAGE CPP #-} -- To remove Storable instances to remove noise when performing analysis of Core {-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-} -- Turn off noise {-# OPTIONS_GHC -Wno-type-defaults #-} -- Turn off noise @@ -33,7 +35,9 @@ module Posit.Internal.PositC (PositC(..), - ES(..) + ES(..), + IntN, + FixedWidthInteger() ) where import Prelude hiding (exponent,significand) @@ -48,13 +52,14 @@ import Data.Int (Int8,Int16,Int32,Int64) -- Import standard Int sizes import Data.DoubleWord (Word128,Int128,Int256,fromHiAndLo,hiWord,loWord) -- Import large Int sizes import Data.Word (Word64) -import Data.Bits ((.|.), shiftL, shift, testBit, (.&.), shiftR) +import Data.Bits (Bits(..), (.|.), shiftL, shift, testBit, (.&.), shiftR) -- Import Naturals and Rationals {-@ embed Natural * as int @-} import GHC.Natural (Natural) -- Import the Natural Numbers ℕ (u+2115) -{-@ embed Ratio * as int @-} -import Data.Ratio (Rational, (%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D) +{-@ embed Ratio * as real @-} +{-@ embed Rational * as real @-} +import Data.Ratio ((%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D) -- | The Exponent Size 'ES' kind, the constructor for the Type is a Roman Numeral. @@ -65,43 +70,103 @@ | IV | V +-- | Type of the Finite Precision Representation, in our case Int8, +-- Int16, Int32, Int64, Int128, Int256. The 'es' of kind 'ES' will +-- determine a result of 'r' such that you can determine the 'es' by the +-- 'r' +{-@ embed IntN * as int @-} +type family IntN (es :: ES) = r | r -> es + where + IntN Z = Int8 + IntN I = Int16 + IntN II = Int32 + IntN III = Int64 + IntN IV = Int128 + IntN V = Int256 --- | The 'Posit' class is an approximation of ℝ, it is like a sampling on the Projective Real line ℙ(ℝ) with Maybe ℚ as the internal type. +-- | 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. +type FixedWidthInteger a = + (Bits a + ,Bounded a + ,Enum a + ,Integral a + ,Eq a + ,Ord a + ,Num a + ,Read a + ,Show a +#ifndef O_NO_STORABLE + ,Storable a +#endif + ) + + +-- | The 'Posit' class is an approximation of ℝ, it is like a sampling +-- on the Projective Real line ℙ(ℝ) with Maybe ℚ as the internal type. -- The 'es' is an index that controlls the log2 word size of the Posit's -- fininte precision representation. -class PositC (es :: ES) where - -- | Type of the Finite Precision Representation, in our case Int8, Int16, Int32, Int64, Int128, Int256. The 'es' of kind 'ES' will determine a result of 'r' such that you can determine the 'es' by the 'r' - type IntN es = r | r -> es - - +class (FixedWidthInteger (IntN es)) => PositC (es :: ES) where + -- | Transform to/from the Infinite Precision Representation encode :: Maybe Rational -> IntN es -- ^ Maybe you have some Rational Number and you want to encode it as some integer with a finite integer log2 word size. + encode Nothing = unReal @es + encode (Just 0) = 0 + encode (Just r) + | r > maxPosRat @es = mostPosVal @es + | r < minNegRat @es = mostNegVal @es + | r > 0 && r < minPosRat @es = leastPosVal @es + | r < 0 && r > maxNegRat @es = leastNegVal @es + | otherwise = buildIntRep @es r + decode :: IntN es -> Maybe Rational -- ^ You have an integer with a finite integer log2 word size decode it and Maybe it is Rational - + decode int + | int == unReal @es = Nothing + | int == 0 = Just 0 + | otherwise = + let sgn = int < 0 + int' = if sgn + then negate int + else int + (regime,nR) = regime2Integer @es int' + exponent = exponent2Nat @es nR int' -- if no e or some bits missing, then they are considered zero + rat = fraction2Posit @es nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 + in tupPosit2Posit @es (sgn,regime,exponent,rat) + + -- | Exponent Size based on the Posit Exponent kind ES exponentSize :: Natural -- ^ The exponent size, 'es' is a Natural number - + -- | Various other size definitions used in the Posit format with their default definitions nBytes :: Natural -- ^ 'nBytes' the number of bytes of the Posit Representation nBytes = 2^(exponentSize @es) - + nBits :: Natural -- ^ 'nBits' the number of bits of the Posit Representation nBits = 8 * (nBytes @es) - + signBitSize :: Natural -- ^ 'signBitSize' the size of the sign bit signBitSize = 1 - + uSeed :: Natural -- ^ 'uSeed' scaling factor for the regime of the Posit Representation uSeed = 2^(nBytes @es) - + -- | Integer Representation of common bounds unReal :: IntN es -- ^ 'unReal' is something that is not Real, the integer value that is not a Real number - + unReal = minBound @(IntN es) + mostPosVal :: IntN es + mostPosVal = maxBound @(IntN es) + leastPosVal :: IntN es + leastPosVal = 1 + leastNegVal :: IntN es + leastNegVal = -1 + mostNegVal :: IntN es - + mostNegVal = negate mostPosVal + -- Rational Value of common bounds maxPosRat :: Rational maxPosRat = 2^((nBytes @es) * ((nBits @es) - 2)) % 1 @@ -111,9 +176,9 @@ maxNegRat = negate (minPosRat @es) minNegRat :: Rational minNegRat = negate (maxPosRat @es) - + -- Functions to support encode and decode - + -- log base uSeed -- After calculating the regime the rational should be in the range [1,uSeed), it starts with (0,rational) log_uSeed :: (Integer, Rational) -> (Integer, Rational) @@ -121,91 +186,47 @@ | r < 1 = log_uSeed @es (regime-1,r * fromRational (toInteger (uSeed @es) % 1)) | r >= fromRational (toInteger (uSeed @es) % 1) = log_uSeed @es (regime+1,r * fromRational (1 % toInteger (uSeed @es))) | otherwise = (regime,r) - + getRegime :: Rational -> (Integer, Rational) getRegime r = log_uSeed @es (0,r) - + posit2TupPosit :: Rational -> (Bool, Integer, Natural, Rational) posit2TupPosit r = let (sgn,r') = getSign r -- returns the sign and a positive rational (regime,r'') = getRegime @es r' -- returns the regime and a rational between uSeed^-1 to uSeed^1 (exponent,significand) = getExponent r'' -- returns the exponent and a rational between [1,2), the significand in (sgn,regime,exponent,significand) - + buildIntRep :: Rational -> IntN es - mkIntRep :: Integer -> Natural -> Rational -> IntN es - formRegime :: Integer -> (IntN es, Integer) - formExponent :: Natural -> Integer -> (IntN es, Integer) - formFraction :: Rational -> Integer -> IntN es - - tupPosit2Posit :: (Bool,Integer,Natural,Rational) -> Maybe Rational - tupPosit2Posit (sgn,regime,exponent,rat) = -- s = isNeg posit == True - let pow2 = toRational (uSeed @es)^^regime * 2^exponent - scale = if sgn - then negate pow2 - else pow2 - in Just $ scale * rat - - regime2Integer :: IntN es -> (Integer, Int) - findRegimeFormat :: IntN es -> Bool - countRegimeBits :: Bool -> IntN es -> Int - exponent2Nat :: Int -> IntN es -> Natural - fraction2Posit :: Int -> IntN es -> Rational - - -- prints out the IntN es value in 0b... format - displayBin :: IntN es -> String - -- decimal Precision - decimalPrec :: Int - decimalPrec = fromIntegral $ 2 * (nBytes @es) + 1 - - - -instance PositC Z where - type IntN Z = Int8 - exponentSize = 0 - - -- Posit Integer Rep of various values - unReal = minBound @Int8 - - mostPosVal = maxBound @Int8 - leastPosVal = 1 - leastNegVal = -1 - mostNegVal = negate mostPosVal - - encode Nothing = unReal @Z - encode (Just 0) = 0 - encode (Just r) - | r > maxPosRat @Z = mostPosVal @Z - | r < minNegRat @Z = mostNegVal @Z - | r > 0 && r < minPosRat @Z = leastPosVal @Z - | r < 0 && r > maxNegRat @Z = leastNegVal @Z - | otherwise = buildIntRep @Z r - buildIntRep r = - let (signBit,regime,exponent,significand) = posit2TupPosit @Z r - intRep = mkIntRep @Z regime exponent significand + let (signBit,regime,exponent,significand) = posit2TupPosit @es r + intRep = mkIntRep @es regime exponent significand in if signBit then negate intRep else intRep - + + mkIntRep :: Integer -> Natural -> Rational -> IntN es mkIntRep regime exponent significand = - let (regime', offset) = formRegime @Z regime -- offset is the number of binary digits remaining after the regime is formed - (exponent', offset') = formExponent @Z exponent offset -- offset' is the number of binary digits remaining after the exponent is formed - fraction = formFraction @Z significand offset' + let (regime', offset) = formRegime @es regime -- offset is the number of binary digits remaining after the regime is formed + (exponent', offset') = formExponent @es exponent offset -- offset' is the number of binary digits remaining after the exponent is formed + fraction = formFraction @es significand offset' in regime' .|. exponent' .|. fraction - + + formRegime :: Integer -> (IntN es, Integer) formRegime power | 0 <= power = - let offset = (fromIntegral (nBits @Z - 1) - power - 1) + let offset = (fromIntegral (nBits @es - 1) - power - 1) in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) | otherwise = - let offset = (fromIntegral (nBits @Z - 1) - abs power - 1) + let offset = (fromIntegral (nBits @es - 1) - abs power - 1) in (1 `shiftL` fromInteger offset, offset) - + + formExponent :: Natural -> Integer -> (IntN es, Integer) formExponent power offset = - let offset' = offset - fromIntegral (exponentSize @Z) + let offset' = offset - fromIntegral (exponentSize @es) in (fromIntegral power `shift` fromInteger offset', offset') - + + formFraction :: Rational -> Integer -> IntN es formFraction r offset = let numFractionBits = offset fractionSize = 2^numFractionBits @@ -213,63 +234,65 @@ in if numFractionBits >= 1 then fromInteger normFraction else 0 - - decode int - | int == unReal @Z = Nothing - | int == 0 = Just 0 - | otherwise = - let sgn = int < 0 - int' = if sgn - then negate int - else int - (regime,nR) = regime2Integer @Z int' - exponent = exponent2Nat @Z nR int' -- if no e or some bits missing, then they are considered zero - rat = fraction2Posit @Z nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 - in tupPosit2Posit @Z (sgn,regime,exponent,rat) - + + tupPosit2Posit :: (Bool,Integer,Natural,Rational) -> Maybe Rational + tupPosit2Posit (sgn,regime,exponent,rat) = -- s = isNeg posit == True + let pow2 = toRational (uSeed @es)^^regime * 2^exponent + scale = if sgn + then negate pow2 + else pow2 + in Just $ scale * rat + + regime2Integer :: IntN es -> (Integer, Int) regime2Integer posit = - let regimeFormat = findRegimeFormat @Z posit - regimeCount = countRegimeBits @Z regimeFormat posit + let regimeFormat = findRegimeFormat @es posit + regimeCount = countRegimeBits @es regimeFormat posit regime = calcRegimeInt regimeFormat regimeCount in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime - + -- will return the format of the regime, either HI or LO; it could get refactored in the future -- True means a 1 is the first bit in the regime - findRegimeFormat posit = testBit posit (fromIntegral (nBits @Z) - 1 - fromIntegral (signBitSize @Z)) - - countRegimeBits format posit = go (fromIntegral (nBits @Z) - 1 - fromIntegral (signBitSize @Z)) 0 + findRegimeFormat :: IntN es -> Bool + findRegimeFormat posit = testBit posit (fromIntegral (nBits @es) - 1 - fromIntegral (signBitSize @es)) + + countRegimeBits :: Bool -> IntN es -> Int + countRegimeBits format posit = go (fromIntegral (nBits @es) - 1 - fromIntegral (signBitSize @es)) 0 where go (-1) acc = acc go index acc | xnor format (testBit posit index) = go (index - 1) (acc + 1) | otherwise = acc - + -- knowing the number of the regime bits, and the sign bit we can extract -- the exponent. We mask to the left of the exponent to remove the sign and regime, and -- then shift to the right to remove the fraction. + exponent2Nat :: Int -> IntN es -> Natural exponent2Nat numBitsRegime posit = - let bitsRemaining = fromIntegral (nBits @Z) - numBitsRegime - fromIntegral (signBitSize @Z) + let bitsRemaining = fromIntegral (nBits @es) - numBitsRegime - fromIntegral (signBitSize @es) signNRegimeMask = 2^bitsRemaining - 1 int = posit .&. signNRegimeMask - nBitsToTheRight = fromIntegral (nBits @Z) - numBitsRegime - fromIntegral (signBitSize @Z) - fromIntegral (exponentSize @Z) + nBitsToTheRight = fromIntegral (nBits @es) - numBitsRegime - fromIntegral (signBitSize @es) - fromIntegral (exponentSize @es) in if bitsRemaining <=0 then 0 else if nBitsToTheRight < 0 then fromIntegral $ int `shiftL` negate nBitsToTheRight else fromIntegral $ int `shiftR` nBitsToTheRight - + -- knowing the number of the regime bits, sign bit, and the number of the -- exponent bits we can extract the fraction. We mask to the left of the fraction to -- remove the sign, regime, and exponent. If there is no fraction then the value is 1. + fraction2Posit :: Int -> IntN es -> Rational fraction2Posit numBitsRegime posit = - let offset = fromIntegral $ (signBitSize @Z) + fromIntegral numBitsRegime + (exponentSize @Z) - fractionSize = fromIntegral (nBits @Z) - offset + let offset = fromIntegral $ (signBitSize @es) + fromIntegral numBitsRegime + (exponentSize @es) + fractionSize = fromIntegral (nBits @es) - offset fractionBits = posit .&. (2^fractionSize - 1) in if fractionSize >= 1 then (2^fractionSize + toInteger fractionBits) % 2^fractionSize else 1 % 1 - - displayBin int = "0b" ++ go (fromIntegral (nBits @Z) - 1) + + -- prints out the IntN es value in 0b... format + displayBin :: IntN es -> String + displayBin int = "0b" ++ go (fromIntegral (nBits @es) - 1) where go :: Int -> String go 0 = if testBit int 0 @@ -277,614 +300,44 @@ else "0" go idx = if testBit int idx then "1" ++ go (idx - 1) - else "0" ++ go (idx -1) + else "0" ++ go (idx - 1) + + -- decimal Precision + decimalPrec :: Int + decimalPrec = fromIntegral $ 2 * (nBytes @es) + 1 + + {-# MINIMAL exponentSize #-} +-- ===================================================================== +-- === PositC Instances === +-- ===================================================================== +instance PositC Z where + exponentSize = 0 + + instance PositC I where - type IntN I = Int16 exponentSize = 1 - - -- Posit Integer Rep of various values - unReal = minBound @Int16 - - mostPosVal = maxBound @Int16 - leastPosVal = 1 - leastNegVal = -1 - mostNegVal = negate mostPosVal - - encode Nothing = unReal @I - encode (Just 0) = 0 - encode (Just r) - | r > maxPosRat @I = mostPosVal @I - | r < minNegRat @I = mostNegVal @I - | r > 0 && r < minPosRat @I = leastPosVal @I - | r < 0 && r > maxNegRat @I = leastNegVal @I - | otherwise = buildIntRep @I r - - buildIntRep r = - let (signBit,regime,exponent,significand) = posit2TupPosit @I r - intRep = mkIntRep @I regime exponent significand - in if signBit - then negate intRep - else intRep - - mkIntRep regime exponent significand = - let (regime', offset) = formRegime @I regime -- offset is the number of binary digits remaining after the regime is formed - (exponent', offset') = formExponent @I exponent offset -- offset' is the number of binary digits remaining after the exponent is formed - fraction = formFraction @I significand offset' - in regime' .|. exponent' .|. fraction - - formRegime power - | 0 <= power = - let offset = (fromIntegral (nBits @I - 1) - power - 1) - in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) - | otherwise = - let offset = (fromIntegral (nBits @I - 1) - abs power - 1) - in (1 `shiftL` fromInteger offset, offset) - - formExponent power offset = - let offset' = offset - fromIntegral (exponentSize @I) - in (fromIntegral power `shift` fromInteger offset', offset') - - formFraction r offset = - let numFractionBits = offset - fractionSize = 2^numFractionBits - normFraction = round $ (r - 1) * fractionSize -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1) - in if numFractionBits >= 1 - then fromInteger normFraction - else 0 - - decode int - | int == unReal @I = Nothing - | int == 0 = Just 0 - | otherwise = - let sgn = int < 0 - int' = if sgn - then negate int - else int - (regime,nR) = regime2Integer @I int' - exponent = exponent2Nat @I nR int' -- if no e or some bits missing, then they are considered zero - rat = fraction2Posit @I nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 - in tupPosit2Posit @I (sgn,regime,exponent,rat) - - regime2Integer posit = - let regimeFormat = findRegimeFormat @I posit - regimeCount = countRegimeBits @I regimeFormat posit - regime = calcRegimeInt regimeFormat regimeCount - in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime - - -- will return the format of the regime, either HI or LO; it could get refactored in the future - -- True means a 1 is the first bit in the regime - findRegimeFormat posit = testBit posit (fromIntegral (nBits @I) - 1 - fromIntegral (signBitSize @I)) - - countRegimeBits format posit = go (fromIntegral (nBits @I) - 1 - fromIntegral (signBitSize @I)) 0 - where - go (-1) acc = acc - go index acc - | xnor format (testBit posit index) = go (index - 1) (acc + 1) - | otherwise = acc - - -- knowing the number of the regime bits, and the sign bit we can extract - -- the exponent. We mask to the left of the exponent to remove the sign and regime, and - -- then shift to the right to remove the fraction. - exponent2Nat numBitsRegime posit = - let bitsRemaining = fromIntegral (nBits @I) - numBitsRegime - fromIntegral (signBitSize @I) - signNRegimeMask = 2^bitsRemaining - 1 - int = posit .&. signNRegimeMask - nBitsToTheRight = fromIntegral (nBits @I) - numBitsRegime - fromIntegral (signBitSize @I) - fromIntegral (exponentSize @I) - in if bitsRemaining <=0 - then 0 - else if nBitsToTheRight < 0 - then fromIntegral $ int `shiftL` negate nBitsToTheRight - else fromIntegral $ int `shiftR` nBitsToTheRight - - -- knowing the number of the regime bits, sign bit, and the number of the - -- exponent bits we can extract the fraction. We mask to the left of the fraction to - -- remove the sign, regime, and exponent. If there is no fraction then the value is 1. - fraction2Posit numBitsRegime posit = - let offset = fromIntegral $ (signBitSize @I) + fromIntegral numBitsRegime + (exponentSize @I) - fractionSize = fromIntegral (nBits @I) - offset - fractionBits = posit .&. (2^fractionSize - 1) - in if fractionSize >= 1 - then (2^fractionSize + toInteger fractionBits) % 2^fractionSize - else 1 % 1 - - displayBin int = "0b" ++ go (fromIntegral (nBits @I) - 1) - where - go :: Int -> String - go 0 = if testBit int 0 - then "1" - else "0" - go idx = if testBit int idx - then "1" ++ go (idx - 1) - else "0" ++ go (idx -1) - instance PositC II where - type IntN II = Int32 exponentSize = 2 - - -- Posit Integer Rep of various values - unReal = minBound @Int32 - - mostPosVal = maxBound @Int32 - leastPosVal = 1 - leastNegVal = -1 - mostNegVal = negate mostPosVal - - encode Nothing = unReal @II - encode (Just 0) = 0 - encode (Just r) - | r > maxPosRat @II = mostPosVal @II - | r < minNegRat @II = mostNegVal @II - | r > 0 && r < minPosRat @II = leastPosVal @II - | r < 0 && r > maxNegRat @II = leastNegVal @II - | otherwise = buildIntRep @II r - - buildIntRep r = - let (signBit,regime,exponent,significand) = posit2TupPosit @II r - intRep = mkIntRep @II regime exponent significand - in if signBit - then negate intRep - else intRep - - mkIntRep regime exponent significand = - let (regime', offset) = formRegime @II regime -- offset is the number of binary digits remaining after the regime is formed - (exponent', offset') = formExponent @II exponent offset -- offset' is the number of binary digits remaining after the exponent is formed - fraction = formFraction @II significand offset' - in regime' .|. exponent' .|. fraction - - formRegime power - | 0 <= power = - let offset = (fromIntegral (nBits @II - 1) - power - 1) - in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) - | otherwise = - let offset = (fromIntegral (nBits @II - 1) - abs power - 1) - in (1 `shiftL` fromInteger offset, offset) - - formExponent power offset = - let offset' = offset - fromIntegral (exponentSize @II) - in (fromIntegral power `shift` fromInteger offset', offset') - - formFraction r offset = - let numFractionBits = offset - fractionSize = 2^numFractionBits - normFraction = round $ (r - 1) * fractionSize -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1) - in if numFractionBits >= 1 - then fromInteger normFraction - else 0 - - decode int - | int == unReal @II = Nothing - | int == 0 = Just 0 - | otherwise = - let sgn = int < 0 - int' = if sgn - then negate int - else int - (regime,nR) = regime2Integer @II int' - exponent = exponent2Nat @II nR int' -- if no e or some bits missing, then they are considered zero - rat = fraction2Posit @II nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 - in tupPosit2Posit @II (sgn,regime,exponent,rat) - - regime2Integer posit = - let regimeFormat = findRegimeFormat @II posit - regimeCount = countRegimeBits @II regimeFormat posit - regime = calcRegimeInt regimeFormat regimeCount - in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime - - -- will return the format of the regime, either HI or LO; it could get refactored in the future - -- True means a 1 is the first bit in the regime - findRegimeFormat posit = testBit posit (fromIntegral (nBits @II) - 1 - fromIntegral (signBitSize @II)) - - countRegimeBits format posit = go (fromIntegral (nBits @II) - 1 - fromIntegral (signBitSize @II)) 0 - where - go (-1) acc = acc - go index acc - | xnor format (testBit posit index) = go (index - 1) (acc + 1) - | otherwise = acc - - -- knowing the number of the regime bits, and the sign bit we can extract - -- the exponent. We mask to the left of the exponent to remove the sign and regime, and - -- then shift to the right to remove the fraction. - exponent2Nat numBitsRegime posit = - let bitsRemaining = fromIntegral (nBits @II) - numBitsRegime - fromIntegral (signBitSize @II) - signNRegimeMask = 2^bitsRemaining - 1 - int = posit .&. signNRegimeMask - nBitsToTheRight = fromIntegral (nBits @II) - numBitsRegime - fromIntegral (signBitSize @II) - fromIntegral (exponentSize @II) - in if bitsRemaining <=0 - then 0 - else if nBitsToTheRight < 0 - then fromIntegral $ int `shiftL` negate nBitsToTheRight - else fromIntegral $ int `shiftR` nBitsToTheRight - - -- knowing the number of the regime bits, sign bit, and the number of the - -- exponent bits we can extract the fraction. We mask to the left of the fraction to - -- remove the sign, regime, and exponent. If there is no fraction then the value is 1. - fraction2Posit numBitsRegime posit = - let offset = fromIntegral $ (signBitSize @II) + fromIntegral numBitsRegime + (exponentSize @II) - fractionSize = fromIntegral (nBits @II) - offset - fractionBits = posit .&. (2^fractionSize - 1) - in if fractionSize >= 1 - then (2^fractionSize + toInteger fractionBits) % 2^fractionSize - else 1 % 1 - - displayBin int = "0b" ++ go (fromIntegral (nBits @II) - 1) - where - go :: Int -> String - go 0 = if testBit int 0 - then "1" - else "0" - go idx = if testBit int idx - then "1" ++ go (idx - 1) - else "0" ++ go (idx -1) - instance PositC III where - type IntN III = Int64 exponentSize = 3 - - -- Posit Integer Rep of various values - unReal = minBound @Int64 - - mostPosVal = maxBound @Int64 - leastPosVal = 1 - leastNegVal = -1 - mostNegVal = negate mostPosVal - - encode Nothing = unReal @III - encode (Just 0) = 0 - encode (Just r) - | r > maxPosRat @III = mostPosVal @III - | r < minNegRat @III = mostNegVal @III - | r > 0 && r < minPosRat @III = leastPosVal @III - | r < 0 && r > maxNegRat @III = leastNegVal @III - | otherwise = buildIntRep @III r - - buildIntRep r = - let (signBit,regime,exponent,significand) = posit2TupPosit @III r - intRep = mkIntRep @III regime exponent significand - in if signBit - then negate intRep - else intRep - - mkIntRep regime exponent significand = - let (regime', offset) = formRegime @III regime -- offset is the number of binary digits remaining after the regime is formed - (exponent', offset') = formExponent @III exponent offset -- offset' is the number of binary digits remaining after the exponent is formed - fraction = formFraction @III significand offset' - in regime' .|. exponent' .|. fraction - - formRegime power - | 0 <= power = - let offset = (fromIntegral (nBits @III - 1) - power - 1) - in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) - | otherwise = - let offset = (fromIntegral (nBits @III - 1) - abs power - 1) - in (1 `shiftL` fromInteger offset, offset) - - formExponent power offset = - let offset' = offset - fromIntegral (exponentSize @III) - in (fromIntegral power `shift` fromInteger offset', offset') - - formFraction r offset = - let numFractionBits = offset - fractionSize = 2^numFractionBits - normFraction = round $ (r - 1) * fractionSize -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1) - in if numFractionBits >= 1 - then fromInteger normFraction - else 0 - - decode int - | int == unReal @III = Nothing - | int == 0 = Just 0 - | otherwise = - let sgn = int < 0 - int' = if sgn - then negate int - else int - (regime,nR) = regime2Integer @III int' - exponent = exponent2Nat @III nR int' -- if no e or some bits missing, then they are considered zero - rat = fraction2Posit @III nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 - in tupPosit2Posit @III (sgn,regime,exponent,rat) - - regime2Integer posit = - let regimeFormat = findRegimeFormat @III posit - regimeCount = countRegimeBits @III regimeFormat posit - regime = calcRegimeInt regimeFormat regimeCount - in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime - - -- will return the format of the regime, either HI or LO; it could get refactored in the future - -- True means a 1 is the first bit in the regime - findRegimeFormat posit = testBit posit (fromIntegral (nBits @III) - 1 - fromIntegral (signBitSize @III)) - - countRegimeBits format posit = go (fromIntegral (nBits @III) - 1 - fromIntegral (signBitSize @III)) 0 - where - go (-1) acc = acc - go index acc - | xnor format (testBit posit index) = go (index - 1) (acc + 1) - | otherwise = acc - - -- knowing the number of the regime bits, and the sign bit we can extract - -- the exponent. We mask to the left of the exponent to remove the sign and regime, and - -- then shift to the right to remove the fraction. - exponent2Nat numBitsRegime posit = - let bitsRemaining = fromIntegral (nBits @III) - numBitsRegime - fromIntegral (signBitSize @III) - signNRegimeMask = 2^bitsRemaining - 1 - int = posit .&. signNRegimeMask - nBitsToTheRight = fromIntegral (nBits @III) - numBitsRegime - fromIntegral (signBitSize @III) - fromIntegral (exponentSize @III) - in if bitsRemaining <=0 - then 0 - else if nBitsToTheRight < 0 - then fromIntegral $ int `shiftL` negate nBitsToTheRight - else fromIntegral $ int `shiftR` nBitsToTheRight - - -- knowing the number of the regime bits, sign bit, and the number of the - -- exponent bits we can extract the fraction. We mask to the left of the fraction to - -- remove the sign, regime, and exponent. If there is no fraction then the value is 1. - fraction2Posit numBitsRegime posit = - let offset = fromIntegral $ (signBitSize @III) + fromIntegral numBitsRegime + (exponentSize @III) - fractionSize = fromIntegral (nBits @III) - offset - fractionBits = posit .&. (2^fractionSize - 1) - in if fractionSize >= 1 - then (2^fractionSize + toInteger fractionBits) % 2^fractionSize - else 1 % 1 - - displayBin int = "0b" ++ go (fromIntegral (nBits @III) - 1) - where - go :: Int -> String - go 0 = if testBit int 0 - then "1" - else "0" - go idx = if testBit int idx - then "1" ++ go (idx - 1) - else "0" ++ go (idx -1) - instance PositC IV where - type IntN IV = Int128 exponentSize = 4 - - -- Posit Integer Rep of various values - unReal = minBound @Int128 - - mostPosVal = maxBound @Int128 - leastPosVal = 1 - leastNegVal = -1 - mostNegVal = negate mostPosVal - - encode Nothing = unReal @IV - encode (Just 0) = 0 - encode (Just r) - | r > maxPosRat @IV = mostPosVal @IV - | r < minNegRat @IV = mostNegVal @IV - | r > 0 && r < minPosRat @IV = leastPosVal @IV - | r < 0 && r > maxNegRat @IV = leastNegVal @IV - | otherwise = buildIntRep @IV r - - buildIntRep r = - let (signBit,regime,exponent,significand) = posit2TupPosit @IV r - intRep = mkIntRep @IV regime exponent significand - in if signBit - then negate intRep - else intRep - - mkIntRep regime exponent significand = - let (regime', offset) = formRegime @IV regime -- offset is the number of binary digits remaining after the regime is formed - (exponent', offset') = formExponent @IV exponent offset -- offset' is the number of binary digits remaining after the exponent is formed - fraction = formFraction @IV significand offset' - in regime' .|. exponent' .|. fraction - - formRegime power - | 0 <= power = - let offset = (fromIntegral (nBits @IV - 1) - power - 1) - in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) - | otherwise = - let offset = (fromIntegral (nBits @IV - 1) - abs power - 1) - in (1 `shiftL` fromInteger offset, offset) - - formExponent power offset = - let offset' = offset - fromIntegral (exponentSize @IV) - in (fromIntegral power `shift` fromInteger offset', offset') - - formFraction r offset = - let numFractionBits = offset - fractionSize = 2^numFractionBits - normFraction = round $ (r - 1) * fractionSize -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1) - in if numFractionBits >= 1 - then fromInteger normFraction - else 0 - - decode int - | int == unReal @IV = Nothing - | int == 0 = Just 0 - | otherwise = - let sgn = int < 0 - int' = if sgn - then negate int - else int - (regime,nR) = regime2Integer @IV int' - exponent = exponent2Nat @IV nR int' -- if no e or some bits missing, then they are considered zero - rat = fraction2Posit @IV nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 - in tupPosit2Posit @IV (sgn,regime,exponent,rat) - - regime2Integer posit = - let regimeFormat = findRegimeFormat @IV posit - regimeCount = countRegimeBits @IV regimeFormat posit - regime = calcRegimeInt regimeFormat regimeCount - in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime - - -- will return the format of the regime, either HI or LO; it could get refactored in the future - -- True means a 1 is the first bit in the regime - findRegimeFormat posit = testBit posit (fromIntegral (nBits @IV) - 1 - fromIntegral (signBitSize @IV)) - - countRegimeBits format posit = go (fromIntegral (nBits @IV) - 1 - fromIntegral (signBitSize @IV)) 0 - where - go (-1) acc = acc - go index acc - | xnor format (testBit posit index) = go (index - 1) (acc + 1) - | otherwise = acc - - -- knowing the number of the regime bits, and the sign bit we can extract - -- the exponent. We mask to the left of the exponent to remove the sign and regime, and - -- then shift to the right to remove the fraction. - exponent2Nat numBitsRegime posit = - let bitsRemaining = fromIntegral (nBits @IV) - numBitsRegime - fromIntegral (signBitSize @IV) - signNRegimeMask = 2^bitsRemaining - 1 - int = posit .&. signNRegimeMask - nBitsToTheRight = fromIntegral (nBits @IV) - numBitsRegime - fromIntegral (signBitSize @IV) - fromIntegral (exponentSize @IV) - in if bitsRemaining <=0 - then 0 - else if nBitsToTheRight < 0 - then fromIntegral $ int `shiftL` negate nBitsToTheRight - else fromIntegral $ int `shiftR` nBitsToTheRight - - -- knowing the number of the regime bits, sign bit, and the number of the - -- exponent bits we can extract the fraction. We mask to the left of the fraction to - -- remove the sign, regime, and exponent. If there is no fraction then the value is 1. - fraction2Posit numBitsRegime posit = - let offset = fromIntegral $ (signBitSize @IV) + fromIntegral numBitsRegime + (exponentSize @IV) - fractionSize = fromIntegral (nBits @IV) - offset - fractionBits = posit .&. (2^fractionSize - 1) - in if fractionSize >= 1 - then (2^fractionSize + toInteger fractionBits) % 2^fractionSize - else 1 % 1 - - displayBin int = "0b" ++ go (fromIntegral (nBits @IV) - 1) - where - go :: Int -> String - go 0 = if testBit int 0 - then "1" - else "0" - go idx = if testBit int idx - then "1" ++ go (idx - 1) - else "0" ++ go (idx -1) - instance PositC V where - type IntN V = Int256 exponentSize = 5 - - -- Posit Integer Rep of various values - unReal = minBound @Int256 - - mostPosVal = maxBound @Int256 - leastPosVal = 1 - leastNegVal = -1 - mostNegVal = negate mostPosVal - - encode Nothing = unReal @V - encode (Just 0) = 0 - encode (Just r) - | r > maxPosRat @V = mostPosVal @V - | r < minNegRat @V = mostNegVal @V - | r > 0 && r < minPosRat @V = leastPosVal @V - | r < 0 && r > maxNegRat @V = leastNegVal @V - | otherwise = buildIntRep @V r - - buildIntRep r = - let (signBit,regime,exponent,significand) = posit2TupPosit @V r - intRep = mkIntRep @V regime exponent significand - in if signBit - then negate intRep - else intRep - - mkIntRep regime exponent significand = - let (regime', offset) = formRegime @V regime -- offset is the number of binary digits remaining after the regime is formed - (exponent', offset') = formExponent @V exponent offset -- offset' is the number of binary digits remaining after the exponent is formed - fraction = formFraction @V significand offset' - in regime' .|. exponent' .|. fraction - - formRegime power - | 0 <= power = - let offset = (fromIntegral (nBits @V - 1) - power - 1) - in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) - | otherwise = - let offset = (fromIntegral (nBits @V - 1) - abs power - 1) - in (1 `shiftL` fromInteger offset, offset) - - formExponent power offset = - let offset' = offset - fromIntegral (exponentSize @V) - in (fromIntegral power `shift` fromInteger offset', offset') - - formFraction r offset = - let numFractionBits = offset - fractionSize = 2^numFractionBits - normFraction = round $ (r - 1) * fractionSize -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1) - in if numFractionBits >= 1 - then fromInteger normFraction - else 0 - - decode int - | int == unReal @V = Nothing - | int == 0 = Just 0 - | otherwise = - let sgn = int < 0 - int' = if sgn - then negate int - else int - (regime,nR) = regime2Integer @V int' - exponent = exponent2Nat @V nR int' -- if no e or some bits missing, then they are considered zero - rat = fraction2Posit @V nR int' -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1 - in tupPosit2Posit @V (sgn,regime,exponent,rat) - - regime2Integer posit = - let regimeFormat = findRegimeFormat @V posit - regimeCount = countRegimeBits @V regimeFormat posit - regime = calcRegimeInt regimeFormat regimeCount - in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime - - -- will return the format of the regime, either HI or LO; it could get refactored in the future - -- True means a 1 is the first bit in the regime - findRegimeFormat posit = testBit posit (fromIntegral (nBits @V) - 1 - fromIntegral (signBitSize @V)) - - countRegimeBits format posit = go (fromIntegral (nBits @V) - 1 - fromIntegral (signBitSize @V)) 0 - where - go (-1) acc = acc - go index acc - | xnor format (testBit posit index) = go (index - 1) (acc + 1) - | otherwise = acc - - -- knowing the number of the regime bits, and the sign bit we can extract - -- the exponent. We mask to the left of the exponent to remove the sign and regime, and - -- then shift to the right to remove the fraction. - exponent2Nat numBitsRegime posit = - let bitsRemaining = fromIntegral (nBits @V) - numBitsRegime - fromIntegral (signBitSize @V) - signNRegimeMask = 2^bitsRemaining - 1 - int = posit .&. signNRegimeMask - nBitsToTheRight = fromIntegral (nBits @V) - numBitsRegime - fromIntegral (signBitSize @V) - fromIntegral (exponentSize @V) - in if bitsRemaining <=0 - then 0 - else if nBitsToTheRight < 0 - then fromIntegral $ int `shiftL` negate nBitsToTheRight - else fromIntegral $ int `shiftR` nBitsToTheRight - - -- knowing the number of the regime bits, sign bit, and the number of the - -- exponent bits we can extract the fraction. We mask to the left of the fraction to - -- remove the sign, regime, and exponent. If there is no fraction then the value is 1. - fraction2Posit numBitsRegime posit = - let offset = fromIntegral $ (signBitSize @V) + fromIntegral numBitsRegime + (exponentSize @V) - fractionSize = fromIntegral (nBits @V) - offset - fractionBits = posit .&. (2^fractionSize - 1) - in if fractionSize >= 1 - then (2^fractionSize + toInteger fractionBits) % 2^fractionSize - else 1 % 1 - - displayBin int = "0b" ++ go (fromIntegral (nBits @V) - 1) - where - go :: Int -> String - go 0 = if testBit int 0 - then "1" - else "0" - go idx = if testBit int idx - then "1" ++ go (idx - 1) - else "0" ++ go (idx -1) + -- ===================================================================== -- === Encode and Decode Helpers === -- ===================================================================== @@ -900,7 +353,7 @@ else r in (s,absPosit) -- pretty much the same as 'abs') --- Exponent should be an integer in the range of [0,uSeed), and also return the posit [1,2) +-- Exponent should be an integer in the range of [0,uSeed), and also return an exponent and a rational in the range of [1,2) getExponent :: Rational -> (Natural, Rational) getExponent r = log_2 (0,r) @@ -930,17 +383,15 @@ sizeOf _ = 16 alignment _ = 16 peek ptr = do - hi <- peek $ offsetInt 0 + hi <- peek $ offsetWord 0 lo <- peek $ offsetWord 1 return $ fromHiAndLo hi lo where - offsetInt i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8) offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8) poke ptr int = do - poke (offsetInt 0) (hiWord int) + poke (offsetWord 0) (hiWord int) poke (offsetWord 1) (loWord int) where - offsetInt i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8) offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8) -- Orphan Instance for Int128 using the DoubleWord type class diff --git a/test/TestPosit.hs b/test/TestPosit.hs --- a/test/TestPosit.hs +++ b/test/TestPosit.hs @@ -19,6 +19,9 @@ main :: IO () main = do -- + print $ "exp(1)**(pi*sqrt 43): " ++ show (exp(1 :: Posit256) ** (pi * sqrt 43)) -- + print $ "exp(1)**(pi*sqrt 67): " ++ show (exp(1 :: Posit256) ** (pi * sqrt 67)) -- + print $ "exp(1)**(pi*sqrt 163): " ++ show (exp(1 :: Posit256) ** (pi * sqrt 163)) -- print $ "Machine Alpha Posit8 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit8)) -- succ (Posit int) = Posit (succ int) print $ "Machine Alpha Posit16 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit16)) -- print $ "Machine Alpha Posit32 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit32)) -- @@ -31,6 +34,16 @@ let sqrtTuma = (funLogDomainReduction funLogTuma).(/2).(funExp2 funExpTuma).(/log 2) print $ "sqrt phi using a Tuma algorithm: " ++ show (sqrtTuma phi) print $ "Tuma is fasta: " ++ show (sqrtTaylor (1/1000000) - sqrtTuma (1/1000000)) + let truth = 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: Posit256 + eval "Standard: gamma(phi): " (gamma (phi)) truth + eval "Fused Gamma: gamma(phi): " (funGammaSeriesFused (phi)) truth + eval "Ramanujan Gamma: gamma(phi): " (funGammaRamanujan (phi)) truth + eval "Calc Gamma: gamma(phi): " (funGammaCalc (phi)) truth + eval "Nemes Gamma: gamma(phi): " (funGammaNemes (phi)) truth + eval "Yang Gamma: gamma(phi): " (funGammaYang (phi)) truth + eval "Chen Gamma: gamma(phi): " (funGammaChen (phi)) truth + eval "Gamma (x - 1): gamma(phi): " (funGammaXminus1 (phi)) truth + eval "Wolfram alpha: gamma(phi): " truth truth let truth = 5.0431656433600286513118821892854247103235901754138463603020001967777869609108929428415187821843384653305404495551887666992776792 :: Posit256 eval "Standard: exp(phi):" (exp (phi)) truth eval "Taylor: exp(phi):" (funExp2 funExpTaylor (phi / log 2)) truth