diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,8 +1,19 @@ # Changelog for Posit Numbers +# posit-2022.1.0.0 + + * Optimized the `Show` instance to properly handle the tapered precision + * Changes to round better + * Improved accuracy of `exp` + * Changed titles of accuracy charts to be more consistant + * Added test to confirm `read.show == id` to perhaps optimise the textual representation, and verify that the implementation is sufficent + * Test round trip with command to run: `stack test posit:test-posit-readShowId` + * Added catagory to .cabal file + # posit-2022.0.1.4 * `atan 1 :: Posit256` is very slow, added a few precice constants to speed it up + * Improved accuracy of `sin` and `cos` at the extreme positive and negative values # posit-2022.0.1.3 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# posit 2022.0.1.4 +# posit 2022.1.0.0 The [Posit Standard 2022](https://posithub.org/docs/posit_standard-2.pdf), and [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf), diff --git a/posit.cabal b/posit.cabal --- a/posit.cabal +++ b/posit.cabal @@ -1,7 +1,8 @@ cabal-version: 1.12 name: posit -version: 2022.0.1.4 +version: 2022.1.0.0 +category: Numeric, Math description: The Posit Number format attempting to conform to the Posit Standard Versions 3.2 and 2022. Where Real numbers are approximated by `Maybe Rational` and sampled in a similar way to the projective real line. homepage: https://github.com/waivio/posit#readme bug-reports: https://github.com/waivio/posit/issues @@ -15,7 +16,7 @@ GHC == 8.10.7, GHC == 9.0.2, GHC == 9.2.8, - GHC == 9.4.5, + GHC == 9.4.6, GHC == 9.6.2 synopsis: Posit Numbers extra-source-files: @@ -113,8 +114,19 @@ main-is: WeighPosit.hs ghc-options: -Wall -O2 build-depends: - base >=4.7 && <5, + base >= 4.7 && < 5, posit, vector, weigh + default-language: Haskell2010 + +-- Test `(show.read == id)` +benchmark test-posit-readShowId + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: TestReadShow.hs + ghc-options: -Wall -O2 + build-depends: + base >= 4.7 && < 5, + posit default-language: Haskell2010 diff --git a/src/Posit.hs b/src/Posit.hs --- a/src/Posit.hs +++ b/src/Posit.hs @@ -201,7 +201,7 @@ -- instance PositC es => Show (Posit es) where show NaR = "NaR" - show (R r) = formatScientific Generic (Just $ decimalPrec @es) (fst.fromRationalRepetendUnlimited $ r) + show p@(Posit int) = formatScientific Generic (Just $ decimalPrec @es int) (fst.fromRationalRepetendUnlimited $ toRational p) -- #endif @@ -710,11 +710,14 @@ approx_exp :: forall es. PositC es => Posit es -> Posit es -- Comment by Abigale Emily: xcddfffff -approx_exp (R bx) = fromIntegral (uSeed @es)^^m * (taylor_approx_expm1 (R c * log_USeed) + 1) +approx_exp NaR = NaR +approx_exp (R x) = 2^^k * funExpTaylor (R r) where - (m,c) = properFraction $ bx / luS - R luS = log_USeed @es + k = floor (x / ln2 + 0.5) -- should be Integer or Int + r = x - fromInteger k * ln2 -- should be Rational +ln2 :: Rational +ln2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202 log_USeed :: forall es. PositC es => Posit es log_USeed = approx_log $ fromIntegral (uSeed @es) @@ -870,7 +873,6 @@ | otherwise = tuma_approx_cos $ 2*approx_pi * x -- - -- Using the CORDIC domain reduction and some approximation function of log funLogDomainReduction :: forall es. PositC es => (Posit es -> Posit es) -> Posit es -> Posit es funLogDomainReduction _ NaR = NaR @@ -881,9 +883,9 @@ where (ex, sig) = (int * fromIntegral (2^(exponentSize @es)) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1 (_,int,nat,rat) = (posit2TupPosit @es).toRational $ x -- sign should always be positive - - + + -- natural log with log phi acurate to 9 ULP funLogTaylor :: forall es. PositC es => Posit es -> Posit es funLogTaylor NaR = NaR @@ -925,9 +927,25 @@ taylor_approx_expm1 x = go 0 [x^n / fromIntegral (fac n) | n <- [1..]] where go :: Posit es -> [Posit es] -> Posit es + go !acc [] = acc go !acc (h:t) | acc == acc + h = acc | otherwise = go (acc + h) t -- need to use a Taylor Series, the `tanh` formulation doesn't work because it requires something that depends on `exp` + +-- +-- calculate exp, its most accurate near zero +-- sum k=0 to k=inf of the terms, iterate until a fixed point is reached +funExpTaylor :: forall es. PositC es => Posit es -> Posit es +funExpTaylor NaR = NaR +funExpTaylor 0 = 1 +funExpTaylor z = go 0 0 + where + go :: Natural -> Posit es -> Posit es + go !k !acc + | acc == (acc + term k) = acc -- if x == x + dx then terminate and return x + | otherwise = go (k+1) (acc + term k) + term :: Natural -> Posit es + term k = (z^k) / (fromIntegral.fac $ k) -- ===================================================================== diff --git a/src/Posit/Internal/PositC.hs b/src/Posit/Internal/PositC.hs --- a/src/Posit/Internal/PositC.hs +++ b/src/Posit/Internal/PositC.hs @@ -25,6 +25,7 @@ {-# LANGUAGE DerivingVia #-} -- To Derive instances for newtypes to eliminate Orphan Instances {-# LANGUAGE UndecidableInstances #-} -- For deriving DoubleWord {-# LANGUAGE CPP #-} -- To remove Storable instances to remove noise when performing analysis of Core +{-# LANGUAGE MultiWayIf #-} {-# OPTIONS_GHC -Wno-unticked-promoted-constructors #-} -- Turn off noise {-# OPTIONS_GHC -Wno-type-defaults #-} -- Turn off noise @@ -32,7 +33,6 @@ -- |Posit Class, implementing: -- -- * PositC --- * Orphan Instances of Storable for Word128, Int128, Int256 -- ---- module Posit.Internal.PositC @@ -61,7 +61,7 @@ import Data.Int (Int8,Int16,Int32,Int64) -- Import standard Int sizes import Data.DoubleWord (Word128,Int128,Int256,fromHiAndLo,hiWord,loWord,DoubleWord,BinaryWord) -- Import large Int sizes import Data.Word (Word64) -import Data.Bits (Bits(..), shiftL, shift, testBit, (.&.), shiftR,FiniteBits) +import Data.Bits (Bits(..), shiftL, shiftR, testBit, (.&.), FiniteBits) -- Import Naturals and Rationals @@ -291,7 +291,7 @@ formRegime :: Integer -> (IntN es, Integer) formRegime power | 0 <= power = - let offset = (fromIntegral (nBits @es - 1) - power - 1) + let offset = (fromIntegral (nBits @es - 1) - power - 1) in (fromIntegral (2^(power + 1) - 1) `shiftL` fromInteger offset, offset - 1) | otherwise = let offset = (fromIntegral (nBits @es - 1) - abs power - 1) @@ -300,16 +300,23 @@ formExponent :: Natural -> Integer -> (IntN es, Integer) formExponent power offset = let offset' = offset - fromIntegral (exponentSize @es) - in (fromIntegral power `shift` fromInteger offset', offset') + result = fromIntegral power `shiftR` negate (fromInteger offset') + in if offset' >= 0 + then (fromIntegral power `shiftL` fromInteger offset', offset') + else if testBit (fromIntegral power :: Int) (pred.fromIntegral.negate $ offset') -- What is "Twilight Zone" Posit Number? + then (succ result, offset') + else (result, offset') -- rounding case formFraction :: Rational -> Integer -> IntN es 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 + in if | numFractionBits >= 1 -> fromInteger normFraction + | numFractionBits == 0 -> if 2 * (r - 1) > 1 -- rounding case -- Previously was 0, this should be "banker's rounding" per Gustafson's ch3 + then 1 + else 0 + | otherwise -> 0 tupPosit2Posit :: (Bool,Integer,Natural,Rational) -> Maybe Rational tupPosit2Posit (sgn,regime,exponent,rat) = -- s = isNeg posit == True @@ -379,9 +386,17 @@ else "0" ++ go (idx - 1) -- decimal Precision - decimalPrec :: Int - decimalPrec = fromIntegral $ 2 * (nBytes @es) + 1 - + decimalPrec :: IntN es -> Int + -- decimalPrec _ = fromIntegral $ 2 * (nBytes @es) + 1 -- The Olden way, 1.37% fail to round trip + -- + -- `1 + ` Gustafson's Decimal Accuracy (from coorispondance with Alessandro) -- Sucessfuly round triped for both Posit32 and P32, and below + decimalPrec posit = + let regimeFormat = findRegimeFormat @es posit + regimeCount = countRegimeBits @es regimeFormat posit + fractionSize = max (fromIntegral (nBits @es) - fromIntegral (signBitSize @es) - regimeCount - fromIntegral (exponentSize @es)) 0 -- fractionSize is at least zero + in succ.ceiling $ (fromIntegral fractionSize + 1) * log10Of2 + halflog10Of2 + -- + {-# MINIMAL exponentSize | nBytes #-} @@ -462,6 +477,14 @@ xnor :: Bool -> Bool -> Bool xnor a b = not $ (a || b) && not (b && a) + + +log10Of2 :: Double +log10Of2 = 0.3010299956639811952137388947 + + +halflog10Of2 :: Double +halflog10Of2 = 0.1505149978319905976068694473 #ifndef O_NO_STORABLE_RANDOM diff --git a/stack.yaml b/stack.yaml --- a/stack.yaml +++ b/stack.yaml @@ -2,7 +2,7 @@ # that coorispond to a specific GHC or Stackage version # resolver: nightly-2023-07-06 # ghc-9.6.2 -resolver: lts-21.1 # ghc-9.4.5 +resolver: lts-21.9 # ghc-9.4.6 # resolver: lts-20.26 # ghc-9.2.8 # resolver: lts-19.33 # ghc-9.0.2 # resolver: lts-18.28 # ghc-8.10.7 diff --git a/test/Test/Algorithms.hs b/test/Test/Algorithms.hs --- a/test/Test/Algorithms.hs +++ b/test/Test/Algorithms.hs @@ -126,11 +126,13 @@ z = 1 / 2^n - - - - - +{- +approx_exp :: forall es. PositC es => Posit es -> Posit es +approx_exp (R bx) = fromIntegral (uSeed @es)^^m * (taylor_approx_expm1 (R c * log_USeed) + 1) + where + (m,c) = properFraction $ bx / luS + R luS = log_USeed @es +-} @@ -301,8 +303,8 @@ go !prevA !prevS !n !a !s | prevA == a = a | prevS == s = a - | abs (prevA - a) <= 2*eps = a -- P256 or Posit128, will not reach a fixed point where `prevA == a` it sort of oscelates until divergence occurs, if we test for less than 2*eps it can stop early - | abs (prevS - s) <= 2*eps = a + | abs (prevA - a) <= 3*eps = a -- P256 or Posit128, will not reach a fixed point where `prevA == a` it sort of oscelates until divergence occurs, if we test for less than 3*eps it can stop early + | abs (prevS - s) <= 3*eps = a | otherwise = let x = 5 / s - 1 y = (x - 1)^2 + 7 diff --git a/test/TestElementaryFunctions.hs b/test/TestElementaryFunctions.hs --- a/test/TestElementaryFunctions.hs +++ b/test/TestElementaryFunctions.hs @@ -32,91 +32,91 @@ expPlotP16Posit16 = toFile def "./test/Results/Bits Accuracy of exp with P16 and Posit16.png" $ do let expP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (exp p) exp) | p <- enumFrom (NaR :: P16)] expPosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (exp p) exp) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 exp" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'exp'" plot (line "exp P16 error" [expP16]) plot (line "exp Posit16 error" [expPosit16]) logPlotP16Posit16 = toFile def "./test/Results/Bits Accuracy of log with P16 and Posit16.png" $ do let lnP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (log p) log) | p <- enumFrom (0 :: P16)] lnPosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (log p) log) | p <- enumFrom (0 :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Log" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'log'" plot (line "log P16 error" [lnP16]) plot (line "log Posit16 error" [lnPosit16]) sqrtPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of sqrt with P16 and Posit16.png" $ do let sqrtP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (sqrt p) sqrt) | p <- enumFrom (0 :: P16)] sqrtPosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (sqrt p) sqrt) | p <- enumFrom (0 :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Square Root" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'sqrt'" plot (line "sqrt P16 error" [sqrtP16]) plot (line "sqrt Posit16 error" [sqrtPosit16]) sinPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of sin with P16 and Posit16.png" $ do let sineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (sin p) sin) | p <- enumFrom (NaR :: P16)] sinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (sin p) sin) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Sine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'sin'" plot (line "sin P16 error" [sineP16]) plot (line "sin Posit16 error" [sinePosit16]) cosPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of cos with P16 and Posit16.png" $ do let cosineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (cos p) cos) | p <- enumFrom (NaR :: P16)] cosinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (cos p) cos) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Cosine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'cos'" plot (line "cos P16 error" [cosineP16]) plot (line "cos Posit16 error" [cosinePosit16]) asinPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of asin with P16 and Posit16.png" $ do let arcsineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (asin p) sin) | p <- enumFromTo (-1 :: P16) 1] arcsinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (asin p) sin) | p <- enumFromTo (-1 :: Posit16) 1] - layout_title .= "Number of Accurate Bits of P16 and Posit16 arcSine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'asin'" plot (line "asin P16 error" [arcsineP16]) plot (line "asin Posit16 error" [arcsinePosit16]) acosPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of acos with P16 and Posit16.png" $ do let arccosineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (acos p) acos) | p <- enumFromTo (-1 :: P16) 1] arccosinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (acos p) acos) | p <- enumFromTo (-1 :: Posit16) 1] - layout_title .= "Number of Accurate Bits of P16 and Posit16 arcCosine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'acos'" plot (line "acos P16 error" [arccosineP16]) plot (line "acos Posit16 error" [arccosinePosit16]) atanPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of atan with P16 and Posit16.png" $ do let arctangentP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (atan p) atan) | p <- enumFrom (NaR :: P16)] arctangentPosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (atan p) atan) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 arcTangent" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'atan'" plot (line "atan P16 error" [arctangentP16]) plot (line "atan Posit16 error" [arctangentPosit16]) sinhPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of sinh with P16 and Posit16.png" $ do let hypsineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (sinh p) sinh) | p <- enumFrom (NaR :: P16)] hypsinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (sinh p) sinh) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Hyperbolic Sine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'sinh'" plot (line "sinh P16 error" [hypsineP16]) plot (line "sinh Posit16 error" [hypsinePosit16]) coshPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of cosh with P16 and Posit16.png" $ do let hypcosineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (cosh p) cosh) | p <- enumFrom (NaR :: P16)] hypcosinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (cosh p) cosh) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Hyperbolic Cosine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'cosh'" plot (line "cosh P16 error" [hypcosineP16]) plot (line "cosh Posit16 error" [hypcosinePosit16]) asinhPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of asinh with P16 and Posit16.png" $ do let archypsineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (asinh p) asinh) | p <- enumFrom (NaR :: P16)] archypsinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (asinh p) asinh) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Hyperbolic Sine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'asinh'" plot (line "asinh P16 error" [archypsineP16]) plot (line "asinh Posit16 error" [archypsinePosit16]) acoshPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of acosh with P16 and Posit16.png" $ do let archypcosineP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (acosh p) acosh) | p <- enumFrom (NaR :: P16)] archypcosinePosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (acosh p) acosh) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Inv Hyperbolic Cosine" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'acosh'" plot (line "acosh P16 error" [archypcosineP16]) plot (line "acosh Posit16 error" [archypcosinePosit16]) atanhPlotsP16Posit16 = toFile def "./test/Results/Bits Accuracy of atanh with P16 and Posit16.png" $ do let archyptangentP16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (atanh p) atanh) | p <- enumFrom (NaR :: P16)] archyptangentPosit16 = filter (\(_,d) -> not $ nanOrInf d) [(read (displayIntegral p) :: Double, err p (atanh p) atanh) | p <- enumFrom (NaR :: Posit16)] - layout_title .= "Number of Accurate Bits of P16 and Posit16 Inv Hyperbolic Tangent" + layout_title .= "Number of Accurate Bits of P16 and Posit16 'atanh'" plot (line "atanh P16 error" [archyptangentP16]) plot (line "atanh Posit16 error" [archyptangentPosit16]) diff --git a/test/TestPosit.hs b/test/TestPosit.hs --- a/test/TestPosit.hs +++ b/test/TestPosit.hs @@ -75,13 +75,21 @@ print $ "Machine epsilon P64 ~1.0: " ++ show (eps :: P64) -- print $ "Machine epsilon P128 ~1.0: " ++ show (eps :: P128) -- print $ "Machine epsilon P256 ~1.0: " ++ show (eps :: P256) -- - -- | Taylor vs. Tuma print $ "Does (1 - 1) == 0 ?: " ++ show ((1 - 1) == (0 :: Posit256)) -- [(1 - 1) == zero | zero = 0 :: Posit es, es <- Z .. V] - let sqrtTaylor = (funLogDomainReduction funLogTaylor).(/2).(funExp2 funExpTaylor).(/log 2) - print $ "sqrt phi using a Taylor algorithm: " ++ show (sqrtTaylor phi) - let sqrtTuma = (funLogDomainReduction funLogTuma).(/2).(funExp2 funExpTuma).(/log 2) - print $ "sqrt phi using a Tuma algorithm: " ++ show (sqrtTuma phi) - print $ "Tuma is fasta: " ++ show (sqrtTaylor (1/1000000) - sqrtTuma (1/1000000)) + print $ "truth :: Posit256: " ++ show (1.2720196495140689642524224617374914917156080418400962486166403825392975755360680118303842149884602585385141476367280265057103381 :: Posit256) + print $ "sqrt phi :: Posit256: " ++ show (sqrt phi :: Posit256) + print $ "sqrt phi :: Posit128: " ++ show (sqrt phi :: Posit128) + print $ "sqrt phi :: Posit64: " ++ show (sqrt phi :: Posit64) + print $ "sqrt phi :: Posit32: " ++ show (sqrt phi :: Posit32) + print $ "sqrt phi :: Posit16: " ++ show (sqrt phi :: Posit16) + print $ "sqrt phi :: Posit8: " ++ show (sqrt phi :: Posit8) + print $ "truth :: P256: " ++ show (1.2720196495140689642524224617374914917156080418400962486166403825392975755360680118303842149884602585385141476367280265057103381 :: P256) + print $ "sqrt phi :: P256: " ++ show (sqrt phi :: P256) + print $ "sqrt phi :: P128: " ++ show (sqrt phi :: P128) + print $ "sqrt phi :: P64: " ++ show (sqrt phi :: P64) + print $ "sqrt phi :: P32: " ++ show (sqrt phi :: P32) + print $ "sqrt phi :: P16: " ++ show (sqrt phi :: P16) + print $ "sqrt phi :: P8: " ++ show (sqrt phi :: P8) {- let truthPosit256 = 0.8956731517052878608869612167009786079379812529831641161347143256 :: Posit256 -- 0.89566032673209158354178209470474131001971567786620187475744721557 :: Posit256 -- 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: Posit256 let truthP256 = 0.8956731517052878608869612167009786079379812529831641161347143256 :: P256 -- 0.89566032673209158354178209470474131001971567786620187475744721557 :: P256 -- 0.8956731517052878608869612167009786079379812529831641161347143256836782657295966290940929214799036260987761959338755143914935872 :: P256 diff --git a/test/TestReadShow.hs b/test/TestReadShow.hs new file mode 100644 --- /dev/null +++ b/test/TestReadShow.hs @@ -0,0 +1,178 @@ + +-------------------------------------------------------------------------------------------- +-- | Posit Numbers +-- Copyright : (C) 2022-2023 Nathan Waivio +-- License : BSD3 +-- Maintainer : Nathan Waivio +-- Stability : Stable +-- Portability : Portable +-- +-- Test Suite for a Library implementing standard Posit Numbers +-- +--------------------------------------------------------------------------------------------- + +-- Hmm... phi is the most irrational number? + +{-# LANGUAGE ScopedTypeVariables #-} -- To reduce some code duplication, this is important +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE TypeApplications #-} + +import Posit +import Posit.Internal.PositC + +import Data.List (partition) + +import Data.Ratio ((%)) -- Import the Rational Numbers ℚ (u+211A), ℚ can get arbitrarily close to Real numbers ℝ (u+211D), used for some of the Transcendental Functions + +main :: IO () +main = do + let (pPosit8, fPosit8) = partition testCase [ x | x <- enumFrom (NaR :: Posit8)] + lpPosit8 = length pPosit8 :: Int -- Int + lfPosit8 = length fPosit8 :: Int -- Int + print $ "Fraction Error (read.show) :: Posit8 : " ++ show (lfPosit8 % (lpPosit8 + lfPosit8)) + print $ "Failure case: " ++ show fPosit8 + let (pP8, fP8) = partition testCase [ x | x <- enumFrom (NaR :: P8)] + lpP8 = length pP8 :: Int + lfP8 = length fP8 :: Int + print $ "Fraction Error (read.show) :: P8 : " ++ show (lfP8 % (lpP8 + lfP8)) + print $ "Failure case: " ++ show fP8 + let (pPosit16, fPosit16) = partition testCase [ x | x <- enumFrom (NaR :: Posit16)] + lpPosit16 = length pPosit16 :: Int -- Int + lfPosit16 = length fPosit16 :: Int -- Int + print $ "Fraction Error (read.show) :: Posit16 : " ++ show (lfPosit16 % (lpPosit16 + lfPosit16)) + print $ "Failure case: " ++ show fPosit16 + let (pP16, fP16) = partition testCase [ x | x <- enumFrom (NaR :: P16)] + lpP16 = length pP16 :: Int + lfP16 = length fP16 :: Int + print $ "Fraction Error (read.show) :: P16 : " ++ show (lfP16 % (lpP16 + lfP16)) + print $ "Failure case: " ++ show fP16 + -- Larger size can only test percentage + -- + let lfPosit32 = length [ x | x <- enumFrom (NaR :: Posit32), not (testCase x)] + print $ "Fraction Error (read.show) :: Posit32 : " ++ show (lfPosit32 % 2^32) + let lfP32 = length [ x | x <- enumFrom (NaR :: P32), not (testCase x)] + print $ "Fraction Error (read.show) :: P32 : " ++ show (lfP32 % 2^32) + -- + + +testCase :: forall es. PositC es => Posit es -> Bool +testCase x = x == ((read @(Posit es)).show $ x) + + +-- * Gustafson's coorispondance with Alessandro +-- +-- Decimal accuracy is not the same as the minimal number of decimals +-- needed to assure conversion back to the original posit. The latter is +-- larger. For example, a 32-bit posit has 28 bits of significance in +-- its highest-accuracy region, which is about: +-- > 28 * log10(2) + 0.5 log10(2) ≈ 8.6 +-- decimals of minimum relative accuracy, but if you only use 9 decimals +-- to express a posit and you convert those 9 decimals back into a +-- posit, occasionally you will get a different posit from what you +-- started with. + +-- * And then responding again to Alessandro + +-- First of all, there are two regime bits (01 or 10), not one. The +-- termination (opposite) bit is part of the regime. That leaves 27 bits +-- for the fraction, which is 28 bits of significand including the +-- hidden value left of the radix point. + +-- Yes, your formula should then work, but remember that you have +-- wobbling relative accuracy (wobble of ±½ log10(2) ≈ ±0.15 decimals) +-- for both posits and floats within any given binade. You should +-- probably be pessimistic and use the lower bound on the sawtooth shape +-- if you are trying to prove something about the accuracy of a +-- computation. + + +-- * During test of "round tripping" a Decimal text string, a rounding +-- issue was found. + +-- After resolving a rounding error that added many failure cases when +-- there should have been rounding for when `2 * remainingSignificand - 1 > 1` +-- for Posits with no fraction. +-- > ... +-- > else if 2 * (r - 1) > 1 -- Previously was 0, this should be "banker's rounding" as per Gustafson (upcoming Ch3), where if `2*(r - 1) > 1` then 1 else 0 +-- > then 1 +-- > else 0 + + +-- The baseline results using `decimalPrec = fromIntegral $ 2 * (nBytes @es) + 1` -- [3,5,9,17,33,65] +-- Benchmark test-posit-readShowId: RUNNING... +-- "Fraction Error (read.show) :: Posit8 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P8 : 1 % 128" +-- "Failure case: [-2.441e-4,2.441e-4]" +-- "Fraction Error (read.show) :: Posit16 : 1 % 8192" +-- "Failure case: [-2.68435e8,-3.35544e7,-5.96046e-8,-2.98023e-8,2.98023e-8,5.96046e-8,3.35544e7,2.68435e8]" +-- "Fraction Error (read.show) :: P16 : 3 % 16384" +-- "Failure case: [-1.40737e14,-7.03687e13,-5.68434e-14,-2.84217e-14,-3.55271e-15,-8.88178e-16,8.88178e-16,3.55271e-15,2.84217e-14,5.68434e-14,7.03687e13,1.40737e14]" +-- "Fraction Error (read.show) :: Posit32 : 37904825 % 2147483648" +-- "Fraction Error (read.show) :: P32 : 37904825 % 2147483648" +-- Benchmark test-posit-readShowId: FINISH +-- +-- After bug fix: +-- Benchmark test-posit-readShowId: RUNNING... +-- "Fraction Error (read.show) :: Posit8 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P8 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: Posit16 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P16 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: Posit32 : 37904819 % 2147483648" +-- "Fraction Error (read.show) :: P32 : 37904819 % 2147483648" +-- Benchmark test-posit-readShowId: FINISH + + +-- Gustafson's Decimal Accuracy: -- [] +-- > decimalPrec posit = +-- > let regimeFormat = findRegimeFormat @es posit +-- > regimeCount = countRegimeBits @es regimeFormat posit +-- > fractionSize = max (fromIntegral (nBits @es) - fromIntegral (signBitSize @es) - regimeCount - fromIntegral (exponentSize @es)) 0 -- fractionSize is at least zero +-- > in ceiling $ (fromIntegral fractionSize + 1) * log10Of2 + halflog10Of2 +-- Benchmark test-posit-readShowId: RUNNING... +-- "Fraction Error (read.show) :: Posit8 : 1 % 128" +-- "Failure case: [-3.1e-2,3.12e-2]" +-- "Fraction Error (read.show) :: P8 : 7 % 128" +-- "Failure case: [-0.12,-2.4e-4,-1.2e-4,-6.1e-5,-1.5e-5,-3.8e-6,-9.5e-7,9.5e-7,3.8e-6,1.5e-5,6.1e-5,1.2e-4,2.44e-4,0.12]" +-- "Fraction Error (read.show) :: Posit16 : 3 % 65536" +-- "Failure case: [-6.7e7,5.96e-8,6.7e7]" +-- "Fraction Error (read.show) :: P16 : 23 % 65536" +-- "Failure case: [-7.2e16,-4.5e15,-1.1e15,-2.8e14,-1.4e14,-7.0e13,-3.5e13,-2.8e-14,-1.4e-14,-7.1e-15,-2.2e-16,2.2e-16,7.1e-15,1.4e-14,2.8e-14,5.68e-14,3.5e13,7.0e13,1.4e14,2.8e14,1.1e15,4.5e15,7.2e16]" +-- "Fraction Error (read.show) :: Posit32 : 75809653 % 4294967296" +-- "Fraction Error (read.show) :: P32 : 75809653 % 4294967296" +-- Benchmark test-posit-readShowId: FINISH +-- +-- After bug fix: +-- Benchmark test-posit-readShowId: RUNNING... +-- "Fraction Error (read.show) :: Posit8 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P8 : 1 % 128" +-- "Failure case: [-0.12,0.12]" +-- "Fraction Error (read.show) :: Posit16 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P16 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: Posit32 : 37904819 % 2147483648" +-- "Fraction Error (read.show) :: P32 : 37904819 % 2147483648" +-- Benchmark test-posit-readShowId: FINISH +-- +-- After bug fix and incrementing by 1 decimal digit +-- Benchmark test-posit-readShowId: RUNNING... +-- "Fraction Error (read.show) :: Posit8 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P8 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: Posit16 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: P16 : 0 % 1" +-- "Failure case: []" +-- "Fraction Error (read.show) :: Posit32 : 0 % 1" +-- "Fraction Error (read.show) :: P32 : 0 % 1" +-- Benchmark test-posit-readShowId: FINISH + + +