content: -- Copyright 2025 Lennart Augustsson
  -- See LICENSE file for full license.
  module Data.Double(Double) where
  import qualified Prelude()              -- do not import Prelude
  import Primitives
  import Control.Error
  import Data.Bits.Base
  import Data.Bool
  import Data.Char
  import Data.Enum
  import Data.Eq
  import Data.Floating
  import Data.Fractional
  import Data.Function
  import Data.Integer
  import Data.Integer.Internal(_integerToDouble, _wordToInteger)
  import Data.Integral
  import Data.List
  import Data.Ord
  import Data.Ratio
  import Data.Real
  import Data.RealFloat
  import Data.RealFrac
  import Data.Num
  import Data.Word.Word64
  import {-# SOURCE #-} Numeric.FormatFloat(showFloat)
  import Numeric.Show(showSignedNeg)
  import Text.Show
  
  --
  -- This is a floating point type that is a wide as the word width of the machine.
  --
  
  instance Num Double where
    (+)  = primDoubleAdd
    (-)  = primDoubleSub
    (*)  = primDoubleMul
    negate = primDoubleNeg
    abs x = if x < 0.0 then - x else x
    signum x =
      case compare x 0.0 of
        LT -> -1.0
        EQ ->  0.0
        GT ->  1.0
    fromInteger = _integerToDouble
  
  instance Fractional Double where
    (/) = primDoubleDiv
    fromRational x =
      let n = numerator x
          d = denominator x
          maxVal = 2 ^ (1023::Int)  -- actually bigger, but this converts correctly
      in  if abs n <= maxVal && d < maxVal then
            fromInteger n / fromInteger d
          else
            let (ne, nm) = scaleToMax maxVal n
                (de, dm) = scaleToMax maxVal d
            in  scaleFloat64 (ne - de) (fromInteger nm / fromInteger dm)
  
  -- XXX Very inefficient scaling
  -- If we get here the number has > 1023 bits, so shifting by 256 is fine.
  scaleToMax :: Integer -> Integer -> (Int, Integer)
  scaleToMax m = f 0
    where f e x = if abs x < m then (e, x) else f (e+256) (x `shiftR` 256)
  
  instance Eq Double where
    (==) = primDoubleEQ
    (/=) = primDoubleNE
  
  instance Ord Double where
    (<)  = primDoubleLT
    (<=) = primDoubleLE
    (>)  = primDoubleGT
    (>=) = primDoubleGE
  
  instance Show Double where
    showsPrec = showSignedNeg (\ x -> x < 0 || isNegZeroFloat64 x) showFloat
  
  {- in Text.Read.Internal
  instance Read Double where
    readsPrec _ = readSigned $ \ r -> [ (primDoubleRead s, t) | (s@(c:_), t) <- lex r, isDigit c ]
  -}
  
  instance Enum Double where
    succ x = x + 1
    pred x = x - 1
    toEnum n = fromIntegral n
    fromEnum = fromInteger . truncate
    enumFrom = numericEnumFrom
    enumFromThen = numericEnumFromThen
    enumFromTo = numericEnumFromTo
    enumFromThenTo = numericEnumFromThenTo
  
  instance Real Double where
    toRational x | isNaN x = rationalNaN
                 | isInfinite x = if x < 0 then -rationalInfinity else rationalInfinity
                 | isNegativeZero x = rationalNegativeZero
                 | otherwise =
      case decodeFloat x of
        (m, e) -> toRational m * 2^^e
  
  instance RealFrac Double where
    properFraction x =
      case decodeFloat x of
        (m, e) ->   -- x = m * 2^e
          let m' | e < 0     = m `quot` 2^(-e)
                 | otherwise = m * 2^e
          in  (fromInteger m', x - fromInteger m')
  
  instance Floating Double where
    pi     = 3.141592653589793
    log  x = clog x
    exp  x = cexp x
    sqrt x = csqrt x
    sin  x = csin x
    cos  x = ccos x
    tan  x = ctan x
    asin x = casin x
    acos x = cacos x
    atan x = catan x
  
  foreign import ccall "log"  clog  :: Double -> Double
  foreign import ccall "exp"  cexp  :: Double -> Double
  foreign import ccall "sqrt" csqrt :: Double -> Double
  foreign import ccall "sin"  csin  :: Double -> Double
  foreign import ccall "cos"  ccos  :: Double -> Double
  foreign import ccall "tan"  ctan  :: Double -> Double
  foreign import ccall "asin" casin :: Double -> Double
  foreign import ccall "acos" cacos :: Double -> Double
  foreign import ccall "atan" catan :: Double -> Double
  foreign import ccall "atan2" catan2 :: Double -> Double -> Double
  foreign import ccall "scalbn" cscalbn :: Double -> Int -> Double
  
  -- Assumes 32/64 bit floats
  instance RealFloat Double where
    floatRadix     _ = 2
    floatDigits    _ = 53
    floatRange     _ = (-1021,1024)
    decodeFloat      = decodeFloat64
    encodeFloat      = encodeFloat64
    isNaN            = isNaNFloat64
    isInfinite       = isInfFloat64
    isDenormalized   = isDenFloat64
    isNegativeZero   = isNegZeroFloat64
    isIEEE         _ = True
    atan2 x y        = catan2 x y
    scaleFloat       = scaleFloat64
  
  decodeFloat64 :: Double -> (Integer, Int)
  decodeFloat64 x =
    let xw   = primWord64FromDoubleRaw x
        sign =  xw .&. 0x8000000000000000
        expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
        mant =  xw .&. 0x000fffffffffffff
        neg  = if sign /= 0 then negate else id
        -- Haskell promises that the mantissa is normalized
        adjDenorm m e | abs m < (1 `shiftL` 52) = adjDenorm (2 * m) (e - 1)
                      | otherwise = (m, e)
    in  if expn == 0 then
          -- denormalized or 0
          if mant == 0 then
            (0, 0)  -- 0
          else
            adjDenorm (neg (_word64ToInteger mant)) (-1074)  -- denormalized
        else if expn == 0x7ff then
          -- infinity or NaN
          --(0, 0)
          error "decodeFloat64: inf or NaN"
        else
          -- ordinary number, add hidden bit
          -- mant is offset-1023, and assumes scaled mantissa (thus -52)
          (neg (_word64ToInteger (mant .|. 0x0010000000000000)),
           primWordToInt (primWord64ToWord expn) - 1023 - 52)
  
  isNaNFloat64 :: Double -> Bool
  isNaNFloat64 x =
    let xw   = primWord64FromDoubleRaw x
        expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
        mant =  xw .&. 0x000fffffffffffff
    in  expn == 0x7ff && mant /= 0
  
  isInfFloat64 :: Double -> Bool
  isInfFloat64 x =
    let xw   = primWord64FromDoubleRaw x
        expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
        mant =  xw .&. 0x000fffffffffffff
    in  expn == 0x7ff && mant == 0
  
  isDenFloat64 :: Double -> Bool
  isDenFloat64 x =
    let xw   = primWord64FromDoubleRaw x
        expn = (xw .&. 0x7fffffffffffffff) `shiftR` 52
        mant =  xw .&. 0x000fffffffffffff
    in  expn == 0 && mant /= 0
  
  isNegZeroFloat64 :: Double -> Bool
  isNegZeroFloat64 x =
    primWord64FromDoubleRaw x == 0x8000000000000000
  
  scaleFloat64 :: Int -> Double -> Double
  scaleFloat64 n x = cscalbn x n
  
  encodeFloat64 :: Integer -> Int -> Double
  encodeFloat64 mant expn = scaleFloat64 expn (fromInteger mant)
  