posit 3.2.0.2 → 3.2.0.3
raw patch · 6 files changed
+87/−36 lines, 6 files
Files
- ChangeLog.md +11/−0
- README.md +2/−2
- posit.cabal +1/−1
- src/Posit.hs +65/−25
- src/Posit/Internal/PositC.hs +5/−7
- test/TestPosit.hs +3/−1
ChangeLog.md view
@@ -1,5 +1,16 @@ # Changelog for Posit Numbers +# posit-3.2.0.3++ * Made the following changes in anticipation of adding the 2022 Posit Standard:+ * Made the `IntN` type family non-Injective, and added more visable type applications to help the compiler select the proper types+ * Corrected some bad uses of `nBytes @es`, with `2^(exponentSize @es)`, in order to be more general+ * Chagned `maxPosRat` to match the more general form as described in "Posit Arithmetic" (John L Gustafson, 10 October 2017)+ * Changed `lnOf2` to be a long decimal value, in order to be more general+ * Changed Borwein's algorithm, with quintic convergence, to check for a fixed point of both `a` and `s`+ * Added Borwein's Quadradic 1985+ * Added Borewein's Cubic+ # posit-3.2.0.2 * Added `FlexableContexts` back in to Posit.hs, a build error occured on GHC-9.2 that didn't occur with GHC-9.0 or GHC-8.10
README.md view
@@ -1,4 +1,4 @@-# posit 3.2.0.2+# posit 3.2.0.3 The [Posit Standard 3.2](https://posithub.org/docs/posit_standard.pdf), where Real numbers are approximated by Maybe Rational. The Posit @@ -25,7 +25,7 @@ * 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.2` es is instantiated as Z, I,+word size are related. In `posit-3.2.0.3` 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
posit.cabal view
@@ -1,7 +1,7 @@ cabal-version: 1.12 name: posit-version: 3.2.0.2+version: 3.2.0.3 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
src/Posit.hs view
@@ -8,7 +8,7 @@ -- Portability : Portable -- -- | Library implementing standard Posit Numbers (Posit Standard version--- 3.2.0.0, with some improvements) a fixed width word size of+-- 3.2, with some improvements) a fixed width word size of -- 2^es bytes. -- ---------------------------------------------------------------------------------------------@@ -111,6 +111,8 @@ funPi2, funPi3, funPi4,+ funPi5,+ funPi6, funPsiSha1, funPsiSha2, funPsiSha3@@ -168,16 +170,16 @@ 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+pattern NaR :: forall es. PositC es => Posit es+pattern NaR <- (Posit (decode @es -> Nothing)) where+ NaR = Posit (unReal @es) -- -- -- |A Real or at least Rational Number, rounded to the nearest Posit Rational representation-pattern R :: PositC es => Rational -> Posit es-pattern R r <- (Posit (decode -> Just r)) where- R r = Posit (encode $ Just r)+pattern R :: forall es. PositC es => Rational -> Posit es+pattern R r <- (Posit (decode @es -> Just r)) where+ R r = Posit (encode @es $ Just r) -- -- Posit functions are complete if the following two patterns are completely defined.@@ -233,7 +235,7 @@ -- '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)+ fromInteger int = R $ fromInteger int -- 'negate', Negates the sign of the directionality. negate of a posit is the same as negate of the integer representation negate = viaIntegral negate --@@ -365,9 +367,9 @@ -- I'm bound to want this definition: instance PositC es => Bounded (Posit es) where -- 'minBound' the most negative number represented- minBound = Posit mostNegVal+ minBound = Posit (mostNegVal @es) -- 'maxBound' the most positive number represented- maxBound = Posit mostPosVal+ maxBound = Posit (mostPosVal @es) -- @@ -437,11 +439,11 @@ -- Fuse Sum of 4 Posits fsum4 = viaRational4 fsum4 -- Fuse Sum of a List- fsumL (toList -> l) = Posit $ encode (Just $ go l 0)+ fsumL (toList -> l) = Posit $ encode @es (Just $ go l 0) where go :: [Posit es] -> Rational -> Rational go [] !acc = acc- go ((Posit int) : xs) !acc = case decode int of+ go ((Posit int) : xs) !acc = case decode @es int of Nothing -> error "Posit List contains NaR" Just r -> go xs (acc + r) -- Fuse Dot Product of a 3-Vector@@ -449,14 +451,14 @@ -- 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)+ fdotL (toList -> l1) (toList -> l2) = Posit $ encode @es (Just $ go l1 l2 0) where go [] [] !acc = acc go [] _ _ = error "Lists not the same length" go _ [] _ = error "Lists not the same length"- go ((Posit int1) : bs) ((Posit int2) : cs) !acc = case decode int1 of+ go ((Posit int1) : bs) ((Posit int2) : cs) !acc = case decode @es int1 of Nothing -> error "First Posit List contains NaR"- Just r1 -> case decode int2 of+ Just r1 -> case decode @es int2 of Nothing -> error "Second Posit List contains NaR" Just r2 -> go bs cs (acc + (r1 * r2)) --@@ -497,7 +499,7 @@ -- instance PositC es => AltShow (Posit es) where- displayBinary (Posit int) = displayBin int+ displayBinary (Posit int) = displayBin @es int displayIntegral (Posit int) = show int @@ -901,7 +903,7 @@ | 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)+ (ex, sig) = (int * fromIntegral (2^(exponentSize @V)) + fromIntegral nat + 1, fromRational rat / 2) (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x' -- sign should always be positive x' = 1 - x t = (2 - sig - x') / (2 + sig - x')@@ -1093,7 +1095,7 @@ -- -- Use the constant, for performance lnOf2 :: Posit256-lnOf2 = Posit 28670435363615573179632300308403400109260626501925370561166468529302554498548+lnOf2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875420014810205706857336855202 -- --@@ -1122,7 +1124,7 @@ | 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)+ (ex, sig) = (int * fromIntegral (2^(exponentSize @V)) + fromIntegral nat, fromRational rat) (_,int,nat,rat) = (posit2TupPosit @V).toRational $ z -- sign should always be positive findSquaring m s | s >= 2 && s < 4 = (m, s/2)@@ -1150,18 +1152,19 @@ -- 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)+funPi2 = recip $ go 0 0 0 0.5 (5 / phi^3) where- go :: Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256- go !prev !n !a !s- | prev == a = a+ go :: Posit256 -> Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256+ go !prevA !prevS !n !a !s+ | prevA == a = a+ | prevS == s = a | otherwise = let x = 5 / s - 1 y = (x - 1)^2 + 7 z = (0.5 * x * (y + sqrt (y^2 - 4 * x^3)))**(1/5) a' = s^2 * a - (5^n * ((s^2 - 5)/2 + sqrt (s * (s^2 - 2*s + 5)))) s' = 25 / ((z + x/z + 1)^2 * s)- in go a (n+1) (trace (show a') a') s'+ in go a s (n+1) (trace ("ΔA: " ++ show (a' - a)) a') (trace ("ΔS: " ++ show (s' - s)) s') -- #endif @@ -1192,7 +1195,44 @@ -- +-- Borwin's Quadradic Alogrithm 1985+funPi5 :: Posit256+funPi5 = recip $ go 0 0 1 (6 - 4 * sqrt 2) (sqrt 2 - 1)+ where+ go :: Posit256 -> Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256+ go !prevA !prevY !n a y+ | prevA == a = a+ | prevY == y = a+ | otherwise =+ let f = (1 - y^4)**(1/4)+ y' = (1 - f) / (1 + f)+ a' = a * (1 + y')^4 - 2^(2 * n + 1) * y' * (1 + y' + y'^2) + in if n == 3+ then a'+ else go a y (n+1) (trace ("A: " ++ show a') a') (trace ("Y: " ++ show y') y')+--+-- 3.14159265358979323846264338327950288419716939937510582097494459231+-- ULP: -97 +-- Borwin's Cubic Algirthm+funPi6 :: Posit256+funPi6 = recip $ go 0 0 1 (1/3) ((sqrt 3 - 1) / 2)+ where+ go :: Posit256 -> Posit256 -> Natural -> Posit256 -> Posit256 -> Posit256+ go !prevA !prevS !n !a !s+ | prevA == a = a+ | prevS == s = a+ | otherwise =+ let r = 3 / (1 + 2 * (1 - s^3)**(1/3))+ s'= (r - 1) / 2+ a'= r^2 * a - 3^(n-1) * (r^2 - 1)+ in if n == 4+ then a'+ else go a s (n+1) a' s'+-- 3.14159265358979323846264338327950288419716939937510582097494459231+-- ULP: 216++ -- -- looks to be about 4 ULP accurate at -100, right on the money at -1000 funExp :: Posit256 -> Posit256@@ -1307,7 +1347,7 @@ | 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+ (ex, sig) = (int * fromIntegral (2^(exponentSize @V)) + fromIntegral nat + 1, fromRational rat / 2) -- move significand range from 1,2 to 0.5,1 (_,int,nat,rat) = (posit2TupPosit @V).toRational $ x -- sign should always be positive
src/Posit/Internal/PositC.hs view
@@ -71,11 +71,9 @@ | V -- | 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'+-- Int16, Int32, Int64, Int128, Int256. {-@ embed IntN * as int @-}-type family IntN (es :: ES) = r | r -> es+type family IntN (es :: ES) where IntN Z = Int8 IntN I = Int16@@ -149,7 +147,7 @@ signBitSize = 1 uSeed :: Natural -- ^ 'uSeed' scaling factor for the regime of the Posit Representation- uSeed = 2^(nBytes @es)+ uSeed = 2^2^(exponentSize @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@@ -165,11 +163,11 @@ leastNegVal = -1 mostNegVal :: IntN es- mostNegVal = negate mostPosVal+ mostNegVal = negate (mostPosVal @es) -- Rational Value of common bounds maxPosRat :: Rational- maxPosRat = 2^((nBytes @es) * ((nBits @es) - 2)) % 1+ maxPosRat = (fromIntegral (uSeed @es)^(nBits @es - 2)) % 1 minPosRat :: Rational minPosRat = recip (maxPosRat @es) maxNegRat :: Rational
test/TestPosit.hs view
@@ -90,9 +90,11 @@ 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 "Borwein's Quintic algorithm: pi:" funPi2 truth eval "Bailey–Borwein–Plouffe (BBP) formula: pi:" funPi3 truth eval "Fabrice Bellard improvement on the BBP: pi:" funPi4 truth+ eval "Borwein's Quadradic 1985 formula: pi:" funPi5 truth+ eval "Borwein Cubic: pi:" funPi6 truth eval "Wolfram Alpha: pi:" truth truth eval "Bailey–Borwein–Plouffe (BBP) formula: but succ pi:" (succ funPi3) truth --