{-# LANGUAGE NoImplicitPrelude #-}
module Algebra.RealRing where
import qualified Algebra.RealRing98 as RealRing98
import qualified Algebra.PrincipalIdealDomain as PID
import qualified Algebra.Field as Field
import qualified Algebra.Ring as Ring
import qualified Algebra.ToRational as ToRational
import qualified Algebra.ToInteger as ToInteger
import qualified Algebra.Absolute as Absolute
import qualified Algebra.OrderDecision as OrdDec
import Algebra.OrderDecision ((<?), (>=?), )
import Algebra.Field (fromRational, )
import Algebra.RealIntegral (quotRem, )
import Algebra.IntegralDomain (divMod, even, )
import Algebra.Ring ((*), fromInteger, one, )
import Algebra.Additive ((+), (-), negate, zero, )
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 Data.List as List
import Data.Tuple.HT (mapFst, mapPair, )
import Prelude (Integer, Float, Double, )
import qualified Prelude as P
import NumericPrelude.Base
class (Absolute.C a, Ord a) => C a where
{-# MINIMAL splitFraction | floor #-}
splitFraction :: (Ring.C b) => a -> (b,a)
fraction :: a -> a
ceiling, floor :: (Ring.C b) => a -> b
truncate :: (Ring.C b) => a -> b
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 =
if x>=0
then floor x
else ceiling x
round x =
let (n,r) = splitFraction x
in case compare (2*r) one of
LT -> n
EQ -> if even n then n else n+1
GT -> n+1
roundSimple :: (C a, Ring.C b) => a -> b
roundSimple x =
let (n,r) = splitFraction x
in case compare (2*r) one of
LT -> n
EQ -> if x<0 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 Integer where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromInteger x, zero)
fraction _ = zero
floor x = fromInteger x
ceiling x = fromInteger x
round x = fromInteger x
truncate x = fromInteger x
instance C Int where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Int8 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Int16 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Int32 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Int64 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Word8 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Word16 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Word32 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Word64 where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction x = (fromIntegral x, zero)
fraction _ = zero
floor x = fromIntegral x
ceiling x = fromIntegral x
round x = fromIntegral x
truncate x = fromIntegral x
instance C Float where
{-# INLINE splitFraction #-}
{-# INLINE fraction #-}
{-# INLINE floor #-}
{-# INLINE ceiling #-}
{-# INLINE round #-}
{-# INLINE truncate #-}
splitFraction = fastSplitFraction GHC.float2Int GHC.int2Float
fraction = RealRing98.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 = RealRing98.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, Absolute.C a, Ring.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)
{-# INLINE fixSplitFraction #-}
fixSplitFraction :: (Ring.C a, Ring.C b, Ord a) => (b,a) -> (b,a)
fixSplitFraction (n,f) =
if f>=0
then (n, f)
else (n-1, f+1)
{-# INLINE fixFraction #-}
fixFraction :: (Ring.C a, Ord a) => a -> a
fixFraction y =
if y>=0 then y else y+1
{-# 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
halfUp = x+half
n = floorInt trunc toFloat halfUp
in if toFloat n == halfUp && P.odd n
then pred n
else n
{-# INLINE roundSimpleInt #-}
roundSimpleInt ::
(Field.C a, Absolute.C a, Ord a) =>
(a -> Int) -> (Int -> a) -> a -> Int
roundSimpleInt trunc _toFloat x =
trunc (x + Absolute.signum x * 0.5)
{-# RULES
"NP.round :: Float -> Int" round = roundInt GHC.float2Int GHC.int2Float;
"NP.roundSimple :: Float -> Int" round = roundSimpleInt 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.roundSimple :: Double -> Int" round = roundSimpleInt 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;
#-}
{-# RULES
"NP.round :: a -> Int8" round = (P.fromIntegral :: Int -> Int8) . round;
"NP.roundSimple :: a -> Int8" roundSimple = (P.fromIntegral :: Int -> Int8) . roundSimple;
"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.roundSimple :: a -> Int16" roundSimple = (P.fromIntegral :: Int -> Int16) . roundSimple;
"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.roundSimple :: a -> Int32" roundSimple = (P.fromIntegral :: Int -> Int32) . roundSimple;
"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.roundSimple :: a -> Int64" roundSimple = (P.fromIntegral :: Int -> Int64) . roundSimple;
"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.roundSimple :: a -> Word" roundSimple = (P.fromIntegral :: Int -> Word) . roundSimple;
"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.roundSimple :: a -> Word8" roundSimple = (P.fromIntegral :: Int -> Word8) . roundSimple;
"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.roundSimple :: a -> Word16" roundSimple = (P.fromIntegral :: Int -> Word16) . roundSimple;
"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.roundSimple :: a -> Word32" roundSimple = (P.fromIntegral :: Int -> Word32) . roundSimple;
"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.roundSimple :: a -> Word64" roundSimple = (P.fromIntegral :: Int -> Word64) . roundSimple;
"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;
#-}
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'
| 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
powersOfTwo :: (Ring.C a) => [a]
powersOfTwo = iterate (2*) one
pairsOfPowersOfTwo :: (Ring.C a, Ring.C b) => [(a,b)]
pairsOfPowersOfTwo =
zip powersOfTwo powersOfTwo
genericFloor :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericFloor a =
if a>=zero
then genericPosFloor a
else negate $ genericPosCeiling $ negate a
genericCeiling :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericCeiling a =
if a>=zero
then genericPosCeiling a
else negate $ genericPosFloor $ negate a
genericTruncate :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericTruncate a =
if a>=zero
then genericPosFloor a
else negate $ genericPosFloor $ negate a
genericRound :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericRound a =
if a>=zero
then genericPosRound a
else negate $ genericPosRound $ negate a
genericFraction :: (Ord a, Ring.C a) => a -> a
genericFraction a =
if a>=zero
then genericPosFraction a
else fixFraction $ negate $ genericPosFraction $ negate a
genericSplitFraction :: (Ord a, Ring.C a, Ring.C b) => a -> (b,a)
genericSplitFraction a =
if a>=zero
then genericPosSplitFraction a
else fixSplitFraction $ mapPair (negate, negate) $
genericPosSplitFraction $ negate a
genericPosFloor :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericPosFloor a =
snd $
foldr
(\(pa,pb) acc@(accA,accB) ->
let newA = accA+pa
in if newA>a then acc else (newA,accB+pb))
(zero,zero) $
takeWhile ((a>=) . fst) $
pairsOfPowersOfTwo
genericPosCeiling :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericPosCeiling a =
snd $
(\(ps,u:_) ->
foldr
(\(pa,pb) acc@(accA,accB) ->
let newA = accA-pa
in if newA>=a then (newA,accB-pb) else acc)
u ps) $
span ((a>) . fst) $
(zero,zero) : pairsOfPowersOfTwo
genericHalfPosFloorDigits :: (Ord a, Ring.C a, Ring.C b) => a -> ((a,b), [Bool])
genericHalfPosFloorDigits a =
List.mapAccumR
(\acc@(accA,accB) (pa,pb) ->
let newA = accA+pa
b = newA<=a
in (if b then (newA,accB+pb) else acc, b))
(zero,zero) $
takeWhile ((a>=) . fst) $
zip powersOfTwo (zero:powersOfTwo)
genericPosRound :: (Ord a, Ring.C a, Ring.C b) => a -> b
genericPosRound a =
let a2 = 2*a
((ai,bi), ds) = genericHalfPosFloorDigits a2
in if ai==a2
then
case ds of
True : True : _ -> bi+one
_ -> bi
else
case ds of
True : _ -> bi+one
_ -> bi
genericPosFraction :: (Ord a, Ring.C a) => a -> a
genericPosFraction a =
foldr
(\p acc ->
if p>acc then acc else acc-p)
a $
takeWhile (a>=) $
powersOfTwo
genericPosSplitFraction :: (Ord a, Ring.C a, Ring.C b) => a -> (b,a)
genericPosSplitFraction a =
foldr
(\(pb,pa) acc@(accB,accA) ->
if pa>accA then acc else (accB+pb,accA-pa))
(zero,a) $
takeWhile ((a>=) . snd) $
pairsOfPowersOfTwo
decisionPosFraction :: (OrdDec.C a, Ring.C a) => a -> a
decisionPosFraction a0 =
(\ps ->
foldr
(\p cont a ->
(a<?one) a $ cont $
(a>=?p) (a-p) a)
(error "decisionPosFraction: end of list should never be reached")
ps a0) $
concatMap (reverse . flip take powersOfTwo) powersOfTwo
decisionPosFractionSqrTime :: (OrdDec.C a, Ring.C a) => a -> a
decisionPosFractionSqrTime a0 =
(\ps ->
foldr
(\p cont a ->
(a<?one) a $ cont $
(a>=?p) (a-p) a)
(error "decisionPosFraction: end of list should never be reached")
ps a0) $
concatMap reverse $
inits powersOfTwo