{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} -- -fglasgow-exts for RULES module Algebra.RealField where import qualified Algebra.Field as Field import qualified Algebra.PrincipalIdealDomain as PID import qualified Algebra.Real as Real import qualified Algebra.Ring as Ring import qualified Algebra.ToRational as ToRational import qualified Algebra.ToInteger as ToInteger import Algebra.Field ((/), fromRational, ) import Algebra.RealIntegral (quotRem, ) import Algebra.IntegralDomain (divMod, even, ) import Algebra.Ring ((*), fromInteger, ) import Algebra.Additive ((+), (-), negate, ) import Algebra.ZeroTestable (isZero, ) import Algebra.ToInteger (fromIntegral, ) import qualified Number.Ratio as Ratio import Number.Ratio (T((:%)), Rational) import Data.Int (Int, Int8, Int16, Int32, Int64, ) import Data.Word (Word, Word8, Word16, Word32, Word64, ) import qualified GHC.Float as GHC import NumericPrelude.Tuple (mapFst, ) import Prelude(Int, Integer, Float, Double) import qualified Prelude as P import PreludeBase {- | Minimal complete definition: 'splitFraction' or 'floor' There are probably more laws, but some laws are > (fromInteger.fst.splitFraction) a + (snd.splitFraction) a === a > ceiling (toRational x) === ceiling x :: Integer > truncate (toRational x) === truncate x :: Integer > floor (toRational x) === floor x :: Integer If there wouldn't be @Real.C a@ and @ToInteger.C b@ constraints, we could also use this class for splitting ratios of polynomials. As an aside, let me note the similarities between @splitFraction x@ and @x divMod 1@ (if that were defined). In particular, it might make sense to unify the rounding modes somehow. IEEEFloat-specific calls are removed here (cf. 'Prelude.RealFloat') so probably nobody will actually use this default definition. Henning: New function 'fraction' doesn't return the integer part of the number. This also removes a type ambiguity if the integer part is not needed. The new methods 'fraction' and 'splitFraction' differ from 'Prelude.properFraction' semantics. They always round to 'floor'. This means that the fraction is always non-negative and is always smaller than 1. This is more useful in practice and can be generalised to more than real numbers. Since every 'Number.Ratio.T' denominator type supports 'Algebra.IntegralDomain.divMod', every 'Number.Ratio.T' can provide 'fraction' and 'splitFraction', e.g. fractions of polynomials. However the ''integral'' part would not be of type class 'ToInteger.C'. Can there be a separate class for 'fraction', 'splitFraction', 'floor' and 'ceiling' since they do not need reals and their ordering? -} class (Real.C a, Field.C a) => C a where splitFraction :: (ToInteger.C b) => a -> (b,a) fraction :: a -> a ceiling, floor :: (ToInteger.C b) => a -> b truncate, round :: (ToInteger.C b) => a -> b splitFraction x = (floor x, fraction x) fraction x = x - fromInteger (floor x) floor x = fromInteger (fst (splitFraction x)) ceiling x = - floor (-x) -- truncate x = signum x * floor (abs x) truncate x = if x>=0 then floor x else ceiling x round x = let (n,r) = splitFraction x in case compare r (1/2) of LT -> n EQ -> if even n then n else n+1 GT -> n+1 instance (ToInteger.C a, PID.C a) => C (Ratio.T a) where splitFraction (x:%y) = (fromIntegral q, r:%y) where (q,r) = divMod x y instance C Float where {-# INLINE splitFraction #-} {-# INLINE fraction #-} {-# INLINE floor #-} {-# INLINE ceiling #-} {-# INLINE round #-} {-# INLINE truncate #-} splitFraction = fastSplitFraction GHC.float2Int GHC.int2Float fraction = fastFraction (GHC.int2Float . GHC.float2Int) floor = fromInteger . P.floor ceiling = fromInteger . P.ceiling round = fromInteger . P.round truncate = fromInteger . P.truncate instance C Double where {-# INLINE splitFraction #-} {-# INLINE fraction #-} {-# INLINE floor #-} {-# INLINE ceiling #-} {-# INLINE round #-} {-# INLINE truncate #-} splitFraction = fastSplitFraction GHC.double2Int GHC.int2Double fraction = fastFraction (GHC.int2Double . GHC.double2Int) floor = fromInteger . P.floor ceiling = fromInteger . P.ceiling round = fromInteger . P.round truncate = fromInteger . P.truncate {-# INLINE fastSplitFraction #-} fastSplitFraction :: (P.RealFrac a, Real.C a, ToInteger.C b) => (a -> Int) -> (Int -> a) -> a -> (b,a) fastSplitFraction trunc toFloat x = fixSplitFraction $ if toFloat minBound <= x && x <= toFloat maxBound then case trunc x of n -> (fromIntegral n, x - toFloat n) else case P.properFraction x of (n,f) -> (fromInteger n, f) {-# fixSplitFraction #-} fixSplitFraction :: (Ring.C a, Ring.C b, Ord a) => (b,a) -> (b,a) fixSplitFraction (n,f) = -- if x>=0 || f==0 if f>=0 then (n, f) else (n-1, f+1) {-# INLINE fastFraction #-} fastFraction :: (P.RealFrac a, Real.C a) => (a -> a) -> a -> a fastFraction trunc x = fixFraction $ if fromIntegral (minBound :: Int) <= x && x <= fromIntegral (maxBound :: Int) then x - trunc x else preludeFraction x {-# INLINE preludeFraction #-} preludeFraction :: (P.RealFrac a, Ring.C a) => a -> a preludeFraction x = let second :: (Integer, a) -> a second = snd in second (P.properFraction x) {-# INLINE fixFraction #-} fixFraction :: (Ring.C a, Ord a) => a -> a fixFraction y = if y>=0 then y else y+1 {- mapM_ (\n -> let x = fromInteger n / 10 in print (x, floorInt GHC.double2Int GHC.int2Double x)) [-20,-19..20] -} {-# INLINE splitFractionInt #-} splitFractionInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> (Int, a) splitFractionInt trunc toFloat x = let n = trunc x in fixSplitFraction (n, x - toFloat n) {-# INLINE floorInt #-} floorInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int floorInt trunc toFloat x = let n = trunc x in if x >= toFloat n then n else pred n {-# INLINE ceilingInt #-} ceilingInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int ceilingInt trunc toFloat x = let n = trunc x in if x <= toFloat n then n else succ n {-# INLINE roundInt #-} roundInt :: (Field.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int roundInt trunc toFloat x = let half = 0.5 -- P.fromRational halfUp = x+half n = floorInt trunc toFloat halfUp in if toFloat n == halfUp && P.odd n then pred n else n {-# RULEZ maybe used, when Prelude implementations become more efficient "NP.round :: Float -> Int" round = P.round :: Float -> Int; "NP.truncate :: Float -> Int" truncate = P.truncate :: Float -> Int; "NP.floor :: Float -> Int" floor = P.floor :: Float -> Int; "NP.ceiling :: Float -> Int" ceiling = P.ceiling :: Float -> Int; "NP.round :: Double -> Int" round = P.round :: Double -> Int; "NP.truncate :: Double -> Int" truncate = P.truncate :: Double -> Int; "NP.floor :: Double -> Int" floor = P.floor :: Double -> Int; "NP.ceiling :: Double -> Int" ceiling = P.ceiling :: Double -> Int; #-} -- these rules will also be needed for Int16 et.al. {-# RULES "NP.round :: Float -> Int" round = roundInt GHC.float2Int GHC.int2Float; "NP.truncate :: Float -> Int" truncate = GHC.float2Int ; "NP.floor :: Float -> Int" floor = floorInt GHC.float2Int GHC.int2Float; "NP.ceiling :: Float -> Int" ceiling = ceilingInt GHC.float2Int GHC.int2Float; "NP.round :: Double -> Int" round = roundInt GHC.double2Int GHC.int2Double; "NP.truncate :: Double -> Int" truncate = GHC.double2Int ; "NP.floor :: Double -> Int" floor = floorInt GHC.double2Int GHC.int2Double; "NP.ceiling :: Double -> Int" ceiling = ceilingInt GHC.double2Int GHC.int2Double; "NP.splitFraction :: Float -> (Int, Float)" splitFraction = splitFractionInt GHC.float2Int GHC.int2Float; "NP.splitFraction :: Double -> (Int, Double)" splitFraction = splitFractionInt GHC.double2Int GHC.int2Double; #-} -- generated by GenerateRules.hs {-# RULES "NP.round :: a -> Int8" round = (P.fromIntegral :: Int -> Int8) . round; "NP.truncate :: a -> Int8" truncate = (P.fromIntegral :: Int -> Int8) . truncate; "NP.floor :: a -> Int8" floor = (P.fromIntegral :: Int -> Int8) . floor; "NP.ceiling :: a -> Int8" ceiling = (P.fromIntegral :: Int -> Int8) . ceiling; "NP.round :: a -> Int16" round = (P.fromIntegral :: Int -> Int16) . round; "NP.truncate :: a -> Int16" truncate = (P.fromIntegral :: Int -> Int16) . truncate; "NP.floor :: a -> Int16" floor = (P.fromIntegral :: Int -> Int16) . floor; "NP.ceiling :: a -> Int16" ceiling = (P.fromIntegral :: Int -> Int16) . ceiling; "NP.round :: a -> Int32" round = (P.fromIntegral :: Int -> Int32) . round; "NP.truncate :: a -> Int32" truncate = (P.fromIntegral :: Int -> Int32) . truncate; "NP.floor :: a -> Int32" floor = (P.fromIntegral :: Int -> Int32) . floor; "NP.ceiling :: a -> Int32" ceiling = (P.fromIntegral :: Int -> Int32) . ceiling; "NP.round :: a -> Int64" round = (P.fromIntegral :: Int -> Int64) . round; "NP.truncate :: a -> Int64" truncate = (P.fromIntegral :: Int -> Int64) . truncate; "NP.floor :: a -> Int64" floor = (P.fromIntegral :: Int -> Int64) . floor; "NP.ceiling :: a -> Int64" ceiling = (P.fromIntegral :: Int -> Int64) . ceiling; "NP.round :: a -> Word" round = (P.fromIntegral :: Int -> Word) . round; "NP.truncate :: a -> Word" truncate = (P.fromIntegral :: Int -> Word) . truncate; "NP.floor :: a -> Word" floor = (P.fromIntegral :: Int -> Word) . floor; "NP.ceiling :: a -> Word" ceiling = (P.fromIntegral :: Int -> Word) . ceiling; "NP.round :: a -> Word8" round = (P.fromIntegral :: Int -> Word8) . round; "NP.truncate :: a -> Word8" truncate = (P.fromIntegral :: Int -> Word8) . truncate; "NP.floor :: a -> Word8" floor = (P.fromIntegral :: Int -> Word8) . floor; "NP.ceiling :: a -> Word8" ceiling = (P.fromIntegral :: Int -> Word8) . ceiling; "NP.round :: a -> Word16" round = (P.fromIntegral :: Int -> Word16) . round; "NP.truncate :: a -> Word16" truncate = (P.fromIntegral :: Int -> Word16) . truncate; "NP.floor :: a -> Word16" floor = (P.fromIntegral :: Int -> Word16) . floor; "NP.ceiling :: a -> Word16" ceiling = (P.fromIntegral :: Int -> Word16) . ceiling; "NP.round :: a -> Word32" round = (P.fromIntegral :: Int -> Word32) . round; "NP.truncate :: a -> Word32" truncate = (P.fromIntegral :: Int -> Word32) . truncate; "NP.floor :: a -> Word32" floor = (P.fromIntegral :: Int -> Word32) . floor; "NP.ceiling :: a -> Word32" ceiling = (P.fromIntegral :: Int -> Word32) . ceiling; "NP.round :: a -> Word64" round = (P.fromIntegral :: Int -> Word64) . round; "NP.truncate :: a -> Word64" truncate = (P.fromIntegral :: Int -> Word64) . truncate; "NP.floor :: a -> Word64" floor = (P.fromIntegral :: Int -> Word64) . floor; "NP.ceiling :: a -> Word64" ceiling = (P.fromIntegral :: Int -> Word64) . ceiling; "NP.splitFraction :: a -> (Int8,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Int8) . splitFraction; "NP.splitFraction :: a -> (Int16,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Int16) . splitFraction; "NP.splitFraction :: a -> (Int32,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Int32) . splitFraction; "NP.splitFraction :: a -> (Int64,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Int64) . splitFraction; "NP.splitFraction :: a -> (Word,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Word) . splitFraction; "NP.splitFraction :: a -> (Word8,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Word8) . splitFraction; "NP.splitFraction :: a -> (Word16,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Word16) . splitFraction; "NP.splitFraction :: a -> (Word32,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Word32) . splitFraction; "NP.splitFraction :: a -> (Word64,a)" splitFraction = mapFst (P.fromIntegral :: Int -> Word64) . splitFraction; #-} {- | TODO: Should be moved to a continued fraction module. -} approxRational :: (ToRational.C a, C a) => a -> a -> Rational approxRational rat eps = simplest (rat-eps) (rat+eps) where simplest x y | y < x = simplest y x | x == y = xr | x > 0 = simplest' n d n' d' | y < 0 = - simplest' (-n') d' (-n) d | otherwise = 0 :% 1 where xr@(n:%d) = ToRational.toRational x (n':%d') = ToRational.toRational y simplest' n d n' d' -- assumes 0 < n%d < n'%d' | isZero r = q :% 1 | q /= q' = (q+1) :% 1 | otherwise = (q*n''+d'') :% n'' where (q,r) = quotRem n d (q',r') = quotRem n' d' (n'':%d'') = simplest' d' r' d r