posit (empty) → 3.2.0.0

raw patch · 8 files changed

+2762/−0 lines, 8 filesdep +basedep +data-dworddep +deepseqsetup-changed

Dependencies added: base, data-dword, deepseq, liquid-base, liquidhaskell, posit, scientific

Files

• ChangeLog.md +9/−0
• LICENSE +30/−0
• README.md +91/−0
• Setup.hs +2/−0
• posit.cabal +96/−0
• src/Posit.hs +1368/−0
• src/Posit/Internal/PositC.hs +983/−0
• test/TestPosit.hs +183/−0

@@ -0,0 +1,9 @@+# Changelog for Posit Numbers++## posit-3.2.0.0++  * Posit Standard 3.2 [Posit Standard] (https://posithub.org/docs/posit_standard.pdf)+  * LiquidHaskell support: stack build --flag posit:do-liquid+  * To run the test suite: stack test --flag posit:do-test+  * To play around: stack repl --flag posit:do-test+

@@ -0,0 +1,30 @@+Copyright Nathan Waivio (c) 2021-2022++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Nathan Waivio nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@@ -0,0 +1,91 @@+# posit 3.2.0.0++The [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf),+where Real numbers are approximated by Maybe Rational.  The Posit type+is mapped to a 2's complement integer type; smoothly and with tapering+precision, in a similar way to the projective real line.  The 'posit'+library implements the following standard classes:++ * Show+ * Eq+ * Ord  -- compare as an integer representation+ * Num  -- Addition, subtraction, multiplication, and other operations+ * Enum  -- Successor and Predecessor+ * Fractional  -- division, divide by zero is Not a Real (NaR) number+ * Real+ * Bounded+ * FusedOps  -- dot product and others+ * Convertable  -- Conversions between different posit formats+ * AltShow+ * Read+ * Storable  -- Formats for binary data, for computation and data interchange+ * RealFrac+ * RealFloat+ * Floating  -- Mathematical functions such as logarithm, exponential, trigonometric, and hyperbolic functions. Warning! May induce trance.++The Posits are indexed by the type (es :: ES) where exponent size and+word size are related.  In `posit-3.2.0.0` es is instantiated as Z, I,+II, III, IV, V.  The word size (in bits) of the value is `= 8 * 2^es`,+that is `2^es` bytes.  The Types: 'Posit8', 'Posit16', 'Posit32',+'Posit64', 'Posit128', and 'Posit256' are implemented and include a+couple of auxiliary classes, like AltShow, AltFloating, and FusedOps.++```+class AltShow a where+  -- Display the Posit in its Binary Representation+  displayBinary :: a -> String+  -- Display the Posit in its Integral Representation+  displayIntegral :: a -> String+  -- Display the Posit as a Rational+  displayRational :: a -> String+  -- Display the Posit as a Decimal until the Repented occurs+  displayDecimal :: a -> String+```++```+class AltFloating p where+  phi :: p+  gamma :: p -> p+  sinc :: p -> p+  expm1 :: p -> p+```++```+class Num a => FusedOps a where+  -- |Fused Multiply Add: (a * b) + c+  fma :: a -> a -> a -> a+  -- |Fused Add Multiply: (a + b) * c+  fam :: a -> a -> a -> a+  -- |Fused Multiply Multiply Subtract: (a * b) - (c * d)+  fmms :: a -> a -> a -> a -> a+  -- |Fused Sum of 3 values: a + b + c+  fsum3 :: a -> a -> a -> a+  -- |Fused Sum of 4 values: a + b + c + d+  fsum4 :: a -> a -> a -> a -> a+  -- |Fused Sum of a List of Posits+  fsumL :: Foldable t => t a -> a+  -- |Fused Dot Product of 3 element vector: (a1 * b1) + (a2 * b2) + (a3 * b3)+  fdot3 :: a -> a -> a -> a -> a -> a -> a+  -- |Fused Dot Product of 4 element veector: (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)+  fdot4 :: a -> a -> a -> a -> a -> a -> a -> a -> a+  -- |Fused Dot Product of Two Lists+  fdotL :: Foldable t => t a -> t a -> a+  -- |Fused Subtract Multiply: a - (b * c)+  fsm :: a -> a -> a -> a+```++The Posit type is 'Convertible' between other Posit lengths.++```+class Convertible a b where+  convert :: a -> b+```++The Posit Library is built on top of two of the most excellent libraries:+[data-dword](https://hackage.haskell.org/package/data-dword), and+[scientific](https://hackage.haskell.org/package/scientific).  The+'data-dword' library provides the underlining machine word+representation, it can provide 2^es word size, 2's complement fixed+length integers.  The 'scientific' library provides 'read' and 'show'+instances.+

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,96 @@+cabal-version: 1.12++name:           posit+version:        3.2.0.0+description:    The Posit Number format.  Please see the README on GitHub at <https://github.com/waivio/posit#readme>+homepage:       https://github.com/waivio/posit#readme+bug-reports:    https://github.com/waivio/posit/issues+author:         Nathan Waivio+maintainer:     nathan.waivio@gmail.com+copyright:      2021-2022 Nathan Waivio+license:        BSD3+license-file:   LICENSE+build-type:     Simple+extra-source-files:+    README.md+    ChangeLog.md++source-repository head+  type: git+  location: https://github.com/waivio/posit++flag do-no-storable+  description: Build without Storable Class support+  manual:      True+  default:     False++flag do-no-orphans+  description: Build without Orphan Instances if data-dword gets updated for Storable Instances+  manual:      True+  default:     False++flag do-liquid+  description: Build with Liquid Haskell checking+  manual:      True+  default:     False++flag do-test+  description: Export additional algorithms for calculating primitive functions for test purposes+  manual:      True+  default:     False++library+  exposed-modules:+      Posit+      Posit.Internal.PositC+  other-modules:+  hs-source-dirs:+      src+  build-depends:+      data-dword,+      scientific+  default-language: Haskell2010++  -- Compiler options+  ghc-options: -Wall -O2+ +  if flag(do-liquid)+    ghc-options: -fplugin=LiquidHaskell -fplugin-opt=LiquidHaskell:--fast -fplugin-opt=LiquidHaskell:--max-case-expand=4 -fplugin-opt=LiquidHaskell:--no-termination -fplugin-opt=LiquidHaskell:--short-names+ +  if flag(do-no-storable)+    cpp-options: -DO_NO_STORABLE+ +  if flag(do-no-orphans)+    cpp-options: -DO_NO_ORPHANS+ +  if flag(do-liquid)+    cpp-options: -DO_LIQUID -DO_NO_STORABLE -DO_NO_READ -DO_NO_SHOW+ +  if flag(do-test)+    cpp-options: -DO_TEST+ +  -- Other library packages from which modules are imported.+  build-depends:+    deepseq >=1.1 && <2+ +  if !flag(do-liquid)+    build-depends:+      base >=4.7 && <5+ +  if flag(do-liquid)+    build-depends:+      liquid-base,+      liquidhaskell >= 0.8.10++-- perhaps one day: -threaded -rtsopts -with-rtsopts=-N+test-suite posit-test+  type: exitcode-stdio-1.0+  main-is: TestPosit.hs+  hs-source-dirs:+      test+  ghc-options: -O2+  cpp-options: -DO_TEST+  build-depends:+      base >=4.7 && <5+    , posit+  default-language: Haskell2010

@@ -0,0 +1,1368 @@++--------------------------------------------------------------------------------------------+--   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.0.0, with some improvements) a fixed width word size of+--   2^es bytes.+-- +---------------------------------------------------------------------------------------------+++{-# LANGUAGE GADTs #-} --   For our main type Posit (es :: ES)+{-# LANGUAGE DataKinds #-}  --   For our ES kind and the constructors Z, I, II, III, IV, V for exponent size type+{-# LANGUAGE KindSignatures #-}  --   For defining the type of kind ES that indexes the GADT+{-# LANGUAGE ViewPatterns #-}  --   To decode the posit in the pattern+{-# LANGUAGE BangPatterns #-}  --   Added Strictness for some fixed point algorithms+{-# LANGUAGE PatternSynonyms #-}  --   for a nice NaR interface+{-# LANGUAGE FlexibleInstances #-} --   To make instances for each specific type [Posit8 .. Posit256]+{-# LANGUAGE FlexibleContexts #-}  --   Allow non-type variables in the constraints+{-# LANGUAGE TypeApplications #-} --   To apply types: @Type, it seems to select the specific class instance, when GHC is not able to reason about things, commenting this out shows an interesting interface+{-# LANGUAGE MultiParamTypeClasses #-}  --   To convert between Posit Types+{-# LANGUAGE ScopedTypeVariables #-} --   To reduce some code duplication+{-# LANGUAGE UndecidableInstances #-}  --   To reduce some code duplication, I think the code is decidable but GHC is not smart enough ;), like there being only 1 instance that is polymorphic and works for all of my types.+{-# LANGUAGE CPP #-} --   To remove Storable instances to remove noise when performing analysis of Core+{-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-}  --   Turn off noise+{-# OPTIONS_GHC -Wno-type-defaults #-}  --   Turn off noise+{-# OPTIONS_GHC -Wno-unused-top-binds #-}  --   Turn off noise++-- ----+--  Posit numbers implementing:+--+--    * Show+--    * Eq+--    * Ord  -- compare as an integer representation+--    * Num  -- Addition, subtraction, multiplication, and other operations+--    * Enum  -- Successor and Predecessor+--    * Fractional  -- division, divide by zero is Not a Real (NaR) number+--    * Real+--    * Bounded+--    * FusedOps  -- dot product and others+--    * Convertible  -- Conversions between different posit formats+--    * AltShow+--    * Read+--    * Storable  -- Formats for binary data, for computation and data interchange+--    * RealFrac+--    * RealFloat+--    * Floating  -- Mathematical functions such as logarithm, exponential, trigonometric, and hyperbolic functions. Warning! May induce trance.+--+-- ----++module Posit+(-- * Main Exported Types+ Posit8, -- |An 8-bit Posit number with 'es' ~ 'Z'+ Posit16, -- |An 16-bit Posit number with 'es' ~ 'I'+ Posit32, -- |An 32-bit Posit number with 'es' ~ 'II'+ Posit64, -- |An 64-bit Posit number with 'es' ~ 'III'+ Posit128, -- |An 128-bit Posit number with 'es' ~ 'IV'+ Posit256, -- |An 256-bit Posit number with 'es' ~ 'V'+ + -- * Patterns for Matching Exported Types+ pattern NaR,  -- |A pattern for Exception handling when a value is Not a Real number (NaR).+ pattern R,  -- |A pattern for the non-Exceptional case, yielding a Rational, will make a total function when paired with NaR, if the Rational implementation is total.+ + -- * Fused Operation Interface defined by the Posit Standard+ FusedOps(..),+ + -- * Posits are Convertable between different Posit representations+ Convertible(..),+ + -- * Additional functions to show the Posit in different formats+ AltShow(..),+ + -- * Additional Special Functions+ AltFloating(..),+ + -- * Functions to lift functions of Integers or Rationals to operate on Posit Types+ viaIntegral,+ viaRational,+ viaRational2,+ viaRational3,+ viaRational4,+ viaRational6,+ viaRational8,+#ifdef O_TEST+ -- * Alternative algorithms for test purposes+ funExp,+ funExp2,+ funExpTaylor,+ funLogTaylor,+ funExpTuma,+ funLogTuma,+ funLogDomainReduction,+ funPi1,+ funPi2,+ funPi3,+ funPi4,+ funPsiSha1,+ funPsiSha2,+ funPsiSha3+#endif+ ) where+++import Prelude hiding (rem)++-- Imports for Show and Read Instances+import Data.Scientific (Scientific+                       ,scientificP+                       ,fromRationalRepetendUnlimited+                       ,formatScientific+                       ,FPFormat(Generic)) -- Used to print/show and read the rational value++import Text.Read (Lexeme(Ident)+                 ,readPrec+                 ,readListPrec+                 ,(+++)+                 ,pfail+                 ,readListPrecDefault+                 ,lexP+                 ,lift+                 ,parens) -- Used to read a Posit value++-- Imports for Vectorization Class Instances+import Data.Foldable (toList)  -- Used for fused operations on foldable/lists++-- Imports for Storable Instance+import Foreign.Storable (Storable, sizeOf, alignment, peek, poke)  -- Used for Storable Instances of Posit+import Foreign.Ptr (Ptr, plusPtr, castPtr)  -- Used for dealing with Pointers for the Posit Storable Instance+++-- would like to:+-- import Posit.Internal.ElementaryFunctions+-- Perhaps on the chopping block if we are moving to ElementaryFunctions+-- Imports for implementing the Transcendental Functions+import GHC.Natural (Natural) -- Import the Natural Numbers ℕ (u+2115) for some of the Transcendental Functions+import Data.Ratio (Rational, (%))  -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions++import Debug.Trace (trace) -- temporary for debug purposes+++-- =====================================================================+-- ===                  Posit Implementation                         ===+-- =====================================================================++-- The machine implementation of the Posit encoding/decoding+import Posit.Internal.PositC (ES(..), PositC(..))  -- The main internal implementation details+++-- |Base GADT rapper type, that uses the Exponent Size kind to index the various implementations+data Posit (es :: ES) where+     Posit :: PositC es => !(IntN es) -> Posit es++-- |Not a Real Number, the Posit is like a Maybe type, it's either a real number or not+pattern NaR :: (PositC es) => Posit es+pattern NaR <- (Posit (decode -> Nothing)) where+  NaR = Posit unReal+--++--+-- |A Real or at least Rational Number, rounded to the nearest Posit Rational representation+pattern R :: PositC es => Rational -> Posit es+pattern R r <- (Posit (decode -> Just r)) where+  R r = Posit (encode $ Just r)+--++-- Posit functions are complete if the following two patterns are completely defined.+{-# COMPLETE NaR, R #-}++-- Concrete types exported for use.+type Posit8 = Posit Z+type Posit16 = Posit I+type Posit32 = Posit II+type Posit64 = Posit III+type Posit128 = Posit IV+type Posit256 = Posit V++#ifndef O_NO_SHOW+-- Show+--+instance forall es. (PositC es) => Show (Posit es) where+  show NaR = "NaR"+  show (R r) = formatScientific Generic (Just $ decimalPrec @es) (fst.fromRationalRepetendUnlimited $ r)+--+#endif++++-- Two Posit Numbers are Equal if their Finite Precision Integer representation is Equal+--+-- All things equal I would rather write it like this:+instance forall es. (Eq (IntN es)) => Eq (Posit es) where+  (Posit int1) == (Posit int2) = int1 == int2+--++++-- Two Posit Numbers are ordered by their Finite Precision Integer representation+--+-- Ordinarily I would only like one instance to cover them all+instance forall es. (Ord (IntN es), PositC es) => Ord (Posit es) where+  compare (Posit int1) (Posit int2) = compare int1 int2+--++++-- Num+--+-- I'm num trying to get this definition:+instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Num (Posit es) where+  -- Addition+  (+) = viaRational2 (+)+  -- Multiplication+  (*) = viaRational2 (*)+  -- 'abs', Absolute Value, it's like a magnitude of sorts, abs of a posit is the same as abs of the integer representation+  abs = viaIntegral abs+  -- 'signum' it is a kind of an representation of directionality, the sign of a number for instance+  signum = viaRational signum+  -- 'fromInteger' rounds the integer into the closest posit number+  fromInteger int = Posit $ encode (Just $ fromInteger int)+  -- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation+  negate = viaIntegral negate+--++-- deriving via Integral Class, for the Integral representation of the posit+viaIntegral :: forall es. PositC es => (IntN es -> IntN es) -> Posit es -> Posit es+viaIntegral f (Posit int) = Posit $ f int+--++++-- Enum-ish, A Posit has a Successor and Predecessor so its an ordinal number, as per Posit standard next, prior+-- The Posit Standard requires 2's complement integer overflow to be ignored+instance forall es. (Num (IntN es), Enum (IntN es), Ord (IntN es), PositC es) => Enum (Posit es) where+  -- succ (Posit int) = Posit (int + 1)+  succ = viaIntegral (+1)+  -- succ = viaIntegral succ  -- Non-compliant, runtime error pred NaR, and worse it is Int64 for types of greater precision, probably because of Preludes gross abomination of toEnum/fromEnum+  -- pred (Posit int) = Posit (int - 1)+  pred = viaIntegral (subtract 1)+  -- pred = viaIntegral pred  -- Non-compliant, runtime error pred NaR, and worse it is Int64 for types of greater precision, probably because of Preludes gross abomination of toEnum/fromEnum+  -- enumFrom :: Posit es -> [Posit es]+  enumFrom n = enumFromTo n maxBound+  enumFromTo n m+    | n == m = [n]+    | n < m = n : enumFromTo (succ n) m+    | otherwise = []+  -- enumFromThen n m :: Posit es -> Posit es -> [Posit es]+  enumFromThen NaR _ = [NaR]+  enumFromThen _ NaR = [NaR]+  enumFromThen n m = n : go n+    where+      step = m - n+      go :: Posit es -> [Posit es]+      go NaR = [NaR]+      go !l = case compare step 0 of+                LT -> let !n' = l + step  -- rounding occurs here, because the next comparison needs it, it wouldn't make sense otherwise...+                      in if n' - l > step+                         then []+                         else n' : go n'+                EQ -> [n, m]+                GT -> let !n' = l + step+                      in if n' - l < step+                         then []  -- with tapered resolution this algorithm can reach a fixed point where the next value is equal to the previous value+                         else n' : go n'+  enumFromThenTo NaR  _   _  = [NaR]+  enumFromThenTo  _  NaR  _  = [NaR]+  enumFromThenTo  _   _  NaR = [NaR]+  enumFromThenTo  e1  e2  e3 = takeWhile predicate (enumFromThen e1 e2)+    where+      mid = (e2 - e1) / 2+      predicate | e2 >= e1  = (<= e3 + mid)+                | otherwise = (>= e3 + mid)+--++++-- Fractional Instances; (Num => Fractional)+--+-- How the Frac do I get this definition:+instance forall es. (Num (IntN es), Ord (IntN es), Eq (IntN es), PositC es) => Fractional (Posit es) where+  fromRational = R+ +  recip 0 = NaR+  recip p = viaRational recip p+--++-- Rational Instances; Num & Ord Instanced => Real+--+-- I for real want this definition:+instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => Real (Posit es) where+  toRational NaR = error "Your input is Not a Real or Rational (NaR) number, please try again!"+  toRational (R r) = r+--++-- Implementing instances via Rational Data Type's instance,+-- The function checks for NaR, to protect against the runtime error 'toRational' would generate if called with a NaR value+-- Unary::Arity NaR guarded pass through with wrapping and unwrapping use of a Rational function+viaRational :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational) -> Posit es -> Posit es+viaRational _ NaR = NaR+viaRational f (R r) = fromRational $ f r++-- Binary NaR guarded pass through with wrapping and unwrapping use of a Rational function+viaRational2 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es+viaRational2 _ NaR  _  = NaR+viaRational2 _  _  NaR = NaR+viaRational2 f (R r1) (R r2) = R $ r1 `f` r2++-- Ternary NaR guarded pass through with wrapping and unwrapping use of a Rational function+viaRational3 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es+viaRational3 _ NaR  _   _  = NaR+viaRational3 _  _  NaR  _  = NaR+viaRational3 _  _   _  NaR = NaR+viaRational3 f (R r1) (R r2) (R r3) = R $ f r1 r2 r3++-- Quaternary NaR guarded pass through with wrapping and unwrapping use of a Rational function+viaRational4 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es+viaRational4 _ NaR  _   _   _  = NaR+viaRational4 _  _  NaR  _   _  = NaR+viaRational4 _  _   _  NaR  _  = NaR+viaRational4 _  _   _   _  NaR = NaR+viaRational4 f (R r0) (R r1) (R r2) (R r3) = R $ f r0 r1 r2 r3++-- Senary NaR guarded pass through with wrapping and unwrapping use of a Rational function+viaRational6 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es+viaRational6 _ NaR  _   _   _   _   _  = NaR+viaRational6 _  _  NaR  _   _   _   _  = NaR+viaRational6 _  _   _  NaR  _   _   _  = NaR+viaRational6 _  _   _   _  NaR  _   _  = NaR+viaRational6 _  _   _   _   _  NaR  _  = NaR+viaRational6 _  _   _   _   _   _  NaR = NaR+viaRational6 f (R a1) (R a2) (R a3) (R b1) (R b2) (R b3) = R $ f a1 a2 a3 b1 b2 b3++-- Octonary NaR guarded pass through with wrapping and unwrapping use of a Rational function+viaRational8 :: (PositC es, Ord (IntN es), Num (IntN es)) => (Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational -> Rational) -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es -> Posit es+viaRational8 _ NaR  _   _   _   _   _   _   _  = NaR+viaRational8 _  _  NaR  _   _   _   _   _   _  = NaR+viaRational8 _  _   _  NaR  _   _   _   _   _  = NaR+viaRational8 _  _   _   _  NaR  _   _   _   _  = NaR+viaRational8 _  _   _   _   _  NaR  _   _   _  = NaR+viaRational8 _  _   _   _   _   _  NaR  _   _  = NaR+viaRational8 _  _   _   _   _   _   _  NaR  _  = NaR+viaRational8 _  _   _   _   _   _   _   _  NaR = NaR+viaRational8 f (R a0) (R a1) (R a2) (R a3) (R b0) (R b1) (R b2) (R b3) = R $ f a0 a1 a2 a3 b0 b1 b2 b3++++-- Bounded, bounded to what?!? To the ℝ! NaR is out of bounds!!!+--+-- I'm bound to want this definition:+instance forall es. PositC es => Bounded (Posit es) where+  -- 'minBound' the most negative number represented+  minBound = Posit mostNegVal+  -- 'maxBound' the most positive number represented+  maxBound = Posit mostPosVal+--+++-- =====================================================================+-- ===                    Fused Operations                           ===+-- =====================================================================++-- |A class that delays the rounding operation until the end for some operations+class Num a => FusedOps a where+  -- |Fused Multiply Add: (a * b) + c+  fma :: a -> a -> a -> a+  -- |Fused Add Multiply: (a + b) * c+  fam :: a -> a -> a -> a+  -- |Fused Multiply Multiply Subtract: (a * b) - (c * d)+  fmms :: a -> a -> a -> a -> a+  -- |Fused Sum of 3 values: a + b + c+  fsum3 :: a -> a -> a -> a+  -- |Fused Sum of 4 values: a + b + c + d+  fsum4 :: a -> a -> a -> a -> a+  -- |Fused Sum of a List of Posits+  fsumL :: Foldable t => t a -> a+  -- |Fused Dot Product of 3 element vector: (a1 * b1) + (a2 * b2) + (a3 * b3)+  fdot3 :: a -> a -> a -> a -> a -> a -> a+  -- |Fused Dot Product of 4 element vector: (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)+  fdot4 :: a -> a -> a -> a -> a -> a -> a -> a -> a+  -- |Fused Dot Product of Two Lists+  fdotL :: Foldable t => t a -> t a -> a+  -- |Fused Subtract Multiply: a - (b * c)+  fsm :: a -> a -> a -> a+ +++-- Rational Instance+instance FusedOps Rational where+  fsm a b c = a - (b * c)+  fma a b c = (a * b) + c+  fam a b c = (a + b) * c+  fmms a b c d = (a * b) - (c * d)+  fsum3 a b c = a + b + c+  fsum4 a b c d = a + b + c + d+  fsumL (toList -> l) = go l 0+    where+      go [] acc = acc+      go (x : xs) acc = go xs (acc + x)+  fdot3 a1 a2 a3 b1 b2 b3 = (a1 * b1) + (a2 * b2) + (a3 * b3)+  fdot4 a0 a1 a2 a3 b0 b1 b2 b3 = (a0 * b0) + (a1 * b1) + (a2 * b2) + (a3 * b3)+  fdotL (toList -> l1) (toList -> l2) = go l1 l2 0+    where+      go [] [] acc = acc+      go []  _  _  = error "Lists not the same length"+      go _  []  _  = error "Lists not the same length"+      go (b : bs) (c : cs) acc = go bs cs (fma b c acc)+--++--+instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => FusedOps (Posit es) where+  -- Fused Subtract Multiply+  fsm = viaRational3 fsm+  -- Fuse Multiply Add+  fma = viaRational3 fma+  -- Fuse Add Multiply+  fam = viaRational3 fam+  -- Fuse Multiply Multiply Subtract+  fmms = viaRational4 fmms+  -- Fuse Sum of 3 Posits+  fsum3 = viaRational3 fsum3+  -- Fuse Sum of 4 Posits+  fsum4 = viaRational4 fsum4+  -- Fuse Sum of a List+  fsumL (toList -> l) = Posit $ encode (Just $ go l 0)+    where+      go :: [Posit es] -> Rational -> Rational+      go [] !acc = acc+      go ((Posit int) : xs) !acc = case decode int of+                                     Nothing -> error "Posit List contains NaR"+                                     Just r -> go xs (acc + r)+  -- Fuse Dot Product of a 3-Vector+  fdot3 = viaRational6 fdot3+  -- Fuse Dot Product of a 4-Vector+  fdot4 = viaRational8 fdot4+  -- Fuse Dot Product of two Lists+  fdotL (toList -> l1) (toList -> l2) = Posit $ encode (Just $ go l1 l2 0)+    where+      go [] [] !acc = acc+      go []  _   _  = error "Lists not the same length"+      go _  []   _  = error "Lists not the same length"+      go ((Posit int1) : bs) ((Posit int2) : cs) !acc = case decode int1 of+                                                          Nothing -> error "First Posit List contains NaR"+                                                          Just r1 -> case decode int2 of+                                                                       Nothing -> error "Second Posit List contains NaR"+                                                                       Just r2 -> go bs cs (acc + (r1 * r2))+--+++++-- =====================================================================+-- ===                  Conversion Between Posits Types              ===+-- =====================================================================++-- |A Convertible class that will cast or 'convert' between two different Posit es types+class Convertible a b where+  convert :: a -> b++instance forall es1 es2. (PositC es1, PositC es2, Ord (IntN es1), Ord (IntN es2), Num (IntN es1), Num (IntN es2)) => Convertible (Posit es1) (Posit es2) where+  convert NaR = NaR+  convert (R r) = R r+--+++#ifndef O_NO_SHOW+-- =====================================================================+-- ===                Alternative Show Formats                       ===+-- =====================================================================++-- |A Alternative to the typical 'Show' class to assist in displaying the Posit es type in different formats+class AltShow a where+  -- |Display the Posit in its Binary Representation+  displayBinary :: a -> String+  -- |Display the Posit in its Integral Representation+  displayIntegral :: a -> String+  -- |Display the Posit as a Rational+  displayRational :: a -> String+  -- |Display the Posit as a Decimal until the Repetend occurs+  displayDecimal :: a -> String+--++--+instance forall es. (Show (IntN es), Ord (IntN es), Num (IntN es), PositC es) => AltShow (Posit es) where+  displayBinary (Posit int) = displayBin int+ +  displayIntegral (Posit int) = show int+ +  displayRational = viaShowable id+ +  displayDecimal = viaShowable (fst.fromRationalRepetendUnlimited)+--++viaShowable :: forall es a. (Show a, Ord (IntN es), Num (IntN es), PositC es) => (Rational -> a) -> Posit es -> String+viaShowable _ NaR = "NaR"+viaShowable f (R r) = show $ f r+#endif++#ifndef O_NO_READ+-- =====================================================================+-- ===                         Read Posit                            ===+-- =====================================================================++--+instance forall es. (PositC es) => Read (Posit es) where+  readPrec =+    parens $ do+      x <- lexP+      case x of+        Ident "NaR" -> return NaR+        _ -> pfail+      ++++      do+        s <- lift scientificP+        return $ R (toRational s)+ +  readListPrec = readListPrecDefault+--+#endif+++-- =====================================================================+-- ===                  Storable Instances                           ===+-- =====================================================================+--+#ifndef O_NO_STORABLE+--+instance forall es. (Storable (IntN es), PositC es) => Storable (Posit es) where+  sizeOf _ = fromIntegral $ nBytes @es+  alignment _ = fromIntegral $ nBytes @es+  peek ptr = do+    int <- peek (castPtr ptr :: Ptr (IntN es))+    return $ Posit int+  poke ptr (Posit int) = do+    poke (castPtr ptr :: Ptr (IntN es)) int+--+#endif+++-- =====================================================================+-- ===                        Real Frac                              ===+-- =====================================================================++--+instance forall es. (Num (IntN es), Ord (IntN es), PositC es) => RealFrac (Posit es) where+  -- properFraction :: Integral b => a -> (b, a)+  properFraction = viaRationalErrTrunkation "NaR value is not a RealFrac" properFraction+--++viaRationalErrTrunkation :: forall es a. (Num (IntN es), (Ord (IntN es)), PositC es, Integral a) => String -> (Rational -> (a, Rational)) -> Posit es -> (a, Posit es)+viaRationalErrTrunkation err _ NaR = error err+viaRationalErrTrunkation _ f (R r) =+  let (int, r') = f r+  in (int, R r')++-- =====================================================================+-- ===                         Real Float                            ===+-- =====================================================================+--+instance forall es. (Eq (IntN es), Ord (IntN es), Num (IntN es), Floating (Posit es), PositC es) => RealFloat (Posit es) where+  isIEEE _ = False+  isDenormalized _ = False+  isNegativeZero _ = False+ +  isNaN NaR = True+  isNaN  _  = False+ +  isInfinite NaR = True+  isInfinite _ = False+ +  -- 'atan2' of y x is the argument "arg function" (also called phase or angle) of the complex number x + i y.+  -- angle from an x basis vector to some other vector+  --+  --     Y+  --     ^+  --     |    ^ (x,y)+  --     |   /+  --     |  / <-  alpha (radians)+  --     | /                      \+  --      /                        |+  --      -----------------------------------> X+  --+  --+  atan2 NaR  _  = NaR+  atan2  _  NaR = NaR+  atan2 y x+    | x == 0 && y == 0 = NaR+    | x > 0             = atan (y/x)+    | x < 0  && y >= 0  = atan (y/x) + pi+    | x < 0  && y  < 0  = atan (y/x) - pi+    | x == 0 && y  > 0  = pi / 2+    | x == 0 && y  < 0  = negate $ pi / 2+    | otherwise = error "What!?!?!" -- The case where x == 0 && y == 0+ +  floatRadix _ = 2+  floatDigits _ = undefined+  floatRange _ = (negate maxExponent, maxExponent)+    where+      maxExponent = fromIntegral $ (nBytes @es) * ((nBits @es) - 2)+  decodeFloat = undefined+  encodeFloat = undefined+--++++-- =====================================================================+-- ===                         Floating                              ===+-- =====================================================================+++instance Floating Posit8 where+  pi = convert (pi :: Posit256) :: Posit8+  exp x = convert (exp (convert x) :: Posit256) :: Posit8+  log x = convert (log (convert x) :: Posit256) :: Posit8+  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit8+  sin x = convert (sin (convert x) :: Posit256) :: Posit8+  cos x = convert (cos (convert x) :: Posit256) :: Posit8+  asin x = convert (asin (convert x) :: Posit256) :: Posit8+  acos x = convert (acos (convert x) :: Posit256) :: Posit8+  atan x = convert (atan (convert x) :: Posit256) :: Posit8+  sinh x = convert (sinh (convert x) :: Posit256) :: Posit8+  cosh x = convert (cosh (convert x) :: Posit256) :: Posit8+  asinh x = convert (asinh (convert x) :: Posit256) :: Posit8+  acosh x = convert (acosh (convert x) :: Posit256) :: Posit8+  atanh x = convert (atanh (convert x) :: Posit256) :: Posit8++instance Floating Posit16 where+  pi = convert (pi :: Posit256) :: Posit16+  exp x = convert (exp (convert x) :: Posit256) :: Posit16+  log x = convert (log (convert x) :: Posit256) :: Posit16+  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit16+  sin x = convert (sin (convert x) :: Posit256) :: Posit16+  cos x = convert (cos (convert x) :: Posit256) :: Posit16+  asin x = convert (asin (convert x) :: Posit256) :: Posit16+  acos x = convert (acos (convert x) :: Posit256) :: Posit16+  atan x = convert (atan (convert x) :: Posit256) :: Posit16+  sinh x = convert (sinh (convert x) :: Posit256) :: Posit16+  cosh x = convert (cosh (convert x) :: Posit256) :: Posit16+  asinh x = convert (asinh (convert x) :: Posit256) :: Posit16+  acosh x = convert (acosh (convert x) :: Posit256) :: Posit16+  atanh x = convert (atanh (convert x) :: Posit256) :: Posit16++instance Floating Posit32 where+  pi = convert (pi :: Posit256) :: Posit32+  exp x = convert (exp (convert x) :: Posit256) :: Posit32+  log x = convert (log (convert x) :: Posit256) :: Posit32+  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit32+  sin x = convert (sin (convert x) :: Posit256) :: Posit32+  cos x = convert (cos (convert x) :: Posit256) :: Posit32+  asin x = convert (asin (convert x) :: Posit256) :: Posit32+  acos x = convert (acos (convert x) :: Posit256) :: Posit32+  atan x = convert (atan (convert x) :: Posit256) :: Posit32+  sinh x = convert (sinh (convert x) :: Posit256) :: Posit32+  cosh x = convert (cosh (convert x) :: Posit256) :: Posit32+  asinh x = convert (asinh (convert x) :: Posit256) :: Posit32+  acosh x = convert (acosh (convert x) :: Posit256) :: Posit32+  atanh x = convert (atanh (convert x) :: Posit256) :: Posit32++instance Floating Posit64 where+  pi = convert (pi :: Posit256) :: Posit64+  exp x = convert (exp (convert x) :: Posit256) :: Posit64+  log x = convert (log (convert x) :: Posit256) :: Posit64+  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit64+  sin x = convert (sin (convert x) :: Posit256) :: Posit64+  cos x = convert (cos (convert x) :: Posit256) :: Posit64+  asin x = convert (asin (convert x) :: Posit256) :: Posit64+  acos x = convert (acos (convert x) :: Posit256) :: Posit64+  atan x = convert (atan (convert x) :: Posit256) :: Posit64+  sinh x = convert (sinh (convert x) :: Posit256) :: Posit64+  cosh x = convert (cosh (convert x) :: Posit256) :: Posit64+  asinh x = convert (asinh (convert x) :: Posit256) :: Posit64+  acosh x = convert (acosh (convert x) :: Posit256) :: Posit64+  atanh x = convert (atanh (convert x) :: Posit256) :: Posit64++instance Floating Posit128 where+  pi = convert (pi :: Posit256) :: Posit128+  exp x = convert (exp (convert x) :: Posit256) :: Posit128+  log x = convert (log (convert x) :: Posit256) :: Posit128+  x ** y = convert $ (convert x :: Posit256) ** (convert y :: Posit256) :: Posit128+  sin x = convert (sin (convert x) :: Posit256) :: Posit128+  cos x = convert (cos (convert x) :: Posit256) :: Posit128+  asin x = convert (asin (convert x) :: Posit256) :: Posit128+  acos x = convert (acos (convert x) :: Posit256) :: Posit128+  atan x = convert (atan (convert x) :: Posit256) :: Posit128+  sinh x = convert (sinh (convert x) :: Posit256) :: Posit128+  cosh x = convert (cosh (convert x) :: Posit256) :: Posit128+  asinh x = convert (asinh (convert x) :: Posit256) :: Posit128+  acosh x = convert (acosh (convert x) :: Posit256) :: Posit128+  atanh x = convert (atanh (convert x) :: Posit256) :: Posit128++instance Floating Posit256 where+  pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286 :: Posit256+  exp = funExp+  log = funLogDomainReduction funLogTaylor+  (**) = funPow+  sin = funSin+  cos = funCos+  asin = funAsin+  acos = funAcos+  atan = funAtan+  sinh = funSinh+  cosh = funCosh+  asinh = funAsinh+  acosh = funAcosh+  atanh = funAtanh++++++class AltFloating p where+  phi :: p+  gamma :: p -> p+  sinc :: p -> p+  expm1 :: p -> p++instance AltFloating Posit8 where+  phi = convert (phi :: Posit256) :: Posit8+  gamma x = convert (gamma (convert x) :: Posit256) :: Posit8+  sinc x = convert (sinc (convert x) :: Posit256) :: Posit8+  expm1 x =+    let b = atanh $ x / 2+    in (2 * b) / (1 - b)++instance AltFloating Posit16 where+  phi = convert (phi :: Posit256) :: Posit16+  gamma x = convert (gamma (convert x) :: Posit256) :: Posit16+  sinc x = convert (sinc (convert x) :: Posit256) :: Posit16+  expm1 x =+    let b = atanh $ x / 2+    in (2 * b) / (1 - b)++instance AltFloating Posit32 where+  phi = convert (phi :: Posit256) :: Posit32+  gamma x = convert (gamma (convert x) :: Posit256) :: Posit32+  sinc x = convert (sinc (convert x) :: Posit256) :: Posit32+  expm1 x =+    let b = atanh $ x / 2+    in (2 * b) / (1 - b)++instance AltFloating Posit64 where+  phi = convert (phi :: Posit256) :: Posit64+  gamma x = convert (gamma (convert x) :: Posit256) :: Posit64+  sinc x = convert (sinc (convert x) :: Posit256) :: Posit64+  expm1 x =+    let b = atanh $ x / 2+    in (2 * b) / (1 - b)++instance AltFloating Posit128 where+  phi = convert (phi :: Posit256) :: Posit128+  gamma x = convert (gamma (convert x) :: Posit256) :: Posit128+  sinc x = convert (sinc (convert x) :: Posit256) :: Posit128+  expm1 x =+    let b = atanh $ x / 2+    in (2 * b) / (1 - b)++instance AltFloating Posit256 where+  phi = funPhi 1.6+  gamma = funGammaSeries+  sinc = funSinc+  expm1 x =+    let b = atanh $ x / 2+    in (2 * b) / (1 - b)+++-- | 'phi' fixed point recursive algorithm,+funPhi :: Posit256 -> Posit256+funPhi  px@(Posit x)+    | x == x' = Posit x+    | otherwise = funPhi (Posit x')+      where+        (Posit x') = (px^2 + 2*px) / (px^2 + 1)+++-- calculate atan(1/2^n)+-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached+funArcTan :: Natural -> Posit256+funArcTan 0 = pi / 4+funArcTan n+  | n <= 122 = go 0 0+  | otherwise = z  -- at small z... (atan z) == z "small angle approximation"+    where+      go !k !acc+        | acc == (acc + term k) = acc+        | otherwise = go (k+1) (acc + term k)+      term :: Integer -> Posit256+      term k = ((-1)^k * z^(2 * k + 1)) / fromIntegral (2 * k + 1)+      z = 1 / 2^n  -- recip $ 2^n :: Posit256 -- inv2PowN++-- seems pretty close to 1 ULP with the input of 0.7813+funAtan :: Posit256 -> Posit256+funAtan NaR = NaR+funAtan x+  | abs x < 1/2^122 = x  -- small angle approximaiton, found emperically+  | x < 0 = negate.funAtan $ negate x  -- if negative turn it positive, it reduces the other domain reductions by half, found from Universal CORDIC+  | x > 1 = pi/2 - funAtan (recip x)  -- if larger than one use the complementary angle, found from Universal CORDIC+  | x > twoMsqrt3 = pi/6 + funAtan ((sqrt 3 * x - 1)/(sqrt 3 + x))  -- another domain reduction, using an identity, found from https://mathonweb.com/help_ebook/html/algorithms.htm+  | otherwise = funArcTanTaylor x+--++twoMsqrt3 :: Posit256+twoMsqrt3 = 2 - sqrt 3++--+funArcTanTaylor :: Posit256 -> Posit256+funArcTanTaylor x = go 0 0+  where+    go !k !acc+      | acc == (acc + term k) = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Integer -> Posit256+    term k = ((-1)^k * x^(2 * k + 1)) / fromIntegral (2 * k + 1)+--++--+funAsin :: Posit256 -> Posit256+funAsin NaR = NaR+funAsin x+  | abs x > 1 = NaR+  | x == 1 = pi/2+  | x == -1 = -pi/2+  | otherwise = funAtan w+    where+      w = x / sqrt (1 - x^2)+--++--+funAcos :: Posit256 -> Posit256+funAcos NaR = NaR+funAcos x+  | abs x > 1 = NaR+  | x < 0 = pi + funAtan invw+  | x == 0 = pi/2+  | x > 0 = funAtan invw+  | otherwise = error "Prove it covers for Rational Numbers."+    where+      invw = sqrt (1 - x^2) / x+--++-- fI2PN = (1 /) . (2 ^)+funInv2PowN :: Natural -> Posit256+funInv2PowN n = 1 / 2^n+++-- calculate atanh(1/2^n)+-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached+funArcHypTan :: Natural -> Posit256+funArcHypTan 0 = NaR+funArcHypTan n+  | n <= 122 = go 0 0+  | otherwise = z  -- at small z... (atan z) == z "small angle approximation"+    where+      go !k !acc+        | acc == (acc + term k) = acc+        | otherwise = go (k+1) (acc + term k)+      term :: Integer -> Posit256+      term k = (z^(2 * k + 1)) / fromIntegral (2 * k + 1)+      z = 1 / 2^n+++fac :: Natural -> Natural+fac 0 = 1+fac n = n * fac (n - 1)++--+funAsinh :: Posit256 -> Posit256+funAsinh NaR = NaR+funAsinh x = log $ x + sqrt (x^2 + 1)+--++--+funAcosh :: Posit256 -> Posit256+funAcosh NaR = NaR+funAcosh x+  | x < 1 = NaR+  | otherwise = log $ x + sqrt (x^2 - 1)+--++--+funAtanh :: Posit256 -> Posit256+funAtanh NaR = NaR+funAtanh x+  | abs x >= 1 = NaR+  | x < 0 = negate.funAtanh.negate $ x  -- make use of odd parity to only calculate the positive part+  | otherwise = 0.5 * log ((1+t) / (1-t)) - (fromIntegral ex / 2) * lnOf2+    where+      (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat + 1, fromRational rat / 2)+      (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x' -- sign should always be positive+      x' = 1 - x+      t = (2 - sig - x') / (2 + sig - x')+--++--+funAtanhTaylor :: Posit256 -> Posit256+funAtanhTaylor NaR = NaR+funAtanhTaylor x+  | abs x >= 1 = NaR+  | abs x < 1/2^122 = x  -- small angle approximaiton, found emperically+  | x < 0 = negate.funAtanhTaylor.negate $ x+  | otherwise = go 0 0+    where+      go !k !acc+        | acc == (acc + term k) = acc+        | otherwise = go (k+1) (acc + term k)+      term :: Integer -> Posit256+      term k = (x^(2 * k + 1)) / fromIntegral (2 * k + 1)+--++--+funSin :: Posit256 -> Posit256+funSin NaR = NaR+funSin 0 = 0+funSin x = funSin' $ x / (2*pi)+--+-- funSin' is sine normalized by 2*pi+funSin' :: Posit256 -> Posit256+funSin' x+  | x == 0 = 0+  | x == 0.25 = 1+  | x == 0.5 = 0+  | x == 0.75 = -1+  | x == 1 = 0+  | x < 0 = negate.funSin'.negate $ x+  | x > 1 =+    let (_,rem) = properFraction x+    in funSin' rem+  | x > 0.75 && x < 1 = negate.funSin' $ 1 - x -- reduce domain by quadrant symmetry+  | x > 0.5 && x < 0.75 = negate.funSin' $ x - 0.5+  | x > 0.25 && x < 0.5 = funSin' $ 0.5 - x+  | x > 0.125 && x < 0.25 = funCosTuma $ 2*pi * (0.25 - x) -- reduce domain and use cofunction+  | otherwise = funSinTuma $ 2*pi * x+--++-- Taylor series expansion and fixed point algorithm, most accurate near zero+funSinTaylor :: Posit256 -> Posit256+funSinTaylor NaR = NaR+funSinTaylor z = go 0 0+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | acc == (acc + term k) = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Natural -> Posit256+    term k = (-1)^k * z^(2*k+1) / (fromIntegral.fac $ 2*k+1)+--++--+funSinTuma :: Posit256 -> Posit256+funSinTuma NaR = NaR+funSinTuma z = go 19 1+  where+    go :: Natural -> Posit256 -> Posit256+    go 1 !acc = z * acc+    go !k !acc = go (k-1) (1 - (z^2 / fromIntegral ((2*k-2)*(2*k-1))) * acc)+--++--+funCos :: Posit256 -> Posit256+funCos NaR = NaR+funCos 0 = 1+funCos x = funCos' $ x / (2*pi)+--+-- funCos' is cosine normalized for 2*pi+funCos' :: Posit256 -> Posit256+funCos' NaR = NaR+funCos' x+  | x == 0 = 1+  | x == 0.25 = 0+  | x == 0.5 = -1+  | x == 0.75 = 0+  | x == 1 = 1+  | x < 0 = funCos'.negate $ x  -- reduce domain by symmetry across 0 to turn x positive+  | x > 1 = -- reduce domain by using perodicity+    let (_,rem) = properFraction x+    in funCos' rem+  | x > 0.75 && x < 1 = funCos' $ 1 - x  -- reduce domain by quadrant symmetry+  | x > 0.5 && x < 0.75 = negate.funCos' $ x - 0.5+  | x > 0.25 && x < 0.5 = negate.funCos' $ 0.5 - x+  | x > 0.125 && x < 0.25 = funSinTuma $ 2*pi * (0.25 - x) -- reduce domain and use cofunction+  | otherwise = funCosTuma $ 2*pi * x --+--++-- Taylor series expansion and fixed point algorithm, most accurate near zero+funCosTaylor :: Posit256 -> Posit256+funCosTaylor NaR = NaR+funCosTaylor z = go 0 0+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | acc == (acc + term k) = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Natural -> Posit256+    term k = (-1)^k * z^(2*k) / (fromIntegral.fac $ 2*k)+--++--+funCosTuma :: Posit256 -> Posit256+funCosTuma NaR = NaR+funCosTuma z = go 19 1+  where+    go :: Natural -> Posit256 -> Posit256+    go 1 !acc = acc+    go !k !acc = go (k-1) (1 - (z^2 / fromIntegral ((2*k-3)*(2*k-2))) * acc)+--++-- ~16 ULP for 42+funSinh :: Posit256 -> Posit256+funSinh NaR = NaR+funSinh x = (exp x - exp (negate x))/2+--++-- ~2 ULP for 42+funSinhTaylor :: Posit256 -> Posit256+funSinhTaylor NaR = NaR+funSinhTaylor z = go 0 0+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | acc == (acc + term k) = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Natural -> Posit256+    term k = z^(2*k+1) / (fromIntegral.fac $ 2*k+1)+--++--+funSinhTuma :: Posit256 -> Posit256+funSinhTuma NaR = NaR+funSinhTuma 0 = 0+funSinhTuma z | z < 0 = negate.funSinhTuma.negate $ z+funSinhTuma z | z > 80 = 0.5 * funExpTuma z+funSinhTuma z = go 256 1+  where+    go :: Natural -> Posit256 -> Posit256+    go 1 !acc = z * acc+    go !k !acc = go (k-1) (1 + (z^2 / fromIntegral ((2*k-2) * (2*k-1))) * acc)+--++-- ~17 ULP for 42+funCosh :: Posit256 -> Posit256+funCosh NaR = NaR+funCosh x = (exp x + exp (negate x))/2+--++-- ~3 ULP for 42+funCoshTaylor :: Posit256 -> Posit256+funCoshTaylor NaR = NaR+funCoshTaylor z = go 0 0+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | acc == (acc + term k) = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Natural -> Posit256+    term k = z^(2*k) / (fromIntegral.fac $ 2*k)+--++--+funCoshTuma :: Posit256 -> Posit256+funCoshTuma NaR = NaR+funCoshTuma 0 = 1+funCoshTuma z | z < 0 = funCoshTuma.negate $ z+funCoshTuma z | z > 3 = 0.5 * (funExpTuma z + funExpTuma (negate z))+funCoshTuma z = go 20 1+  where+    go :: Natural -> Posit256 -> Posit256+    go 1 !acc = acc+    go !k !acc = go (k-1) (1 + (z^2 / fromIntegral ((2*k-3)*(2*k-2)))*acc)+--+++--+funLog :: Posit256 -> Posit256+funLog x = funLog2 x * lnOf2+--++--+-- Use the constant, for performance+lnOf2 :: Posit256+lnOf2 = Posit 28670435363615573179632300308403400109260626501925370561166468529302554498548+--++--+-- Some series don't converge reliably, this one does+funLnOf2 :: Posit256+funLnOf2 = go 1 0+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | acc == (acc + term k) = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Natural -> Posit256+    term k = 1 / fromIntegral (2^k * k)+--++--+funLog2 :: Posit256 -> Posit256+funLog2 NaR = NaR+funLog2 z+  | z <= 0 = NaR -- includes the NaR case+  | otherwise = go (fromInteger ex) 1 sig  -- domain reduction+    where+      go :: Posit256 -> Posit256 -> Posit256 -> Posit256+      go !acc !mak !sig' -- fixed point iteration, y is [1,2) :: Posit256+        | sig == 1 = acc+        | acc == (acc + mak * 2^^(negate.fst.term $ sig')) = acc  -- stop when fixed point is reached+        | otherwise = go (acc + mak * 2^^(negate.fst.term $ sig')) (mak * 2^^(negate.fst.term $ sig')) (snd.term $ sig')+      term = findSquaring 0  -- returns (m,s') m the number of times to square, and the new significand+      (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat, fromRational rat)+      (_,int,nat,rat) = (posit2TupPosit @V).toRational $ z -- sign should always be positive+      findSquaring m s+        | s >= 2 && s < 4 = (m, s/2)+        | otherwise = findSquaring (m+1) (s^2)+--+++--  Gauss–Legendre algorithm, Seems only accurate to 2-3 ULP, but really slow+funPi1 :: Posit256+funPi1 = go 0 3 1 (recip.sqrt $ 2) (recip 4) 1+  where+    go :: Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256 -> Posit256+    go !prev !next !a !b !t !p+      | prev == next = next+      | otherwise =+        let a' = (a + b) / 2+            b' = sqrt $ a * b+            t' = t - p * (a - ((a + b) / 2))^2+            p' = 2 * p+        in go next ((a' + b')^2 / (4 * t')) a' b' t' p'+--++--  Borwein's algorithm, with quintic convergence,+--  gets to 7 ULP in 4 iterations, but really slow due to expensive function evaluations+--  quite unstable and will not converge if sqrt is not accurate, which means log must be accurate+funPi2 :: Posit256+funPi2 = recip $ go 0 0 0.5 (5 / phi^3)+  where+    go :: Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256+    go !prev !n !a !s+      | prev == a = a+      | otherwise =+        let x = 5 / s - 1+            y = (x - 1)^2 + 7+            z = (0.5 * x * (y + sqrt (y^2 - 4 * x^3)))**(1/5)+            a' = s^2 * a - (5^n * ((s^2 - 5)/2 + sqrt (s * (s^2 - 2*s + 5))))+            s' = 25 / ((z + x/z + 1)^2 * s)+        in go a (n+1) (trace (show a') a') s'+--+++-- Bailey–Borwein–Plouffe (BBP) formula, to 1-2 ULP, and blazing fast, converges in 60 iterations+funPi3 :: Posit256+funPi3 = go 0 0+  where+    go :: Integer -> Posit256 -> Posit256+    go !k !acc+      | acc == acc + term k = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Integer -> Posit256+    term k = fromRational $ (1 % 16^k) * ((120 * k^2 + 151 * k + 47) % (512 * k^4 + 1024 * k^3 + 712 * k^2 + 194 * k + 15))+--+++-- Fabrice Bellard improvement on the BBP, 2-3 ULP, even faster, converges in 25 iterations, really fast+funPi4 :: Posit256+funPi4 = (1/2^6) * go 0 0+  where+    go :: Integer -> Posit256 -> Posit256+    go !k !acc+      | acc == acc + term k = acc+      | otherwise = go (k+1) (acc + term k)+    term :: Integer -> Posit256+    term k = fromRational $ ((-1)^k % (2^(10*k))) * ((1 % (10 * k + 9)) - (2^2 % (10 * k + 7)) - (2^2 % (10 * k + 5)) - (2^6 % (10 * k + 3)) + (2^8 % (10 * k + 1)) - (1 % (4 * k + 3)) - (2^5 % (4 * k + 1)))+--++++--+-- looks to be about 4 ULP accurate at -100, right on the money at -1000+funExp :: Posit256 -> Posit256+funExp x = funExp2 funExpTaylor (x / lnOf2)+--++--+--+funExp2 :: (Posit256 -> Posit256) -> Posit256 -> Posit256+funExp2 _ NaR = NaR+funExp2 _ 0 = 1+funExp2 f x+  | x < 0 = recip.funExp2 f.negate $ x  -- always calculate the positive method+  | otherwise = case properFraction x of+                  (int,rem) -> fromIntegral (2^int) * f (lnOf2 * rem)++++--+-- calculate exp, its most accurate near zero+-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached+funExpTaylor :: Posit256 -> Posit256+funExpTaylor NaR = NaR+funExpTaylor 0 = 1+funExpTaylor z = go 0 0+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | acc == (acc + term k) = acc  -- if x == x + dx then terminate and return x+      | otherwise = go (k+1) (acc + term k)+    term :: Natural -> Posit256+    term k = (z^k) / (fromIntegral.fac $ k)+--++--+-- calculate exp, its most accurate near zero+-- use the Nested Series of Jan J Tuma+funExpTuma :: Posit256 -> Posit256+funExpTuma NaR = NaR+funExpTuma 0 = 1+funExpTuma z = go 57 1 -- was 66+  where+    go :: Natural -> Posit256 -> Posit256+    go !k !acc+      | k == 0 = acc+      | otherwise = go (k-1) (1 + (z / fromIntegral k) * acc)+--++--+--+funPow :: Posit256 -> Posit256 -> Posit256+NaR `funPow` _ = NaR+_ `funPow` NaR = NaR+funPow 0 y+  | y < 0 = NaR -- NaR: Divide by Zero+  | y == 0 = NaR -- NaR: Indeterminate+  | y > 0 = 0+funPow x y+  | y < 0 = recip $ funPow x (negate y)+  | x < 0 = -- NaR if y is not an integer+    let (int,rem) = properFraction y+    in if rem == 0+       then x^^int+       else NaR -- NaR: Imaginary Number+  | otherwise = exp $ y * log x+--++-- Looks like 1 ULP for 0.7813+funSinc :: Posit256 -> Posit256+funSinc NaR = NaR+funSinc 0 = 1  -- Why the hell not!+funSinc theta = sin theta / theta+--++-- Interestingly enough, wikipedia defines two alternative solutions+-- for the Shannon Wavelet, eventhough there are infinite solutions+-- where the functions are equal, they are not equal+-- Shannon wavelet+funPsiSha1 :: Posit256 -> Posit256+funPsiSha1 NaR = NaR+funPsiSha1 t = 2 * funSinc (2 * t) - funSinc t+--++-- Shannon wavelet+funPsiSha2 :: Posit256 -> Posit256+funPsiSha2 NaR = NaR+funPsiSha2 t = funSinc (t/2) * cos (3*pi*t/2)+--++-- Shannon wavelet+funPsiSha3 :: Posit256 -> Posit256+funPsiSha3 NaR = NaR+funPsiSha3 0 = 1  -- Why the hell not!+funPsiSha3 t =+  let pit = pi * t+      invpit = recip $ pit +  in invpit * (sin (2 * pit) - sin pit)+--++-- funPsiSha1.(pi*) === funPsiSha3++-- Using the CORDIC domain reduction and some approximation function+funLogDomainReduction :: (Posit256 -> Posit256) -> Posit256 -> Posit256+funLogDomainReduction _ NaR = NaR+funLogDomainReduction _ 1 = 0+funLogDomainReduction f x+  | x <= 0 = NaR+  | otherwise = f sig + (fromIntegral ex * lnOf2)+    where+      (ex, sig) = (int * fromIntegral (nBytes @V) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1+      (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x -- sign should always be positive+     + ++-- natural log with log phi acurate to 9 ULP+funLogTaylor :: Posit256 -> Posit256+funLogTaylor NaR = NaR+funLogTaylor 1 = 0+funLogTaylor x | x <= 0 = NaR+funLogTaylor x+  | x <= 2 = go 1 0+  | otherwise = error "The funLogTaylor algorithm is being used improperly"+    where+      go :: Natural -> Posit256 -> Posit256+      go !k !acc+        | acc == (acc + term k) = acc+        | otherwise = go (k + 1) (acc + term k)+      term :: Natural -> Posit256+      term k = (-1)^(k+1) * (x - 1)^k / fromIntegral k+     ++++-- natural log the Jan J Tuma way+funLogTuma :: Posit256 -> Posit256+funLogTuma NaR = NaR+funLogTuma 1 = 0  -- domain reduced input is [0.5,1) and/or , where funLogTuma 1 = 0+funLogTuma x | x <= 0 = NaR  -- zero and less than zero is NaR+funLogTuma x+  = go 242 1+    where+      xM1 = x - 1  -- now [-0.5, 0)+      go :: Natural -> Posit256 -> Posit256+      go !k !acc+        | k == 0 = xM1 * acc+        | otherwise = go (k-1) (recip (fromIntegral k) - xM1 * acc)+++funGammaRamanujan :: Posit256 -> Posit256+funGammaRamanujan z = sqrt pi * (x / exp 1)**x * (8*x^3 + 4*x^2 + x + (1/30))**(1/6)+  where+    x = z - 1++--+a001163 :: [Integer] -- Numerator+a001163 = [1, 1, -139, -571, 163879, 5246819, -534703531, -4483131259, 432261921612371, 6232523202521089, -25834629665134204969, -1579029138854919086429, 746590869962651602203151, 1511513601028097903631961, -8849272268392873147705987190261, -142801712490607530608130701097701]+a001164 :: [Integer]  -- Denominator+a001164 = [12, 288, 51840, 2488320, 209018880, 75246796800, 902961561600, 86684309913600, 514904800886784000, 86504006548979712000, 13494625021640835072000, 9716130015581401251840000, 116593560186976815022080000, 2798245444487443560529920000, 299692087104605205332754432000000, 57540880724084199423888850944000000]++funGammaSeries :: Posit256 -> Posit256+funGammaSeries z = sqrt(2 * pi) * (z**(z - 0.5)) * exp (negate z) * (1 + series)+  where+    series :: Posit256+    series = sum $ zipWith (*) [fromRational (a % b) | (a,b) <- zip a001163 a001164] [recip $ z^n |  n <- [1..len]]  -- zipWith (\x y -> ) a001163 a001164+    lenA = length a001163+    lenB = length a001164+    len = if lenA == lenB+            then lenA+            else error "Seiries Numerator and Denominator do not have the same length."+--++funGammaCalc :: Posit256 -> Posit256+funGammaCalc z = sqrt (2*pi / z) * ((z / exp 1) * sqrt (z * sinh (recip z) + recip (810 * z^6)))**z+++funGammaNemes :: Posit256 -> Posit256+funGammaNemes z = sqrt (2*pi / z) * (recip (exp 1) * (z + recip (12 * z - recip (10 * z))))**z++funGammaYang :: Posit256 -> Posit256+funGammaYang z = sqrt (2 * pi * x) * (x / exp 1)**x * (x * sinh (recip x))**(x/2) * exp (fromRational (7 % 324) * recip (x^3 * (35 * x^2 + 33)))+  where+    x = z - 1++funGammaChen :: Posit256 -> Posit256+funGammaChen z = sqrt (2 * pi * x) * (x / exp 1)**x * (1 + recip (12*x^3 + (24/7)*x - 0.5))**(x^2 + fromRational (53 % 210))+  where+    x = z - 1++

@@ -0,0 +1,983 @@++--------------------------------------------------------------------------------------------+--+--   Copyright   :  (C) 2022 Nathan Waivio+--   License     :  BSD3+--   Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>+--   Stability   :  Stable+--   Portability :  Portable+--+-- | Library implementing standard 'Posit-3.2' numbers, as defined by+--   the Posit Working Group 23 June 2018.+-- +-- +---------------------------------------------------------------------------------------------+++{-# LANGUAGE TypeFamilyDependencies #-} -- For the associated bidirectional type family that the Posit library is based on+{-# LANGUAGE DataKinds #-}  -- For our ES kind and the constructors Z, I, II, III, IV, V, for exponent size type+{-# LANGUAGE TypeApplications #-}  -- The most excellent syntax @Int256+{-# LANGUAGE AllowAmbiguousTypes #-} -- The Haskell/GHC Type checker seems to have trouble things in the PositC class+{-# LANGUAGE ScopedTypeVariables #-} -- To reduce some code duplication+{-# LANGUAGE FlexibleContexts #-} -- To reduce some code duplication by claiming the type family provides some constraints, that GHC can't do without fully evaluating the type family+{-# LANGUAGE CPP #-} -- To remove Storable instances to remove noise when performing analysis of Core+{-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-}  -- Turn off noise+{-# OPTIONS_GHC -Wno-type-defaults #-}  -- Turn off noise++-- ----+--  |Posit Class, implementing:+--+--   * PositC+--   * Orphan Instances of Storable for Word128, Int128, Int256+-- ----++module Posit.Internal.PositC+(PositC(..),+ ES(..)+ ) where++import Prelude hiding (exponent,significand)++-- Imports for Storable Instance of Data.DoubleWord+import Foreign.Storable (Storable, sizeOf, alignment, peek, poke)  -- Used for Storable Instances of Data.DoubleWord+import Foreign.Ptr (Ptr, plusPtr, castPtr)  -- Used for dealing with Pointers for the Data.DoubleWord Storable Instance++-- Machine Integers and Operations+{-@ embed Int128 * as int @-}+{-@ embed Int256 * as int @-}+import Data.Int (Int8,Int16,Int32,Int64)  -- Import standard Int sizes+import Data.DoubleWord (Word128,Int128,Int256,fromHiAndLo,hiWord,loWord) -- Import large Int sizes+import Data.Word (Word64)+import Data.Bits ((.|.), shiftL, shift, testBit, (.&.), shiftR)++-- Import Naturals and Rationals+{-@ embed Natural * as int @-}+import GHC.Natural (Natural) -- Import the Natural Numbers ℕ (u+2115)+{-@ embed Ratio * as int @-}+import Data.Ratio (Rational, (%))  -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D)+++-- | The Exponent Size 'ES' kind, the constructor for the Type is a Roman Numeral.+data ES = Z+        | I+        | II+        | III+        | IV+        | V+++-- | The 'Posit' class is an approximation of ℝ, it is like a sampling on the Projective Real line ℙ(ℝ) with Maybe ℚ as the internal type.+-- The 'es' is an index that controlls the log2 word size of the Posit's+-- fininte precision representation.+class PositC (es :: ES) where+  -- | Type of the Finite Precision Representation, in our case Int8, Int16, Int32, Int64, Int128, Int256. The 'es' of kind 'ES' will determine a result of 'r' such that you can determine the 'es' by the 'r'+  type IntN es = r | r -> es+ + +  -- | Transform to/from the Infinite Precision Representation+  encode :: Maybe Rational -> IntN es  -- ^ Maybe you have some Rational Number and you want to encode it as some integer with a finite integer log2 word size.+  decode :: IntN es -> Maybe Rational  -- ^ You have an integer with a finite integer log2 word size decode it and Maybe it is Rational+ +  -- | Exponent Size based on the Posit Exponent kind ES+  exponentSize :: Natural  -- ^ The exponent size, 'es' is a Natural number+ +  -- | Various other size definitions used in the Posit format with their default definitions+  nBytes :: Natural  -- ^ 'nBytes' the number of bytes of the Posit Representation+  nBytes = 2^(exponentSize @es)+ +  nBits :: Natural  -- ^ 'nBits' the number of bits of the Posit Representation+  nBits = 8 * (nBytes @es)+ +  signBitSize :: Natural  -- ^ 'signBitSize' the size of the sign bit+  signBitSize = 1+ +  uSeed :: Natural  -- ^ 'uSeed' scaling factor for the regime of the Posit Representation+  uSeed = 2^(nBytes @es)+ +  -- | Integer Representation of common bounds+  unReal :: IntN es  -- ^ 'unReal' is something that is not Real, the integer value that is not a Real number+ +  mostPosVal :: IntN es+  leastPosVal :: IntN es+  leastNegVal :: IntN es+  mostNegVal :: IntN es+ +  -- Rational Value of common bounds+  maxPosRat :: Rational+  maxPosRat = 2^((nBytes @es) * ((nBits @es) - 2)) % 1+  minPosRat :: Rational+  minPosRat = recip (maxPosRat @es)+  maxNegRat :: Rational+  maxNegRat = negate (minPosRat @es)+  minNegRat :: Rational+  minNegRat = negate (maxPosRat @es)+ +  -- Functions to support encode and decode+ +  -- log base uSeed+  -- After calculating the regime the rational should be in the range [1,uSeed), it starts with (0,rational)+  log_uSeed :: (Integer, Rational) -> (Integer, Rational)+  log_uSeed (regime,r)+    | r < 1 = log_uSeed @es (regime-1,r * fromRational (toInteger (uSeed @es) % 1))+    | r >= fromRational (toInteger (uSeed @es) % 1) = log_uSeed @es (regime+1,r * fromRational (1 % toInteger (uSeed @es)))+    | otherwise = (regime,r)+ +  getRegime :: Rational -> (Integer, Rational)+  getRegime r = log_uSeed @es (0,r)+ +  posit2TupPosit :: Rational -> (Bool, Integer, Natural, Rational)+  posit2TupPosit r =+    let (sgn,r') = getSign r -- returns the sign and a positive rational+        (regime,r'') = getRegime @es r' -- returns the regime and a rational between uSeed^-1 to uSeed^1+        (exponent,significand) = getExponent r'' -- returns the exponent and a rational between [1,2), the significand+    in (sgn,regime,exponent,significand)+ +  buildIntRep :: Rational -> IntN es+  mkIntRep :: Integer -> Natural -> Rational -> IntN es+  formRegime :: Integer -> (IntN es, Integer)+  formExponent :: Natural -> Integer -> (IntN es, Integer)+  formFraction :: Rational -> Integer -> IntN es+ +  tupPosit2Posit :: (Bool,Integer,Natural,Rational) -> Maybe Rational+  tupPosit2Posit (sgn,regime,exponent,rat) = -- s = isNeg posit == True+    let pow2 = toRational (uSeed @es)^^regime * 2^exponent+        scale = if sgn+                then negate pow2+                else pow2+    in Just $ scale * rat+ +  regime2Integer :: IntN es -> (Integer, Int)+  findRegimeFormat :: IntN es -> Bool+  countRegimeBits :: Bool -> IntN es -> Int+  exponent2Nat :: Int -> IntN es -> Natural+  fraction2Posit :: Int -> IntN es -> Rational+ +  -- prints out the IntN es value in 0b... format+  displayBin :: IntN es -> String+  -- decimal Precision+  decimalPrec :: Int+  decimalPrec = fromIntegral $ 2 * (nBytes @es) + 1++++instance PositC Z where+  type IntN Z = Int8+  exponentSize = 0+ +  -- Posit Integer Rep of various values+  unReal = minBound @Int8+ +  mostPosVal = maxBound @Int8+  leastPosVal = 1+  leastNegVal = -1+  mostNegVal = negate mostPosVal+ +  encode Nothing = unReal @Z+  encode (Just 0) = 0+  encode (Just r)+    | r > maxPosRat @Z = mostPosVal @Z+    | r < minNegRat @Z = mostNegVal @Z+    | r > 0 && r < minPosRat @Z = leastPosVal @Z+    | r < 0 && r > maxNegRat @Z = leastNegVal @Z+    | otherwise = buildIntRep @Z r+ +  buildIntRep r =+    let (signBit,regime,exponent,significand) = posit2TupPosit @Z r+        intRep = mkIntRep @Z regime exponent significand+    in if signBit+       then negate intRep+       else intRep+ +  mkIntRep regime exponent significand =+    let (regime', offset) = formRegime @Z regime  -- offset is the number of binary digits remaining after the regime is formed+        (exponent', offset') = formExponent @Z exponent offset  -- offset' is the number of binary digits remaining after the exponent is formed+        fraction = formFraction @Z significand offset'+    in regime' .|. exponent' .|. fraction+ +  formRegime power+    | 0 <= power =+      let offset = (fromIntegral (nBits @Z - 1) -     power - 1)+      in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1)+    | otherwise =+      let offset = (fromIntegral (nBits @Z - 1) - abs power - 1)+      in (1 `shiftL` fromInteger offset, offset)+ +  formExponent power offset =+    let offset' = offset - fromIntegral (exponentSize @Z)+    in (fromIntegral power `shift` fromInteger offset', offset')+ +  formFraction r offset =+    let numFractionBits = offset+        fractionSize = 2^numFractionBits+        normFraction = round $ (r - 1) * fractionSize  -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1)+    in if numFractionBits >= 1+       then fromInteger normFraction+       else 0+ +  decode int+    | int == unReal @Z = Nothing+    | int == 0 = Just 0+    | otherwise =+      let sgn = int < 0+          int' = if sgn+                 then negate int+                 else int+          (regime,nR) = regime2Integer @Z int'+          exponent = exponent2Nat @Z nR int'  -- if no e or some bits missing, then they are considered zero+          rat = fraction2Posit @Z nR int'  -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1+      in tupPosit2Posit @Z (sgn,regime,exponent,rat)+ +  regime2Integer posit =+    let regimeFormat = findRegimeFormat @Z posit+        regimeCount = countRegimeBits @Z regimeFormat posit+        regime = calcRegimeInt regimeFormat regimeCount+    in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime+ +  -- will return the format of the regime, either HI or LO; it could get refactored in the future+  -- True means a 1 is the first bit in the regime+  findRegimeFormat posit = testBit posit (fromIntegral (nBits @Z) - 1 - fromIntegral (signBitSize @Z))+ +  countRegimeBits format posit = go (fromIntegral (nBits @Z) - 1 - fromIntegral (signBitSize @Z)) 0+    where+      go (-1) acc = acc+      go index acc+        | xnor format (testBit posit index)  = go (index - 1) (acc + 1)+        | otherwise = acc+ +  -- knowing the number of the regime bits, and the sign bit we can extract+  -- the exponent.  We mask to the left of the exponent to remove the sign and regime, and+  -- then shift to the right to remove the fraction.+  exponent2Nat numBitsRegime posit =+    let bitsRemaining = fromIntegral (nBits @Z) - numBitsRegime - fromIntegral (signBitSize @Z)+        signNRegimeMask = 2^bitsRemaining - 1+        int = posit .&. signNRegimeMask+        nBitsToTheRight = fromIntegral (nBits @Z) - numBitsRegime - fromIntegral (signBitSize @Z) - fromIntegral (exponentSize @Z)+    in if bitsRemaining <=0+       then 0+       else if nBitsToTheRight < 0+            then fromIntegral $ int `shiftL` negate nBitsToTheRight+            else fromIntegral $ int `shiftR` nBitsToTheRight+ +  -- knowing the number of the regime bits, sign bit, and the number of the+  -- exponent bits we can extract the fraction.  We mask to the left of the fraction to+  -- remove the sign, regime, and exponent. If there is no fraction then the value is 1.+  fraction2Posit numBitsRegime posit =+    let offset = fromIntegral $ (signBitSize @Z) + fromIntegral numBitsRegime + (exponentSize @Z)+        fractionSize = fromIntegral (nBits @Z) - offset+        fractionBits = posit .&. (2^fractionSize - 1)+    in if fractionSize >= 1+       then (2^fractionSize + toInteger fractionBits) % 2^fractionSize+       else 1 % 1+ +  displayBin int = "0b" ++ go (fromIntegral (nBits @Z) - 1)+    where+      go :: Int -> String+      go 0 = if testBit int 0+             then "1"+             else "0"+      go idx = if testBit int idx+               then "1" ++ go (idx - 1)+               else "0" ++ go (idx -1)++++instance PositC I where+  type IntN I = Int16+  exponentSize = 1+ +  -- Posit Integer Rep of various values+  unReal = minBound @Int16+ +  mostPosVal = maxBound @Int16+  leastPosVal = 1+  leastNegVal = -1+  mostNegVal = negate mostPosVal+ +  encode Nothing = unReal @I+  encode (Just 0) = 0+  encode (Just r)+    | r > maxPosRat @I = mostPosVal @I+    | r < minNegRat @I = mostNegVal @I+    | r > 0 && r < minPosRat @I = leastPosVal @I+    | r < 0 && r > maxNegRat @I = leastNegVal @I+    | otherwise = buildIntRep @I r+ +  buildIntRep r =+    let (signBit,regime,exponent,significand) = posit2TupPosit @I r+        intRep = mkIntRep @I regime exponent significand+    in if signBit+       then negate intRep+       else intRep+ +  mkIntRep regime exponent significand =+    let (regime', offset) = formRegime @I regime  -- offset is the number of binary digits remaining after the regime is formed+        (exponent', offset') = formExponent @I exponent offset  -- offset' is the number of binary digits remaining after the exponent is formed+        fraction = formFraction @I significand offset'+    in regime' .|. exponent' .|. fraction+ +  formRegime power+    | 0 <= power =+      let offset = (fromIntegral (nBits @I - 1) -     power - 1)+      in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1)+    | otherwise =+      let offset = (fromIntegral (nBits @I - 1) - abs power - 1)+      in (1 `shiftL` fromInteger offset, offset)+ +  formExponent power offset =+    let offset' = offset - fromIntegral (exponentSize @I)+    in (fromIntegral power `shift` fromInteger offset', offset')+ +  formFraction r offset =+    let numFractionBits = offset+        fractionSize = 2^numFractionBits+        normFraction = round $ (r - 1) * fractionSize  -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1)+    in if numFractionBits >= 1+       then fromInteger normFraction+       else 0+ +  decode int+    | int == unReal @I = Nothing+    | int == 0 = Just 0+    | otherwise =+      let sgn = int < 0+          int' = if sgn+                 then negate int+                 else int+          (regime,nR) = regime2Integer @I int'+          exponent = exponent2Nat @I nR int'  -- if no e or some bits missing, then they are considered zero+          rat = fraction2Posit @I nR int'  -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1+      in tupPosit2Posit @I (sgn,regime,exponent,rat)+ +  regime2Integer posit =+    let regimeFormat = findRegimeFormat @I posit+        regimeCount = countRegimeBits @I regimeFormat posit+        regime = calcRegimeInt regimeFormat regimeCount+    in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime+ +  -- will return the format of the regime, either HI or LO; it could get refactored in the future+  -- True means a 1 is the first bit in the regime+  findRegimeFormat posit = testBit posit (fromIntegral (nBits @I) - 1 - fromIntegral (signBitSize @I))+ +  countRegimeBits format posit = go (fromIntegral (nBits @I) - 1 - fromIntegral (signBitSize @I)) 0+    where+      go (-1) acc = acc+      go index acc+        | xnor format (testBit posit index)  = go (index - 1) (acc + 1)+        | otherwise = acc+ +  -- knowing the number of the regime bits, and the sign bit we can extract+  -- the exponent.  We mask to the left of the exponent to remove the sign and regime, and+  -- then shift to the right to remove the fraction.+  exponent2Nat numBitsRegime posit =+    let bitsRemaining = fromIntegral (nBits @I) - numBitsRegime - fromIntegral (signBitSize @I)+        signNRegimeMask = 2^bitsRemaining - 1+        int = posit .&. signNRegimeMask+        nBitsToTheRight = fromIntegral (nBits @I) - numBitsRegime - fromIntegral (signBitSize @I) - fromIntegral (exponentSize @I)+    in if bitsRemaining <=0+       then 0+       else if nBitsToTheRight < 0+            then fromIntegral $ int `shiftL` negate nBitsToTheRight+            else fromIntegral $ int `shiftR` nBitsToTheRight+ +  -- knowing the number of the regime bits, sign bit, and the number of the+  -- exponent bits we can extract the fraction.  We mask to the left of the fraction to+  -- remove the sign, regime, and exponent. If there is no fraction then the value is 1.+  fraction2Posit numBitsRegime posit =+    let offset = fromIntegral $ (signBitSize @I) + fromIntegral numBitsRegime + (exponentSize @I)+        fractionSize = fromIntegral (nBits @I) - offset+        fractionBits = posit .&. (2^fractionSize - 1)+    in if fractionSize >= 1+       then (2^fractionSize + toInteger fractionBits) % 2^fractionSize+       else 1 % 1+ +  displayBin int = "0b" ++ go (fromIntegral (nBits @I) - 1)+    where+      go :: Int -> String+      go 0 = if testBit int 0+             then "1"+             else "0"+      go idx = if testBit int idx+               then "1" ++ go (idx - 1)+               else "0" ++ go (idx -1)++++instance PositC II where+  type IntN II = Int32+  exponentSize = 2+ +  -- Posit Integer Rep of various values+  unReal = minBound @Int32+ +  mostPosVal = maxBound @Int32+  leastPosVal = 1+  leastNegVal = -1+  mostNegVal = negate mostPosVal+ +  encode Nothing = unReal @II+  encode (Just 0) = 0+  encode (Just r)+    | r > maxPosRat @II = mostPosVal @II+    | r < minNegRat @II = mostNegVal @II+    | r > 0 && r < minPosRat @II = leastPosVal @II+    | r < 0 && r > maxNegRat @II = leastNegVal @II+    | otherwise = buildIntRep @II r+ +  buildIntRep r =+    let (signBit,regime,exponent,significand) = posit2TupPosit @II r+        intRep = mkIntRep @II regime exponent significand+    in if signBit+       then negate intRep+       else intRep+ +  mkIntRep regime exponent significand =+    let (regime', offset) = formRegime @II regime  -- offset is the number of binary digits remaining after the regime is formed+        (exponent', offset') = formExponent @II exponent offset  -- offset' is the number of binary digits remaining after the exponent is formed+        fraction = formFraction @II significand offset'+    in regime' .|. exponent' .|. fraction+ +  formRegime power+    | 0 <= power =+      let offset = (fromIntegral (nBits @II - 1) -     power - 1)+      in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1)+    | otherwise =+      let offset = (fromIntegral (nBits @II - 1) - abs power - 1)+      in (1 `shiftL` fromInteger offset, offset)+ +  formExponent power offset =+    let offset' = offset - fromIntegral (exponentSize @II)+    in (fromIntegral power `shift` fromInteger offset', offset')+ +  formFraction r offset =+    let numFractionBits = offset+        fractionSize = 2^numFractionBits+        normFraction = round $ (r - 1) * fractionSize  -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1)+    in if numFractionBits >= 1+       then fromInteger normFraction+       else 0+ +  decode int+    | int == unReal @II = Nothing+    | int == 0 = Just 0+    | otherwise =+      let sgn = int < 0+          int' = if sgn+                 then negate int+                 else int+          (regime,nR) = regime2Integer @II int'+          exponent = exponent2Nat @II nR int'  -- if no e or some bits missing, then they are considered zero+          rat = fraction2Posit @II nR int'  -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1+      in tupPosit2Posit @II (sgn,regime,exponent,rat)+ +  regime2Integer posit =+    let regimeFormat = findRegimeFormat @II posit+        regimeCount = countRegimeBits @II regimeFormat posit+        regime = calcRegimeInt regimeFormat regimeCount+    in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime+ +  -- will return the format of the regime, either HI or LO; it could get refactored in the future+  -- True means a 1 is the first bit in the regime+  findRegimeFormat posit = testBit posit (fromIntegral (nBits @II) - 1 - fromIntegral (signBitSize @II))+ +  countRegimeBits format posit = go (fromIntegral (nBits @II) - 1 - fromIntegral (signBitSize @II)) 0+    where+      go (-1) acc = acc+      go index acc+        | xnor format (testBit posit index)  = go (index - 1) (acc + 1)+        | otherwise = acc+ +  -- knowing the number of the regime bits, and the sign bit we can extract+  -- the exponent.  We mask to the left of the exponent to remove the sign and regime, and+  -- then shift to the right to remove the fraction.+  exponent2Nat numBitsRegime posit =+    let bitsRemaining = fromIntegral (nBits @II) - numBitsRegime - fromIntegral (signBitSize @II)+        signNRegimeMask = 2^bitsRemaining - 1+        int = posit .&. signNRegimeMask+        nBitsToTheRight = fromIntegral (nBits @II) - numBitsRegime - fromIntegral (signBitSize @II) - fromIntegral (exponentSize @II)+    in if bitsRemaining <=0+       then 0+       else if nBitsToTheRight < 0+            then fromIntegral $ int `shiftL` negate nBitsToTheRight+            else fromIntegral $ int `shiftR` nBitsToTheRight+ +  -- knowing the number of the regime bits, sign bit, and the number of the+  -- exponent bits we can extract the fraction.  We mask to the left of the fraction to+  -- remove the sign, regime, and exponent. If there is no fraction then the value is 1.+  fraction2Posit numBitsRegime posit =+    let offset = fromIntegral $ (signBitSize @II) + fromIntegral numBitsRegime + (exponentSize @II)+        fractionSize = fromIntegral (nBits @II) - offset+        fractionBits = posit .&. (2^fractionSize - 1)+    in if fractionSize >= 1+       then (2^fractionSize + toInteger fractionBits) % 2^fractionSize+       else 1 % 1+ +  displayBin int = "0b" ++ go (fromIntegral (nBits @II) - 1)+    where+      go :: Int -> String+      go 0 = if testBit int 0+             then "1"+             else "0"+      go idx = if testBit int idx+               then "1" ++ go (idx - 1)+               else "0" ++ go (idx -1)++++instance PositC III where+  type IntN III = Int64+  exponentSize = 3+ +  -- Posit Integer Rep of various values+  unReal = minBound @Int64+ +  mostPosVal = maxBound @Int64+  leastPosVal = 1+  leastNegVal = -1+  mostNegVal = negate mostPosVal+ +  encode Nothing = unReal @III+  encode (Just 0) = 0+  encode (Just r)+    | r > maxPosRat @III = mostPosVal @III+    | r < minNegRat @III = mostNegVal @III+    | r > 0 && r < minPosRat @III = leastPosVal @III+    | r < 0 && r > maxNegRat @III = leastNegVal @III+    | otherwise = buildIntRep @III r+ +  buildIntRep r =+    let (signBit,regime,exponent,significand) = posit2TupPosit @III r+        intRep = mkIntRep @III regime exponent significand+    in if signBit+       then negate intRep+       else intRep+ +  mkIntRep regime exponent significand =+    let (regime', offset) = formRegime @III regime  -- offset is the number of binary digits remaining after the regime is formed+        (exponent', offset') = formExponent @III exponent offset  -- offset' is the number of binary digits remaining after the exponent is formed+        fraction = formFraction @III significand offset'+    in regime' .|. exponent' .|. fraction+ +  formRegime power+    | 0 <= power =+      let offset = (fromIntegral (nBits @III - 1) -     power - 1)+      in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1)+    | otherwise =+      let offset = (fromIntegral (nBits @III - 1) - abs power - 1)+      in (1 `shiftL` fromInteger offset, offset)+ +  formExponent power offset =+    let offset' = offset - fromIntegral (exponentSize @III)+    in (fromIntegral power `shift` fromInteger offset', offset')+ +  formFraction r offset =+    let numFractionBits = offset+        fractionSize = 2^numFractionBits+        normFraction = round $ (r - 1) * fractionSize  -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1)+    in if numFractionBits >= 1+       then fromInteger normFraction+       else 0+ +  decode int+    | int == unReal @III = Nothing+    | int == 0 = Just 0+    | otherwise =+      let sgn = int < 0+          int' = if sgn+                 then negate int+                 else int+          (regime,nR) = regime2Integer @III int'+          exponent = exponent2Nat @III nR int'  -- if no e or some bits missing, then they are considered zero+          rat = fraction2Posit @III nR int'  -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1+      in tupPosit2Posit @III (sgn,regime,exponent,rat)+ +  regime2Integer posit =+    let regimeFormat = findRegimeFormat @III posit+        regimeCount = countRegimeBits @III regimeFormat posit+        regime = calcRegimeInt regimeFormat regimeCount+    in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime+ +  -- will return the format of the regime, either HI or LO; it could get refactored in the future+  -- True means a 1 is the first bit in the regime+  findRegimeFormat posit = testBit posit (fromIntegral (nBits @III) - 1 - fromIntegral (signBitSize @III))+ +  countRegimeBits format posit = go (fromIntegral (nBits @III) - 1 - fromIntegral (signBitSize @III)) 0+    where+      go (-1) acc = acc+      go index acc+        | xnor format (testBit posit index)  = go (index - 1) (acc + 1)+        | otherwise = acc+ +  -- knowing the number of the regime bits, and the sign bit we can extract+  -- the exponent.  We mask to the left of the exponent to remove the sign and regime, and+  -- then shift to the right to remove the fraction.+  exponent2Nat numBitsRegime posit =+    let bitsRemaining = fromIntegral (nBits @III) - numBitsRegime - fromIntegral (signBitSize @III)+        signNRegimeMask = 2^bitsRemaining - 1+        int = posit .&. signNRegimeMask+        nBitsToTheRight = fromIntegral (nBits @III) - numBitsRegime - fromIntegral (signBitSize @III) - fromIntegral (exponentSize @III)+    in if bitsRemaining <=0+       then 0+       else if nBitsToTheRight < 0+            then fromIntegral $ int `shiftL` negate nBitsToTheRight+            else fromIntegral $ int `shiftR` nBitsToTheRight+ +  -- knowing the number of the regime bits, sign bit, and the number of the+  -- exponent bits we can extract the fraction.  We mask to the left of the fraction to+  -- remove the sign, regime, and exponent. If there is no fraction then the value is 1.+  fraction2Posit numBitsRegime posit =+    let offset = fromIntegral $ (signBitSize @III) + fromIntegral numBitsRegime + (exponentSize @III)+        fractionSize = fromIntegral (nBits @III) - offset+        fractionBits = posit .&. (2^fractionSize - 1)+    in if fractionSize >= 1+       then (2^fractionSize + toInteger fractionBits) % 2^fractionSize+       else 1 % 1+ +  displayBin int = "0b" ++ go (fromIntegral (nBits @III) - 1)+    where+      go :: Int -> String+      go 0 = if testBit int 0+             then "1"+             else "0"+      go idx = if testBit int idx+               then "1" ++ go (idx - 1)+               else "0" ++ go (idx -1)++++instance PositC IV where+  type IntN IV = Int128+  exponentSize = 4+ +  -- Posit Integer Rep of various values+  unReal = minBound @Int128+ +  mostPosVal = maxBound @Int128+  leastPosVal = 1+  leastNegVal = -1+  mostNegVal = negate mostPosVal+ +  encode Nothing = unReal @IV+  encode (Just 0) = 0+  encode (Just r)+    | r > maxPosRat @IV = mostPosVal @IV+    | r < minNegRat @IV = mostNegVal @IV+    | r > 0 && r < minPosRat @IV = leastPosVal @IV+    | r < 0 && r > maxNegRat @IV = leastNegVal @IV+    | otherwise = buildIntRep @IV r+ +  buildIntRep r =+    let (signBit,regime,exponent,significand) = posit2TupPosit @IV r+        intRep = mkIntRep @IV regime exponent significand+    in if signBit+       then negate intRep+       else intRep+ +  mkIntRep regime exponent significand =+    let (regime', offset) = formRegime @IV regime  -- offset is the number of binary digits remaining after the regime is formed+        (exponent', offset') = formExponent @IV exponent offset  -- offset' is the number of binary digits remaining after the exponent is formed+        fraction = formFraction @IV significand offset'+    in regime' .|. exponent' .|. fraction+ +  formRegime power+    | 0 <= power =+      let offset = (fromIntegral (nBits @IV - 1) -     power - 1)+      in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1)+    | otherwise =+      let offset = (fromIntegral (nBits @IV - 1) - abs power - 1)+      in (1 `shiftL` fromInteger offset, offset)+ +  formExponent power offset =+    let offset' = offset - fromIntegral (exponentSize @IV)+    in (fromIntegral power `shift` fromInteger offset', offset')+ +  formFraction r offset =+    let numFractionBits = offset+        fractionSize = 2^numFractionBits+        normFraction = round $ (r - 1) * fractionSize  -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1)+    in if numFractionBits >= 1+       then fromInteger normFraction+       else 0+ +  decode int+    | int == unReal @IV = Nothing+    | int == 0 = Just 0+    | otherwise =+      let sgn = int < 0+          int' = if sgn+                 then negate int+                 else int+          (regime,nR) = regime2Integer @IV int'+          exponent = exponent2Nat @IV nR int'  -- if no e or some bits missing, then they are considered zero+          rat = fraction2Posit @IV nR int'  -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1+      in tupPosit2Posit @IV (sgn,regime,exponent,rat)+ +  regime2Integer posit =+    let regimeFormat = findRegimeFormat @IV posit+        regimeCount = countRegimeBits @IV regimeFormat posit+        regime = calcRegimeInt regimeFormat regimeCount+    in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime+ +  -- will return the format of the regime, either HI or LO; it could get refactored in the future+  -- True means a 1 is the first bit in the regime+  findRegimeFormat posit = testBit posit (fromIntegral (nBits @IV) - 1 - fromIntegral (signBitSize @IV))+ +  countRegimeBits format posit = go (fromIntegral (nBits @IV) - 1 - fromIntegral (signBitSize @IV)) 0+    where+      go (-1) acc = acc+      go index acc+        | xnor format (testBit posit index)  = go (index - 1) (acc + 1)+        | otherwise = acc+ +  -- knowing the number of the regime bits, and the sign bit we can extract+  -- the exponent.  We mask to the left of the exponent to remove the sign and regime, and+  -- then shift to the right to remove the fraction.+  exponent2Nat numBitsRegime posit =+    let bitsRemaining = fromIntegral (nBits @IV) - numBitsRegime - fromIntegral (signBitSize @IV)+        signNRegimeMask = 2^bitsRemaining - 1+        int = posit .&. signNRegimeMask+        nBitsToTheRight = fromIntegral (nBits @IV) - numBitsRegime - fromIntegral (signBitSize @IV) - fromIntegral (exponentSize @IV)+    in if bitsRemaining <=0+       then 0+       else if nBitsToTheRight < 0+            then fromIntegral $ int `shiftL` negate nBitsToTheRight+            else fromIntegral $ int `shiftR` nBitsToTheRight+ +  -- knowing the number of the regime bits, sign bit, and the number of the+  -- exponent bits we can extract the fraction.  We mask to the left of the fraction to+  -- remove the sign, regime, and exponent. If there is no fraction then the value is 1.+  fraction2Posit numBitsRegime posit =+    let offset = fromIntegral $ (signBitSize @IV) + fromIntegral numBitsRegime + (exponentSize @IV)+        fractionSize = fromIntegral (nBits @IV) - offset+        fractionBits = posit .&. (2^fractionSize - 1)+    in if fractionSize >= 1+       then (2^fractionSize + toInteger fractionBits) % 2^fractionSize+       else 1 % 1+ +  displayBin int = "0b" ++ go (fromIntegral (nBits @IV) - 1)+    where+      go :: Int -> String+      go 0 = if testBit int 0+             then "1"+             else "0"+      go idx = if testBit int idx+               then "1" ++ go (idx - 1)+               else "0" ++ go (idx -1)++++instance PositC V where+  type IntN V = Int256+  exponentSize = 5+ +  -- Posit Integer Rep of various values+  unReal = minBound @Int256+ +  mostPosVal = maxBound @Int256+  leastPosVal = 1+  leastNegVal = -1+  mostNegVal = negate mostPosVal+ +  encode Nothing = unReal @V+  encode (Just 0) = 0+  encode (Just r)+    | r > maxPosRat @V = mostPosVal @V+    | r < minNegRat @V = mostNegVal @V+    | r > 0 && r < minPosRat @V = leastPosVal @V+    | r < 0 && r > maxNegRat @V = leastNegVal @V+    | otherwise = buildIntRep @V r+ +  buildIntRep r =+    let (signBit,regime,exponent,significand) = posit2TupPosit @V r+        intRep = mkIntRep @V regime exponent significand+    in if signBit+       then negate intRep+       else intRep+ +  mkIntRep regime exponent significand =+    let (regime', offset) = formRegime @V regime  -- offset is the number of binary digits remaining after the regime is formed+        (exponent', offset') = formExponent @V exponent offset  -- offset' is the number of binary digits remaining after the exponent is formed+        fraction = formFraction @V significand offset'+    in regime' .|. exponent' .|. fraction+ +  formRegime power+    | 0 <= power =+      let offset = (fromIntegral (nBits @V - 1) -     power - 1)+      in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1)+    | otherwise =+      let offset = (fromIntegral (nBits @V - 1) - abs power - 1)+      in (1 `shiftL` fromInteger offset, offset)+ +  formExponent power offset =+    let offset' = offset - fromIntegral (exponentSize @V)+    in (fromIntegral power `shift` fromInteger offset', offset')+ +  formFraction r offset =+    let numFractionBits = offset+        fractionSize = 2^numFractionBits+        normFraction = round $ (r - 1) * fractionSize  -- "posit - 1" is changing it from the significand to the fraction: [1,2) -> [0,1)+    in if numFractionBits >= 1+       then fromInteger normFraction+       else 0+ +  decode int+    | int == unReal @V = Nothing+    | int == 0 = Just 0+    | otherwise =+      let sgn = int < 0+          int' = if sgn+                 then negate int+                 else int+          (regime,nR) = regime2Integer @V int'+          exponent = exponent2Nat @V nR int'  -- if no e or some bits missing, then they are considered zero+          rat = fraction2Posit @V nR int'  -- if no fraction or some bits missing, then the missing bits are zero, making the significand p=1+      in tupPosit2Posit @V (sgn,regime,exponent,rat)+ +  regime2Integer posit =+    let regimeFormat = findRegimeFormat @V posit+        regimeCount = countRegimeBits @V regimeFormat posit+        regime = calcRegimeInt regimeFormat regimeCount+    in (regime, regimeCount + 1) -- a rational representation of the regime and the regimeCount plus rBar which is the numBitsRegime+ +  -- will return the format of the regime, either HI or LO; it could get refactored in the future+  -- True means a 1 is the first bit in the regime+  findRegimeFormat posit = testBit posit (fromIntegral (nBits @V) - 1 - fromIntegral (signBitSize @V))+ +  countRegimeBits format posit = go (fromIntegral (nBits @V) - 1 - fromIntegral (signBitSize @V)) 0+    where+      go (-1) acc = acc+      go index acc+        | xnor format (testBit posit index)  = go (index - 1) (acc + 1)+        | otherwise = acc+ +  -- knowing the number of the regime bits, and the sign bit we can extract+  -- the exponent.  We mask to the left of the exponent to remove the sign and regime, and+  -- then shift to the right to remove the fraction.+  exponent2Nat numBitsRegime posit =+    let bitsRemaining = fromIntegral (nBits @V) - numBitsRegime - fromIntegral (signBitSize @V)+        signNRegimeMask = 2^bitsRemaining - 1+        int = posit .&. signNRegimeMask+        nBitsToTheRight = fromIntegral (nBits @V) - numBitsRegime - fromIntegral (signBitSize @V) - fromIntegral (exponentSize @V)+    in if bitsRemaining <=0+       then 0+       else if nBitsToTheRight < 0+            then fromIntegral $ int `shiftL` negate nBitsToTheRight+            else fromIntegral $ int `shiftR` nBitsToTheRight+ +  -- knowing the number of the regime bits, sign bit, and the number of the+  -- exponent bits we can extract the fraction.  We mask to the left of the fraction to+  -- remove the sign, regime, and exponent. If there is no fraction then the value is 1.+  fraction2Posit numBitsRegime posit =+    let offset = fromIntegral $ (signBitSize @V) + fromIntegral numBitsRegime + (exponentSize @V)+        fractionSize = fromIntegral (nBits @V) - offset+        fractionBits = posit .&. (2^fractionSize - 1)+    in if fractionSize >= 1+       then (2^fractionSize + toInteger fractionBits) % 2^fractionSize+       else 1 % 1+ +  displayBin int = "0b" ++ go (fromIntegral (nBits @V) - 1)+    where+      go :: Int -> String+      go 0 = if testBit int 0+             then "1"+             else "0"+      go idx = if testBit int idx+               then "1" ++ go (idx - 1)+               else "0" ++ go (idx -1)+++-- =====================================================================+-- ===                Encode and Decode Helpers                      ===+-- =====================================================================+++-- getSign finds the sign value and then returns the absolute value of the Posit+getSign :: Rational -> (Bool, Rational)+getSign r =+  let s = r <= 0+      absPosit =+        if s+        then negate r+        else r+  in (s,absPosit)  -- pretty much the same as 'abs')++-- Exponent should be an integer in the range of [0,uSeed), and also return the posit [1,2)+getExponent :: Rational -> (Natural, Rational)+getExponent r = log_2 (0,r)++log_2 :: (Natural, Rational) -> (Natural, Rational)+log_2 (exponent,r) | r <  1 = error "Should never happen, exponent should be a natural number, i.e. positive integer."+                   | r >= (2 % 1) = log_2 (exponent+1,r * (1 % 2))+                   | otherwise = (exponent,r)+++calcRegimeInt :: Bool -> Int -> Integer+calcRegimeInt format count | format = fromIntegral (count - 1)+                           | otherwise = fromIntegral $ negate count+++xnor :: Bool -> Bool -> Bool+xnor a b = not $ (a || b) && not (b && a)+++#ifndef O_NO_ORPHANS+#ifndef O_NO_STORABLE+-- =====================================================================+-- ===                  Storable Instances                           ===+-- =====================================================================+--+-- Orphan Instance for Word128 using the DoubleWord type class+instance Storable Word128 where+  sizeOf _ = 16+  alignment _ = 16+  peek ptr = do+    hi <- peek $ offsetInt 0+    lo <- peek $ offsetWord 1+    return $ fromHiAndLo hi lo+      where+        offsetInt i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8)+        offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8)+  poke ptr int = do+    poke (offsetInt 0) (hiWord int)+    poke (offsetWord 1) (loWord int)+      where+        offsetInt i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8)+        offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8)++-- Orphan Instance for Int128 using the DoubleWord type class+instance Storable Int128 where+  sizeOf _ = 16+  alignment _ = 16+  peek ptr = do+    hi <- peek $ offsetInt 0+    lo <- peek $ offsetWord 1+    return $ fromHiAndLo hi lo+      where+        offsetInt i = (castPtr ptr :: Ptr Int64) `plusPtr` (i*8)+        offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8)+  poke ptr int = do+    poke (offsetInt 0) (hiWord int)+    poke (offsetWord 1) (loWord int)+      where+        offsetInt i = (castPtr ptr :: Ptr Int64) `plusPtr` (i*8)+        offsetWord i = (castPtr ptr :: Ptr Word64) `plusPtr` (i*8)++-- Orphan Instance for Int256 using the DoubleWord type class+instance Storable Int256 where+  sizeOf _ = 32+  alignment _ = 32+  peek ptr = do+    hi <- peek $ offsetInt 0+    lo <- peek $ offsetWord 1+    return $ fromHiAndLo hi lo+      where+        offsetInt i = (castPtr ptr :: Ptr Int128) `plusPtr` (i*16)+        offsetWord i = (castPtr ptr :: Ptr Word128) `plusPtr` (i*16)+  poke ptr int = do+    poke (offsetInt 0) (hiWord int)+    poke (offsetWord 1) (loWord int)+      where+        offsetInt i = (castPtr ptr :: Ptr Int128) `plusPtr` (i*16)+        offsetWord i = (castPtr ptr :: Ptr Word128) `plusPtr` (i*16)+--+#endif+#endif

@@ -0,0 +1,183 @@++--------------------------------------------------------------------------------------------+-- | Posit Numbers+--   Copyright   :  (C) 2022 Nathan Waivio+--   License     :  BSD3+--   Maintainer  :  Nathan Waivio <nathan.waivio@gmail.com>+--   Stability   :  Stable+--   Portability :  Portable+--+--   Test Suite for a Library implementing standard Posit Numbers+-- +---------------------------------------------------------------------------------------------++import Posit+import Posit.Internal.PositC++++main :: IO ()+main = do+--+  print $ "Machine Alpha Posit8 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit8)) -- succ (Posit int) = Posit (succ int)+  print $ "Machine Alpha Posit16 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit16)) -- +  print $ "Machine Alpha Posit32 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit32)) -- +  print $ "Machine Alpha Posit64 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit64)) -- +  print $ "Machine Alpha Posit128 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit128)) -- +  print $ "Machine Alpha Posit256 ~1.0: " ++ show (1.0 - succ (1.0 :: Posit256)) -- +  print $ "Does (1 - 1) == 0 ?: " ++ show ((1 - 1) == (0 :: Posit256)) -- [(1 - 1) == zero | zero = 0 :: Posit es, es <- Z .. V]+  let sqrtTaylor = (funLogDomainReduction funLogTaylor).(/2).(funExp2 funExpTaylor).(/log 2)+  print $ "sqrt phi using a Taylor algorithm: " ++ show (sqrtTaylor phi)+  let sqrtTuma = (funLogDomainReduction funLogTuma).(/2).(funExp2 funExpTuma).(/log 2)+  print $ "sqrt phi using a Tuma algorithm: " ++ show (sqrtTuma phi)+  print $ "Tuma is fasta: " ++ show (sqrtTaylor (1/1000000) - sqrtTuma (1/1000000))+  let truth = 5.0431656433600286513118821892854247103235901754138463603020001967777869609108929428415187821843384653305404495551887666992776792 :: Posit256+  eval "Standard: exp(phi):" (exp (phi)) truth+  eval "Taylor: exp(phi):" (funExp2 funExpTaylor (phi / log 2)) truth+  eval "Tuma: exp(phi):" (funExp2 funExpTuma (phi / log 2)) truth+  eval "Wolfram Alpha: exp(phi):" truth truth+  let truth = 2.6881171418161354484126255515800135873611118773741922415191608615280287034909564914158871097219845710811670879190576068697e43 :: Posit256+  eval "Standard: exp(100):" (exp (100)) truth+  eval "Taylor: exp(100):" (funExp2 funExpTaylor (100 / log 2)) truth+  eval "Tuma: exp(100):" (funExp2 funExpTuma (100 / log 2))  truth+  eval "Wolfram Alpha: exp(100):" truth truth+  let truth = 3.7200759760208359629596958038631183373588922923767819671206138766632904758958157181571187786422814966019356176423110698002e-44 :: Posit256+  eval "Standard: exp(-100):" (exp (-100)) truth+  eval "Taylor: exp(-100):" (funExp2 funExpTaylor (-100 / log 2)) truth+  eval "Tuma: exp(-100):" (funExp2 funExpTuma (-100 / log 2)) truth+  eval "Wolfram Alpha: exp(-100):" truth truth+  let truth = 1.9700711140170469938888793522433231253169379853238457899528029913850638507824411934749780765630268899309638179875202269359e434 :: Posit256+  eval "Standard: exp(1000):" (exp (1000)) truth+  eval "Taylor: exp(1000):" (funExp2 funExpTaylor (1000 / log 2)) truth+  eval "Tuma: exp(1000):" (funExp2 funExpTuma (1000 / log 2)) truth+  eval "Wolfram Alpha: exp(1000):" truth truth+  let truth = 5.075958897549456765291809479574336919305599282892837361832393845410540542974819175679662169046542867863667106831065285113e-435 :: Posit256+  eval "Standard: exp(-1000):" (exp (-1000)) truth+  eval "Taylor: exp(-1000):" (funExp2 funExpTaylor (-1000 / log 2)) truth+  eval "Tuma: exp(-1000):" (funExp2 funExpTuma (-1000 / log 2)) truth+  eval "Wolfram Alpha: exp(-1000):" truth truth+  let truth = 0.4812118250596034474977589134243684231351843343856605196610181688401638676082217744120094291227234749972318399582936564112725683 :: Posit256+  eval "Standard: log(phi):" (log (phi)) truth+  eval "Taylor: log(phi):" (funLogDomainReduction funLogTaylor (phi)) truth+  eval "Tuma: log(phi):" (funLogDomainReduction funLogTuma (phi)) truth+  eval "Wolfram Alpha: log(phi):" truth truth+  let truth = -4.6051701859880913680359829093687284152022029772575459520666558019351452193547049604719944101791965966839355680845724972668190 :: Posit256+  eval "Standard: log(1/100):" (log (1/100)) truth+  eval "Taylor: log(1/100):" (funLogDomainReduction funLogTaylor (1/100)) truth+  eval "Tuma: log(1/100):" (funLogDomainReduction funLogTuma (1/100)) truth+  eval "Wolfram Alpha: log(1/100):" truth truth+  let truth = -6.9077552789821370520539743640530926228033044658863189280999837029027178290320574407079916152687948950259033521268587459002285 :: Posit256+  eval "Standard: log(1/1000):" (log (1/1000)) truth+  eval "Taylor: log(1/1000):" (funLogDomainReduction funLogTaylor (1/1000)) truth+  eval "Tuma: log(1/1000):" (funLogDomainReduction funLogTuma (1/1000)) truth+  eval "Wolfram Alpha: log(1/1000):" truth truth+  let truth = 4.5347571611551792889915884948567915637887680293971326427244942079650289300980475282698882636812383679690567084677326507550787791 :: Posit256+  eval "Standard: phi^pi:" ((phi) ** pi) truth+  eval "Wolfram Alpha: phi^pi:" truth truth+  let truth = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446 :: Posit256+  eval "Standard pi:" pi truth+  eval "Gauss–Legendre algorithm: pi:" funPi1 truth+  eval "Borwein's algorithm: pi:" funPi2 truth+  eval "Bailey–Borwein–Plouffe (BBP) formula: pi:" funPi3 truth+  eval "Fabrice Bellard improvement on the BBP: pi:" funPi4 truth+  eval "Wolfram Alpha: pi:" truth truth+  eval "Bailey–Borwein–Plouffe (BBP) formula: but succ pi:" (succ funPi3) truth+--+  -- print $ "Does (1 - 1) == 0 ?: " ++ (1 - 1) == (0 :: Posit256) -- [(1 - 1) == zero | zero = 0 :: Posit es, es <- Z .. V]+  print "Now for Property testing of Posit8... (This should generalize for all other Posit types)"+  print $ "Does associtivity of (+) hold?: " ++ (show assoc8)+  print $ "Does commutitivity of (+) hold?: " ++ (show commutative8)+  print $ "Is `fromInteger 0` the additive identity?: " ++ (show additiveIdent8)+  print $ "Does `negate` give the additive inverse? (excluding NaR): " ++ (show additiveInv8)+  print $ "Does `negate.negate == id`?: " ++ (show nn8)+  print $ "Does associtivity of (*) hold?: " ++ (show assocMult8)+  print $ "Is `fromInteger 1` the multiplicitive identity?: " ++ (show multIdent8)+  print $ "Does Reflexivity of Eq hold?: " ++ (show reflEq8)+  print $ "Does Symmetry of Eq hold?: " ++ (show symEq8)+  print $ "Does Transitivity of Eq hold?: " ++ (show transEq8)+  print $ "Does Extensionality of Eq hold?: " ++ (show extEq8)+  print $ "Does Negation of Eq hold?: " ++ (show negEq8)+  print $ "Does Comparability of Ord hold?: " ++ (show comp8)+  print $ "Does Transitivity of Ord hold?: " ++ (show trans8)+  print $ "Does Reflexivity of Ord hold?: " ++ (show refl8)+  print $ "Does Antisymmetry of Ord hold?: " ++ (show anti8)+  print $ "Does the `abs x * signum x == x` law hold?: " ++ (show absSignumLaw)+  print $ "Is recip a multiplicative inverse?: " ++ (show recipInv8)+  print $ "Are there any `recip.recip == id` values: " ++ (show rr8)+  print $ "Are there any `recip.recip /= id` values: " ++ (show rrne8)+  print $ "Does the distributive property hold with posits all the time?: " ++ (show doesItDistribute)+  print $ "Exaustive Proof... for fused ops recovering the distributeive property... and it turns out to be true."+  print $ "Can fused ops recover the distributive property for `fmms a b (negate a) c == fam b c a` ?: " ++ (show fusedDistribute)++++eval :: String -> Posit256 -> Posit256 -> IO ()+eval msg val tru = putStr $ msg ++ "\n" ++ (show val) ++ "\n" ++ "ULP: " ++ (show $ valInt - truInt) ++ "\n"+  where+    valInt = read (displayIntegral val) :: Integer+    truInt = read (displayIntegral tru) :: Integer++-- exaustive testing, enum from to+assoc8 :: Bool+assoc8 = and [(x + y) + z == x + (y + z) | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), z <- enumFrom (NaR :: Posit8)]++commutative8 :: Bool+commutative8 = and [x + y == y + x | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8)]++additiveIdent8 :: Bool+additiveIdent8 = and [x + fromInteger 0 == x | x <- enumFrom (NaR :: Posit8)]++additiveInv8 :: Bool+additiveInv8 = and [x + negate x == fromInteger 0 | x <- enumFrom (minBound :: Posit8)]++assocMult8 :: Bool+assocMult8 = and [(x * y) * z == x * (y * z) | x <- enumFrom (minBound :: Posit8), y <- enumFrom (minBound :: Posit8), z <- enumFrom (minBound :: Posit8)]++multIdent8 :: Bool+multIdent8 = and [x * fromInteger 1 == x && fromInteger 1 * x == x | x <- enumFrom (NaR :: Posit8)]++reflEq8 = and [(x == x) | x <- enumFrom (NaR :: Posit8)]++symEq8 = and [y == x | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), (x == y)]++transEq8 = and [x == z | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), z <- enumFrom (NaR :: Posit8), (x == y) && (y == z)]++extEq8 = and [sin x == sin y | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), x == y]++negEq8 = and [not (x == y) | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8), x /= y]++comp8 :: Bool+comp8 = and [(x <= y || y <= x) == True | x <- enumFrom (NaR :: Posit8), y <- enumFrom (NaR :: Posit8)]++trans8 :: Bool+trans8 = and [(x <= z) == True | x <- enumFrom (minBound :: Posit8), y <- enumFrom (minBound :: Posit8), z <- enumFrom (minBound :: Posit8), (x <= y && y <= z) == True]++refl8 :: Bool+refl8 = and [(x <= x) == True | x <- enumFrom (NaR :: Posit8)]++anti8 :: Bool+anti8 = and [(x == y) == True | x <- enumFrom (minBound :: Posit8), y <- enumFrom (minBound :: Posit8), (x <= y && y <= x) == True]++nn8 :: Bool+nn8 = and [(negate.negate $ x) == x | x <- enumFrom (NaR :: Posit8)]++-- recip.recip == id +rr8 :: [Posit8]+rr8 = [x| x <- enumFrom (NaR :: Posit8), (recip.recip $ x) == x]++-- recip.recip /= id+rrne8 :: [Posit8]+rrne8 = [x| x <- enumFrom (NaR :: Posit8), (recip.recip $ x) /= x]++doesItDistribute :: Bool+doesItDistribute = and [a*b + a*c == a*(b+c) | a <- enumFrom (NaR :: Posit8), b <- enumFrom (NaR :: Posit8), c <- enumFrom (NaR :: Posit8)]++fusedDistribute :: Bool+fusedDistribute = and [fmms a b (negate a) c == fam b c a | a <- enumFrom (NaR :: Posit8), b <- enumFrom (NaR :: Posit8), c <- enumFrom (NaR :: Posit8)]++absSignumLaw :: Bool+absSignumLaw = and [abs x * signum x == x | x <- enumFrom (NaR :: Posit8)]++recipInv8 :: Bool+recipInv8 = and [((x * recip x) == fromInteger 1) && ((recip x * x) == fromInteger 1)  | x <- enumFrom (NaR :: Posit8)]+