posit 3.2.0.5 → 2022.0.0.0

raw patch · 9 files changed

+1859/−1558 lines, 9 filesdep ~liquidhaskellPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: liquidhaskell

API changes (from Hackage documentation)

Files

• ChangeLog.md +10/−0
• README.md +12/−7
• posit.cabal +11/−15
• src/Posit.hs +1000/−1447
• src/Posit/Internal/PositC.hs +98/−33
• stack.yaml +40/−0
• test/Test/Algorithms.hs +605/−0
• test/TestPosit.hs +67/−31
• test/WeighPosit.hs +16/−25

@@ -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

@@ -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

@@ -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 <https://github.com/waivio/posit#readme>+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

@@ -1,1450 +1,1003 @@  -------------------------------------------------------------------------------------------- --   Posit Numbers---   Copyright   :  (C) 2022 Nathan Waivio---   License     :  BSD3---   Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>---   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 <nathan.waivio@gmail.com>+--   Stability   :  Stable+--   Portability :  Portable+--+-- | Library implementing standard Posit Numbers both Posit Standard version+--   3.2 and 2022, with some improvements.  Posit is the interface, PositC +--   provides the implemetation.  2's Complement Fixed Point Integers,+--   and Rational numbers, are used throughout, as well as Integers & Naturals.+--   Encode and Decode are indexed through a Type Family.+-- +---------------------------------------------------------------------------------------------+++{-# LANGUAGE GADTs #-} --   For our main type Posit (es :: ES)+{-# LANGUAGE DataKinds #-}  --   For our ES kind and the constructors Z, I, II, III, IV, V for exponent size type, post-pended with the version.+{-# LANGUAGE KindSignatures #-}  --   For defining the type of kind ES that indexes the GADT+{-# LANGUAGE ViewPatterns #-}  --   To decode the posit in the pattern+{-# LANGUAGE BangPatterns #-}  --   Added Strictness for some fixed point algorithms+{-# LANGUAGE PatternSynonyms #-}  --   for a nice NaR interface+{-# LANGUAGE FlexibleInstances #-} --   To make instances for each specific type [Posit8 .. Posit256], and [P8 .. P256]+{-# LANGUAGE FlexibleContexts #-} --   If anybody knows what's this for let me know...+{-# LANGUAGE TypeApplications #-} --   To apply types: @Type, it seems to select the specific class instance, when GHC is not able to reason about things, commenting this out shows an interesting interface+{-# LANGUAGE MultiParamTypeClasses #-}  --   To convert between Posit Types, via Rational+{-# LANGUAGE ScopedTypeVariables #-} --   To reduce some code duplication, this is important+{-# LANGUAGE UndecidableInstances #-}  --   To reduce some code duplication, I think the code is decidable but GHC is not smart enough ;), like there being only 1 instance that is polymorphic and works for all of my types.+{-# LANGUAGE CPP #-} --   To remove Storable instances to remove noise when performing analysis of Core+{-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-}  --   Turn off noise+{-# OPTIONS_GHC -Wno-type-defaults #-}  --   Turn off noise+{-# OPTIONS_GHC -Wno-unused-top-binds #-}  --   Turn off noise+++-- ----+--  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+++

@@ -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   

@@ -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: <something> # 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

@@ -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)+--++++++++++

@@ -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)]+ 

@@ -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-