content: {-# OPTIONS_JHC -fno-prelude -fffi -fm4  -funboxed-tuples #-}
  
  module Prelude.Float(readDouble,doubleToDigits,doubleToRational) where
  
  import Jhc.Basics
  import Jhc.Class.Real
  import Jhc.Float
  import Jhc.List(length,notElem,take,elem)
  import Jhc.Num
  import Jhc.Numeric((^),(^^))
  import Jhc.Order
  import Jhc.Type.C
  import Numeric
  import Prelude.Text
  
  m4_define(INST,{{
  
  foreign import primitive "FDiv" divide$2 ::  $2 -> $2 -> $2
  foreign import primitive "FPwr" exponent$2 ::  $2 -> $2 -> $2
  foreign import primitive "FAtan2" atan2$1 ::  $1 -> $1 -> $1
  foreign import primitive "F2I"  toInteger$1 :: $1 -> Integer
  foreign import primitive "const.M_PI" c_pi$1 :: $1
  
  instance Fractional $1 where
      $1 x / $1 y = $1 (divide$2 x y)
      fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
  
  m4_define(FI,foreign import primitive "{{$}}1"  {{$}}2$2 :: $2 -> $2)
  
  FI(Sqrt,sqrt)
  FI(Exp,exp)
  FI(Log,log)
  FI(Sin,sin)
  FI(Cos,cos)
  FI(Tan,tan)
  FI(Sinh,sinh)
  FI(Cosh,cosh)
  FI(Tanh,tanh)
  FI(Asin,asin)
  FI(Acos,acos)
  FI(Atan,atan)
  
  m4_undefine({{FI}})
  
  instance Floating $1 where
      pi = c_pi$1
      sqrt ($1 x) = $1 (sqrt$2 x)
      exp ($1 x) = $1 (exp$2 x)
      log ($1 x) = $1 (log$2 x)
      sin ($1 x) = $1 (sin$2 x)
      cos ($1 x) = $1 (cos$2 x)
      tan ($1 x) = $1 (tan$2 x)
      asin ($1 x) = $1 (asin$2 x)
      acos ($1 x) = $1 (acos$2 x)
      atan ($1 x) = $1 (atan$2 x)
      sinh ($1 x) = $1 (sinh$2 x)
      cosh ($1 x) = $1 (cosh$2 x)
      tanh ($1 x) = $1 (tanh$2 x)
      $1 x ** $1 y = $1 (exponent$2 x y)
  
      asinh = c_asinh$1
      acosh = c_acosh$1
      atanh = c_atanh$1
  
  instance RealFrac $1 where
      properFraction x
        = case (decodeFloat x)      of { (m,n) ->
      	let  b = floatRadix x     in
      	if n >= 0 then
  	    (fromInteger m * fromInteger b ^ n, 0.0)
      	else
  	    case (quotRem m (b^(negate n))) of { (w,r) ->
  	    (fromInteger w, encodeFloat r n)
  	    }
          }
  
      truncate x = fromInteger (toInteger$1 x)
      round x = fromInteger (toInteger$1 (roundf x))
      ceiling x = fromInteger (toInteger$1 (ceilingf x))
      floor x = fromInteger (toInteger$1 (floorf x))
  
      properFractionf x = (c_trunc$1 x,x - c_trunc$1 x)
      truncatef x = c_trunc$1 x
      roundf x = c_nearbyint$1 x
      ceilingf x = c_ceil$1 x
      floorf x = c_floor$1 x
  
  foreign import ccall "-lm math.h asinh$3" c_asinh$1 :: $1 -> $1
  foreign import ccall "-lm math.h acosh$3" c_acosh$1 :: $1 -> $1
  foreign import ccall "-lm math.h atanh$3" c_atanh$1 :: $1 -> $1
  foreign import ccall "-lm math.h trunc$3" c_trunc$1 :: $1 -> $1
  foreign import ccall "-lm math.h ceil$3" c_ceil$1 :: $1 -> $1
  foreign import ccall "-lm math.h floor$3" c_floor$1 :: $1 -> $1
  foreign import ccall "-lm math.h nearbyint$3" c_nearbyint$1 :: $1 -> $1
  
  foreign import ccall "math.h isnan" c_isnan$3 :: $1 -> CInt
  foreign import ccall "math.h isinf" c_isinfinite$3 :: $1 -> CInt
  foreign import ccall "math.h signbit" c_signbit$3 :: $1 -> CInt
  
  foreign import ccall "math.h ldexp$3"  c_ldexp$3 :: $1 -> CInt -> $1
  foreign import ccall "math.h frexp$3"  c_frexp$3 :: $1 -> (# $1, CInt #)
  
  }})
  
  INST(Float,Float32_,f)
  INST(Double,Float64_)
  
  instance Real Float where
      toRational x	=  (m:%1)*(b:%1)^^n
  			   where (m,n) = decodeFloat x
  				 b     = floatRadix  x
      toDouble x = floatToDouble x
  
  instance Real Double where
      toRational x = doubleToRational x
      toDouble x = x
  
  instance RealFloat Float where
      floatRadix _ = 2
      floatDigits _ = 24
      floatRange _ = (-125,128)
  
      exponent x		= case decodeFloatf x of (_,n) -> n
      significand x	= case decodeFloatf x of (m,_) -> m
  
      isNaN x = c_isnanf x /= 0
      isInfinite x = c_isinfinitef x /= 0
      isDenormalized _ = False
      isNegativeZero x = x == 0 && c_signbitf x /= 0
      isIEEE _ = True
  
      scaleFloat k x = c_ldexpf x (fromInt k)
      decodeFloatf x = case c_frexpf x of
          (# x', exp #) -> (x', fromIntegral exp)
  
      encodeFloat i e = c_ldexpf (fromInteger i) (fromInt e)
      decodeFloat x = case c_frexpf x of
          (# x', exp #) -> let x'' =  c_ldexp (floatToDouble x') (fromInt $ floatDigits x) in
              (double2integer x'', fromIntegral exp  - floatDigits x)
  
      atan2 = atan2Float
  
  instance RealFloat Double where
      floatRadix _ = 2
      floatDigits _ = 53
      floatRange _ = (-1021,1024)
  
      exponent x		= case decodeFloatf x of (_,n) -> n
      significand x	= case decodeFloatf x of (m,_) -> m
  
      isNaN x = c_isnan x /= 0
      isInfinite x = c_isinfinite x /= 0
      isDenormalized _ = False
      isNegativeZero x = x == 0 && c_signbit x /= 0
      isIEEE _ = True
      scaleFloat k x = c_ldexp x (fromInt k)
      decodeFloatf x = case c_frexp x of
          (# x', exp #) -> (x', fromIntegral exp)
  
      encodeFloat i e =  c_ldexp (integer2double i) (fromInt e)
      decodeFloat x = case c_frexp x of
          (# x', exp #) -> let x'' = c_ldexp x' (fromInt $ floatDigits x) in
              (double2integer x'', fromIntegral exp  - floatDigits x)
  
      atan2 = atan2Double
  
  foreign import primitive "I2F" integer2float :: Integer -> Float
  foreign import primitive "I2F" integer2double :: Integer -> Double
  foreign import primitive "F2I" double2integer :: Double -> Integer
  
  readDouble :: ReadS Double
  readDouble r    = [((fromInteger n * (10^^(k-d))),t) | (n,d,s) <- readFix r,(k,t)   <- readExp s] ++
                   [ (0/0, t) | ("NaN",t)      <- lex r] ++
                   [ (1/0, t) | ("Infinity",t) <- lex r]
                 where
                   readFix r = [(read (take 15 $ ds++ds'), length ds', t)
                               | (ds,d) <- lexDigits r,
                                 (ds',t) <- lexFrac d ]
  
                   lexFrac ('.':ds) = lexDigits ds
                   lexFrac s        = [("",s)]
  
                   readExp (e:s) | e `elem` "eE" = readExp' s
                   readExp s                     = [(0,s)]
  
                   readExp' ('-':s) = [(-k,t) | (k,t) <- readDec s]
                   readExp' ('+':s) = readDec s
                   readExp' s       = readDec s
  
  doubleToDigits :: Integer -> Double -> ([Int], Int)
  doubleToDigits n d | n `seq` d `seq` d == 0 = ([], 0)
  doubleToDigits base' x =
      let (f0', e0) = decodeFloat x
          base, f0 :: WordMax
          base = fromInteger base'
          f0 = fromIntegral f0'
          (minExp0, _) = floatRange x
          p = floatDigits x
          b :: WordMax
          b = fromInteger $ 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 :: WordMax
          e :: Int
          (f, e) = let n = minExp - e0
                   in  if n > 0 then (f0 `div` (b^n), e0+n) else (f0, e0)
  
          (r, s, mUp, mDn) =
             if e >= 0 then
                 let be = b^e in
                 if f == b^(p-1) then
                     (f*be*b*2, 2*b, be*b, b)
                 else
                     (f*be*2, 2, be, be)
             else
                 if e > minExp && f == b^(p-1) then
                     (f*b*2, b^(-e+1)*2, b, 1)
                 else
                     (f*2, b^(-e)*2, 1, 1)
          k =
              let k0 =
                      if b==2 && base==10 then
                          -- logBase 10 2 is slightly bigger than 3/10 so
                          -- the following will err on the low side.  Ignoring
                          -- the fraction will make it err even more.
                          -- Haskell promises that p-1 <= logBase b f < p.
                          (p - 1 + e0) * 3 `div` 10
                      else
                          ceiling ((log ((fromIntegral (f+1))::Double) +
                                   fromIntegral e * log (fromIntegral b)) /
                                    log (fromIntegral base))
                  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::Int)
  
          gen ds rn sN mUpN mDnN | rn `seq` sN `seq` mUpN `seq` mDnN `seq` True =
              let (dn, rn') = (rn * base) `divMod` sN
                  mUpN' = mUpN * base
                  mDnN' = mDnN * base
              in  case (rn' < mDnN', rn' + mUpN' > sN) of
                  (True,  False) -> toInt dn : ds
                  (False, True)  -> toInt (dn+1) : ds
                  (True,  True)  -> if rn' * 2 < sN then toInt dn : ds else toInt (dn+1) : ds
                  (False, False) -> gen (toInt dn:ds) rn' sN mUpN' mDnN'
          rds,rrds :: [Int]
          rrds = reverse rds
          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)
          expt :: WordMax -> Int -> WordMax
          expt base n = base^n
      in  k `seq` f `seq` e `seq` b `seq` rrds `seq` (rrds, k)
  
  doubleToRational :: Double -> Rational
  doubleToRational x  =  (m:%1)*(b:%1)^^n where
      (m,n) = decodeFloat x
      b     = floatRadix  x
  