ref: 51f905138babfc67f3291a40c8b5ddd7ec85b73d
dir: /lib/Numeric/FormatFloat.hs/
module Numeric.FormatFloat(
showFloat,
formatRealFloat, formatRealFloatAlt,
showEFloat,
showFFloat,
showGFloat,
showFFloatAlt,
showGFloatAlt,
showHFloat,
) where
import Data.Char
import Numeric.Show
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