diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,9 @@ +# Changelog for Posit Numbers + +## posit-3.2.0.0 + + * Posit Standard 3.2 [Posit Standard] (https://posithub.org/docs/posit_standard.pdf) + * LiquidHaskell support: stack build --flag posit:do-liquid + * To run the test suite: stack test --flag posit:do-test + * To play around: stack repl --flag posit:do-test + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Nathan Waivio (c) 2021-2022 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Nathan Waivio nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,91 @@ +# posit 3.2.0.0 + +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: + + * Show + * Eq + * Ord -- compare as an integer representation + * Num -- Addition, subtraction, multiplication, and other operations + * Enum -- Successor and Predecessor + * Fractional -- division, divide by zero is Not a Real (NaR) number + * Real + * Bounded + * FusedOps -- dot product and others + * Convertable -- Conversions between different posit formats + * AltShow + * Read + * Storable -- Formats for binary data, for computation and data interchange + * RealFrac + * RealFloat + * 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, +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 +couple of auxiliary classes, like AltShow, AltFloating, and FusedOps. + +``` +class AltShow a where + -- Display the Posit in its Binary Representation + displayBinary :: a -> String + -- Display the Posit in its Integral Representation + displayIntegral :: a -> String + -- Display the Posit as a Rational + displayRational :: a -> String + -- Display the Posit as a Decimal until the Repented occurs + displayDecimal :: a -> String +``` + +``` +class AltFloating p where + phi :: p + gamma :: p -> p + sinc :: p -> p + expm1 :: p -> p +``` + +``` +class Num a => FusedOps a where + -- |Fused Multiply Add: (a * b) + c + fma :: a -> a -> a -> a + -- |Fused Add Multiply: (a + b) * c + fam :: a -> a -> a -> a + -- |Fused Multiply Multiply Subtract: (a * b) - (c * d) + fmms :: a -> a -> a -> a -> a + -- |Fused Sum of 3 values: a + b + c + fsum3 :: a -> a -> a -> a + -- |Fused Sum of 4 values: a + b + c + d + fsum4 :: a -> a -> a -> a -> a + -- |Fused Sum of a List of Posits + fsumL :: Foldable t => t a -> a + -- |Fused Dot Product of 3 element vector: (a1 * b1) + (a2 * b2) + (a3 * b3) + fdot3 :: a -> a -> a -> a -> a -> a -> a + -- |Fused Dot Product of 4 element veector: (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3) + fdot4 :: a -> a -> a -> a -> a -> a -> a -> a -> a + -- |Fused Dot Product of Two Lists + fdotL :: Foldable t => t a -> t a -> a + -- |Fused Subtract Multiply: a - (b * c) + fsm :: a -> a -> a -> a +``` + +The Posit type is 'Convertible' between other Posit lengths. + +``` +class Convertible a b where + convert :: a -> b +``` + +The Posit Library is built on top of two of the most excellent libraries: +[data-dword](https://hackage.haskell.org/package/data-dword), and +[scientific](https://hackage.haskell.org/package/scientific). The +'data-dword' library provides the underlining machine word +representation, it can provide 2^es word size, 2's complement fixed +length integers. The 'scientific' library provides 'read' and 'show' +instances. + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/posit.cabal b/posit.cabal new file mode 100644 --- /dev/null +++ b/posit.cabal @@ -0,0 +1,96 @@ +cabal-version: 1.12 + +name: posit +version: 3.2.0.0 +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 +author: Nathan Waivio +maintainer: nathan.waivio@gmail.com +copyright: 2021-2022 Nathan Waivio +license: BSD3 +license-file: LICENSE +build-type: Simple +extra-source-files: + README.md + ChangeLog.md + +source-repository head + type: git + location: https://github.com/waivio/posit + +flag do-no-storable + description: Build without Storable Class support + manual: True + default: False + +flag do-no-orphans + description: Build without Orphan Instances if data-dword gets updated for Storable Instances + manual: True + default: False + +flag do-liquid + description: Build with Liquid Haskell checking + manual: True + default: False + +flag do-test + description: Export additional algorithms for calculating primitive functions for test purposes + manual: True + default: False + +library + exposed-modules: + Posit + Posit.Internal.PositC + other-modules: + hs-source-dirs: + src + build-depends: + data-dword, + scientific + default-language: Haskell2010 + + -- Compiler options + 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 + + if flag(do-no-storable) + cpp-options: -DO_NO_STORABLE + + if flag(do-no-orphans) + cpp-options: -DO_NO_ORPHANS + + if flag(do-liquid) + cpp-options: -DO_LIQUID -DO_NO_STORABLE -DO_NO_READ -DO_NO_SHOW + + if flag(do-test) + cpp-options: -DO_TEST + + -- Other library packages from which modules are imported. + build-depends: + deepseq >=1.1 && <2 + + if !flag(do-liquid) + build-depends: + base >=4.7 && <5 + + if flag(do-liquid) + build-depends: + liquid-base, + liquidhaskell >= 0.8.10 + +-- perhaps one day: -threaded -rtsopts -with-rtsopts=-N +test-suite posit-test + type: exitcode-stdio-1.0 + main-is: TestPosit.hs + hs-source-dirs: + test + ghc-options: -O2 + cpp-options: -DO_TEST + build-depends: + base >=4.7 && <5 + , posit + default-language: Haskell2010 diff --git a/src/Posit.hs b/src/Posit.hs new file mode 100644 --- /dev/null +++ b/src/Posit.hs @@ -0,0 +1,1368 @@ + +-------------------------------------------------------------------------------------------- +-- Posit Numbers +-- Copyright : (C) 2022 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- Stability : Stable +-- Portability : Portable +-- +-- | Library implementing standard Posit Numbers (Posit Standard version +-- 3.2.0.0, with some improvements) a fixed width word size of +-- 2^es bytes. +-- +--------------------------------------------------------------------------------------------- + + +{-# LANGUAGE GADTs #-} -- For our main type Posit (es :: ES) +{-# LANGUAGE DataKinds #-} -- For our ES kind and the constructors Z, I, II, III, IV, V for exponent size type +{-# LANGUAGE KindSignatures #-} -- For defining the type of kind ES that indexes the GADT +{-# LANGUAGE ViewPatterns #-} -- To decode the posit in the pattern +{-# 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 +{-# LANGUAGE UndecidableInstances #-} -- To reduce some code duplication, I think the code is decidable but GHC is not smart enough ;), like there being only 1 instance that is polymorphic and works for all of my types. +{-# 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 +{-# OPTIONS_GHC -Wno-unused-top-binds #-} -- Turn off noise + +-- ---- +-- Posit numbers implementing: +-- +-- * Show +-- * Eq +-- * Ord -- compare as an integer representation +-- * Num -- Addition, subtraction, multiplication, and other operations +-- * Enum -- Successor and Predecessor +-- * Fractional -- division, divide by zero is Not a Real (NaR) number +-- * Real +-- * Bounded +-- * FusedOps -- dot product and others +-- * Convertible -- Conversions between different posit formats +-- * AltShow +-- * Read +-- * Storable -- Formats for binary data, for computation and data interchange +-- * RealFrac +-- * RealFloat +-- * Floating -- Mathematical functions such as logarithm, exponential, trigonometric, and hyperbolic functions. Warning! May induce trance. +-- +-- ---- + +module 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' + Posit64, -- |An 64-bit Posit number with 'es' ~ 'III' + Posit128, -- |An 128-bit Posit number with 'es' ~ 'IV' + Posit256, -- |An 256-bit Posit number with 'es' ~ 'V' + + -- * Patterns for Matching Exported Types + pattern NaR, -- |A pattern for Exception handling when a value is Not a Real number (NaR). + pattern R, -- |A pattern for the non-Exceptional case, yielding a Rational, will make a total function when paired with NaR, if the Rational implementation is total. + + -- * Fused Operation Interface defined by the Posit Standard + FusedOps(..), + + -- * Posits are Convertable between different Posit representations + Convertible(..), + + -- * Additional functions to show the Posit in different formats + AltShow(..), + + -- * Additional Special Functions + AltFloating(..), + + -- * Functions to lift functions of Integers or Rationals to operate on Posit Types + viaIntegral, + viaRational, + viaRational2, + viaRational3, + viaRational4, + viaRational6, + viaRational8, +#ifdef O_TEST + -- * Alternative algorithms for test purposes + funExp, + funExp2, + funExpTaylor, + funLogTaylor, + funExpTuma, + funLogTuma, + funLogDomainReduction, + funPi1, + funPi2, + funPi3, + funPi4, + funPsiSha1, + funPsiSha2, + funPsiSha3 +#endif + ) where + + +import Prelude hiding (rem) + +-- Imports for Show and Read Instances +import Data.Scientific (Scientific + ,scientificP + ,fromRationalRepetendUnlimited + ,formatScientific + ,FPFormat(Generic)) -- Used to print/show and read the rational value + +import Text.Read (Lexeme(Ident) + ,readPrec + ,readListPrec + ,(+++) + ,pfail + ,readListPrecDefault + ,lexP + ,lift + ,parens) -- Used to read a Posit value + +-- Imports for Vectorization Class Instances +import Data.Foldable (toList) -- Used for fused operations on foldable/lists + +-- 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 + + +-- would like to: +-- import Posit.Internal.ElementaryFunctions +-- 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 Debug.Trace (trace) -- temporary for debug purposes + + +-- ===================================================================== +-- === Posit Implementation === +-- ===================================================================== + +-- The machine implementation of the Posit encoding/decoding +import Posit.Internal.PositC (ES(..), PositC(..)) -- The main internal implementation details + + +-- |Base GADT rapper type, that uses the Exponent Size kind to index the various implementations +data Posit (es :: ES) where + 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 <- (Posit (decode -> Nothing)) where + NaR = Posit unReal +-- + +-- +-- |A Real or at least Rational Number, rounded to the nearest Posit Rational representation +pattern R :: PositC es => Rational -> Posit es +pattern R r <- (Posit (decode -> Just r)) where + R r = Posit (encode $ Just r) +-- + +-- Posit functions are complete if the following two patterns are completely defined. +{-# COMPLETE NaR, R #-} + +-- Concrete types exported for use. +type Posit8 = Posit Z +type Posit16 = Posit I +type Posit32 = Posit II +type Posit64 = Posit III +type Posit128 = Posit IV +type Posit256 = Posit V + +#ifndef O_NO_SHOW +-- Show +-- +instance forall es. (PositC es) => Show (Posit es) where + show NaR = "NaR" + show (R r) = formatScientific Generic (Just $ decimalPrec @es) (fst.fromRationalRepetendUnlimited $ r) +-- +#endif + + + +-- 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 + (Posit int1) == (Posit int2) = int1 == int2 +-- + + + +-- 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 + compare (Posit int1) (Posit int2) = compare int1 int2 +-- + + + +-- Num +-- +-- I'm num trying to get this definition: +instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Num (Posit es) where + -- Addition + (+) = viaRational2 (+) + -- Multiplication + (*) = viaRational2 (*) + -- '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 + -- 'signum' it is a kind of an representation of directionality, the sign of a number for instance + signum = viaRational signum + -- 'fromInteger' rounds the integer into the closest posit number + fromInteger int = Posit $ encode (Just $ 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 +-- + +-- 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 f (Posit int) = Posit $ f int +-- + + + +-- 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 + -- 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 + -- pred (Posit int) = Posit (int - 1) + pred = viaIntegral (subtract 1) + -- pred = viaIntegral pred -- 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 + -- enumFrom :: Posit es -> [Posit es] + enumFrom n = enumFromTo n maxBound + enumFromTo n m + | n == m = [n] + | n < m = n : enumFromTo (succ n) m + | otherwise = [] + -- enumFromThen n m :: Posit es -> Posit es -> [Posit es] + enumFromThen NaR _ = [NaR] + enumFromThen _ NaR = [NaR] + enumFromThen n m = n : go n + where + step = m - n + go :: Posit es -> [Posit es] + go NaR = [NaR] + go !l = case compare step 0 of + LT -> let !n' = l + step -- rounding occurs here, because the next comparison needs it, it wouldn't make sense otherwise... + in if n' - l > step + then [] + else n' : go n' + EQ -> [n, m] + GT -> let !n' = l + step + in if n' - l < step + then [] -- with tapered resolution this algorithm can reach a fixed point where the next value is equal to the previous value + else n' : go n' + enumFromThenTo NaR _ _ = [NaR] + enumFromThenTo _ NaR _ = [NaR] + enumFromThenTo _ _ NaR = [NaR] + enumFromThenTo e1 e2 e3 = takeWhile predicate (enumFromThen e1 e2) + where + mid = (e2 - e1) / 2 + predicate | e2 >= e1 = (<= e3 + mid) + | otherwise = (>= e3 + mid) +-- + + + +-- 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 + fromRational = R + + recip 0 = NaR + recip p = viaRational recip p +-- + +-- 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 + toRational NaR = error "Your input is Not a Real or Rational (NaR) number, please try again!" + toRational (R r) = r +-- + +-- 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 _ 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 _ 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 _ 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 _ NaR _ _ _ = NaR +viaRational4 _ _ NaR _ _ = NaR +viaRational4 _ _ _ NaR _ = NaR +viaRational4 _ _ _ _ NaR = NaR +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 _ NaR _ _ _ _ _ = NaR +viaRational6 _ _ NaR _ _ _ _ = NaR +viaRational6 _ _ _ NaR _ _ _ = NaR +viaRational6 _ _ _ _ NaR _ _ = NaR +viaRational6 _ _ _ _ _ NaR _ = NaR +viaRational6 _ _ _ _ _ _ NaR = NaR +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 _ NaR _ _ _ _ _ _ _ = NaR +viaRational8 _ _ NaR _ _ _ _ _ _ = NaR +viaRational8 _ _ _ NaR _ _ _ _ _ = NaR +viaRational8 _ _ _ _ NaR _ _ _ _ = NaR +viaRational8 _ _ _ _ _ NaR _ _ _ = NaR +viaRational8 _ _ _ _ _ _ NaR _ _ = NaR +viaRational8 _ _ _ _ _ _ _ NaR _ = NaR +viaRational8 _ _ _ _ _ _ _ _ NaR = NaR +viaRational8 f (R a0) (R a1) (R a2) (R a3) (R b0) (R b1) (R b2) (R b3) = R $ f a0 a1 a2 a3 b0 b1 b2 b3 + + + +-- 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 + -- 'minBound' the most negative number represented + minBound = Posit mostNegVal + -- 'maxBound' the most positive number represented + maxBound = Posit mostPosVal +-- + + +-- ===================================================================== +-- === Fused Operations === +-- ===================================================================== + +-- |A class that delays the rounding operation until the end for some operations +class Num a => FusedOps a where + -- |Fused Multiply Add: (a * b) + c + fma :: a -> a -> a -> a + -- |Fused Add Multiply: (a + b) * c + fam :: a -> a -> a -> a + -- |Fused Multiply Multiply Subtract: (a * b) - (c * d) + fmms :: a -> a -> a -> a -> a + -- |Fused Sum of 3 values: a + b + c + fsum3 :: a -> a -> a -> a + -- |Fused Sum of 4 values: a + b + c + d + fsum4 :: a -> a -> a -> a -> a + -- |Fused Sum of a List of Posits + fsumL :: Foldable t => t a -> a + -- |Fused Dot Product of 3 element vector: (a1 * b1) + (a2 * b2) + (a3 * b3) + fdot3 :: a -> a -> a -> a -> a -> a -> a + -- |Fused Dot Product of 4 element vector: (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3) + fdot4 :: a -> a -> a -> a -> a -> a -> a -> a -> a + -- |Fused Dot Product of Two Lists + fdotL :: Foldable t => t a -> t a -> a + -- |Fused Subtract Multiply: a - (b * c) + fsm :: a -> a -> a -> a + + + +-- Rational Instance +instance FusedOps Rational where + fsm a b c = a - (b * c) + fma a b c = (a * b) + c + fam a b c = (a + b) * c + fmms a b c d = (a * b) - (c * d) + fsum3 a b c = a + b + c + fsum4 a b c d = a + b + c + d + fsumL (toList -> l) = go l 0 + where + go [] acc = acc + go (x : xs) acc = go xs (acc + x) + fdot3 a1 a2 a3 b1 b2 b3 = (a1 * b1) + (a2 * b2) + (a3 * b3) + fdot4 a0 a1 a2 a3 b0 b1 b2 b3 = (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3) + fdotL (toList -> l1) (toList -> l2) = go l1 l2 0 + where + go [] [] acc = acc + go [] _ _ = error "Lists not the same length" + go _ [] _ = error "Lists not the same length" + go (b : bs) (c : cs) acc = go bs cs (fma b c acc) +-- + +-- +instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => FusedOps (Posit es) where + -- Fused Subtract Multiply + fsm = viaRational3 fsm + -- Fuse Multiply Add + fma = viaRational3 fma + -- Fuse Add Multiply + fam = viaRational3 fam + -- Fuse Multiply Multiply Subtract + fmms = viaRational4 fmms + -- Fuse Sum of 3 Posits + fsum3 = viaRational3 fsum3 + -- Fuse Sum of 4 Posits + fsum4 = viaRational4 fsum4 + -- Fuse Sum of a List + fsumL (toList -> l) = Posit $ encode (Just $ go l 0) + where + go :: [Posit es] -> Rational -> Rational + go [] !acc = acc + go ((Posit int) : xs) !acc = case decode int of + Nothing -> error "Posit List contains NaR" + Just r -> go xs (acc + r) + -- Fuse Dot Product of a 3-Vector + fdot3 = viaRational6 fdot3 + -- Fuse Dot Product of a 4-Vector + fdot4 = viaRational8 fdot4 + -- Fuse Dot Product of two Lists + fdotL (toList -> l1) (toList -> l2) = Posit $ encode (Just $ go l1 l2 0) + where + go [] [] !acc = acc + go [] _ _ = error "Lists not the same length" + go _ [] _ = error "Lists not the same length" + go ((Posit int1) : bs) ((Posit int2) : cs) !acc = case decode int1 of + Nothing -> error "First Posit List contains NaR" + Just r1 -> case decode int2 of + Nothing -> error "Second Posit List contains NaR" + Just r2 -> go bs cs (acc + (r1 * r2)) +-- + + + + +-- ===================================================================== +-- === Conversion Between Posits Types === +-- ===================================================================== + +-- |A Convertible class that will cast or 'convert' between two different Posit es types +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 + convert NaR = NaR + convert (R r) = R r +-- + + +#ifndef O_NO_SHOW +-- ===================================================================== +-- === Alternative Show Formats === +-- ===================================================================== + +-- |A Alternative to the typical 'Show' class to assist in displaying the Posit es type in different formats +class AltShow a where + -- |Display the Posit in its Binary Representation + displayBinary :: a -> String + -- |Display the Posit in its Integral Representation + displayIntegral :: a -> String + -- |Display the Posit as a Rational + displayRational :: a -> String + -- |Display the Posit as a Decimal until the Repetend occurs + displayDecimal :: a -> String +-- + +-- +instance forall es. (Show (IntN es), Ord (IntN es), Num (IntN es), PositC es) => AltShow (Posit es) where + displayBinary (Posit int) = displayBin int + + displayIntegral (Posit int) = show int + + displayRational = viaShowable id + + displayDecimal = viaShowable (fst.fromRationalRepetendUnlimited) +-- + +viaShowable :: forall es a. (Show a, Ord (IntN es), Num (IntN es), PositC es) => (Rational -> a) -> Posit es -> String +viaShowable _ NaR = "NaR" +viaShowable f (R r) = show $ f r +#endif + +#ifndef O_NO_READ +-- ===================================================================== +-- === Read Posit === +-- ===================================================================== + +-- +instance forall es. (PositC es) => Read (Posit es) where + readPrec = + parens $ do + x <- lexP + case x of + Ident "NaR" -> return NaR + _ -> pfail + +++ + do + s <- lift scientificP + return $ R (toRational s) + + readListPrec = readListPrecDefault +-- +#endif + + +-- ===================================================================== +-- === Storable Instances === +-- ===================================================================== +-- +#ifndef O_NO_STORABLE +-- +instance forall es. (Storable (IntN es), PositC es) => Storable (Posit es) where + sizeOf _ = fromIntegral $ nBytes @es + alignment _ = fromIntegral $ nBytes @es + peek ptr = do + int <- peek (castPtr ptr :: Ptr (IntN es)) + return $ Posit int + poke ptr (Posit int) = do + poke (castPtr ptr :: Ptr (IntN es)) int +-- +#endif + + +-- ===================================================================== +-- === Real Frac === +-- ===================================================================== + +-- +instance forall es. (Num (IntN es), Ord (IntN es), 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 err _ NaR = error err +viaRationalErrTrunkation _ f (R r) = + let (int, r') = f r + in (int, R r') + +-- ===================================================================== +-- === Real Float === +-- ===================================================================== +-- +instance forall es. (Eq (IntN es), Ord (IntN es), Num (IntN es), Floating (Posit es), PositC es) => RealFloat (Posit es) where + isIEEE _ = False + isDenormalized _ = False + isNegativeZero _ = False + + isNaN NaR = True + isNaN _ = False + + isInfinite NaR = True + isInfinite _ = False + + -- 'atan2' of y x is the argument "arg function" (also called phase or angle) of the complex number x + i y. + -- angle from an x basis vector to some other vector + -- + -- Y + -- ^ + -- | ^ (x,y) + -- | / + -- | / <- alpha (radians) + -- | / \ + -- / | + -- -----------------------------------> X + -- + -- + atan2 NaR _ = NaR + atan2 _ NaR = NaR + atan2 y x + | x == 0 && y == 0 = NaR + | x > 0 = atan (y/x) + | x < 0 && y >= 0 = atan (y/x) + pi + | x < 0 && y < 0 = atan (y/x) - pi + | x == 0 && y > 0 = pi / 2 + | x == 0 && y < 0 = negate $ pi / 2 + | otherwise = error "What!?!?!" -- The case where x == 0 && y == 0 + + floatRadix _ = 2 + floatDigits _ = undefined + floatRange _ = (negate maxExponent, maxExponent) + where + maxExponent = fromIntegral $ (nBytes @es) * ((nBits @es) - 2) + decodeFloat = undefined + encodeFloat = undefined +-- + + + +-- ===================================================================== +-- === Floating === +-- ===================================================================== + + +instance Floating Posit8 where + pi = convert (pi :: Posit256) :: Posit8 + exp x = convert (exp (convert x) :: Posit256) :: Posit8 + log x = convert (log (convert x) :: Posit256) :: Posit8 + x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit8 + sin x = convert (sin (convert x) :: Posit256) :: Posit8 + cos x = convert (cos (convert x) :: Posit256) :: Posit8 + asin x = convert (asin (convert x) :: Posit256) :: Posit8 + acos x = convert (acos (convert x) :: Posit256) :: Posit8 + atan x = convert (atan (convert x) :: Posit256) :: Posit8 + sinh x = convert (sinh (convert x) :: Posit256) :: Posit8 + cosh x = convert (cosh (convert x) :: Posit256) :: Posit8 + asinh x = convert (asinh (convert x) :: Posit256) :: Posit8 + acosh x = convert (acosh (convert x) :: Posit256) :: Posit8 + atanh x = convert (atanh (convert x) :: Posit256) :: Posit8 + +instance Floating Posit16 where + pi = convert (pi :: Posit256) :: Posit16 + exp x = convert (exp (convert x) :: Posit256) :: Posit16 + log x = convert (log (convert x) :: Posit256) :: Posit16 + x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit16 + sin x = convert (sin (convert x) :: Posit256) :: Posit16 + cos x = convert (cos (convert x) :: Posit256) :: Posit16 + asin x = convert (asin (convert x) :: Posit256) :: Posit16 + acos x = convert (acos (convert x) :: Posit256) :: Posit16 + atan x = convert (atan (convert x) :: Posit256) :: Posit16 + sinh x = convert (sinh (convert x) :: Posit256) :: Posit16 + cosh x = convert (cosh (convert x) :: Posit256) :: Posit16 + asinh x = convert (asinh (convert x) :: Posit256) :: Posit16 + acosh x = convert (acosh (convert x) :: Posit256) :: Posit16 + atanh x = convert (atanh (convert x) :: Posit256) :: Posit16 + +instance Floating Posit32 where + pi = convert (pi :: Posit256) :: Posit32 + exp x = convert (exp (convert x) :: Posit256) :: Posit32 + log x = convert (log (convert x) :: Posit256) :: Posit32 + x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit32 + sin x = convert (sin (convert x) :: Posit256) :: Posit32 + cos x = convert (cos (convert x) :: Posit256) :: Posit32 + asin x = convert (asin (convert x) :: Posit256) :: Posit32 + acos x = convert (acos (convert x) :: Posit256) :: Posit32 + atan x = convert (atan (convert x) :: Posit256) :: Posit32 + sinh x = convert (sinh (convert x) :: Posit256) :: Posit32 + cosh x = convert (cosh (convert x) :: Posit256) :: Posit32 + asinh x = convert (asinh (convert x) :: Posit256) :: Posit32 + acosh x = convert (acosh (convert x) :: Posit256) :: Posit32 + atanh x = convert (atanh (convert x) :: Posit256) :: Posit32 + +instance Floating Posit64 where + pi = convert (pi :: Posit256) :: Posit64 + exp x = convert (exp (convert x) :: Posit256) :: Posit64 + log x = convert (log (convert x) :: Posit256) :: Posit64 + x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit64 + sin x = convert (sin (convert x) :: Posit256) :: Posit64 + cos x = convert (cos (convert x) :: Posit256) :: Posit64 + asin x = convert (asin (convert x) :: Posit256) :: Posit64 + acos x = convert (acos (convert x) :: Posit256) :: Posit64 + atan x = convert (atan (convert x) :: Posit256) :: Posit64 + sinh x = convert (sinh (convert x) :: Posit256) :: Posit64 + cosh x = convert (cosh (convert x) :: Posit256) :: Posit64 + asinh x = convert (asinh (convert x) :: Posit256) :: Posit64 + acosh x = convert (acosh (convert x) :: Posit256) :: Posit64 + atanh x = convert (atanh (convert x) :: Posit256) :: Posit64 + +instance Floating Posit128 where + pi = convert (pi :: Posit256) :: Posit128 + exp x = convert (exp (convert x) :: Posit256) :: Posit128 + log x = convert (log (convert x) :: Posit256) :: Posit128 + x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit128 + sin x = convert (sin (convert x) :: Posit256) :: Posit128 + cos x = convert (cos (convert x) :: Posit256) :: Posit128 + asin x = convert (asin (convert x) :: Posit256) :: Posit128 + acos x = convert (acos (convert x) :: Posit256) :: Posit128 + atan x = convert (atan (convert x) :: Posit256) :: Posit128 + sinh x = convert (sinh (convert x) :: Posit256) :: Posit128 + cosh x = convert (cosh (convert x) :: Posit256) :: Posit128 + asinh x = convert (asinh (convert x) :: Posit256) :: Posit128 + acosh x = convert (acosh (convert x) :: Posit256) :: Posit128 + atanh x = convert (atanh (convert x) :: Posit256) :: Posit128 + +instance Floating Posit256 where + pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286 :: Posit256 + exp = funExp + log = funLogDomainReduction funLogTaylor + (**) = funPow + sin = funSin + cos = funCos + asin = funAsin + acos = funAcos + atan = funAtan + sinh = funSinh + cosh = funCosh + asinh = funAsinh + acosh = funAcosh + atanh = funAtanh + + + + + +class AltFloating p where + phi :: p + gamma :: p -> p + sinc :: p -> p + expm1 :: p -> p + +instance AltFloating Posit8 where + phi = convert (phi :: Posit256) :: Posit8 + gamma x = convert (gamma (convert x) :: Posit256) :: Posit8 + sinc x = convert (sinc (convert x) :: Posit256) :: Posit8 + expm1 x = + let b = atanh $ x / 2 + in (2 * b) / (1 - b) + +instance AltFloating Posit16 where + phi = convert (phi :: Posit256) :: Posit16 + gamma x = convert (gamma (convert x) :: Posit256) :: Posit16 + sinc x = convert (sinc (convert x) :: Posit256) :: Posit16 + expm1 x = + let b = atanh $ x / 2 + in (2 * b) / (1 - b) + +instance AltFloating Posit32 where + phi = convert (phi :: Posit256) :: Posit32 + gamma x = convert (gamma (convert x) :: Posit256) :: Posit32 + sinc x = convert (sinc (convert x) :: Posit256) :: Posit32 + expm1 x = + let b = atanh $ x / 2 + in (2 * b) / (1 - b) + +instance AltFloating Posit64 where + phi = convert (phi :: Posit256) :: Posit64 + gamma x = convert (gamma (convert x) :: Posit256) :: Posit64 + sinc x = convert (sinc (convert x) :: Posit256) :: Posit64 + expm1 x = + let b = atanh $ x / 2 + in (2 * b) / (1 - b) + +instance AltFloating Posit128 where + phi = convert (phi :: Posit256) :: Posit128 + gamma x = convert (gamma (convert x) :: Posit256) :: Posit128 + sinc x = convert (sinc (convert x) :: Posit256) :: Posit128 + expm1 x = + let b = atanh $ x / 2 + in (2 * b) / (1 - b) + +instance AltFloating Posit256 where + phi = funPhi 1.6 + gamma = funGammaSeries + sinc = funSinc + expm1 x = + let b = atanh $ x / 2 + in (2 * b) / (1 - b) + + +-- | 'phi' fixed point recursive algorithm, +funPhi :: Posit256 -> Posit256 +funPhi px@(Posit x) + | x == x' = Posit x + | otherwise = funPhi (Posit x') + where + (Posit x') = (px^2 + 2*px) / (px^2 + 1) + + +-- calculate atan(1/2^n) +-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached +funArcTan :: Natural -> Posit256 +funArcTan 0 = pi / 4 +funArcTan n + | n <= 122 = go 0 0 + | otherwise = z -- at small z... (atan z) == z "small angle approximation" + where + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit256 + term k = ((-1)^k * z^(2 * k + 1)) / fromIntegral (2 * k + 1) + z = 1 / 2^n -- recip $ 2^n :: Posit256 -- inv2PowN + +-- seems pretty close to 1 ULP with the input of 0.7813 +funAtan :: Posit256 -> Posit256 +funAtan NaR = NaR +funAtan x + | abs x < 1/2^122 = x -- small angle approximaiton, found emperically + | x < 0 = negate.funAtan $ negate x -- if negative turn it positive, it reduces the other domain reductions by half, found from Universal CORDIC + | x > 1 = pi/2 - funAtan (recip x) -- if larger than one use the complementary angle, found from Universal CORDIC + | x > twoMsqrt3 = pi/6 + funAtan ((sqrt 3 * x - 1)/(sqrt 3 + x)) -- another domain reduction, using an identity, found from https://mathonweb.com/help_ebook/html/algorithms.htm + | otherwise = funArcTanTaylor x +-- + +twoMsqrt3 :: Posit256 +twoMsqrt3 = 2 - sqrt 3 + +-- +funArcTanTaylor :: Posit256 -> Posit256 +funArcTanTaylor x = go 0 0 + where + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit256 + term k = ((-1)^k * x^(2 * k + 1)) / fromIntegral (2 * k + 1) +-- + +-- +funAsin :: Posit256 -> Posit256 +funAsin NaR = NaR +funAsin x + | abs x > 1 = NaR + | x == 1 = pi/2 + | x == -1 = -pi/2 + | otherwise = funAtan w + where + w = x / sqrt (1 - x^2) +-- + +-- +funAcos :: Posit256 -> Posit256 +funAcos NaR = NaR +funAcos x + | abs x > 1 = NaR + | x < 0 = pi + funAtan invw + | x == 0 = pi/2 + | x > 0 = funAtan invw + | otherwise = error "Prove it covers for Rational Numbers." + where + invw = sqrt (1 - x^2) / x +-- + +-- fI2PN = (1 /) . (2 ^) +funInv2PowN :: Natural -> Posit256 +funInv2PowN n = 1 / 2^n + + +-- calculate atanh(1/2^n) +-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached +funArcHypTan :: Natural -> Posit256 +funArcHypTan 0 = NaR +funArcHypTan n + | n <= 122 = go 0 0 + | otherwise = z -- at small z... (atan z) == z "small angle approximation" + where + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit256 + term k = (z^(2 * k + 1)) / fromIntegral (2 * k + 1) + z = 1 / 2^n + + +fac :: Natural -> Natural +fac 0 = 1 +fac n = n * fac (n - 1) + +-- +funAsinh :: Posit256 -> Posit256 +funAsinh NaR = NaR +funAsinh x = log $ x + sqrt (x^2 + 1) +-- + +-- +funAcosh :: Posit256 -> Posit256 +funAcosh NaR = NaR +funAcosh x + | x < 1 = NaR + | otherwise = log $ x + sqrt (x^2 - 1) +-- + +-- +funAtanh :: Posit256 -> Posit256 +funAtanh NaR = NaR +funAtanh x + | abs x >= 1 = NaR + | x < 0 = negate.funAtanh.negate $ x -- make use of odd parity to only calculate the positive part + | otherwise = 0.5 * log ((1+t) / (1-t)) - (fromIntegral ex / 2) * lnOf2 + where + (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat + 1, fromRational rat / 2) + (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x' -- sign should always be positive + x' = 1 - x + t = (2 - sig - x') / (2 + sig - x') +-- + +-- +funAtanhTaylor :: Posit256 -> Posit256 +funAtanhTaylor NaR = NaR +funAtanhTaylor x + | abs x >= 1 = NaR + | abs x < 1/2^122 = x -- small angle approximaiton, found emperically + | x < 0 = negate.funAtanhTaylor.negate $ x + | otherwise = go 0 0 + where + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit256 + term k = (x^(2 * k + 1)) / fromIntegral (2 * k + 1) +-- + +-- +funSin :: Posit256 -> Posit256 +funSin NaR = NaR +funSin 0 = 0 +funSin x = funSin' $ x / (2*pi) +-- +-- funSin' is sine normalized by 2*pi +funSin' :: Posit256 -> Posit256 +funSin' x + | x == 0 = 0 + | x == 0.25 = 1 + | x == 0.5 = 0 + | x == 0.75 = -1 + | x == 1 = 0 + | x < 0 = negate.funSin'.negate $ x + | x > 1 = + let (_,rem) = properFraction x + in funSin' rem + | x > 0.75 && x < 1 = negate.funSin' $ 1 - x -- reduce domain by quadrant symmetry + | x > 0.5 && x < 0.75 = negate.funSin' $ x - 0.5 + | x > 0.25 && x < 0.5 = funSin' $ 0.5 - x + | x > 0.125 && x < 0.25 = funCosTuma $ 2*pi * (0.25 - x) -- reduce domain and use cofunction + | otherwise = funSinTuma $ 2*pi * x +-- + +-- Taylor series expansion and fixed point algorithm, most accurate near zero +funSinTaylor :: Posit256 -> Posit256 +funSinTaylor NaR = NaR +funSinTaylor z = go 0 0 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit256 + term k = (-1)^k * z^(2*k+1) / (fromIntegral.fac $ 2*k+1) +-- + +-- +funSinTuma :: Posit256 -> Posit256 +funSinTuma NaR = NaR +funSinTuma z = go 19 1 + where + go :: Natural -> Posit256 -> Posit256 + go 1 !acc = z * acc + go !k !acc = go (k-1) (1 - (z^2 / fromIntegral ((2*k-2)*(2*k-1))) * acc) +-- + +-- +funCos :: Posit256 -> Posit256 +funCos NaR = NaR +funCos 0 = 1 +funCos x = funCos' $ x / (2*pi) +-- +-- funCos' is cosine normalized for 2*pi +funCos' :: Posit256 -> Posit256 +funCos' NaR = NaR +funCos' x + | x == 0 = 1 + | x == 0.25 = 0 + | x == 0.5 = -1 + | x == 0.75 = 0 + | x == 1 = 1 + | x < 0 = funCos'.negate $ x -- reduce domain by symmetry across 0 to turn x positive + | x > 1 = -- reduce domain by using perodicity + let (_,rem) = properFraction x + in funCos' rem + | x > 0.75 && x < 1 = funCos' $ 1 - x -- reduce domain by quadrant symmetry + | x > 0.5 && x < 0.75 = negate.funCos' $ x - 0.5 + | x > 0.25 && x < 0.5 = negate.funCos' $ 0.5 - x + | x > 0.125 && x < 0.25 = funSinTuma $ 2*pi * (0.25 - x) -- reduce domain and use cofunction + | otherwise = funCosTuma $ 2*pi * x -- +-- + +-- Taylor series expansion and fixed point algorithm, most accurate near zero +funCosTaylor :: Posit256 -> Posit256 +funCosTaylor NaR = NaR +funCosTaylor z = go 0 0 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit256 + term k = (-1)^k * z^(2*k) / (fromIntegral.fac $ 2*k) +-- + +-- +funCosTuma :: Posit256 -> Posit256 +funCosTuma NaR = NaR +funCosTuma z = go 19 1 + where + go :: Natural -> Posit256 -> Posit256 + go 1 !acc = acc + go !k !acc = go (k-1) (1 - (z^2 / fromIntegral ((2*k-3)*(2*k-2))) * acc) +-- + +-- ~16 ULP for 42 +funSinh :: Posit256 -> Posit256 +funSinh NaR = NaR +funSinh x = (exp x - exp (negate x))/2 +-- + +-- ~2 ULP for 42 +funSinhTaylor :: Posit256 -> Posit256 +funSinhTaylor NaR = NaR +funSinhTaylor z = go 0 0 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit256 + term k = z^(2*k+1) / (fromIntegral.fac $ 2*k+1) +-- + +-- +funSinhTuma :: Posit256 -> Posit256 +funSinhTuma NaR = NaR +funSinhTuma 0 = 0 +funSinhTuma z | z < 0 = negate.funSinhTuma.negate $ z +funSinhTuma z | z > 80 = 0.5 * funExpTuma z +funSinhTuma z = go 256 1 + where + go :: Natural -> Posit256 -> Posit256 + go 1 !acc = z * acc + go !k !acc = go (k-1) (1 + (z^2 / fromIntegral ((2*k-2) * (2*k-1))) * acc) +-- + +-- ~17 ULP for 42 +funCosh :: Posit256 -> Posit256 +funCosh NaR = NaR +funCosh x = (exp x + exp (negate x))/2 +-- + +-- ~3 ULP for 42 +funCoshTaylor :: Posit256 -> Posit256 +funCoshTaylor NaR = NaR +funCoshTaylor z = go 0 0 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit256 + term k = z^(2*k) / (fromIntegral.fac $ 2*k) +-- + +-- +funCoshTuma :: Posit256 -> Posit256 +funCoshTuma NaR = NaR +funCoshTuma 0 = 1 +funCoshTuma z | z < 0 = funCoshTuma.negate $ z +funCoshTuma z | z > 3 = 0.5 * (funExpTuma z + funExpTuma (negate z)) +funCoshTuma z = go 20 1 + where + go :: Natural -> Posit256 -> Posit256 + go 1 !acc = acc + go !k !acc = go (k-1) (1 + (z^2 / fromIntegral ((2*k-3)*(2*k-2)))*acc) +-- + + +-- +funLog :: Posit256 -> Posit256 +funLog x = funLog2 x * lnOf2 +-- + +-- +-- Use the constant, for performance +lnOf2 :: Posit256 +lnOf2 = Posit 28670435363615573179632300308403400109260626501925370561166468529302554498548 +-- + +-- +-- Some series don't converge reliably, this one does +funLnOf2 :: Posit256 +funLnOf2 = go 1 0 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit256 + term k = 1 / fromIntegral (2^k * k) +-- + +-- +funLog2 :: Posit256 -> Posit256 +funLog2 NaR = NaR +funLog2 z + | z <= 0 = NaR -- includes the NaR case + | otherwise = go (fromInteger ex) 1 sig -- domain reduction + where + go :: Posit256 -> Posit256 -> Posit256 -> Posit256 + go !acc !mak !sig' -- fixed point iteration, y is [1,2) :: Posit256 + | sig == 1 = acc + | acc == (acc + mak * 2^^(negate.fst.term $ sig')) = acc -- stop when fixed point is reached + | otherwise = go (acc + mak * 2^^(negate.fst.term $ sig')) (mak * 2^^(negate.fst.term $ sig')) (snd.term $ sig') + term = findSquaring 0 -- returns (m,s') m the number of times to square, and the new significand + (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat, fromRational rat) + (_,int,nat,rat) = (posit2TupPosit @V).toRational $ z -- sign should always be positive + findSquaring m s + | s >= 2 && s < 4 = (m, s/2) + | otherwise = findSquaring (m+1) (s^2) +-- + + +-- Gauss–Legendre algorithm, Seems only accurate to 2-3 ULP, but really slow +funPi1 :: Posit256 +funPi1 = go 0 3 1 (recip.sqrt $ 2) (recip 4) 1 + where + go :: Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 + go !prev !next !a !b !t !p + | prev == next = next + | otherwise = + let a' = (a + b) / 2 + b' = sqrt $ a * b + t' = t - p * (a - ((a + b) / 2))^2 + p' = 2 * p + in go next ((a' + b')^2 / (4 * t')) a' b' t' p' +-- + +-- 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 :: Posit256 +funPi2 = recip $ go 0 0 0.5 (5 / phi^3) + where + go :: Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256 + go !prev !n !a !s + | prev == a = a + | otherwise = + let x = 5 / s - 1 + y = (x - 1)^2 + 7 + z = (0.5 * x * (y + sqrt (y^2 - 4 * x^3)))**(1/5) + a' = s^2 * a - (5^n * ((s^2 - 5)/2 + sqrt (s * (s^2 - 2*s + 5)))) + s' = 25 / ((z + x/z + 1)^2 * s) + in go a (n+1) (trace (show a') a') s' +-- + + +-- Bailey–Borwein–Plouffe (BBP) formula, to 1-2 ULP, and blazing fast, converges in 60 iterations +funPi3 :: Posit256 +funPi3 = go 0 0 + where + go :: Integer -> Posit256 -> Posit256 + go !k !acc + | acc == acc + term k = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit256 + term k = fromRational $ (1 % 16^k) * ((120 * k^2 + 151 * k + 47) % (512 * k^4 + 1024 * k^3 + 712 * k^2 + 194 * k + 15)) +-- + + +-- Fabrice Bellard improvement on the BBP, 2-3 ULP, even faster, converges in 25 iterations, really fast +funPi4 :: Posit256 +funPi4 = (1/2^6) * go 0 0 + where + go :: Integer -> Posit256 -> Posit256 + go !k !acc + | acc == acc + term k = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit256 + term k = fromRational $ ((-1)^k % (2^(10*k))) * ((1 % (10 * k + 9)) - (2^2 % (10 * k + 7)) - (2^2 % (10 * k + 5)) - (2^6 % (10 * k + 3)) + (2^8 % (10 * k + 1)) - (1 % (4 * k + 3)) - (2^5 % (4 * k + 1))) +-- + + + +-- +-- looks to be about 4 ULP accurate at -100, right on the money at -1000 +funExp :: Posit256 -> Posit256 +funExp x = funExp2 funExpTaylor (x / lnOf2) +-- + +-- +-- +funExp2 :: (Posit256 -> Posit256) -> Posit256 -> Posit256 +funExp2 _ NaR = NaR +funExp2 _ 0 = 1 +funExp2 f x + | x < 0 = recip.funExp2 f.negate $ x -- always calculate the positive method + | otherwise = case properFraction x of + (int,rem) -> fromIntegral (2^int) * f (lnOf2 * rem) + + + +-- +-- calculate exp, its most accurate near zero +-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached +funExpTaylor :: Posit256 -> Posit256 +funExpTaylor NaR = NaR +funExpTaylor 0 = 1 +funExpTaylor z = go 0 0 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc -- if x == x + dx then terminate and return x + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit256 + term k = (z^k) / (fromIntegral.fac $ k) +-- + +-- +-- calculate exp, its most accurate near zero +-- use the Nested Series of Jan J Tuma +funExpTuma :: Posit256 -> Posit256 +funExpTuma NaR = NaR +funExpTuma 0 = 1 +funExpTuma z = go 57 1 -- was 66 + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | k == 0 = acc + | otherwise = go (k-1) (1 + (z / fromIntegral k) * acc) +-- + +-- +-- +funPow :: Posit256 -> Posit256 -> Posit256 +NaR `funPow` _ = NaR +_ `funPow` NaR = NaR +funPow 0 y + | y < 0 = NaR -- NaR: Divide by Zero + | y == 0 = NaR -- NaR: Indeterminate + | y > 0 = 0 +funPow x y + | y < 0 = recip $ funPow x (negate y) + | x < 0 = -- NaR if y is not an integer + let (int,rem) = properFraction y + in if rem == 0 + then x^^int + else NaR -- NaR: Imaginary Number + | otherwise = exp $ y * log x +-- + +-- Looks like 1 ULP for 0.7813 +funSinc :: Posit256 -> Posit256 +funSinc NaR = NaR +funSinc 0 = 1 -- Why the hell not! +funSinc theta = sin theta / theta +-- + +-- 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 +-- Shannon wavelet +funPsiSha1 :: Posit256 -> Posit256 +funPsiSha1 NaR = NaR +funPsiSha1 t = 2 * funSinc (2 * t) - funSinc t +-- + +-- Shannon wavelet +funPsiSha2 :: Posit256 -> Posit256 +funPsiSha2 NaR = NaR +funPsiSha2 t = funSinc (t/2) * cos (3*pi*t/2) +-- + +-- Shannon wavelet +funPsiSha3 :: Posit256 -> Posit256 +funPsiSha3 NaR = NaR +funPsiSha3 0 = 1 -- Why the hell not! +funPsiSha3 t = + let pit = pi * t + 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 +funLogDomainReduction _ 1 = 0 +funLogDomainReduction f x + | x <= 0 = NaR + | otherwise = f sig + (fromIntegral ex * lnOf2) + where + (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1 + (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x -- sign should always be positive + + + +-- natural log with log phi acurate to 9 ULP +funLogTaylor :: Posit256 -> Posit256 +funLogTaylor NaR = NaR +funLogTaylor 1 = 0 +funLogTaylor x | x <= 0 = NaR +funLogTaylor x + | x <= 2 = go 1 0 + | otherwise = error "The funLogTaylor algorithm is being used improperly" + where + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | acc == (acc + term k) = acc + | otherwise = go (k + 1) (acc + term k) + term :: Natural -> Posit256 + term k = (-1)^(k+1) * (x - 1)^k / fromIntegral k + + + + +-- natural log the Jan J Tuma way +funLogTuma :: Posit256 -> Posit256 +funLogTuma NaR = NaR +funLogTuma 1 = 0 -- domain reduced input is [0.5,1) and/or , where funLogTuma 1 = 0 +funLogTuma x | x <= 0 = NaR -- zero and less than zero is NaR +funLogTuma x + = go 242 1 + where + xM1 = x - 1 -- now [-0.5, 0) + go :: Natural -> Posit256 -> Posit256 + go !k !acc + | k == 0 = xM1 * acc + | otherwise = go (k-1) (recip (fromIntegral k) - xM1 * acc) + + +funGammaRamanujan :: Posit256 -> Posit256 +funGammaRamanujan z = sqrt pi * (x / exp 1)**x * (8*x^3 + 4*x^2 + x + (1/30))**(1/6) + where + x = z - 1 + +-- +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] + +funGammaSeries :: Posit256 -> Posit256 +funGammaSeries z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series) + where + series :: Posit256 + 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." +-- + +funGammaCalc :: Posit256 -> Posit256 +funGammaCalc z = sqrt (2*pi / z) * ((z / exp 1) * sqrt (z * sinh (recip z) + recip (810 * z^6)))**z + + +funGammaNemes :: Posit256 -> Posit256 +funGammaNemes z = sqrt (2*pi / z) * (recip (exp 1) * (z + recip (12 * z - recip (10 * z))))**z + +funGammaYang :: Posit256 -> Posit256 +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 :: Posit256 -> Posit256 +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 + + diff --git a/src/Posit/Internal/PositC.hs b/src/Posit/Internal/PositC.hs new file mode 100644 --- /dev/null +++ b/src/Posit/Internal/PositC.hs @@ -0,0 +1,983 @@ + +-------------------------------------------------------------------------------------------- +-- +-- Copyright : (C) 2022 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- Stability : Stable +-- Portability : Portable +-- +-- | Library implementing standard 'Posit-3.2' numbers, as defined by +-- the Posit Working Group 23 June 2018. +-- +-- +--------------------------------------------------------------------------------------------- + + +{-# LANGUAGE TypeFamilyDependencies #-} -- For the associated bidirectional type family that the Posit library is based on +{-# LANGUAGE DataKinds #-} -- For our ES kind and the constructors Z, I, II, III, IV, V, for exponent size type +{-# LANGUAGE TypeApplications #-} -- The most excellent syntax @Int256 +{-# 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 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 + +-- ---- +-- |Posit Class, implementing: +-- +-- * PositC +-- * Orphan Instances of Storable for Word128, Int128, Int256 +-- ---- + +module Posit.Internal.PositC +(PositC(..), + ES(..) + ) where + +import Prelude hiding (exponent,significand) + +-- Imports for Storable Instance of Data.DoubleWord +import Foreign.Storable (Storable, sizeOf, alignment, peek, poke) -- Used for Storable Instances of Data.DoubleWord +import Foreign.Ptr (Ptr, plusPtr, castPtr) -- Used for dealing with Pointers for the Data.DoubleWord Storable Instance + +-- Machine Integers and Operations +{-@ embed Int128 * as int @-} +{-@ embed Int256 * as int @-} +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 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) + + +-- | The Exponent Size 'ES' kind, the constructor for the Type is a Roman Numeral. +data ES = Z + | I + | II + | III + | IV + | V + + +-- | 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 + + + -- | 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. + decode :: IntN es -> Maybe Rational -- ^ You have an integer with a finite integer log2 word size decode it and Maybe it is Rational + + -- | 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 + + mostPosVal :: IntN es + leastPosVal :: IntN es + leastNegVal :: IntN es + mostNegVal :: IntN es + + -- Rational Value of common bounds + maxPosRat :: Rational + maxPosRat = 2^((nBytes @es) * ((nBits @es) - 2)) % 1 + minPosRat :: Rational + minPosRat = recip (maxPosRat @es) + maxNegRat :: Rational + 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) + log_uSeed (regime,r) + | 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 + in if signBit + then negate intRep + else intRep + + 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' + in regime' .|. exponent' .|. fraction + + formRegime power + | 0 <= power = + let offset = (fromIntegral (nBits @Z - 1) - power - 1) + in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) + | otherwise = + let offset = (fromIntegral (nBits @Z - 1) - abs power - 1) + in (1 `shiftL` fromInteger offset, offset) + + formExponent power offset = + let offset' = offset - fromIntegral (exponentSize @Z) + 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 @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) + + regime2Integer posit = + let regimeFormat = findRegimeFormat @Z posit + regimeCount = countRegimeBits @Z 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 + 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 @Z) - numBitsRegime - fromIntegral (signBitSize @Z) + signNRegimeMask = 2^bitsRemaining - 1 + int = posit .&. signNRegimeMask + nBitsToTheRight = fromIntegral (nBits @Z) - numBitsRegime - fromIntegral (signBitSize @Z) - fromIntegral (exponentSize @Z) + 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 @Z) + fromIntegral numBitsRegime + (exponentSize @Z) + fractionSize = fromIntegral (nBits @Z) - 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) + 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 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 === +-- ===================================================================== + + +-- getSign finds the sign value and then returns the absolute value of the Posit +getSign :: Rational -> (Bool, Rational) +getSign r = + let s = r <= 0 + absPosit = + if s + then negate r + 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) +getExponent :: Rational -> (Natural, Rational) +getExponent r = log_2 (0,r) + +log_2 :: (Natural, Rational) -> (Natural, Rational) +log_2 (exponent,r) | r < 1 = error "Should never happen, exponent should be a natural number, i.e. positive integer." + | r >= (2 % 1) = log_2 (exponent+1,r * (1 % 2)) + | otherwise = (exponent,r) + + +calcRegimeInt :: Bool -> Int -> Integer +calcRegimeInt format count | format = fromIntegral (count - 1) + | otherwise = fromIntegral $ negate count + + +xnor :: Bool -> Bool -> Bool +xnor a b = not $ (a || b) && not (b && a) + + +#ifndef O_NO_ORPHANS +#ifndef O_NO_STORABLE +-- ===================================================================== +-- === Storable Instances === +-- ===================================================================== +-- +-- Orphan Instance for Word128 using the DoubleWord type class +instance Storable Word128 where + sizeOf _ = 16 + alignment _ = 16 + peek ptr = do + hi <- peek $ offsetInt 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 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 +instance Storable Int128 where + sizeOf _ = 16 + alignment _ = 16 + peek ptr = do + hi <- peek $ offsetInt 0 + lo <- peek $ offsetWord 1 + return $ fromHiAndLo hi lo + where + offsetInt i = (castPtr ptr :: Ptr Int64) `plusPtr` (i*8) + offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8) + poke ptr int = do + poke (offsetInt 0) (hiWord int) + poke (offsetWord 1) (loWord int) + where + offsetInt i = (castPtr ptr :: Ptr Int64) `plusPtr` (i*8) + offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8) + +-- Orphan Instance for Int256 using the DoubleWord type class +instance Storable Int256 where + sizeOf _ = 32 + alignment _ = 32 + peek ptr = do + hi <- peek $ offsetInt 0 + lo <- peek $ offsetWord 1 + return $ fromHiAndLo hi lo + where + offsetInt i = (castPtr ptr :: Ptr Int128) `plusPtr` (i*16) + offsetWord i = (castPtr ptr :: Ptr Word128) `plusPtr` (i*16) + poke ptr int = do + poke (offsetInt 0) (hiWord int) + poke (offsetWord 1) (loWord int) + where + offsetInt i = (castPtr ptr :: Ptr Int128) `plusPtr` (i*16) + offsetWord i = (castPtr ptr :: Ptr Word128) `plusPtr` (i*16) +-- +#endif +#endif diff --git a/test/TestPosit.hs b/test/TestPosit.hs new file mode 100644 --- /dev/null +++ b/test/TestPosit.hs @@ -0,0 +1,183 @@ + +-------------------------------------------------------------------------------------------- +-- | Posit Numbers +-- Copyright : (C) 2022 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- Stability : Stable +-- Portability : Portable +-- +-- Test Suite for a Library implementing standard Posit Numbers +-- +--------------------------------------------------------------------------------------------- + +import Posit +import Posit.Internal.PositC + + + +main :: IO () +main = do +-- + 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)) -- + print $ "Machine Alpha Posit64 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit64)) -- + print $ "Machine Alpha Posit128 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit128)) -- + print $ "Machine Alpha Posit256 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit256)) -- + print $ "Does (1 - 1) == 0 ?: " ++ show ((1 - 1) == (0 :: Posit256)) -- [(1 - 1) == zero | zero = 0 :: Posit es, es <- Z .. V] + let sqrtTaylor = (funLogDomainReduction funLogTaylor).(/2).(funExp2 funExpTaylor).(/log 2) + print $ "sqrt phi using a Taylor algorithm: " ++ show (sqrtTaylor phi) + 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 = 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 + let truth = 2.6881171418161354484126255515800135873611118773741922415191608615280287034909564914158871097219845710811670879190576068697e43 :: Posit256 + eval "Standard: exp(100):" (exp (100)) truth + eval "Taylor: exp(100):" (funExp2 funExpTaylor (100 / log 2)) truth + eval "Tuma: exp(100):" (funExp2 funExpTuma (100 / log 2)) truth + eval "Wolfram Alpha: exp(100):" truth truth + let truth = 3.7200759760208359629596958038631183373588922923767819671206138766632904758958157181571187786422814966019356176423110698002e-44 :: Posit256 + eval "Standard: exp(-100):" (exp (-100)) truth + eval "Taylor: exp(-100):" (funExp2 funExpTaylor (-100 / log 2)) truth + eval "Tuma: exp(-100):" (funExp2 funExpTuma (-100 / log 2)) truth + eval "Wolfram Alpha: exp(-100):" truth truth + let truth = 1.9700711140170469938888793522433231253169379853238457899528029913850638507824411934749780765630268899309638179875202269359e434 :: Posit256 + eval "Standard: exp(1000):" (exp (1000)) truth + eval "Taylor: exp(1000):" (funExp2 funExpTaylor (1000 / log 2)) truth + eval "Tuma: exp(1000):" (funExp2 funExpTuma (1000 / log 2)) truth + eval "Wolfram Alpha: exp(1000):" truth truth + let truth = 5.075958897549456765291809479574336919305599282892837361832393845410540542974819175679662169046542867863667106831065285113e-435 :: Posit256 + eval "Standard: exp(-1000):" (exp (-1000)) truth + eval "Taylor: exp(-1000):" (funExp2 funExpTaylor (-1000 / log 2)) truth + 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 + let truth = -4.6051701859880913680359829093687284152022029772575459520666558019351452193547049604719944101791965966839355680845724972668190 :: Posit256 + eval "Standard: log(1/100):" (log (1/100)) truth + eval "Taylor: log(1/100):" (funLogDomainReduction funLogTaylor (1/100)) truth + eval "Tuma: log(1/100):" (funLogDomainReduction funLogTuma (1/100)) truth + eval "Wolfram Alpha: log(1/100):" truth truth + let truth = -6.9077552789821370520539743640530926228033044658863189280999837029027178290320574407079916152687948950259033521268587459002285 :: Posit256 + eval "Standard: log(1/1000):" (log (1/1000)) truth + eval "Taylor: log(1/1000):" (funLogDomainReduction funLogTaylor (1/1000)) truth + 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 + let truth = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: Posit256 + eval "Standard pi:" pi truth + eval "Gauss–Legendre algorithm: pi:" funPi1 truth + eval "Borwein's algorithm: pi:" funPi2 truth + eval "Bailey–Borwein–Plouffe (BBP) formula: pi:" funPi3 truth + eval "Fabrice Bellard improvement on the BBP: pi:" funPi4 truth + eval "Wolfram Alpha: pi:" truth truth + eval "Bailey–Borwein–Plouffe (BBP) formula: but succ pi:" (succ funPi3) truth +-- + -- print $ "Does (1 - 1) == 0 ?: " ++ (1 - 1) == (0 :: Posit256) -- [(1 - 1) == zero | zero = 0 :: Posit es, es <- Z .. V] + print "Now for Property testing of Posit8... (This should generalize for all other Posit types)" + print $ "Does associtivity of (+) hold?: " ++ (show assoc8) + print $ "Does commutitivity of (+) hold?: " ++ (show commutative8) + print $ "Is `fromInteger 0` the additive identity?: " ++ (show additiveIdent8) + print $ "Does `negate` give the additive inverse? (excluding NaR): " ++ (show additiveInv8) + print $ "Does `negate.negate == id`?: " ++ (show nn8) + print $ "Does associtivity of (*) hold?: " ++ (show assocMult8) + print $ "Is `fromInteger 1` the multiplicitive identity?: " ++ (show multIdent8) + print $ "Does Reflexivity of Eq hold?: " ++ (show reflEq8) + print $ "Does Symmetry of Eq hold?: " ++ (show symEq8) + print $ "Does Transitivity of Eq hold?: " ++ (show transEq8) + print $ "Does Extensionality of Eq hold?: " ++ (show extEq8) + print $ "Does Negation of Eq hold?: " ++ (show negEq8) + print $ "Does Comparability of Ord hold?: " ++ (show comp8) + print $ "Does Transitivity of Ord hold?: " ++ (show trans8) + print $ "Does Reflexivity of Ord hold?: " ++ (show refl8) + print $ "Does Antisymmetry of Ord hold?: " ++ (show anti8) + print $ "Does the `abs x * signum x == x` law hold?: " ++ (show absSignumLaw) + print $ "Is recip a multiplicative inverse?: " ++ (show recipInv8) + print $ "Are there any `recip.recip == id` values: " ++ (show rr8) + print $ "Are there any `recip.recip /= id` values: " ++ (show rrne8) + print $ "Does the distributive property hold with posits all the time?: " ++ (show doesItDistribute) + print $ "Exaustive Proof... for fused ops recovering the distributeive property... and it turns out to be true." + print $ "Can fused ops recover the distributive property for `fmms a b (negate a) c == fam b c a` ?: " ++ (show fusedDistribute) + + + +eval :: String -> Posit256 -> Posit256 -> IO () +eval msg val tru = putStr $ msg ++ "\n" ++ (show val) ++ "\n" ++ "ULP: " ++ (show $ valInt - truInt) ++ "\n" + where + valInt = read (displayIntegral val) :: Integer + truInt = read (displayIntegral tru) :: Integer + +-- exaustive testing, enum from to +assoc8 :: Bool +assoc8 = and [(x + y) + z == x + (y + z) | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), z <- enumFrom (NaR :: Posit8)] + +commutative8 :: Bool +commutative8 = and [x + y == y + x | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8)] + +additiveIdent8 :: Bool +additiveIdent8 = and [x + fromInteger 0 == x | x <- enumFrom (NaR :: Posit8)] + +additiveInv8 :: Bool +additiveInv8 = and [x + negate x == fromInteger 0 | x <- enumFrom (minBound :: Posit8)] + +assocMult8 :: Bool +assocMult8 = and [(x * y) * z == x * (y * z) | x <- enumFrom (minBound :: Posit8), y <- enumFrom (minBound :: Posit8), z <- enumFrom (minBound :: Posit8)] + +multIdent8 :: Bool +multIdent8 = and [x * fromInteger 1 == x && fromInteger 1 * x == x | x <- enumFrom (NaR :: Posit8)] + +reflEq8 = and [(x == x) | x <- enumFrom (NaR :: Posit8)] + +symEq8 = and [y == x | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), (x == y)] + +transEq8 = and [x == z | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), z <- enumFrom (NaR :: Posit8), (x == y) && (y == z)] + +extEq8 = and [sin x == sin y | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), x == y] + +negEq8 = and [not (x == y) | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), x /= y] + +comp8 :: Bool +comp8 = and [(x <= y || y <= x) == True | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8)] + +trans8 :: Bool +trans8 = and [(x <= z) == True | x <- enumFrom (minBound :: Posit8), y <- enumFrom (minBound :: Posit8), z <- enumFrom (minBound :: Posit8), (x <= y && y <= z) == True] + +refl8 :: Bool +refl8 = and [(x <= x) == True | x <- enumFrom (NaR :: Posit8)] + +anti8 :: Bool +anti8 = and [(x == y) == True | x <- enumFrom (minBound :: Posit8), y <- enumFrom (minBound :: Posit8), (x <= y && y <= x) == True] + +nn8 :: Bool +nn8 = and [(negate.negate $ x) == x | x <- enumFrom (NaR :: Posit8)] + +-- recip.recip == id +rr8 :: [Posit8] +rr8 = [x| x <- enumFrom (NaR :: Posit8), (recip.recip $ x) == x] + +-- recip.recip /= id +rrne8 :: [Posit8] +rrne8 = [x| x <- enumFrom (NaR :: Posit8), (recip.recip $ x) /= x] + +doesItDistribute :: Bool +doesItDistribute = and [a*b + a*c == a*(b+c) | a <- enumFrom (NaR :: Posit8), b <- enumFrom (NaR :: Posit8), c <- enumFrom (NaR :: Posit8)] + +fusedDistribute :: Bool +fusedDistribute = and [fmms a b (negate a) c == fam b c a | a <- enumFrom (NaR :: Posit8), b <- enumFrom (NaR :: Posit8), c <- enumFrom (NaR :: Posit8)] + +absSignumLaw :: Bool +absSignumLaw = and [abs x * signum x == x | x <- enumFrom (NaR :: Posit8)] + +recipInv8 :: Bool +recipInv8 = and [((x * recip x) == fromInteger 1) && ((recip x * x) == fromInteger 1) | x <- enumFrom (NaR :: Posit8)] +