numeric-prelude 0.0.2 → 0.0.4
raw patch · 34 files changed
+1792/−319 lines, 34 filesdep ~non-negative
Dependency ranges changed: non-negative
Files
- docs/NOTES +61/−26
- numeric-prelude.cabal +14/−2
- src/Algebra/Additive.hs +152/−17
- src/Algebra/DimensionTerm.hs +19/−0
- src/Algebra/Field.hs +16/−1
- src/Algebra/GenerateRules.hs +86/−0
- src/Algebra/IntegralDomain.hs +98/−3
- src/Algebra/Module.hs +11/−0
- src/Algebra/NonNegative.hs +6/−0
- src/Algebra/Real.hs +84/−4
- src/Algebra/RealField.hs +165/−20
- src/Algebra/RealIntegral.hs +96/−2
- src/Algebra/Ring.hs +111/−12
- src/Algebra/ToInteger.hs +68/−5
- src/Algebra/ToRational.hs +18/−3
- src/Algebra/Transcendental.hs +61/−1
- src/Algebra/VectorSpace.hs +7/−3
- src/Algebra/ZeroTestable.hs +20/−11
- src/MathObj/LaurentPolynomial.hs +2/−2
- src/MathObj/Polynomial.hs +78/−24
- src/MathObj/PowerSeries.hs +64/−7
- src/MathObj/RootSet.hs +2/−2
- src/Number/Complex.hs +184/−133
- src/Number/DimensionTerm.hs +18/−4
- src/Number/NonNegative.hs +4/−19
- src/Number/NonNegativeChunky.hs +266/−0
- src/Number/Positional/Check.hs +0/−7
- src/Number/Ratio.hs +40/−10
- src/NumericPrelude.hs +1/−1
- src/NumericPrelude/Condition.hs +4/−0
- src/NumericPrelude/List.hs +20/−0
- src/NumericPrelude/Monad.hs +1/−0
- src/NumericPrelude/Text.hs +2/−0
- src/NumericPrelude/Tuple.hs +13/−0
docs/NOTES view
@@ -1,3 +1,28 @@+** -> ^?++* non-negative++NonNegative could require ZeroTestable+because Eq provides (==) and Additive provides zero+and thus ZeroTestable can be implemented for all these instances.++NonNegative class on top of Additive is the wrong way around+since non-negative types like Peano+do not provide full subtraction.+There should be AdditiveMonoid and Additive on top of this,+and NonNegative should be a distinct sub-class of AdditiveMonoid.++Unfortunately we often have a type which provides also negative numbers+(Integer, Int, Float, Double)+without a counterpart type for non-negative numbers.+That's the reason why the wrapper Number.NonNegative exists.+But it is clearly the wrong way around.+There should be a multi-parameter type class with functional dependency+for the relation of a full additive/subtractive with its non-negative counterpart+in order to write algorithms that work with+say Card and Int, Peano and signed Peano.++ * Proper place of abs and signum After reflection, perhaps 'abs' and 'signum' should be names for canonically multiplying an element by a unit, and need not necessarily@@ -24,13 +49,14 @@ used; e.g., in the definition of 'x % y'. This factorization seems useful in somewhat surprising generality.-However, there are useful spaces where it's not defined; e.g.,-computable reals. (abs x is computable, but signum x is not-continuous so not computable.)+However, there are useful spaces where it's not defined,+e.g. computable reals.+@abs x@ is computable, but @signum x@ is not continuous so not computable. -* Names of floating point classes-The name 'Float' does seem to imply that the number can be represented-as an integer times an appropriate power of a base.+Also note, that @abs@ can be defined within @Additive@ and @Ord@+(namely with the methods @negate@ and @max@)+whereas @signum@ needs a @one@ or another unit+which requires at least @Ring@ capabilities. * GHC bugs -fno-implicit-prelude is happy to use locally defined 'fromInteger',@@ -241,7 +267,7 @@ In contrast to multi-parameter VectorSpace,-we cannot force that 'v a' is also a method of Additive.+we cannot force that 'v a' is also a member of Additive. We cannot restrict the vector element types by a class constraint, but the routines acting on Vector containers can have these restrictions.@@ -252,7 +278,12 @@ hornerScalar :: Ring a => [a] -> a -> a hornerVector :: (Ring a, Vector v) => [v a] -> a -> v a +Actually we could again setup a multi-parameter type class+ VectorSpace v a+where 'v' is not a type but a type constructor+much like in the MArray class. + Advantages: - scale :: (Complex Double) -> [Complex Double] -> [Complex Double] is possible@@ -294,8 +325,9 @@ * ToDo: Classes - - Hilbert space (scalar product)- - Affine space+ - HilbertSpace (scalar product)+ - AffineSpace (affine combination a*>x + (1-a)*>y) as superclass of module,+ contains ConvexSpace * ToDo: Types @@ -355,29 +387,32 @@ Module is named Algebra.Module since there might be many people who want to define some type named Module.- the prop_* routines in NumExtras could be rewritten as simplification rules for GHC,+ the properties in Algebra.Laws could be rewritten as simplification rules for GHC, though they should be disabled by default, because the rules doesn't always apply due to overflows and rounding. How can one handle errors in a computation?- say, vectors mismatch,- there is an overflow,- a sum of two physical values with different units fails- (I have already implemented modules for dealing with units) etc.- Making the operations undefined for these cases- is ok if the programmer has control over the operands.- But if the values are given by the user- the programmer might want to obtain something- from which he can build a user friendly error message,- say, "the values 1m and 2s can't be added:+ On the one hand there are errors+ that can be avoided by respecting restrictions for function arguments.+ E.g. vectors dimension mismatch, unit mismatch, index out of array, division by zero.+ These should be prevented by proper types+ and if that is not possible then by 'undefined' values.+ There is no need to recover from them+ or handle them otherwise because they can be avoided by proper call of the functions,+ and in several cases are fulfilled automatically.+ However there might be some support for detecting programming errors+ e.g. by reporting+ "the values 1m and 2s can't be added: expression 1m+2s, sub-expression of ..."+ instad of just "unit mismatch".++ On the other hand there are errors+ that cannot be avoided easily, e.g. overflow.+ Overflow can be considered as using the wrong type,+ e.g. Int instead of Integer.+ But in more complicated cases you should return a Maybe. Examples for implicit configuration residue classes: modulus matrix computation: matrix size positional numbers: base- fixed point numbers: position of the dot-----+ fixed point numbers: position of the dot, i.e. denominator
numeric-prelude.cabal view
@@ -1,5 +1,5 @@ Name: numeric-prelude-Version: 0.0.2+Version: 0.0.4 License: GPL License-File: LICENSE Author: Dylan Thurston <dpt@math.harvard.edu>, Henning Thielemann <numericprelude@henning-thielemann.de>, Mikael Johansson@@ -120,12 +120,18 @@ Makefile docs/NOTES docs/README+ src/Algebra/GenerateRules.hs Flag splitBase description: Choose the new smaller, split-up base package. +Flag buildTests+ description: Build test executables+ default: False+ Library- Build-Depends: parsec >= 1, HUnit >=1 && <2, QuickCheck >=1 && <2, non-negative >=0.0.1 && <0.1+ Build-Depends: parsec >= 1, HUnit >=1 && <2, QuickCheck >=1 && <2+ Build-Depends: non-negative>=0.0.2 && <0.1 If flag(splitBase) Build-Depends: base >= 2, array, containers, random Else@@ -190,6 +196,7 @@ Number.FixedPoint Number.FixedPoint.Check Number.NonNegative+ Number.NonNegativeChunky Number.PartiallyTranscendental Number.Peano Number.Positional@@ -214,11 +221,14 @@ NumericPrelude.List NumericPrelude.Monad NumericPrelude.Text+ NumericPrelude.Tuple PreludeBase Executable test Hs-Source-Dirs: src, test Main-Is: Test.hs+ If !flag(buildTests)+ Buildable: False Executable testsuite Hs-Source-Dirs: src, test@@ -230,3 +240,5 @@ Test.MathObj.Polynomial Test.MathObj.PowerSeries Main-Is: Test/Run.hs+ If !flag(buildTests)+ Buildable: False
src/Algebra/Additive.hs view
@@ -18,9 +18,12 @@ import qualified Algebra.Laws as Laws +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import qualified Data.Ratio as Ratio98 import qualified Prelude as P-import Prelude(fromInteger)+import Prelude(Int, Integer, Float, Double, fromInteger, ) import PreludeBase @@ -50,7 +53,9 @@ -- | inverse with respect to '+' negate :: a -> a + {-# INLINE negate #-} negate a = zero - a+ {-# INLINE (-) #-} a - b = a + negate b {- |@@ -84,36 +89,158 @@ {-* Instances for atomic types -} -instance C P.Integer where- (+) = (P.+)- zero = P.fromInteger 0- negate = P.negate+instance C Integer where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-) -instance C P.Int where- (+) = (P.+)- zero = P.fromInteger 0- negate = P.negate+instance C Float where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-) -instance C P.Float where- (+) = (P.+)- zero = P.fromInteger 0- negate = P.negate+instance C Double where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-) -instance C P.Double where- (+) = (P.+)- zero = P.fromInteger 0- negate = P.negate +instance C Int where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-) +instance C Int8 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Int16 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Int32 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Int64 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)+++instance C Word where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Word8 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Word16 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Word32 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++instance C Word64 where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-}+ zero = P.fromInteger 0+ negate = P.negate+ (+) = (P.+)+ (-) = (P.-)++++ {-* Instances for composed types -} instance (C v0, C v1) => C (v0, v1) where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-} zero = (zero, zero) (+) (x0,x1) (y0,y1) = ((+) x0 y0, (+) x1 y1) (-) (x0,x1) (y0,y1) = ((-) x0 y0, (-) x1 y1) negate (x0,x1) = (negate x0, negate x1) instance (C v0, C v1, C v2) => C (v0, v1, v2) where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-} zero = (zero, zero, zero) (+) (x0,x1,x2) (y0,y1,y2) = ((+) x0 y0, (+) x1 y1, (+) x2 y2) (-) (x0,x1,x2) (y0,y1,y2) = ((-) x0 y0, (-) x1 y1, (-) x2 y2)@@ -132,6 +259,10 @@ instance (C v) => C (b -> v) where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-} zero _ = zero (+) f g x = (+) (f x) (g x) (-) f g x = (-) (f x) (g x)@@ -154,6 +285,10 @@ -- legacy instance (P.Integral a) => C (Ratio98.Ratio a) where+ {-# INLINE zero #-}+ {-# INLINE negate #-}+ {-# INLINE (+) #-}+ {-# INLINE (-) #-} zero = 0 (+) = (P.+) (-) = (P.-)
src/Algebra/DimensionTerm.hs view
@@ -127,12 +127,27 @@ cancelRight _ = noValue invertRecip :: C u => Recip (Recip u) -> u invertRecip _ = noValue+doubleRecip :: C u => u -> Recip (Recip u)+doubleRecip _ = noValue recipScalar :: Recip Scalar -> Scalar recipScalar _ = noValue {- * Example dimensions -} +{- ** Scalar -}++{- |+This class allows defining instances that are exclusively for 'Scalar' dimension.+You won't want to define instances by yourself.+-}+class C dim => IsScalar dim where+ toScalar :: dim -> Scalar++instance IsScalar Scalar where+ toScalar = id++ {- ** Basis dimensions -} data Length = Length@@ -184,6 +199,10 @@ {- ** Derived dimensions -} type Frequency = Recip Time++frequency :: Frequency+frequency = noValue+ data Voltage = Voltage
src/Algebra/Field.hs view
@@ -75,9 +75,13 @@ fromRational' :: Rational -> a (^-) :: a -> Integer -> a + {-# INLINE recip #-} recip a = one / a+ {-# INLINE (/) #-} a / b = a * recip b+ {-# INLINE fromRational' #-} fromRational' r = fromInteger (numerator r) / fromInteger (denominator r)+ {-# INLINE (^-) #-} a ^- n = if n < zero then recip (a^(-n)) else a^n@@ -88,6 +92,7 @@ -- | Needed to work around shortcomings in GHC. +{-# INLINE fromRational #-} fromRational :: (C a) => P.Rational -> a fromRational x = fromRational' (Ratio98.numerator x :% Ratio98.denominator x) @@ -103,14 +108,21 @@ -} instance C Float where+ {-# INLINE (/) #-}+ {-# INLINE recip #-} (/) = (P./) recip = (P.recip) instance C Double where+ {-# INLINE (/) #-}+ {-# INLINE recip #-} (/) = (P./) recip = (P.recip) instance (PID.C a) => C (Ratio.T a) where+ {-# INLINE (/) #-}+ {-# INLINE recip #-}+ {-# INLINE fromRational' #-} -- (/) = Ratio.liftOrd (%) (x:%y) / (x':%y') = (x*y') % (y*x') recip (x:%y) = (y:%x)@@ -129,4 +141,7 @@ -- legacy instance (P.Integral a) => C (Ratio98.Ratio a) where- (/) = (P./)+ {-# INLINE (/) #-}+ {-# INLINE recip #-}+ (/) = (P./)+ recip = (P.recip)
+ src/Algebra/GenerateRules.hs view
@@ -0,0 +1,86 @@+{- |+Poor man's Template Haskell:+Generate RULES for handling of primitive number types.+-}+module Main where++import Data.Maybe (fromMaybe, )++import Prelude hiding (fromIntegral, )+++pad :: Int -> String -> String+pad n str =+ zipWith fromMaybe+ (replicate n ' ')+ (map Just str ++ repeat Nothing)+++machineIntegerTypes :: [String]+machineIntegerTypes =+ do typeSign <- "Int" : "Word" : []+ typeSize <- "" : "8" : "16" : "32" : "64" : []+ return $ typeSign ++ typeSize++functionSignature :: String -> String -> String -> String+functionSignature functionName sourceType targetType =+ functionName ++ " :: " ++ sourceType ++ " -> " ++ targetType++{-+Simply replace NumericPrelude.roundFunc by Prelude98.roundFunc.+This is only sensible where Prelude functions are optimized.+Unfortunately there seems to be no optimization for target type Int8 et.al.+-}+realField :: [String]+realField =+ do sourceType <- "Float" : "Double" : []+ targetType <- machineIntegerTypes+ method <- "round" : "truncate" : "floor" : "ceiling" : []+ let methodPad = pad 8 method+ let signature = functionSignature methodPad sourceType targetType+ return $ " " +++ pad 40 ("\"NP." ++ signature ++ "\"") +++ methodPad ++ " = P." ++ signature ++ ";"++realFieldIndirect :: [String]+realFieldIndirect =+ do targetType <- tail machineIntegerTypes+ method <- "round" : "truncate" : "floor" : "ceiling" : []+ let methodPad = pad 8 method+ let signature = functionSignature methodPad "a" targetType+ return $ " " +++ pad 30 ("\"NP." ++ signature ++ "\"") +++ methodPad ++ " = (" ++ functionSignature "P.fromIntegral" "Int" targetType ++ ") . "+ ++ method ++ ";"++splitFractionIndirect :: [String]+splitFractionIndirect =+ do targetType <- tail machineIntegerTypes+ method <- "splitFraction" : []+ let methodPad = pad 13 method+ let signature = functionSignature methodPad "a" ("("++targetType++",a)")+ return $ " " +++ pad 40 ("\"NP." ++ signature ++ "\"") +++ methodPad ++ " = mapFst (" ++ functionSignature "P.fromIntegral" "Int" targetType ++ ") . "+ ++ method ++ ";"+++fromIntegral :: [String]+fromIntegral =+ do sourceType <- "Integer" : machineIntegerTypes+ targetType <- "Int" : "Integer" : "Float" : "Double" : []+ let function = "fromIntegral"+ let signature = functionSignature function sourceType targetType+ return $ " " +++ pad 40 ("\"NP." ++ signature ++ "\"") +++ function ++ " = P." ++ signature ++ ";"+++main :: IO ()+main =+ putStrLn "module Algebra.RealField" >>+ mapM_ putStrLn realFieldIndirect >>+ mapM_ putStrLn splitFractionIndirect >>++ putStrLn "module Algebra.ToInteger" >>+ mapM_ putStrLn fromIntegral
src/Algebra/IntegralDomain.hs view
@@ -40,6 +40,9 @@ import Test.QuickCheck ((==>), Property) +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import PreludeBase import Prelude (Integer, Int) import qualified Prelude as P@@ -91,14 +94,19 @@ div, mod :: a -> a -> a divMod :: a -> a -> (a,a) + {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-} div a b = fst (divMod a b) mod a b = snd (divMod a b) divMod a b = (div a b, mod a b) +{-# INLINE divides #-} divides :: (C a, ZeroTestable.C a) => a -> a -> Bool divides y x = isZero (mod x y) +{-# INLINE sameResidueClass #-} sameResidueClass :: (C a, ZeroTestable.C a) => a -> a -> a -> Bool sameResidueClass m x y = divides m (x-y) @@ -134,6 +142,7 @@ Returns the result of the division, if divisible. Otherwise undefined. -}+{-# INLINE safeDiv #-} safeDiv :: (ZeroTestable.C a, C a) => a -> a -> a safeDiv a b = let (q,r) = divMod a b@@ -145,6 +154,7 @@ Allows division by zero. If the divisor is zero, then the divident is returned as remainder. -}+{-# INLINE divModZero #-} divModZero :: (C a, ZeroTestable.C a) => a -> a -> (a,a) divModZero x y = if isZero y@@ -153,6 +163,8 @@ +{-# INLINE even #-}+{-# INLINE odd #-} even, odd :: (C a, ZeroTestable.C a) => a -> Bool even n = divides 2 n odd = not . even@@ -161,10 +173,93 @@ {- * Instances for atomic types -} instance C Integer where- divMod = P.divMod+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod -instance C Int where- divMod = P.divMod+instance C Int where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Int8 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Int16 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Int32 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Int64 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod+++instance C Word where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Word8 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Word16 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Word32 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod++instance C Word64 where+ {-# INLINE div #-}+ {-# INLINE mod #-}+ {-# INLINE divMod #-}+ div = P.div+ mod = P.mod+ divMod = P.divMod
src/Algebra/Module.hs view
@@ -52,21 +52,27 @@ {-* Instances for atomic types -} instance C Float Float where+ {-# INLINE (*>) #-} (*>) = (*) instance C Double Double where+ {-# INLINE (*>) #-} (*>) = (*) instance C Int Int where+ {-# INLINE (*>) #-} (*>) = (*) instance C Integer Integer where+ {-# INLINE (*>) #-} (*>) = (*) instance (PID.C a) => C (Ratio.T a) (Ratio.T a) where+ {-# INLINE (*>) #-} (*>) = (*) instance (PID.C a) => C Integer (Ratio.T a) where+ {-# INLINE (*>) #-} x *> y = fromInteger x * y @@ -74,15 +80,19 @@ {-* Instances for composed types -} instance (C a b0, C a b1) => C a (b0, b1) where+ {-# INLINE (*>) #-} s *> (x0,x1) = (s *> x0, s *> x1) instance (C a b0, C a b1, C a b2) => C a (b0, b1, b2) where+ {-# INLINE (*>) #-} s *> (x0,x1,x2) = (s *> x0, s *> x1, s *> x2) instance (C a b) => C a [b] where+ {-# INLINE (*>) #-} (*>) = map . (*>) instance (C a b) => C a (c -> b) where+ {-# INLINE (*>) #-} (*>) s f = (*>) s . f @@ -103,6 +113,7 @@ Better move to "Algebra.Additive"? -}+{-# INLINE integerMultiply #-} integerMultiply :: (ToInteger.C a, Additive.C b) => a -> b -> b integerMultiply a b = reduceRepeated (+) zero b (ToInteger.toInteger a)
src/Algebra/NonNegative.hs view
@@ -16,7 +16,9 @@ module Algebra.NonNegative (C(..)) where import qualified Algebra.Additive as Additive+import qualified Algebra.Real as Real +infixl 6 -|, -? {- | Instances of this class must ensure non-negative values.@@ -35,3 +37,7 @@ -} (-|) :: a -> a -> a x -| y = if x >= y then x Additive.- y else Additive.zero+++(-?) :: (Real.C a) => a -> a -> (Bool, a)+(-?) x y = (x >= y, Real.abs (x Additive.- y))
src/Algebra/Real.hs view
@@ -10,7 +10,11 @@ import Algebra.Ring (one, ) -- fromInteger import Algebra.Additive (zero, negate,) +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import PreludeBase+import qualified Prelude as P import Prelude(Int,Integer,Float,Double) @@ -42,7 +46,83 @@ LT -> negate one -instance C Integer-instance C Int-instance C Float-instance C Double+instance C Integer where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Float where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Double where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum+++instance C Int where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Int8 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Int16 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Int32 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Int64 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum+++instance C Word where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Word8 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Word16 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Word32 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum++instance C Word64 where+ {-# INLINE abs #-}+ {-# INLINE signum #-}+ abs = P.abs+ signum = P.signum+
src/Algebra/RealField.hs view
@@ -1,4 +1,5 @@-{-# OPTIONS -fno-implicit-prelude #-}+{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-}+-- -fglasgow-exts for RULES module Algebra.RealField where import qualified Algebra.Field as Field@@ -9,7 +10,7 @@ import qualified Algebra.ToInteger as ToInteger -import Algebra.Field ((/))+import Algebra.Field ((/), fromRational, ) import Algebra.RealIntegral (quotRem, ) import Algebra.IntegralDomain (divMod, even, ) import Algebra.Ring ((*), fromInteger, )@@ -20,8 +21,12 @@ 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 Prelude(Int,Float,Double)+import NumericPrelude.Tuple (mapFst, )+import Prelude(Int, Integer, Float, Double) import qualified Prelude as P import PreludeBase @@ -100,43 +105,183 @@ where (q,r) = divMod x y instance C Float where- splitFraction = preludeSplitFraction- fraction = fractionTrunc (GHC.int2Float . GHC.float2Int)- -- preludeFraction+ {-# 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- splitFraction = preludeSplitFraction- fraction = fractionTrunc (GHC.int2Double . GHC.double2Int)+ {-# 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 -preludeSplitFraction :: (P.RealFrac a, Ring.C a, ToInteger.C b) => a -> (b,a)-preludeSplitFraction x =- let (n,f) = P.properFraction x++{-# 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- in if f>=0- then (fromInteger n, f)- else (fromInteger n-1, f+1)+ 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 :: (Int, a) -> a+ let second :: (Integer, a) -> a second = snd- in fixFraction (second (P.properFraction x))--fractionTrunc :: (Ring.C a, Ord a) => (a -> a) -> a -> a-fractionTrunc trunc x =- fixFraction (x - trunc x)+ 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. -}
src/Algebra/RealIntegral.hs view
@@ -26,7 +26,11 @@ import Algebra.Ring (one, ) -- fromInteger import Algebra.Additive (zero, (+), (-), ) +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import PreludeBase+import qualified Prelude as P import Prelude (Int, Integer, ) @@ -46,6 +50,9 @@ quot, rem :: a -> a -> a quotRem :: a -> a -> (a,a) + {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-} quot a b = fst (quotRem a b) rem a b = snd (quotRem a b) quotRem a b = let (d,m) = divMod a b in@@ -53,5 +60,92 @@ (d+one,m-b) else (d,m) -instance C Integer-instance C Int+instance C Integer where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Int where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Int8 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Int16 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Int32 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Int64 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem+++instance C Word where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Word8 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Word16 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Word32 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem++instance C Word64 where+ {-# INLINE quot #-}+ {-# INLINE rem #-}+ {-# INLINE quotRem #-}+ quot = P.quot+ rem = P.rem+ quotRem = P.quotRem+
src/Algebra/Ring.hs view
@@ -32,6 +32,9 @@ import Test.QuickCheck ((==>), Property) +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import PreludeBase import Prelude(Integer,Int,Float,Double) import qualified Data.Ratio as Ratio98@@ -78,12 +81,15 @@ -} (^) :: a -> Integer -> a + {-# INLINE fromInteger #-} fromInteger n = if n < 0 then reduceRepeated (+) zero (negate one) (negate n) else reduceRepeated (+) zero one n+ {-# INLINE (^) #-} a ^ n = if n >= zero then reduceRepeated (*) one a n else error "(^): Illegal negative exponent"+ {-# INLINE one #-} one = fromInteger 1 @@ -104,25 +110,115 @@ {- * Instances for atomic types -} instance C Integer where- (*) = (P.*)- one = P.fromInteger 1- fromInteger = P.fromInteger+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*) -instance C Int where- (*) = (P.*)- one = P.fromInteger 1+instance C Float where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*) -instance C Float where- (*) = (P.*)- one = P.fromInteger 1+instance C Double where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*) -instance C Double where- (*) = (P.*)- one = P.fromInteger 1 +instance C Int where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*) +instance C Int8 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*) +instance C Int16 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*) +instance C Int32 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)++instance C Int64 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)+++instance C Word where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)++instance C Word8 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)++instance C Word16 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)++instance C Word32 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)++instance C Word64 where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-}+ one = P.fromInteger 1+ fromInteger = P.fromInteger+ (*) = (P.*)+++++ propAssociative :: (Eq a, C a) => a -> a -> a -> Bool propLeftDistributive :: (Eq a, C a) => a -> a -> a -> Bool propRightDistributive :: (Eq a, C a) => a -> a -> a -> Bool@@ -152,6 +248,9 @@ -- legacy instance (P.Integral a) => C (Ratio98.Ratio a) where+ {-# INLINE one #-}+ {-# INLINE fromInteger #-}+ {-# INLINE (*) #-} one = 1 fromInteger = P.fromInteger (*) = (P.*)
src/Algebra/ToInteger.hs view
@@ -1,3 +1,5 @@+{-# OPTIONS -fglasgow-exts #-}+-- -fglasgow-exts for RULES module Algebra.ToInteger where import qualified Number.Ratio as Ratio@@ -13,9 +15,12 @@ import Algebra.Field ((^-), ) import Algebra.Ring ((^), fromInteger, ) +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import qualified Prelude as P import PreludeBase-import Prelude(Int,Integer)+import Prelude (Int, Integer, Float, Double, ) {- |@@ -43,11 +48,69 @@ fromIntegral = fromInteger . toInteger -instance C Integer where- toInteger = id+-- generated by GenerateRules.hs+{-# RULES+ "NP.fromIntegral :: Integer -> Int" fromIntegral = P.fromIntegral :: Integer -> Int;+ "NP.fromIntegral :: Integer -> Integer" fromIntegral = P.fromIntegral :: Integer -> Integer;+ "NP.fromIntegral :: Integer -> Float" fromIntegral = P.fromIntegral :: Integer -> Float;+ "NP.fromIntegral :: Integer -> Double" fromIntegral = P.fromIntegral :: Integer -> Double;+ "NP.fromIntegral :: Int -> Int" fromIntegral = P.fromIntegral :: Int -> Int;+ "NP.fromIntegral :: Int -> Integer" fromIntegral = P.fromIntegral :: Int -> Integer;+ "NP.fromIntegral :: Int -> Float" fromIntegral = P.fromIntegral :: Int -> Float;+ "NP.fromIntegral :: Int -> Double" fromIntegral = P.fromIntegral :: Int -> Double;+ "NP.fromIntegral :: Int8 -> Int" fromIntegral = P.fromIntegral :: Int8 -> Int;+ "NP.fromIntegral :: Int8 -> Integer" fromIntegral = P.fromIntegral :: Int8 -> Integer;+ "NP.fromIntegral :: Int8 -> Float" fromIntegral = P.fromIntegral :: Int8 -> Float;+ "NP.fromIntegral :: Int8 -> Double" fromIntegral = P.fromIntegral :: Int8 -> Double;+ "NP.fromIntegral :: Int16 -> Int" fromIntegral = P.fromIntegral :: Int16 -> Int;+ "NP.fromIntegral :: Int16 -> Integer" fromIntegral = P.fromIntegral :: Int16 -> Integer;+ "NP.fromIntegral :: Int16 -> Float" fromIntegral = P.fromIntegral :: Int16 -> Float;+ "NP.fromIntegral :: Int16 -> Double" fromIntegral = P.fromIntegral :: Int16 -> Double;+ "NP.fromIntegral :: Int32 -> Int" fromIntegral = P.fromIntegral :: Int32 -> Int;+ "NP.fromIntegral :: Int32 -> Integer" fromIntegral = P.fromIntegral :: Int32 -> Integer;+ "NP.fromIntegral :: Int32 -> Float" fromIntegral = P.fromIntegral :: Int32 -> Float;+ "NP.fromIntegral :: Int32 -> Double" fromIntegral = P.fromIntegral :: Int32 -> Double;+ "NP.fromIntegral :: Int64 -> Int" fromIntegral = P.fromIntegral :: Int64 -> Int;+ "NP.fromIntegral :: Int64 -> Integer" fromIntegral = P.fromIntegral :: Int64 -> Integer;+ "NP.fromIntegral :: Int64 -> Float" fromIntegral = P.fromIntegral :: Int64 -> Float;+ "NP.fromIntegral :: Int64 -> Double" fromIntegral = P.fromIntegral :: Int64 -> Double;+ "NP.fromIntegral :: Word -> Int" fromIntegral = P.fromIntegral :: Word -> Int;+ "NP.fromIntegral :: Word -> Integer" fromIntegral = P.fromIntegral :: Word -> Integer;+ "NP.fromIntegral :: Word -> Float" fromIntegral = P.fromIntegral :: Word -> Float;+ "NP.fromIntegral :: Word -> Double" fromIntegral = P.fromIntegral :: Word -> Double;+ "NP.fromIntegral :: Word8 -> Int" fromIntegral = P.fromIntegral :: Word8 -> Int;+ "NP.fromIntegral :: Word8 -> Integer" fromIntegral = P.fromIntegral :: Word8 -> Integer;+ "NP.fromIntegral :: Word8 -> Float" fromIntegral = P.fromIntegral :: Word8 -> Float;+ "NP.fromIntegral :: Word8 -> Double" fromIntegral = P.fromIntegral :: Word8 -> Double;+ "NP.fromIntegral :: Word16 -> Int" fromIntegral = P.fromIntegral :: Word16 -> Int;+ "NP.fromIntegral :: Word16 -> Integer" fromIntegral = P.fromIntegral :: Word16 -> Integer;+ "NP.fromIntegral :: Word16 -> Float" fromIntegral = P.fromIntegral :: Word16 -> Float;+ "NP.fromIntegral :: Word16 -> Double" fromIntegral = P.fromIntegral :: Word16 -> Double;+ "NP.fromIntegral :: Word32 -> Int" fromIntegral = P.fromIntegral :: Word32 -> Int;+ "NP.fromIntegral :: Word32 -> Integer" fromIntegral = P.fromIntegral :: Word32 -> Integer;+ "NP.fromIntegral :: Word32 -> Float" fromIntegral = P.fromIntegral :: Word32 -> Float;+ "NP.fromIntegral :: Word32 -> Double" fromIntegral = P.fromIntegral :: Word32 -> Double;+ "NP.fromIntegral :: Word64 -> Int" fromIntegral = P.fromIntegral :: Word64 -> Int;+ "NP.fromIntegral :: Word64 -> Integer" fromIntegral = P.fromIntegral :: Word64 -> Integer;+ "NP.fromIntegral :: Word64 -> Float" fromIntegral = P.fromIntegral :: Word64 -> Float;+ "NP.fromIntegral :: Word64 -> Double" fromIntegral = P.fromIntegral :: Word64 -> Double;+ #-} -instance C Int where- toInteger = P.toInteger++instance C Integer where {-#INLINE toInteger #-}; toInteger = id++instance C Int where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Int8 where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Int16 where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Int32 where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Int64 where {-#INLINE toInteger #-}; toInteger = P.toInteger++instance C Word where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Word8 where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Word16 where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Word32 where {-#INLINE toInteger #-}; toInteger = P.toInteger+instance C Word64 where {-#INLINE toInteger #-}; toInteger = P.toInteger+ instance (C a, PID.C a) => ToRational.C (Ratio.T a) where toRational (x:%y) = toInteger x :% toInteger y
src/Algebra/ToRational.hs view
@@ -7,6 +7,9 @@ import Number.Ratio (Rational, ) +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )+ import qualified Prelude as P import PreludeBase import Prelude(Int,Integer,Float,Double)@@ -30,13 +33,25 @@ toRational :: a -> Rational instance C Integer where+ {-#INLINE toRational #-} toRational = fromInteger -instance C Int where- toRational = toRational . P.toInteger- instance C Float where+ {-#INLINE toRational #-} toRational = fromRational . P.toRational instance C Double where+ {-#INLINE toRational #-} toRational = fromRational . P.toRational++instance C Int where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int8 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int16 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int32 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int64 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger++instance C Word where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word8 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word16 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word32 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word64 where {-#INLINE toRational #-}; toRational = toRational . P.toInteger
src/Algebra/Transcendental.hs view
@@ -16,7 +16,7 @@ import PreludeBase -infixr 8 **+infixr 8 **, ^? {-| Transcendental is the type of numbers supporting the elementary@@ -42,6 +42,24 @@ sinh, cosh, tanh :: a -> a asinh, acosh, atanh :: a -> a + {-# INLINE pi #-}+ {-# INLINE exp #-}+ {-# INLINE log #-}+ {-# INLINE logBase #-}+ {-# INLINE (**) #-}+ {-# INLINE sin #-}+ {-# INLINE tan #-}+ {-# INLINE cos #-}+ {-# INLINE asin #-}+ {-# INLINE atan #-}+ {-# INLINE acos #-}+ {-# INLINE sinh #-}+ {-# INLINE tanh #-}+ {-# INLINE cosh #-}+ {-# INLINE asinh #-}+ {-# INLINE atanh #-}+ {-# INLINE acosh #-}+ x ** y = exp (log x * y) logBase x y = log y / log x @@ -62,6 +80,24 @@ instance C P.Float where+ {-# INLINE pi #-}+ {-# INLINE exp #-}+ {-# INLINE log #-}+ {-# INLINE logBase #-}+ {-# INLINE (**) #-}+ {-# INLINE sin #-}+ {-# INLINE tan #-}+ {-# INLINE cos #-}+ {-# INLINE asin #-}+ {-# INLINE atan #-}+ {-# INLINE acos #-}+ {-# INLINE sinh #-}+ {-# INLINE tanh #-}+ {-# INLINE cosh #-}+ {-# INLINE asinh #-}+ {-# INLINE atanh #-}+ {-# INLINE acosh #-}+ (**) = (P.**) exp = P.exp; log = P.log; logBase = P.logBase pi = P.pi;@@ -71,6 +107,24 @@ asinh = P.asinh; acosh = P.acosh; atanh = P.atanh instance C P.Double where+ {-# INLINE pi #-}+ {-# INLINE exp #-}+ {-# INLINE log #-}+ {-# INLINE logBase #-}+ {-# INLINE (**) #-}+ {-# INLINE sin #-}+ {-# INLINE tan #-}+ {-# INLINE cos #-}+ {-# INLINE asin #-}+ {-# INLINE atan #-}+ {-# INLINE acos #-}+ {-# INLINE sinh #-}+ {-# INLINE tanh #-}+ {-# INLINE cosh #-}+ {-# INLINE asinh #-}+ {-# INLINE atanh #-}+ {-# INLINE acosh #-}+ (**) = (P.**) exp = P.exp; log = P.log; logBase = P.logBase pi = P.pi;@@ -78,6 +132,12 @@ asin = P.asin; acos = P.acos; atan = P.atan sinh = P.sinh; cosh = P.cosh; tanh = P.tanh asinh = P.asinh; acosh = P.acosh; atanh = P.atanh++++{-# INLINE (^?) #-}+(^?) :: C a => a -> a -> a+(^?) = (**) {-* Transcendental laws, will only hold approximately on floating point numbers -}
src/Algebra/VectorSpace.hs view
@@ -1,14 +1,16 @@ {-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} module Algebra.VectorSpace where -import qualified Algebra.Module-import qualified Algebra.Field as Field+import qualified Algebra.Module as Module+import qualified Algebra.Field as Field+import qualified Algebra.PrincipalIdealDomain as PID+import qualified Number.Ratio as Ratio -- import NumericPrelude import qualified Prelude as P -class (Field.C a, Algebra.Module.C a b) => C a b+class (Field.C a, Module.C a b) => C a b {-* Instances for atomic types -}@@ -18,6 +20,8 @@ instance C P.Double P.Double {-* Instances for composed types -}++instance (PID.C a) => C (Ratio.T a) (Ratio.T a) instance (C a b0, C a b1) => C a (b0, b1)
src/Algebra/ZeroTestable.hs view
@@ -1,11 +1,15 @@ {-# OPTIONS -fno-implicit-prelude #-} module Algebra.ZeroTestable where -import qualified Prelude as P-import PreludeBase- import qualified Algebra.Additive as Additive +import Data.Int (Int, Int8, Int16, Int32, Int64, )+import Data.Word (Word, Word8, Word16, Word32, Word64, )++-- import qualified Prelude as P+import Prelude(Int,Integer,Float,Double)+import PreludeBase+ {- | Maybe the naming should be according to Algebra.Unit: Algebra.Zero as module name, and @query@ as method name.@@ -30,17 +34,22 @@ {-* Instances for atomic types -} -instance C P.Integer where- isZero = defltIsZero+instance C Integer where isZero = defltIsZero+instance C Float where isZero = defltIsZero+instance C Double where isZero = defltIsZero -instance C P.Int where- isZero = defltIsZero+instance C Int where isZero = defltIsZero+instance C Int8 where isZero = defltIsZero+instance C Int16 where isZero = defltIsZero+instance C Int32 where isZero = defltIsZero+instance C Int64 where isZero = defltIsZero -instance C P.Float where- isZero = defltIsZero+instance C Word where isZero = defltIsZero+instance C Word8 where isZero = defltIsZero+instance C Word16 where isZero = defltIsZero+instance C Word32 where isZero = defltIsZero+instance C Word64 where isZero = defltIsZero -instance C P.Double where- isZero = defltIsZero {-* Instances for composed types -}
src/MathObj/LaurentPolynomial.hs view
@@ -59,10 +59,10 @@ fromShiftCoeffs = Cons fromPolynomial :: Poly.T a -> T a-fromPolynomial (Poly.Cons xs) = fromCoeffs xs+fromPolynomial = fromCoeffs . Poly.coeffs fromPowerSeries :: PS.T a -> T a-fromPowerSeries (PS.Cons xs) = fromCoeffs xs+fromPowerSeries = fromCoeffs . PS.coeffs bounds :: T a -> (Int, Int) bounds (Cons xt x) = (xt, xt + length x - 1)
src/MathObj/Polynomial.hs view
@@ -41,15 +41,18 @@ may not be representable by radicals. -} -module MathObj.Polynomial(T(..), fromCoeffs, showsExpressionPrec, const,- eval, compose, equal, add, sub, negate,- shift, unShift,- mul, scale, divMod,- tensorProduct, tensorProduct',- mulShear, mulShearTranspose,- horner, horner',- progression, differentiate, integrate, integrateInt,- fromRoots, alternate)+module MathObj.Polynomial+ (T, fromCoeffs, coeffs,+ showsExpressionPrec, const,+ evaluate, evaluateCoeffVector, evaluateArgVector,+ compose, equal, add, sub, negate,+ horner, hornerCoeffVector, hornerArgVector,+ shift, unShift,+ mul, scale, divMod,+ tensorProduct, tensorProductAlt,+ mulShear, mulShearTranspose,+ progression, differentiate, integrate, integrateInt,+ fromRoots, alternate) where import qualified Algebra.Differential as Differential@@ -82,17 +85,23 @@ import PreludeBase hiding (const) import NumericPrelude hiding (divMod, negate, stdUnit) + newtype T a = Cons {coeffs :: [a]} ++{-# INLINE fromCoeffs #-} fromCoeffs :: [a] -> T a fromCoeffs = lift0 +{-# INLINE lift0 #-} lift0 :: [a] -> T a lift0 = Cons +{-# INLINE lift1 #-} lift1 :: ([a] -> [a]) -> (T a -> T a) lift1 f (Cons x0) = Cons (f x0) +{-# INLINE lift2 #-} lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a) lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1) @@ -104,6 +113,8 @@ instance Functor T where fmap f (Cons xs) = Cons (map f xs) +{-# INLINE plusPrec #-}+{-# INLINE appPrec #-} plusPrec, appPrec :: Int plusPrec = 6 appPrec = 10@@ -112,6 +123,7 @@ showsPrec p (Cons xs) = showParen (p >= appPrec) (showString "Polynomial.fromCoeffs " . shows xs) +{-# INLINE showsExpressionPrec #-} showsExpressionPrec :: (Show a, ZeroTestable.C a, Additive.C a) => Int -> String -> T a -> String -> String showsExpressionPrec p var poly =@@ -131,25 +143,46 @@ map (uncurry showsTerm) terms) {- |-Horner's scheme for evaluating an polynomial+Horner's scheme for evaluating a polynomial in a ring. -}--horner' :: Ring.C a => a -> [a] -> a-horner' x = foldr (\c val -> c+x*val) zero+{-# INLINE horner #-}+horner :: Ring.C a => a -> [a] -> a+horner x = foldr (\c val -> c+x*val) zero -{--***** Module types are more general,-but most times this flexibility is not needed and-let type inference fail.+{- |+Horner's scheme for evaluating a polynomial in a module. -}+{-# INLINE hornerCoeffVector #-}+hornerCoeffVector :: Module.C a v => a -> [v] -> v+hornerCoeffVector x = foldr (\c val -> c+x*>val) zero -horner :: Module.C a b => a -> [b] -> b-horner x = foldr (\c val -> c+x*>val) zero+{-# INLINE hornerArgVector #-}+hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v+hornerArgVector x = foldr (\c val -> c*>one+val*x) zero -eval :: Module.C a b => T b -> a -> b-eval (Cons y) x = horner x y +{-# INLINE evaluate #-}+evaluate :: Ring.C a => T a -> a -> a+evaluate (Cons y) x = horner x y+ {- |+Here the coefficients are vectors,+for example the coefficients are real and the coefficents are real vectors.+-}+{-# INLINE evaluateCoeffVector #-}+evaluateCoeffVector :: Module.C a v => T v -> a -> v+evaluateCoeffVector (Cons y) x = hornerCoeffVector x y++{- |+Here the argument is a vector,+for example the coefficients are complex numbers or square matrices+and the coefficents are reals.+-}+{-# INLINE evaluateArgVector #-}+evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v+evaluateArgVector (Cons y) x = hornerArgVector x y++{- | 'compose' is the functional composition of polynomials. It fulfills@@ -158,14 +191,16 @@ -- compose :: Module.C a b => T b -> T a -> T a -- compose (Cons x) y = horner y (map const x)+{-# INLINE compose #-} compose :: (Ring.C a) => T a -> T a -> T a-compose (Cons x) y = horner' y (map const x)+compose (Cons x) y = horner y (map const x) {- | It's also helpful to put a polynomial in canonical form. 'normalize' strips leading coefficients that are zero. -} +{-# INLINE normalize #-} normalize :: (ZeroTestable.C a) => [a] -> [a] normalize = dropWhileRev isZero @@ -173,17 +208,21 @@ Multiply by the variable, used internally. -} +{-# INLINE shift #-} shift :: (Additive.C a) => [a] -> [a] shift [] = [] shift l = zero : l +{-# INLINE unShift #-} unShift :: [a] -> [a] unShift [] = [] unShift (_:xs) = xs +{-# INLINE const #-} const :: a -> T a const x = lift0 [x] +{-# INLINE equal #-} equal :: (Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool equal x y = and (zipWithOverlap isZero isZero (==) x y) @@ -201,6 +240,7 @@ add = (+) sub = (-) +{-# INLINE negate #-} negate :: (Additive.C a) => [a] -> [a] negate = map NP.negate @@ -211,6 +251,7 @@ negate = lift1 negate +{-# INLINE scale #-} scale :: Ring.C a => a -> [a] -> [a] scale s = map (s*) @@ -226,17 +267,19 @@ instance (Field.C a, Module.C a b) => VectorSpace.C a (T b) +{-# INLINE tensorProduct #-} tensorProduct :: Ring.C a => [a] -> [a] -> [[a]] tensorProduct = outerProduct (*) -tensorProduct' :: Ring.C a => [a] -> [a] -> [[a]]-tensorProduct' xs ys = map (flip scale ys) xs+tensorProductAlt :: Ring.C a => [a] -> [a] -> [[a]]+tensorProductAlt xs ys = map (flip scale ys) xs {- | 'mul' is fast if the second argument is a short polynomial, 'MathObj.PowerSeries.**' relies on that fact. -} +{-# INLINE mul #-} mul :: Ring.C a => [a] -> [a] -> [a] {- prevent from generation of many zeros if the first operand is the empty list -}@@ -245,9 +288,11 @@ -- this one fails on infinite lists -- mul xs = foldr (\y zs -> add (scale y xs) (shift zs)) [] +{-# INLINE mulShear #-} mulShear :: Ring.C a => [a] -> [a] -> [a] mulShear xs ys = map sum (shear (tensorProduct xs ys)) +{-# INLINE mulShearTranspose #-} mulShearTranspose :: Ring.C a => [a] -> [a] -> [a] mulShearTranspose xs ys = map sum (shearTranspose (tensorProduct xs ys)) @@ -278,6 +323,7 @@ let (d,m) = divMod x y in (Cons d, Cons m) +{-# INLINE stdUnit #-} stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a stdUnit x = case normalize x of [] -> one@@ -296,12 +342,15 @@ instance (ZeroTestable.C a, Field.C a) => PID.C (T a) +{-# INLINE progression #-} progression :: Ring.C a => [a] progression = iterate (one+) one +{-# INLINE differentiate #-} differentiate :: (Ring.C a) => [a] -> [a] differentiate = zipWith (*) progression . tail +{-# INLINE integrate #-} integrate :: (Field.C a) => a -> [a] -> [a] integrate c x = c : zipWith (/) x progression @@ -310,6 +359,7 @@ in the given ring. Otherwise undefined coefficients occur. -}+{-# INLINE integrateInt #-} integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a] integrateInt c x = c : zipWith Integral.safeDiv x progression@@ -319,13 +369,16 @@ differentiate = lift1 differentiate +{-# INLINE fromRoots #-} fromRoots :: (Ring.C a) => [a] -> T a fromRoots = Cons . foldl (flip mulLinearFactor) [1] +{-# INLINE mulLinearFactor #-} mulLinearFactor :: Ring.C a => a -> [a] -> [a] mulLinearFactor x yt@(y:ys) = Additive.negate (x*y) : yt - scale x ys mulLinearFactor _ [] = [] +{-# INLINE alternate #-} alternate :: Additive.C a => [a] -> [a] alternate = zipWith ($) (cycle [id, Additive.negate]) @@ -357,6 +410,7 @@ in (x^2+x+1)*(x-1) However the output looks much different. -}+{-# INLINE legacyInstance #-} legacyInstance :: a legacyInstance = error "legacy Ring.C instance for simple input of numeric literals"
src/MathObj/PowerSeries.hs view
@@ -33,20 +33,26 @@ sqrt, exp, log, sin, cos, tan, asin, acos, atan) + newtype T a = Cons {coeffs :: [a]} deriving (Ord) +{-# INLINE fromCoeffs #-} fromCoeffs :: [a] -> T a fromCoeffs = lift0 +{-# INLINE lift0 #-} lift0 :: [a] -> T a lift0 = Cons +{-# INLINE lift1 #-} lift1 :: ([a] -> [a]) -> (T a -> T a) lift1 f (Cons x0) = Cons (f x0) +{-# INLINE lift2 #-} lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a) lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1) +{-# INLINE const #-} const :: a -> T a const x = lift0 [x] @@ -58,6 +64,7 @@ instance Functor T where fmap f (Cons xs) = Cons (map f xs) +{-# INLINE appPrec #-} appPrec :: Int appPrec = 10 @@ -66,29 +73,79 @@ showParen (p >= appPrec) (showString "PowerSeries.fromCoeffs " . shows xs) +{-# INLINE truncate #-} truncate :: Int -> T a -> T a truncate n = lift1 (take n) {- | Evaluate (truncated) power series. -}--eval :: Module.C a b => [b] -> a -> b+{-# INLINE eval #-}+eval :: Ring.C a => [a] -> a -> a eval = flip Poly.horner -evaluate :: Module.C a b => T b -> a -> b+{-# INLINE evaluate #-}+evaluate :: Ring.C a => T a -> a -> a evaluate (Cons y) = eval y {- |-Evaluate approximations that is evaluated all truncations of the series.+Evaluate (truncated) power series. -}+{-# INLINE evalCoeffVector #-}+evalCoeffVector :: Module.C a v => [v] -> a -> v+evalCoeffVector = flip Poly.hornerCoeffVector -approx :: Module.C a b => [b] -> a -> [b]+{-# INLINE evaluateCoeffVector #-}+evaluateCoeffVector :: Module.C a v => T v -> a -> v+evaluateCoeffVector (Cons y) = evalCoeffVector y+++{-# INLINE evalArgVector #-}+evalArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> v+evalArgVector = flip Poly.hornerArgVector++{-# INLINE evaluateArgVector #-}+evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v+evaluateArgVector (Cons y) = evalArgVector y++{- |+Evaluate approximations that is evaluate all truncations of the series.+-}+{-# INLINE approx #-}+approx :: Ring.C a => [a] -> a -> [a] approx y x =- scanl (+) zero (zipWith (*>) (iterate (x*) 1) y)+ scanl (+) zero (zipWith (*) (iterate (x*) 1) y) -approximate :: Module.C a b => T b -> a -> [b]+{-# INLINE approximate #-}+approximate :: Ring.C a => T a -> a -> [a] approximate (Cons y) = approx y+++{- |+Evaluate approximations that is evaluate all truncations of the series.+-}+{-# INLINE approxCoeffVector #-}+approxCoeffVector :: Module.C a v => [v] -> a -> [v]+approxCoeffVector y x =+ scanl (+) zero (zipWith (*>) (iterate (x*) 1) y)++{-# INLINE approximateCoeffVector #-}+approximateCoeffVector :: Module.C a v => T v -> a -> [v]+approximateCoeffVector (Cons y) = approxCoeffVector y+++{- |+Evaluate approximations that is evaluate all truncations of the series.+-}+{-# INLINE approxArgVector #-}+approxArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> [v]+approxArgVector y x =+ scanl (+) zero (zipWith (*>) y (iterate (x*) 1))++{-# INLINE approximateArgVector #-}+approximateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> [v]+approximateArgVector (Cons y) = approxArgVector y+ {- * Simple series manipulation -}
src/MathObj/RootSet.hs view
@@ -52,10 +52,10 @@ toPolynomial :: Poly.T a -> T a-toPolynomial (Poly.Cons xs) = Cons (reverse xs)+toPolynomial xs = Cons (reverse (Poly.coeffs xs)) fromPolynomial :: T a -> Poly.T a-fromPolynomial (Cons xs) = Poly.Cons (reverse xs)+fromPolynomial (Cons xs) = Poly.fromCoeffs (reverse xs)
src/Number/Complex.hs view
@@ -1,4 +1,6 @@-{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-}+{-# OPTIONS_GHC -fno-implicit-prelude -fglasgow-exts #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# LANGUAGE_HOW_CAN_WE_ENABLE Rules ? #-} {- | Module : Number.Complex Copyright : (c) The University of Glasgow 2001@@ -12,37 +14,33 @@ -} module Number.Complex- (- -- * Cartesian form- T(real,imag)- , fromReal+ (+ -- * Cartesian form+ T(real,imag),+ imaginaryUnit,+ fromReal, - , (+:)- , (-:)- -- * Polar form- , fromPolar- , cis- , signum- , toPolar- , magnitude- , phase- , Polar- , defltMagnitude- , defltPhase- -- * Conjugate- , conjugate+ (+:),+ (-:),+ -- * Polar form+ fromPolar,+ cis,+ signum,+ toPolar,+ magnitude,+ phase,+ -- * Conjugate+ conjugate, -- * Properties- , propPolar+ propPolar, -- * Auxiliary classes- , Divisible(divide)- , defltDiv- , Power(power)- , defltPow+ Power(power),+ defltPow, ) where -import qualified Number.Ratio as Ratio+-- import qualified Number.Ratio as Ratio import qualified Algebra.NormedSpace.Euclidean as NormedEuc import qualified Algebra.NormedSpace.Sum as NormedSum@@ -68,6 +66,9 @@ import Algebra.Module((*>)) import Algebra.Algebraic((^/)) +import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary)+import Control.Monad (liftM2)+ import qualified Prelude as P import PreludeBase import NumericPrelude hiding (signum)@@ -87,10 +88,16 @@ } deriving (Eq) +{-# INLINE imaginaryUnit #-}+imaginaryUnit :: Ring.C a => T a+imaginaryUnit = zero +: one++{-# INLINE fromReal #-} fromReal :: Additive.C a => a -> T a fromReal x = Cons x zero +{-# INLINE plusPrec #-} plusPrec :: Int plusPrec = 6 @@ -100,29 +107,40 @@ instance (Read a) => Read (T a) where readsPrec prec = readsInfixPrec "+:" plusPrec prec (+:) +instance (Arbitrary a) => Arbitrary (T a) where+ {-# INLINE arbitrary #-}+ arbitrary = liftM2 Cons arbitrary arbitrary+ {-# INLINE coarbitrary #-}+ coarbitrary = undefined {- * Functions -} -- | Construct a complex number from real and imaginary part.+{-# INLINE (+:) #-} (+:) :: a -> a -> T a (+:) = Cons -- | Construct a complex number with negated imaginary part.+{-# INLINE (-:) #-} (-:) :: Additive.C a => a -> a -> T a (-:) x y = Cons x (-y) -- | The conjugate of a complex number.-{-# SPECIALISE conjugate :: T Double -> T Double #-}-conjugate :: (Additive.C a) => T a -> T a+{-# SPECULATE 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.-{-# SPECIALISE scale :: Double -> T Double -> T Double #-}-scale :: (Ring.C a) => a -> T a -> T a+{-# SPECULATE 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) -- | Turn the point one quarter to the right.+{-# INLINE orthoRight #-}+{-# INLINE orthoLeft #-} orthoRight, orthoLeft :: (Additive.C a) => T a -> T a orthoRight (Cons x y) = Cons y (-x) orthoLeft (Cons x y) = Cons (-y) x@@ -133,7 +151,8 @@ but oriented in the positive real direction, whereas @'signum' z@ has the phase of @z@, but unit magnitude. -}-{-# SPECIALISE signum :: T Double -> T Double #-}+{-# SPECULATE signum :: T Double -> T Double #-}+{-# INLINE signum #-} signum :: (Algebraic.C a, NormedEuc.C a a, ZeroTestable.C a) => T a -> T a signum z = if isZero z@@ -141,84 +160,76 @@ else scale (recip (NormedEuc.norm z)) z -- | Form a complex number from polar components of magnitude and phase.-{-# SPECIALISE fromPolar :: Double -> Double -> T Double #-}-fromPolar :: (Trans.C a) => a -> a -> T a-fromPolar r theta = scale r (cis theta)+{-# SPECULATE 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'@).-{-# SPECIALISE cis :: Double -> T Double #-}-cis :: (Trans.C a) => a -> T a-cis theta = Cons (cos theta) (sin theta)+{-# SPECULATE cis :: Double -> T Double #-}+{-# INLINE cis #-}+cis :: (Trans.C a) => a -> T a+cis theta = Cons (cos theta) (sin theta) -propPolar :: (Polar a, RealTrans.C a) => T a -> Bool-propPolar z = uncurry fromPolar (toPolar z) == z+propPolar :: (RealTrans.C a) => T a -> Bool+propPolar z = uncurry fromPolar (toPolar z) == z --- | The nonnegative magnitude of a complex number.+{- |+The nonnegative magnitude of a complex number.+This implementation respects the limited range of floating point numbers.+The trivial implementation 'magnitude'+would overflow for floating point exponents greater than+the half of the maximum admissible exponent.+We automatically drop in this implementation for 'Float' and 'Double'+by optimizer rules.+You should do so for your custom floating point types.+-}+{-# INLINE floatMagnitude #-} floatMagnitude :: (P.RealFloat a, Algebraic.C a) => T a -> a-floatMagnitude (Cons x y) = P.scaleFloat k- (sqrt (P.scaleFloat mk x ^ 2 +- P.scaleFloat mk y ^ 2))- where k = max (P.exponent x) (P.exponent y)- mk = - k+floatMagnitude (Cons x y) =+ let k = max (P.exponent x) (P.exponent y)+ mk = - k+ in P.scaleFloat k+ (sqrt (P.scaleFloat mk x ^ 2 ++ P.scaleFloat mk y ^ 2)) -defltMagnitude :: (Algebraic.C a) => T a -> a-defltMagnitude = sqrt . defltMagnitudeSqr+{-# INLINE [1] magnitude #-}+magnitude :: (Algebraic.C a) => T a -> a+magnitude = sqrt . magnitudeSqr +{-# RULES+ "Complex.magnitude :: Double"+ magnitude = floatMagnitude :: T Double -> Double;++ "Complex.magnitude :: Float"+ magnitude = floatMagnitude :: T Float -> Float;+ #-}+ -- like NormedEuc.normSqr with lifted class constraints-defltMagnitudeSqr :: (Ring.C a) => T a -> a-defltMagnitudeSqr (Cons x y) = x^2 + y^2+{-# INLINE magnitudeSqr #-}+magnitudeSqr :: (Ring.C a) => T a -> a+magnitudeSqr (Cons x y) = x^2 + y^2 -- | The phase of a complex number, in the range @(-'pi', 'pi']@. -- If the magnitude is zero, then so is the phase.-defltPhase :: (RealTrans.C a, ZeroTestable.C a) => T a -> a-defltPhase z =+{-# INLINE phase #-}+phase :: (RealTrans.C a, ZeroTestable.C a) => T a -> a+phase z = if isZero z then zero -- SLPJ July 97 from John Peterson else case z of (Cons x y) -> atan2 y x {- |-Minimal implementation: toPolar or (magnitude and phase),-usually the instance definition--@-magnitude = defltMagnitude-phase = defltPhase-@--is enough.--This class requires transcendent number types-although 'magnitude' can be computed algebraically.+The function 'toPolar' takes a complex number and+returns a (magnitude, phase) pair in canonical form:+the magnitude is nonnegative, and the phase in the range @(-'pi', 'pi']@;+if the magnitude is zero, then so is the phase. -}-class RealTrans.C a => Polar a where- {- |- The function 'toPolar' takes a complex number and- returns a (magnitude, phase) pair in canonical form:- the magnitude is nonnegative, and the phase in the range @(-'pi', 'pi']@;- if the magnitude is zero, then so is the phase.- -}- {--# SPECIALISE toPolar :: T Double -> (Double,Double) #--}- toPolar :: T a -> (a,a)- toPolar z = (magnitude z, phase z)-- {--# SPECIALISE magnitude :: T Double -> Double #--}- magnitude :: T a -> a- magnitude = fst . toPolar-- {--# SPECIALISE phase :: T Double -> Double #--}- phase :: T a -> a- phase = snd . toPolar--instance Polar Float where- magnitude = floatMagnitude- phase = defltPhase--instance Polar Double where- magnitude = floatMagnitude- phase = defltPhase+toPolar :: (RealTrans.C a) => T a -> (a,a)+toPolar z = (magnitude z, phase z) @@ -231,48 +242,65 @@ -} instance (Indexable.C a) => Indexable.C (T a) where+ {-# INLINE compare #-} compare (Cons x y) (Cons x' y') = Indexable.compare (x,y) (x',y') instance (ZeroTestable.C a) => ZeroTestable.C (T a) where+ {-# INLINE isZero #-} isZero (Cons x y) = isZero x && isZero y instance (Additive.C a) => Additive.C (T a) where- {-# SPECIALISE instance Additive.C (T Float) #-}- {-# SPECIALISE instance Additive.C (T Double) #-}- zero = Cons zero zero- (Cons x y) + (Cons x' y') = Cons (x+x') (y+y')- (Cons x y) - (Cons x' y') = Cons (x-x') (y-y')- negate (Cons x y) = Cons (negate x) (negate y)+ {-# SPECULATE instance Additive.C (T Float) #-}+ {-# SPECULATE instance Additive.C (T Double) #-}+ {-# INLINE zero #-}+ zero = Cons zero zero+ {-# INLINE (+) #-}+ (Cons x y) + (Cons x' y') = Cons (x+x') (y+y')+ {-# INLINE (-) #-}+ (Cons x y) - (Cons x' y') = Cons (x-x') (y-y')+ {-# INLINE negate #-}+ negate (Cons x y) = Cons (negate x) (negate y) instance (Ring.C a) => Ring.C (T a) where- {-# SPECIALISE instance Ring.C (T Float) #-}- {-# SPECIALISE instance Ring.C (T Double) #-}- one = Cons one zero- (Cons x y) * (Cons x' y') = Cons (x*x'-y*y') (x*y'+y*x')- fromInteger = fromReal . fromInteger+ {-# SPECULATE instance Ring.C (T Float) #-}+ {-# SPECULATE instance Ring.C (T Double) #-}+ {-# INLINE one #-}+ one = Cons one zero+ {-# INLINE (*) #-}+ (Cons x y) * (Cons x' y') = Cons (x*x'-y*y') (x*y'+y*x')+ {-# INLINE fromInteger #-}+ fromInteger = fromReal . fromInteger instance Vector.C T where+ {-# INLINE zero #-} zero = zero+ {-# INLINE (<+>) #-} (<+>) = (+)+ {-# INLINE (*>) #-} (*>) = scale -- | The '(*>)' method can't replace 'scale' -- because it requires the Algebra.Module constraint instance (Module.C a b) => Module.C a (T b) where+ {-# INLINE (*>) #-} s *> (Cons x y) = Cons (s *> x) (s *> y) instance (VectorSpace.C a b) => VectorSpace.C a (T b) instance (Additive.C a, NormedSum.C a v) => NormedSum.C a (T v) where+ {-# INLINE norm #-} norm x = NormedSum.norm (real x) + NormedSum.norm (imag x) instance (NormedEuc.Sqr a b) => NormedEuc.Sqr a (T b) where+ {-# INLINE normSqr #-} normSqr x = NormedEuc.normSqr (real x) + NormedEuc.normSqr (imag x) instance (Algebraic.C a, NormedEuc.Sqr a b) => NormedEuc.C a (T b) where+ {-# INLINE norm #-} norm = NormedEuc.defltNorm instance (Ord a, NormedMax.C a v) => NormedMax.C a (T v) where+ {-# INLINE norm #-} norm x = max (NormedMax.norm (real x)) (NormedMax.norm (imag x)) @@ -287,7 +315,7 @@ instance (Integral.C a) => Integral.C (T a) where divMod z z' =- let denom = defltMagnitudeSqr z'+ let denom = magnitudeSqr z' zBig = z * conjugate z' re = divMod (real zBig) denom im = divMod (imag zBig) denom@@ -300,9 +328,10 @@ Thus the remainder has smaller magnitude than the divisor. This variant of divModCent can be used for Euclidean's algorithm. -}+{-# INLINE divModCent #-} divModCent :: (Ord a, Integral.C a) => T a -> T a -> (T a, T a) divModCent z z' =- let denom = defltMagnitudeSqr z'+ let denom = magnitudeSqr z' zBig = z * conjugate z' re = divMod (real zBig) denom im = divMod (imag zBig) denom@@ -313,16 +342,20 @@ (imag q + if 2 * imag r > denom then one else zero) in (q', z-q'*z') +{-# INLINE modCent #-} modCent :: (Ord a, Integral.C a) => T a -> T a -> T a modCent z z' = snd (divModCent z z') instance (Ord a, Units.C a) => Units.C (T a) where+ {-# INLINE isUnit #-} isUnit (Cons x y) = isUnit x && y==zero || isUnit y && x==zero+ {-# INLINE stdAssociate #-} stdAssociate z@(Cons x y) = let z' = if y<0 || y==0 && x<0 then negate z else z in if real z'<=0 then orthoRight z' else z'+ {-# INLINE stdUnit #-} stdUnit z@(Cons x y) = if z==zero then 1@@ -334,51 +367,45 @@ instance (Ord a, ZeroTestable.C a, Units.C a) => PID.C (T a) where+ {-# INLINE gcd #-} gcd = euclid modCent+ {-# INLINE extendedGCD #-} extendedGCD = extendedEuclid divModCent -defltDiv :: (Field.C a) => T a -> T a -> T a-defltDiv (Cons x y) z'@(Cons x' y') =- let d = defltMagnitudeSqr z'+{-# INLINE [1] divide #-}+divide :: (Field.C a) => T a -> T a -> T a+divide (Cons x y) z'@(Cons x' y') =+ let d = magnitudeSqr z' in Cons ((x*x'+y*y') / d) ((y*x'-x*y') / d) -- | Special implementation of @(\/)@ for floating point numbers -- which prevent intermediate overflows.-floatDiv :: (P.RealFloat a, Field.C a) => T a -> T a -> T a-floatDiv (Cons x y) (Cons x' y') =+{-# INLINE floatDivide #-}+floatDivide :: (P.RealFloat a, Field.C a) => T a -> T a -> T a+floatDivide (Cons x y) (Cons x' y') = let k = - max (P.exponent x') (P.exponent y') x'' = P.scaleFloat k x' y'' = P.scaleFloat k y' d = x'*x'' + y'*y'' in Cons ((x*x''+y*y'') / d) ((y*x''-x*y'') / d) -{-|- In order to have an efficient implementation- for both complex floats and exact complex numbers,- we define the intermediate class Complex.Divisible- which in fact implements the complex division.- This way we avoid overlapping and undecidable instances.- In most cases it should suffice to define- an instance of Complex.Divisible with no method implementation- for each instance of Fractional.--}-class (Field.C a) => Divisible a where- divide :: T a -> T a -> T a- divide = defltDiv+{-# RULES+ "Complex.divide :: Double"+ divide = floatDivide :: T Double -> T Double -> T Double; -instance Divisible Float where- divide = floatDiv+ "Complex.divide :: Float"+ divide = floatDivide :: T Float -> T Float -> T Float;+ #-} -instance Divisible Double where- divide = floatDiv -instance (PID.C a) => Divisible (Ratio.T a) -instance (Divisible a) => Field.C (T a) where- (/) = divide- fromRational' = fromReal . fromRational'+instance (Field.C a) => Field.C (T a) where+ {-# INLINE (/) #-}+ (/) = divide+ {-# INLINE fromRational' #-}+ fromRational' = fromReal . fromRational' {-| We like to build the Complex Algebraic instance@@ -394,7 +421,8 @@ power :: Rational -> T a -> T a -defltPow :: (Polar a, RealTrans.C a) =>+{-# INLINE defltPow #-}+defltPow :: (RealTrans.C a) => Rational -> T a -> T a defltPow r x = let (mag,arg) = toPolar x@@ -403,14 +431,17 @@ instance Power Float where+ {-# INLINE power #-} power = defltPow instance Power Double where+ {-# INLINE power #-} power = defltPow -instance (Polar a, Real.C a, Algebraic.C a, Divisible a, Power a) =>+instance (Real.C a, Algebraic.C a, Power a) => Algebraic.C (T a) where+ {-# INLINE sqrt #-} sqrt z@(Cons x y) = if z == zero then zero else@@ -418,27 +449,38 @@ u' = sqrt ((magnitude z + abs x) / 2) (u,v) = if x < 0 then (v',u') else (u',v') in Cons u (if y < 0 then -v else v)+ {-# INLINE (^/) #-} (^/) = flip power -instance (Polar a, Real.C a, RealTrans.C a, Divisible a, Power a) =>+instance (Real.C a, RealTrans.C a, Power a) => Trans.C (T a) where- {-# SPECIALISE instance Trans.C (T Float) #-}- {-# SPECIALISE instance Trans.C (T Double) #-}+ {-# SPECULATE instance Trans.C (T Float) #-}+ {-# SPECULATE instance Trans.C (T Double) #-}+ {-# INLINE pi #-} pi = fromReal pi+ {-# INLINE exp #-} exp (Cons x y) = scale (exp x) (cis y)+ {-# INLINE log #-} log z = let (m,p) = toPolar z in Cons (log m) p -- use defaults for tan, tanh + {-# INLINE sin #-} sin (Cons x y) = Cons (sin x * cosh y) ( cos x * sinh y)+ {-# INLINE cos #-} cos (Cons x y) = Cons (cos x * cosh y) (- sin x * sinh y) + {-# INLINE sinh #-} sinh (Cons x y) = Cons (cos y * sinh x) (sin y * cosh x)+ {-# INLINE cosh #-} cosh (Cons x y) = Cons (cos y * cosh x) (sin y * sinh x) + {-# INLINE asin #-} asin z = orthoRight (log (orthoLeft z + sqrt (1 - z^2)))+ {-# INLINE acos #-} acos z = orthoRight (log (z + orthoLeft (sqrt (1 - z^2))))+ {-# INLINE atan #-} atan z@(Cons x y) = orthoRight (log (Cons (1-y) x / sqrt (1+z^2))) {- use the default implementation@@ -450,18 +492,27 @@ {- * legacy instances -} +{-# INLINE legacyInstance #-} legacyInstance :: a legacyInstance = error "legacy Ring.C instance for simple input of numeric literals" instance (Ring.C a, Eq a, Show a) => P.Num (T a) where+ {-# INLINE fromInteger #-} fromInteger = fromReal . fromInteger+ {-# INLINE negate #-} negate = negate -- for unary minus+ {-# INLINE (+) #-} (+) = legacyInstance+ {-# INLINE (*) #-} (*) = legacyInstance+ {-# INLINE abs #-} abs = legacyInstance+ {-# INLINE signum #-} signum = legacyInstance -instance (Ring.C a, Eq a, Show a, Divisible a) => P.Fractional (T a) where+instance (Field.C a, Eq a, Show a) => P.Fractional (T a) where+ {-# INLINE fromRational #-} fromRational = fromRational+ {-# INLINE (/) #-} (/) = legacyInstance
src/Number/DimensionTerm.hs view
@@ -1,4 +1,5 @@ {-# OPTIONS -fglasgow-exts #-}+-- glasgow-exts for multi-parameter type class instances {- | Copyright : (c) Henning Thielemann 2008 License : GPL@@ -23,11 +24,14 @@ import qualified Algebra.Ring as Ring import qualified Algebra.Additive as Additive +import Algebra.Field ((/), fromRational', )+import Algebra.Ring ((*), one, fromInteger, ) import Algebra.Additive ((+), (-), zero, negate, )-import Algebra.Module ((*>), )+import Algebra.Module ((*>), ) import System.Random (Random, randomR, random) +import NumericPrelude.Tuple (mapFst, ) import PreludeBase import Prelude () @@ -70,14 +74,21 @@ instance (Dim.C u, Module.C a b) => Module.C a (T u b) where a *> (Cons b) = Cons (a *> b) +instance (Dim.IsScalar u, Ring.C a) => Ring.C (T u a) where+ one = Cons one+ (Cons a) * (Cons b) = Cons (a*b)+ fromInteger a = Cons (fromInteger a)++instance (Dim.IsScalar u, Field.C a) => Field.C (T u a) where+ (Cons a) / (Cons b) = Cons (a/b)+ recip (Cons a) = Cons (Field.recip a)+ fromRational' a = Cons (fromRational' a)+ instance (OccScalar.C a b) => OccScalar.C a (Scalar b) where toScalar = OccScalar.toScalar . toNumber toMaybeScalar = OccScalar.toMaybeScalar . toNumber fromScalar = fromNumber . OccScalar.fromScalar -mapFst :: (a -> c) -> (a,b) -> (c,b)-mapFst f ~(x,y) = (f x, y)- instance (Dim.C u, Random a) => Random (T u a) where randomR (Cons l, Cons u) = mapFst Cons . randomR (l,u) random = mapFst Cons . random@@ -169,6 +180,9 @@ type Frequency a = T Dim.Frequency a type Voltage a = T Dim.Voltage a ++scalar :: a -> Scalar a+scalar = fromNumber length :: a -> Length a length = Cons
src/Number/NonNegative.hs view
@@ -17,7 +17,7 @@ Ratio, Rational) where import Numeric.NonNegative.Wrapper- (T, fromNumberUnsafe, toNumber)+ (T, fromNumberUnsafe, toNumber, ) import qualified Numeric.NonNegative.Wrapper as NonNegW import qualified Algebra.NonNegative as NonNeg@@ -41,7 +41,8 @@ import qualified Prelude as P import PreludeBase-import NumericPrelude hiding (Int, Integer, Float, Double, Rational)+import NumericPrelude.Tuple (mapSnd, mapPair, )+import NumericPrelude hiding (Int, Integer, Float, Double, Rational, ) {- |@@ -193,23 +194,7 @@ type Ratio a = T (R.T a)-type Rational = T P.Rational----- auxiliary functions--mapPair :: (a -> c, b -> d) -> (a,b) -> (c,d)-mapPair ~(f,g) ~(x,y) = (f x, g y)--{--mapFst :: (a -> c) -> (a,b) -> (c,b)-mapFst f ~(x,y) = (f x, y)--}--mapSnd :: (b -> d) -> (a,b) -> (a,d)-mapSnd g ~(x,y) = (x, g y)--+type Rational = T R.Rational {- legacy instances already defined in non-negative package -}
+ src/Number/NonNegativeChunky.hs view
@@ -0,0 +1,266 @@+{- |+Copyright : (c) Henning Thielemann 2007++Maintainer : haskell@henning-thielemann.de+Stability : stable+Portability : Haskell 98++A lazy number type, which is a generalization of lazy Peano numbers.+Comparisons can be made lazy and+thus computations are possible which are impossible with strict number types,+e.g. you can compute @let y = min (1+y) 2 in y@.+You can even work with infinite values.+However, depending on the granularity,+the memory consumption is higher than that for strict number types.+This number type is of interest for the merge operation of event lists,+which allows for co-recursive merges.+-}+module Number.NonNegativeChunky+ (T, fromChunks, fromNumber, toNumber, fromChunky98, toChunky98,+ normalize, isNull, isPositive) where++import qualified Numeric.NonNegative.ChunkyPrivate as Chunky98+import qualified Numeric.NonNegative.Class as NonNeg98++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+import qualified Algebra.RealIntegral as RealIntegral+import qualified Algebra.ZeroTestable as ZeroTestable+import Algebra.ZeroTestable (isZero, )++import Control.Monad (liftM, liftM2, )++import Test.QuickCheck (Arbitrary(..))++import NumericPrelude+import NumericPrelude.Tuple (mapFst, mapPair, )+import PreludeBase+import qualified Prelude as P98+++{- |+A chunky non-negative number is a list of non-negative numbers.+It represents the sum of the list elements.+It is possible to represent a finite number with infinitely many chunks+by using an infinite number of zeros.++Note the following problems:++Addition is commutative only for finite representations.+E.g. @let y = min (1+y) 2 in y@ is defined,+@let y = min (y+1) 2 in y@ is not.++The type is equivalent to 'Numeric.NonNegative.Chunky'.+-}+newtype T a = Cons {decons :: [a]}+++fromChunks :: NonNeg.C a => [a] -> T a+fromChunks = Cons++toChunks :: NonNeg.C a => T a -> [a]+toChunks = decons++fromChunky98 :: (NonNeg.C a, NonNeg98.C a) => Chunky98.T a -> T a+fromChunky98 = fromChunks . Chunky98.toChunksUnsafe++toChunky98 :: (NonNeg.C a, NonNeg98.C a) => T a -> Chunky98.T a+toChunky98 = Chunky98.fromChunks . toChunks++fromNumber :: NonNeg.C a => a -> T a+fromNumber = fromChunks . (:[])++toNumber :: NonNeg.C a => T a -> a+toNumber = sum . toChunks++++lift2 :: NonNeg.C a => ([a] -> [a] -> [a]) -> (T a -> T a -> T a)+lift2 f x y =+ fromChunks $ f (toChunks x) (toChunks y)++{- |+Remove zero chunks.+-}+normalize :: NonNeg.C a => T a -> T a+normalize = fromChunks . filter (>zero) . toChunks++isNullList :: NonNeg.C a => [a] -> Bool+isNullList = null . filter (>zero)++isNull :: NonNeg.C a => T a -> Bool+isNull = isNullList . toChunks+ -- null . toChunks . normalize++isPositive :: NonNeg.C a => T a -> Bool+isPositive = not . isNull++++{-+normalizeZT :: ZeroTestable.C a => T a -> T a+normalizeZT = fromChunks . filter (not . isZero) . toChunks+-}++isNullListZT :: ZeroTestable.C a => [a] -> Bool+isNullListZT = null . filter (not . isZero)++isNullZT :: ZeroTestable.C a => T a -> Bool+isNullZT = isNullListZT . decons+ -- null . toChunks . normalize+{-+isPositiveZT :: ZeroTestable.C a => T a -> Bool+isPositiveZT = not . isNull+-}+++check :: String -> Bool -> a -> a+check funcName b x =+ if b+ then x+ 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<y ==> 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)+glue (x:xs) (y:ys) =+ let (z,(zs,rs,b)) =+ case compare x y of+ LT -> (x, glue xs ((y-x):ys))+ GT -> (y, glue ((x-y):xs) ys)+ EQ -> (x, glue xs ys)+ in (z:zs,rs,b)++equalList :: (NonNeg.C a) => [a] -> [a] -> Bool+equalList x y =+ let (_,r,_) = glue x y+ in isNullList r++compareList :: (NonNeg.C a) => [a] -> [a] -> Ordering+compareList x y =+ let (_,r,b) = glue x y+ in if isNullList r+ then EQ+ else if b then LT else GT++minList :: (NonNeg.C a) => [a] -> [a] -> [a]+minList x y =+ let (z,_,_) = glue x y in z++maxList :: (NonNeg.C a) => [a] -> [a] -> [a]+maxList x y =+ let (z,r,_) = glue x y in z++r+++instance (NonNeg.C a) => Eq (T a) where+ (Cons x) == (Cons y) = equalList x y++instance (NonNeg.C a) => Ord (T a) where+ compare (Cons x) (Cons y) = compareList x y+ min = lift2 minList+ max = lift2 maxList+++instance (NonNeg.C a) => NonNeg.C (T a) where+ (-|) =+ lift2 (\x w ->+ let sub _ [] = []+ sub z (y:ys) =+ if z<y then (y-z):ys else sub (z-y) ys+ in foldr sub x w)++instance (ZeroTestable.C a) => ZeroTestable.C (T a) where+ isZero = isNullZT++instance (NonNeg.C a) => Additive.C (T a) where+ zero = Cons []+ (+) = lift2 (++)+ x - y =+ let (_,d,b) = glue (toChunks x) (toChunks y)+ d' = fromChunks d+ in check "-" (not b || isNull d') d'+ negate x = check "negate" (isNull x) x+{-+ x0 - y0 =+ let d' = lift2 (\x y -> let (_,d,b) = glue x y in d) x0 y0+ in check "-" (not b || isNull d') d'+-}++instance (Ring.C a, NonNeg.C a) => Ring.C (T a) where+ one = fromNumber one+ (*) = lift2 (liftM2 (*))+ fromInteger = fromNumber . fromInteger++instance (Ring.C a, ZeroTestable.C a, NonNeg.C a) => Real.C (T a) where+ abs = id+ signum = fromNumber . (\b -> if b then one else zero) . isPositive++instance (ToInteger.C a, NonNeg.C a) => ToInteger.C (T a) where+ toInteger = sum . map toInteger . toChunks++instance (ToRational.C a, NonNeg.C a) => ToRational.C (T a) where+ toRational = sum . map toRational . toChunks+++instance (RealIntegral.C a, NonNeg.C a) => RealIntegral.C (T a) where+ quot = div+ rem = mod+ quotRem = divMod++instance (Ord a, Integral.C a, NonNeg.C a) => Integral.C (T a) where+ divMod x0 y0 =+ let y = toChunks y0+ recurse x =+ let (r,d,b) = glue y x+ in if not b+ then ([], r)+ else mapFst (one:) (recurse d)+ in mapPair+ (fromChunks, fromChunks)+ (recurse (toChunks x0))+++++instance (Show a) => Show (T a) where+ showsPrec p x =+ showParen (p>10)+ (showString "Chunky.fromChunks " . showsPrec 10 (decons x))+++instance (NonNeg.C a, Arbitrary a) => Arbitrary (T a) where+ arbitrary = liftM Cons arbitrary+ coarbitrary = undefined++++{- * legacy instances -}++legacyInstance :: a+legacyInstance =+ error "legacy Ring.C instance for simple input of numeric literals"++instance (Ring.C a, Eq a, Show a, NonNeg.C a) => P98.Num (T a) where+ fromInteger = fromNumber . fromInteger+ negate = Additive.negate -- for unary minus+ (+) = legacyInstance+ (*) = legacyInstance+ abs = legacyInstance+ signum = legacyInstance++instance (Field.C a, Eq a, Show a, NonNeg.C a) => P98.Fractional (T a) where+ fromRational = fromNumber . fromRational+ (/) = legacyInstance
src/Number/Positional/Check.hs view
@@ -205,15 +205,8 @@ -- for complex numbers -instance Complex.Polar T where- magnitude = Complex.defltMagnitude- phase = Complex.defltPhase- instance Complex.Power T where power = Complex.defltPow--instance Complex.Divisible T where- divide = Complex.defltDiv
src/Number/Ratio.hs view
@@ -35,13 +35,14 @@ import Algebra.Units (stdUnitInv, stdAssociate) import Algebra.IntegralDomain (div, divMod) import Algebra.Ring (one, (*), fromInteger)-import Algebra.Additive (zero, (+), negate)+import Algebra.Additive (zero, (+), (-), negate) import Algebra.ZeroTestable (isZero) -- import NumericPrelude.Monad(untilM) import Control.Monad(liftM, liftM2) import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary))+import System.Random (Random(..), RandomGen, ) import qualified Data.Ratio as Ratio98 @@ -158,8 +159,34 @@ (liftM (\x -> if isZero x then one else x) arbitrary) coarbitrary = undefined +{-+This instance may not be appropriate for mathematical objects other than numbers.+If we encounter such a type of object+we should define an intermediate class+which provides the necessary functions.+I should remark that methods of Random like 'randomR'+cannot sensibly be defined for ratios of polynomials.+-}+instance (Random a, PID.C a, ZeroTestable.C a) => Random (T a) where+ random g0 =+ let (numer, g1) = random g0+ (denom, g2) = random g1+ in (numer % if isZero denom then one else denom, g2)+ randomR (lower,upper) g0 =+ let (k, g1) = randomR01 g0+ in (lower + k*(upper-lower), g1) +randomR01 ::+ (Random a, PID.C a, RandomGen g) =>+ g -> (T a, g)+randomR01 g0 =+ let (denom0, g1) = random g0+ denom = if isZero denom0 then one else denom0+ (numer, g2) = randomR (zero,denom) g1+ in (numer % denom, g2)++ -- * Legacy Instances @@ -169,21 +196,24 @@ toRational98 x = numerator x Ratio98.% denominator x -legacyInstance :: a-legacyInstance = error "legacy Ring instance for simple input of numeric literals"+legacyInstance :: String -> a+legacyInstance op =+ error ("Ratio." ++ op ++ ": legacy Ring instance for simple input of numeric literals") -instance (P.Num a, PID.C a) => P.Num (T a) where+-- instance (P.Num a, PID.C a) => P.Num (T a) where+instance (P.Num a, PID.C a, Real.C a) => P.Num (T a) where fromInteger n = P.fromInteger n % 1 negate = negate -- for unary minus- (+) = legacyInstance- (*) = legacyInstance- abs = legacyInstance- signum = legacyInstance+ (+) = legacyInstance "(+)"+ (*) = legacyInstance "(*)"+ abs = Real.abs -- needed for Arbitrary instance of NonNegative.Ratio+ signum = legacyInstance "signum" -instance (P.Num a, PID.C a) => P.Fractional (T a) where+-- instance (P.Num a, PID.C a) => P.Fractional (T a) where+instance (P.Num a, PID.C a, Real.C a) => P.Fractional (T a) where -- fromRational = Field.fromRational fromRational x = fromInteger (Ratio98.numerator x) :% fromInteger (Ratio98.denominator x)- (/) = legacyInstance+ (/) = legacyInstance "(/)"
src/NumericPrelude.hs view
@@ -7,7 +7,7 @@ {- Field -} (/), recip, fromRational', (^-), fieldPower, fromRational, {- Algebraic -} (^/), sqrt, {- Transcendental -}- pi, exp, log, logBase, (**), sin, cos, tan,+ pi, exp, log, logBase, (**), (^?), sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, {- Real -} abs, signum, {- RealIntegral -} quot, rem, quotRem,
src/NumericPrelude/Condition.hs view
@@ -5,6 +5,7 @@ {- | Returns 'Just' if the precondition is fulfilled. -}+{-# INLINE toMaybe #-} toMaybe :: Bool -> a -> Maybe a toMaybe False _ = Nothing toMaybe True x = Just x@@ -13,6 +14,7 @@ A purely functional implementation of @if@. Very useful in connection with 'zipWith3'. -}+{-# INLINE if' #-} if' :: Bool {-^ condition -} -> a {-^ then -}@@ -29,6 +31,7 @@ > [(x>0, "positive"), > (x<0, "negative")] -}+{-# INLINE select #-} select :: a -> [(Bool, a)] -> a select = foldr (uncurry if') @@ -41,6 +44,7 @@ Funnily because of the ordering of 'Bool' it holds @implies == (<=)@. -}+{-# INLINE implies #-} implies :: Bool -> Bool -> Bool implies prerequisite conclusion = not prerequisite || conclusion
src/NumericPrelude/List.hs view
@@ -44,9 +44,16 @@ {- * Use lists as counters -} {- | Make a list as long as another one -}+{-# INLINE takeMatch #-} takeMatch :: [b] -> [a] -> [a] takeMatch = flip (zipWith const) +{-# INLINE dropMatch #-}+dropMatch :: [b] -> [a] -> [a]+dropMatch (_:xs) (_:ys) = dropMatch xs ys+dropMatch _ ys = ys++{-# INLINE splitAtMatch #-} splitAtMatch :: [b] -> [a] -> ([a],[a]) splitAtMatch (_:ns) (x:xs) = let (as,bs) = splitAtMatch ns xs@@ -54,6 +61,7 @@ splitAtMatch _ [] = ([],[]) splitAtMatch [] xs = ([],xs) +{-# INLINE replicateMatch #-} replicateMatch :: [a] -> b -> [b] replicateMatch xs y = takeMatch xs (repeat y)@@ -63,6 +71,7 @@ For finite lists it is equivalent to (compare (length xs) (length ys)) but more efficient. -}+{-# INLINE compareLength #-} compareLength :: [a] -> [b] -> Ordering compareLength (_:xs) (_:ys) = compareLength xs ys compareLength [] [] = EQ@@ -76,6 +85,7 @@ {- * Zip lists -} {- | zip two lists using an arbitrary function, the shorter list is padded -}+{-# INLINE zipWithPad #-} zipWithPad :: a {-^ padding value -} -> (a -> a -> b) {-^ function applied to corresponding elements of the lists -} -> [a]@@ -87,6 +97,7 @@ aux (x:xs) (y:ys) = f x y : aux xs ys in aux +{-# INLINE zipWithOverlap #-} zipWithOverlap :: (a -> c) -> (b -> c) -> (a -> b -> c) -> [a] -> [b] -> [c] zipWithOverlap fa fb fab = let aux (x:xs) (y:ys) = fab x y : aux xs ys@@ -110,6 +121,7 @@ aux _ _ = error "zipWith: lists must have the same length" in aux +{-# INLINE zipNeighborsWith #-} zipNeighborsWith :: (a -> a -> a) -> [a] -> [a] zipNeighborsWith f xs = zipWith f xs (drop 1 xs) @@ -182,6 +194,7 @@ {- * Various -} +{-# INLINE partitionMaybe #-} partitionMaybe :: (a -> Maybe b) -> [a] -> ([b], [a]) partitionMaybe f = foldr (\x ~(y,z) -> case f x of@@ -194,6 +207,7 @@ if the last element is accessed after the initial ones, because it avoids memoizing list. -}+{-# INLINE splitLast #-} splitLast :: [a] -> ([a], a) splitLast [] = error "splitLast: empty list" splitLast [x] = ([], x)@@ -209,6 +223,7 @@ In contrast to 'reverse . dropWhile p . reverse' this works for infinite lists, too. -}+{-# INLINE dropWhileRev #-} dropWhileRev :: (a -> Bool) -> [a] -> [a] dropWhileRev p = foldr (\x xs -> if p x && null xs then [] else x:xs) []@@ -218,6 +233,7 @@ Apply a function to the last element of a list. If the list is empty, nothing changes. -}+{-# INLINE mapLast #-} mapLast :: (a -> a) -> [a] -> [a] mapLast f = let recurse [] = [] -- behaviour as needed in powerBasis@@ -226,10 +242,12 @@ recurse (x:xs) = x : recurse xs in recurse +{-# INLINE padLeft #-} padLeft :: a -> Int -> [a] -> [a] padLeft c n xs = replicate (n - length xs) c ++ xs +{-# INLINE padRight #-} padRight :: a -> Int -> [a] -> [a] padRight c n xs = xs ++ replicate (n - length xs) c @@ -244,6 +262,8 @@ but applies "op" O(log n) times and works for large n. -} +{-# INLINE reduceRepeated #-}+{-# INLINE reduceRepeatedSlow #-} reduceRepeated, reduceRepeatedSlow :: (a -> a -> a) -> a -> a -> Integer -> a reduceRepeated _ a0 _ 0 = a0
src/NumericPrelude/Monad.hs view
@@ -1,6 +1,7 @@ module NumericPrelude.Monad where {- | repeat action until result fulfills condition -}+{-# INLINE untilM #-} untilM :: (Monad m) => (a -> Bool) -> m a -> m a untilM p m = do x <- m
src/NumericPrelude/Text.hs view
@@ -3,6 +3,7 @@ {-* Formatting and parsing. -} {-| Show a value using an infix operator. -}+{-# INLINE showsInfixPrec #-} showsInfixPrec :: (Show a, Show b) => String -> Int -> Int -> a -> b -> ShowS showsInfixPrec opStr opPrec prec x y =@@ -13,6 +14,7 @@ showsPrec opPrec y) {-| Parse a string containing an infix operator. -}+{-# INLINE readsInfixPrec #-} readsInfixPrec :: (Read a, Read b) => String -> Int -> Int -> (a -> b -> c) -> ReadS c readsInfixPrec opStr opPrec prec cons =
+ src/NumericPrelude/Tuple.hs view
@@ -0,0 +1,13 @@+module NumericPrelude.Tuple where++{-# INLINE mapPair #-}+mapPair :: (a -> c, b -> d) -> (a,b) -> (c,d)+mapPair ~(f,g) ~(x,y) = (f x, g y)++{-# INLINE mapFst #-}+mapFst :: (a -> c) -> (a,b) -> (c,b)+mapFst f ~(x,y) = (f x, y)++{-# INLINE mapSnd #-}+mapSnd :: (b -> d) -> (a,b) -> (a,d)+mapSnd g ~(x,y) = (x, g y)