ref: de1a985cb74165f0b8689e5ba52391db288e4369
parent: f63e6ed90adf61b213ff5ffe3ef5bf16560cdf2f
author: Lennart Augustsson <lennart@augustsson.net>
date: Sun Nov 5 12:39:03 EST 2023
Add RealFloat class.
--- a/TODO
+++ b/TODO
@@ -47,3 +47,5 @@
Bugs
* Removing [] from prim table
* :reload doesn't show error message
+ * default methods not exported
+
\ No newline at end of file
--- a/lib/Data/Double.hs
+++ b/lib/Data/Double.hs
@@ -1,17 +1,21 @@
-- Copyright 2023 Lennart Augustsson
-- See LICENSE file for full license.
-module Data.Double(module Data.Double, Double) where
+module Data.Double(Double) where
import Primitives
import Control.Error
-import Data.Bool_Type
+import Data.Bits
+import Data.Bool
import Data.Eq
import Data.Floating
import Data.Fractional
+import Data.Function
import Data.Integer
import Data.Ord
import Data.Ratio
import Data.Real
+import Data.RealFloat
import Data.Num
+import Data.Word
import Text.Show
instance Num Double where
@@ -43,14 +47,14 @@
(>) = primDoubleGT
(>=) = primDoubleGE
--- | this primitive will print doubles with up to 6 decimal points
--- it turns out that doubles are extremely tricky, and just printing them is a
--- herculean task of its own...
+-- For now, cheat and call C
instance Show Double where
show = primDoubleShow
instance Real Double where
- toRational _ = error "Double.toRational not implemented"
+ toRational x =
+ let (m, e) = decodeDouble x
+ in toRational m * 2^^e
instance Floating Double where
pi = 3.141592653589793
@@ -73,3 +77,67 @@
foreign import ccall "asin" casin :: Double -> IO Double
foreign import ccall "acos" cacos :: Double -> IO Double
foreign import ccall "atan" catan :: Double -> IO Double
+
+-- Assumes 64 bit floats
+instance RealFloat Double where
+ floatRadix _ = 2
+ floatDigits _ = 53
+ floatRange _ = (-1021,1024)
+ decodeFloat = decodeDouble
+ encodeFloat = encodeDouble
+ isNaN = isNaNDouble
+ isInfinite = isInfDouble
+ isDenormalized = isDenDouble
+ isNegativeZero = isNegZeroDouble
+ isIEEE _ = True
+
+decodeDouble :: Double -> (Integer, Int)
+decodeDouble x =
+ let xw = primWordFromDoubleRaw x
+ sign = xw .&. 0x8000000000000000
+ expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
+ mant = xw .&. 0x000fffffffffffff
+ neg = if sign /= 0 then negate else id
+ in if expn == 0 then
+ -- subnormal or 0
+ (neg (_wordToInteger mant), 0)
+ else if expn == 0x7ff then
+ -- infinity or NaN
+ (0, 0)
+ else
+ -- ordinary number, add hidden bit
+ -- mant is offset-1023, and assumes scaled mantissa (thus -52)
+ (neg (_wordToInteger (mant .|. 0x0010000000000000)),
+ primWordToInt expn - 1023 - 52)
+
+isNaNDouble :: Double -> Bool
+isNaNDouble x =
+ let xw = primWordFromDoubleRaw x
+ expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
+ mant = xw .&. 0x000fffffffffffff
+ in expn == 0x7ff && mant /= 0
+
+isInfDouble :: Double -> Bool
+isInfDouble x =
+ let xw = primWordFromDoubleRaw x
+ expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
+ mant = xw .&. 0x000fffffffffffff
+ in expn == 0x7ff && mant == 0
+
+isDenDouble :: Double -> Bool
+isDenDouble x =
+ let xw = primWordFromDoubleRaw x
+ expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
+ mant = xw .&. 0x000fffffffffffff
+ in expn == 0 && mant /= 0
+
+isNegZeroDouble :: Double -> Bool
+isNegZeroDouble x =
+ let xw = primWordFromDoubleRaw x
+ sign = xw .&. 0x8000000000000000
+ rest = xw .&. 0x7fffffffffffffff
+ in sign /= 0 && rest == 0
+
+-- Simple (and sometimes wrong) encoder
+encodeDouble :: Integer -> Int -> Double
+encodeDouble mant expn = fromInteger mant * 2^^expn
--- /dev/null
+++ b/lib/Data/RealFloat.hs
@@ -1,0 +1,69 @@
+module Data.RealFloat(
+ --RealFloat(..), BUG
+ module Data.RealFloat,
+ ) where
+import Primitives
+import Data.Bool
+import Data.Eq
+import Data.Floating
+import Data.Fractional
+import Data.Int
+import Data.Integer
+import Data.Num
+import Data.Ord
+
+class (Fractional a, Ord a, Floating a) => RealFloat a where
+ floatRadix :: a -> Integer
+ floatDigits :: a -> Int
+ floatRange :: a -> (Int,Int)
+ decodeFloat :: a -> (Integer,Int)
+ encodeFloat :: Integer -> Int -> a
+ exponent :: a -> Int
+ significand :: a -> a
+ scaleFloat :: Int -> a -> a
+ isNaN :: a -> Bool
+ isInfinite :: a -> Bool
+ isDenormalized :: a -> Bool
+ isNegativeZero :: a -> Bool
+ isIEEE :: a -> Bool
+ atan2 :: a -> a -> a
+
+ exponent x = if m == 0 then 0 else n + floatDigits x
+ where (m,n) = decodeFloat x
+
+ significand x = encodeFloat m (negate (floatDigits x))
+ where (m,_) = decodeFloat x
+
+ scaleFloat 0 x = x
+ scaleFloat k x
+ | isFix = x
+ | otherwise = encodeFloat m (n + clamp b k)
+ where (m,n) = decodeFloat x
+ (l,h) = floatRange x
+ d = floatDigits x
+ b = h - l + 4*d
+ -- n+k may overflow, which would lead
+ -- to wrong results, hence we clamp the
+ -- scaling parameter.
+ -- If n + k would be larger than h,
+ -- n + clamp b k must be too, similar
+ -- for smaller than l - d.
+ -- Add a little extra to keep clear
+ -- from the boundary cases.
+ isFix = x == 0 || isNaN x || isInfinite x
+
+ atan2 y x
+ | x > 0 = atan (y/x)
+ | x == 0 && y > 0 = pi/2
+ | x < 0 && y > 0 = pi + atan (y/x)
+ |(x <= 0 && y < 0) ||
+ (x < 0 && isNegativeZero y) ||
+ (isNegativeZero x && isNegativeZero y)
+ = negate (atan2 (negate y) x)
+ | y == 0 && (x < 0 || isNegativeZero x)
+ = pi -- must be after the previous test on zero y
+ | x==0 && y==0 = y -- must be after the other double zero tests
+ | otherwise = x + y -- x or y is a NaN, return a NaN (via +)
+
+clamp :: Int -> Int -> Int
+clamp bd k = max (negate bd) (min bd k)
--- a/lib/Prelude.hs
+++ b/lib/Prelude.hs
@@ -24,6 +24,7 @@
module Data.Ord,
module Data.Ratio,
module Data.Real,
+ module Data.RealFloat,
module Data.Tuple,
module System.IO,
module Text.Show,
@@ -54,6 +55,7 @@
import Data.Ord
import Data.Ratio(Rational)
import Data.Real
+import Data.RealFloat
import Data.Tuple
import System.IO
import Text.Show
--
⑨