diff --git a/docs/NOTES b/docs/NOTES --- a/docs/NOTES +++ b/docs/NOTES @@ -65,6 +65,11 @@ * people probably interested in NumPrelude: + Date: Wed, 28 Oct 2009 21:06:40 -0700 + From: Paul McJones + + Andrew nochwas aus Wien oder so ... + Mike Thomas http://www.haskell.org/pipermail/haskell-cafe/2002-February/002660.html diff --git a/numeric-prelude.cabal b/numeric-prelude.cabal --- a/numeric-prelude.cabal +++ b/numeric-prelude.cabal @@ -1,5 +1,5 @@ Name: numeric-prelude -Version: 0.1.2 +Version: 0.1.3 License: GPL License-File: LICENSE Author: Dylan Thurston , Henning Thielemann , Mikael Johansson diff --git a/src/Algebra/Monoid.hs b/src/Algebra/Monoid.hs --- a/src/Algebra/Monoid.hs +++ b/src/Algebra/Monoid.hs @@ -1,7 +1,7 @@ {-# LANGUAGE NoImplicitPrelude #-} {- | -Copyright : (c) Mikael Johansson 2006 -Maintainer : mik@math.uni-jena.de +Copyright : (c) Henning Thielemann 2009, Mikael Johansson 2006 +Maintainer : numericprelude@henning-thielemann.de Stability : provisional Portability : @@ -28,7 +28,7 @@ idt :: a (<*>) :: a -> a -> a -instance C Mn.All where +instance C All where idt = mempty (<*>) = mappend @@ -60,4 +60,3 @@ instance Additive.C a => C (Sum a) where idt = Mn.Sum Additive.zero (Mn.Sum x) <*> (Mn.Sum y) = Mn.Sum (x Additive.+ y) - diff --git a/src/Algebra/Real.hs b/src/Algebra/Real.hs --- a/src/Algebra/Real.hs +++ b/src/Algebra/Real.hs @@ -33,6 +33,11 @@ is not defined. To me, this seems to have the wrong type anyway; Complex.magnitude has the correct type. + +Note: +The abs function can be defined for Additive and Ord, +the Ring constraint is not needed. +We may move 'signum' to the new 'RealRing' class. -} class (Ring.C a, ZeroTestable.C a, Ord a) => C a where abs :: a -> a diff --git a/src/Algebra/RealField.hs b/src/Algebra/RealField.hs --- a/src/Algebra/RealField.hs +++ b/src/Algebra/RealField.hs @@ -14,8 +14,8 @@ import Algebra.Field ((/), fromRational, ) import Algebra.RealIntegral (quotRem, ) import Algebra.IntegralDomain (divMod, even, ) -import Algebra.Ring ((*), fromInteger, ) -import Algebra.Additive ((+), (-), negate, ) +import Algebra.Ring ((*), fromInteger, one, ) +import Algebra.Additive ((+), (-), negate, zero, ) import Algebra.ZeroTestable (isZero, ) import Algebra.ToInteger (fromIntegral, ) @@ -26,7 +26,8 @@ import Data.Word (Word, Word8, Word16, Word32, Word64, ) import qualified GHC.Float as GHC -import Data.Tuple.HT (mapFst, ) +import Data.List as List +import Data.Tuple.HT (mapFst, mapPair, ) import Prelude(Int, Integer, Float, Double) import qualified Prelude as P import PreludeBase @@ -72,11 +73,28 @@ Can there be a separate class for 'fraction', 'splitFraction', 'floor' and 'ceiling' since they do not need reals and their ordering? + +Note: +All of these methods can be defined +exclusively with functions from Ord and Ring. +We could write a power-of-two-algorithm +like the one for finding the number of digits of an Integer +in FixedPoint-fractions module. +This would even be reasonably efficient. +I think the module should be renamed to RealRing, +and the superclass constraint should be lifted from Field to Ring. + +We might also add a round method, +that rounds 0.5 always up or always down. +This is much more efficient in inner loops +and is acceptable or even preferable for many applications. + +The ToInteger constraint can be lifted to Ring. -} class (Real.C a, Field.C a) => C a where splitFraction :: (ToInteger.C b) => a -> (b,a) - fraction :: a -> a + fraction :: a -> a ceiling, floor :: (ToInteger.C b) => a -> b truncate, round :: (ToInteger.C b) => a -> b @@ -143,7 +161,7 @@ then case trunc x of n -> (fromIntegral n, x - toFloat n) else case P.properFraction x of (n,f) -> (fromInteger n, f) -{-# fixSplitFraction #-} +{-# INLINE fixSplitFraction #-} fixSplitFraction :: (Ring.C a, Ring.C b, Ord a) => (b,a) -> (b,a) fixSplitFraction (n,f) = -- if x>=0 || f==0 @@ -208,7 +226,7 @@ else n -{-# RULEZ maybe used, when Prelude implementations become more efficient +{- RULES 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; @@ -217,7 +235,7 @@ "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 @@ -304,4 +322,138 @@ where (q,r) = quotRem n d (q',r') = quotRem n' d' (n'':%d'') = simplest' d' r' d r + + +-- * generic implementation of round functions + +powersOfTwo :: (Ring.C a) => [a] +powersOfTwo = iterate (2*) one + +pairsOfPowersOfTwo :: (Ring.C a, Ring.C b) => [(a,b)] +pairsOfPowersOfTwo = + zip powersOfTwo powersOfTwo + +{- | +The generic rounding functions need a number of operations +proportional to the number of binary digits of the integer portion. +If operations like multiplication with two and comparison +need time proportional to the number of binary digits, +then the overall rounding requires quadratic time. +-} +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 + +{- +genericPosFloorDigits :: (Ord a, Ring.C a, Ring.C b) => a -> ((a,b), [Bool]) +genericPosFloorDigits 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) $ + 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 diff --git a/src/Number/Complex.hs b/src/Number/Complex.hs --- a/src/Number/Complex.hs +++ b/src/Number/Complex.hs @@ -158,13 +158,13 @@ (-:) x y = Cons x (-y) -- | The conjugate of a complex number. -{-# SPECULATE conjugate :: T Double -> T Double #-} +{- SPECIALISE conjugate :: T Double -> T Double -} {-# INLINE conjugate #-} conjugate :: (Additive.C a) => T a -> T a conjugate (Cons x y) = Cons x (-y) -- | Scale a complex number by a real number. -{-# SPECULATE scale :: Double -> T Double -> T Double #-} +{- SPECIALISE scale :: Double -> T Double -> T Double -} {-# INLINE scale #-} scale :: (Ring.C a) => a -> T a -> T a scale r (Cons x y) = Cons (r * x) (r * y) @@ -187,7 +187,7 @@ but oriented in the positive real direction, whereas @'signum' z@ has the phase of @z@, but unit magnitude. -} -{-# SPECULATE signum :: T Double -> T Double #-} +{- SPECIALISE signum :: T Double -> T Double -} {-# INLINE signum #-} signum :: (Algebraic.C a, NormedEuc.C a a, ZeroTestable.C a) => T a -> T a signum z = @@ -196,14 +196,14 @@ else scale (recip (NormedEuc.norm z)) z -- | Form a complex number from polar components of magnitude and phase. -{-# SPECULATE fromPolar :: Double -> Double -> T Double #-} +{- SPECIALISE fromPolar :: Double -> Double -> T Double -} {-# INLINE fromPolar #-} fromPolar :: (Trans.C a) => a -> a -> T a fromPolar r theta = scale r (cis theta) -- | @'cis' t@ is a complex value with magnitude @1@ -- and phase @t@ (modulo @2*'pi'@). -{-# SPECULATE cis :: Double -> T Double #-} +{- SPECIALISE cis :: Double -> T Double -} {-# INLINE cis #-} cis :: (Trans.C a) => a -> T a cis theta = Cons (cos theta) (sin theta) @@ -286,8 +286,8 @@ isZero (Cons x y) = isZero x && isZero y instance (Additive.C a) => Additive.C (T a) where - {-# SPECULATE instance Additive.C (T Float) #-} - {-# SPECULATE instance Additive.C (T Double) #-} + {- SPECIALISE instance Additive.C (T Float) -} + {- SPECIALISE instance Additive.C (T Double) -} {-# INLINE zero #-} {-# INLINE negate #-} {-# INLINE (+) #-} @@ -298,8 +298,8 @@ negate = Elem.run $ Elem.with Cons <*>.-$ real <*>.-$ imag instance (Ring.C a) => Ring.C (T a) where - {-# SPECULATE instance Ring.C (T Float) #-} - {-# SPECULATE instance Ring.C (T Double) #-} + {- SPECIALISE instance Ring.C (T Float) -} + {- SPECIALISE instance Ring.C (T Double) -} {-# INLINE one #-} one = Cons one zero {-# INLINE (*) #-} @@ -492,8 +492,8 @@ instance (Real.C a, RealTrans.C a, Power a) => Trans.C (T a) where - {-# SPECULATE instance Trans.C (T Float) #-} - {-# SPECULATE instance Trans.C (T Double) #-} + {- SPECIALISE instance Trans.C (T Float) -} + {- SPECIALISE instance Trans.C (T Double) -} {-# INLINE pi #-} pi = fromReal pi {-# INLINE exp #-} diff --git a/src/Number/NonNegativeChunky.hs b/src/Number/NonNegativeChunky.hs --- a/src/Number/NonNegativeChunky.hs +++ b/src/Number/NonNegativeChunky.hs @@ -16,17 +16,17 @@ which allows for co-recursive merges. -} module Number.NonNegativeChunky - (T, fromChunks, fromNumber, toNumber, fromChunky98, toChunky98, - normalize, isNull, isPositive) where + (T, fromChunks, toChunks, fromNumber, toNumber, fromChunky98, toChunky98, + minMaxDiff, normalize, isNull, isPositive) where import qualified Numeric.NonNegative.ChunkyPrivate as Chunky98 import qualified Numeric.NonNegative.Class as NonNeg98 +import qualified Algebra.NonNegative as NonNeg import qualified Algebra.Field as Field import qualified Algebra.Real as Real import qualified Algebra.Ring as Ring import qualified Algebra.Additive as Additive -import qualified Algebra.NonNegative as NonNeg import qualified Algebra.ToInteger as ToInteger import qualified Algebra.ToRational as ToRational import qualified Algebra.IntegralDomain as Integral @@ -34,6 +34,9 @@ import qualified Algebra.ZeroTestable as ZeroTestable import Algebra.ZeroTestable (isZero, ) +import qualified Algebra.Monoid as Monoid +import qualified Data.Monoid as Mn98 + import Control.Monad (liftM, liftM2, ) import Test.QuickCheck (Arbitrary(..)) @@ -127,12 +130,6 @@ else error ("Numeric.NonNegative.Chunky."++funcName++": negative number") -{- | -In @glue x y == (z,r,b)@ -@z@ represents @min x y@, -@r@ represents @max x y - min x y@, -and @x b@ or @x>y ==> not b@, for @x==y@ the value of b is arbitrary. --} glue :: (NonNeg.C a) => [a] -> [a] -> ([a], [a], Bool) glue [] ys = ([], ys, True) glue xs [] = ([], xs, False) @@ -144,6 +141,18 @@ EQ -> (x, glue xs ys) in (z:zs,rs,b) +{- | +In @minMaxDiff x y == (z,r,b)@ +@z@ represents @min x y@, +@r@ represents @max x y - min x y@, +and @x b@ or @x>y ==> not b@, + for @x==y@ the value of b is arbitrary. +-} +minMaxDiff :: (NonNeg.C a) => T a -> T a -> (T a, T a, Bool) +minMaxDiff (Cons xs) (Cons ys) = + let (zs, rs, b) = glue xs ys + in (Cons zs, Cons rs, b) + equalList :: (NonNeg.C a) => [a] -> [a] -> Bool equalList x y = let (_,r,_) = glue x y @@ -186,8 +195,8 @@ isZero = isNullZT instance (NonNeg.C a) => Additive.C (T a) where - zero = Cons [] - (+) = lift2 (++) + zero = Monoid.idt + (+) = (Monoid.<*>) x - y = let (_,d,b) = glue (toChunks x) (toChunks y) d' = fromChunks d @@ -264,3 +273,11 @@ instance (Field.C a, Eq a, Show a, NonNeg.C a) => P98.Fractional (T a) where fromRational = fromNumber . fromRational (/) = legacyInstance + +instance (NonNeg.C a) => Mn98.Monoid (T a) where + mempty = Monoid.idt + mappend = (Monoid.<*>) + +instance (NonNeg.C a) => Monoid.C (T a) where + idt = Cons [] + (<*>) = lift2 (++) diff --git a/test/Test/Run.hs b/test/Test/Run.hs --- a/test/Test/Run.hs +++ b/test/Test/Run.hs @@ -1,5 +1,6 @@ module Main where +import qualified Test.Algebra.RealRing as RealRing import qualified Test.MathObj.Gaussian.Polynomial as GaussPoly import qualified Test.MathObj.Gaussian.Variance as GaussVariance import qualified Test.MathObj.Gaussian.Bell as GaussBell @@ -14,6 +15,7 @@ main :: IO () main = do HUnitText.runTestTT (HUnit.TestList $ + RealRing.tests : GaussVariance.tests : GaussBell.tests : GaussPoly.tests :