diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,5 +1,15 @@ # Changelog for Posit Numbers +# posit-2022 + + * Added Types (P8, P16, P32, P64, P128, P256) for the Posit Standard 2022 encoding, exponent size = 2, and with nBytes = 2^es + * Refactored `Floating` to step up in resolution and then calculate a function, and then round it down to the the lower resolution + * Added polymorphic `Posit es` approximations for the `Floating` class + * Moved functions used in the test suite to the Test.Algorithms module, to eliminate the `do-test` flag + * Since the test flag has been removed the test can be run by: stack test + * Please forgive the lack of camelCase in some of the Floating functions... I think it reads better this time + * The Weigh test can be run as a benchmark: stack bench + # posit-3.2.0.5 * Bug fix for `mkIntRep` to resolve an overflow issue with the fractional part when it rounds up, in anticipation of the 2022 Standard release diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ -# posit 3.2.0.5 +# posit 2022.0.0.0 -The [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf), +The [Posit Standard 2022](https://posithub.org/docs/posit_standard-2.pdf), +and [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf), where Real numbers are approximated by Maybe Rational. The Posit Numbers are a drop in replacement for `Float` or `Double` mapped to a 2's complement integer type; smoothly and with tapering precision, in a @@ -25,11 +26,14 @@ * 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.4` 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. +word size are related. In `posit-3.2` es is instantiated as Z, I, +II, III, IV, V. In `posit-2022` es is instantiated as Z_2022, I_2022, +II_2022, III_2022, IV_2022, V_2022. 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' as well as, +'P8', 'P16', 'P32', 'P64', 'P128', and 'P256' are implemented and +include a couple of auxiliary classes, like AltShow, AltFloating, and +FusedOps. ``` class AltShow a where @@ -45,6 +49,7 @@ ``` class AltFloating p where + eps :: p -- Machine Epsilon near 1.0 phi :: p gamma :: p -> p sinc :: p -> p diff --git a/posit.cabal b/posit.cabal --- a/posit.cabal +++ b/posit.cabal @@ -1,8 +1,8 @@ cabal-version: 1.12 name: posit -version: 3.2.0.5 -description: The Posit Number format. Please see the README on GitHub at +version: 2022.0.0.0 +description: The Posit Number format attempting to conform to the Posit Standard Versions 3.2 and 2022. Where Real numbers are approximated by `Maybe Rational` and sampled in a similar way to the projective real line. homepage: https://github.com/waivio/posit#readme bug-reports: https://github.com/waivio/posit/issues author: Nathan Waivio @@ -14,11 +14,13 @@ tested-with: GHC == 8.10.4, GHC == 8.10.7, GHC == 9.0.2, - GHC == 9.2.5, + GHC == 9.2.7, GHC == 9.4.4 +synopsis: Posit Numbers extra-source-files: README.md ChangeLog.md + stack.yaml source-repository head type: git @@ -34,16 +36,11 @@ 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: @@ -63,8 +60,6 @@ if flag(do-liquid) cpp-options: -DO_LIQUID -DO_NO_STORABLE - if flag(do-test) - cpp-options: -DO_TEST -- Other library packages from which modules are imported. build-depends: @@ -77,19 +72,20 @@ if flag(do-liquid) build-depends: liquid-base, - liquidhaskell >= 0.8.10 + liquidhaskell -- perhaps one day: -threaded -rtsopts -with-rtsopts=-N test-suite posit-test type: exitcode-stdio-1.0 main-is: TestPosit.hs + other-modules: + Test.Algorithms hs-source-dirs: test ghc-options: -O2 - cpp-options: -DO_TEST build-depends: - base >=4.7 && <5 - , posit + base >=4.7 && <5, + posit default-language: Haskell2010 -- Weigh based benchmark for Vector @@ -99,8 +95,8 @@ main-is: WeighPosit.hs ghc-options: -Wall -O2 build-depends: - posit, base >=4.7 && <5, + posit, vector, weigh default-language: Haskell2010 diff --git a/src/Posit.hs b/src/Posit.hs --- a/src/Posit.hs +++ b/src/Posit.hs @@ -1,1450 +1,1003 @@ -------------------------------------------------------------------------------------------- -- 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, 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 #-} -{-# 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 -(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(..), - -#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, - -#ifdef O_TEST - -- * Alternative algorithms for test purposes - funExp, - funExp2, - funExpTaylor, - funLogTaylor, - funExpTuma, - funGammaSeriesFused, - funGammaRamanujan, - funGammaCalc, - funGammaNemes, - funGammaYang, - funGammaChen, - funGammaXminus1, - funLogTuma, - funLogDomainReduction, - funPi1, - funPi2, - funPi3, - funPi4, - funPi5, - funPi6, - funPsiSha1, - funPsiSha2, - funPsiSha3 -#endif - - ) 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 - --- 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 - - --- 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 ((%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions - --- 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 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 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 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 - (+) = 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 = R $ fromInteger int - -- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation - negate = viaIntegral negate --- - --- 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) - 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 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 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 - - - --- 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 - - --- ===================================================================== --- === Storable Instances === --- ===================================================================== --- -#ifndef O_NO_STORABLE --- -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 --- -#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 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) - -- LiquidHaskell is telling me this is unsafe if px is imaginary - -- lucky for us Posit256 is not imaginary - - --- 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 (2^(exponentSize @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 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202 --- - --- --- 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 (2^(exponentSize @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' --- - -#ifndef O_NO_SHOW --- Borwein's algorithm, with quintic convergence, --- gets to 7 ULP in 4 iterations, but really slow due to expensive function evaluations --- quite unstable and will not converge if sqrt is not accurate, which means log must be accurate -funPi2 :: Posit256 -funPi2 = recip $ go 0 0 0 0.5 (5 / phi^3) - where - go :: Posit256 -> Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256 - go !prevA !prevS !n !a !s - | prevA == a = a - | prevS == s = 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 s (n+1) (trace ("ΔA: " ++ show (a' - a)) a') (trace ("ΔS: " ++ show (s' - s)) s') --- -#endif - - --- 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))) --- - - --- Borwin's Quadradic Alogrithm 1985 -funPi5 :: Posit256 -funPi5 = recip $ go 0 0 1 (6 - 4 * sqrt 2) (sqrt 2 - 1) - where - go :: Posit256 -> Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256 - go !prevA !prevY !n a y - | prevA == a = a - | prevY == y = a - | otherwise = - let f = (1 - y^4)**(1/4) - y' = (1 - f) / (1 + f) - a' = a * (1 + y')^4 - 2^(2 * n + 1) * y' * (1 + y' + y'^2) - in if n == 3 - then a' - else go a y (n+1) (trace ("A: " ++ show a') a') (trace ("Y: " ++ show y') y') --- --- 3.14159265358979323846264338327950288419716939937510582097494459231 --- ULP: -97 - --- Borwin's Cubic Algirthm -funPi6 :: Posit256 -funPi6 = recip $ go 0 0 1 (1/3) ((sqrt 3 - 1) / 2) - where - go :: Posit256 -> Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256 - go !prevA !prevS !n !a !s - | prevA == a = a - | prevS == s = a - | otherwise = - let r = 3 / (1 + 2 * (1 - s^3)**(1/3)) - s'= (r - 1) / 2 - a'= r^2 * a - 3^(n-1) * (r^2 - 1) - in if n == 4 - then a' - else go a s (n+1) a' s' --- 3.14159265358979323846264338327950288419716939937510582097494459231 --- ULP: 216 - - --- --- 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. It a class of --- functions with the charicteristic of being a band pass filter in the --- frequency space. --- Shannon wavelet -funPsiSha1 :: Posit256 -> Posit256 -funPsiSha1 NaR = NaR -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, same as funPsiSha1 but with a factor of pi, with the --- Law: funPsiSha1.(pi*) === funPsiSha3 --- or : funPsiSha1 === funpsiSha3.(/pi) --- Posit256 seems to hold to a few ULP -funPsiSha3 :: Posit256 -> Posit256 -funPsiSha3 NaR = NaR -funPsiSha3 0 = 1 -- Why the hell not! -funPsiSha3 t = - let pit = pi * t - invpit = recip pit - in invpit * (sin (2 * pit) - sin pit) --- - - - --- 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 (2^(exponentSize @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." - -funGammaSeriesFused :: Posit256 -> Posit256 -funGammaSeriesFused z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series) - where - series :: Posit256 - series = fsumL $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n | n <- [1..len]] -- zipWith (\x y -> ) a001163 a001164 - lenA = length a001163 - lenB = length a001164 - len = if lenA == lenB - then lenA - else error "Seiries Numerator and Denominator do not have the same length." --- - -funGammaCalc :: Posit256 -> Posit256 -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 - -funGammaXminus1 :: Posit256 -> Posit256 -funGammaXminus1 x = go (x - 1) - where - go z = sqrt (2 * pi) * exp z ** (negate z) * z ** (z + 0.5) +-- Copyright : (C) 2022-2023 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- 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 + + +-- ---- +-- 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 + +-- 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 + + +-- 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 ((%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions + +-- 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 (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 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 + (+) = 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 = R $ fromInteger int + -- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation + negate = viaIntegral negate +-- + +-- 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) +-- + + + +-- Fractional Instances; (Num => Fractional) +-- +-- How the Frac do I get this definition: +instance 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 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 + + + +-- 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 + + +-- ===================================================================== +-- === Storable Instances === +-- ===================================================================== +-- +#ifndef O_NO_STORABLE +-- +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 +-- +#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 (PositC es, PositC (Next 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. (PositC es, PositC (Next 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. (PositC es, PositC (Next 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 + + +approx_exp :: PositC es => Posit es -> Posit es -- Comment by Abigale Emily: xcddfffff +approx_exp x = approx_2exp taylor_approx_exp (x / lnOf2) + + +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 x = normalizedSine $ x / (2*approx_pi) + + +approx_cos :: PositC es => Posit es -> Posit es +approx_cos NaR = NaR +approx_cos 0 = 1 +approx_cos x = normalizedCosine $ x / (2*approx_pi) + + +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 (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 (1 - x^2) / x + + +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 emperically + | 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 ((approx_sqrt 3 * x - 1)/(approx_sqrt 3 + 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 (x^2 + 1) + + +approx_acosh :: PositC es => Posit es -> Posit es +approx_acosh NaR = NaR +approx_acosh x + | x < 1 = NaR + | otherwise = approx_log $ x + approx_sqrt (x^2 - 1) + + +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 -- + + +-- Approximation of 2^x Domain Reduction +approx_2exp :: PositC es => (Posit es -> Posit es) -> Posit es -> Posit es +approx_2exp _ NaR = NaR +approx_2exp _ 0 = 1 +approx_2exp f x + | x < 0 = recip.approx_2exp f.negate $ x -- always calculate the positive method + | otherwise = case properFraction x of + (int,rem) -> fromIntegral (2^int) * f (lnOf2 * rem) + + + + +-- 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) +-- + + +-- calculate exp, its most accurate near zero +-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached +taylor_approx_exp :: forall es. PositC es => Posit es -> Posit es +taylor_approx_exp NaR = NaR +taylor_approx_exp 0 = 1 +taylor_approx_exp 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 + eps :: p + phi :: p + gamma :: p -> p + sinc :: p -> p + expm1 :: p -> p + +-- +instance PositC es => AltFloating (Posit es) where + phi = 1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374847540880753868917521266338 -- approx_phi 1.6 + eps = succ 1.0 - 1.0 + gamma = approx_gamma + sinc = approx_sinc + expm1 x = + let b = approx_atanh $ x / 2 + in (2 * b) / (1 - b) + + + + + + +approx_gamma :: forall es. PositC es => Posit es -> Posit es +approx_gamma z = approx_sqrt(2 * approx_pi) * (z `approx_pow` (z - 0.5)) * approx_exp (negate z) * (1 + series) + where + series :: Posit es + series = sum $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n | n <- [1..len]] -- zipWith (\x y -> ) a001163 a001164 + lenA = length a001163 + lenB = length a001164 + len = if lenA == lenB + then lenA + else error "Seiries Numerator and Denominator do not have the same length." +-- + + +-- Looks like 1 ULP for 0.7813 +approx_sinc :: PositC es => Posit es -> Posit es +approx_sinc NaR = NaR +approx_sinc 0 = 1 -- Why the hell not! +approx_sinc theta = approx_sin theta / theta +-- + + + +-- ===================================================================== +-- Useful Constants +-- ===================================================================== + +-- +-- Use the constant, for performance +lnOf2 :: PositC es => Posit es +lnOf2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202 +-- + +-- +a001163 :: [Integer] -- Numerator +a001163 = [1, 1, -139, -571, 163879, 5246819, -534703531, -4483131259, 432261921612371, 6232523202521089, -25834629665134204969, -1579029138854919086429, 746590869962651602203151, 1511513601028097903631961, -8849272268392873147705987190261, -142801712490607530608130701097701] +a001164 :: [Integer] -- Denominator +a001164 = [12, 288, 51840, 2488320, 209018880, 75246796800, 902961561600, 86684309913600, 514904800886784000, 86504006548979712000, 13494625021640835072000, 9716130015581401251840000, 116593560186976815022080000, 2798245444487443560529920000, 299692087104605205332754432000000, 57540880724084199423888850944000000] +-- + +twoMsqrt3 :: PositC es => Posit es +twoMsqrt3 = 2 - approx_sqrt 3 + + + +-- ===================================================================== +-- 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 + + + diff --git a/src/Posit/Internal/PositC.hs b/src/Posit/Internal/PositC.hs --- a/src/Posit/Internal/PositC.hs +++ b/src/Posit/Internal/PositC.hs @@ -7,8 +7,8 @@ -- Stability : Stable -- Portability : Portable -- --- | Library implementing standard 'Posit-3.2' numbers, as defined by --- the Posit Working Group 23 June 2018. +-- | Library implementing standard 'Posit-3.2', and 'Posit-2022' numbers, as defined by +-- the Posit Working Group 23 June 2018, and in 2022 respectively. -- -- --------------------------------------------------------------------------------------------- @@ -39,7 +39,9 @@ (PositC(..), ES(..), IntN, - FixedWidthInteger() + FixedWidthInteger(), + Max, + Next ) where import Prelude hiding (exponent,significand) @@ -65,29 +67,45 @@ -- | The Exponent Size 'ES' kind, the constructor for the Type is a Roman Numeral. -data ES = Z - | I - | II - | III - | IV - | V +data ES = Z_3_2 + | I_3_2 + | II_3_2 + | III_3_2 + | IV_3_2 + | V_3_2 + | Z_2022 + | I_2022 + | II_2022 + | III_2022 + | IV_2022 + | V_2022 -- | Type of the Finite Precision Representation, in our case Int8, -- Int16, Int32, Int64, Int128, Int256. {-@ embed IntN * as int @-} type family IntN (es :: ES) where - IntN Z = Int8 - IntN I = Int16 - IntN II = Int32 - IntN III = Int64 + IntN Z_3_2 = Int8 + IntN I_3_2 = Int16 + IntN II_3_2 = Int32 + IntN III_3_2 = Int64 #ifdef O_NO_STORABLE - IntN IV = Int128 - IntN V = Int256 + IntN IV_3_2 = Int128 + IntN V_3_2 = Int256 +#else + IntN IV_3_2 = Int128_Storable + IntN V_3_2 = Int256_Storable #endif -#ifndef O_NO_STORABLE - IntN IV = Int128_Storable - IntN V = Int256_Storable + IntN Z_2022 = Int8 + IntN I_2022 = Int16 + IntN II_2022 = Int32 + IntN III_2022 = Int64 +#ifdef O_NO_STORABLE + IntN IV_2022 = Int128 + IntN V_2022 = Int256 +#else + IntN IV_2022 = Int128_Storable + IntN V_2022 = Int256_Storable -- | New Type Wrappers to resolve Orphan Instance Issue newtype Int128_Storable = Int128_Storable Int128 @@ -101,6 +119,38 @@ via Word128 #endif + +-- | Type Max of Kind ES +type family Max (es :: ES) + where + Max Z_3_2 = V_3_2 + Max I_3_2 = V_3_2 + Max II_3_2 = V_3_2 + Max III_3_2 = V_3_2 + Max IV_3_2 = V_3_2 + Max V_3_2 = V_3_2 + Max Z_2022 = V_2022 + Max I_2022 = V_2022 + Max II_2022 = V_2022 + Max III_2022 = V_2022 + Max IV_2022 = V_2022 + Max V_2022 = V_2022 + +type family Next (es :: ES) + where + Next Z_3_2 = I_3_2 + Next I_3_2 = II_3_2 + Next II_3_2 = III_3_2 + Next III_3_2 = IV_3_2 + Next IV_3_2 = V_3_2 + Next V_3_2 = V_3_2 + Next Z_2022 = I_2022 + Next I_2022 = II_2022 + Next II_2022 = III_2022 + Next III_2022 = IV_2022 + Next IV_2022 = V_2022 + Next V_2022 = V_2022 + -- | The 'FixedWidthInteger' is a Constraint Synonym that contains all -- of the constraints provided by the 'IntN' Type Family. It is a super -- class for the Posit Class. @@ -152,8 +202,9 @@ in tupPosit2Posit @es (sgn,regime,exponent,rat) - -- | Exponent Size based on the Posit Exponent kind ES + -- | Exponent Size based on the Posit Exponent kind ES, Posit-2022 sets the default to 2. exponentSize :: Natural -- ^ The exponent size, 'es' is a Natural number + exponentSize = 2 -- | Various other size definitions used in the Posit format with their default definitions nBytes :: Natural -- ^ 'nBytes' the number of bytes of the Posit Representation @@ -229,7 +280,7 @@ fraction = formFraction @es significand offset' in regime' + exponent' + fraction -- Previously bad code... -- Was previously Bitwise OR'd (regime' .|. exponent' .|. fraction), but that failed when an overflow occurs in the fraction: - -- (R @es (6546781215792283740026379393655198304433284092086129578966582736192267592809066457889108741457440782093636999212155773298525238592782299216095867171579 % 6546781215792283740026379393655198304433284092086129578966582736192267592809349109766540184651808314301773368255120142018434513091770786106657055178752)) + -- (R @V_3_2 (6546781215792283740026379393655198304433284092086129578966582736192267592809066457889108741457440782093636999212155773298525238592782299216095867171579 % 6546781215792283740026379393655198304433284092086129578966582736192267592809349109766540184651808314301773368255120142018434513091770786106657055178752)) formRegime :: Integer -> (IntN es, Integer) formRegime power @@ -325,35 +376,49 @@ decimalPrec :: Int decimalPrec = fromIntegral $ 2 * (nBytes @es) + 1 - {-# MINIMAL exponentSize #-} + {-# MINIMAL exponentSize | nBytes #-} -- ===================================================================== -- === PositC Instances === -- ===================================================================== - -instance PositC Z where +-- | Standard 3.2 +instance PositC Z_3_2 where exponentSize = 0 - -instance PositC I where +instance PositC I_3_2 where exponentSize = 1 - -instance PositC II where +instance PositC II_3_2 where exponentSize = 2 - -instance PositC III where +instance PositC III_3_2 where exponentSize = 3 - -instance PositC IV where +instance PositC IV_3_2 where exponentSize = 4 - -instance PositC V where +instance PositC V_3_2 where exponentSize = 5 + +-- | Standard 2022 +instance PositC Z_2022 where + nBytes = 2^0 + +instance PositC I_2022 where + nBytes = 2^1 + +instance PositC II_2022 where + nBytes = 2^2 + +instance PositC III_2022 where + nBytes = 2^3 + +instance PositC IV_2022 where + nBytes = 2^4 + +instance PositC V_2022 where + nBytes = 2^5 diff --git a/stack.yaml b/stack.yaml new file mode 100644 --- /dev/null +++ b/stack.yaml @@ -0,0 +1,40 @@ +# This file is attempting to maintain the working Liquid Haskell versions +# that coorispond to a specific GHC or Stackage version + +# resolver: nightly-2023-03-30 # nightly-2023-02-20 # ghc-9.4.4 +resolver: lts-20.16 # ghc-9.2.7 # Currently the only version that seems to work with LiquidHaskell +# resolver: lts-19.33 # ghc-9.0.2 +# resolver: lts-18.28 # ghc-8.10.7 +# resolver: lts-18.6 # ghc-8.10.4 +# resolver: lts-16.31 # ghc-8.8.4 # Fails To Build! ghc: panic! (the 'impossible' happened) +# resolver: lts-14.27 # ghc-8.6.5 # Fails To Build! ghc: panic! (the 'impossible' happened) +packages: + - . +allow-newer: true +extra-deps: + # For LiquidHaskell: + - hashable-1.3.5.0 # lts-20.16 and below + # - hashable-1.4.2.0 # ghc-9.4.4 + - text-format-0.3.2 + - Diff-0.3.4 + - optparse-applicative-0.16.1.0 + # - rest-rewrite-0.3.0 # ye olde reliable + - rest-rewrite-0.4.1 # latest + - smtlib-backends-0.3 # ghc-9.2.7 + - smtlib-backends-process-0.3 # ghc-9.2.7 + - git: https://github.com/ucsd-progsys/liquidhaskell + # commit: # ghc-9.4.4 "Generically" errors out! Ambiguous occurrence ‘Generically’: It could refer to... ‘GHC.Generics.Generically’ or 'Language.Haskell.Liquid.Types.Generics.Generically' + commit: 63337d432b47c1ba1ec9925db0994fc5cdce3eaf # ghc-9.2.7 + # commit: b8780ee8d73d123adb39675ef87a2883f8aa1ecd # ghc-9.0.2 + # commit: f917323a1f9db1677e592d6ffc81467d53588d70 # ghc-8.10.7 + subdirs: + - . + - liquid-base + - liquid-vector + - liquid-bytestring + - liquid-containers + - liquid-ghc-prim + - git: https://github.com/ucsd-progsys/liquid-fixpoint + commit: 0e1a4725793740f495c26957044c56488d6e1efc # ghc-9.2.7 + # commit: 5aed39ec3210b9093ed635693d01bf351e25392f # ghc-9.0.2 + # commit: 544f8b0ba6d03b060701961250cce012412039c4 # ghc-8.10.7 diff --git a/test/Test/Algorithms.hs b/test/Test/Algorithms.hs new file mode 100644 --- /dev/null +++ b/test/Test/Algorithms.hs @@ -0,0 +1,605 @@ + + +{-# 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 ScopedTypeVariables #-} -- To reduce some code duplication, this is important +{-# LANGUAGE FlexibleContexts #-} -- to talk about class constraints like: (PositC es, PositC (Next es)) => +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE DataKinds #-} -- For our ES kind and the constructors Z, I, II, III, IV, V for exponent size type, post-pended with the version. +{-# OPTIONS_GHC -Wno-type-defaults #-} -- Turn off noise +{-# OPTIONS_GHC -Wno-unused-top-binds #-} -- Turn off noise + + +module Test.Algorithms + ( funLogDomainReduction + , funLogTaylor + , funExp2 + , funExpTaylor + , funLogTuma + , funExpTuma + , funGammaSeriesFused + , funGammaRamanujan + , funGammaCalc + , funGammaNemes + , funGammaYang + , funGammaChen + , funGammaXminus1 + -- , funGammaViaLngamma + , funPi1 + , funPi2 + , funPi3 + , funPi4 + , funPi5 + , funPi6 + ) where + +import Posit -- run with -O_TEST CPP directive + +import Prelude hiding (rem) + +-- 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 ((%)) -- 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 + + + +-- The machine implementation of the Posit encoding/decoding +import Posit.Internal.PositC -- The main internal implementation details + + -- Algorithms in Type: `Posit es` + + +-- ============================================================== +-- Other functions: +-- ============================================================== + + + +-- Approximation of log2 "Log Base 2" +approx_log2 :: forall es. PositC es => Posit es -> Posit es +approx_log2 NaR = NaR +approx_log2 z + | z <= 0 = NaR -- includes the NaR case + | otherwise = go (fromInteger ex) 1 sig -- domain reduction + where + go :: Posit es -> Posit es -> Posit es -> Posit es + 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 (2^(exponentSize @es)) + fromIntegral nat, fromRational rat) + (_,int,nat,rat) = (posit2TupPosit @es).toRational $ z -- sign should always be positive + findSquaring m s + | s >= 2 && s < 4 = (m, s/2) + | otherwise = findSquaring (m+1) (s^2) + + + +-- 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 + + + + + + + + + + +-- 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 + + + + + + + + + + +-- +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) +-- + + + + +-- 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) +-- + + + + + +-- 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) +-- + + +-- ~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) +-- + + +{- +-- | 'phi' fixed point recursive algorithm, +approx_phi :: (PositC es) => Posit es -> Posit es +approx_phi px@(Posit x) + | x == x' = Posit x + | otherwise = approx_phi (Posit x') + where + (Posit x') = (px^2 + 2*px) / (px^2 + 1) + -- LiquidHaskell is telling me this is unsafe if px is imaginary + -- lucky for us Posit256 is not imaginary +-} + + +-- +-- 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) +-- + + + +-- +-- Gauss–Legendre algorithm, Seems only accurate to 2-3 ULP, but really slow +funPi1 :: forall es. (PositC es, PositC (Next es)) => Posit es +funPi1 = go 0 3 1 (recip.sqrt $ 2) (recip 4) 1 + where + go :: Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es + 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 :: forall es. (PositC es, PositC (Next es)) => Posit es +funPi2 = recip $ go 0 0 0 0.5 (5 / phi^3) + where + go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es + go !prevA !prevS !n !a !s + | prevA == a = a + | prevS == s = a + | abs (prevA - a) < eps = a -- P256 will not reach a fixed point where `prevA == a` it sort of oscelates until divergence occurs, if we test for less than eps it can stop early + | 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 s (n+1) (trace ("ΔA: " ++ show (a' - a)) a') (trace ("ΔS: " ++ show (s' - s)) s') +-- + + + +-- Bailey–Borwein–Plouffe (BBP) formula, to 1-2 ULP, and blazing fast, converges in 60 iterations +funPi3 :: forall es. (PositC es) => Posit es +funPi3 = go 0 0 + where + go :: Integer -> Posit es -> Posit es + go !k !acc + | acc == acc + term k = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit es + 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 :: forall es. (PositC es) => Posit es +funPi4 = (1/2^6) * go 0 0 + where + go :: Integer -> Posit es -> Posit es + go !k !acc + | acc == acc + term k = acc + | otherwise = go (k+1) (acc + term k) + term :: Integer -> Posit es + 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))) +-- + + +-- Borwin's Quadradic Alogrithm 1985 +funPi5 :: forall es. (PositC es, PositC (Next es)) => Posit es +funPi5 = recip $ go 0 0 1 (6 - 4 * sqrt 2) (sqrt 2 - 1) + where + go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es + go !prevA !prevY !n a y + | prevA == a = a + | prevY == y = a + | otherwise = + let f = (1 - y^4)**(1/4) + y' = (1 - f) / (1 + f) + a' = a * (1 + y')^4 - 2^(2 * n + 1) * y' * (1 + y' + y'^2) + in if n == 3 + then a' + else go a y (n+1) (trace ("A: " ++ show a') a') (trace ("Y: " ++ show y') y') +-- +-- 3.14159265358979323846264338327950288419716939937510582097494459231 +-- ULP: -97 + +-- Borwin's Cubic Algirthm +funPi6 :: forall es. (PositC es, PositC (Next es)) => Posit es +funPi6 = recip $ go 0 0 1 (1/3) ((sqrt 3 - 1) / 2) + where + go :: Posit es -> Posit es -> Natural -> Posit es -> Posit es -> Posit es + go !prevA !prevS !n !a !s + | prevA == a = a + | prevS == s = a + | otherwise = + let r = 3 / (1 + 2 * (1 - s^3)**(1/3)) + s'= (r - 1) / 2 + a'= r^2 * a - 3^(n-1) * (r^2 - 1) + in if n == 4 + then a' + else go a s (n+1) a' s' +-- 3.14159265358979323846264338327950288419716939937510582097494459231 +-- ULP: 216 +-- +-- + + + + + + + +-- +-- 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) +-- + + + +-- +-- 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. It a class of +-- functions with the charicteristic of being a band pass filter in the +-- frequency space. +-- Shannon wavelet +funPsiSha1 :: Posit256 -> Posit256 +funPsiSha1 NaR = NaR +funPsiSha1 t = 2 * sinc (2 * t) - sinc t +-- + +-- Shannon wavelet +funPsiSha2 :: Posit256 -> Posit256 +funPsiSha2 NaR = NaR +funPsiSha2 t = sinc (t/2) * cos (3*pi*t/2) +-- + +-- Shannon wavelet, same as funPsiSha1 but with a factor of pi, with the +-- Law: funPsiSha1.(pi*) === funPsiSha3 +-- or : funPsiSha1 === funpsiSha3.(/pi) +-- Posit256 seems to hold to a few ULP +funPsiSha3 :: Posit256 -> Posit256 +funPsiSha3 NaR = NaR +funPsiSha3 0 = 1 -- Why the hell not! +funPsiSha3 t = + let pit = pi * t + invpit = recip pit + in invpit * (sin (2 * pit) - sin pit) +-- + + +-- +-- 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 (2^(exponentSize @V_3_2)) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1 + (_,int,nat,rat) = (posit2TupPosit @V_3_2).toRational $ x -- sign should always be positive +-- + +-- Use the constant, for performance +lnOf2 :: PositC es => Posit es +lnOf2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202 +-- + +-- +-- 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) +-- + +-- +-- +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) + + + +funGammaSeriesFused :: forall es. (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaSeriesFused z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series) + where + series :: Posit es + series = fsumL $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n | n <- [1..len]] -- zipWith (\x y -> ) a001163 a001164 + lenA = length a001163 + lenB = length a001164 + len = if lenA == lenB + then lenA + else error "Seiries Numerator and Denominator do not have the same length." +-- + +-- +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] + + +-- +-- 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 :: (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaRamanujan z = sqrt pi * (x / exp 1)**x * (8*x^3 + 4*x^2 + x + (1/30))**(1/6) + where + x = z - 1 +-- + + +-- +funGammaCalc :: (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaCalc z = sqrt (2*pi / z) * ((z / exp 1) * sqrt (z * sinh (recip z) + recip (810 * z^6)))**z + + +funGammaNemes :: (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaNemes z = sqrt (2*pi / z) * (recip (exp 1) * (z + recip (12 * z - recip (10 * z))))**z + +funGammaYang :: (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaYang z = sqrt (2 * pi * x) * (x / exp 1)**x * (x * sinh (recip x))**(x/2) * exp (fromRational (7 % 324) * recip (x^3 * (35 * x^2 + 33))) + where + x = z - 1 + +funGammaChen :: (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaChen z = sqrt (2 * pi * x) * (x / exp 1)**x * (1 + recip (12*x^3 + (24/7)*x - 0.5))**(x^2 + fromRational (53 % 210)) + where + x = z - 1 + +funGammaXminus1 :: (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaXminus1 x = go (x - 1) + where + go z = sqrt (2 * pi) * exp z ** (negate z) * z ** (z + 0.5) + +{- +funGammaInfProd :: Posit es -> Posit es +funGammaInfProd + + +funGammaViaInv :: Posit es -> Posit es +funGammaViaInv +-} +{- +funGammaViaLngamma :: forall es. (PositC es, PositC (Next es)) => Posit es -> Posit es +funGammaViaLngamma z = exp $ lngamma + where + lngamma :: Posit es + lngamma = negate eulersConstant * z - log z + go 0 1 + go :: Posit es -> Integer -> Posit es + go NaR _ = NaR + go prev k | prev == prev + next k = prev + | otherwise = go (trace ("Next: " ++ show (prev + next k)) (prev + next k)) (k+1) + next :: Integer -> Posit es + next k = z / fromIntegral k - (log $ 1 + z / fromIntegral k) + +eulersConstant :: PositC es => Posit es +eulersConstant = 0.57721566490153286060651209008240243104215933593992 +-} + +fac :: Natural -> Natural +fac 0 = 1 +fac n = n * fac (n - 1) +-- + + + + + + + + + + diff --git a/test/TestPosit.hs b/test/TestPosit.hs --- a/test/TestPosit.hs +++ b/test/TestPosit.hs @@ -16,6 +16,7 @@ import Posit import Posit.Internal.PositC +import Test.Algorithms 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 @@ -23,32 +24,58 @@ main :: IO () main = do -- - print $ "bitwise OR causes problem when fraction overflows Posit256: should be close to 1.0 not 0.5 ==> " ++ show (R @V (6546781215792283740026379393655198304433284092086129578966582736192267592809066457889108741457440782093636999212155773298525238592782299216095867171579 % 6546781215792283740026379393655198304433284092086129578966582736192267592809349109766540184651808314301773368255120142018434513091770786106657055178752)) + print $ "bitwise OR causes problem when fraction overflows Posit256: should be close to 1.0 not 0.5 ==> " ++ show (R @V_3_2 (6546781215792283740026379393655198304433284092086129578966582736192267592809066457889108741457440782093636999212155773298525238592782299216095867171579 % 6546781215792283740026379393655198304433284092086129578966582736192267592809349109766540184651808314301773368255120142018434513091770786106657055178752)) + print $ "bitwise OR causes problem when fraction overflows P256: should be close to 1.0 not 0.5 ==> " ++ show (R @V_2022 (6546781215792283740026379393655198304433284092086129578966582736192267592809066457889108741457440782093636999212155773298525238592782299216095867171579 % 6546781215792283740026379393655198304433284092086129578966582736192267592809349109766540184651808314301773368255120142018434513091770786106657055178752)) print $ "exp(1)**(pi*sqrt 43) :: Posit256 " ++ show (exp(1 :: Posit256) ** (pi * sqrt 43)) -- + print $ "exp(1)**(pi*sqrt 43) :: P256 " ++ show (exp(1 :: P256) ** (pi * sqrt 43)) -- print $ "exp(1)**(pi*sqrt 67) :: Posit256 " ++ show (exp(1 :: Posit256) ** (pi * sqrt 67)) -- - print $ "exp(1)**(pi*sqrt 163) :: Posit256 " ++ show (exp(1 :: Posit256) ** (pi * sqrt 163)) -- - print $ "Machine epsilon Posit8 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit8)) -- succ (Posit int) = Posit (succ int) - print $ "Machine epsilon Posit16 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit16)) -- - print $ "Machine epsilon Posit32 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit32)) -- - print $ "Machine epsilon Posit64 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit64)) -- - print $ "Machine epsilon Posit128 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit128)) -- - print $ "Machine epsilon Posit256 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit256)) -- + print $ "exp(1)**(pi*sqrt 67) :: P256 " ++ show (exp(1 :: P256) ** (pi * sqrt 67)) -- + print $ "exp(1)**(pi*sqrt 163):: Posit256 " ++ show (exp(1 :: Posit256) ** (pi * sqrt 163)) -- + print $ "exp(1)**(pi*sqrt 163):: P256 " ++ show (exp(1 :: P256) ** (pi * sqrt 163)) -- +-- | 'EPS' + print $ "Machine epsilon Posit8 ~1.0: " ++ show (eps :: Posit8) -- succ (Posit int) = Posit (succ int) + print $ "Machine epsilon Posit16 ~1.0: " ++ show (eps :: Posit16) -- + print $ "Machine epsilon Posit32 ~1.0: " ++ show (eps :: Posit32) -- + print $ "Machine epsilon Posit64 ~1.0: " ++ show (eps :: Posit64) -- + print $ "Machine epsilon Posit128 ~1.0: " ++ show (eps :: Posit128) -- + print $ "Machine epsilon Posit256 ~1.0: " ++ show (eps :: Posit256) -- + print $ "Machine epsilon P8 ~1.0: " ++ show (eps :: P8) -- succ (Posit int) = Posit (succ int) + print $ "Machine epsilon P16 ~1.0: " ++ show (eps :: P16) -- + print $ "Machine epsilon P32 ~1.0: " ++ show (eps :: P32) -- + print $ "Machine epsilon P64 ~1.0: " ++ show (eps :: P64) -- + print $ "Machine epsilon P128 ~1.0: " ++ show (eps :: P128) -- + print $ "Machine epsilon P256 ~1.0: " ++ show (eps :: P256) -- + -- | Taylor vs. Tuma 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 = 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: Posit256 - eval "Standard: gamma(phi): " (gamma (phi)) truth - eval "Fused Gamma: gamma(phi): " (funGammaSeriesFused (phi)) truth - eval "Ramanujan Gamma: gamma(phi): " (funGammaRamanujan (phi)) truth - eval "Calc Gamma: gamma(phi): " (funGammaCalc (phi)) truth - eval "Nemes Gamma: gamma(phi): " (funGammaNemes (phi)) truth - eval "Yang Gamma: gamma(phi): " (funGammaYang (phi)) truth - eval "Chen Gamma: gamma(phi): " (funGammaChen (phi)) truth - eval "Gamma (x - 1): gamma(phi): " (funGammaXminus1 (phi)) truth - eval "Wolfram alpha: gamma(phi): " truth truth + {- + let truthPosit256 = 0.8956731517052878608869612167009786079379812529831641161347143256 :: Posit256 -- 0.89566032673209158354178209470474131001971567786620187475744721557 :: Posit256 -- 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: Posit256 + let truthP256 = 0.8956731517052878608869612167009786079379812529831641161347143256 :: P256 -- 0.89566032673209158354178209470474131001971567786620187475744721557 :: P256 -- 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: P256 + eval "Standard: gamma(phi) :: Posit256 " (gamma (phi)) truthPosit256 + eval "Standard: gamma(phi) :: P256 " (gamma (phi)) truthP256 + eval "Fused Gamma: gamma(phi) :: Posit256 " (funGammaSeriesFused (phi)) truthPosit256 + eval "Fused Gamma: gamma(phi) :: P256 " (funGammaSeriesFused (phi)) truthP256 + eval "Ramanujan Gamma: gamma(phi) :: Posit256 " (funGammaRamanujan (phi)) truthPosit256 + eval "Ramanujan Gamma: gamma(phi) :: P256 " (funGammaRamanujan (phi)) truthP256 + eval "Calc Gamma: gamma(phi) :: Posit256 " (funGammaCalc (phi)) truthPosit256 + eval "Calc Gamma: gamma(phi) :: P256 " (funGammaCalc (phi)) truthP256 + eval "Nemes Gamma: gamma(phi) :: Posit256 " (funGammaNemes (phi)) truthPosit256 + eval "Nemes Gamma: gamma(phi) :: P256 " (funGammaNemes (phi)) truthP256 + eval "Yang Gamma: gamma(phi) :: Posit256 " (funGammaYang (phi)) truthPosit256 + eval "Yang Gamma: gamma(phi) :: P256 " (funGammaYang (phi)) truthP256 + eval "Chen Gamma: gamma(phi) :: Posit256 " (funGammaChen (phi)) truthPosit256 + eval "Chen Gamma: gamma(phi) :: P256 " (funGammaChen (phi)) truthP256 + eval "Gamma (x - 1): gamma(phi) :: Posit256 " (funGammaXminus1 (phi)) truthPosit256 + eval "Gamma (x - 1): gamma(phi) :: P256 " (funGammaXminus1 (phi)) truthP256 + eval "Calcuation of gamma(phi) using lngamma :: Posit256" (funGammaViaLngamma (phi)) truthPosit256 + eval "Calcuation of gamma(phi) using lngamma :: P256" (funGammaViaLngamma (phi)) truthP256 + eval "Wolfram alpha: gamma(phi) :: Posit256 " truthPosit256 truthPosit256 + eval "Wolfram alpha: gamma(phi) :: P256 " truthP256 truthP256 + -} let truth = 5.0431656433600286513118821892854247103235901754138463603020001967777869609108929428415187821843384653305404495551887666992776792 :: Posit256 eval "Standard: exp(phi):" (exp (phi)) truth eval "Taylor: exp(phi):" (funExp2 funExpTaylor (phi / log 2)) truth @@ -90,18 +117,26 @@ 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 Quintic algorithm: pi:" funPi2 truth - eval "Bailey–Borwein–Plouffe (BBP) formula: pi:" funPi3 truth - eval "Fabrice Bellard improvement on the BBP: pi:" funPi4 truth - eval "Borwein's Quadradic 1985 formula: pi:" funPi5 truth - eval "Borwein Cubic: pi:" funPi6 truth - eval "Wolfram Alpha: pi:" truth truth - eval "Bailey–Borwein–Plouffe (BBP) formula: but succ pi:" (succ funPi3) truth + eval "Standard: phi**pi:" ((phi) ** pi) truth + eval "Wolfram Alpha: phi**pi:" truth truth + let tPiPosit256 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: Posit256 + let tPiP256 = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: P256 + eval "Standard pi :: Posit256" pi tPiPosit256 + eval "Standard pi :: P256" pi tPiP256 + eval "Gauss–Legendre algorithm: pi :: Posit256" funPi1 tPiPosit256 + eval "Gauss–Legendre algorithm: pi :: P256" funPi1 tPiP256 + eval "Borwein's Quintic algorithm: pi :: Posit256" funPi2 tPiPosit256 + eval "Borwein's Quintic algorithm: pi :: P256" funPi2 tPiP256 + eval "Bailey–Borwein–Plouffe (BBP) formula: pi :: Posit256" funPi3 tPiPosit256 + eval "Bailey–Borwein–Plouffe (BBP) formula: pi :: P256" funPi3 tPiP256 + eval "Fabrice Bellard improvement on the BBP: pi :: Posit256" funPi4 tPiPosit256 + eval "Fabrice Bellard improvement on the BBP: pi :: P256" funPi4 tPiP256 + eval "Borwein's Quadradic 1985 formula: pi :: Posit256" funPi5 tPiPosit256 + eval "Borwein's Quadradic 1985 formula: pi :: P256" funPi5 tPiP256 + eval "Borwein Cubic: pi :: Posit256" funPi6 tPiPosit256 + eval "Borwein Cubic: pi :: P256" funPi6 tPiP256 + eval "Wolfram Alpha: pi :: Posit256" tPiPosit256 tPiPosit256 + eval "Wolfram Alpha: pi :: P256" tPiP256 tPiP256 -- -- 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)" @@ -131,7 +166,7 @@ -eval :: String -> Posit256 -> Posit256 -> IO () +eval :: (PositC es) => String -> Posit es -> Posit es -> IO () eval msg val tru = putStr $ msg ++ "\n" ++ (show val) ++ "\n" ++ "ULP: " ++ (show $ valInt - truInt) ++ "\n" where valInt = read (displayIntegral val) :: Integer @@ -200,4 +235,5 @@ recipInv8 :: Bool recipInv8 = and [((x * recip x) == fromInteger 1) && ((recip x * x) == fromInteger 1) | x <- enumFrom (NaR :: Posit8)] + diff --git a/test/WeighPosit.hs b/test/WeighPosit.hs --- a/test/WeighPosit.hs +++ b/test/WeighPosit.hs @@ -1,3 +1,7 @@ + +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE DataKinds #-} + import Weigh import Data.Vector.Storable as V @@ -6,34 +10,21 @@ main :: IO () main = mainWith $ do - func' "Posit8 in 1M Vector" vecOf unitPosit8 - func' "Posit16 in 1M Vector" vecOf unitPosit16 - func' "Posit32 in 1M Vector" vecOf unitPosit32 - func' "Posit64 in 1M Vector" vecOf unitPosit64 - func' "Posit128 in 1M Vector" vecOf unitPosit128 - func' "Posit256 in 1M Vector" vecOf unitPosit256 + func' "Posit8 in 1M Vector" vecOf (1.0 :: Posit8) + func' "Posit16 in 1M Vector" vecOf (1.0 :: Posit16) + func' "Posit32 in 1M Vector" vecOf (1.0 :: Posit32) + func' "Posit64 in 1M Vector" vecOf (1.0 :: Posit64) + func' "Posit128 in 1M Vector" vecOf (1.0 :: Posit128) + func' "Posit256 in 1M Vector" vecOf (1.0 :: Posit256) + func' "P8 in 1M Vector" vecOf (1.0 :: P8) + func' "P16 in 1M Vector" vecOf (1.0 :: P16) + func' "P32 in 1M Vector" vecOf (1.0 :: P32) + func' "P64 in 1M Vector" vecOf (1.0 :: P64) + func' "P128 in 1M Vector" vecOf (1.0 :: P128) + func' "P256 in 1M Vector" vecOf (1.0 :: P256) vecOf :: PositC es => Posit es -> V.Vector (Posit es) vecOf x = V.replicate (1024*1024) x - -unitPosit8 :: Posit8 -unitPosit8 = 1 - -unitPosit16 :: Posit16 -unitPosit16 = 1 - -unitPosit32 :: Posit32 -unitPosit32 = 1 - -unitPosit64 :: Posit64 -unitPosit64 = 1 - -unitPosit128 :: Posit128 -unitPosit128 = 1 - -unitPosit256 :: Posit256 -unitPosit256 = 1 -