posit-2022.2.0.3: src/Posit.hs
--------------------------------------------------------------------------------------------
-- Posit Numbers
-- Copyright : (C) 2022-2025 Nathan Waivio
-- License : BSD3
-- Maintainer : Nathan Waivio <nathan.waivio@gmail.com>
-- Stability : Stable
-- Portability : Portable
--
-- | Library implementing standard Posit Numbers both Posit Standard version
-- 3.2 and 2022, with some improvements. Posit is the interface, PositC
-- provides the implemetation. 2's Complement Fixed Point Integers,
-- and Rational numbers, are used throughout, as well as Integers & Naturals.
-- Encode and Decode are indexed through a Type Family.
--
---------------------------------------------------------------------------------------------
{-# 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, post-pended with the version.
{-# 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], and [P8 .. P256]
{-# LANGUAGE FlexibleContexts #-} -- If anybody knows what's this for let me know...
{-# 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, via Rational
{-# LANGUAGE ScopedTypeVariables #-} -- To reduce some code duplication, this is important
{-# 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
{-# LANGUAGE AllowAmbiguousTypes #-}
-- ----
-- Posit numbers implementing:
--
-- * Show
-- * Eq -- equality via an integer representation
-- * Ord -- compare via an integer representation
-- * Num -- Addition, subtraction, multiplication, and other operations most via Rational, negate is via an integer representation
-- * 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
(Posit(),
-- * Main Exported Types
Posit8, -- |A Posit-3.2 8-bit Posit number with 'exponentSize' = '0', and 1 byte wide
Posit16, -- |A Posit-3.2 16-bit Posit number with 'exponentSize' = '1', and 2 bytes wide
Posit32, -- |A Posit-3.2 32-bit Posit number with 'exponentSize' = '2', and 4 bytes wide
Posit64, -- |A Posit-3.2 64-bit Posit number with 'exponentSize' = '3', and 8 bytes wide
Posit128, -- |A Posit-3.2 128-bit Posit number with 'exponentSize' = '4', and 16 bytes wide
Posit256, -- |A Posit-3.2 256-bit Posit number with 'exponentSize' = '5', and 32 bytes wide
P8, -- |A Posit-2022 8-bit Posit number with 'exponentSize' = '2', and 1 byte wide
P16, -- |A Posit-2022 16-bit Posit number with 'exponentSize' = '2', and 2 bytes wide
P32, -- |A Posit-2022 32-bit Posit number with 'exponentSize' = '2', and 4 bytes wide
P64, -- |A Posit-2022 64-bit Posit number with 'exponentSize' = '2', and 8 bytes wide
P128, -- |A Posit-2022 128-bit Posit number with 'exponentSize' = '2', and 16 bytes wide
P256, -- |A Posit-2022 256-bit Posit number with 'exponentSize' = '2', and 32 bytes wide
-- * A Complete Pair of 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(..),
#ifndef O_NO_SHOW
-- * Additional functions to show the Posit in different formats
AltShow(..),
#endif
-- * Additional Special Functions
AltFloating(..),
-- * Functions to lift functions of Integers or Rationals to operate on Posit Types
viaIntegral,
viaRational,
viaRational2,
viaRational3,
viaRational4,
viaRational6,
viaRational8
) where
import Prelude hiding (rem)
-- Imports for Show and Read Instances
import Data.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
#ifndef O_NO_STORABLE_RANDOM
-- Imports for Storable Instance
import Foreign.Storable (Storable, sizeOf, alignment, peek, poke) -- Used for Storable Instances of Posit
import Foreign.Ptr (Ptr, castPtr) -- Used for dealing with Pointers for the Posit Storable Instance
import System.Random (Random(random,randomR))
import System.Random.Stateful (Uniform, uniform, uniformM)
import Data.Bits (shiftL, (.&.), (.|.))
#endif
-- 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 Numeric.Natural (Natural) -- Import the Natural Numbers ℕ (u+2115) for some of the Transcendental Functions
import Data.Ratio () -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions, no more (%) now.
-- for NFData instance
import Control.DeepSeq (NFData, rnf)
-- import Debug.Trace (trace) -- temporary for debug purposes
-- =====================================================================
-- === Posit Implementation ===
-- =====================================================================
-- The machine implementation of the Posit encoding/decoding
import Posit.Internal.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
-- |NFData Instance
instance NFData (Posit es) where
rnf (Posit _) = ()
-- |Not a Real Number, the Posit is like a Maybe type, it's either a real number or not
pattern NaR :: forall es. PositC es => Posit es
pattern NaR <- (Posit (decode @es -> Nothing)) where
NaR = Posit (unReal @es)
--
--
-- |A Real or at least Rational Number, rounded to the nearest Posit Rational representation
pattern R :: forall es. PositC es => Rational -> Posit es
pattern R r <- (Posit (decode @es -> Just r)) where
R r = Posit (encode @es $ Just r)
--
-- Posit functions are complete if the following two patterns are completely defined.
{-# COMPLETE NaR, R #-}
-- Concrete 3.2 types exported for use.
type Posit8 = Posit Z_3_2
type Posit16 = Posit I_3_2
type Posit32 = Posit II_3_2
type Posit64 = Posit III_3_2
type Posit128 = Posit IV_3_2
type Posit256 = Posit V_3_2
-- Concrete 2022 types exported for use.
type P8 = Posit Z_2022
type P16 = Posit I_2022
type P32 = Posit II_2022
type P64 = Posit III_2022
type P128 = Posit IV_2022
type P256 = Posit V_2022
#ifndef O_NO_SHOW
-- Show
--
instance PositC es => Show (Posit es) where
show NaR = "NaR"
show p@(Posit int) = formatScientific Generic (Just $ decimalPrec @es int) (fst.fromRationalRepetendUnlimited $ toRational p)
--
#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 PositC 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 PositC es => Ord (Posit es) where
compare (Posit int1) (Posit int2) = compare int1 int2
--
-- Num
--
-- I'm num trying to get this definition:
instance PositC es => Num (Posit es) where
-- Addition
(+) = positADD
-- Multiplication
(*) = positMULT
-- 'abs', Absolute Value, it's like a magnitude of sorts, abs of a posit is the same as abs of the integer representation
abs = positABS
-- 'signum' it is a kind of an representation of directionality, the sign of a number for instance
signum = positSIGNUM
-- 'fromInteger' rounds the integer into the closest posit number
fromInteger int = R $ fromInteger int
-- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation
negate = positNEG
-- '(-)', minus, but explisit with an implementaiton that will fuse
(-) = positSUB
--
--
-- To be able to have rewrite rules to function the instance needs to be some un-inlined function
{-# NOINLINE [1] positADD #-}
positADD :: forall es. PositC es => Posit es -> Posit es -> Posit es
positADD = viaRational2 (+)
{-# NOINLINE [1] positMULT #-}
positMULT :: forall es. PositC es => Posit es -> Posit es -> Posit es
positMULT = viaRational2 (*)
{-# NOINLINE [1] positNEG #-}
positNEG :: forall es. PositC es => Posit es -> Posit es
positNEG = viaIntegral negate
{-# NOINLINE [1] positABS #-}
positABS :: forall es. PositC es => Posit es -> Posit es
positABS = viaIntegral abs
{-# NOINLINE [1] positSIGNUM #-}
positSIGNUM :: forall es. PositC es => Posit es -> Posit es
positSIGNUM = viaRational signum
{-# NOINLINE [1] positSUB #-}
positSUB :: forall es. PositC es => Posit es -> Posit es -> Posit es
positSUB = \ p1 p2 -> positADD p1 (positNEG p2)
-- deriving via Integral Class, for the Integral representation of the posit
viaIntegral :: PositC es => (IntN es -> IntN es) -> Posit es -> Posit es
viaIntegral f (Posit int) = Posit $ f int
--
-- 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 PositC es => Enum (Posit es) where
-- succ (Posit int) = Posit (int + 1) -- Successor
succ = viaIntegral (+1) -- Posit Standard `next`
-- 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) -- Predicessor
pred = viaIntegral (subtract 1) -- Posit Standard `prior`
-- 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)
fromEnum _ = error "Please do not use 'fromEnum', it is size limited to Int, and can be machine dependant. Please advocate for the function to be size polymorphic of a FixedWidthInteger."
toEnum _ = error "Please do not use 'toEnum', it is size limited to Int, and can be machine dependant. Please advocate for the function to be size polymorphic of a FixedWidthInteger."
--
-- Fractional Instances; (Num => Fractional)
--
-- How the Frac do I get this definition:
instance PositC es => Fractional (Posit es) where
fromRational = R
_ / 0 = NaR
a / b = viaRational2 (/) a b
recip 0 = NaR
recip p = positRECIP p
--
{-# NOINLINE [1] positRECIP #-}
positRECIP :: forall es. PositC es => Posit es -> Posit es
positRECIP = viaRational recip
-- Rational Instances; Num & Ord Instanced => Real
--
-- I for real want this definition:
instance PositC es => Real (Posit es) where
toRational NaR = error "Your input is Not a Real or Rational (NaR) number, please try again!"
toRational (R r) = r
--
-- 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 => (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 => (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 => (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 => (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 => (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 => (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 PositC es => Bounded (Posit es) where
-- 'minBound' the most negative number represented
minBound = Posit (mostNegVal @es)
-- 'maxBound' the most positive number represented
maxBound = Posit (mostPosVal @es)
--
-- =====================================================================
-- === 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
#ifdef O_REWRITE
{-# RULES
"posit/fdot4" forall a0 a1 a2 a3 b0 b1 b2 b3. positADD (positADD (positADD (positMULT a0 b0) (positMULT a1 b1)) (positMULT a2 b2)) (positMULT a3 b3) = fdot4 a0 a1 a2 a3 b0 b1 b2 b3 -- (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)
"posit/fsum4" forall a b c d. positADD a (positADD b (positADD c d)) = fsum4 a b c d -- a + b + c + d
"posit/fdot3" forall a1 a2 a3 b1 b2 b3. positADD (positADD (positMULT a1 b1) (positMULT a2 b2)) (positMULT a3 b3) = fdot3 a1 a2 a3 b1 b2 b3 -- (a1 * b1) + (a2 * b2) + (a3 * b3)
"posit/fsum3" forall a b c. positADD a (positADD b c) = fsum3 a b c -- a + b + c
"posit/fsmSub" forall a b c. positSUB a (positMULT b c) = fsm a b c -- a - (b * c)
"posit/fsm" forall a b c. positADD a (positNEG (positMULT b c)) = fsm a b c -- a - (b * c)
"posit/fsmSwaped" forall a b c. positADD (positNEG (positMULT b c)) a = fsm a b c -- negate (b * c) + a
"posit/fma" forall a b c. positADD (positMULT a b) c = fma a b c -- (a * b) + c
"posit/fmaSwaped" forall a b c. positADD c (positMULT a b) = fma a b c -- c + (a * b)
"posit/fam" forall a b c. positMULT (positADD a b) c = fam a b c -- (a + b) * c
"posit/famSwaped" forall a b c. positMULT c (positADD a b) = fam a b c -- c * (a + b)
"posit/fmmsSub" forall a b c d. positSUB (positMULT a b) (positMULT c d) = fmms a b c d -- (a * b) - (c * d)
"posit/fmms" forall a b c d. positADD (positMULT a b) (positNEG (positMULT c d)) = fmms a b c d -- (a * b) - (c * d)
"posit/fmmsSwapped" forall a b c d. positADD (positNEG (positMULT c d)) (positMULT a b) = fmms a b c d -- negate (c * d) + (a * b)
#-}
#endif
-- Rational Instance
instance FusedOps Rational where
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 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 @es (Just $ go l 0)
where
go :: [Posit es] -> Rational -> Rational
go [] !acc = acc
go ((Posit int) : xs) !acc = case decode @es 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 @es (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 @es int1 of
Nothing -> error "First Posit List contains NaR"
Just r1 -> case decode @es 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 (PositC es1, PositC 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 PositC es => AltShow (Posit es) where
displayBinary (Posit int) = displayBin @es int
displayIntegral (Posit int) = show int
displayRational = viaShowable id
displayDecimal = viaShowable (fst.fromRationalRepetendUnlimited)
--
viaShowable :: (Show a, 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 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
#ifndef O_NO_STORABLE_RANDOM
-- =====================================================================
-- === Storable Instances ===
-- =====================================================================
--
--
instance 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
--
-- | Random instance for the Posit Sampling of R [0,1), this is for the
-- real numbers, not on the projective real numbers, for projective
-- real numbers, use Uniform.
instance forall es. PositC es => Random (Posit es) where
-- First we take a uniform distributed random posit, then we mask out
-- the sign, 2 bits of regime, and the exponent, then we write in the
-- sign, regime and exponent of 1.0, to get a posit [1,2) then subtract
-- 1.0 to adjust the range to [0,1). This approach is credited to a
-- coorispondance between Shin Yee Chung and John L. Gustafson titled:
-- "random number generators for posit" in the Unum Computing Google Group
random g = case uniform g of
(Posit int :: Posit es, g') -> (Posit ((int .&. maskFraction @es) .|. patt) - 1.0, g')
where
(Posit patt) = 1.0 :: Posit es
randomR (lo,hi) g
| lo > hi = randomR (hi,lo) g
| otherwise = case random g of
(p,g') -> let scaled_p = (hi - lo) * p + lo
in (scaled_p, g')
-- | Uniform instance for the Posit Sampling of the projective real line
instance PositC es => Uniform (Posit es) where
uniformM g = do
int <- uniformM g
return $ Posit int
maskFraction :: forall es. PositC es => IntN es
maskFraction =
let twoRegimeBits = 2 -- regimeBitSize set to 2, the good range is [1,2)
sreSize = signBitSize @es + twoRegimeBits + exponentSize @es
in (1 `shiftL` fromIntegral (nBits @es - sreSize) - 1)
#endif
-- =====================================================================
-- === Real Frac ===
-- =====================================================================
--
instance PositC es => RealFrac (Posit es) where
-- properFraction :: Integral b => a -> (b, a)
properFraction = viaRationalErrTrunkation "NaR value is not a RealFrac" properFraction
--
viaRationalErrTrunkation :: PositC es => String -> (Rational -> (a, Rational)) -> Posit es -> (a, Posit es)
viaRationalErrTrunkation err _ NaR = error err
viaRationalErrTrunkation _ f (R r) =
let (int, r') = f r
in (int, R r')
-- =====================================================================
-- === Real Float ===
-- =====================================================================
--
instance (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 PositF es => Floating (Posit es) where
pi = approx_pi
exp = hiRezNext approx_exp
log = hiRezNext approx_log
x ** y = hiRezNext2 approx_pow x y
sin = hiRezNext approx_sin
cos = hiRezNext approx_cos
asin = hiRezNext approx_asin
acos = hiRezNext approx_acos
atan = hiRezNext approx_atan
sinh = hiRezNext approx_sinh
cosh = hiRezNext approx_cosh
asinh = hiRezNext approx_asinh
acosh = hiRezNext approx_acosh
atanh = hiRezNext approx_atanh
-- Functions to step up and down in Resolution of the trancendental
-- functions so that we get properly rounded results upto 128-bits
-- Note: 256-bit resolution will not have ulp accuracy
hiRezNext :: forall es. PositF es => (Posit (Next es) -> Posit (Next es)) -> Posit es -> Posit es
hiRezNext f x = convert (f (convert x) :: Posit (Next es)) :: Posit es
hiRezMax :: forall es. (PositC es, PositC (Max es)) => (Posit (Max es) -> Posit (Max es)) -> Posit es -> Posit es
hiRezMax f x = convert (f (convert x) :: Posit (Max es)) :: Posit es
hiRezNext2 :: forall es. PositF es => (Posit (Next es) -> Posit (Next es) -> Posit (Next es)) -> Posit es -> Posit es -> Posit es
hiRezNext2 f x y = convert (f (convert x :: Posit (Next es)) (convert y :: Posit (Next es)) ) :: Posit es
hiRezMax2 :: forall es. (PositC es, PositC (Max es)) => (Posit (Max es) -> Posit (Max es) -> Posit (Max es)) -> Posit es -> Posit es -> Posit es
hiRezMax2 f x y = convert (f (convert x :: Posit (Max es)) (convert y :: Posit (Max es)) ) :: Posit es
-- =====================================================================
-- Approximations of Trancendental Funcitons
-- =====================================================================
approx_pi :: PositC es => Posit es
approx_pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446
-- euler's constant
approx_e :: PositC es => Posit es
approx_e = 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919320030599218174135
approx_exp :: forall es. PositC es => Posit es -> Posit es -- Comment by Abigale Emily: xcddfffff
approx_exp NaR = NaR
approx_exp (R x) = 2^^k * funExpTaylor (R r)
where
k = floor (x / ln2 + 0.5) -- should be Integer or Int
r = x - fromInteger k * ln2 -- should be Rational
ln2 :: Rational
ln2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202
log_USeed :: forall es. PositC es => Posit es
log_USeed = approx_log $ fromIntegral (uSeed @es)
approx_log :: PositC es => Posit es -> Posit es
approx_log = funLogDomainReduction funLogTaylor -- lnOf2 * approx_log2 x -- the commented out was slightly less accurate
approx_pow :: (PositC es) => Posit es -> Posit es -> Posit es
NaR `approx_pow` _ = NaR
_ `approx_pow` NaR = NaR
approx_pow 0 y
| y < 0 = NaR -- NaR: Divide by Zero
| y == 0 = NaR -- NaR: Indeterminate
| y > 0 = 0
approx_pow x y
| y < 0 = recip $ approx_pow 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 = approx_exp $ y * approx_log x
approx_sin :: forall es. PositC es => Posit es -> Posit es
approx_sin NaR = NaR
approx_sin 0 = 0
approx_sin (R x) = normalizedSine (R x')
where
(_, x') = properFraction $ x / twoPi
approx_cos :: forall es. PositC es => Posit es -> Posit es
approx_cos NaR = NaR
approx_cos 0 = 1
approx_cos (R x) = normalizedCosine (R x')
where
(_, x') = properFraction $ x / twoPi
approx_asin :: PositC es => Posit es -> Posit es
approx_asin NaR = NaR
approx_asin x
| abs x > 1 = NaR
| x == 1 = approx_pi/2
| x == -1 = -approx_pi/2
| otherwise = approx_atan w
where
w = x / approx_sqrt (fsm 1 x x) -- (1 - x^2)
approx_acos :: PositC es => Posit es -> Posit es
approx_acos NaR = NaR
approx_acos x
| abs x > 1 = NaR
| x < 0 = approx_pi + approx_atan invw
| x == 0 = approx_pi/2
| x > 0 = approx_atan invw
| otherwise = error "Prove it covers for Rational Numbers."
where
invw = approx_sqrt (fsm 1 x x) / x -- (1 - x^2)
approx_atan :: PositC es => Posit es -> Posit es
approx_atan NaR = NaR
approx_atan x
| abs x < 1/2^122 = x -- small angle approximaiton, found empirically
| x < 0 = negate.approx_atan $ negate x -- if negative turn it positive, it reduces the other domain reductions by half, found from Universal CORDIC
| x > 1 = approx_pi/2 - approx_atan (recip x) -- if larger than one use the complementary angle, found from Universal CORDIC
| x > twoMsqrt3 = approx_pi/6 + approx_atan ((sqrt3 * x - 1)/(sqrt3 + x)) -- another domain reduction, using an identity, found from https://mathonweb.com/help_ebook/html/algorithms.htm
| otherwise = taylor_approx_atan x
approx_sinh :: PositC es => Posit es -> Posit es
approx_sinh NaR = NaR
approx_sinh x = (approx_exp x - approx_exp (negate x))/2
approx_cosh :: PositC es => Posit es -> Posit es
approx_cosh NaR = NaR
approx_cosh x = (approx_exp x + approx_exp (negate x))/2
approx_asinh :: PositC es => Posit es -> Posit es
approx_asinh NaR = NaR
approx_asinh x = approx_log $ x + approx_sqrt (fma x x 1) -- (x^2 + 1) --
approx_acosh :: PositC es => Posit es -> Posit es
approx_acosh NaR = NaR
approx_acosh x
| x < 1 = NaR
| otherwise = approx_log $ x + approx_sqrt (fma x x (-1)) -- (x^2 - 1) --
approx_atanh :: forall es. PositC es => Posit es -> Posit es
approx_atanh NaR = NaR
approx_atanh x
| abs x >= 1 = NaR
| x < 0 = negate.approx_atanh.negate $ x -- make use of odd parity to only calculate the positive part
| otherwise = 0.5 * approx_log ((1+t) / (1-t)) - (fromIntegral ex / 2) * lnOf2
where
(ex, sig) = (int * fromIntegral (2^(exponentSize @es)) + fromIntegral nat + 1, fromRational rat / 2)
(_,int,nat,rat) = (posit2TupPosit @es).toRational $ x' -- sign should always be positive
x' = 1 - x
t = (2 - sig - x') / (2 + sig - x')
-- =====================================================================
-- Normalized Functions or Alternative Bases
-- =====================================================================
-- normalizedSine is sine normalized by 2*pi
normalizedSine :: PositC es => Posit es -> Posit es
normalizedSine NaR = NaR
normalizedSine x
| x == 0 = 0
| x == 0.25 = 1
| x == 0.5 = 0
| x == 0.75 = -1
| x == 1 = 0
| x < 0 = negate.normalizedSine.negate $ x
| x > 1 =
let (_,rem) = properFraction x
in normalizedSine rem
| x > 0.75 && x < 1 = negate.normalizedSine $ 1 - x -- reduce domain by quadrant symmetry
| x > 0.5 && x < 0.75 = negate.normalizedSine $ x - 0.5
| x > 0.25 && x < 0.5 = normalizedSine $ 0.5 - x
| x > 0.125 && x < 0.25 = tuma_approx_cos $ 2*approx_pi * (0.25 - x) -- reduce domain and use cofunction
| otherwise = tuma_approx_sin $ 2*approx_pi * x
-- normalizedCosine is cosine normalized for 2*pi
normalizedCosine :: PositC es => Posit es -> Posit es
normalizedCosine NaR = NaR
normalizedCosine x
| x == 0 = 1
| x == 0.25 = 0
| x == 0.5 = -1
| x == 0.75 = 0
| x == 1 = 1
| x < 0 = normalizedCosine.negate $ x -- reduce domain by symmetry across 0 to turn x positive
| x > 1 = -- reduce domain by using perodicity
let (_,rem) = properFraction x
in normalizedCosine rem
| x > 0.75 && x < 1 = normalizedCosine $ 1 - x -- reduce domain by quadrant symmetry
| x > 0.5 && x < 0.75 = negate.normalizedCosine $ x - 0.5
| x > 0.25 && x < 0.5 = negate.normalizedCosine $ 0.5 - x
| x > 0.125 && x < 0.25 = tuma_approx_sin $ 2*approx_pi * (0.25 - x) -- reduce domain and use cofunction
| otherwise = tuma_approx_cos $ 2*approx_pi * x --
-- Using the CORDIC domain reduction and some approximation function of log
funLogDomainReduction :: forall es. PositC es => (Posit es -> Posit es) -> Posit es -> Posit es
funLogDomainReduction _ NaR = NaR
funLogDomainReduction _ 1 = 0
funLogDomainReduction f x
| x <= 0 = NaR
| otherwise = f sig + (fromIntegral ex * lnOf2)
where
(ex, sig) = (int * fromIntegral (2^(exponentSize @es)) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1
(_,int,nat,rat) = (posit2TupPosit @es).toRational $ x -- sign should always be positive
-- natural log with log phi acurate to 9 ULP
funLogTaylor :: forall es. PositC es => Posit es -> Posit es
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 -> Posit es -> Posit es
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k + 1) (acc + term k)
term :: Natural -> Posit es
term k = (-1)^(k+1) * (x - 1)^k / fromIntegral k
-- =====================================================================
-- Taylor Series Fixed Point Approximations
-- =====================================================================
--
taylor_approx_atan :: forall es. PositC es => Posit es -> Posit es
taylor_approx_atan NaR = NaR
taylor_approx_atan x = go 0 0
where
go !k !acc
| acc == (acc + term k) = acc
| otherwise = go (k+1) (acc + term k)
term :: Integer -> Posit es
term k = ((-1)^k * x^(2 * k + 1)) / fromIntegral (2 * k + 1)
--
-- exp x - 1, has a decent Taylor Series
taylor_approx_expm1 :: forall es. PositC es => Posit es -> Posit es
taylor_approx_expm1 x = go 0 [x^n / fromIntegral (fac n) | n <- [1..]]
where
go :: Posit es -> [Posit es] -> Posit es
go !acc [] = acc
go !acc (h:t) | acc == acc + h = acc
| otherwise = go (acc + h) t
-- need to use a Taylor Series, the `tanh` formulation doesn't work because it requires something that depends on `exp`
--
-- calculate exp, its most accurate near zero
-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached
funExpTaylor :: forall es. PositC es => Posit es -> Posit es
funExpTaylor NaR = NaR
funExpTaylor 0 = 1
funExpTaylor z = go 0 0
where
go :: Natural -> Posit es -> Posit es
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 -> Posit es
term k = (z^k) / (fromIntegral.fac $ k)
-- =====================================================================
-- High Order Taylor Series transformed to Horner's Method
-- from Jan J Tuma's "Handbook of Numerical Calculations in Engineering"
-- =====================================================================
--
tuma_approx_cos :: forall es. PositC es => Posit es -> Posit es
tuma_approx_cos NaR = NaR
tuma_approx_cos z = go 19 1 -- TODO can the order be selected based on the word size?
where
go :: Natural -> Posit es -> Posit es
go 1 !acc = acc
go !k !acc = go (k-1) (1 - (z^2 / fromIntegral ((2*k-3)*(2*k-2))) * acc)
--
--
tuma_approx_sin :: forall es. PositC es => Posit es -> Posit es
tuma_approx_sin NaR = NaR
tuma_approx_sin z = go 19 1 -- TODO can the order be selected based on the word size?
where
go :: Natural -> Posit es -> Posit es
go 1 !acc = z * acc
go !k !acc = go (k-1) (1 - (z^2 / fromIntegral ((2*k-2)*(2*k-1))) * acc)
--
-- =========================================================
-- Alternate Floating of a Posit es
-- =========================================================
class AltFloating p where
machEps :: p
approxEq :: p -> p -> Bool
goldenRatio :: p
hypot2 :: p -> p -> p
hypot3 :: p -> p -> p -> p
hypot4 :: p -> p -> p -> p -> p
--
instance forall es. PositF es => AltFloating (Posit es) where
machEps = succ 1.0 - 1.0
approxEq a b =
let a' = convert a :: Posit (Prev es)
b' = convert b :: Posit (Prev es)
in a' == b'
goldenRatio = 1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374847540880753868917521266338 -- approx_phi 1.6
hypot2 a b = let a' :: Posit (Next es) = convert a
b' :: Posit (Next es) = convert b
in convert (approx_sqrt $ a'^2 + b'^2) :: Posit es
hypot3 a b c = let a' :: Posit (Next es) = convert a
b' :: Posit (Next es) = convert b
c' :: Posit (Next es) = convert c
in convert (approx_sqrt $ fsum3 (a'^2) (b'^2) (c'^2)) :: Posit es
hypot4 a b c d = let a' :: Posit (Next es) = convert a
b' :: Posit (Next es) = convert b
c' :: Posit (Next es) = convert c
d' :: Posit (Next es) = convert d
in convert (approx_sqrt $ fsum4 (a'^2) (b'^2) (c'^2) (d'^2)) :: Posit es
-- =====================================================================
-- Useful Constants used internally for Elementary Functions
-- =====================================================================
--
-- Use the constant, for performance
lnOf2 :: PositC es => Posit es
lnOf2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202
--
twoMsqrt3 :: PositC es => Posit es
twoMsqrt3 = 0.2679491924311227064725536584941276330571947461896193719441930205480669830911999629188538132427514243243738585845932969700300549
sqrt3 :: PositC es => Posit es
sqrt3 = 1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756756261414154067030299699450
twoPi :: Rational
twoPi = 6.2831853071795864769252867665590057683943387987502116419498891846156328125724179972560696506842341359642961730265646132941876892
-- =====================================================================
-- Helper Funcitons
-- =====================================================================
-- Factorial Function of type Natural
fac :: Natural -> Natural
fac 0 = 1
fac n = n * fac (n - 1)
--
approx_sqrt :: PositC es => Posit es -> Posit es
approx_sqrt x = approx_pow x 0.5