diff --git a/Makefile b/Makefile --- a/Makefile +++ b/Makefile @@ -54,7 +54,7 @@ ./testsuite ghci: - $(HCI) -Wall -i:src +RTS -M256m -c30 -RTS test/Test.hs + $(HCI) -Wall -i:src:test +RTS -M256m -c30 -RTS test/Test.hs build: -mkdir $(OBJECT_DIR) diff --git a/numeric-prelude.cabal b/numeric-prelude.cabal --- a/numeric-prelude.cabal +++ b/numeric-prelude.cabal @@ -1,5 +1,5 @@ Name: numeric-prelude -Version: 0.0.5 +Version: 0.1 License: GPL License-File: LICENSE Author: Dylan Thurston , Henning Thielemann , Mikael Johansson @@ -8,12 +8,13 @@ Package-URL: http://darcs.haskell.org/numericprelude/ Category: Math Stability: Experimental +Tested-With: GHC==6.4.1, GHC==6.8.2 Synopsis: An experimental alternative hierarchy of numeric type classes Description: Revisiting the Numeric Classes . - The Prelude for Haskell 98 offers a well-considered set of numeric - classes which cover the standard numeric types + The Prelude for Haskell 98 offers a well-considered set of numeric classes + which covers the standard numeric types ('Integer', 'Int', 'Rational', 'Float', 'Double', 'Complex') quite well. But they offer limited extensibility and have a few other flaws. In this proposal we will revisit these classes, addressing the following concerns: @@ -64,7 +65,7 @@ . For inexact number types like floating point types, thus these laws should be interpreted as guidelines rather than absolute rules. - In particular, the compiler is not allowed to use them. + In particular, the compiler is not allowed to use them for optimization. Unless stated otherwise, default definitions should also be taken as laws. . Thanks to Brian Boutel, Joe English, William Lee Irwin II, Marcin @@ -76,7 +77,7 @@ (a total ordering). This is not addressed here. . - * In some cases, this hierarchy may not be fine-grained enough. + * In some cases, this hierarchy may not yet be fine-grained enough. For instance, time spans (\"5 minutes\") can be added to times (\"12:34\"), but two times are not addable. (\"12:34 + 8:23\") As it stands, @@ -130,10 +131,15 @@ default: False Library - Build-Depends: parsec >= 1, HUnit >=1 && <2, QuickCheck >=1 && <2 - Build-Depends: non-negative>=0.0.2 && <0.1 + Build-Depends: parsec >=1 && <3, HUnit >=1 && <2, QuickCheck >=1 && <2 + Build-Depends: non-negative >=0.0.2 && <0.1 + Build-Depends: utility-ht >=0.0.4 && <0.1 If flag(splitBase) - Build-Depends: base >= 2, array, containers, random + Build-Depends: + base >= 2 && <5, + array >=0.1 && <0.3, + containers >=0.1 && <0.3, + random >=1.0 && <1.1 Else Build-Depends: base >= 1.0 && < 2 @@ -217,12 +223,12 @@ Number.Physical.Read Number.Physical.Show NumericPrelude - NumericPrelude.Condition - NumericPrelude.List - NumericPrelude.Monad - NumericPrelude.Text - NumericPrelude.Tuple PreludeBase + Other-modules: + NumericPrelude.List + MathObj.Gaussian.Variance + MathObj.Gaussian.Bell + MathObj.Gaussian.Polynomial Executable test Hs-Source-Dirs: src, test @@ -235,10 +241,11 @@ GHC-Options: -Wall Other-modules: Test.NumericPrelude.Utility - Test.NumericPrelude.List Test.MathObj.PartialFraction Test.MathObj.Polynomial Test.MathObj.PowerSeries + Test.MathObj.Gaussian.Variance + Test.MathObj.Gaussian.Bell Main-Is: Test/Run.hs If !flag(buildTests) Buildable: False diff --git a/src/Algebra/Additive.hs b/src/Algebra/Additive.hs --- a/src/Algebra/Additive.hs +++ b/src/Algebra/Additive.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Additive ( {- * Class -} C, diff --git a/src/Algebra/Algebraic.hs b/src/Algebra/Algebraic.hs --- a/src/Algebra/Algebraic.hs +++ b/src/Algebra/Algebraic.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Algebraic where import qualified Algebra.Field as Field diff --git a/src/Algebra/Differential.hs b/src/Algebra/Differential.hs --- a/src/Algebra/Differential.hs +++ b/src/Algebra/Differential.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Differential where import qualified Algebra.Ring as Ring diff --git a/src/Algebra/DimensionTerm.hs b/src/Algebra/DimensionTerm.hs --- a/src/Algebra/DimensionTerm.hs +++ b/src/Algebra/DimensionTerm.hs @@ -143,9 +143,11 @@ -} class C dim => IsScalar dim where toScalar :: dim -> Scalar + fromScalar :: Scalar -> dim instance IsScalar Scalar where toScalar = id + fromScalar = id {- ** Basis dimensions -} diff --git a/src/Algebra/DivisibleSpace.hs b/src/Algebra/DivisibleSpace.hs --- a/src/Algebra/DivisibleSpace.hs +++ b/src/Algebra/DivisibleSpace.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} module Algebra.DivisibleSpace where import qualified Prelude import qualified Algebra.VectorSpace as VectorSpace diff --git a/src/Algebra/Field.hs b/src/Algebra/Field.hs --- a/src/Algebra/Field.hs +++ b/src/Algebra/Field.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Field ( {- * Class -} C, @@ -26,8 +26,6 @@ import Algebra.Ring ((*), (^), one, fromInteger) import Algebra.Additive (zero, negate) import Algebra.ZeroTestable (isZero) - --- import NumericPrelude.List(reduceRepeated) import PreludeBase import Prelude (Integer, Float, Double) diff --git a/src/Algebra/IntegralDomain.hs b/src/Algebra/IntegralDomain.hs --- a/src/Algebra/IntegralDomain.hs +++ b/src/Algebra/IntegralDomain.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.IntegralDomain ( {- * Class -} C, @@ -35,8 +35,8 @@ import Algebra.Additive (zero, (+), (-)) import Algebra.ZeroTestable (isZero) -import NumericPrelude.Condition(implies) -import Data.List(mapAccumL) +import Data.Bool.HT (implies, ) +import Data.List (mapAccumL, ) import Test.QuickCheck ((==>), Property) diff --git a/src/Algebra/Lattice.hs b/src/Algebra/Lattice.hs --- a/src/Algebra/Lattice.hs +++ b/src/Algebra/Lattice.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Lattice ( C(up, dn) , max, min, abs diff --git a/src/Algebra/Module.hs b/src/Algebra/Module.hs --- a/src/Algebra/Module.hs +++ b/src/Algebra/Module.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Dylan Thurston, Henning Thielemann 2004-2005 @@ -20,11 +22,11 @@ import qualified Algebra.Laws as Laws -import Algebra.Ring ((*), fromInteger) -import Algebra.Additive ((+), zero) +import Algebra.Ring ((*), fromInteger, ) +import Algebra.Additive ((+), zero, ) -import NumericPrelude.List (reduceRepeated) -import Data.List (map, zipWith, foldl) +import Data.Function.HT (powerAssociative, ) +import Data.List (map, zipWith, foldl, ) import Prelude((.), Eq, Bool, Int, Integer, Float, Double) -- import qualified Prelude as P @@ -116,7 +118,7 @@ {-# INLINE integerMultiply #-} integerMultiply :: (ToInteger.C a, Additive.C b) => a -> b -> b integerMultiply a b = - reduceRepeated (+) zero b (ToInteger.toInteger a) + powerAssociative (+) zero b (ToInteger.toInteger a) {- * Properties -} diff --git a/src/Algebra/ModuleBasis.hs b/src/Algebra/ModuleBasis.hs --- a/src/Algebra/ModuleBasis.hs +++ b/src/Algebra/ModuleBasis.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Maintainer : numericprelude@henning-thielemann.de Stability : provisional @@ -9,8 +11,9 @@ module Algebra.ModuleBasis where -import Number.Ratio (Rational) +import qualified Number.Ratio as Ratio +import qualified Algebra.PrincipalIdealDomain as PID import qualified Algebra.Module as Module import qualified Algebra.Additive as Additive import Algebra.Ring (one, fromInteger) @@ -59,7 +62,7 @@ flatten = (:[]) dimension _ _ = 1 -instance C Rational Rational where +instance (PID.C a) => C (Ratio.T a) (Ratio.T a) where basis _ = [one] flatten = (:[]) dimension _ _ = 1 diff --git a/src/Algebra/Monoid.hs b/src/Algebra/Monoid.hs --- a/src/Algebra/Monoid.hs +++ b/src/Algebra/Monoid.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Mikael Johansson 2006 Maintainer : mik@math.uni-jena.de diff --git a/src/Algebra/NormedSpace/Euclidean.hs b/src/Algebra/NormedSpace/Euclidean.hs --- a/src/Algebra/NormedSpace/Euclidean.hs +++ b/src/Algebra/NormedSpace/Euclidean.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2005 diff --git a/src/Algebra/NormedSpace/Maximum.hs b/src/Algebra/NormedSpace/Maximum.hs --- a/src/Algebra/NormedSpace/Maximum.hs +++ b/src/Algebra/NormedSpace/Maximum.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2005 diff --git a/src/Algebra/NormedSpace/Sum.hs b/src/Algebra/NormedSpace/Sum.hs --- a/src/Algebra/NormedSpace/Sum.hs +++ b/src/Algebra/NormedSpace/Sum.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2005 diff --git a/src/Algebra/OccasionallyScalar.hs b/src/Algebra/OccasionallyScalar.hs --- a/src/Algebra/OccasionallyScalar.hs +++ b/src/Algebra/OccasionallyScalar.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fglasgow-exts -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | diff --git a/src/Algebra/PrincipalIdealDomain.hs b/src/Algebra/PrincipalIdealDomain.hs --- a/src/Algebra/PrincipalIdealDomain.hs +++ b/src/Algebra/PrincipalIdealDomain.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.PrincipalIdealDomain ( {- * Class -} C, @@ -50,7 +50,7 @@ import Algebra.Additive (zero, (+), (-)) import Algebra.ZeroTestable (isZero) -import NumericPrelude.Condition (toMaybe) +import Data.Maybe.HT (toMaybe) import Control.Monad (foldM, liftM) import Data.List (mapAccumL, mapAccumR, unfoldr) diff --git a/src/Algebra/Real.hs b/src/Algebra/Real.hs --- a/src/Algebra/Real.hs +++ b/src/Algebra/Real.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Real ( C(abs, signum), ) where diff --git a/src/Algebra/RealField.hs b/src/Algebra/RealField.hs --- a/src/Algebra/RealField.hs +++ b/src/Algebra/RealField.hs @@ -1,4 +1,5 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# OPTIONS_GHC -fglasgow-exts #-} -- -fglasgow-exts for RULES module Algebra.RealField where @@ -25,7 +26,7 @@ import Data.Word (Word, Word8, Word16, Word32, Word64, ) import qualified GHC.Float as GHC -import NumericPrelude.Tuple (mapFst, ) +import Data.Tuple.HT (mapFst, ) import Prelude(Int, Integer, Float, Double) import qualified Prelude as P import PreludeBase diff --git a/src/Algebra/RealIntegral.hs b/src/Algebra/RealIntegral.hs --- a/src/Algebra/RealIntegral.hs +++ b/src/Algebra/RealIntegral.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Generally before using 'quot' and 'rem', think twice. In most cases 'divMod' and friends are the right choice, diff --git a/src/Algebra/RealTranscendental.hs b/src/Algebra/RealTranscendental.hs --- a/src/Algebra/RealTranscendental.hs +++ b/src/Algebra/RealTranscendental.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.RealTranscendental where import qualified Algebra.Transcendental as Trans @@ -9,7 +9,7 @@ import Algebra.Ring (fromInteger) import Algebra.Additive ((+), negate) -import NumericPrelude.Condition (select) +import Data.Bool.HT (select, ) import qualified Prelude as P import PreludeBase diff --git a/src/Algebra/RightModule.hs b/src/Algebra/RightModule.hs --- a/src/Algebra/RightModule.hs +++ b/src/Algebra/RightModule.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} module Algebra.RightModule where import qualified Algebra.Ring as Ring diff --git a/src/Algebra/Ring.hs b/src/Algebra/Ring.hs --- a/src/Algebra/Ring.hs +++ b/src/Algebra/Ring.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Ring ( {- * Class -} C, @@ -28,7 +28,8 @@ import Algebra.Additive(zero, (+), negate, sum) -import NumericPrelude.List(reduceRepeated, zipWithMatch) +import Data.Function.HT (powerAssociative, ) +import NumericPrelude.List (zipWithMatch, ) import Test.QuickCheck ((==>), Property) @@ -83,11 +84,11 @@ {-# INLINE fromInteger #-} fromInteger n = if n < 0 - then reduceRepeated (+) zero (negate one) (negate n) - else reduceRepeated (+) zero one n + then powerAssociative (+) zero (negate one) (negate n) + else powerAssociative (+) zero one n {-# INLINE (^) #-} a ^ n = if n >= zero - then reduceRepeated (*) one a n + then powerAssociative (*) one a n else error "(^): Illegal negative exponent" {-# INLINE one #-} one = fromInteger 1 diff --git a/src/Algebra/ToInteger.hs b/src/Algebra/ToInteger.hs --- a/src/Algebra/ToInteger.hs +++ b/src/Algebra/ToInteger.hs @@ -1,5 +1,12 @@ -{-# OPTIONS -fglasgow-exts #-} --- -fglasgow-exts for RULES +{-# OPTIONS_GHC -fglasgow-exts -fno-warn-orphans #-} +{- +-fglasgow-exts for RULES + +The orphan instance could be fixed +by making this module mutually recursive with ToRational.hs, +but that's not worth the complication. +-} + module Algebra.ToInteger where import qualified Number.Ratio as Ratio diff --git a/src/Algebra/ToRational.hs b/src/Algebra/ToRational.hs --- a/src/Algebra/ToRational.hs +++ b/src/Algebra/ToRational.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.ToRational where import qualified Algebra.Real as Real diff --git a/src/Algebra/Transcendental.hs b/src/Algebra/Transcendental.hs --- a/src/Algebra/Transcendental.hs +++ b/src/Algebra/Transcendental.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Transcendental where import qualified Algebra.Algebraic as Algebraic diff --git a/src/Algebra/Units.hs b/src/Algebra/Units.hs --- a/src/Algebra/Units.hs +++ b/src/Algebra/Units.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.Units ( {- * Class -} C, diff --git a/src/Algebra/Vector.hs b/src/Algebra/Vector.hs --- a/src/Algebra/Vector.hs +++ b/src/Algebra/Vector.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2004-2005 diff --git a/src/Algebra/VectorSpace.hs b/src/Algebra/VectorSpace.hs --- a/src/Algebra/VectorSpace.hs +++ b/src/Algebra/VectorSpace.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} module Algebra.VectorSpace where import qualified Algebra.Module as Module diff --git a/src/Algebra/ZeroTestable.hs b/src/Algebra/ZeroTestable.hs --- a/src/Algebra/ZeroTestable.hs +++ b/src/Algebra/ZeroTestable.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Algebra.ZeroTestable where import qualified Algebra.Additive as Additive diff --git a/src/MathObj/Algebra.hs b/src/MathObj/Algebra.hs --- a/src/MathObj/Algebra.hs +++ b/src/MathObj/Algebra.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Mikael Johansson 2006 Maintainer : mik@math.uni-jena.de diff --git a/src/MathObj/DiscreteMap.hs b/src/MathObj/DiscreteMap.hs --- a/src/MathObj/DiscreteMap.hs +++ b/src/MathObj/DiscreteMap.hs @@ -1,18 +1,30 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} + +{- FIXME: +Rationale for -fno-warn-orphans: + * The orphan instances can't be put into Numeric.NonNegative.Wrapper + since that's in another package. + * We had to spread the instance declarations + over the modules defining the typeclasses instantiated. + Do we want that? + * We could define the DiscreteMap as newtype. +-} + {- | -DiscreteMap is a class that unifies -Map and Array, -thus one can simply choose between +DiscreteMap was originally intended as a type class +that unifies Map and Array. +One should be able to simply choose between - Map for sparse arrays - Array for full arrays. -Ok, forget it, -the Edison package provides the class AssocX -which probably will do it. +However, the Edison package provides the class AssocX +which already exists for that purpose. -So long I use this module for some numeric instances for FiniteMaps +Currently I use this module for some numeric instances of Data.Map. -} - module MathObj.DiscreteMap where import qualified Algebra.NormedSpace.Sum as NormedSum @@ -32,7 +44,7 @@ import qualified Prelude as P import PreludeBase --- *** Should this be implemented by isZero? +-- FIXME: Should this be implemented by isZero? -- | Remove all zero values from the map. strip :: (Ord i, Eq v, Additive.C v) => Map i v -> Map i v strip = Map.filter (zero /=) diff --git a/src/MathObj/Gaussian/Bell.hs b/src/MathObj/Gaussian/Bell.hs new file mode 100644 --- /dev/null +++ b/src/MathObj/Gaussian/Bell.hs @@ -0,0 +1,256 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{- +Complex translated Gaussian bell curve +with amplitude abstracted away. +-} +module MathObj.Gaussian.Bell where + +import qualified MathObj.Polynomial as Poly +import qualified Number.Complex as Complex + +import qualified Algebra.Transcendental as Trans +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 Number.Complex ((+:), ) +import Algebra.Transcendental (pi, ) +import Algebra.Ring ((*), (^), ) +import Algebra.Additive ((+), ) + +import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, ) +import Control.Monad (liftM3, ) + +-- import Prelude (($)) +import NumericPrelude +import PreludeBase hiding (reverse, ) + + +data T a = Cons {c0, c1 :: Complex.T a, c2 :: a} + deriving (Eq, Show) + +instance (Real.C a, Arbitrary a) => Arbitrary (T a) where + arbitrary = + liftM3 + (\a b c -> Cons a b (1 + abs c)) + arbitrary arbitrary arbitrary + coarbitrary = undefined + + +constant :: Additive.C a => T a +constant = Cons zero zero zero + +{-# INLINE evaluate #-} +evaluate :: (Trans.C a) => + T a -> a -> Complex.T a +evaluate f x = + Complex.exp $ Complex.scale (-pi) $ + c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2) + +evaluateSqRt :: (Trans.C a) => + T a -> a -> Complex.T a +evaluateSqRt f x0 = + let x = sqrt pi * x0 + in Complex.exp $ negate $ + c0 f + Complex.scale x (c1 f) + Complex.fromReal (c2 f * x^2) + +exponentPolynomial :: (Additive.C a) => + T a -> Poly.T (Complex.T a) +exponentPolynomial f = + Poly.fromCoeffs [c0 f, c1 f, Complex.fromReal (c2 f)] + + +multiply :: (Additive.C a) => + T a -> T a -> T a +multiply f g = + Cons (c0 f + c0 g) (c1 f + c1 g) (c2 f + c2 g) + + +{- +let x=Cons (1+:3) (4+:5) (7::Rational); y=Cons (1+:4) (3+:2) (5::Rational) +-} +convolve :: (Field.C a) => + T a -> T a -> T a +convolve f g = + let s = c2 f + c2 g + {- + fd = f1/(2*f2) + gd = g1/(2*g2) + c = f2*g2/(f2+g2) + + c*(fd+gd) = (f1*g2+f2*g1)/(2*(f2+g2)) = b/2 + + c*(fd+gd)^2 - fd^2*f2 - gd^2*g2 + = f2*g2*(fd+gd)^2/(f2 + g2) - (fd^2*f2 + gd^2*g2) + = (f2*g2*(fd+gd)^2 - (f2+g2)*(fd^2*f2+gd^2*g2)) / (f2 + g2) + = (2*f2*g2*fd*gd - (fd^2*f2^2+gd^2*g2^2)) / (f2 + g2) + = (2*f1*g1 - (f1^2+g1^2)) / (4*(f2 + g2)) + = -(f1 - g1)^2/(4*(f2 + g2)) + -} + in Cons + (c0 f + c0 g + - Complex.scale (recip (4*s)) ((c1 f - c1 g)^2)) + (Complex.scale (c2 g / s) (c1 f) + + Complex.scale (c2 f / s) (c1 g)) + (c2 f * c2 g / s) + -- recip $ recip (c2 f) + recip (c2 g) +{- + Cons + (c0 f + c0 g) (c1 f + c1 g) + (recip $ recip (c2 f) + recip (c2 g)) +-} + +convolveByTranslation :: (Field.C a) => + T a -> T a -> T a +convolveByTranslation f0 g0 = + let fd = Complex.scale (recip (2 * c2 f0)) $ c1 f0 + gd = Complex.scale (recip (2 * c2 g0)) $ c1 g0 + f1 = translateComplex fd f0 + g1 = translateComplex gd g0 + in translateComplex (negate $ fd + gd) $ + Cons + (c0 f1 + c0 g1) zero + (recip $ recip (c2 f1) + recip (c2 g1)) + +convolveByFourier :: (Field.C a) => + T a -> T a -> T a +convolveByFourier f g = + reverse $ fourier $ multiply (fourier f) (fourier g) + +fourier :: (Field.C a) => + T a -> T a +fourier f = + let a = c0 f + b = c1 f + c = c2 f + rc = recip c + in Cons + (Complex.scale (rc/4) (-b^2) + a) + (Complex.scale rc $ Complex.quarterRight b) + rc + +fourierByTranslation :: (Field.C a) => + T a -> T a +fourierByTranslation f = + translateComplex (Complex.scale (1/2) $ Complex.quarterLeft $ c1 f) $ + Cons (c0 f) zero (recip $ c2 f) + +{- +a + b*x + c*x^2 + = c*(a/c + b/c*x + x^2) + = c*((x-b/(2*c))^2 + (4*a*c+b^2)/(4*c^2)) + = c*(x-b/(2*c))^2 + (4*a*c+b^2)/(4*c) + +fourier -> + x^2/c - i*b/c*x + (4*a*c+b^2)/(4*c) + +fourier (x -> exp(-pi*c*(x-t)^2)) + = fourier $ translate t $ shrink (sqrt c) $ x -> exp(-pi*x^2) + = modulate t $ dilate (sqrt c) $ fourier $ x -> exp(-pi*x^2) + = modulate t $ dilate (sqrt c) $ x -> exp(-pi*x^2) + = modulate t $ x -> exp(-pi*x^2/c) + = x -> exp(-pi*x^2/c) * exp(-2*pi*i*x*t) + = x -> exp(-pi*(x^2/c - 2*i*x*t)) +-} + +{- +b*x + c*x^2 + = c*(b/c*x + x^2) + = c*((x-br/(2*c))^2 + i*x*bi/c - br^2/(4*c^2)) + = c*(x-br/(2*c))^2 + i*x*bi - br^2/(4*c) + +fourier -> + (x+bi/2)^2/c - i*br/c*(x+bi/2) - br^2/(4*c) + = (1/c) * ((x+bi/2)^2 - i*br*(x+bi/2) + (br/2)^2) + = (1/c) * (x^2 - i*b*x + -(br/2)^2 + (bi/2)^2 - i*br*bi/2) + = (1/c) * (x^2 - i*b*x - (br^2-bi^2+2*br*bi*i)^2 /4) + = (1/c) * (x^2 - i*b*x - b^2 / 4) + = (1/c) * (x^2 - i*b*x + (i*b/2)^2) + = (1/c) * (x - i*b/2)^2 + +Example: + (x-b)^2 = b^2 - 2*b*x + x^2 + -> (- i*2*b*x + x^2) + + +fourier (x -> exp(-pi*(c*(x-t)^2 + 2*i*m*x))) + = fourier $ modulate m $ translate t $ shrink (sqrt c) $ x -> exp(-pi*x^2) + = translate (-m) $ modulate t $ dilate (sqrt c) $ fourier $ x -> exp(-pi*x^2) + = translate (-m) $ modulate t $ dilate (sqrt c) $ x -> exp(-pi*x^2) + = translate (-m) $ modulate t $ x -> exp(-pi*x^2/c) + = translate (-m) $ x -> exp(-pi*x^2/c) * exp(-2*pi*i*x*t) + = x -> exp(-pi*(x+m)^2/c) * exp(-2*pi*i*(x+m)*t) + = x -> exp(-pi*((x+m)^2/c - 2*i*(x+m)*t)) +-} + +{- +fourier (Cons a 0 0) = + Cons a 0 infinity + +fourier (Cons 0 0 c) = + Cons 0 0 (recip c) + +fourier (Cons 0 b 1) = + Cons 0 (i*b) 1 +-} + +translate :: Ring.C a => a -> T a -> T a +translate d f = + let a = c0 f + b = c1 f + c = c2 f + in Cons + (Complex.fromReal (c*d^2) - Complex.scale d b + a) + (Complex.fromReal (-2*c*d) + b) + c + +translateComplex :: Ring.C a => Complex.T a -> T a -> T a +translateComplex d f = + let a = c0 f + b = c1 f + c = c2 f + in Cons + (Complex.scale c (d^2) - b*d + a) + (Complex.scale (-2*c) d + b) + c + +modulate :: Ring.C a => a -> T a -> T a +modulate d f = + Cons + (c0 f) + (c1 f + (zero +: 2*d)) + (c2 f) + +turn :: Ring.C a => a -> T a -> T a +turn d f = + Cons + (c0 f + (zero +: 2*d)) + (c1 f) + (c2 f) + +reverse :: Additive.C a => T a -> T a +reverse f = + f{c1 = negate $ c1 f} + + +dilate :: Field.C a => a -> T a -> T a +dilate k f = + Cons + (c0 f) + (Complex.scale (recip k) $ c1 f) + (c2 f / k^2) + +shrink :: Ring.C a => a -> T a -> T a +shrink k f = + Cons + (c0 f) + (Complex.scale k $ c1 f) + (k^2 * c2 f) + + +{- laws +fourier (convolve f g) = fourier f * fourier g + +fourier (fourier f) = reverse f +-} diff --git a/src/MathObj/Gaussian/Polynomial.hs b/src/MathObj/Gaussian/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/MathObj/Gaussian/Polynomial.hs @@ -0,0 +1,106 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{- +Complex Gaussian bell multiplied with a polynomial. + +In order to make this free of @pi@ factors, +we have to choose @recip (sqrt pi)@ +as unit for translations and modulations, +for linear factors and in the differentiation. +-} +module MathObj.Gaussian.Polynomial where + +import qualified MathObj.Gaussian.Bell as Bell + +import qualified MathObj.Polynomial as Poly +import qualified Number.Complex as Complex + +import qualified Algebra.Differential as Differential +import qualified Algebra.Transcendental as Trans +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 Algebra.Transcendental (pi, ) +import Algebra.Ring ((*), ) +-- import Algebra.Additive ((+)) + +import NumericPrelude +import PreludeBase hiding (reverse, ) +-- import Prelude () + + +data T a = Cons {bell :: Bell.T a, polynomial :: Poly.T (Complex.T a)} + deriving (Eq, Show) + + +{-# INLINE evaluate #-} +evaluate :: (Trans.C a) => + T a -> a -> Complex.T a +evaluate f x = + Bell.evaluateSqRt (bell f) x * + Poly.evaluate (polynomial f) (Complex.fromReal $ sqrt pi * x) +{- ToDo: evaluating a complex polynomial for a real argument can be optimized -} + + +multiply :: (Ring.C a) => + T a -> T a -> T a +multiply x y = + Cons + (Bell.multiply (bell x) (bell y)) + (polynomial x * polynomial y) + +convolve :: (Field.C a) => + T a -> T a -> T a +convolve f g = + reverse $ fourier $ multiply (fourier f) (fourier g) + +reverse :: Additive.C a => T a -> T a +reverse x = + Cons + (Bell.reverse $ bell x) + (Poly.reverse $ polynomial x) + +{- +We use a Horner like scheme +in order to translate multiplications with @id@ +to differentations on the Fourier side. +Quadratic runtime. + +fourier (Cons bell (Poly.const a + Poly.shift x)) + = fourier (Cons bell (Poly.const a)) + fourier (Cons bell (Poly.shift x)) + = fourier (Cons bell (Poly.const a)) + differentiate (fourier (Cons bell x)) + +untested +-} +fourier :: (Field.C a) => + T a -> T a +fourier x = + foldr + (\c p -> + let q = differentiate p + in q{polynomial = + Poly.const c + + fmap Complex.quarterLeft (polynomial q)}) + (Cons (Bell.fourier $ bell x) zero) $ + Poly.coeffs $ polynomial x + +{- +Differentiate and divide by @sqrt pi@ in order to stay in a ring. +This way, we do not need to fiddle with pi factors. +-} +differentiate :: (Ring.C a) => T a -> T a +differentiate f = + f{polynomial = + Bell.exponentPolynomial (bell f) * polynomial f + + Differential.differentiate (polynomial f)} + +{- No Ring instance for Gaussians +instance (Ring.C a) => Differential.C (T a) where + differentiate = differentiate +-} + +{- laws +differentiate (f*g) = + (differentiate f) * g + f * (differentiate g) +-} diff --git a/src/MathObj/Gaussian/Variance.hs b/src/MathObj/Gaussian/Variance.hs new file mode 100644 --- /dev/null +++ b/src/MathObj/Gaussian/Variance.hs @@ -0,0 +1,123 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{- +Reciprocal of variance of a Gaussian bell curve. +We describe the curve only in terms of its variance +thus we represent a bell curve at the coordinate origin +neglecting its amplitude. + +We could also define the amplitude as @root 4 c@, +but then @dilate@ and @shrink@ also include an amplification. + +We could do some projective geometry in the exponent +in order to also have zero variance, +which corresponds to the dirac impulse. +-} +module MathObj.Gaussian.Variance where + +import qualified MathObj.Polynomial as Poly + +import qualified Algebra.Transcendental as Trans +import qualified Algebra.Algebraic as Algebraic +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 Algebra.Transcendental (pi, ) +import Algebra.Ring ((*), (^), ) +import Algebra.Additive ((+)) + +import Test.QuickCheck (Arbitrary, arbitrary, coarbitrary, ) + + +-- import Prelude (($)) +import NumericPrelude +import PreludeBase + + +data T a = Cons {c :: a} + deriving (Eq, Show) + +instance (Real.C a, Arbitrary a) => Arbitrary (T a) where + arbitrary = fmap (Cons . (1+) . abs) arbitrary + coarbitrary = undefined + + +constant :: Additive.C a => T a +constant = Cons zero + +{-# INLINE evaluate #-} +evaluate :: (Trans.C a) => + T a -> a -> a +evaluate f x = + exp $ (-pi * c f * x^2) + +exponentPolynomial :: (Additive.C a) => + T a -> Poly.T a +exponentPolynomial f = + Poly.fromCoeffs [zero, zero, c f] + + +norm1 :: (Algebraic.C a) => T a -> a +norm1 f = + recip $ sqrt $ c f + +norm2 :: (Algebraic.C a) => T a -> a +norm2 f = + recip $ sqrt $ sqrt $ 2 * c f + +normP :: (Trans.C a) => a -> T a -> a +normP p f = + (p * c f) ^? (- recip (2*p)) + + +variance :: (Trans.C a) => + T a -> a +variance f = + recip $ c f * 2*pi + +multiply :: (Additive.C a) => + T a -> T a -> T a +multiply f g = + Cons $ c f + c g + +{- | +> convolve x y t = +> integrate $ \s -> x s * y(t-s) +-} +convolve :: (Field.C a) => + T a -> T a -> T a +convolve f g = + Cons $ recip $ recip (c f) + recip (c g) + +{- | +> fourier x f = +> integrate $ \t -> x t * cis (-2*pi*t*f) +-} +fourier :: (Field.C a) => + T a -> T a +fourier f = + Cons $ recip $ c f +{- +fourier (t -> exp(-(a*t)^2)) +-} + +dilate :: (Field.C a) => a -> T a -> T a +dilate k f = + Cons $ c f / k^2 + +shrink :: (Ring.C a) => a -> T a -> T a +shrink k f = + Cons $ c f * k^2 + +{- laws +fourier (convolve f g) = multiply (fourier f) (fourier g) + +dilate k (dilate m f) = dilate (k*m) f + +dilate k (shrink k f) = f + +variance (dilate k f) = k^2 * variance f + +variance (convolve f g) = variance f + variance g +-} diff --git a/src/MathObj/LaurentPolynomial.hs b/src/MathObj/LaurentPolynomial.hs --- a/src/MathObj/LaurentPolynomial.hs +++ b/src/MathObj/LaurentPolynomial.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2004-2006 @@ -33,7 +35,7 @@ import NumericPrelude hiding (div, negate, ) import qualified Data.List as List -import NumericPrelude.List (zipNeighborsWith) +import Data.List.HT (mapAdjacent) {- | Polynomial including negative exponents -} @@ -109,7 +111,7 @@ series ps = let es = map expon ps xs = map coeffs ps - ds = zipNeighborsWith subtract es + ds = mapAdjacent subtract es in Cons (head es) (addShiftedMany ds xs) {- | diff --git a/src/MathObj/Matrix.hs b/src/MathObj/Matrix.hs --- a/src/MathObj/Matrix.hs +++ b/src/MathObj/Matrix.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Mikael Johansson 2006 Maintainer : mik@math.uni-jena.de @@ -16,19 +18,20 @@ import qualified Algebra.Ring as Ring import qualified Algebra.Additive as Additive -import Algebra.Module((*>)) -import Algebra.Ring((*), fromInteger, scalarProduct) -import Algebra.Additive((+), (-), zero, subtract) +import Algebra.Module((*>), ) +import Algebra.Ring((*), fromInteger, scalarProduct, ) +import Algebra.Additive((+), (-), zero, subtract, ) -import Data.Array (Array, listArray, elems, bounds, (!), ixmap, range) +import Data.Array (Array, listArray, elems, bounds, (!), ixmap, range, ) import qualified Data.List as List -import Control.Monad (liftM2) -import Control.Exception (assert) +import Control.Monad (liftM2, ) +import Control.Exception (assert, ) -import NumericPrelude.List (outerProduct) -import NumericPrelude(Integer) -import PreludeBase hiding (zipWith) +import Data.Tuple.HT (swap, ) +import Data.List.HT (outerProduct, ) +import NumericPrelude (Integer, ) +import PreludeBase hiding (zipWith, ) {- | A matrix is a twodimensional array of ring elements, indexed by integers. @@ -43,14 +46,10 @@ -} --- candidate for Utility -twist :: (Integer,Integer) -> (Integer,Integer) -twist (x,y) = (y,x) - transpose :: T a -> T a transpose (Cons m) = let (lower,upper) = bounds m - in Cons (ixmap (twist lower, twist upper) twist m) + in Cons (ixmap (swap lower, swap upper) swap m) rows :: T a -> [[a]] rows (Cons m) = diff --git a/src/MathObj/PartialFraction.hs b/src/MathObj/PartialFraction.hs --- a/src/MathObj/PartialFraction.hs +++ b/src/MathObj/PartialFraction.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2007 Maintainer : numericprelude@henning-thielemann.de @@ -34,9 +34,10 @@ import Data.Map(Map) import qualified Data.Map as Map -import Data.Maybe(fromMaybe) -import NumericPrelude.List(replicateMatch, dropWhileRev) -import Data.List(group, sortBy, mapAccumR) +import Data.Maybe(fromMaybe, ) +import qualified Data.List.Match as Match +import Data.List.HT (dropWhileRev, ) +import Data.List (group, sortBy, mapAccumR, ) import PreludeBase hiding (zipWith) @@ -92,7 +93,7 @@ toFactoredFraction :: (PID.C a) => T a -> ([a], a) toFactoredFraction x@(Cons _ m) = let r = toFraction x - denoms = concat $ Map.elems $ indexMapMapWithKey (flip replicateMatch) m + denoms = concat $ Map.elems $ indexMapMapWithKey (flip Match.replicate) m numer = foldl (flip Ratio.scale) r denoms {- From the theory it must be Ratio.denominator denom==1. We could check this dynamically, if there would be an Eq instance. @@ -337,7 +338,7 @@ in (sc, (dis, rc))) . Map.elems . indexMapMapWithKey - (\d l -> (replicateMatch l d, multiToFraction d l)) + (\d l -> (Match.replicate l d, multiToFraction d l)) in removeZeros $ reduceHeads $ Cons z (mapApplySplit dsOrd (+) (uncurry (:) . carryRipple ds . map (ns*)) diff --git a/src/MathObj/Permutation.hs b/src/MathObj/Permutation.hs --- a/src/MathObj/Permutation.hs +++ b/src/MathObj/Permutation.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 Maintainer : numericprelude@henning-thielemann.de diff --git a/src/MathObj/Permutation/CycleList.hs b/src/MathObj/Permutation/CycleList.hs --- a/src/MathObj/Permutation/CycleList.hs +++ b/src/MathObj/Permutation/CycleList.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Mikael Johansson 2006 Maintainer : mik@math.uni-jena.de @@ -17,8 +17,8 @@ import Data.Array(Ix) import qualified Data.Array as Array -import NumericPrelude.List (takeMatch) -import NumericPrelude.Condition (toMaybe) +import qualified Data.List.Match as Match +import Data.Maybe.HT (toMaybe) import NumericPrelude (fromInteger) import PreludeBase @@ -71,18 +71,18 @@ orbit :: (Ord i) => (i -> i) -> i -> [i] orbit op x0 = takeUntilRepetition (iterate op x0) --- | candidates for NumericPrelude.List ? +-- | candidates for Utility ? takeUntilRepetition :: Ord a => [a] -> [a] takeUntilRepetition xs = let accs = scanl (flip Set.insert) Set.empty xs lenlist = takeWhile not (zipWith Set.member xs accs) - in takeMatch lenlist xs + in Match.take lenlist xs takeUntilRepetitionSlow :: Eq a => [a] -> [a] takeUntilRepetitionSlow xs = let accs = scanl (flip (:)) [] xs lenlist = takeWhile not (zipWith elem xs accs) - in takeMatch lenlist xs + in Match.take lenlist xs {- diff --git a/src/MathObj/Permutation/CycleList/Check.hs b/src/MathObj/Permutation/CycleList/Check.hs --- a/src/MathObj/Permutation/CycleList/Check.hs +++ b/src/MathObj/Permutation/CycleList/Check.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 Maintainer : numericprelude@henning-thielemann.de diff --git a/src/MathObj/Permutation/Table.hs b/src/MathObj/Permutation/Table.hs --- a/src/MathObj/Permutation/Table.hs +++ b/src/MathObj/Permutation/Table.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 Maintainer : numericprelude@henning-thielemann.de @@ -18,9 +18,10 @@ import Data.Array(Array,(!),(//),Ix) import qualified Data.Array as Array -import Data.List ((\\), nub, unfoldr) +import Data.List ((\\), nub, unfoldr, ) -import NumericPrelude.Condition (toMaybe) +import Data.Tuple.HT (swap, ) +import Data.Maybe.HT (toMaybe, ) -- import NumericPrelude (Integer) import PreludeBase hiding (cycle) @@ -70,10 +71,6 @@ then fmap (p!) q else error "compose: ranges differ" -- ++ show pRng ++ " /= " ++ show qRng) - --- | candidate for Utility -swap :: (a,b) -> (b,a) -swap (x,y) = (y,x) {- | diff --git a/src/MathObj/Polynomial.hs b/src/MathObj/Polynomial.hs --- a/src/MathObj/Polynomial.hs +++ b/src/MathObj/Polynomial.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Polynomials and rational functions in a single indeterminate. @@ -52,7 +54,7 @@ tensorProduct, tensorProductAlt, mulShear, mulShearTranspose, progression, differentiate, integrate, integrateInt, - fromRoots, alternate) + fromRoots, alternate, reverse, ) where import qualified Algebra.Differential as Differential @@ -71,10 +73,11 @@ import Algebra.Module((*>)) import Algebra.ZeroTestable(isZero) -import Control.Monad (liftM) +import Control.Monad (liftM, ) import qualified Data.List as List -import NumericPrelude.List - (zipWithOverlap, dropWhileRev, shear, shearTranspose, outerProduct) +import NumericPrelude.List (zipWithOverlap, ) +import Data.List.HT + (dropWhileRev, shear, shearTranspose, outerProduct, ) import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary)) @@ -82,8 +85,8 @@ import qualified PreludeBase as P import qualified NumericPrelude as NP -import PreludeBase hiding (const) -import NumericPrelude hiding (divMod, negate, stdUnit) +import PreludeBase hiding (const, reverse, ) +import NumericPrelude hiding (divMod, negate, stdUnit, ) newtype T a = Cons {coeffs :: [a]} @@ -304,7 +307,7 @@ divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a]) divMod x y = - let (y0:ys) = dropWhile isZero (reverse y) + let (y0:ys) = dropWhile isZero (List.reverse y) aux l xs' = if l < 0 then ([], xs') @@ -313,10 +316,10 @@ q0 = x0/y0 (d',m') = aux (l-1) (sub xs (scale q0 ys)) in (q0:d',m') - (d, m) = aux (length x - length y) (reverse x) + (d, m) = aux (length x - length y) (List.reverse x) in if isZero y then error "MathObj.Polynomial: division by zero" - else (reverse d, reverse m) + else (List.reverse d, List.reverse m) instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where divMod (Cons x) (Cons y) = @@ -382,13 +385,17 @@ alternate :: Additive.C a => [a] -> [a] alternate = zipWith ($) (cycle [id, Additive.negate]) +{-# INLINE reverse #-} +reverse :: Additive.C a => T a -> T a +reverse = lift1 alternate + {- see htam: Wavelet/DyadicResultant resultant :: Ring.C a => [a] -> [a] -> [a] resultant xs ys = -discriminant :: Ring.C a => [a] -> [a] +discriminant :: Ring.C a => [a] -> a discriminant xs = let degree = genericLength xs in parityFlip (safeDiv (degree*(degree-1)) 2) diff --git a/src/MathObj/PowerSeries.hs b/src/MathObj/PowerSeries.hs --- a/src/MathObj/PowerSeries.hs +++ b/src/MathObj/PowerSeries.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Power series, either finite or unbounded. (zipWith does exactly the @@ -24,7 +26,7 @@ import Algebra.Module((*>)) import Algebra.ZeroTestable(isZero) -import NumericPrelude.List(splitAtMatch) +import qualified Data.List.Match as Match import qualified NumericPrelude as NP import qualified PreludeBase as P @@ -256,7 +258,7 @@ divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a],[a]) divMod xs ys = let (yZero,yRem) = span isZero ys - (xMod, xRem) = splitAtMatch yZero xs + (xMod, xRem) = Match.splitAt yZero xs in (divide xRem yRem, xMod) instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where diff --git a/src/MathObj/PowerSeries/DifferentialEquation.hs b/src/MathObj/PowerSeries/DifferentialEquation.hs --- a/src/MathObj/PowerSeries/DifferentialEquation.hs +++ b/src/MathObj/PowerSeries/DifferentialEquation.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Lazy evaluation allows for the solution of differential equations in terms of power series. diff --git a/src/MathObj/PowerSeries/Example.hs b/src/MathObj/PowerSeries/Example.hs --- a/src/MathObj/PowerSeries/Example.hs +++ b/src/MathObj/PowerSeries/Example.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module MathObj.PowerSeries.Example where import qualified MathObj.PowerSeries as PS @@ -12,7 +12,7 @@ import Algebra.Additive (zero, subtract, negate) import Data.List (map, tail, cycle, zipWith, scanl, intersperse) -import NumericPrelude.List (sieve) +import Data.List.HT (sieve) import NumericPrelude (one, (*), (/), fromInteger, {-fromRational,-} pi) diff --git a/src/MathObj/PowerSeries/Mean.hs b/src/MathObj/PowerSeries/Mean.hs --- a/src/MathObj/PowerSeries/Mean.hs +++ b/src/MathObj/PowerSeries/Mean.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | This module computes power series for representing some means as generalized $f$-means. @@ -12,7 +12,7 @@ import qualified Algebra.Field as Field import qualified Algebra.Ring as Ring -import NumericPrelude.List (shearTranspose) +import Data.List.HT (shearTranspose) import NumericPrelude import PreludeBase diff --git a/src/MathObj/PowerSeries2.hs b/src/MathObj/PowerSeries2.hs --- a/src/MathObj/PowerSeries2.hs +++ b/src/MathObj/PowerSeries2.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Two-variate power series. @@ -20,8 +22,8 @@ import qualified NumericPrelude as NP import qualified PreludeBase as P -import Data.List (isPrefixOf) -import NumericPrelude.List (compareLength) +import Data.List (isPrefixOf, ) +import qualified Data.List.Match as Match import PreludeBase hiding (const) import NumericPrelude hiding (negate, stdUnit, @@ -49,7 +51,7 @@ check :: [[a]] -> [[a]] check xs = zipWith (\n x -> - if compareLength n x == EQ + if Match.compareLength n x == EQ then x else error "PowerSeries2.check: invalid length of sub-list") (iterate (():) [()]) xs diff --git a/src/MathObj/PowerSum.hs b/src/MathObj/PowerSum.hs --- a/src/MathObj/PowerSum.hs +++ b/src/MathObj/PowerSum.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2004-2005 @@ -30,7 +32,7 @@ import Control.Monad(liftM2) import qualified Data.List as List -import NumericPrelude.List (shearTranspose, sieve) +import Data.List.HT (shearTranspose, sieve) import PreludeBase as P hiding (const) import NumericPrelude as NP diff --git a/src/MathObj/RootSet.hs b/src/MathObj/RootSet.hs --- a/src/MathObj/RootSet.hs +++ b/src/MathObj/RootSet.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2004-2005 @@ -25,7 +25,7 @@ import qualified Algebra.Additive as Additive import qualified Algebra.ZeroTestable as ZeroTestable -import NumericPrelude.List (takeMatch) +import qualified Data.List.Match as Match import Control.Monad (liftM2) import PreludeBase as P hiding (const) @@ -81,12 +81,12 @@ liftPowerSum1Gen :: ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a]) liftPowerSum1Gen fromPS toPS op x = - takeMatch x (fromPS (op (toPS x))) + Match.take x (fromPS (op (toPS x))) liftPowerSum2Gen :: ([a] -> [a]) -> ([a] -> [a]) -> ([a] -> [a] -> [a]) -> ([a] -> [a] -> [a]) liftPowerSum2Gen fromPS toPS op x y = - takeMatch (undefined : liftM2 (,) (tail x) (tail y)) + Match.take (undefined : liftM2 (,) (tail x) (tail y)) (fromPS (op (toPS x) (toPS y))) diff --git a/src/MyPrelude.hs b/src/MyPrelude.hs --- a/src/MyPrelude.hs +++ b/src/MyPrelude.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module MyPrelude(module NumericPrelude, module PreludeBase, max, min, abs) where import NumericPrelude hiding (abs) import PreludeBase hiding (max,min) diff --git a/src/Number/Complex.hs b/src/Number/Complex.hs --- a/src/Number/Complex.hs +++ b/src/Number/Complex.hs @@ -1,5 +1,6 @@ -{-# OPTIONS_GHC -fno-implicit-prelude -fglasgow-exts #-} {-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- Rules should be processed -} {- | Module : Number.Complex @@ -24,6 +25,8 @@ (-:), scale, exp, + quarterLeft, + quarterRight, -- * Polar form fromPolar, @@ -75,7 +78,8 @@ import qualified Prelude as P import PreludeBase import NumericPrelude hiding (signum, exp, ) -import NumericPrelude.Text (showsInfixPrec, readsInfixPrec) +import Text.Show.HT (showsInfixPrec, ) +import Text.Read.HT (readsInfixPrec, ) -- import qualified Data.Typeable as Ty @@ -147,11 +151,11 @@ exp (Cons x y) = scale (Trans.exp x) (cis 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 +{-# INLINE quarterRight #-} +{-# INLINE quarterLeft #-} +quarterRight, quarterLeft :: (Additive.C a) => T a -> T a +quarterRight (Cons x y) = Cons y (-x) +quarterLeft (Cons x y) = Cons (-y) x {- | Scale a complex number to magnitude 1. @@ -362,7 +366,7 @@ {-# 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' + in if real z'<=0 then quarterRight z' else z' {-# INLINE stdUnit #-} stdUnit z@(Cons x y) = if z==zero @@ -371,7 +375,7 @@ let (x',sgn') = if y<0 || y==0 && x<0 then (negate x, -1) else (x, 1) - in if x'<=0 then orthoLeft sgn' else sgn' + in if x'<=0 then quarterLeft sgn' else sgn' instance (Ord a, ZeroTestable.C a, Units.C a) => PID.C (T a) where @@ -485,11 +489,11 @@ 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))) + asin z = quarterRight (log (quarterLeft z + sqrt (1 - z^2))) {-# INLINE acos #-} - acos z = orthoRight (log (z + orthoLeft (sqrt (1 - z^2)))) + acos z = quarterRight (log (z + quarterLeft (sqrt (1 - z^2)))) {-# INLINE atan #-} - atan z@(Cons x y) = orthoRight (log (Cons (1-y) x / sqrt (1+z^2))) + atan z@(Cons x y) = quarterRight (log (Cons (1-y) x / sqrt (1+z^2))) {- use the default implementation asinh z = log (z + sqrt (1+z^2)) diff --git a/src/Number/DimensionTerm.hs b/src/Number/DimensionTerm.hs --- a/src/Number/DimensionTerm.hs +++ b/src/Number/DimensionTerm.hs @@ -1,5 +1,5 @@ -{-# OPTIONS -fglasgow-exts #-} --- glasgow-exts for multi-parameter type class instances +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2008 License : GPL @@ -31,7 +31,7 @@ import System.Random (Random, randomR, random) -import NumericPrelude.Tuple (mapFst, ) +import Data.Tuple.HT (mapFst, ) import PreludeBase import Prelude () @@ -84,10 +84,13 @@ 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 +instance (Dim.IsScalar u, OccScalar.C a b) => OccScalar.C a (T u b) where + toScalar = + OccScalar.toScalar . toNumber . rewriteDimension Dim.toScalar + toMaybeScalar = + OccScalar.toMaybeScalar . toNumber . rewriteDimension Dim.toScalar + fromScalar = + rewriteDimension Dim.fromScalar . fromNumber . OccScalar.fromScalar instance (Dim.C u, Random a) => Random (T u a) where randomR (Cons l, Cons u) = mapFst Cons . randomR (l,u) diff --git a/src/Number/DimensionTerm/SI.hs b/src/Number/DimensionTerm/SI.hs --- a/src/Number/DimensionTerm/SI.hs +++ b/src/Number/DimensionTerm/SI.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2003 License : GPL diff --git a/src/Number/FixedPoint.hs b/src/Number/FixedPoint.hs --- a/src/Number/FixedPoint.hs +++ b/src/Number/FixedPoint.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 @@ -22,10 +22,12 @@ import qualified Algebra.Transcendental as Trans import qualified MathObj.PowerSeries.Example as PSE -import NumericPrelude.List (dropWhileRev, mapLast, padLeft) -import NumericPrelude.Condition (toMaybe) -import Data.List (transpose, unfoldr) -import Data.Char (intToDigit) +import NumericPrelude.List (mapLast, ) +import Data.Function.HT (powerAssociative, ) +import Data.List.HT (dropWhileRev, padLeft, ) +import Data.Maybe.HT (toMaybe, ) +import Data.List (transpose, unfoldr, ) +import Data.Char (intToDigit, ) import PreludeBase import NumericPrelude hiding (recip, sqrt, exp, sin, cos, tan, @@ -200,8 +202,8 @@ expFrac = expSmall den (frac-den2) in case compare int 0 of EQ -> expFrac - GT -> reduceRepeated (mul den) expFrac (eConst den) int - LT -> reduceRepeated (mul den) expFrac (recipEConst den) (-int) + GT -> powerAssociative (mul den) expFrac (eConst den) int + LT -> powerAssociative (mul den) expFrac (recipEConst den) (-int) -- LT -> nest (-int) (divide den e) expFrac diff --git a/src/Number/FixedPoint/Check.hs b/src/Number/FixedPoint/Check.hs --- a/src/Number/FixedPoint/Check.hs +++ b/src/Number/FixedPoint/Check.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Number.FixedPoint.Check where import qualified Number.FixedPoint as FP diff --git a/src/Number/NonNegative.hs b/src/Number/NonNegative.hs --- a/src/Number/NonNegative.hs +++ b/src/Number/NonNegative.hs @@ -1,4 +1,14 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# OPTIONS -XNoImplicitPrelude -fno-warn-orphans #-} + +{- +Rationale for -fno-warn-orphans: + * The orphan instances can't be put into Numeric.NonNegative.Wrapper + since that's in another package. + * We had to spread the instance declarations + over the modules defining the typeclasses instantiated. + Do we want that? +-} + {- | Copyright : (c) Henning Thielemann 2007 @@ -41,7 +51,7 @@ import qualified Prelude as P import PreludeBase -import NumericPrelude.Tuple (mapSnd, mapPair, ) +import Data.Tuple.HT (mapSnd, mapPair, ) import NumericPrelude hiding (Int, Integer, Float, Double, Rational, ) diff --git a/src/Number/NonNegativeChunky.hs b/src/Number/NonNegativeChunky.hs --- a/src/Number/NonNegativeChunky.hs +++ b/src/Number/NonNegativeChunky.hs @@ -39,7 +39,7 @@ import Test.QuickCheck (Arbitrary(..)) import NumericPrelude -import NumericPrelude.Tuple (mapFst, mapPair, ) +import Data.Tuple.HT (mapFst, mapPair, ) import PreludeBase import qualified Prelude as P98 diff --git a/src/Number/OccasionallyScalarExpression.hs b/src/Number/OccasionallyScalarExpression.hs --- a/src/Number/OccasionallyScalarExpression.hs +++ b/src/Number/OccasionallyScalarExpression.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2004 License : GPL @@ -23,7 +25,7 @@ import qualified Algebra.ZeroTestable as ZeroTestable import Algebra.Algebraic (sqrt, (^/)) -import Algebra.OccasionallyScalar as OccScalar +import qualified Algebra.OccasionallyScalar as OccScalar import Data.Maybe(fromMaybe) import Data.Array(listArray,(!)) @@ -81,13 +83,17 @@ lift :: (v -> v) -> (T a v -> T a v) lift f (Cons xe x) = Cons xe (f x) +fromScalar :: (Show v, OccScalar.C a v) => + a -> T a v +fromScalar = OccScalar.fromScalar + scalarMap :: (Show v, OccScalar.C a v) => (a -> a) -> (T a v -> T a v) -scalarMap f x = fromScalar (f (toScalar x)) +scalarMap f x = OccScalar.fromScalar (f (OccScalar.toScalar x)) scalarMap2 :: (Show v, OccScalar.C a v) => (a -> a -> a) -> (T a v -> T a v -> T a v) -scalarMap2 f x y = fromScalar (f (toScalar x) (toScalar y)) +scalarMap2 f x y = OccScalar.fromScalar (f (OccScalar.toScalar x) (OccScalar.toScalar y)) instance (Show v) => Show (T a v) where @@ -133,7 +139,7 @@ instance (Trans.C a, Field.C v, Show v, OccScalar.C a v) => Trans.C (T a v) where - pi = fromScalar (pi::a) + pi = fromScalar pi log = scalarMap log exp = scalarMap exp logBase = scalarMap2 logBase @@ -158,9 +164,9 @@ fromMaybe (error (show xe ++ " is not a scalar value.\n" ++ showUnitError True 0 x xe)) - (toMaybeScalar x) - toMaybeScalar (Cons _ x) = toMaybeScalar x - fromScalar = fromValue . fromScalar + (OccScalar.toMaybeScalar x) + toMaybeScalar (Cons _ x) = OccScalar.toMaybeScalar x + fromScalar = fromValue . OccScalar.fromScalar {- @@ -178,7 +184,7 @@ instance Functor (T i) where fmap f (Cons xu x) = if Unit.isScalar xu - then fromScalar (f x) + then OccScalar.fromScalar (f x) else error "Physics.Quantity.Value.fmap: function for scalars, only" instance Monad (T i) where @@ -186,5 +192,5 @@ if Unit.isScalar xu then f x else error "Physics.Quantity.Value.(>>=): function for scalars, only" - return = fromScalar + return = OccScalar.fromScalar -} diff --git a/src/Number/PartiallyTranscendental.hs b/src/Number/PartiallyTranscendental.hs --- a/src/Number/PartiallyTranscendental.hs +++ b/src/Number/PartiallyTranscendental.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Define Transcendental functions on arbitrary fields. These functions are defined for only a few (in most cases only one) arguments, diff --git a/src/Number/Peano.hs b/src/Number/Peano.hs --- a/src/Number/Peano.hs +++ b/src/Number/Peano.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2007 Maintainer : numericprelude@henning-thielemann.de @@ -25,11 +25,13 @@ import qualified Algebra.ToRational as ToRational import qualified Algebra.NonNegative as NonNeg +import Data.Maybe (catMaybes, ) import Data.Array(Ix(..)) import qualified Prelude as P98 import qualified PreludeBase as P import qualified NumericPrelude as NP +import Data.List.HT (mapAdjacent, shearTranspose, ) import PreludeBase import NumericPrelude @@ -66,17 +68,10 @@ subNeg x Zero = (True, x) subNeg (Succ x) (Succ y) = subNeg x y -{- efficient implementation of x0 <= x1 && x1 <= x2 ... -} -isAscending :: [T] -> Bool -isAscending [] = True -isAscending (x:xs) = - if isZero x - then isAscending xs - else not (any isZero xs) && - isAscending (map pred xs) mul :: T -> T -> T mul Zero _ = Zero +mul _ Zero = Zero mul (Succ x) y = add y (mul x y) fromPosEnum :: (ZeroTestable.C a, Enum a) => a -> T @@ -141,6 +136,21 @@ max Zero y = y max x Zero = x + {- + This special implementation works also for undefined < Zero. + Thanks to Peter Divianszky for the hint. + -} + _ < Zero = False + Zero < _ = True + Succ n < Succ m = n < m + + x > y = y < x + + x <= y = not (y < x) + + x >= y = not (x < y) + + {- | cf. To how to find the shortest list in a list of lists efficiently, this means, also in the presence of infinite lists. @@ -183,6 +193,125 @@ argMaximum :: [(T,a)] -> a argMaximum = snd . foldl1 argMaxFull + + +-- isAscending - naive implementations + +{- | +@x0 <= x1 && x1 <= x2 ... @ +for possibly infinite numbers in finite lists. +-} +isAscendingFiniteList :: [T] -> Bool +isAscendingFiniteList [] = True +isAscendingFiniteList (x:xs) = + let decrement (Succ y) = Just y + decrement _ = Nothing + in case x of + Zero -> isAscendingFiniteList xs + Succ xd -> + case mapM decrement xs of + Nothing -> False + Just xsd -> isAscendingFiniteList (xd : xsd) + +isAscendingFiniteNumbers :: [T] -> Bool +isAscendingFiniteNumbers = and . mapAdjacent (<=) + + +-- isAscending - sophisticated implementations - explicit + +toListMaybe :: a -> T -> [Maybe a] +toListMaybe a = + let recurse Zero = [Just a] + recurse (Succ x) = Nothing : recurse x + in recurse + +{- | +In @glue x y == (z,r,b)@ +@z@ represents @min x y@, +@r@ represents @max x y - min x y@, +and @x<=y == b@. + +Cf. Numeric.NonNegative.Chunky +-} +glue :: T -> T -> (T, T, Bool) +glue Zero ys = (Zero, ys, True) +glue xs Zero = (Zero, xs, False) +glue (Succ xs) (Succ ys) = + let (common, difference, sign) = glue xs ys + in (Succ common, difference, sign) + +{- +Implementation notes: +We check all pairs of adjacent numbers for correct order. +We obtain a set of booleans, which must all be True. +The order of checking these booleans is crucial. +Pairs of numbers that are infinitely big or infinitely far in the list +must be checked \"last\". +Thus we order the booleans according to their computation costs +(list position + magnitude of number) +using 'shearTranspose'. +-} +isAscending :: [T] -> Bool +isAscending = + and . catMaybes . concat . + shearTranspose . + mapAdjacent (\x y -> + let (costs0,_,le) = glue x y + in toListMaybe le costs0) + + +-- isAscending - use a cost measuring data type (could generalized to a monad, when considered as Writer monad, see htam and unique-logic packages + +-- following an idea of vixy http://moonpatio.com:8080/fastcgi/hpaste.fcgi/view?id=562 + +data Valuable a = Valuable {costs :: T, value :: a} + deriving (Show, Eq, Ord) + + +increaseCosts :: T -> Valuable a -> Valuable a +increaseCosts inc ~(Valuable c x) = Valuable (inc+c) x + +{- | +Compute '(&&)' with minimal costs. +-} +infixr 3 &&~ +(&&~) :: Valuable Bool -> Valuable Bool -> Valuable Bool +(&&~) (Valuable xc xb) (Valuable yc yb) = + let (minc,difc,le) = glue xc yc + (bCheap,bExpensive) = + if le + then (xb,yb) + else (yb,xb) + in increaseCosts minc $ + uncurry Valuable $ + if bCheap + then (difc, bExpensive) + else (Zero, False) + +andW :: [Valuable Bool] -> Valuable Bool +andW = + foldr + (\b acc -> b &&~ increaseCosts one acc) + (Valuable Zero True) + +leW :: T -> T -> Valuable Bool +leW x y = + let (minc,_difc,le) = glue x y + in Valuable minc le + +isAscendingW :: [T] -> Valuable Bool +isAscendingW = + andW . mapAdjacent leW + +{- +test with + +*Number.Peano> isAscendingW [0,infinity,infinity,5] +False +-} + + +-- instances instance Real.C T where signum Zero = zero diff --git a/src/Number/Physical.hs b/src/Number/Physical.hs --- a/src/Number/Physical.hs +++ b/src/Number/Physical.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2003-2006 License : GPL @@ -34,7 +36,7 @@ import Control.Monad(guard,liftM,liftM2) -import NumericPrelude.Condition(toMaybe) +import Data.Maybe.HT(toMaybe) import Data.Maybe(fromMaybe) import NumericPrelude diff --git a/src/Number/Physical/Read.hs b/src/Number/Physical/Read.hs --- a/src/Number/Physical/Read.hs +++ b/src/Number/Physical/Read.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2004 License : GPL diff --git a/src/Number/Physical/Show.hs b/src/Number/Physical/Show.hs --- a/src/Number/Physical/Show.hs +++ b/src/Number/Physical/Show.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2004 License : GPL diff --git a/src/Number/Physical/Unit.hs b/src/Number/Physical/Unit.hs --- a/src/Number/Physical/Unit.hs +++ b/src/Number/Physical/Unit.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2003-2006 License : GPL @@ -19,7 +19,7 @@ import qualified Number.Ratio as Ratio -import NumericPrelude.Condition(toMaybe) +import Data.Maybe.HT(toMaybe) import PreludeBase import NumericPrelude diff --git a/src/Number/Physical/UnitDatabase.hs b/src/Number/Physical/UnitDatabase.hs --- a/src/Number/Physical/UnitDatabase.hs +++ b/src/Number/Physical/UnitDatabase.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2003 License : GPL @@ -19,7 +19,7 @@ -- import Algebra.Module((*>)) import Algebra.NormedSpace.Sum(norm) -import NumericPrelude.Condition (toMaybe) +import Data.Maybe.HT (toMaybe) import Data.List (findIndices, partition, unfoldr, find, minimumBy) import PreludeBase diff --git a/src/Number/Positional.hs b/src/Number/Positional.hs --- a/src/Number/Positional.hs +++ b/src/Number/Positional.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 License : GPL @@ -31,9 +31,15 @@ import qualified Data.List as List import Data.Char (intToDigit) -import NumericPrelude.Condition (toMaybe, select, if') -import NumericPrelude.List (replicateMatch, sliceVert, zipNeighborsWith, - padLeft, padRight, mapLast) +import qualified Data.List.Match as Match +import Data.Function.HT (powerAssociative, nest, ) +import Data.Tuple.HT (swap, ) +import Data.Maybe.HT (toMaybe, ) +import Data.Bool.HT (select, if', ) +import NumericPrelude.List (mapLast, ) +import Data.List.HT + (sliceVertical, mapAdjacent, + padLeft, padRight, ) {- @@ -106,8 +112,8 @@ [] -> (x:xs) -- keep trailing zeros, because they show precision in 'show' functions (y:ys) -> if y>=0 -- equivalent to y>0 - then x : replicateMatch xZeros 0 ++ recurse xRem - else (x-1) : replicateMatch xZeros (b-1) ++ recurse ((y+b) : ys) + then x : Match.replicate xZeros 0 ++ recurse xRem + else (x-1) : Match.replicate xZeros (b-1) ++ recurse ((y+b) : ys) in recurse @@ -381,7 +387,7 @@ powerBasis b e (xe,xm) = let (ye,r) = divMod xe e (y0,y1) = splitAtPadZero (r+1) xm - y1pad = mapLast (padRight 0 e) (sliceVert e y1) + y1pad = mapLast (padRight 0 e) (sliceVertical e y1) in (ye, map (mantissaToNum b) (y0 : y1pad)) {- | @@ -564,7 +570,7 @@ addMany _ [] = zero addMany b ys = let recurse xs = - case map (addSome b) (sliceVert b xs) of + case map (addSome b) (sliceVertical b xs) of [s] -> s sums -> recurse sums in recurse ys @@ -594,9 +600,9 @@ seriesPlain _ [] = error "empty series: don't know a good exponent" seriesPlain b summands = let (es,m:ms) = unzip (map (uncurry (align b)) summands) - eDifs = zipNeighborsWith (-) es - eDifLists = sliceVert (pred b) eDifs - mLists = sliceVert (pred b) ms + eDifs = mapAdjacent (-) es + eDifLists = sliceVertical (pred b) eDifs + mLists = sliceVertical (pred b) ms accum sumM (eDifList,mList) = let subM = LPoly.addShiftedMany eDifList (sumM:mList) -- lazy unary sum @@ -640,7 +646,7 @@ -- must get a case for negative index splitAtMatchPadZero :: [()] -> Mantissa -> (Mantissa, Mantissa) -splitAtMatchPadZero n [] = (replicateMatch n 0, []) +splitAtMatchPadZero n [] = (Match.replicate n 0, []) splitAtMatchPadZero [] xs = ([], xs) splitAtMatchPadZero n (x:xs) = let (ys, zs) = splitAtMatchPadZero (tail n) xs @@ -876,7 +882,7 @@ let b = fromIntegral bInt chunkLengths = iterate (2*) 1 xChunks = map (mantissaToNum bInt) $ snd $ - List.mapAccumL (\x cl -> flipPair (splitAtPadZero cl x)) + List.mapAccumL (\x cl -> swap (splitAtPadZero cl x)) xs chunkLengths basisPowers = iterate (^2) b truncXs = scanl (\acc (bp,frac) -> acc*bp+frac) x0 @@ -897,18 +903,24 @@ @inits = foldr (\x ys -> [] : map (x:) ys) [[]]@ This is too strict for our application. -@ -Prelude> List.inits (0:1:2:undefined) -[[],[0],[0,1]*** Exception: Prelude.undefined -@ -The following routine is more lazy -but restricted to infinite lists. +> Prelude> List.inits (0:1:2:undefined) +> [[],[0],[0,1]*** Exception: Prelude.undefined + +The following routine is more lazy than 'List.inits' +and even lazier than 'Data.List.HT.inits' from @utility-ht@ package, +but it is restricted to infinite lists. +This degree of laziness is needed for @sqrtFP@. + +> Prelude> lazyInits (0:1:2:undefined) +> [[],[0],[0,1],[0,1,2],[0,1,2,*** Exception: Prelude.undefined -} lazyInits :: [a] -> [[a]] lazyInits ~(x:xs) = [] : map (x:) (lazyInits xs) --- this pattern is never reached, GHC does not complain about it -lazyInits [] = [[]] +{- +The lazy match above is irrefutable, +so the pattern @[]@ would never be reached. +-} @@ -969,7 +981,7 @@ cardPower :: Basis -> Integer -> T -> T -> T cardPower b expon neutral x = if expon >= 0 - then reduceRepeated (mul b) neutral x expon + then powerAssociative (mul b) neutral x expon else error "negative exponent - use intPower" @@ -1182,8 +1194,8 @@ ei = fromIntegral e y = trim $ if e<0 - then reduceRepeated (mul b) x (eConst b) (-ei) - else reduceRepeated (mul b) x (recipEConst b) ei + then powerAssociative (mul b) x (eConst b) (-ei) + else powerAssociative (mul b) x (recipEConst b) ei estimate = liftDoubleApprox b log y residue = mul b (expSmall b (neg b estimate)) y in addSome b [(0,[e]), estimate, lnSmall b residue] @@ -1330,18 +1342,6 @@ {- * auxilary functions -} - -{- | -Candidate for a Utility module. --} -nest :: Int -> (a -> a) -> a -> a -nest 0 _ x = x -nest n f x = f (nest (n-1) f x) - - -flipPair :: (a,b) -> (b,a) -flipPair ~(x,y) = (y,x) - sliceVertPair :: [a] -> [(a,a)] sliceVertPair (x0:x1:xs) = (x0,x1) : sliceVertPair xs diff --git a/src/Number/Positional/Check.hs b/src/Number/Positional/Check.hs --- a/src/Number/Positional/Check.hs +++ b/src/Number/Positional/Check.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2006 License : GPL diff --git a/src/Number/Quaternion.hs b/src/Number/Quaternion.hs --- a/src/Number/Quaternion.hs +++ b/src/Number/Quaternion.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Maintainer : numericprelude@henning-thielemann.de Stability : provisional @@ -59,7 +61,8 @@ import qualified Prelude as P import PreludeBase import NumericPrelude hiding (signum) -import NumericPrelude.Text (showsInfixPrec, readsInfixPrec) +import Text.Show.HT (showsInfixPrec, ) +import Text.Read.HT (readsInfixPrec, ) {- TODO: diff --git a/src/Number/Ratio.hs b/src/Number/Ratio.hs --- a/src/Number/Ratio.hs +++ b/src/Number/Ratio.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Module : Number.Ratio Copyright : (c) Henning Thielemann, Dylan Thurston 2006 @@ -38,8 +38,7 @@ import Algebra.Additive (zero, (+), (-), negate) import Algebra.ZeroTestable (isZero) --- import NumericPrelude.Monad(untilM) -import Control.Monad(liftM, liftM2) +import Control.Monad(liftM, liftM2, ) import Test.QuickCheck (Arbitrary(arbitrary,coarbitrary)) import System.Random (Random(..), RandomGen, ) diff --git a/src/Number/ResidueClass.hs b/src/Number/ResidueClass.hs --- a/src/Number/ResidueClass.hs +++ b/src/Number/ResidueClass.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Number.ResidueClass where import qualified Algebra.PrincipalIdealDomain as PID @@ -10,7 +10,7 @@ import PreludeBase import NumericPrelude hiding (recip) -import NumericPrelude.Condition (toMaybe) +import Data.Maybe.HT (toMaybe) import Data.Maybe (fromMaybe) diff --git a/src/Number/ResidueClass/Check.hs b/src/Number/ResidueClass/Check.hs --- a/src/Number/ResidueClass/Check.hs +++ b/src/Number/ResidueClass/Check.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Number.ResidueClass.Check where import qualified Number.ResidueClass as Res @@ -14,10 +14,11 @@ import PreludeBase import NumericPrelude (Int, Integer, mod, ) -import qualified NumericPrelude -import NumericPrelude.Condition (toMaybe) -import NumericPrelude.Text (showsInfixPrec, readsInfixPrec) +import Data.Maybe.HT (toMaybe, ) +import Text.Show.HT (showsInfixPrec, ) +import Text.Read.HT (readsInfixPrec, ) + infix 7 /:, `Cons` {- | @@ -81,10 +82,10 @@ zero m = Cons m Additive.zero one :: (Ring.C a) => a -> T a -one m = Cons m NumericPrelude.one +one m = Cons m Ring.one fromInteger :: (Integral.C a) => a -> Integer -> T a -fromInteger m x = fromRepresentative m (NumericPrelude.fromInteger x) +fromInteger m x = fromRepresentative m (Ring.fromInteger x) @@ -110,5 +111,5 @@ instance (Eq a, PID.C a) => Field.C (T a) where (/) = lift2 Res.divide - recip = lift1 (flip Res.divide NumericPrelude.one) + recip = lift1 (flip Res.divide Ring.one) fromRational' = error "no conversion from rational to residue class" diff --git a/src/Number/ResidueClass/Func.hs b/src/Number/ResidueClass/Func.hs --- a/src/Number/ResidueClass/Func.hs +++ b/src/Number/ResidueClass/Func.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Number.ResidueClass.Func where import qualified Number.ResidueClass as Res diff --git a/src/Number/ResidueClass/Maybe.hs b/src/Number/ResidueClass/Maybe.hs --- a/src/Number/ResidueClass/Maybe.hs +++ b/src/Number/ResidueClass/Maybe.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Number.ResidueClass.Maybe where import qualified Number.ResidueClass as Res diff --git a/src/Number/ResidueClass/Reader.hs b/src/Number/ResidueClass/Reader.hs --- a/src/Number/ResidueClass/Reader.hs +++ b/src/Number/ResidueClass/Reader.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Number.ResidueClass.Reader where import qualified Number.ResidueClass as Res diff --git a/src/Number/SI.hs b/src/Number/SI.hs --- a/src/Number/SI.hs +++ b/src/Number/SI.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2003-2006 License : GPL @@ -36,8 +38,10 @@ import qualified Algebra.Additive as Additive import qualified Algebra.ZeroTestable as ZeroTestable -import Algebra.Algebraic (sqrt, (^/)) +import Algebra.Algebraic (sqrt, (^/), ) +import Data.Tuple.HT (mapFst, ) + import qualified Prelude as P import NumericPrelude @@ -45,7 +49,7 @@ newtype T a v = Cons (PValue v) -{- -fglasgow-exts allow even this +{- LANGUAGE GeneralizedNewtypeDeriving allows even this deriving (Monad, Functor) -} @@ -102,18 +106,25 @@ instance Eq v => Eq (T a v) where (==) = lift2Gen (==) +showNat :: (Show v, Field.C a, Ord a, NormedMax.C a v) => + UnitDatabase.T Dimension a -> T a v -> String +showNat db = + liftGen (PVShow.showNat db) + instance (Show v, Ord a, Trans.C a, NormedMax.C a v) => Show (T a v) where showsPrec prec x = showParen (prec > PVShow.mulPrec) - (liftGen (PVShow.showNat - (SIUnit.databaseShow :: UnitDatabase.T Dimension a)) x ++) + (showNat SIUnit.databaseShow x ++) +readsNat :: (Read v, VectorSpace.C a v) => + UnitDatabase.T Dimension a -> Int -> ReadS (T a v) +readsNat db prec = + map (mapFst Cons) . PVRead.readsNat db prec + instance (Read v, Ord a, Trans.C a, VectorSpace.C a v) => Read (T a v) where - readsPrec prec str = - map (\(x,s) -> (Cons x, s)) - (PVRead.readsNat (SIUnit.databaseRead :: UnitDatabase.T Dimension a) prec str) + readsPrec = readsNat SIUnit.databaseRead instance (Additive.C v) => Additive.C (T a v) where zero = Cons zero diff --git a/src/Number/SI/Unit.hs b/src/Number/SI/Unit.hs --- a/src/Number/SI/Unit.hs +++ b/src/Number/SI/Unit.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} {- | Copyright : (c) Henning Thielemann 2003 License : GPL diff --git a/src/NumericPrelude.hs b/src/NumericPrelude.hs --- a/src/NumericPrelude.hs +++ b/src/NumericPrelude.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module NumericPrelude ( {- Additive -} (+), (-), negate, zero, subtract, sum, sum1, {- ZeroTestable -} isZero, @@ -15,7 +15,6 @@ {- RealTrans -} atan2, {- ToRational -} toRational, {- ToInteger -} toInteger, fromIntegral, - reduceRepeated, {- Units -} isUnit, stdAssociate, stdUnit, stdUnitInv, {- PID -} extendedGCD, gcd, lcm, euclid, extendedEuclid, {- Ratio -} Rational, (%), numerator, denominator, @@ -43,4 +42,3 @@ import Algebra.ToRational (toRational, ) import Prelude (Int, Integer, Float, Double) -import NumericPrelude.List (reduceRepeated) diff --git a/src/NumericPrelude/Condition.hs b/src/NumericPrelude/Condition.hs deleted file mode 100644 --- a/src/NumericPrelude/Condition.hs +++ /dev/null @@ -1,50 +0,0 @@ -module NumericPrelude.Condition where - -{- some routines that are copied from Henning's Useful.hs -} - -{- | -Returns 'Just' if the precondition is fulfilled. --} -{-# INLINE toMaybe #-} -toMaybe :: Bool -> a -> Maybe a -toMaybe False _ = Nothing -toMaybe True x = Just x - -{- | -A purely functional implementation of @if@. -Very useful in connection with 'zipWith3'. --} -{-# INLINE if' #-} -if' :: - Bool {-^ condition -} - -> a {-^ then -} - -> a {-^ else -} - -> a -if' True x _ = x -if' False _ y = y - -{- | -From a list of expressions choose the one, -whose condition is true. - -> select "zero" -> [(x>0, "positive"), -> (x<0, "negative")] --} -{-# INLINE select #-} -select :: a -> [(Bool, a)] -> a -select = foldr (uncurry if') - - --- precedence below (||) and (&&) -infix 1 `implies` - -{- | -Logical operator for implication. - -Funnily because of the ordering of 'Bool' it holds @implies == (<=)@. --} -{-# INLINE implies #-} -implies :: Bool -> Bool -> Bool -implies prerequisite conclusion = - not prerequisite || conclusion diff --git a/src/NumericPrelude/List.hs b/src/NumericPrelude/List.hs --- a/src/NumericPrelude/List.hs +++ b/src/NumericPrelude/List.hs @@ -1,85 +1,6 @@ module NumericPrelude.List where -import Data.List (unfoldr, genericReplicate) -import NumericPrelude.Condition (toMaybe) - -{- * Slice lists -} - - -{-| keep every k-th value from the list - - Since these implementations check for the end of lists, - they may fail in fixpoint computations on infinite lists. -} -sieve, sieve', sieve'', sieve''' :: Int -> [a] -> [a] -sieve k = - unfoldr (\xs -> toMaybe (not (null xs)) (head xs, drop k xs)) - -sieve' k = map head . sliceVert k - --- this one works only on finite lists -sieve'' k x = map (x!!) [0,k..(length x-1)] - -sieve''' k = map head . takeWhile (not . null) . iterate (drop k) - -{- sliceHoriz is faster than sliceHoriz' but consumes slightly more memory - (although it needs no swapping) -} -sliceHoriz, sliceHoriz' :: Int -> [a] -> [[a]] -sliceHoriz n = - map (sieve n) . take n . iterate (drop 1) - -sliceHoriz' n = - foldr (\x ys -> let y = last ys in takeMatch ys ((x:y):ys)) (replicate n []) - - -sliceVert, sliceVert' :: Int -> [a] -> [[a]] -sliceVert n = - map (take n) . takeWhile (not . null) . iterate (drop n) - {- takeWhile must be performed before (map take) - in order to handle (n==0) correctly -} - -sliceVert' n = - unfoldr (\x -> toMaybe (not (null x)) (splitAt n x)) - - -{- * 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 - in (x:as,bs) -splitAtMatch _ [] = ([],[]) -splitAtMatch [] xs = ([],xs) - -{-# INLINE replicateMatch #-} -replicateMatch :: [a] -> b -> [b] -replicateMatch xs y = - takeMatch xs (repeat y) - -{- | -Compare the length of two lists over different types. -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 -compareLength (_:_) [] = GT -compareLength [] (_:_) = LT - - - +import Data.List.HT (switchL, switchR, ) {- * Zip lists -} @@ -121,192 +42,31 @@ 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) - - -{- * Lists of lists -} - {- | -Transform - -@ - [[00,01,02,...], [[00], - [10,11,12,...], --> [10,01], - [20,21,22,...], [20,11,02], - ...] ...] -@ - -With @concat . shear@ you can perform a Cantor diagonalization, -that is an enumeration of all elements of the sub-lists -where each element is reachable within a finite number of steps. -It is also useful for polynomial multiplication (convolution). --} -shear :: [[a]] -> [[a]] -shear xs@(_:_) = - let (y:ys,zs) = unzip (map (splitAt 1) xs) - zipConc (a:as) (b:bs) = (a++b) : zipConc as bs - zipConc [] bs = bs - zipConc as [] = as - in y : zipConc ys (shear (dropWhileRev null zs)) - {- Dropping trailing empty lists is necessary, - otherwise finite lists are filled with empty lists. -} -shear [] = [] - -{- | -Transform - -@ - [[00,01,02,...], [[00], - [10,11,12,...], --> [01,10], - [20,21,22,...], [02,11,20], - ...] ...] -@ - -It's like 'shear' but the order of elements in the sub list is reversed. -Its implementation seems to be more efficient than that of 'shear'. -If the order does not matter, better choose 'shearTranspose'. --} -shearTranspose :: [[a]] -> [[a]] -shearTranspose = - let -- zipCons is like zipWith (:) keep lists which are too long - zipCons (x:xs) (y:ys) = (x:y) : zipCons xs ys - zipCons [] ys = ys - zipCons xs [] = map (:[]) xs - aux (x:xs) yss = [x] : zipCons xs yss - aux [] yss = []:yss - in foldr aux [] - - -{- | -Operate on each combination of elements of the first and the second list. -In contrast to the list instance of 'Monad.liftM2' -in holds the results in a list of lists. -It holds -@concat (outerProduct f xs ys) == liftM2 f xs ys@ --} -outerProduct :: (a -> b -> c) -> [a] -> [b] -> [[c]] -outerProduct f xs ys = map (flip map ys . f) xs - - - -{- * Various -} - -{-# INLINE partitionMaybe #-} -partitionMaybe :: (a -> Maybe b) -> [a] -> ([b], [a]) -partitionMaybe f = - foldr (\x ~(y,z) -> case f x of - Just x' -> (x' : y, z) - Nothing -> (y, x : z)) ([],[]) - -{- | -It holds @splitLast xs == (init xs, last xs)@, -but 'splitLast' is more efficient -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) -splitLast (x:xs) = - let (xs', lastx) = splitLast xs in (x:xs', lastx) - -propSplitLast :: Eq a => [a] -> Bool -propSplitLast xs = - splitLast xs == (init xs, last xs) - -{- | -Remove the longest suffix of elements satisfying p. -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) [] - - -{- | 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 - -- error "mapLast: empty list" - recurse (x:[]) = f x : [] - 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 - - -{- | -@reduceRepeated@ is an auxiliary function that, -for an associative operation @op@, -computes the same value as - - @reduceRepeated op a0 a n = foldr op a0 (genericReplicate n a)@ - -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 -reduceRepeated op a0 a n = - if even n - then reduceRepeated op a0 (op a a) (div n 2) - else reduceRepeated op (op a0 a) (op a a) (div n 2) - -reduceRepeatedSlow op a0 a n = - foldr op a0 (genericReplicate n a) - - -{- | -For an associative operation @op@ this computes - @iterateAssoc op a = iterate (op a) a@ -but it is even faster than @map (reduceRepeated op a a) [0..]@ -since it shares temporary results. + switchL [] + (\x xs -> + uncurry (:) $ + foldr (\x1 k x0 -> (x0, uncurry (:) (k x1))) + (\x0 -> (f x0, [])) xs x) -The idea is: -From the list @map (reduceRepeated op a a) [0,(2*n)..]@ -we compute the list @map (reduceRepeated op a a) [0,n..]@, -and iterate that until @n==1@. --} -iterateAssoc, iterateLeaky :: (a -> a -> a) -> a -> [a] -iterateAssoc op a = - foldr (\pow xs -> pow : concatMap (\x -> [x, op x pow]) xs) - undefined (iterate (\x -> op x x) a) +mapLast' :: (a -> a) -> [a] -> [a] +mapLast' f = + let recourse [] = [] -- behaviour as needed in powerBasis + -- otherwise: error "mapLast: empty list" + recourse (x:xs) = + uncurry (:) $ + if null xs + then (f x, []) + else (x, recourse xs) + in recourse -{- | -This is equal to 'iterateAssoc'. -The idea is the following: -The list we search is the fixpoint of the function: -"Square all elements of the list, -then spread it and fill the holes with successive numbers -of their left neighbour." -This also preserves log n applications per value. -However it has a space leak, -because for the value with index @n@ -all elements starting at @div n 2@ must be kept. --} -iterateLeaky op x = - let merge (a:as) b = a : merge b as - merge _ _ = error "iterateLeaky: an empty list cannot occur" - sqrs = map (\y -> op y y) z - z = x : merge sqrs (map (op x) sqrs) - in z +mapLast'' :: (a -> a) -> [a] -> [a] +mapLast'' f = + switchR [] (\xs x -> xs ++ [f x]) diff --git a/src/NumericPrelude/Monad.hs b/src/NumericPrelude/Monad.hs deleted file mode 100644 --- a/src/NumericPrelude/Monad.hs +++ /dev/null @@ -1,10 +0,0 @@ -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 - if p x - then return x - else untilM p m diff --git a/src/NumericPrelude/Text.hs b/src/NumericPrelude/Text.hs deleted file mode 100644 --- a/src/NumericPrelude/Text.hs +++ /dev/null @@ -1,32 +0,0 @@ -module NumericPrelude.Text where - -{-* 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 = - showParen - (prec >= opPrec) - (showsPrec opPrec x . showString " " . - showString opStr . showString " " . - 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 = - readParen - (prec >= opPrec) - ((\s -> [(const . cons, s)]) .> - readsPrec opPrec .> - (filter ((opStr==).fst) . lex) .> - readsPrec opPrec) - -{-| Compose two parsers sequentially. -} -infixl 9 .> -(.>) :: ReadS (b->c) -> ReadS b -> ReadS c -(.>) ra rb = - concatMap (\(f,rest) -> map (\(b, rest') -> (f b, rest')) (rb rest)) . ra diff --git a/src/NumericPrelude/Tuple.hs b/src/NumericPrelude/Tuple.hs deleted file mode 100644 --- a/src/NumericPrelude/Tuple.hs +++ /dev/null @@ -1,13 +0,0 @@ -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) diff --git a/test/Test.hs b/test/Test.hs --- a/test/Test.hs +++ b/test/Test.hs @@ -1,4 +1,4 @@ -{-# OPTIONS -fno-implicit-prelude #-} +{-# LANGUAGE NoImplicitPrelude #-} module Main where import Number.Complex((+:), (-:), ) diff --git a/test/Test/MathObj/Gaussian/Bell.hs b/test/Test/MathObj/Gaussian/Bell.hs new file mode 100644 --- /dev/null +++ b/test/Test/MathObj/Gaussian/Bell.hs @@ -0,0 +1,94 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +module Test.MathObj.Gaussian.Bell where + +import qualified MathObj.Gaussian.Bell as G + +-- import qualified Algebra.Ring as Ring + +import qualified Algebra.Laws as Laws + +import qualified Number.Complex as Complex + +import Test.NumericPrelude.Utility (testUnit) +import Test.QuickCheck (Testable, quickCheck, (==>)) +import qualified Test.HUnit as HUnit + +import Data.Function.HT (nest, ) + +import PreludeBase as P +import NumericPrelude as NP + + +simple :: + (Testable t) => + (G.T Rational -> t) -> IO () +simple f = + quickCheck (\x -> f (x :: G.T Rational)) + +tests :: HUnit.Test +tests = + HUnit.TestLabel "polynomial" $ + HUnit.TestList $ + map testUnit $ +{- + ("convolution, dirac", + simple $ Laws.identity (+) zero) : +-} + ("convolution, commutative", + simple $ Laws.commutative G.convolve) : + ("convolution, associative", + simple $ Laws.associative G.convolve) : + ("multiplication, one", + simple $ Laws.identity G.multiply G.constant) : + ("multiplication, commutative", + simple $ Laws.commutative G.multiply) : + ("multiplication, associative", + simple $ Laws.associative G.multiply) : + ("convolution, multplication, fourier", + simple $ \x y -> + G.fourier (G.convolve x y) + == G.multiply (G.fourier x) (G.fourier y)) : + ("convolution via translation", + simple $ \x y -> + G.convolve x y + == G.convolveByTranslation x y) : + ("convolution via fourier", + simple $ \x y -> + G.convolve x y + == G.convolveByFourier x y) : + ("fourier reverse", + simple $ \x -> nest 2 G.fourier x == G.reverse x) : + ("reverse identity", + simple $ \x -> nest 2 G.reverse x == x) : + ("translate additive", + simple $ \x a b -> + G.translate a (G.translate b x) == G.translate (a+b) x) : + ("translateComplex additive", + simple $ \x a b -> + G.translateComplex a (G.translateComplex b x) == G.translateComplex (a+b) x) : + ("translateComplex real", + simple $ \x a -> + G.translateComplex (Complex.fromReal a) x == G.translate a x) : + ("modulate additive", + simple $ \x a b -> + G.modulate a (G.modulate b x) == G.modulate (a+b) x) : + ("commute translate modulate", + simple $ \x a b -> + G.modulate b (G.translate a x) + == G.turn (a*b) (G.translate a (G.modulate b x))) : + ("fourier translate", + simple $ \x a -> + G.fourier (G.translate a x) + == G.modulate a (G.fourier x)) : + ("dilate multiplicative", + simple $ \x a b -> a>0 && b>0 ==> + G.dilate a (G.dilate b x) == G.dilate (a*b) x) : + ("dilate by shrink", + simple $ \x a -> a>0 ==> + G.shrink a x == G.dilate (recip a) x) : + ("fourier dilate", + simple $ \x a -> a>0 ==> + G.fourier (G.dilate a x) == G.shrink a (G.fourier x)) : + [] diff --git a/test/Test/MathObj/Gaussian/Variance.hs b/test/Test/MathObj/Gaussian/Variance.hs new file mode 100644 --- /dev/null +++ b/test/Test/MathObj/Gaussian/Variance.hs @@ -0,0 +1,62 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +module Test.MathObj.Gaussian.Variance where + +import qualified MathObj.Gaussian.Variance as G + +-- import qualified Algebra.Ring as Ring + +import qualified Algebra.Laws as Laws + +import Test.NumericPrelude.Utility (testUnit) +import Test.QuickCheck (Testable, quickCheck, (==>)) +import qualified Test.HUnit as HUnit + +import Data.Function.HT (nest, ) + +import PreludeBase as P +import NumericPrelude as NP + + +simple :: + (Testable t) => + (G.T Rational -> t) -> IO () +simple f = + quickCheck (\x -> f (x :: G.T Rational)) + +tests :: HUnit.Test +tests = + HUnit.TestLabel "polynomial" $ + HUnit.TestList $ + map testUnit $ +{- + ("convolution, dirac", + simple $ Laws.identity (+) zero) : +-} + ("convolution, commutative", + simple $ Laws.commutative G.convolve) : + ("convolution, associative", + simple $ Laws.associative G.convolve) : + ("multiplication, one", + simple $ Laws.identity G.multiply G.constant) : + ("multiplication, commutative", + simple $ Laws.commutative G.multiply) : + ("multiplication, associative", + simple $ Laws.associative G.multiply) : + ("convolution via fourier", + simple $ \x y -> + G.fourier (G.convolve x y) + == G.multiply (G.fourier x) (G.fourier y)) : + ("fourier identity", + simple $ \x -> nest 4 G.fourier x == x) : + ("dilate multiplicative", + simple $ \x a b -> a>0 && b>0 ==> + G.dilate a (G.dilate b x) == G.dilate (a*b) x) : + ("dilate by shrink", + simple $ \x a -> a>0 ==> + G.shrink a x == G.dilate (recip a) x) : + ("fourier dilate", + simple $ \x a -> a>0 ==> + G.fourier (G.dilate a x) == G.shrink a (G.fourier x)) : + [] diff --git a/test/Test/MathObj/PartialFraction.hs b/test/Test/MathObj/PartialFraction.hs --- a/test/Test/MathObj/PartialFraction.hs +++ b/test/Test/MathObj/PartialFraction.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} module Test.MathObj.PartialFraction where import qualified MathObj.PartialFraction as PartialFraction @@ -14,7 +16,7 @@ import qualified Algebra.Laws as Laws import qualified Test.QuickCheck as QC -import NumericPrelude.Monad (untilM) +import Control.Monad.HT as M import Test.NumericPrelude.Utility (testUnit) import Test.QuickCheck (quickCheck) import qualified Test.HUnit as HUnit @@ -134,7 +136,7 @@ instance QC.Arbitrary IrredPoly where arbitrary = do poly <- QC.elements (map Poly.fromCoeffs [[2,3],[2,0,1],[3,0,1],[1,-3,0,1]]) - unit <- untilM (not. isZero) QC.arbitrary + unit <- M.until (not. isZero) QC.arbitrary return (IrredPoly (unit Vector.*> poly)) coarbitrary = undefined diff --git a/test/Test/MathObj/Polynomial.hs b/test/Test/MathObj/Polynomial.hs --- a/test/Test/MathObj/Polynomial.hs +++ b/test/Test/MathObj/Polynomial.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} module Test.MathObj.Polynomial where import qualified MathObj.Polynomial as Poly diff --git a/test/Test/MathObj/PowerSeries.hs b/test/Test/MathObj/PowerSeries.hs --- a/test/Test/MathObj/PowerSeries.hs +++ b/test/Test/MathObj/PowerSeries.hs @@ -1,4 +1,6 @@ -{-# OPTIONS -fno-implicit-prelude -fglasgow-exts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} module Test.MathObj.PowerSeries where import qualified MathObj.PowerSeries as PS diff --git a/test/Test/NumericPrelude/List.hs b/test/Test/NumericPrelude/List.hs deleted file mode 100644 --- a/test/Test/NumericPrelude/List.hs +++ /dev/null @@ -1,85 +0,0 @@ -module Test.NumericPrelude.List where - -import qualified NumericPrelude.List as NList -import qualified Data.List as List -import Control.Monad (liftM2) - -import Test.NumericPrelude.Utility (equalLists, equalInfLists, testUnit) -import Test.QuickCheck (Property, quickCheck, (==>)) -import qualified Test.HUnit as HUnit - - - -sieve :: Eq a => Int -> [a] -> Property -sieve n x = - n>0 ==> - equalLists [NList.sieve n x, - NList.sieve' n x, - NList.sieve'' n x, - NList.sieve''' n x] - - -sliceHoriz :: Eq a => Int -> [a] -> Property -sliceHoriz n x = - n>0 ==> - NList.sliceHoriz n x == NList.sliceHoriz' n x - - -sliceVert :: Eq a => Int -> [a] -> Property -sliceVert n x = - n>0 ==> - NList.sliceVert n x == NList.sliceVert' n x - -slice :: Eq a => Int -> [a] -> Property -slice n x = - 0 - -- problems: NList.sliceHoriz 4 [] == [[],[],[],[]] - NList.sliceHoriz n x == List.transpose (NList.sliceVert n x) && - NList.sliceVert n x == List.transpose (NList.sliceHoriz n x) - - - - -shear :: Eq a => [[a]] -> Bool -shear xs = - NList.shearTranspose xs == map reverse (NList.shear xs) - - - -outerProduct :: (Eq a, Eq b) => [a] -> [b] -> Bool -outerProduct xs ys = - equalLists [concat (NList.outerProduct (,) xs ys), liftM2 (,) xs ys] - - - -reduceRepeated :: Eq a => - (a -> a -> a) -> a -> a -> Integer -> Property -reduceRepeated op a0 a n = - n>0 ==> - NList.reduceRepeated op a0 a n == NList.reduceRepeatedSlow op a0 a n - - -iterate' :: Eq a => (a -> a -> a) -> a -> Bool -iterate' op a = - let xs = List.iterate (op a) a - ys = NList.iterateAssoc op a - zs = NList.iterateLeaky op a - in equalInfLists 1000 [xs, ys, zs] - - - - -tests :: HUnit.Test -tests = - HUnit.TestLabel "list" $ - HUnit.TestList $ - map testUnit $ - ("sieve", quickCheck (sieve :: Int -> [Integer] -> Property)) : - ("sliceHoriz", quickCheck (sliceHoriz :: Int -> [Integer] -> Property)) : - ("sliceVert", quickCheck (sliceVert :: Int -> [Integer] -> Property)) : - ("slice", quickCheck (slice :: Int -> [Integer] -> Property)) : - ("shear", quickCheck (shear :: [[Integer]] -> Bool)) : - ("outerProduct", quickCheck (outerProduct :: [Integer] -> [Int] -> Bool)) : - ("reduceRepeated", quickCheck (reduceRepeated (+) :: Integer -> Integer -> Integer -> Property)) : - ("iterate", quickCheck (iterate' (+) :: Integer -> Bool)) : - [] diff --git a/test/Test/NumericPrelude/Utility.hs b/test/Test/NumericPrelude/Utility.hs --- a/test/Test/NumericPrelude/Utility.hs +++ b/test/Test/NumericPrelude/Utility.hs @@ -1,5 +1,7 @@ +-- cf. utility-ht Test.Utility module Test.NumericPrelude.Utility where +import Data.List.HT (mapAdjacent, ) import qualified Data.List as List import qualified Test.HUnit as HUnit @@ -8,11 +10,11 @@ testUnit (label, check) = HUnit.TestLabel label (HUnit.TestCase check) - +-- compare the lists simultaneously equalLists :: Eq a => [[a]] -> Bool equalLists xs = let equalElems ys = - and (zipWith (==) ys (tail ys)) && length xs == length ys + and (mapAdjacent (==) ys) && length xs == length ys in all equalElems (List.transpose xs) equalInfLists :: Eq a => Int -> [[a]] -> Bool diff --git a/test/Test/Run.hs b/test/Test/Run.hs --- a/test/Test/Run.hs +++ b/test/Test/Run.hs @@ -1,6 +1,7 @@ module Main where -import qualified Test.NumericPrelude.List as NList +import qualified Test.MathObj.Gaussian.Variance as GaussVariance +import qualified Test.MathObj.Gaussian.Bell as GaussBell import qualified Test.MathObj.PartialFraction as PartialFraction import qualified Test.MathObj.Polynomial as Polynomial import qualified Test.MathObj.PowerSeries as PowerSeries @@ -10,7 +11,8 @@ main :: IO () main = do HUnitText.runTestTT (HUnit.TestList $ - NList.tests : + GaussVariance.tests : + GaussBell.tests : PartialFraction.tests : Polynomial.tests : PowerSeries.tests :