diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -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 diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -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 diff --git a/posit.cabal b/posit.cabal --- a/posit.cabal +++ b/posit.cabal @@ -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 homepage: https://github.com/waivio/posit#readme bug-reports: https://github.com/waivio/posit/issues diff --git a/src/Posit.hs b/src/Posit.hs --- a/src/Posit.hs +++ b/src/Posit.hs @@ -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 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 @@ -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 diff --git a/test/TestPosit.hs b/test/TestPosit.hs --- a/test/TestPosit.hs +++ b/test/TestPosit.hs @@ -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 --