ref: 41c2cd888668f697dbcae2b031c7afacde168987
parent: f7236224b0e0d091ab0bd9f40017e527b96f87a3
author: Lennart Augustsson <lennart@augustsson.net>
date: Wed Apr 3 05:38:44 EDT 2024
Add float formatting
--- a/lib/AllOfLib.hs
+++ b/lib/AllOfLib.hs
@@ -80,6 +80,7 @@
import GHC.Stack
import GHC.Types
import Numeric
+import Numeric.FormatFloat
import Numeric.Natural
import Prelude
import Primitives
--- /dev/null
+++ b/lib/Numeric/FormatFloat.hs
@@ -1,0 +1,304 @@
+module Numeric.FormatFloat(
+ formatRealFloat, formatRealFloatAlt,
+
+ showEFloat,
+ showFFloat,
+ showGFloat,
+ showFFloatAlt,
+ showGFloatAlt,
+ showHFloat,
+
+ ) where
+import Data.Char
+import Numeric
+
+--showFloat :: (RealFloat a) => a -> ShowS
+--showFloat x = showString (formatRealFloat FFGeneric Nothing x)
+
+-- These are the format types. This type is not exported.
+data FFFormat = FFExponent | FFFixed | FFGeneric
+
+formatRealFloat :: (RealFloat a) => FFFormat -> Maybe Int -> a -> String
+formatRealFloat fmt decs x = formatRealFloatAlt fmt decs False x
+
+formatRealFloatAlt :: (RealFloat a) => FFFormat -> Maybe Int -> Bool -> a
+ -> String
+formatRealFloatAlt fmt decs alt x
+ | isNaN x = "NaN"
+ | isInfinite x = if x < 0 then "-Infinity" else "Infinity"
+ | x < 0 || isNegativeZero x = '-':doFmt fmt (floatToDigits (toInteger base) (-x))
+ | otherwise = doFmt fmt (floatToDigits (toInteger base) x)
+ where
+ base = 10
+
+ doFmt format (is, e) =
+ let ds = map intToDigit is in
+ case format of
+ FFGeneric ->
+ doFmt (if e < 0 || e > 7 then FFExponent else FFFixed)
+ (is,e)
+ FFExponent ->
+ case decs of
+ Nothing ->
+ let show_e' = show (e-1) in
+ case ds of
+ "0" -> "0.0e0"
+ [d] -> d : ".0e" ++ show_e'
+ (d:ds') -> d : '.' : ds' ++ "e" ++ show_e'
+ [] -> error "formatRealFloat/doFmt/FFExponent: []"
+ Just d | d <= 0 ->
+ -- handle this case specifically since we need to omit the
+ -- decimal point as well (#15115).
+ -- Note that this handles negative precisions as well for consistency
+ -- (see #15509).
+ case is of
+ [0] -> "0e0"
+ _ ->
+ let
+ (ei,is') = roundTo base 1 is
+ n:_ = map intToDigit (if ei > 0 then init is' else is')
+ in n : 'e' : show (e-1+ei)
+ Just dec ->
+ let dec' = max dec 1 in
+ case is of
+ [0] -> '0' :'.' : take dec' (repeat '0') ++ "e0"
+ _ ->
+ let
+ (ei,is') = roundTo base (dec'+1) is
+ (d:ds') = map intToDigit (if ei > 0 then init is' else is')
+ in
+ d:'.':ds' ++ 'e':show (e-1+ei)
+ FFFixed ->
+ let
+ mk0 ls = case ls of { "" -> "0" ; _ -> ls}+ in
+ case decs of
+ Nothing
+ | e <= 0 -> "0." ++ replicate (-e) '0' ++ ds
+ | otherwise ->
+ let
+ f 0 s rs = mk0 (reverse s) ++ '.':mk0 rs
+ f n s "" = f (n-1) ('0':s) ""+ f n s (r:rs) = f (n-1) (r:s) rs
+ in
+ f e "" ds
+ Just dec ->
+ let dec' = max dec 0 in
+ if e >= 0 then
+ let
+ (ei,is') = roundTo base (dec' + e) is
+ (ls,rs) = splitAt (e+ei) (map intToDigit is')
+ in
+ mk0 ls ++ (if null rs && not alt then "" else '.':rs)
+ else
+ let
+ (ei,is') = roundTo base dec' (replicate (-e) 0 ++ is)
+ d:ds' = map intToDigit (if ei > 0 then is' else 0:is')
+ in
+ d : (if null ds' && not alt then "" else '.':ds')
+
+
+roundTo :: Int -> Int -> [Int] -> (Int,[Int])
+roundTo base d is =
+ case f d True is of
+ x@(0,_) -> x
+ (1,xs) -> (1, 1:xs)
+ _ -> error "roundTo: bad Value"
+ where
+ b2 = base `quot` 2
+
+ f n _ [] = (0, replicate n 0)
+ f 0 e (x:xs) | x == b2 && e && all (== 0) xs = (0, []) -- Round to even when at exactly half the base
+ | otherwise = (if x >= b2 then 1 else 0, [])
+ f n _ (i:xs)
+ | i' == base = (1,0:ds)
+ | otherwise = (0,i':ds)
+ where
+ (c,ds) = f (n-1) (even i) xs
+ i' = c + i
+
+-- Based on "Printing Floating-Point Numbers Quickly and Accurately"
+-- by R.G. Burger and R.K. Dybvig in PLDI 96.
+-- This version uses a much slower logarithm estimator. It should be improved.
+
+-- | 'floatToDigits' takes a base and a non-negative 'RealFloat' number,
+-- and returns a list of digits and an exponent.
+-- In particular, if @x>=0@, and
+--
+-- > floatToDigits base x = ([d1,d2,...,dn], e)
+--
+-- then
+--
+-- (1) @n >= 1@
+--
+-- (2) @x = 0.d1d2...dn * (base**e)@
+--
+-- (3) @0 <= di <= base-1@
+
+floatToDigits :: (RealFloat a) => Integer -> a -> ([Int], Int)
+floatToDigits _ 0 = ([0], 0)
+floatToDigits base x =
+ let
+ (f0, e0) = decodeFloat x
+ (minExp0, _) = floatRange x
+ p = floatDigits x
+ b = floatRadix x
+ minExp = minExp0 - p -- the real minimum exponent
+ -- Haskell requires that f be adjusted so denormalized numbers
+ -- will have an impossibly low exponent. Adjust for this.
+ (f, e) =
+ let n = minExp - e0 in
+ if n > 0 then (f0 `quot` (expt b n), e0+n) else (f0, e0)
+ (r, s, mUp, mDn) =
+ if e >= 0 then
+ let be = expt b e in
+ if f == expt b (p-1) then
+ (f*be*b*2, 2*b, be*b, be) -- according to Burger and Dybvig
+ else
+ (f*be*2, 2, be, be)
+ else
+ if e > minExp && f == expt b (p-1) then
+ (f*b*2, expt b (-e+1)*2, b, 1)
+ else
+ (f*2, expt b (-e)*2, 1, 1)
+ k :: Int
+ k =
+ let
+ k0 :: Int
+ k0 =
+ if b == 2 && base == 10 then
+ -- logBase 10 2 is very slightly larger than 8651/28738
+ -- (about 5.3558e-10), so if log x >= 0, the approximation
+ -- k1 is too small, hence we add one and need one fixup step less.
+ -- If log x < 0, the approximation errs rather on the high side.
+ -- That is usually more than compensated for by ignoring the
+ -- fractional part of logBase 2 x, but when x is a power of 1/2
+ -- or slightly larger and the exponent is a multiple of the
+ -- denominator of the rational approximation to logBase 10 2,
+ -- k1 is larger than logBase 10 x. If k1 > 1 + logBase 10 x,
+ -- we get a leading zero-digit we don't want.
+ -- With the approximation 3/10, this happened for
+ -- 0.5^1030, 0.5^1040, ..., 0.5^1070 and values close above.
+ -- The approximation 8651/28738 guarantees k1 < 1 + logBase 10 x
+ -- for IEEE-ish floating point types with exponent fields
+ -- <= 17 bits and mantissae of several thousand bits, earlier
+ -- convergents to logBase 10 2 would fail for long double.
+ -- Using quot instead of div is a little faster and requires
+ -- fewer fixup steps for negative lx.
+ let lx = p - 1 + e0
+ k1 = (lx * 8651) `quot` 28738
+ in if lx >= 0 then k1 + 1 else k1
+ else
+ -- f :: Integer, log :: Float -> Float,
+ -- ceiling :: Float -> Int
+ ceiling ((log (fromInteger (f+1) :: Float) +
+ fromIntegral e * log (fromInteger b)) /
+ log (fromInteger base))
+--WAS: fromInt e * log (fromInteger b))
+
+ fixup n =
+ if n >= 0 then
+ if r + mUp <= expt base n * s then n else fixup (n+1)
+ else
+ if expt base (-n) * (r + mUp) <= s then n else fixup (n+1)
+ in
+ fixup k0
+
+ gen ds rn sN mUpN mDnN =
+ let
+ (dn, rn') = (rn * base) `quotRem` sN
+ mUpN' = mUpN * base
+ mDnN' = mDnN * base
+ in
+ case (rn' < mDnN', rn' + mUpN' > sN) of
+ (True, False) -> dn : ds
+ (False, True) -> dn+1 : ds
+ (True, True) -> if rn' * 2 < sN then dn : ds else dn+1 : ds
+ (False, False) -> gen (dn:ds) rn' sN mUpN' mDnN'
+
+ rds =
+ if k >= 0 then
+ gen [] r (s * expt base k) mUp mDn
+ else
+ let bk = expt base (-k) in
+ gen [] (r * bk) s (mUp * bk) (mDn * bk)
+ in
+ (map fromIntegral (reverse rds), k)
+
+-----------
+
+showEFloat :: (RealFloat a) => Maybe Int -> a -> ShowS
+showFFloat :: (RealFloat a) => Maybe Int -> a -> ShowS
+showGFloat :: (RealFloat a) => Maybe Int -> a -> ShowS
+showFFloatAlt :: (RealFloat a) => Maybe Int -> a -> ShowS
+showGFloatAlt :: (RealFloat a) => Maybe Int -> a -> ShowS
+
+showEFloat d x = showString (formatRealFloat FFExponent d x)
+showFFloat d x = showString (formatRealFloat FFFixed d x)
+showGFloat d x = showString (formatRealFloat FFGeneric d x)
+showFFloatAlt d x = showString (formatRealFloatAlt FFFixed d True x)
+showGFloatAlt d x = showString (formatRealFloatAlt FFGeneric d True x)
+
+showHFloat :: RealFloat a => a -> ShowS
+showHFloat = showString . fmt
+ where
+ fmt x
+ | isNaN x = "NaN"
+ | isInfinite x = (if x < 0 then "-" else "") ++ "Infinity"
+ | x < 0 || isNegativeZero x = '-' : cvt (-x)
+ | otherwise = cvt x
+
+ cvt x
+ | x == 0 = "0x0p+0"
+ | otherwise =
+ case floatToDigits 2 x of
+ r@([], _) -> error $ "Impossible happened: showHFloat: " ++ show r
+ (d:ds, e) -> "0x" ++ show d ++ frac ds ++ "p" ++ show (e-1)
+
+ -- Given binary digits, convert them to hex in blocks of 4
+ -- Special case: If all 0's, just drop it.
+ frac digits
+ | allZ digits = ""
+ | otherwise = "." ++ hex digits
+ where
+ hex ds =
+ case ds of
+ [] -> ""
+ [a] -> hexDigit a 0 0 0 ""
+ [a,b] -> hexDigit a b 0 0 ""
+ [a,b,c] -> hexDigit a b c 0 ""
+ a : b : c : d : r -> hexDigit a b c d (hex r)
+
+ hexDigit a b c d = showHex (8*a + 4*b + 2*c + d)
+
+ allZ xs = case xs of
+ x : more -> x == 0 && allZ more
+ [] -> True
+
+{-+-- Exponentiation with a cache for the most common numbers.
+minExpt, maxExpt :: Int
+minExpt = 0
+maxExpt = 1100
+
+expt :: Integer -> Int -> Integer
+expt base n =
+ if base == 2 && n >= minExpt && n <= maxExpt then
+ expts!n
+ else
+ if base == 10 && n <= maxExpt10 then
+ expts10!n
+ else
+ base^n
+
+expts :: Array Int Integer
+expts = array (minExpt,maxExpt) [(n,2^n) | n <- [minExpt .. maxExpt]]
+
+maxExpt10 :: Int
+maxExpt10 = 324
+
+expts10 :: Array Int Integer
+expts10 = array (minExpt,maxExpt10) [(n,10^n) | n <- [minExpt .. maxExpt10]]
+-}
+expt :: Integer -> Int -> Integer
+expt base n = base^n
--- a/lib/Prelude.hs
+++ b/lib/Prelude.hs
@@ -52,7 +52,7 @@
import Data.Eq(Eq(..))
import Data.Floating(Floating(..))
import Data.Fractional(Fractional(..), (^^))
-import Data.Function(id, const, (.), flip, ($), seq, ($!), until, curry, uncurry)
+import Data.Function(id, const, (.), flip, ($), seq, ($!), until, curry, uncurry, asTypeOf)
import Data.Functor(Functor(..), (<$>))
import Data.Int(Int)
import Data.Int.Instances
--
⑨