packages feed

numeric-prelude 0.1.3.4 → 0.2

raw patch · 106 files changed

+3250/−1698 lines, 106 filesdep ~gnuplotdep ~non-negative

Dependency ranges changed: gnuplot, non-negative

Files

docs/NOTES view
@@ -1,5 +1,16 @@-** -> ^?+* sum (and mconcat) +How to provide a 'sum' function that works optimal for the strict and lazy types?+It must sum strict types from the left and lazy types from the right,+and pairs of both kind of types must be added in a mixed manner.+This seems to be impossible to achieve with a Haskell 98 type class or optimizer rules.++http://projects.haskell.org/pipermail/numeric-prelude/2010-July/000016.html++Shall we leave left-biased sum as the default+and provide an advanced type class for unbounded summation?+It could be used to define a generic 'List.length' function.+ * non-negative  NonNegative could require ZeroTestable@@ -63,26 +74,6 @@ but not a locally defined 'fromRational'.  -* people probably interested in NumPrelude:-- Mike Thomas <miketh@brisbane.paradigmgeo.com>-   http://www.haskell.org/pipermail/haskell-cafe/2002-February/002660.html-- jan.skibinski@sympatico.ca-   indexless linear algebra-- blaetterrascheln@web.de- Christian Sievers <sievers@math2.nat.tu-bs.de>- Remi Turk <buran@xs4all.nl>, rturk@science.uva.nl- Ronny Wichers Schreur <R.WichersSchreur@science.ru.nl>-   floorSqrt- - William Lee Irwin III <wli@holomorphy.com>-   ContFrac, continued fractions-- Juergen Bokowski <bokowski@mathematik.tu-darmstadt.de>-   DMV-Nachrichten 2004/3- * RealFloat Defines the properties of a Floating type, thus should be named 'Floating'.@@ -149,7 +140,7 @@ If they can't assert that (I assume that will only rarely be the case), they must do this check by themselve. -* Numeric type classes for DSLs+* Numeric type classes for EDSLs  It is very common to define instances of Numeric type classes for wrapping operations of a foreign programming language.@@ -179,6 +170,21 @@ although usage of Eq is discouraged, and Ord is of restricted use. (For similar values, the rounding errors might be greater than the difference of the values.)++The EDSL classes also have its use for CPU vectors+   EqualityDecision: (==)+   OrderDecision: compare+   Choice: ifThenElse+   Share: Verallgemeinerung von 'let'++And there could be more classes for vector computing:+   Shift+   Access: insert, extract++Using this class we could write synthesizer:Wave functions+or vectorised frequency filters,+such that we can test it with exact arithmetic+and run it with LLVM on CPU vectors.  * Implicit configuration
numeric-prelude.cabal view
@@ -1,5 +1,5 @@ Name:           numeric-prelude-Version:        0.1.3.4+Version:        0.2 License:        GPL License-File:   LICENSE Author:         Dylan Thurston <dpt@math.harvard.edu>, Henning Thielemann <numericprelude@henning-thielemann.de>, Mikael Johansson@@ -71,7 +71,30 @@   Kowalczyk, Ketil Malde, Tom Schrijvers, Ken Shan, and Henning   Thielemann for helpful comments.   .+  .+  Usage:+  .+  Write modules in the following style:+  .+  > [-# NoImplicitPrelude #-]+  > module MyModule where+  >+  > ... various specific imports ...+  >+  > import NumericPrelude+  .+  Importing @NumericPrelude@ is almost the same as+  .+  > import NumericPrelude.Numeric+  > import NumericPrelude.Base   .+  .+  Instead of the @NoImplicitPrelude@ pragma+  you could also write @import Prelude ()@+  but this will yield problems with numeric literals.+  .+  .   Scope & Limitations\/TODO:+  .   * It might be desireable to split Ord up into Poset and Ord     (a total ordering).     This is not addressed here.@@ -97,16 +120,13 @@   * I stuck with the Haskell 98 names.     In some cases I find them lacking.     Neglecting backwards compatibility, we have renamed classes as follows:-      Num           --> Ring,-      Fractional    --> Field,-      Floating      --> Algebraic + Transcendental,-      RealFloat     --> RealTranscendental,-  .-  * It's slightly unfortunate that 'abs' can no longer be used for complex numbers,-    since it is standard mathematically.-    'magnitude' or more generally 'Algebra.NormedSpace.Euclidean.norm' can be used.-    But it had the wrong type before,-    and I couldn't see how to fit it in without complicating the hierarchy.+      Num           --> Additive, Ring, Absolute+      Integral      --> ToInteger, IntegralDomain, RealIntegral+      Fractional    --> Field+      Floating      --> Algebraic, Transcendental+      Real          --> ToRational+      RealFrac      --> RealRing, RealField+      RealFloat     --> RealTranscendental   .   .   Additional standard libraries might include Enum, IEEEFloat (including@@ -130,7 +150,7 @@   default:     False  Source-Repository this-  Tag:         0.1.3.4+  Tag:         0.2   Type:        darcs   Location:    http://code.haskell.org/numeric-prelude/ @@ -143,7 +163,7 @@     parsec >=1 && <4,     QuickCheck >=1 && <3,     storable-record >=0.0.1 && <0.1,-    non-negative >=0.0.2 && <0.1,+    non-negative >=0.0.5 && <0.2,     utility-ht >=0.0.4 && <0.1   If flag(splitBase)     Build-Depends:@@ -157,6 +177,7 @@   GHC-Options:    -Wall   Hs-source-dirs: src   Exposed-modules:+    Algebra.Absolute     Algebra.Additive     Algebra.Algebraic     Algebra.Differential@@ -176,9 +197,9 @@     Algebra.NormedSpace.Sum     Algebra.OccasionallyScalar     Algebra.PrincipalIdealDomain-    Algebra.Real     Algebra.RealField     Algebra.RealIntegral+    Algebra.RealRing     Algebra.RealTranscendental     Algebra.RightModule     Algebra.Ring@@ -200,14 +221,17 @@     MathObj.Permutation.CycleList.Check     MathObj.Permutation.Table     MathObj.Polynomial+    MathObj.Polynomial.Core     MathObj.PowerSeries+    MathObj.PowerSeries.Core     MathObj.PowerSeries.DifferentialEquation     MathObj.PowerSeries.Example     MathObj.PowerSeries.Mean     MathObj.PowerSeries2+    MathObj.PowerSeries2.Core     MathObj.PowerSum+    MathObj.RefinementMask2     MathObj.RootSet-    MyPrelude     Number.Complex     Number.DimensionTerm     Number.DimensionTerm.SI@@ -235,15 +259,22 @@     Number.Physical     Number.Physical.Read     Number.Physical.Show+    NumericPrelude.List.Checked+    NumericPrelude.List.Generic     NumericPrelude.Elementwise+    NumericPrelude.Numeric+    NumericPrelude.Base     NumericPrelude-    PreludeBase   Other-modules:     NumericPrelude.List     Algebra.AffineSpace     MathObj.Gaussian.Variance     MathObj.Gaussian.Bell     MathObj.Gaussian.Polynomial+    -- I think I won't add them this way.+    -- It is certainly better to split the class into comparison and selection.+    Algebra.EqualityDecision+    Algebra.OrderDecision  Executable test   Hs-Source-Dirs: src, test@@ -257,6 +288,7 @@   Other-modules:     Test.NumericPrelude.Utility     Test.Number.GaloisField2p32m5+    Test.MathObj.RefinementMask2     Test.MathObj.PartialFraction     Test.MathObj.Matrix     Test.MathObj.Polynomial@@ -277,7 +309,7 @@     MathObj.Gaussian.Example   If flag(buildTests)     Build-Depends:-      gnuplot >=0.3 && <0.4,+      gnuplot >=0.3 && <0.5,       HTam >=0.0.2 && <0.1   Else     Buildable: False
+ src/Algebra/Absolute.hs view
@@ -0,0 +1,151 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Algebra.Absolute (+   C(abs, signum),+   absOrd, signumOrd,+   ) where++import qualified Algebra.Ring         as Ring+import qualified Algebra.Additive     as Additive+import qualified Algebra.ZeroTestable as ZeroTestable++import Algebra.Ring (one, ) -- fromInteger+import Algebra.Additive (zero, negate,)++import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )+import Data.Word (Word, Word8, Word16, Word32, Word64, )++import NumericPrelude.Base+import qualified Prelude as P+import Prelude(Int,Integer,Float,Double)+++{- |+This is the type class of a ring with a notion of an absolute value,+satisfying the laws++>                        a * b === b * a+>   a /= 0  =>  abs (signum a) === 1+>             abs a * signum a === a++Minimal definition: 'abs', 'signum'.++If the type is in the 'Ord' class+we expect 'abs' = 'absOrd' and 'signum' = 'signumOrd'+and we expect the following laws to hold:++>      a + (max b c) === max (a+b) (a+c)+>   negate (max b c) === min (negate b) (negate c)+>      a * (max b c) === max (a*b) (a*c) where a >= 0+>           absOrd a === max a (-a)++We do not require 'Ord' as superclass+since we also want to have "Number.Complex" as instance.+'abs' for complex numbers alone may have an inappropriate type,+because it does not reflect that the absolute value is a real number.+You might prefer 'Number.Complex.magnitude'.+This type class is intended for unifying algorithms+that work for both real and complex numbers.+Note the similarity to "Algebra.Units":+'abs' plays the role of @stdAssociate@+and 'signum' plays the role of @stdUnit@.++Actually, since 'abs' can be defined using 'max' and 'negate'+we could relax the superclasses to @Additive@ and 'Ord'+if his class would only contain 'signum'.+-}+class (Ring.C a, ZeroTestable.C a) => C a where+    abs    :: a -> a+    signum :: a -> a+++absOrd :: (Additive.C a, Ord a) => a -> a+absOrd x = max x (negate x)++signumOrd :: (Ring.C a, Ord a) => a -> a+signumOrd x =+   case compare x zero of+      GT ->        one+      EQ ->        zero+      LT -> negate one+++instance C Integer where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Float   where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Double  where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum+++instance C Int     where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Int8    where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Int16   where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Int32   where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Int64   where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum+++instance C Word    where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Word8   where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Word16  where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Word32  where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum++instance C Word64  where+   {-# INLINE abs #-}+   {-# INLINE signum #-}+   abs = P.abs+   signum = P.signum+
src/Algebra/Additive.hs view
@@ -32,7 +32,7 @@ import qualified Data.Ratio as Ratio98 import qualified Prelude as P import Prelude(Int, Integer, Float, Double, fromInteger, )-import PreludeBase+import NumericPrelude.Base   infixl 6  +, -@@ -80,6 +80,10 @@ {- | Sum up all elements of a list. An empty list yields zero.++This function is inappropriate for number types like Peano.+Maybe we should make 'sum' a method of Additive.+This would also make 'lengthLeft' and 'lengthRight' superfluous. -} sum :: (C a) => [a] -> a sum = foldl (+) zero@@ -130,18 +134,21 @@ > addPair :: (Additive.C a, Additive.C b) => (a,b) -> (a,b) -> (a,b) > addPair = Elem.run2 $ Elem.with (,) <*>.+  fst <*>.+  snd -}+{-# INLINE (<*>.+) #-} (<*>.+) ::    (C x) =>    Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a (<*>.+) f acc =    f <*> elementAdd acc +{-# INLINE (<*>.-) #-} (<*>.-) ::    (C x) =>    Elem.T (v,v) (x -> a) -> (v -> x) -> Elem.T (v,v) a (<*>.-) f acc =    f <*> elementSub acc +{-# INLINE (<*>.-$) #-} (<*>.-$) ::    (C x) =>    Elem.T v (x -> a) -> (v -> x) -> Elem.T v a
src/Algebra/AffineSpace.hs view
@@ -33,8 +33,8 @@  import Control.Applicative (Applicative(pure, (<*>)), ) -import NumericPrelude hiding (zero, )-import PreludeBase+import NumericPrelude.Numeric hiding (zero, )+import NumericPrelude.Base import Prelude ()  {- |
src/Algebra/Algebraic.hs view
@@ -12,7 +12,7 @@ import Algebra.Ring ((*), (^), fromInteger) import Algebra.Additive((+)) -import PreludeBase+import NumericPrelude.Base import qualified Prelude as P
src/Algebra/Differential.hs view
@@ -3,7 +3,7 @@  import qualified Algebra.Ring as Ring --- import NumericPrelude+-- import NumericPrelude.Numeric import qualified Prelude  {- |
+ src/Algebra/EqualityDecision.hs view
@@ -0,0 +1,110 @@+{- |+Combination of @(==)@ and @if then else@+that can be instantiated for more types than @Eq@+or can be instantiated in a way+that allows more defined results (\"more total\" functions):++* Reader like types for representing a context+  like 'Number.ResidueClass.Reader'++* Expressions in an EDSL++* More generally every type based on an applicative functor++* Tuples and Vector types++* Positional and Peano numbers,+  a common prefix of two numbers can be emitted+  before the comparison is done.+  (This could also be done with an overloaded 'if',+   what we also do not have.)+-}+module Algebra.EqualityDecision where++import qualified NumericPrelude.Elementwise as Elem+import Control.Applicative (Applicative(pure, (<*>)), )+import Data.Tuple.HT (fst3, snd3, thd3, )+import Data.List (zipWith4, )+++{- |+For atomic types this could be a superclass of 'Eq'.+However for composed types like tuples, lists, functions+we do elementwise comparison+which is incompatible with the complete comparison performed by '(==)'.+-}+class C a where+   {- |+   It holds++   > (a ==? b) eq noteq  ==  if a==b then eq else noteq++   for atomic types where the right hand side can be defined.+   -}+   (==?) :: a -> a -> a -> a -> a++++{-# INLINE deflt #-}+deflt :: Eq a => a -> a -> a -> a -> a+deflt a b eq noteq =+   if a==b then eq else noteq++++instance C Int where+   {-# INLINE (==?) #-}+   (==?) = deflt++instance C Integer where+   {-# INLINE (==?) #-}+   (==?) = deflt++instance C Float where+   {-# INLINE (==?) #-}+   (==?) = deflt++instance C Double where+   {-# INLINE (==?) #-}+   (==?) = deflt++instance C Bool where+   {-# INLINE (==?) #-}+   (==?) = deflt++instance C Ordering where+   {-# INLINE (==?) #-}+   (==?) = deflt++++{-# INLINE element #-}+element ::+   (C x) =>+   (v -> x) -> Elem.T (v,v,v,v) x+element f =+   Elem.element (\(x,y,eq,noteq) -> (f x ==? f y) (f eq) (f noteq))++{-# INLINE (<*>.==?) #-}+(<*>.==?) ::+   (C x) =>+   Elem.T (v,v,v,v) (x -> a) -> (v -> x) -> Elem.T (v,v,v,v) a+(<*>.==?) f acc =+   f <*> element acc+++instance (C a, C b) => C (a,b) where+   {-# INLINE (==?) #-}+   (==?) = Elem.run4 $ pure (,) <*>.==?  fst <*>.==?  snd++instance (C a, C b, C c) => C (a,b,c) where+   {-# INLINE (==?) #-}+   (==?) = Elem.run4 $ pure (,,) <*>.==?  fst3 <*>.==?  snd3 <*>.==?  thd3++instance C a => C [a] where+   {-# INLINE (==?) #-}+   (==?) = zipWith4 (==?)++instance (C a) => C (b -> a) where+   {-# INLINE (==?) #-}+   (==?) x y eq noteq c  =  (x c ==? y c) (eq c) (noteq c)
src/Algebra/Field.hs view
@@ -27,7 +27,7 @@ import Algebra.Additive (zero, negate) import Algebra.ZeroTestable (isZero) -import PreludeBase+import NumericPrelude.Base import Prelude (Integer, Float, Double) import qualified Prelude as P import Test.QuickCheck ((==>), Property)
src/Algebra/GenerateRules.hs view
@@ -45,11 +45,11 @@ realFieldIndirect :: [String] realFieldIndirect =    do targetType <- tail machineIntegerTypes-      method <- "round" : "truncate" : "floor" : "ceiling" : []-      let methodPad = pad 8 method+      method <- "round" : "roundSimple" : "truncate" : "floor" : "ceiling" : []+      let methodPad = pad 11 method       let signature = functionSignature methodPad "a" targetType       return $ "     " ++-         pad 30 ("\"NP." ++ signature ++ "\"") +++         pad 33 ("\"NP." ++ signature ++ "\"") ++          methodPad ++ " = (" ++ functionSignature "P.fromIntegral" "Int" targetType ++ ") . "              ++ method ++ ";" @@ -78,7 +78,7 @@  main :: IO () main =-   putStrLn "module Algebra.RealField" >>+   putStrLn "module Algebra.RealRing" >>    mapM_ putStrLn realFieldIndirect >>    mapM_ putStrLn splitFractionIndirect >>
src/Algebra/IntegralDomain.hs view
@@ -43,7 +43,7 @@ import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  ) import Data.Word (Word, Word8, Word16, Word32, Word64, ) -import PreludeBase+import NumericPrelude.Base import Prelude (Integer, Int) import qualified Prelude as P
src/Algebra/Lattice.hs view
@@ -13,8 +13,8 @@  import qualified Algebra.Laws as Laws -import NumericPrelude hiding (abs)-import PreludeBase hiding (max, min)+import NumericPrelude.Numeric hiding (abs)+import NumericPrelude.Base hiding (max, min) import qualified Prelude as P  infixl 5 `up`, `dn`
src/Algebra/Module.hs view
@@ -60,6 +60,7 @@     (*>) :: a -> v -> v  +{-# INLINE (<*>.*>) #-} (<*>.*>) ::    (C a x) =>    Elem.T (a,v) (x -> c) -> (v -> x) -> Elem.T (a,v) c
src/Algebra/Monoid.hs view
@@ -1,6 +1,5 @@-{-# LANGUAGE NoImplicitPrelude #-} {- |-Copyright    :   (c) Henning Thielemann 2009, Mikael Johansson 2006+Copyright    :   (c) Henning Thielemann 2009-2010, Mikael Johansson 2006 Maintainer   :   numericprelude@henning-thielemann.de Stability    :   provisional Portability  :@@ -27,36 +26,47 @@ class C a where   idt   :: a   (<*>) :: a -> a -> a+  cumulate :: [a] -> a+  cumulate = foldr (<*>) idt + instance C All where   idt = mempty   (<*>) = mappend+  cumulate = mconcat  instance C Any where   idt = mempty   (<*>) = mappend+  cumulate = mconcat  instance C a => C (Dual a) where   idt = Mn.Dual idt   (Mn.Dual x) <*> (Mn.Dual y) = Mn.Dual (y <*> x)+  cumulate = Mn.Dual . cumulate . reverse . map Mn.getDual  instance C (Endo a) where   idt = mempty   (<*>) = mappend+  cumulate = mconcat  instance C (First a) where   idt = mempty   (<*>) = mappend+  cumulate = mconcat  instance C (Last a) where   idt = mempty   (<*>) = mappend+  cumulate = mconcat   instance Ring.C a => C (Product a) where   idt = Mn.Product Ring.one   (Mn.Product x) <*> (Mn.Product y) = Mn.Product (x Ring.* y)+  cumulate = Mn.Product . Ring.product . map Mn.getProduct  instance Additive.C a => C (Sum a) where   idt = Mn.Sum Additive.zero   (Mn.Sum x) <*> (Mn.Sum y) = Mn.Sum (x Additive.+ y)+  cumulate = Mn.Sum . Additive.sum . map Mn.getSum
src/Algebra/NonNegative.hs view
@@ -1,43 +1,130 @@ {- |-Copyright   :  (c) Henning Thielemann 2007+Copyright   :  (c) Henning Thielemann 2007-2010  Maintainer  :  haskell@henning-thielemann.de Stability   :  stable Portability :  Haskell 98  A type class for non-negative numbers.-Prominent instances are 'Numeric.NonNegative.Wrapper.T' and peano numbers.+Prominent instances are 'Number.NonNegative.T' and 'Number.Peano.T' numbers. This class cannot do any checks, but it let you show to the user what arguments your function expects.+Thus you must define class instances with care. In fact many standard functions ('take', '(!!)', ...) should have this type class constraint.-Thus you must define class instances with care. -}-module Algebra.NonNegative (C(..)) where+module Algebra.NonNegative (+   C(..),+   splitDefault, +   (-|),+--   (-?),+   zero,+   add,+   sum,+   ) where+ import qualified Algebra.Additive as Additive-import qualified Algebra.Real     as Real+-- import qualified Algebra.RealRing as RealRing -infixl 6 -|, -?+import qualified Algebra.Monoid as Monoid +-- import Algebra.Absolute (abs, )+import Algebra.Additive ((-), )++import Prelude hiding (sum, (-), abs, )+++infixl 6 -|  -- , -?++ {- | Instances of this class must ensure non-negative values. We cannot enforce this by types, but the type class constraint @NonNegative.C@ avoids accidental usage of types which allow for negative numbers.++The Monoid superclass contributes a zero and an addition. -}-class (Ord a, Additive.C a) => C a where+class (Ord a, Monoid.C a) => C a where    {- |-   @x -| y == max 0 (x-y)@+   @split x y == (m,(b,d))@ means that+   @b == (x<=y)@,+   @m == min x y@,+   @d == max x y - min x y@, that is @d == abs(x-y)@. -   The default implementation is not efficient,-   because it compares the values and then subtracts, again, if safe.-   @max 0 (x-y)@ is more elegant and efficient-   but not possible in the general case,-   since @x-y@ may already yield a negative number.+   We have chosen this function as base function,+   since it provides comparison and subtraction in one go,+   which is important for replacing common structures like++   > if x<=y+   >   then f(x-y)+   >   else g(y-x)++   that lead to a memory leak for peano numbers.+   We have choosen the simple check @x<=y@+   instead of a full-blown @compare@,+   since we want @Zero <= undefined@ for peano numbers.+   Because of undefined values 'split' is in general+   not commutative in the sense++   > let (m0,(b0,d0)) = split x y+   >     (m1,(b1,d1)) = split y x+   > in  m0==m1 && d0==d1++   The result values are in the order+   in which they are generated for Peano numbers.+   We have chosen the nested pair instead of a triple+   in order to prevent a memory leak+   that occurs if you only use @b@ and @d@ and ignore @m@.+   This is demonstrated by test cases+   Chunky.splitSpaceLeak3 and Chunky.splitSpaceLeak4.    -}-   (-|) :: a -> a -> a-   x -| y  =  if x >= y then x Additive.- y else Additive.zero+   split :: a -> a -> (a, (Bool, a))  -(-?) :: (Real.C a) => a -> a -> (Bool, a)-(-?) x y  =  (x >= y, Real.abs (x Additive.- y))+{- |+Default implementation for wrapped types of 'Ord' and 'Num' class.+-}+{-# INLINE splitDefault #-}+splitDefault ::+   (Ord b, Additive.C b) =>+   (a -> b) -> (b -> a) -> a -> a -> (a, (Bool, a))+splitDefault unpack pack px py =+   let x = unpack px+       y = unpack py+   in  if x<=y+         then (pack x, (True,  pack (y-x)))+         else (pack y, (False, pack (x-y)))+++zero :: C a => a+zero = Monoid.idt++-- like (+)+infixl 6 `add`++add :: C a => a -> a -> a+add = (Monoid.<*>)++sum :: C a => [a] -> a+sum = Monoid.cumulate+++{- |+@x -| y == max 0 (x-y)@++The default implementation is not efficient,+because it compares the values and then subtracts, again, if safe.+@max 0 (x-y)@ is more elegant and efficient+but not possible in the general case,+since @x-y@ may already yield a negative number.+-}+(-|) :: C a => a -> a -> a+x -| y  =+   let (b,d) = snd $ split y x+   in  if b then d else zero++{-+(-?) :: (RealRing.C a) => a -> a -> (Bool, a)+(-?) x y  =  snd $ split y x+-}
src/Algebra/NormedSpace/Euclidean.hs view
@@ -3,7 +3,7 @@ {-# LANGUAGE FlexibleInstances #-}  {- |-Copyright   :  (c) Henning Thielemann 2005+Copyright   :  (c) Henning Thielemann 2005-2010 License     :  GPL  Maintainer  :  numericprelude@henning-thielemann.de@@ -15,31 +15,60 @@  module Algebra.NormedSpace.Euclidean where -import PreludeBase-import NumericPrelude (sqr, abs, (+), sum, Float, Double, Int, Integer, )+import NumericPrelude.Base+import NumericPrelude.Numeric (sqr, abs, zero, (+), sum, Float, Double, Int, Integer, )  import qualified Number.Ratio as Ratio  import qualified Algebra.PrincipalIdealDomain as PID import qualified Algebra.Algebraic as Algebraic-import qualified Algebra.Real      as Real+import qualified Algebra.Absolute      as Absolute import qualified Algebra.Module    as Module +import qualified Data.Foldable as Fold++ {-|-A vector space equipped with an Euclidean or a Hilbert norm.+Helper class for 'C' that does not need an algebraic type @a@.  Minimal definition: 'normSqr' -}-class (Real.C a, Module.C a v) => Sqr a v where+class (Absolute.C a, Module.C a v) => Sqr a v where   {-| Square of the Euclidean norm of a vector.       This is sometimes easier to implement. -}   normSqr :: v -> a --  normSqr = sqr . norm +{- |+Default definition for 'normSqr' that is based on 'Fold.Foldable' class.+-}+{-# INLINE normSqrFoldable #-}+normSqrFoldable ::+   (Sqr a v, Fold.Foldable f) => f v -> a+normSqrFoldable =+   Fold.foldl (\a v -> a + normSqr v) zero++{- |+Default definition for 'normSqr' that is based on 'Fold.Foldable' class+and the argument vector has at least one component.+-}+{-# INLINE normSqrFoldable1 #-}+normSqrFoldable1 ::+   (Sqr a v, Fold.Foldable f, Functor f) => f v -> a+normSqrFoldable1 =+   Fold.foldl1 (+) . fmap normSqr+++{-|+A vector space equipped with an Euclidean or a Hilbert norm.++Minimal definition:+'norm'+-} class (Sqr a v) => C a v where   {-| Euclidean norm of a vector. -}-  norm    :: v -> a+  norm :: v -> a   defltNorm :: (Algebraic.C a, Sqr a v) => v -> a@@ -75,7 +104,7 @@  {-* Instances for composed types -} -instance (Real.C a, PID.C a) => Sqr (Ratio.T a) (Ratio.T a) where+instance (Absolute.C a, PID.C a) => Sqr (Ratio.T a) (Ratio.T a) where   normSqr = sqr  instance (Sqr a v0, Sqr a v1) => Sqr a (v0, v1) where
src/Algebra/NormedSpace/Maximum.hs view
@@ -3,7 +3,7 @@ {-# LANGUAGE FlexibleInstances #-}  {- |-Copyright   :  (c) Henning Thielemann 2005+Copyright   :  (c) Henning Thielemann 2005-2010 License     :  GPL  Maintainer  :  numericprelude@henning-thielemann.de@@ -15,18 +15,42 @@  module Algebra.NormedSpace.Maximum where -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  import qualified Number.Ratio as Ratio  import qualified Algebra.PrincipalIdealDomain as PID-import qualified Algebra.Real     as Real+import qualified Algebra.ToInteger as ToInteger+import qualified Algebra.RealRing as RealRing+import qualified Algebra.Absolute as Absolute import qualified Algebra.Module   as Module -class (Real.C a, Module.C a v) => C a v where+import qualified Data.Foldable as Fold+++class (RealRing.C a, Module.C a v) => C a v where   norm :: v -> a +{- |+Default definition for 'norm' that is based on 'Fold.Foldable' class.+-}+{-# INLINE normFoldable #-}+normFoldable ::+   (C a v, Fold.Foldable f) => f v -> a+normFoldable =+   Fold.foldl (\a v -> max a (norm v)) zero++{- |+Default definition for 'norm' that is based on 'Fold.Foldable' class+and the argument vector has at least one component.+-}+{-# INLINE normFoldable1 #-}+normFoldable1 ::+   (C a v, Fold.Foldable f, Functor f) => f v -> a+normFoldable1 =+   Fold.foldl1 max . fmap norm+ {- instance (Ring.C a, Algebra.Module a a) => C a a where   norm = abs@@ -44,7 +68,7 @@   norm = abs  -instance (Real.C a, PID.C a) => C (Ratio.T a) (Ratio.T a) where+instance (RealRing.C a, ToInteger.C a, PID.C a) => C (Ratio.T a) (Ratio.T a) where   norm = abs  instance (Ord a, C a v0, C a v1) => C a (v0, v1) where@@ -55,4 +79,8 @@  instance (Ord a, C a v) => C a [v] where   norm = foldl max zero . map norm---  norm = maximum . map norm+{-+Since the norm is always non-negative,+we can use zero as identity element.+  norm = maximum . map norm+-}
src/Algebra/NormedSpace/Sum.hs view
@@ -3,7 +3,7 @@ {-# LANGUAGE FlexibleInstances #-}  {- |-Copyright   :  (c) Henning Thielemann 2005+Copyright   :  (c) Henning Thielemann 2005-2010 License     :  GPL  Maintainer  :  numericprelude@henning-thielemann.de@@ -15,16 +15,19 @@  module Algebra.NormedSpace.Sum where -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  import qualified Number.Ratio as Ratio  import qualified Algebra.PrincipalIdealDomain as PID-import qualified Algebra.Real     as Real+import qualified Algebra.Absolute     as Absolute import qualified Algebra.Additive as Additive import qualified Algebra.Module   as Module +import qualified Data.Foldable as Fold++ {-|   The super class is only needed to state the laws   @@@ -33,9 +36,29 @@      norm (u+v)       <=   norm u + norm v   @ -}-class (Real.C a, Module.C a v) => C a v where+class (Absolute.C a, Module.C a v) => C a v where   norm :: v -> a +{- |+Default definition for 'norm' that is based on 'Fold.Foldable' class.+-}+{-# INLINE normFoldable #-}+normFoldable ::+   (C a v, Fold.Foldable f) => f v -> a+normFoldable =+   Fold.foldl (\a v -> a + norm v) zero++{- |+Default definition for 'norm' that is based on 'Fold.Foldable' class+and the argument vector has at least one component.+-}+{-# INLINE normFoldable1 #-}+normFoldable1 ::+   (C a v, Fold.Foldable f, Functor f) => f v -> a+normFoldable1 =+   Fold.foldl1 (+) . fmap norm++ {- instance (Ring.C a, Algebra.Module a a) => C a a where   norm = abs@@ -54,7 +77,7 @@   norm = abs  -instance (Real.C a, PID.C a) => C (Ratio.T a) (Ratio.T a) where+instance (Absolute.C a, PID.C a) => C (Ratio.T a) (Ratio.T a) where   norm = abs  instance (Additive.C a, C a v0, C a v1) => C a (v0, v1) where
src/Algebra/OccasionallyScalar.hs view
@@ -27,7 +27,7 @@  module Algebra.OccasionallyScalar where --- import qualified Algebra.RealField    as RealField+-- import qualified Algebra.RealRing    as RealRing import qualified Algebra.ZeroTestable as ZeroTestable import qualified Algebra.Additive     as Additive import qualified Number.Complex       as Complex@@ -36,8 +36,8 @@  import Number.Complex((+:)) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric   -- this is somehow similar to Normalized classes@@ -76,7 +76,7 @@    fromScalar x    = fromScalar x +: zero  {- converting values automatically to integers is a bad idea-instance (Integral b, RealField.C a)+instance (Integral b, RealRing.C a)       => C b a where    toScalar        = toScalarDefault    toMaybeScalar x = mapMaybe round (toMaybeScalar x)
+ src/Algebra/OrderDecision.hs view
@@ -0,0 +1,244 @@+{- |+Combination of @compare@ and @if then else@+that can be instantiated for more types than @Ord@+or can be instantiated in a way+that allows more defined results (\"more total\" functions):++* Reader like types for representing a context+  like 'Number.ResidueClass.Reader'++* Expressions in an EDSL++* More generally every type based on an applicative functor++* Tuples and Vector types++* Positional and Peano numbers,+  a common prefix of two numbers can be emitted+  before the comparison is done.+  (This could also be done with an overloaded 'if',+   what we also do not have.)+-}+module Algebra.OrderDecision where++import qualified Algebra.EqualityDecision as Equality+import Algebra.EqualityDecision ((==?), )++import qualified NumericPrelude.Elementwise as Elem+import Control.Applicative (Applicative(pure, (<*>)), )+import Data.Tuple.HT (fst3, snd3, thd3, )+import Data.List (zipWith4, zipWith5, )++import Prelude hiding (compare, min, max, minimum, maximum, )+import qualified Prelude as P++++{- |+For atomic types this could be a superclass of 'Ord'.+However for composed types like tuples, lists, functions+we do elementwise comparison+which is incompatible with the complete comparison performed by 'P.compare'.+-}+class Equality.C a => C a where+   {- |+   It holds++   > (compare a b) lt eq gt  ==+   >    case Prelude.compare a b of+   >       LT -> lt+   >       EQ -> eq+   >       GT -> gt++   for atomic types where the right hand side can be defined.++   Minimal complete definition:+   'compare' or '(<=?)'.+   -}+   compare :: a -> a -> a -> a -> a -> a+   compare x y lt eq gt =+      (x ==? y) eq ((x <=? y) lt gt)++   {-# INLINE (<=?) #-}+   (<=?) :: a -> a -> a -> a -> a+   (<=?) x y le gt =+      compare x y le le gt++   {-# INLINE (>=?) #-}+   (>=?) :: a -> a -> a -> a -> a+   (>=?) = flip (<=?)++   (<?) :: a -> a -> a -> a -> a+   (<?) x y = flip (x >=? y)++   {-# INLINE (>?) #-}+   (>?) :: a -> a -> a -> a -> a+   (>?) = flip (<?)++{-+   (<?) :: a -> a -> a -> a -> a+   (<?) x y lt ge =+      compare x y lt ge ge++   (>?) :: a -> a -> a -> a -> a+   (>?) x y gt le =+      compare x y le le gt++   (<=?) :: a -> a -> a -> a -> a+   (<=?) x y le gt =+      compare x y le le gt++   (>=?) :: a -> a -> a -> a -> a+   (>=?) x y ge lt =+      compare x y lt ge ge+-}+++max :: C a => a -> a -> a+max x y = (x>?y) x y++min :: C a => a -> a -> a+min x y = (x<?y) x y++maximum :: C a => a -> [a] -> a+maximum x xs = foldl max x xs++minimum :: C a => a -> [a] -> a+minimum x xs = foldl min x xs++++{-# INLINE compareOrd #-}+compareOrd :: Ord a => a -> a -> a -> a -> a -> a+compareOrd a b lt eq gt =+   case P.compare a b of+      LT -> lt+      EQ -> eq+      GT -> gt++instance C Int where+   {-# INLINE compare #-}+   compare = compareOrd++instance C Integer where+   {-# INLINE compare #-}+   compare = compareOrd++instance C Float where+   {-# INLINE compare #-}+   compare = compareOrd++instance C Double where+   {-# INLINE compare #-}+   compare = compareOrd++instance C Bool where+   {-# INLINE compare #-}+   compare = compareOrd++instance C Ordering where+   {-# INLINE compare #-}+   compare = compareOrd++++{-# INLINE elementCompare #-}+elementCompare ::+   (C x) =>+   (v -> x) -> Elem.T (v,v,v,v,v) x+elementCompare f =+   Elem.element (\(x,y,lt,eq,gt) ->+      compare (f x) (f y) (f lt) (f eq) (f gt))++{-# INLINE (<*>.<=>?) #-}+(<*>.<=>?) ::+   (C x) =>+   Elem.T (v,v,v,v,v) (x -> a) -> (v -> x) -> Elem.T (v,v,v,v,v) a+(<*>.<=>?) f acc =+   f <*> elementCompare acc+++{-# INLINE element #-}+element ::+   (C x) =>+   (x -> x -> x -> x -> x) ->+   (v -> x) -> Elem.T (v,v,v,v) x+element cmp f =+   Elem.element (\(x,y,true,false) -> cmp (f x) (f y) (f true) (f false))++{-# INLINE (<*>.<=?) #-}+(<*>.<=?) ::+   (C x) =>+   Elem.T (v,v,v,v) (x -> a) -> (v -> x) -> Elem.T (v,v,v,v) a+(<*>.<=?) f acc =+   f <*> element (<=?) acc++{-# INLINE (<*>.>=?) #-}+(<*>.>=?) ::+   (C x) =>+   Elem.T (v,v,v,v) (x -> a) -> (v -> x) -> Elem.T (v,v,v,v) a+(<*>.>=?) f acc =+   f <*> element (>=?) acc++{-# INLINE (<*>.<?) #-}+(<*>.<?) ::+   (C x) =>+   Elem.T (v,v,v,v) (x -> a) -> (v -> x) -> Elem.T (v,v,v,v) a+(<*>.<?) f acc =+   f <*> element (<?) acc++{-# INLINE (<*>.>?) #-}+(<*>.>?) ::+   (C x) =>+   Elem.T (v,v,v,v) (x -> a) -> (v -> x) -> Elem.T (v,v,v,v) a+(<*>.>?) f acc =+   f <*> element (>?) acc+++instance (C a, C b) => C (a,b) where+   {-# INLINE compare #-}+   compare = Elem.run5 $ pure (,) <*>.<=>? fst <*>.<=>? snd+   {-# INLINE (<=?) #-}+   (<=?)   = Elem.run4 $ pure (,) <*>.<=?  fst <*>.<=?  snd+   {-# INLINE (>=?) #-}+   (>=?)   = Elem.run4 $ pure (,) <*>.>=?  fst <*>.>=?  snd+   {-# INLINE (<?) #-}+   (<?)    = Elem.run4 $ pure (,) <*>.<?   fst <*>.<?   snd+   {-# INLINE (>?) #-}+   (>?)    = Elem.run4 $ pure (,) <*>.>?   fst <*>.>?   snd++instance (C a, C b, C c) => C (a,b,c) where+   {-# INLINE compare #-}+   compare = Elem.run5 $ pure (,,) <*>.<=>? fst3 <*>.<=>? snd3 <*>.<=>? thd3+   {-# INLINE (<=?) #-}+   (<=?)   = Elem.run4 $ pure (,,) <*>.<=?  fst3 <*>.<=?  snd3 <*>.<=?  thd3+   {-# INLINE (>=?) #-}+   (>=?)   = Elem.run4 $ pure (,,) <*>.>=?  fst3 <*>.>=?  snd3 <*>.>=?  thd3+   {-# INLINE (<?) #-}+   (<?)    = Elem.run4 $ pure (,,) <*>.<?   fst3 <*>.<?   snd3 <*>.<?   thd3+   {-# INLINE (>?) #-}+   (>?)    = Elem.run4 $ pure (,,) <*>.>?   fst3 <*>.>?   snd3 <*>.>?   thd3++instance C a => C [a] where+   {-# INLINE compare #-}+   compare = zipWith5 compare+   {-# INLINE (<=?) #-}+   (<=?) = zipWith4 (<=?)+   {-# INLINE (>=?) #-}+   (>=?) = zipWith4 (>=?)+   {-# INLINE (<?) #-}+   (<?)  = zipWith4 (<?)+   {-# INLINE (>?) #-}+   (>?)  = zipWith4 (>?)++instance (C a) => C (b -> a) where+   {-# INLINE compare #-}+   compare x y lt eq gt c  =  compare (x c) (y c) (lt c) (eq c) (gt c)+   {-# INLINE (<=?) #-}+   (<=?) x y true false c  =  (x c <=? y c) (true c) (false c)+   {-# INLINE (>=?) #-}+   (>=?) x y true false c  =  (x c >=? y c) (true c) (false c)+   {-# INLINE (<?) #-}+   (<?)  x y true false c  =  (x c <?  y c) (true c) (false c)+   {-# INLINE (>?) #-}+   (>?)  x y true false c  =  (x c >?  y c) (true c) (false c)
src/Algebra/PrincipalIdealDomain.hs view
@@ -57,7 +57,7 @@  import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  ) -import PreludeBase+import NumericPrelude.Base import Prelude (Integer, Int) import qualified Prelude as P import Test.QuickCheck ((==>), Property)
− src/Algebra/Real.hs 
@@ -1,128 +0,0 @@-{-# LANGUAGE NoImplicitPrelude #-}-module Algebra.Real (-   C(abs, signum),-   ) where--import qualified Algebra.Ring         as Ring-import qualified Algebra.Additive     as Additive-import qualified Algebra.ZeroTestable as ZeroTestable--import Algebra.Ring (one, ) -- fromInteger-import Algebra.Additive (zero, negate,)--import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )-import Data.Word (Word, Word8, Word16, Word32, Word64, )--import PreludeBase-import qualified Prelude as P-import Prelude(Int,Integer,Float,Double)---{- |-This is the type class of an ordered ring, satisfying the laws-->              a * b === b * a->      a + (max b c) === max (a+b) (a+c)->   negate (max b c) === min (negate b) (negate c)->      a * (max b c) === max (a*b) (a*c) where a >= 0--Note that abs is in a rather different place than it is in the Haskell-98 Prelude.  In particular,-->   abs :: Complex -> Complex--is not defined.  To me, this seems to have the wrong type anyway;-Complex.magnitude has the correct type.--}-class (Ring.C a, ZeroTestable.C a, Ord a) => C a where-    abs    :: a -> a-    signum :: a -> a--      -- Minimal definition: nothing-    abs x    = max x (negate x)-    signum x = case compare x zero of-                 GT ->        one-                 EQ ->        zero-                 LT -> negate one---instance C Integer where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Float   where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Double  where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum---instance C Int     where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Int8    where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Int16   where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Int32   where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Int64   where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum---instance C Word    where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Word8   where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Word16  where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Word32  where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum--instance C Word64  where-   {-# INLINE abs #-}-   {-# INLINE signum #-}-   abs = P.abs-   signum = P.signum-
src/Algebra/RealField.hs view
@@ -1,459 +1,26 @@ {-# LANGUAGE NoImplicitPrelude #-}-{-# OPTIONS_GHC -fglasgow-exts #-}--- -fglasgow-exts for RULES-module Algebra.RealField where+module Algebra.RealField (+   C,+   ) where -import qualified Algebra.Field              as Field+import qualified Algebra.Field as Field+import qualified Algebra.RealRing as RealRing import qualified Algebra.PrincipalIdealDomain as PID-import qualified Algebra.Real           as Real-import qualified Algebra.Ring           as Ring-import qualified Algebra.ToRational     as ToRational-import qualified Algebra.ToInteger      as ToInteger---import Algebra.Field          ((/), fromRational, )-import Algebra.RealIntegral   (quotRem, )-import Algebra.IntegralDomain (divMod, even, )-import Algebra.Ring           ((*), fromInteger, one, )-import Algebra.Additive       ((+), (-), negate, zero, )-import Algebra.ZeroTestable   (isZero, )-import Algebra.ToInteger      (fromIntegral, )+import qualified Algebra.ToInteger as ToInteger  import qualified Number.Ratio as Ratio-import Number.Ratio (T((:%)), Rational) -import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )-import Data.Word (Word, Word8, Word16, Word32, Word64, )--import qualified GHC.Float as GHC-import Data.List as List-import Data.Tuple.HT (mapFst, mapPair, )-import Prelude(Int, Integer, Float, Double)+-- import NumericPrelude.Base import qualified Prelude as P-import PreludeBase-+import Prelude (Float, Double, )  {- |-Minimal complete definition:-     'splitFraction' or 'floor'--There are probably more laws, but some laws are--> (fromInteger.fst.splitFraction) a + (snd.splitFraction) a === a->    ceiling (toRational x) === ceiling x :: Integer->   truncate (toRational x) === truncate x :: Integer->      floor (toRational x) === floor x :: Integer--If there wouldn't be @Real.C a@ and @ToInteger.C b@ constraints,-we could also use this class for splitting ratios of polynomials.--As an aside, let me note the similarities-between @splitFraction x@ and @x divMod 1@ (if that were defined).-In particular, it might make sense to unify the rounding modes somehow.--IEEEFloat-specific calls are removed here (cf. 'Prelude.RealFloat')-so probably nobody will actually use this default definition.--Henning:-New function 'fraction' doesn't return the integer part of the number.-This also removes a type ambiguity if the integer part is not needed.--The new methods 'fraction' and 'splitFraction'-differ from 'Prelude.properFraction' semantics.-They always round to 'floor'.-This means that the fraction is always non-negative and-is always smaller than 1.-This is more useful in practice and-can be generalised to more than real numbers.-Since every 'Number.Ratio.T' denominator type supports 'Algebra.IntegralDomain.divMod',-every 'Number.Ratio.T' can provide 'fraction' and 'splitFraction',-e.g. fractions of polynomials.-However the ''integral'' part would not be of type class 'ToInteger.C'.--Can there be a separate class for-'fraction', 'splitFraction', 'floor' and 'ceiling'-since they do not need reals and their ordering?--Note:-All of these methods can be defined-exclusively with functions from Ord and Ring.-We could write a power-of-two-algorithm-like the one for finding the number of digits of an Integer-in FixedPoint-fractions module.-This would even be reasonably efficient.-I think the module should be renamed to RealRing,-and the superclass constraint should be lifted from Field to Ring.--We might also add a round method,-that rounds 0.5 always up or always down.-This is much more efficient in inner loops-and is acceptable or even preferable for many applications.--The ToInteger constraint can be lifted to Ring.+This is a convenient class for common types+that both form a field and have a notion of ordering by magnitude. -}--class (Real.C a, Field.C a) => C a where-    splitFraction    :: (ToInteger.C b) => a -> (b,a)-    fraction         ::                    a -> a-    ceiling, floor   :: (ToInteger.C b) => a -> b-    truncate, round  :: (ToInteger.C b) => a -> b---    splitFraction x   =  (floor x, fraction x)--    fraction x   =  x - fromInteger (floor x)--    floor x      =  fromInteger (fst (splitFraction x))--    ceiling x    =  - floor (-x)----    truncate x   =  signum x * floor (abs x)-    truncate x   =  if x>=0-                      then floor x-                      else ceiling x--    round x      =  let (n,r) = splitFraction x-                    in  case compare r (1/2) of-                           LT -> n-                           EQ -> if even n then n else n+1-                           GT -> n+1---instance (ToInteger.C a, PID.C a) => C (Ratio.T a) where-    splitFraction (x:%y) = (fromIntegral q, r:%y)-                               where (q,r) = divMod x y+class (RealRing.C a, Field.C a) => C a where  instance C Float where-    {-# INLINE splitFraction #-}-    {-# INLINE fraction #-}-    {-# INLINE floor #-}-    {-# INLINE ceiling #-}-    {-# INLINE round #-}-    {-# INLINE truncate #-}-    splitFraction = fastSplitFraction GHC.float2Int GHC.int2Float-    fraction      = fastFraction (GHC.int2Float . GHC.float2Int)-    floor         = fromInteger . P.floor-    ceiling       = fromInteger . P.ceiling-    round         = fromInteger . P.round-    truncate      = fromInteger . P.truncate- instance C Double where-    {-# INLINE splitFraction #-}-    {-# INLINE fraction #-}-    {-# INLINE floor #-}-    {-# INLINE ceiling #-}-    {-# INLINE round #-}-    {-# INLINE truncate #-}-    splitFraction = fastSplitFraction GHC.double2Int GHC.int2Double-    fraction      = fastFraction (GHC.int2Double . GHC.double2Int)-    floor         = fromInteger . P.floor-    ceiling       = fromInteger . P.ceiling-    round         = fromInteger . P.round-    truncate      = fromInteger . P.truncate --{-# INLINE fastSplitFraction #-}-fastSplitFraction :: (P.RealFrac a, Real.C a, ToInteger.C b) =>-   (a -> Int) -> (Int -> a) -> a -> (b,a)-fastSplitFraction trunc toFloat x =-   fixSplitFraction $-   if toFloat minBound <= x && x <= toFloat maxBound-     then case trunc x of n -> (fromIntegral n, x - toFloat n)-     else case P.properFraction x of (n,f) -> (fromInteger n, f)--{-# INLINE fixSplitFraction #-}-fixSplitFraction :: (Ring.C a, Ring.C b, Ord a) => (b,a) -> (b,a)-fixSplitFraction (n,f) =-   --  if x>=0 || f==0-   if f>=0-     then (n,   f)-     else (n-1, f+1)--{-# INLINE fastFraction #-}-fastFraction :: (P.RealFrac a, Real.C a) => (a -> a) -> a -> a-fastFraction trunc x =-   fixFraction $-   if fromIntegral (minBound :: Int) <= x && x <= fromIntegral (maxBound :: Int)-     then x - trunc x-     else preludeFraction x--{-# INLINE preludeFraction #-}-preludeFraction :: (P.RealFrac a, Ring.C a) => a -> a-preludeFraction x =-   let second :: (Integer, a) -> a-       second = snd-   in  second (P.properFraction x)--{-# INLINE fixFraction #-}-fixFraction :: (Ring.C a, Ord a) => a -> a-fixFraction y =-   if y>=0 then y else y+1--{--mapM_ (\n -> let x = fromInteger n / 10 in print (x, floorInt GHC.double2Int GHC.int2Double x)) [-20,-19..20]--}--{-# INLINE splitFractionInt #-}-splitFractionInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> (Int, a)-splitFractionInt trunc toFloat x =-   let n = trunc x-   in  fixSplitFraction (n, x - toFloat n)--{-# INLINE floorInt #-}-floorInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int-floorInt trunc toFloat x =-   let n = trunc x-   in  if x >= toFloat n-         then n-         else pred n--{-# INLINE ceilingInt #-}-ceilingInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int-ceilingInt trunc toFloat x =-   let n = trunc x-   in  if x <= toFloat n-         then n-         else succ n--{-# INLINE roundInt #-}-roundInt :: (Field.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int-roundInt trunc toFloat x =-   let half = 0.5 -- P.fromRational-       halfUp = x+half-       n = floorInt trunc toFloat halfUp-   in  if toFloat n == halfUp  &&  P.odd n-         then pred n-         else n---{- RULES maybe used, when Prelude implementations become more efficient-     "NP.round    :: Float -> Int"    round    = P.round    :: Float -> Int;-     "NP.truncate :: Float -> Int"    truncate = P.truncate :: Float -> Int;-     "NP.floor    :: Float -> Int"    floor    = P.floor    :: Float -> Int;-     "NP.ceiling  :: Float -> Int"    ceiling  = P.ceiling  :: Float -> Int;-     "NP.round    :: Double -> Int"   round    = P.round    :: Double -> Int;-     "NP.truncate :: Double -> Int"   truncate = P.truncate :: Double -> Int;-     "NP.floor    :: Double -> Int"   floor    = P.floor    :: Double -> Int;-     "NP.ceiling  :: Double -> Int"   ceiling  = P.ceiling  :: Double -> Int;-  -}---- these rules will also be needed for Int16 et.al.-{-# RULES-     "NP.round    :: Float -> Int"    round    = roundInt    GHC.float2Int  GHC.int2Float;-     "NP.truncate :: Float -> Int"    truncate =             GHC.float2Int               ;-     "NP.floor    :: Float -> Int"    floor    = floorInt    GHC.float2Int  GHC.int2Float;-     "NP.ceiling  :: Float -> Int"    ceiling  = ceilingInt  GHC.float2Int  GHC.int2Float;-     "NP.round    :: Double -> Int"   round    = roundInt    GHC.double2Int GHC.int2Double;-     "NP.truncate :: Double -> Int"   truncate =             GHC.double2Int               ;-     "NP.floor    :: Double -> Int"   floor    = floorInt    GHC.double2Int GHC.int2Double;-     "NP.ceiling  :: Double -> Int"   ceiling  = ceilingInt  GHC.double2Int GHC.int2Double;--     "NP.splitFraction :: Float ->  (Int, Float)"  splitFraction = splitFractionInt GHC.float2Int GHC.int2Float;-     "NP.splitFraction :: Double -> (Int, Double)" splitFraction = splitFractionInt GHC.double2Int GHC.int2Double;-  #-}---- generated by GenerateRules.hs-{-# RULES-     "NP.round    :: a -> Int8"    round    = (P.fromIntegral :: Int -> Int8) . round;-     "NP.truncate :: a -> Int8"    truncate = (P.fromIntegral :: Int -> Int8) . truncate;-     "NP.floor    :: a -> Int8"    floor    = (P.fromIntegral :: Int -> Int8) . floor;-     "NP.ceiling  :: a -> Int8"    ceiling  = (P.fromIntegral :: Int -> Int8) . ceiling;-     "NP.round    :: a -> Int16"   round    = (P.fromIntegral :: Int -> Int16) . round;-     "NP.truncate :: a -> Int16"   truncate = (P.fromIntegral :: Int -> Int16) . truncate;-     "NP.floor    :: a -> Int16"   floor    = (P.fromIntegral :: Int -> Int16) . floor;-     "NP.ceiling  :: a -> Int16"   ceiling  = (P.fromIntegral :: Int -> Int16) . ceiling;-     "NP.round    :: a -> Int32"   round    = (P.fromIntegral :: Int -> Int32) . round;-     "NP.truncate :: a -> Int32"   truncate = (P.fromIntegral :: Int -> Int32) . truncate;-     "NP.floor    :: a -> Int32"   floor    = (P.fromIntegral :: Int -> Int32) . floor;-     "NP.ceiling  :: a -> Int32"   ceiling  = (P.fromIntegral :: Int -> Int32) . ceiling;-     "NP.round    :: a -> Int64"   round    = (P.fromIntegral :: Int -> Int64) . round;-     "NP.truncate :: a -> Int64"   truncate = (P.fromIntegral :: Int -> Int64) . truncate;-     "NP.floor    :: a -> Int64"   floor    = (P.fromIntegral :: Int -> Int64) . floor;-     "NP.ceiling  :: a -> Int64"   ceiling  = (P.fromIntegral :: Int -> Int64) . ceiling;-     "NP.round    :: a -> Word"    round    = (P.fromIntegral :: Int -> Word) . round;-     "NP.truncate :: a -> Word"    truncate = (P.fromIntegral :: Int -> Word) . truncate;-     "NP.floor    :: a -> Word"    floor    = (P.fromIntegral :: Int -> Word) . floor;-     "NP.ceiling  :: a -> Word"    ceiling  = (P.fromIntegral :: Int -> Word) . ceiling;-     "NP.round    :: a -> Word8"   round    = (P.fromIntegral :: Int -> Word8) . round;-     "NP.truncate :: a -> Word8"   truncate = (P.fromIntegral :: Int -> Word8) . truncate;-     "NP.floor    :: a -> Word8"   floor    = (P.fromIntegral :: Int -> Word8) . floor;-     "NP.ceiling  :: a -> Word8"   ceiling  = (P.fromIntegral :: Int -> Word8) . ceiling;-     "NP.round    :: a -> Word16"  round    = (P.fromIntegral :: Int -> Word16) . round;-     "NP.truncate :: a -> Word16"  truncate = (P.fromIntegral :: Int -> Word16) . truncate;-     "NP.floor    :: a -> Word16"  floor    = (P.fromIntegral :: Int -> Word16) . floor;-     "NP.ceiling  :: a -> Word16"  ceiling  = (P.fromIntegral :: Int -> Word16) . ceiling;-     "NP.round    :: a -> Word32"  round    = (P.fromIntegral :: Int -> Word32) . round;-     "NP.truncate :: a -> Word32"  truncate = (P.fromIntegral :: Int -> Word32) . truncate;-     "NP.floor    :: a -> Word32"  floor    = (P.fromIntegral :: Int -> Word32) . floor;-     "NP.ceiling  :: a -> Word32"  ceiling  = (P.fromIntegral :: Int -> Word32) . ceiling;-     "NP.round    :: a -> Word64"  round    = (P.fromIntegral :: Int -> Word64) . round;-     "NP.truncate :: a -> Word64"  truncate = (P.fromIntegral :: Int -> Word64) . truncate;-     "NP.floor    :: a -> Word64"  floor    = (P.fromIntegral :: Int -> Word64) . floor;-     "NP.ceiling  :: a -> Word64"  ceiling  = (P.fromIntegral :: Int -> Word64) . ceiling;--     "NP.splitFraction :: a -> (Int8,a)"     splitFraction = mapFst (P.fromIntegral :: Int -> Int8) . splitFraction;-     "NP.splitFraction :: a -> (Int16,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Int16) . splitFraction;-     "NP.splitFraction :: a -> (Int32,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Int32) . splitFraction;-     "NP.splitFraction :: a -> (Int64,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Int64) . splitFraction;-     "NP.splitFraction :: a -> (Word,a)"     splitFraction = mapFst (P.fromIntegral :: Int -> Word) . splitFraction;-     "NP.splitFraction :: a -> (Word8,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Word8) . splitFraction;-     "NP.splitFraction :: a -> (Word16,a)"   splitFraction = mapFst (P.fromIntegral :: Int -> Word16) . splitFraction;-     "NP.splitFraction :: a -> (Word32,a)"   splitFraction = mapFst (P.fromIntegral :: Int -> Word32) . splitFraction;-     "NP.splitFraction :: a -> (Word64,a)"   splitFraction = mapFst (P.fromIntegral :: Int -> Word64) . splitFraction;-  #-}---{- | TODO: Should be moved to a continued fraction module. -}--approxRational :: (ToRational.C a, C a) => a -> a -> Rational-approxRational rat eps    =  simplest (rat-eps) (rat+eps)-        where simplest x y | y < x      =  simplest y x-                           | x == y     =  xr-                           | x > 0      =  simplest' n d n' d'-                           | y < 0      =  - simplest' (-n') d' (-n) d-                           | otherwise  =  0 :% 1-                                        where xr@(n:%d) = ToRational.toRational x-                                              (n':%d')  = ToRational.toRational y--              simplest' n d n' d'       -- assumes 0 < n%d < n'%d'-                        | isZero r   =  q :% 1-                        | q /= q'    =  (q+1) :% 1-                        | otherwise  =  (q*n''+d'') :% n''-                                     where (q,r)      =  quotRem n d-                                           (q',r')    =  quotRem n' d'-                                           (n'':%d'') =  simplest' d' r' d r----- * generic implementation of round functions--powersOfTwo :: (Ring.C a) => [a]-powersOfTwo = iterate (2*) one--pairsOfPowersOfTwo :: (Ring.C a, Ring.C b) => [(a,b)]-pairsOfPowersOfTwo =-   zip powersOfTwo powersOfTwo--{- |-The generic rounding functions need a number of operations-proportional to the number of binary digits of the integer portion.-If operations like multiplication with two and comparison-need time proportional to the number of binary digits,-then the overall rounding requires quadratic time.--}-genericFloor :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericFloor a =-   if a>=zero-     then genericPosFloor a-     else negate $ genericPosCeiling $ negate a--genericCeiling :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericCeiling a =-   if a>=zero-     then genericPosCeiling a-     else negate $ genericPosFloor $ negate a--genericTruncate :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericTruncate a =-   if a>=zero-     then genericPosFloor a-     else negate $ genericPosFloor $ negate a--genericRound :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericRound a =-   if a>=zero-     then genericPosRound a-     else negate $ genericPosRound $ negate a--genericFraction :: (Ord a, Ring.C a) => a -> a-genericFraction a =-   if a>=zero-     then genericPosFraction a-     else fixFraction $ negate $ genericPosFraction $ negate a--genericSplitFraction :: (Ord a, Ring.C a, Ring.C b) => a -> (b,a)-genericSplitFraction a =-   if a>=zero-     then genericPosSplitFraction a-     else fixSplitFraction $ mapPair (negate, negate) $-          genericPosSplitFraction $ negate a---genericPosFloor :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericPosFloor a =-   snd $-   foldr-      (\(pa,pb) acc@(accA,accB) ->-         let newA = accA+pa-         in  if newA>a then acc else (newA,accB+pb))-      (zero,zero) $-   takeWhile ((a>=) . fst) $-   pairsOfPowersOfTwo--genericPosCeiling :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericPosCeiling a =-   snd $-   (\(ps,u:_) ->-      foldr-         (\(pa,pb) acc@(accA,accB) ->-            let newA = accA-pa-            in  if newA>=a then (newA,accB-pb) else acc)-         u ps) $-   span ((a>) . fst) $-   (zero,zero) : pairsOfPowersOfTwo--{--genericPosFloorDigits :: (Ord a, Ring.C a, Ring.C b) => a -> ((a,b), [Bool])-genericPosFloorDigits a =-   List.mapAccumR-      (\acc@(accA,accB) (pa,pb) ->-         let newA = accA+pa-             b = newA<=a-         in  (if b then (newA,accB+pb) else acc, b))-      (zero,zero) $-   takeWhile ((a>=) . fst) $-   pairsOfPowersOfTwo--}--genericHalfPosFloorDigits :: (Ord a, Ring.C a, Ring.C b) => a -> ((a,b), [Bool])-genericHalfPosFloorDigits a =-   List.mapAccumR-      (\acc@(accA,accB) (pa,pb) ->-         let newA = accA+pa-             b = newA<=a-         in  (if b then (newA,accB+pb) else acc, b))-      (zero,zero) $-   takeWhile ((a>=) . fst) $-   zip powersOfTwo (zero:powersOfTwo)--genericPosRound :: (Ord a, Ring.C a, Ring.C b) => a -> b-genericPosRound a =-   let a2 = 2*a-       ((ai,bi), ds) = genericHalfPosFloorDigits a2-   in  if ai==a2-         then-           case ds of-             True : True : _ -> bi+one-             _ -> bi-         else-           case ds of-             True : _ -> bi+one-             _ -> bi--genericPosFraction :: (Ord a, Ring.C a) => a -> a-genericPosFraction a =-   foldr-      (\p acc ->-         if p>acc then acc else acc-p)-      a $-   takeWhile (a>=) $-   powersOfTwo--genericPosSplitFraction :: (Ord a, Ring.C a, Ring.C b) => a -> (b,a)-genericPosSplitFraction a =-   foldr-      (\(pb,pa) acc@(accB,accA) ->-         if pa>accA then acc else (accB+pb,accA-pa))-      (zero,a) $-   takeWhile ((a>=) . snd) $-   pairsOfPowersOfTwo-+instance (ToInteger.C a, PID.C a) => C (Ratio.T a) where
src/Algebra/RealIntegral.hs view
@@ -17,11 +17,11 @@    ) where  import qualified Algebra.IntegralDomain as Integral-import qualified Algebra.Real           as Real+import qualified Algebra.Absolute           as Absolute import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive -import Algebra.Real (signum, )+import Algebra.Absolute (signum, ) import Algebra.IntegralDomain (divMod, ) import Algebra.Ring (one, ) -- fromInteger import Algebra.Additive (zero, (+), (-), )@@ -29,7 +29,7 @@ import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  ) import Data.Word (Word, Word8, Word16, Word32, Word64, ) -import PreludeBase+import NumericPrelude.Base import qualified Prelude as P import Prelude (Int, Integer, ) @@ -46,7 +46,7 @@ Minimal definition: nothing required -} -class (Real.C a, Integral.C a) => C a where+class (Absolute.C a, Ord a, Integral.C a) => C a where     quot, rem        :: a -> a -> a     quotRem          :: a -> a -> (a,a)
+ src/Algebra/RealRing.hs view
@@ -0,0 +1,586 @@+{-# LANGUAGE NoImplicitPrelude #-}+{-# OPTIONS_GHC -fglasgow-exts #-}+-- -fglasgow-exts for RULES+module Algebra.RealRing where++import qualified Algebra.Field              as Field+import qualified Algebra.PrincipalIdealDomain as PID+import qualified Algebra.Absolute           as Absolute+import qualified Algebra.Ring           as Ring+import qualified Algebra.ToRational     as ToRational+import qualified Algebra.ToInteger      as ToInteger++import qualified Algebra.OrderDecision as OrdDec+import Algebra.OrderDecision ((<?), (>=?), )++import Algebra.Field          (fromRational, )+import Algebra.RealIntegral   (quotRem, )+import Algebra.IntegralDomain (divMod, even, )+import Algebra.Ring           ((*), fromInteger, one, )+import Algebra.Additive       ((+), (-), negate, zero, )+import Algebra.ZeroTestable   (isZero, )+import Algebra.ToInteger      (fromIntegral, )++import qualified Number.Ratio as Ratio+import Number.Ratio (T((:%)), Rational)++import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )+import Data.Word (Word, Word8, Word16, Word32, Word64, )++import qualified GHC.Float as GHC+import Data.List as List+import Data.Tuple.HT (mapFst, mapPair, )+import Prelude(Int, Integer, Float, Double)+import qualified Prelude as P+import NumericPrelude.Base+++{- |+Minimal complete definition:+     'splitFraction' or 'floor'++There are probably more laws, but some laws are++> splitFraction x === (fromInteger (floor x), fraction x)+> fromInteger (floor x) + fraction x === x+> floor x       <= x       x <  floor x + 1+> ceiling x - 1 <  x       x <= ceiling x+> 0 <= fraction x          fraction x < 1++>               - ceiling x === floor (-x)+>                truncate x === signum x * floor (abs x)+>    ceiling (toRational x) === ceiling x :: Integer+>   truncate (toRational x) === truncate x :: Integer+>      floor (toRational x) === floor x :: Integer++The new function 'fraction' doesn't return the integer part of the number.+This also removes a type ambiguity if the integer part is not needed.++Many people will associate rounding with fractional numbers,+and thus they are surprised about the superclass being @Ring@ not @Field@.+The reason is that all of these methods can be defined+exclusively with functions from @Ord@ and @Ring@.+The implementations of 'genericFloor' and other functions demonstrate that.+They implement power-of-two-algorithms+like the one for finding the number of digits of an 'Integer'+in FixedPoint-fractions module.+They are even reasonably efficient.++I am still uncertain whether it was a good idea+to add instances for @Integer@ and friends,+since calling @floor@ or @fraction@ on an integer may well indicate a bug.+The rounding functions are just the identity function+and 'fraction' is constant zero.+However, I decided to associate our class with @Ring@ rather than @Field@,+after I found myself using repeated subtraction and testing+rather than just calling @fraction@,+just in order to get the constraint @(Ring a, Ord a)@+that was more general than @(RealField a)@.++For the results of the rounding functions+we have chosen the constraint @Ring@ instead of @ToInteger@,+since this is more flexible to use,+but it still signals to the user that only integral numbers can be returned.+This is so, because the plain @Ring@ class only provides+@zero@, @one@ and operations that allow to reach all natural numbers but not more.+++As an aside, let me note the similarities+between @splitFraction x@ and @divMod x 1@ (if that were defined).+In particular, it might make sense to unify the rounding modes somehow.++The new methods 'fraction' and 'splitFraction'+differ from 'Prelude.properFraction' semantics.+They always round to 'floor'.+This means that the fraction is always non-negative and+is always smaller than 1.+This is more useful in practice and+can be generalised to more than real numbers.+Since every 'Number.Ratio.T' denominator type+supports 'Algebra.IntegralDomain.divMod',+every 'Number.Ratio.T' can provide 'fraction' and 'splitFraction',+e.g. fractions of polynomials.+However the @Ring@ constraint for the ''integral'' part of 'splitFraction'+is too weak in order to generate polynomials.+After all, I am uncertain whether this would be useful or not.++Can there be a separate class for+'fraction', 'splitFraction', 'floor' and 'ceiling'+since they do not need reals and their ordering?++We might also add a round method,+that rounds 0.5 always up or always down.+This is much more efficient in inner loops+and is acceptable or even preferable for many applications.+-}++class (Absolute.C a, Ord a) => C a where+    splitFraction    :: (Ring.C b) => a -> (b,a)+    fraction         ::               a -> a+    ceiling, floor   :: (Ring.C b) => a -> b+    truncate         :: (Ring.C b) => a -> b+    round            :: (ToInteger.C b) => a -> b+++    splitFraction x   =  (floor x, fraction x)++    fraction x   =  x - fromInteger (floor x)++    floor x      =  fromInteger (fst (splitFraction x))++    ceiling x    =  - floor (-x)++--    truncate x   =  signum x * floor (abs x)+    truncate x =+       if x>=0+         then floor x+         else ceiling x++    {-+    The ToInteger constraint can be lifted to Ring+    if use Integer temporarily.+    I expect this would not be efficient in many cases.+    -}+    round x =+       let (n,r) = splitFraction x+       in  case compare (2*r) one of+              LT -> n+              EQ -> if even n then n else n+1+              GT -> n+1+++{- |+This function rounds to the closest integer.+For @fraction x == 0.5@ it rounds away from zero.+This function is not the result of an ingenious mathematical insight,+but is simply a kind of rounding that is the fastest+on IEEE floating point architectures.+-}+roundSimple :: (C a, Ring.C b) => a -> b+roundSimple x =+   let (n,r) = splitFraction x+   in  case compare (2*r) one of+          LT -> n+          EQ -> if x<0 then n else n+1+          GT -> n+1+++instance (ToInteger.C a, PID.C a) => C (Ratio.T a) where+    splitFraction (x:%y) = (fromIntegral q, r:%y)+                               where (q,r) = divMod x y++instance C Int where+    {-# INLINE splitFraction #-}+    {-# INLINE fraction #-}+    {-# INLINE floor #-}+    {-# INLINE ceiling #-}+    {-# INLINE round #-}+    {-# INLINE truncate #-}+    splitFraction x = (fromIntegral x, zero)+    fraction      _ = zero+    floor         x = fromIntegral x+    ceiling       x = fromIntegral x+    round         x = fromIntegral x+    truncate      x = fromIntegral x++instance C Integer where+    {-# INLINE splitFraction #-}+    {-# INLINE fraction #-}+    {-# INLINE floor #-}+    {-# INLINE ceiling #-}+    {-# INLINE round #-}+    {-# INLINE truncate #-}+    splitFraction x = (fromInteger x, zero)+    fraction      _ = zero+    floor         x = fromInteger x+    ceiling       x = fromInteger x+    round         x = fromInteger x+    truncate      x = fromInteger x++instance C Float where+    {-# INLINE splitFraction #-}+    {-# INLINE fraction #-}+    {-# INLINE floor #-}+    {-# INLINE ceiling #-}+    {-# INLINE round #-}+    {-# INLINE truncate #-}+    splitFraction = fastSplitFraction GHC.float2Int GHC.int2Float+    fraction      = fastFraction (GHC.int2Float . GHC.float2Int)+    floor         = fromInteger . P.floor+    ceiling       = fromInteger . P.ceiling+    round         = fromInteger . P.round+    truncate      = fromInteger . P.truncate++instance C Double where+    {-# INLINE splitFraction #-}+    {-# INLINE fraction #-}+    {-# INLINE floor #-}+    {-# INLINE ceiling #-}+    {-# INLINE round #-}+    {-# INLINE truncate #-}+    splitFraction = fastSplitFraction GHC.double2Int GHC.int2Double+    fraction      = fastFraction (GHC.int2Double . GHC.double2Int)+    floor         = fromInteger . P.floor+    ceiling       = fromInteger . P.ceiling+    round         = fromInteger . P.round+    truncate      = fromInteger . P.truncate+++{-# INLINE fastSplitFraction #-}+fastSplitFraction :: (P.RealFrac a, Absolute.C a, Ring.C b) =>+   (a -> Int) -> (Int -> a) -> a -> (b,a)+fastSplitFraction trunc toFloat x =+   fixSplitFraction $+   if toFloat minBound <= x && x <= toFloat maxBound+     then case trunc x of n -> (fromIntegral n, x - toFloat n)+     else case P.properFraction x of (n,f) -> (fromInteger n, f)++{-# INLINE fixSplitFraction #-}+fixSplitFraction :: (Ring.C a, Ring.C b, Ord a) => (b,a) -> (b,a)+fixSplitFraction (n,f) =+   --  if x>=0 || f==0+   if f>=0+     then (n,   f)+     else (n-1, f+1)++{-# INLINE fastFraction #-}+fastFraction :: (P.RealFrac a, Absolute.C a) => (a -> a) -> a -> a+fastFraction trunc x =+   fixFraction $+   if fromIntegral (minBound :: Int) <= x && x <= fromIntegral (maxBound :: Int)+     then x - trunc x+     else preludeFraction x++{-# INLINE preludeFraction #-}+preludeFraction :: (P.RealFrac a, Ring.C a) => a -> a+preludeFraction x =+   let second :: (Integer, a) -> a+       second = snd+   in  second (P.properFraction x)++{-# INLINE fixFraction #-}+fixFraction :: (Ring.C a, Ord a) => a -> a+fixFraction y =+   if y>=0 then y else y+1++{-+mapM_ (\n -> let x = fromInteger n / 10 in print (x, floorInt GHC.double2Int GHC.int2Double x)) [-20,-19..20]+-}++{-# INLINE splitFractionInt #-}+splitFractionInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> (Int, a)+splitFractionInt trunc toFloat x =+   let n = trunc x+   in  fixSplitFraction (n, x - toFloat n)++{-# INLINE floorInt #-}+floorInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int+floorInt trunc toFloat x =+   let n = trunc x+   in  if x >= toFloat n+         then n+         else pred n++{-# INLINE ceilingInt #-}+ceilingInt :: (Ring.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int+ceilingInt trunc toFloat x =+   let n = trunc x+   in  if x <= toFloat n+         then n+         else succ n++{-# INLINE roundInt #-}+roundInt :: (Field.C a, Ord a) => (a -> Int) -> (Int -> a) -> a -> Int+roundInt trunc toFloat x =+   let half = 0.5 -- P.fromRational+       halfUp = x+half+       n = floorInt trunc toFloat halfUp+   in  if toFloat n == halfUp  &&  P.odd n+         then pred n+         else n++{-# INLINE roundSimpleInt #-}+roundSimpleInt ::+   (Field.C a, Absolute.C a, Ord a) =>+   (a -> Int) -> (Int -> a) -> a -> Int+roundSimpleInt trunc _toFloat x =+   trunc (x + Absolute.signum x * 0.5)++++{- RULES maybe used, when Prelude implementations become more efficient+     "NP.round    :: Float -> Int"    round    = P.round    :: Float -> Int;+     "NP.truncate :: Float -> Int"    truncate = P.truncate :: Float -> Int;+     "NP.floor    :: Float -> Int"    floor    = P.floor    :: Float -> Int;+     "NP.ceiling  :: Float -> Int"    ceiling  = P.ceiling  :: Float -> Int;+     "NP.round    :: Double -> Int"   round    = P.round    :: Double -> Int;+     "NP.truncate :: Double -> Int"   truncate = P.truncate :: Double -> Int;+     "NP.floor    :: Double -> Int"   floor    = P.floor    :: Double -> Int;+     "NP.ceiling  :: Double -> Int"   ceiling  = P.ceiling  :: Double -> Int;+  -}++-- these rules will also be needed for Int16 et.al.+{-# RULES+     "NP.round       :: Float -> Int"    round    = roundInt       GHC.float2Int  GHC.int2Float;+     "NP.roundSimple :: Float -> Int"    round    = roundSimpleInt GHC.float2Int  GHC.int2Float;+     "NP.truncate    :: Float -> Int"    truncate =                GHC.float2Int               ;+     "NP.floor       :: Float -> Int"    floor    = floorInt       GHC.float2Int  GHC.int2Float;+     "NP.ceiling     :: Float -> Int"    ceiling  = ceilingInt     GHC.float2Int  GHC.int2Float;+     "NP.round       :: Double -> Int"   round    = roundInt       GHC.double2Int GHC.int2Double;+     "NP.roundSimple :: Double -> Int"   round    = roundSimpleInt GHC.double2Int GHC.int2Double;+     "NP.truncate    :: Double -> Int"   truncate =                GHC.double2Int               ;+     "NP.floor       :: Double -> Int"   floor    = floorInt       GHC.double2Int GHC.int2Double;+     "NP.ceiling     :: Double -> Int"   ceiling  = ceilingInt     GHC.double2Int GHC.int2Double;++     "NP.splitFraction :: Float ->  (Int, Float)"  splitFraction = splitFractionInt GHC.float2Int GHC.int2Float;+     "NP.splitFraction :: Double -> (Int, Double)" splitFraction = splitFractionInt GHC.double2Int GHC.int2Double;+  #-}++-- generated by GenerateRules.hs+{-# RULES+     "NP.round       :: a -> Int8"    round       = (P.fromIntegral :: Int -> Int8) . round;+     "NP.roundSimple :: a -> Int8"    roundSimple = (P.fromIntegral :: Int -> Int8) . roundSimple;+     "NP.truncate    :: a -> Int8"    truncate    = (P.fromIntegral :: Int -> Int8) . truncate;+     "NP.floor       :: a -> Int8"    floor       = (P.fromIntegral :: Int -> Int8) . floor;+     "NP.ceiling     :: a -> Int8"    ceiling     = (P.fromIntegral :: Int -> Int8) . ceiling;+     "NP.round       :: a -> Int16"   round       = (P.fromIntegral :: Int -> Int16) . round;+     "NP.roundSimple :: a -> Int16"   roundSimple = (P.fromIntegral :: Int -> Int16) . roundSimple;+     "NP.truncate    :: a -> Int16"   truncate    = (P.fromIntegral :: Int -> Int16) . truncate;+     "NP.floor       :: a -> Int16"   floor       = (P.fromIntegral :: Int -> Int16) . floor;+     "NP.ceiling     :: a -> Int16"   ceiling     = (P.fromIntegral :: Int -> Int16) . ceiling;+     "NP.round       :: a -> Int32"   round       = (P.fromIntegral :: Int -> Int32) . round;+     "NP.roundSimple :: a -> Int32"   roundSimple = (P.fromIntegral :: Int -> Int32) . roundSimple;+     "NP.truncate    :: a -> Int32"   truncate    = (P.fromIntegral :: Int -> Int32) . truncate;+     "NP.floor       :: a -> Int32"   floor       = (P.fromIntegral :: Int -> Int32) . floor;+     "NP.ceiling     :: a -> Int32"   ceiling     = (P.fromIntegral :: Int -> Int32) . ceiling;+     "NP.round       :: a -> Int64"   round       = (P.fromIntegral :: Int -> Int64) . round;+     "NP.roundSimple :: a -> Int64"   roundSimple = (P.fromIntegral :: Int -> Int64) . roundSimple;+     "NP.truncate    :: a -> Int64"   truncate    = (P.fromIntegral :: Int -> Int64) . truncate;+     "NP.floor       :: a -> Int64"   floor       = (P.fromIntegral :: Int -> Int64) . floor;+     "NP.ceiling     :: a -> Int64"   ceiling     = (P.fromIntegral :: Int -> Int64) . ceiling;+     "NP.round       :: a -> Word"    round       = (P.fromIntegral :: Int -> Word) . round;+     "NP.roundSimple :: a -> Word"    roundSimple = (P.fromIntegral :: Int -> Word) . roundSimple;+     "NP.truncate    :: a -> Word"    truncate    = (P.fromIntegral :: Int -> Word) . truncate;+     "NP.floor       :: a -> Word"    floor       = (P.fromIntegral :: Int -> Word) . floor;+     "NP.ceiling     :: a -> Word"    ceiling     = (P.fromIntegral :: Int -> Word) . ceiling;+     "NP.round       :: a -> Word8"   round       = (P.fromIntegral :: Int -> Word8) . round;+     "NP.roundSimple :: a -> Word8"   roundSimple = (P.fromIntegral :: Int -> Word8) . roundSimple;+     "NP.truncate    :: a -> Word8"   truncate    = (P.fromIntegral :: Int -> Word8) . truncate;+     "NP.floor       :: a -> Word8"   floor       = (P.fromIntegral :: Int -> Word8) . floor;+     "NP.ceiling     :: a -> Word8"   ceiling     = (P.fromIntegral :: Int -> Word8) . ceiling;+     "NP.round       :: a -> Word16"  round       = (P.fromIntegral :: Int -> Word16) . round;+     "NP.roundSimple :: a -> Word16"  roundSimple = (P.fromIntegral :: Int -> Word16) . roundSimple;+     "NP.truncate    :: a -> Word16"  truncate    = (P.fromIntegral :: Int -> Word16) . truncate;+     "NP.floor       :: a -> Word16"  floor       = (P.fromIntegral :: Int -> Word16) . floor;+     "NP.ceiling     :: a -> Word16"  ceiling     = (P.fromIntegral :: Int -> Word16) . ceiling;+     "NP.round       :: a -> Word32"  round       = (P.fromIntegral :: Int -> Word32) . round;+     "NP.roundSimple :: a -> Word32"  roundSimple = (P.fromIntegral :: Int -> Word32) . roundSimple;+     "NP.truncate    :: a -> Word32"  truncate    = (P.fromIntegral :: Int -> Word32) . truncate;+     "NP.floor       :: a -> Word32"  floor       = (P.fromIntegral :: Int -> Word32) . floor;+     "NP.ceiling     :: a -> Word32"  ceiling     = (P.fromIntegral :: Int -> Word32) . ceiling;+     "NP.round       :: a -> Word64"  round       = (P.fromIntegral :: Int -> Word64) . round;+     "NP.roundSimple :: a -> Word64"  roundSimple = (P.fromIntegral :: Int -> Word64) . roundSimple;+     "NP.truncate    :: a -> Word64"  truncate    = (P.fromIntegral :: Int -> Word64) . truncate;+     "NP.floor       :: a -> Word64"  floor       = (P.fromIntegral :: Int -> Word64) . floor;+     "NP.ceiling     :: a -> Word64"  ceiling     = (P.fromIntegral :: Int -> Word64) . ceiling;++     "NP.splitFraction :: a -> (Int8,a)"     splitFraction = mapFst (P.fromIntegral :: Int -> Int8) . splitFraction;+     "NP.splitFraction :: a -> (Int16,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Int16) . splitFraction;+     "NP.splitFraction :: a -> (Int32,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Int32) . splitFraction;+     "NP.splitFraction :: a -> (Int64,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Int64) . splitFraction;+     "NP.splitFraction :: a -> (Word,a)"     splitFraction = mapFst (P.fromIntegral :: Int -> Word) . splitFraction;+     "NP.splitFraction :: a -> (Word8,a)"    splitFraction = mapFst (P.fromIntegral :: Int -> Word8) . splitFraction;+     "NP.splitFraction :: a -> (Word16,a)"   splitFraction = mapFst (P.fromIntegral :: Int -> Word16) . splitFraction;+     "NP.splitFraction :: a -> (Word32,a)"   splitFraction = mapFst (P.fromIntegral :: Int -> Word32) . splitFraction;+     "NP.splitFraction :: a -> (Word64,a)"   splitFraction = mapFst (P.fromIntegral :: Int -> Word64) . splitFraction;+  #-}+++{- | TODO: Should be moved to a continued fraction module. -}++approxRational :: (ToRational.C a, C a) => a -> a -> Rational+approxRational rat eps    =  simplest (rat-eps) (rat+eps)+        where simplest x y | y < x      =  simplest y x+                           | x == y     =  xr+                           | x > 0      =  simplest' n d n' d'+                           | y < 0      =  - simplest' (-n') d' (-n) d+                           | otherwise  =  0 :% 1+                                        where xr@(n:%d) = ToRational.toRational x+                                              (n':%d')  = ToRational.toRational y++              simplest' n d n' d'       -- assumes 0 < n%d < n'%d'+                        | isZero r   =  q :% 1+                        | q /= q'    =  (q+1) :% 1+                        | otherwise  =  (q*n''+d'') :% n''+                                     where (q,r)      =  quotRem n d+                                           (q',r')    =  quotRem n' d'+                                           (n'':%d'') =  simplest' d' r' d r+++-- * generic implementation of round functions++powersOfTwo :: (Ring.C a) => [a]+powersOfTwo = iterate (2*) one++pairsOfPowersOfTwo :: (Ring.C a, Ring.C b) => [(a,b)]+pairsOfPowersOfTwo =+   zip powersOfTwo powersOfTwo++{- |+The generic rounding functions need a number of operations+proportional to the number of binary digits of the integer portion.+If operations like multiplication with two and comparison+need time proportional to the number of binary digits,+then the overall rounding requires quadratic time.+-}+genericFloor :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericFloor a =+   if a>=zero+     then genericPosFloor a+     else negate $ genericPosCeiling $ negate a++genericCeiling :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericCeiling a =+   if a>=zero+     then genericPosCeiling a+     else negate $ genericPosFloor $ negate a++genericTruncate :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericTruncate a =+   if a>=zero+     then genericPosFloor a+     else negate $ genericPosFloor $ negate a++genericRound :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericRound a =+   if a>=zero+     then genericPosRound a+     else negate $ genericPosRound $ negate a++genericFraction :: (Ord a, Ring.C a) => a -> a+genericFraction a =+   if a>=zero+     then genericPosFraction a+     else fixFraction $ negate $ genericPosFraction $ negate a++genericSplitFraction :: (Ord a, Ring.C a, Ring.C b) => a -> (b,a)+genericSplitFraction a =+   if a>=zero+     then genericPosSplitFraction a+     else fixSplitFraction $ mapPair (negate, negate) $+          genericPosSplitFraction $ negate a+++genericPosFloor :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericPosFloor a =+   snd $+   foldr+      (\(pa,pb) acc@(accA,accB) ->+         let newA = accA+pa+         in  if newA>a then acc else (newA,accB+pb))+      (zero,zero) $+   takeWhile ((a>=) . fst) $+   pairsOfPowersOfTwo++genericPosCeiling :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericPosCeiling a =+   snd $+   (\(ps,u:_) ->+      foldr+         (\(pa,pb) acc@(accA,accB) ->+            let newA = accA-pa+            in  if newA>=a then (newA,accB-pb) else acc)+         u ps) $+   span ((a>) . fst) $+   (zero,zero) : pairsOfPowersOfTwo++{-+genericPosFloorDigits :: (Ord a, Ring.C a, Ring.C b) => a -> ((a,b), [Bool])+genericPosFloorDigits a =+   List.mapAccumR+      (\acc@(accA,accB) (pa,pb) ->+         let newA = accA+pa+             b = newA<=a+         in  (if b then (newA,accB+pb) else acc, b))+      (zero,zero) $+   takeWhile ((a>=) . fst) $+   pairsOfPowersOfTwo+-}++genericHalfPosFloorDigits :: (Ord a, Ring.C a, Ring.C b) => a -> ((a,b), [Bool])+genericHalfPosFloorDigits a =+   List.mapAccumR+      (\acc@(accA,accB) (pa,pb) ->+         let newA = accA+pa+             b = newA<=a+         in  (if b then (newA,accB+pb) else acc, b))+      (zero,zero) $+   takeWhile ((a>=) . fst) $+   zip powersOfTwo (zero:powersOfTwo)++genericPosRound :: (Ord a, Ring.C a, Ring.C b) => a -> b+genericPosRound a =+   let a2 = 2*a+       ((ai,bi), ds) = genericHalfPosFloorDigits a2+   in  if ai==a2+         then+           case ds of+             True : True : _ -> bi+one+             _ -> bi+         else+           case ds of+             True : _ -> bi+one+             _ -> bi++genericPosFraction :: (Ord a, Ring.C a) => a -> a+genericPosFraction a =+   foldr+      (\p acc ->+         if p>acc then acc else acc-p)+      a $+   takeWhile (a>=) $+   powersOfTwo++genericPosSplitFraction :: (Ord a, Ring.C a, Ring.C b) => a -> (b,a)+genericPosSplitFraction a =+   foldr+      (\(pb,pa) acc@(accB,accA) ->+         if pa>accA then acc else (accB+pb,accA-pa))+      (zero,a) $+   takeWhile ((a>=) . snd) $+   pairsOfPowersOfTwo+++{- |+Needs linear time with respect to the number of digits.++This and other functions using OrderDecision+like @floor@ where argument and result are the same+may be moved to a new module.+-}+decisionPosFraction :: (OrdDec.C a, Ring.C a) => a -> a+decisionPosFraction a0 =+   (\ps ->+      foldr+         (\p cont a ->+            (a<?one) a $ cont $+            (a>=?p) (a-p) a)+         (error "decisionPosFraction: end of list should never be reached")+         ps a0) $+   concatMap (reverse . flip take powersOfTwo) powersOfTwo++{-+Works but needs quadratic time with respect to the number of digits.+I feel that there must be something more efficient.+-}+decisionPosFractionSqrTime :: (OrdDec.C a, Ring.C a) => a -> a+decisionPosFractionSqrTime a0 =+   (\ps ->+      foldr+         (\p cont a ->+            (a<?one) a $ cont $+            (a>=?p) (a-p) a)+         (error "decisionPosFraction: end of list should never be reached")+         ps a0) $+   concatMap reverse $+   inits powersOfTwo
src/Algebra/RealTranscendental.hs view
@@ -12,7 +12,7 @@ import Data.Bool.HT (select, )  import qualified Prelude as P-import PreludeBase+import NumericPrelude.Base
src/Algebra/RightModule.hs view
@@ -6,7 +6,7 @@ import qualified Algebra.Ring     as Ring import qualified Algebra.Additive as Additive --- import NumericPrelude+-- import NumericPrelude.Numeric import qualified Prelude
src/Algebra/Ring.hs view
@@ -36,7 +36,7 @@ import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  ) import Data.Word (Word, Word8, Word16, Word32, Word64, ) -import PreludeBase+import NumericPrelude.Base import Prelude(Integer,Int,Float,Double) import qualified Data.Ratio as Ratio98 import qualified Prelude as P
src/Algebra/ToInteger.hs view
@@ -26,7 +26,7 @@ import Data.Word (Word, Word8, Word16, Word32, Word64, )  import qualified Prelude as P-import PreludeBase+import NumericPrelude.Base import Prelude (Int, Integer, Float, Double, )  @@ -41,7 +41,7 @@  Conversions must be lossless, that is, they do not round in any way.-For rounding see "Algebra.RealField".+For rounding see "Algebra.RealRing". With the instances for 'Prelude.Float' and 'Prelude.Double' we acknowledge that these types actually represent rationals rather than (approximated) real numbers.
src/Algebra/ToRational.hs view
@@ -1,7 +1,8 @@ {-# LANGUAGE NoImplicitPrelude #-} module Algebra.ToRational where -import qualified Algebra.Real           as Real+import qualified Algebra.Field    as Field+import qualified Algebra.Absolute as Absolute import Algebra.Field (fromRational, ) import Algebra.Ring (fromInteger, ) @@ -11,47 +12,89 @@ import Data.Word (Word, Word8, Word16, Word32, Word64, )  import qualified Prelude as P-import PreludeBase+import NumericPrelude.Base import Prelude(Int,Integer,Float,Double)  {- | This class allows lossless conversion from any representation of a rational to the fixed 'Rational' type. \"Lossless\" means - don't do any rounding.-For rounding see "Algebra.RealField".+For rounding see "Algebra.RealRing". With the instances for 'Float' and 'Double' we acknowledge that these types actually represent rationals rather than (approximated) real numbers.-However, this contradicts to the 'Algebra.Transcendental'+However, this contradicts to the 'Algebra.Transcendental' class.  Laws that must be satisfied by instances:  >  fromRational' . toRational === id -}-class (Real.C a) => C a where+class (Absolute.C a) => C a where    -- | Lossless conversion from any representation of a rational to 'Rational'    toRational :: a -> Rational  instance C Integer where-   {-#INLINE toRational #-}+   {-# INLINE toRational #-}    toRational = fromInteger  instance C Float where-   {-#INLINE toRational #-}+   {-# INLINE toRational #-}    toRational = fromRational . P.toRational  instance C Double where-   {-#INLINE toRational #-}+   {-# INLINE toRational #-}    toRational = fromRational . P.toRational -instance C Int     where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Int8    where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Int16   where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Int32   where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Int64   where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int    where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int8   where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int16  where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int32  where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Int64  where {-# INLINE toRational #-}; toRational = toRational . P.toInteger -instance C Word    where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Word8   where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Word16  where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Word32  where {-#INLINE toRational #-}; toRational = toRational . P.toInteger-instance C Word64  where {-#INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word   where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word8  where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word16 where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word32 where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+instance C Word64 where {-# INLINE toRational #-}; toRational = toRational . P.toInteger+++{- |+It should hold++> realToField = fromRational' . toRational++but it should be much more efficient for particular pairs of types,+such as converting 'Float' to 'Double'.+This achieved by optimizer rules.+-}+realToField :: (C a, Field.C b) => a -> b+realToField = Field.fromRational' . toRational++{-# RULES+     "NP.realToField :: Integer  -> Float "  realToField = P.realToFrac :: Integer  -> Float ;+     "NP.realToField :: Int      -> Float "  realToField = P.realToFrac :: Int      -> Float ;+     "NP.realToField :: Int8     -> Float "  realToField = P.realToFrac :: Int8     -> Float ;+     "NP.realToField :: Int16    -> Float "  realToField = P.realToFrac :: Int16    -> Float ;+     "NP.realToField :: Int32    -> Float "  realToField = P.realToFrac :: Int32    -> Float ;+     "NP.realToField :: Int64    -> Float "  realToField = P.realToFrac :: Int64    -> Float ;+     "NP.realToField :: Word     -> Float "  realToField = P.realToFrac :: Word     -> Float ;+     "NP.realToField :: Word8    -> Float "  realToField = P.realToFrac :: Word8    -> Float ;+     "NP.realToField :: Word16   -> Float "  realToField = P.realToFrac :: Word16   -> Float ;+     "NP.realToField :: Word32   -> Float "  realToField = P.realToFrac :: Word32   -> Float ;+     "NP.realToField :: Word64   -> Float "  realToField = P.realToFrac :: Word64   -> Float ;+     "NP.realToField :: Float    -> Float "  realToField = P.realToFrac :: Float    -> Float ;+     "NP.realToField :: Double   -> Float "  realToField = P.realToFrac :: Double   -> Float ;+     "NP.realToField :: Integer  -> Double"  realToField = P.realToFrac :: Integer  -> Double;+     "NP.realToField :: Int      -> Double"  realToField = P.realToFrac :: Int      -> Double;+     "NP.realToField :: Int8     -> Double"  realToField = P.realToFrac :: Int8     -> Double;+     "NP.realToField :: Int16    -> Double"  realToField = P.realToFrac :: Int16    -> Double;+     "NP.realToField :: Int32    -> Double"  realToField = P.realToFrac :: Int32    -> Double;+     "NP.realToField :: Int64    -> Double"  realToField = P.realToFrac :: Int64    -> Double;+     "NP.realToField :: Word     -> Double"  realToField = P.realToFrac :: Word     -> Double;+     "NP.realToField :: Word8    -> Double"  realToField = P.realToFrac :: Word8    -> Double;+     "NP.realToField :: Word16   -> Double"  realToField = P.realToFrac :: Word16   -> Double;+     "NP.realToField :: Word32   -> Double"  realToField = P.realToFrac :: Word32   -> Double;+     "NP.realToField :: Word64   -> Double"  realToField = P.realToFrac :: Word64   -> Double;+     "NP.realToField :: Float    -> Double"  realToField = P.realToFrac :: Float    -> Double;+     "NP.realToField :: Double   -> Double"  realToField = P.realToFrac :: Double   -> Double;+  #-}
src/Algebra/Transcendental.hs view
@@ -13,7 +13,7 @@ import Algebra.Additive  ((+), (-), negate)  import qualified Prelude as P-import PreludeBase+import NumericPrelude.Base   infixr 8  **, ^?
src/Algebra/Units.hs view
@@ -34,7 +34,7 @@  import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  ) -import PreludeBase+import NumericPrelude.Base import Prelude (Integer, Int) import qualified Prelude as P import Test.QuickCheck ((==>), Property)@@ -84,7 +84,7 @@  intQuery :: (P.Integral a, Ring.C a) => a -> Bool intQuery = flip elem [one, negate one]-{- constraint must be replaced by NumericPrelude.Real -}+{- constraint must be replaced by NumericPrelude.Absolute -} intAssociate, intStandard, intStandardInverse ::    (P.Integral a, Ring.C a, ZeroTestable.C a) => a -> a intAssociate = P.abs
src/Algebra/VectorSpace.hs view
@@ -8,7 +8,7 @@ import qualified Algebra.PrincipalIdealDomain as PID import qualified Number.Ratio   as Ratio --- import NumericPrelude+-- import NumericPrelude.Numeric import qualified Prelude as P
src/Algebra/ZeroTestable.hs view
@@ -8,7 +8,7 @@  -- import qualified Prelude as P import Prelude(Int,Integer,Float,Double)-import PreludeBase+import NumericPrelude.Base  {- | Maybe the naming should be according to Algebra.Unit:
src/MathObj/Algebra.hs view
@@ -24,7 +24,7 @@ import qualified Data.Map as Map import Data.List(intersperse) -import PreludeBase(Ord,Eq,{-Read,-}Show,(++),($),+import NumericPrelude.Base(Ord,Eq,{-Read,-}Show,(++),($),                    concat,map,show)
src/MathObj/DiscreteMap.hs view
@@ -42,7 +42,7 @@ import Data.Map (Map)  import qualified Prelude as P-import PreludeBase+import NumericPrelude.Base  -- FIXME: Should this be implemented by isZero? -- | Remove all zero values from the map.
src/MathObj/Gaussian/Bell.hs view
@@ -10,7 +10,7 @@  import qualified Algebra.Transcendental as Trans import qualified Algebra.Field          as Field-import qualified Algebra.Real           as Real+import qualified Algebra.Absolute           as Absolute import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive @@ -23,14 +23,14 @@ import Control.Monad (liftM4, )  -- import Prelude (($))-import NumericPrelude-import PreludeBase hiding (reverse, )+import NumericPrelude.Numeric+import NumericPrelude.Base hiding (reverse, )   data T a = Cons {amp :: a, c0, c1 :: Complex.T a, c2 :: a}    deriving (Eq, Show) -instance (Real.C a, Arbitrary a) => Arbitrary (T a) where+instance (Absolute.C a, Arbitrary a) => Arbitrary (T a) where    arbitrary =       liftM4          (\k a b c -> Cons k a b (1 + abs c))
src/MathObj/Gaussian/Example.hs view
@@ -24,7 +24,7 @@ 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.Absolute           as Absolute import qualified Algebra.Ring           as Ring -- import qualified Algebra.Additive       as Additive @@ -48,8 +48,8 @@ -- import System.Exit (ExitCode, )  -- import Prelude (($))-import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base import qualified Prelude as P
src/MathObj/Gaussian/Polynomial.hs view
@@ -23,14 +23,15 @@ import qualified MathObj.Gaussian.Bell as Bell  import qualified MathObj.LaurentPolynomial as LPoly-import qualified MathObj.Polynomial as Poly+import qualified MathObj.Polynomial.Core   as PolyCore+import qualified MathObj.Polynomial        as Poly import qualified Number.Complex     as Complex  import qualified Algebra.ZeroTestable   as ZeroTestable 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.Absolute           as Absolute import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive @@ -48,15 +49,15 @@ import Test.QuickCheck (Arbitrary, arbitrary, ) import Control.Monad (liftM2, ) -import NumericPrelude-import PreludeBase hiding (reverse, )+import NumericPrelude.Numeric+import NumericPrelude.Base hiding (reverse, ) -- import Prelude ()   data T a = Cons {bell :: Bell.T a, polynomial :: Poly.T (Complex.T a)}    deriving (Show) -instance Real.C a => Eq (T a) where+instance (Ring.C a, Ord a) => Eq (T a) where    (==) = equal  {-@@ -64,7 +65,7 @@ We have to combine the amplitude of the bell with the polynomial, respecting signs and the square root of the bell amplitude. -}-equal :: Real.C a => T a -> T a -> Bool+equal :: (Ring.C a, Ord a) => T a -> T a -> Bool equal x y =    let bx = bell x        by = bell y@@ -84,7 +85,7 @@        scaleSqr by x == scaleSqr bx y  -instance (Real.C a, Arbitrary a) => Arbitrary (T a) where+instance (Absolute.C a, Arbitrary a) => Arbitrary (T a) where    arbitrary = --      liftM2 Cons arbitrary arbitrary       liftM2 Cons@@ -145,7 +146,7 @@    nest n (scale (-1/4) . differentiate) $    Cons (Bell.Cons one zero zero 2) one -eigenfunctionIterative :: (Field.C a, Real.C a) => Int -> T a+eigenfunctionIterative :: (Field.C a, Ord a) => Int -> T a eigenfunctionIterative n =    fst . head . dropWhile (uncurry (/=)) . mapAdjacent (,) $    eigenfunctionIteration $@@ -217,7 +218,7 @@ integrate f =    let fs = Poly.coeffs $ polynomial f        (ys,~[r]) =-          Poly.divModRev+          PolyCore.divModRev              {-              We need the shortening convention of 'zipWith'              in order to limit the result list,
src/MathObj/Gaussian/Variance.hs view
@@ -14,7 +14,7 @@ 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.Absolute           as Absolute import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive @@ -27,14 +27,14 @@ import Control.Monad (liftM2, )  -- import Prelude (($))-import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base   data T a = Cons {amp, c :: a}    deriving (Eq, Show) -instance (Real.C a, Arbitrary a) => Arbitrary (T a) where+instance (Absolute.C a, Arbitrary a) => Arbitrary (T a) where    arbitrary =       liftM2 Cons          arbitrary@@ -56,15 +56,15 @@    Poly.fromCoeffs [zero, zero, c f]  -norm1 :: (Algebraic.C a, Real.C a) => T a -> a+norm1 :: (Algebraic.C a, Absolute.C a) => T a -> a norm1 f =    sqrt $ abs (amp f) / c f -norm2 :: (Algebraic.C a, Real.C a) => T a -> a+norm2 :: (Algebraic.C a, Absolute.C a) => T a -> a norm2 f =    sqrt $ abs (amp f) / (sqrt $ 2 * c f) -normP :: (Trans.C a, Real.C a) => a -> T a -> a+normP :: (Trans.C a, Absolute.C a) => a -> T a -> a normP p f =    sqrt (abs (amp f)) * (p * c f) ^? (- recip (2*p))
src/MathObj/LaurentPolynomial.hs view
@@ -14,6 +14,7 @@  import qualified MathObj.Polynomial  as Poly import qualified MathObj.PowerSeries as PS+import qualified MathObj.PowerSeries.Core as PSCore  import qualified Algebra.VectorSpace  as VectorSpace import qualified Algebra.Module       as Module@@ -28,11 +29,11 @@ import Algebra.ZeroTestable(isZero) import Algebra.Module((*>)) -import qualified PreludeBase as P-import qualified NumericPrelude as NP+import qualified NumericPrelude.Base as P+import qualified NumericPrelude.Numeric as NP -import PreludeBase    hiding (const, reverse, )-import NumericPrelude hiding (div, negate, )+import NumericPrelude.Base    hiding (const, reverse, )+import NumericPrelude.Numeric hiding (div, negate, )  import qualified Data.List as List import Data.List.HT (mapAdjacent)@@ -163,7 +164,7 @@  addShifted :: Additive.C a => Int -> [a] -> [a] -> [a] addShifted del px py =-   let recurse 0 x      = PS.add x py+   let recurse 0 x      = PSCore.add x py        recurse d []     = replicate d zero ++ py        recurse d (x:xs) = x : recurse (d-1) xs    in  if del >= 0@@ -186,9 +187,6 @@  {- * Module -} -scale :: Ring.C a => a -> [a] -> [a]-scale = Poly.scale- instance Vector.C T where    zero  = zero    (<+>) = (+)@@ -203,7 +201,7 @@ {- * Ring -}  mul :: Ring.C a => T a -> T a -> T a-mul (Cons xt x) (Cons yt y) = Cons (xt+yt) (PS.mul x y)+mul (Cons xt x) (Cons yt y) = Cons (xt+yt) (PSCore.mul x y)  instance (Ring.C a) => Ring.C (T a) where     one           = const one@@ -218,7 +216,7 @@    let (xzero,x) = span isZero xs        (yzero,y) = span isZero ys    in  Cons (xt - yt + length xzero - length yzero)-            (PS.divide x y)+            (PSCore.divide x y)  instance (Field.C a, ZeroTestable.C a) => Field.C (T a) where    (/) = div
src/MathObj/Matrix.hs view
@@ -61,8 +61,8 @@ import Data.Tuple.HT (swap, mapFst, ) import Data.List.HT (outerProduct, ) -import NumericPrelude (Int, )-import PreludeBase hiding (zipWith, )+import NumericPrelude.Numeric (Int, )+import NumericPrelude.Base hiding (zipWith, )   {- |
src/MathObj/Monoid.hs view
@@ -7,7 +7,7 @@ import Algebra.Additive (zero, ) import Algebra.Monoid (C, idt, (<*>), ) -import PreludeBase+import NumericPrelude.Base  {- | It is only a monoid for non-negative numbers.
src/MathObj/PartialFraction.hs view
@@ -39,9 +39,9 @@ import Data.List.HT (dropWhileRev, ) import Data.List (group, sortBy, mapAccumR, ) -import PreludeBase hiding (zipWith)+import NumericPrelude.Base hiding (zipWith) -import NumericPrelude(Int, fromInteger)+import NumericPrelude.Numeric(Int, fromInteger)
src/MathObj/Permutation.hs view
@@ -16,8 +16,8 @@  import Data.Array(Ix) --- import NumericPrelude (Integer)-import PreludeBase+-- import NumericPrelude.Numeric (Integer)+import NumericPrelude.Base   {- |
src/MathObj/Permutation/CycleList.hs view
@@ -19,8 +19,8 @@  import qualified Data.List.Match as Match import Data.Maybe.HT (toMaybe)-import NumericPrelude (fromInteger)-import PreludeBase+import NumericPrelude.Numeric (fromInteger)+import NumericPrelude.Base   type Cycle i = [i]
src/MathObj/Permutation/CycleList/Check.hs view
@@ -24,8 +24,8 @@ import Data.Array((!), Ix) import qualified Data.Array as Array --- import NumericPrelude (Integer)-import PreludeBase hiding (cycle)+-- import NumericPrelude.Numeric (Integer)+import NumericPrelude.Base hiding (cycle)  {- | We shall make a little bit of a hack here, enabling us to use additive
src/MathObj/Permutation/Table.hs view
@@ -23,8 +23,8 @@ import Data.Tuple.HT (swap, ) import Data.Maybe.HT (toMaybe, ) --- import NumericPrelude (Integer)-import PreludeBase hiding (cycle)+-- import NumericPrelude.Numeric (Integer)+import NumericPrelude.Base hiding (cycle)   type T i = Array i i
src/MathObj/Polynomial.hs view
@@ -44,20 +44,17 @@ -}  module MathObj.Polynomial-   (T, fromCoeffs, coeffs,+   (T, fromCoeffs, coeffs, degree,     showsExpressionPrec, const,     evaluate, evaluateCoeffVector, evaluateArgVector,-    compose, equal, add, sub, negate,-    horner, hornerCoeffVector, hornerArgVector,-    shift, unShift,-    mul, scale, divMod, divModRev,-    tensorProduct, tensorProductAlt,-    mulShear, mulShearTranspose,-    progression, differentiate, integrate, integrateInt,-    fromRoots, alternate, reverse,+    collinear,+    integrate,+    compose, fromRoots, reverse,     translate, dilate, shrink, ) where +import qualified MathObj.Polynomial.Core as Core+ import qualified Algebra.Differential         as Differential import qualified Algebra.VectorSpace          as VectorSpace import qualified Algebra.Module               as Module@@ -76,19 +73,13 @@  import Control.Monad (liftM, ) import qualified Data.List as List-import NumericPrelude.List (zipWithOverlap, )-import Data.Tuple.HT (mapPair, mapFst, forcePair, )-import Data.List.HT-          (dropWhileRev, switchL, shear, shearTranspose, outerProduct, )  import Test.QuickCheck (Arbitrary(arbitrary)) -import qualified Prelude     as P98-import qualified PreludeBase as P-import qualified NumericPrelude as NP+import NumericPrelude.Base    hiding (const, reverse, )+import NumericPrelude.Numeric -import PreludeBase    hiding (const, reverse, )-import NumericPrelude hiding (divMod, negate, stdUnit, )+import qualified Prelude as P98   newtype T a = Cons {coeffs :: [a]}@@ -110,13 +101,19 @@ lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a) lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1) +degree :: (ZeroTestable.C a) => T a -> Maybe Int+degree x =+   case Core.normalize (coeffs x) of+      [] -> Nothing+      (_:xs) -> Just $ length xs+ {- Functor instance is e.g. useful for showing polynomials in residue rings. @fmap (ResidueClass.concrete 7) (polynomial [1,4,4::ResidueClass.T Integer] * polynomial [1,5,6])@ -}  instance Functor T where-  fmap f (Cons xs) = Cons (map f xs)+   fmap f (Cons xs) = Cons (map f xs)  {-# INLINE plusPrec #-} {-# INLINE appPrec #-}@@ -125,8 +122,8 @@ appPrec  = 10  instance (Show a) => Show (T a) where-  showsPrec p (Cons xs) =-    showParen (p >= appPrec) (showString "Polynomial.fromCoeffs " . shows xs)+   showsPrec p (Cons xs) =+      showParen (p >= appPrec) (showString "Polynomial.fromCoeffs " . shows xs)  {-# INLINE showsExpressionPrec #-} showsExpressionPrec :: (Show a, ZeroTestable.C a, Additive.C a) =>@@ -147,28 +144,10 @@            (foldl (.) id $ List.intersperse (showString " + ") $             map (uncurry showsTerm) terms) -{- |-Horner's scheme for evaluating a polynomial in a ring.--}-{-# INLINE horner #-}-horner :: Ring.C a => a -> [a] -> a-horner x = foldr (\c val -> c+x*val) zero -{- |-Horner's scheme for evaluating a polynomial in a module.--}-{-# INLINE hornerCoeffVector #-}-hornerCoeffVector :: Module.C a v => a -> [v] -> v-hornerCoeffVector x = foldr (\c val -> c+x*>val) zero--{-# INLINE hornerArgVector #-}-hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v-hornerArgVector x = foldr (\c val -> c*>one+val*x) zero-- {-# INLINE evaluate #-} evaluate :: Ring.C a => T a -> a -> a-evaluate (Cons y) x = horner x y+evaluate (Cons y) x = Core.horner x y  {- | Here the coefficients are vectors,@@ -176,7 +155,7 @@ -} {-# INLINE evaluateCoeffVector #-} evaluateCoeffVector :: Module.C a v => T v -> a -> v-evaluateCoeffVector (Cons y) x = hornerCoeffVector x y+evaluateCoeffVector (Cons y) x = Core.hornerCoeffVector x y  {- | Here the argument is a vector,@@ -185,7 +164,7 @@ -} {-# INLINE evaluateArgVector #-} evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v-evaluateArgVector (Cons y) x = hornerArgVector x y+evaluateArgVector (Cons y) x = Core.hornerArgVector x y  {- | 'compose' is the functional composition of polynomials.@@ -195,70 +174,35 @@ -}  -- compose :: Module.C a b => T b -> T a -> T a--- compose (Cons x) y = horner y (map const x)+-- compose (Cons x) y = Core.horner y (map const x) {-# INLINE compose #-} compose :: (Ring.C a) => T a -> T a -> T a-compose (Cons x) y = horner y (map const x)--{- |-It's also helpful to put a polynomial in canonical form.-'normalize' strips leading coefficients that are zero.--}--{-# INLINE normalize #-}-normalize :: (ZeroTestable.C a) => [a] -> [a]-normalize = dropWhileRev isZero--{- |-Multiply by the variable, used internally.--}--{-# INLINE shift #-}-shift :: (Additive.C a) => [a] -> [a]-shift [] = []-shift l  = zero : l--{-# INLINE unShift #-}-unShift :: [a] -> [a]-unShift []     = []-unShift (_:xs) = xs+compose (Cons x) y = Core.horner y (map const x)  {-# INLINE const #-} const :: a -> T a const x = lift0 [x] -{-# INLINE equal #-}-equal :: (Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool-equal x y = and (zipWithOverlap isZero isZero (==) x y) +collinear :: (Eq a, Ring.C a) => T a -> T a -> Bool+collinear (Cons x) (Cons y) = Core.collinear x y++ instance (Eq a, ZeroTestable.C a) => Eq (T a) where-  (Cons x) == (Cons y) = equal x y+   (Cons x) == (Cons y) = Core.equal x y  instance (Indexable.C a, ZeroTestable.C a) => Indexable.C (T a) where-  compare = Indexable.liftCompare coeffs+   compare = Indexable.liftCompare coeffs  instance (ZeroTestable.C a) => ZeroTestable.C (T a) where-  isZero (Cons x) = isZero x+   isZero (Cons x) = isZero x  -add, sub :: (Additive.C a) => [a] -> [a] -> [a]-add = (+)-sub = (-)--{-# INLINE negate #-}-negate :: (Additive.C a) => [a] -> [a]-negate = map NP.negate- instance (Additive.C a) => Additive.C (T a) where-  (+)    = lift2 add-  (-)    = lift2 sub-  zero   = lift0 []-  negate = lift1 negate---{-# INLINE scale #-}-scale :: Ring.C a => a -> [a] -> [a]-scale s = map (s*)+   (+)    = lift2 Core.add+   (-)    = lift2 Core.sub+   zero   = lift0 []+   negate = lift1 Core.negate   instance Vector.C T where@@ -272,81 +216,29 @@ instance (Field.C a, Module.C a b) => VectorSpace.C a (T b)  -{-# INLINE tensorProduct #-}-tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]-tensorProduct = outerProduct (*)--tensorProductAlt :: Ring.C a => [a] -> [a] -> [[a]]-tensorProductAlt xs ys = map (flip scale ys) xs--{- |-'mul' is fast if the second argument is a short polynomial,-'MathObj.PowerSeries.**' relies on that fact.--}--{-# INLINE mul #-}-mul :: Ring.C a => [a] -> [a] -> [a]-{- prevent from generation of many zeros-   if the first operand is the empty list -}-mul [] = P.const []-mul xs = foldr (\y zs -> let (v:vs) = scale y xs in v : add vs zs) []--- this one fails on infinite lists---    mul xs = foldr (\y zs -> add (scale y xs) (shift zs)) []--{-# INLINE mulShear #-}-mulShear :: Ring.C a => [a] -> [a] -> [a]-mulShear xs ys = map sum (shear (tensorProduct xs ys))--{-# INLINE mulShearTranspose #-}-mulShearTranspose :: Ring.C a => [a] -> [a] -> [a]-mulShearTranspose xs ys = map sum (shearTranspose (tensorProduct xs ys))- instance (Ring.C a) => Ring.C (T a) where-  one         = const one-  fromInteger = const . fromInteger-  (*)         = lift2 mul+   one         = const one+   fromInteger = const . fromInteger+   (*)         = lift2 Core.mul  -divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])-divMod x y =-   mapPair (List.reverse, List.reverse) $-   divModRev (List.reverse x) (List.reverse y)--{--snd $ Poly.divMod (repeat (1::Double)) [1,1]+{- |+The 'Integral.C' instance is intensionally built+from the 'Field.C' structure of the polynomial coefficients.+If we would use @Integral.C a@ superclass,+then the Euclidean algorithm could not determine+the greatest common divisor of e.g. @[1,1]@ and @[2]@. -}-divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])-divModRev x y =-   let (y0:ys) = dropWhile isZero y-       -- the second parameter represents lazily (length x - length y)-       aux xs' =-         forcePair .-         switchL-           ([], xs')-           (P.const $-              let (x0:xs) = xs'-                  q0      = x0/y0-              in  mapFst (q0:) . aux (sub xs (scale q0 ys)))-   in  if isZero y-         then error "MathObj.Polynomial: division by zero"-         else aux x (drop (length y - 1) x)- instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where-  divMod (Cons x) (Cons y) =-     let (d,m) = divMod x y-     in  (Cons d, Cons m)--{-# INLINE stdUnit #-}-stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a-stdUnit x = case normalize x of-    [] -> one-    l  -> last l+   divMod (Cons x) (Cons y) =+      let (d,m) = Core.divMod x y+      in  (Cons d, Cons m)  instance (ZeroTestable.C a, Field.C a) => Units.C (T a) where-  isUnit (Cons []) = False-  isUnit (Cons (x0:xs)) = not (isZero x0) && all isZero xs-  stdUnit    (Cons x) = const        (stdUnit x)-  stdUnitInv (Cons x) = const (recip (stdUnit x))+   isUnit (Cons []) = False+   isUnit (Cons (x0:xs)) = not (isZero x0) && all isZero xs+   stdUnit    (Cons x) = const        (Core.stdUnit x)+   stdUnitInv (Cons x) = const (recip (Core.stdUnit x))  {- Polynomials are a Euclidean domain, so no instance is necessary@@ -355,53 +247,28 @@  instance (ZeroTestable.C a, Field.C a) => PID.C (T a) -{-# INLINE progression #-}-progression :: Ring.C a => [a]-progression = iterate (one+) one -{-# INLINE differentiate #-}-differentiate :: (Ring.C a) => [a] -> [a]-differentiate = zipWith (*) progression . drop 1+instance (Ring.C a) => Differential.C (T a) where+   differentiate = lift1 Core.differentiate + {-# INLINE integrate #-}-integrate :: (Field.C a) => a -> [a] -> [a]-integrate c x = c : zipWith (/) x progression+integrate :: (Field.C a) => a -> T a -> T a+integrate = lift1 . Core.integrate -{- |-Integrates if it is possible to represent the integrated polynomial-in the given ring.-Otherwise undefined coefficients occur.--}-{-# INLINE integrateInt #-}-integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a]-integrateInt c x =-   c : zipWith Integral.safeDiv x progression  -instance (Ring.C a) => Differential.C (T a) where-  differentiate = lift1 differentiate-- {-# INLINE fromRoots #-} fromRoots :: (Ring.C a) => [a] -> T a-fromRoots = Cons . foldl (flip mulLinearFactor) [1]--{-# INLINE mulLinearFactor #-}-mulLinearFactor :: Ring.C a => a -> [a] -> [a]-mulLinearFactor x yt@(y:ys) = Additive.negate (x*y) : yt - scale x ys-mulLinearFactor _ [] = []--{-# INLINE alternate #-}-alternate :: Additive.C a => [a] -> [a]-alternate = zipWith ($) (cycle [id, Additive.negate])+fromRoots = Cons . foldl (flip Core.mulLinearFactor) [one]  {-# INLINE reverse #-} reverse :: Additive.C a => T a -> T a-reverse = lift1 alternate+reverse = lift1 Core.alternate  translate :: Ring.C a => a -> T a -> T a translate d =-   lift1 $ foldr (\c p -> [c] + mulLinearFactor d p) []+   lift1 $ foldr (\c p -> [c] + Core.mulLinearFactor d p) []  shrink :: Ring.C a => a -> T a -> T a shrink k =@@ -410,22 +277,9 @@ dilate :: Field.C a => a -> T a -> T a dilate = shrink . Field.recip -{--see htam: Wavelet/DyadicResultant -resultant :: Ring.C a => [a] -> [a] -> [a]-resultant xs ys =--discriminant :: Ring.C a => [a] -> a-discriminant xs =-   let degree = genericLength xs-   in  parityFlip (safeDiv (degree*(degree-1)) 2)-                  (resultant xs (differentiate xs))-          `safeDiv` last xs--}- instance (Arbitrary a, ZeroTestable.C a) => Arbitrary (T a) where-   arbitrary = liftM (fromCoeffs . normalize) arbitrary+   arbitrary = liftM (fromCoeffs . Core.normalize) arbitrary   {- * legacy instances -}
+ src/MathObj/Polynomial/Core.hs view
@@ -0,0 +1,228 @@+{-# LANGUAGE NoImplicitPrelude #-}+{- |+This module implements polynomial functions on plain lists.+We use such functions in order to implement methods of other datatypes.++The module organization differs from that of @ResidueClass@:+Here the @Polynomial@ module exports the type+that fits to the NumericPrelude type classes,+whereas in @ResidueClass@ the sub-modules export various flavors of them.+-}+module MathObj.Polynomial.Core (+   horner, hornerCoeffVector, hornerArgVector,+   normalize,+   shift, unShift,+   equal,+   add, sub, negate,+   scale, collinear,+   tensorProduct, tensorProductAlt,+   mul, mulShear, mulShearTranspose,+   divMod, divModRev,+   stdUnit,+   progression, differentiate, integrate, integrateInt,+   mulLinearFactor,+   alternate,+   ) where++import qualified Algebra.Module               as Module+import qualified Algebra.Field                as Field+import qualified Algebra.IntegralDomain       as Integral+import qualified Algebra.Ring                 as Ring+import qualified Algebra.Additive             as Additive+import qualified Algebra.ZeroTestable         as ZeroTestable++import Algebra.Module((*>))+import Algebra.ZeroTestable(isZero)++import qualified Data.List as List+import NumericPrelude.List (zipWithOverlap, )+import Data.Tuple.HT (mapPair, mapFst, forcePair, )+import Data.List.HT+          (dropWhileRev, switchL, shear, shearTranspose, outerProduct, )++import qualified Prelude     as P98+import qualified NumericPrelude.Base as P+import qualified NumericPrelude.Numeric as NP++import NumericPrelude.Base    hiding (const, reverse, )+import NumericPrelude.Numeric hiding (divMod, negate, stdUnit, )+++{- |+Horner's scheme for evaluating a polynomial in a ring.+-}+{-# INLINE horner #-}+horner :: Ring.C a => a -> [a] -> a+horner x = foldr (\c val -> c+x*val) zero++{- |+Horner's scheme for evaluating a polynomial in a module.+-}+{-# INLINE hornerCoeffVector #-}+hornerCoeffVector :: Module.C a v => a -> [v] -> v+hornerCoeffVector x = foldr (\c val -> c+x*>val) zero++{-# INLINE hornerArgVector #-}+hornerArgVector :: (Module.C a v, Ring.C v) => v -> [a] -> v+hornerArgVector x = foldr (\c val -> c*>one+val*x) zero+++{- |+It's also helpful to put a polynomial in canonical form.+'normalize' strips leading coefficients that are zero.+-}+{-# INLINE normalize #-}+normalize :: (ZeroTestable.C a) => [a] -> [a]+normalize = dropWhileRev isZero++{- |+Multiply by the variable, used internally.+-}+{-# INLINE shift #-}+shift :: (Additive.C a) => [a] -> [a]+shift [] = []+shift l  = zero : l++{-# INLINE unShift #-}+unShift :: [a] -> [a]+unShift []     = []+unShift (_:xs) = xs++{-# INLINE equal #-}+equal :: (Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool+equal x y = and (zipWithOverlap isZero isZero (==) x y)+++add, sub :: (Additive.C a) => [a] -> [a] -> [a]+add = (+)+sub = (-)++{-# INLINE negate #-}+negate :: (Additive.C a) => [a] -> [a]+negate = map NP.negate+++{-# INLINE scale #-}+scale :: Ring.C a => a -> [a] -> [a]+scale s = map (s*)+++collinear :: (Eq a, Ring.C a) => [a] -> [a] -> Bool+collinear (x:xs) (y:ys) =+   if x==zero && y==zero+     then collinear xs ys+     else scale x ys == scale y xs+-- here at least one of xs and ys is empty+collinear xs ys =+   all (==zero) xs && all (==zero) ys+++{-# INLINE tensorProduct #-}+tensorProduct :: Ring.C a => [a] -> [a] -> [[a]]+tensorProduct = outerProduct (*)++tensorProductAlt :: Ring.C a => [a] -> [a] -> [[a]]+tensorProductAlt xs ys = map (flip scale ys) xs+++{- |+'mul' is fast if the second argument is a short polynomial,+'MathObj.PowerSeries.**' relies on that fact.+-}++{-# INLINE mul #-}+mul :: Ring.C a => [a] -> [a] -> [a]+{- prevent from generation of many zeros+   if the first operand is the empty list -}+mul [] = P.const []+mul xs = foldr (\y zs -> let (v:vs) = scale y xs in v : add vs zs) []+-- this one fails on infinite lists+--    mul xs = foldr (\y zs -> add (scale y xs) (shift zs)) []++{-# INLINE mulShear #-}+mulShear :: Ring.C a => [a] -> [a] -> [a]+mulShear xs ys = map sum (shear (tensorProduct xs ys))++{-# INLINE mulShearTranspose #-}+mulShearTranspose :: Ring.C a => [a] -> [a] -> [a]+mulShearTranspose xs ys = map sum (shearTranspose (tensorProduct xs ys))+++divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])+divMod x y =+   mapPair (List.reverse, List.reverse) $+   divModRev (List.reverse x) (List.reverse y)++{-+snd $ Poly.divMod (repeat (1::Double)) [1,1]+-}+divModRev :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a], [a])+divModRev x y =+   let (y0:ys) = dropWhile isZero y+       -- the second parameter represents lazily (length x - length y)+       aux xs' =+         forcePair .+         switchL+           ([], xs')+           (P.const $+              let (x0:xs) = xs'+                  q0      = x0/y0+              in  mapFst (q0:) . aux (sub xs (scale q0 ys)))+   in  if isZero y+         then error "MathObj.Polynomial: division by zero"+         else aux x (drop (length y - 1) x)++{-# INLINE stdUnit #-}+stdUnit :: (ZeroTestable.C a, Ring.C a) => [a] -> a+stdUnit x = case normalize x of+    [] -> one+    l  -> last l+++{-# INLINE progression #-}+progression :: Ring.C a => [a]+progression = iterate (one+) one++{-# INLINE differentiate #-}+differentiate :: (Ring.C a) => [a] -> [a]+differentiate = zipWith (*) progression . drop 1++{-# INLINE integrate #-}+integrate :: (Field.C a) => a -> [a] -> [a]+integrate c x = c : zipWith (/) x progression++{- |+Integrates if it is possible to represent the integrated polynomial+in the given ring.+Otherwise undefined coefficients occur.+-}+{-# INLINE integrateInt #-}+integrateInt :: (ZeroTestable.C a, Integral.C a) => a -> [a] -> [a]+integrateInt c x =+   c : zipWith Integral.safeDiv x progression+++{-# INLINE mulLinearFactor #-}+mulLinearFactor :: Ring.C a => a -> [a] -> [a]+mulLinearFactor x yt@(y:ys) = Additive.negate (x*y) : yt - scale x ys+mulLinearFactor _ [] = []++{-# INLINE alternate #-}+alternate :: Additive.C a => [a] -> [a]+alternate = zipWith ($) (cycle [id, Additive.negate])+++{-+see htam: Wavelet/DyadicResultant++resultant :: Ring.C a => [a] -> [a] -> [a]+resultant xs ys =++discriminant :: Ring.C a => [a] -> a+discriminant xs =+   let degree = genericLength xs+   in  parityFlip (safeDiv (degree*(degree-1)) 2)+                  (resultant xs (differentiate xs))+          `safeDiv` last xs+-}+
src/MathObj/PowerSeries.hs view
@@ -3,13 +3,13 @@ {-# LANGUAGE FlexibleInstances #-}  {- |-Power series, either finite or unbounded.  (zipWith does exactly the-right thing to make it work almost transparently.)+Power series, either finite or unbounded.+(zipWith does exactly the right thing to make it work almost transparently.) -}- module MathObj.PowerSeries where -import qualified MathObj.Polynomial     as Poly+import qualified MathObj.PowerSeries.Core as Core+import qualified MathObj.Polynomial.Core as Poly  import qualified Algebra.Differential   as Differential import qualified Algebra.IntegralDomain as Integral@@ -26,14 +26,8 @@ import Algebra.Module((*>)) import Algebra.ZeroTestable(isZero) -import qualified Data.List.Match as Match-import qualified NumericPrelude as NP-import qualified PreludeBase as P--import PreludeBase    hiding (const)-import NumericPrelude hiding (negate, stdUnit, divMod,-                              sqrt, exp, log,-                              sin, cos, tan, asin, acos, atan)+import NumericPrelude.Base    hiding (const)+import NumericPrelude.Numeric   newtype T a = Cons {coeffs :: [a]} deriving (Ord)@@ -64,15 +58,15 @@ -}  instance Functor T where-  fmap f (Cons xs) = Cons (map f xs)+   fmap f (Cons xs) = Cons (map f xs)  {-# INLINE appPrec #-} appPrec :: Int appPrec  = 10  instance (Show a) => Show (T a) where-  showsPrec p (Cons xs) =-    showParen (p >= appPrec) (showString "PowerSeries.fromCoeffs " . shows xs)+   showsPrec p (Cons xs) =+     showParen (p >= appPrec) (showString "PowerSeries.fromCoeffs " . shows xs)   {-# INLINE truncate #-}@@ -82,132 +76,63 @@ {- | Evaluate (truncated) power series. -}-{-# INLINE eval #-}-eval :: Ring.C a => [a] -> a -> a-eval = flip Poly.horner- {-# INLINE evaluate #-} evaluate :: Ring.C a => T a -> a -> a-evaluate (Cons y) = eval y+evaluate (Cons y) = Core.evaluate y  {- | Evaluate (truncated) power series. -}-{-# INLINE evalCoeffVector #-}-evalCoeffVector :: Module.C a v => [v] -> a -> v-evalCoeffVector = flip Poly.hornerCoeffVector- {-# INLINE evaluateCoeffVector #-} evaluateCoeffVector :: Module.C a v => T v -> a -> v-evaluateCoeffVector (Cons y) = evalCoeffVector y+evaluateCoeffVector (Cons y) = Core.evaluateCoeffVector y  -{-# INLINE evalArgVector #-}-evalArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> v-evalArgVector = flip Poly.hornerArgVector- {-# INLINE evaluateArgVector #-} evaluateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> v-evaluateArgVector (Cons y) = evalArgVector y+evaluateArgVector (Cons y) = Core.evaluateArgVector y  {- | Evaluate approximations that is evaluate all truncations of the series. -}-{-# INLINE approx #-}-approx :: Ring.C a => [a] -> a -> [a]-approx y x =-   scanl (+) zero (zipWith (*) (iterate (x*) 1) y)- {-# INLINE approximate #-} approximate :: Ring.C a => T a -> a -> [a]-approximate (Cons y) = approx y+approximate (Cons y) = Core.approximate y   {- | Evaluate approximations that is evaluate all truncations of the series. -}-{-# INLINE approxCoeffVector #-}-approxCoeffVector :: Module.C a v => [v] -> a -> [v]-approxCoeffVector y x =-   scanl (+) zero (zipWith (*>) (iterate (x*) 1) y)- {-# INLINE approximateCoeffVector #-} approximateCoeffVector :: Module.C a v => T v -> a -> [v]-approximateCoeffVector (Cons y) = approxCoeffVector y+approximateCoeffVector (Cons y) = Core.approximateCoeffVector y   {- | Evaluate approximations that is evaluate all truncations of the series. -}-{-# INLINE approxArgVector #-}-approxArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> [v]-approxArgVector y x =-   scanl (+) zero (zipWith (*>) y (iterate (x*) 1))- {-# INLINE approximateArgVector #-} approximateArgVector :: (Module.C a v, Ring.C v) => T a -> v -> [v]-approximateArgVector (Cons y) = approxArgVector y---{- * Simple series manipulation -}--{- |-For the series of a real function @f@-compute the series for @\x -> f (-x)@--}--alternate :: Additive.C a => [a] -> [a]-alternate = zipWith id (cycle [id, NP.negate])--{- |-For the series of a real function @f@-compute the series for @\x -> (f x + f (-x)) \/ 2@--}--holes2 :: Additive.C a => [a] -> [a]-holes2 = zipWith id (cycle [id, P.const zero])--{- |-For the series of a real function @f@-compute the real series for @\x -> (f (i*x) + f (-i*x)) \/ 2@--}-holes2alternate :: Additive.C a => [a] -> [a]-holes2alternate =-   zipWith id (cycle [id, P.const zero, NP.negate, P.const zero])+approximateArgVector (Cons y) = Core.approximateArgVector y  -{- * Series arithmetic -}--add, sub :: (Additive.C a) => [a] -> [a] -> [a]-add = Poly.add-sub = Poly.sub--negate :: (Additive.C a) => [a] -> [a]-negate = Poly.negate--scale :: Ring.C a => a -> [a] -> [a]-scale = Poly.scale--mul :: Ring.C a => [a] -> [a] -> [a]-mul = Poly.mul- {- Note that the derived instances only make sense for finite series. -}  instance (Eq a, ZeroTestable.C a) => Eq (T a) where-    (Cons x) == (Cons y) = Poly.equal x y+   (Cons x) == (Cons y) = Poly.equal x y  instance (Additive.C a) => Additive.C (T a) where-    negate = lift1 Poly.negate-    (+)    = lift2 Poly.add-    (-)    = lift2 Poly.sub-    zero   = lift0 []+   negate = lift1 Poly.negate+   (+)    = lift2 Poly.add+   (-)    = lift2 Poly.sub+   zero   = lift0 []  instance (Ring.C a) => Ring.C (T a) where-    one           = const one-    fromInteger n = const (fromInteger n)-    (*)           = lift2 mul+   one           = const one+   fromInteger n = const (fromInteger n)+   (*)           = lift2 Core.mul  instance Vector.C T where    zero  = zero@@ -215,190 +140,45 @@    (*>)  = Vector.functorScale  instance (Module.C a b) => Module.C a (T b) where-    (*>) x = lift1 (x *>)+   (*>) x = lift1 (x *>)  instance (Field.C a, Module.C a b) => VectorSpace.C a (T b) -stripLeadZero :: (ZeroTestable.C a) => [a] -> [a] -> ([a],[a])-stripLeadZero (x:xs) (y:ys) =-  if isZero x && isZero y-    then stripLeadZero xs ys-    else (x:xs,y:ys)-stripLeadZero xs ys = (xs,ys) -{- |-Divide two series where the absolute term of the divisor is non-zero.-That is, power series with leading non-zero terms are the units-in the ring of power series.--Knuth: Seminumerical algorithms--}-divide :: (Field.C a) => [a] -> [a] -> [a]-divide (x:xs) (y:ys) =-   let zs = map (/y) (x : sub xs (mul zs ys))-   in  zs-divide [] _ = []-divide _ [] = error "PowerSeries.divide: division by empty series"--{- |-Divide two series also if the divisor has leading zeros.--}-divideStripZero :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> [a]-divideStripZero x' y' =-   let (x0,y0) = stripLeadZero x' y'-   in  if null y0 || isZero (head y0)-         then error "PowerSeries.divideStripZero: Division by zero."-         else divide x0 y0-- instance (Field.C a) => Field.C (T a) where-  (/) = lift2 divide+   (/) = lift2 Core.divide  -divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a],[a])-divMod xs ys =-   let (yZero,yRem) = span isZero ys-       (xMod, xRem) = Match.splitAt yZero xs-   in  (divide xRem yRem, xMod)- instance (ZeroTestable.C a, Field.C a) => Integral.C (T a) where-  divMod (Cons x) (Cons y) =-     let (d,m) = divMod x y-     in  (Cons d, Cons m)+   divMod (Cons x) (Cons y) =+      let (d,m) = Core.divMod x y+      in  (Cons d, Cons m)  -progression :: Ring.C a => [a]-progression = Poly.progression--recipProgression :: (Field.C a) => [a]-recipProgression = map recip progression--differentiate :: (Ring.C a) => [a] -> [a]-differentiate = Poly.differentiate--integrate :: (Field.C a) => a -> [a] -> [a]-integrate = Poly.integrate- instance (Ring.C a) => Differential.C (T a) where-  differentiate = lift1 differentiate+   differentiate = lift1 Core.differentiate  -{- |-We need to compute the square root only of the first term.-That is, if the first term is rational,-then all terms of the series are rational.--}--sqrt :: Field.C a => (a -> a) -> [a] -> [a]-sqrt _ [] = []-sqrt f0 (x:xs) =-   let y  = f0 x-       ys = map (/(y+y)) (xs - (0 : mul ys ys))-   in  y:ys--{--pow alpha t = t^alpha-(pow alpha . x)' = alpha * (pow (alpha-1) . x) * x'-alpha * (pow alpha . x) = x * x' * (pow alpha . x)'-y = pow alpha . x-alpha * y = x * x' * y'--}--{- |-Input series must start with non-zero term.--}-pow :: (Field.C a) => (a -> a) -> a -> [a] -> [a]-pow f0 expon x =-   let y  = integrate (f0 (head x)) y'-       y' = scale expon (divide y (mul x (differentiate x)))-   in  y- instance (Algebraic.C a) => Algebraic.C (T a) where-   sqrt   = lift1 (sqrt Algebraic.sqrt)-   x ^/ y = lift1 (pow (Algebraic.^/ y)+   sqrt   = lift1 (Core.sqrt Algebraic.sqrt)+   x ^/ y = lift1 (Core.pow (Algebraic.^/ y)                        (fromRational' y)) x -{- |-The first term needs a transcendent computation but the others do not.-That's why we accept a function which computes the first term. -> (exp . x)' =   (exp . x) * x'-> (sin . x)' =   (cos . x) * x'-> (cos . x)' = - (sin . x) * x'--}--exp :: Field.C a => (a -> a) -> [a] -> [a]-exp f0 x =-   let x' = differentiate x-       y  = integrate (f0 (head x)) (mul y x')-   in  y--sinCos :: Field.C a => (a -> (a,a)) -> [a] -> ([a],[a])-sinCos f0 x =-   let (y0Sin, y0Cos) = f0 (head x)-       x'   = differentiate x-       ySin = integrate y0Sin         (mul yCos x')-       yCos = integrate y0Cos (negate (mul ySin x'))-   in  (ySin, yCos)--sinCosScalar :: Transcendental.C a => a -> (a,a)-sinCosScalar x = (Transcendental.sin x, Transcendental.cos x)--sin, cos :: Field.C a => (a -> (a,a)) -> [a] -> [a]-sin f0 = fst . sinCos f0-cos f0 = snd . sinCos f0--tan :: (Field.C a) => (a -> (a,a)) -> [a] -> [a]-tan f0 = uncurry divide . sinCos f0--{--(log x)' == x'/x-(asin x)' == (acos x) == x'/sqrt(1-x^2)-(atan x)' == x'/(1+x^2)--}--{- |-Input series must start with non-zero term.--}-log :: (Field.C a) => (a -> a) -> [a] -> [a]-log f0 x = integrate (f0 (head x)) (derivedLog x)--{- |-Computes @(log x)'@, that is @x'\/x@--}-derivedLog :: (Field.C a) => [a] -> [a]-derivedLog x = divide (differentiate x) x--atan :: (Field.C a) => (a -> a) -> [a] -> [a]-atan f0 x =-   let x' = differentiate x-   in  integrate (f0 (head x)) (divide x' ([1] + mul x x))--asin, acos :: (Field.C a) =>-   (a -> a) -> (a -> a) -> [a] -> [a]-asin sqrt0 f0 x =-   let x' = differentiate x-   in  integrate (f0 (head x))-                 (divide x' (sqrt sqrt0 ([1] - mul x x)))-acos = asin---- instance (Transcendental.C a) =>              Transcendental.C (T a) where-   pi = const NP.pi-   exp = lift1 (exp Transcendental.exp)-   sin = lift1 (sin sinCosScalar)-   cos = lift1 (cos sinCosScalar)-   tan = lift1 (tan sinCosScalar)+   pi = const Transcendental.pi+   exp = lift1 (Core.exp Transcendental.exp)+   sin = lift1 (Core.sin Core.sinCosScalar)+   cos = lift1 (Core.cos Core.sinCosScalar)+   tan = lift1 (Core.tan Core.sinCosScalar)    x ** y = Transcendental.exp (Transcendental.log x * y)                 {- This order of multiplication is especially fast                    when y is a singleton. -}-   log  = lift1 (log  Transcendental.log)-   asin = lift1 (asin Algebraic.sqrt Transcendental.asin)-   acos = lift1 (acos Algebraic.sqrt Transcendental.acos)-   atan = lift1 (atan Transcendental.atan)+   log  = lift1 (Core.log  Transcendental.log)+   asin = lift1 (Core.asin Algebraic.sqrt Transcendental.asin)+   acos = lift1 (Core.acos Algebraic.sqrt Transcendental.acos)+   atan = lift1 (Core.atan Transcendental.atan)  {- | It fulfills@@ -410,46 +190,5 @@ compose (Cons (x:_)) (Cons []) = Cons [x] compose (Cons x) (Cons (y:ys)) =    if isZero y-     then Cons (comp x ys)+     then Cons (Core.compose x ys)      else error "PowerSeries.compose: inner series must not have an absolute term."--{- |-Since the inner series must start with a zero,-the first term is omitted in y.--}-comp :: (Ring.C a) => [a] -> [a] -> [a]-comp xs y = foldr (\x acc -> x : mul y acc) [] xs---{- |-Compose two power series where the outer series-can be developed for any expansion point.-To be more precise:-The outer series must be expanded with respect to the leading term-of the inner series.--}-composeTaylor :: Ring.C a => (a -> [a]) -> [a] -> [a]-composeTaylor x (y:ys) = comp (x y) ys-composeTaylor x []     = x 0----{--(x . y) = id-(x' . y) * y' = 1-y' = 1 / (x' . y)--}--{- |-This function returns the series of the function in the form:-(point of the expansion, power series)--This is exceptionally slow and needs cubic run-time.--}--inv :: (Field.C a) => [a] -> (a, [a])-inv x =-   let y' = divide [1] (comp (differentiate x) (tail y))-       y  = integrate 0 y'-            -- the first term is zero, which is required for composition-   in  (head x, y)
+ src/MathObj/PowerSeries/Core.hs view
@@ -0,0 +1,282 @@+{-# LANGUAGE NoImplicitPrelude #-}+module MathObj.PowerSeries.Core where++import qualified MathObj.Polynomial.Core as Poly++import qualified Algebra.Module         as Module+import qualified Algebra.Transcendental as Transcendental+import qualified Algebra.Field          as Field+import qualified Algebra.Ring           as Ring+import qualified Algebra.Additive       as Additive+import qualified Algebra.ZeroTestable   as ZeroTestable++import Algebra.Module((*>))+import Algebra.ZeroTestable(isZero)++import qualified Data.List.Match as Match+import qualified NumericPrelude.Numeric as NP+import qualified NumericPrelude.Base as P++import NumericPrelude.Base    hiding (const)+import NumericPrelude.Numeric hiding (negate, stdUnit, divMod,+                              sqrt, exp, log,+                              sin, cos, tan, asin, acos, atan)+++{-# INLINE evaluate #-}+evaluate :: Ring.C a => [a] -> a -> a+evaluate = flip Poly.horner++{-# INLINE evaluateCoeffVector #-}+evaluateCoeffVector :: Module.C a v => [v] -> a -> v+evaluateCoeffVector = flip Poly.hornerCoeffVector++{-# INLINE evaluateArgVector #-}+evaluateArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> v+evaluateArgVector = flip Poly.hornerArgVector+++{-# INLINE approximate #-}+approximate :: Ring.C a => [a] -> a -> [a]+approximate y x =+   scanl (+) zero (zipWith (*) (iterate (x*) 1) y)++{-# INLINE approximateCoeffVector #-}+approximateCoeffVector :: Module.C a v => [v] -> a -> [v]+approximateCoeffVector y x =+   scanl (+) zero (zipWith (*>) (iterate (x*) 1) y)++{-# INLINE approximateArgVector #-}+approximateArgVector :: (Module.C a v, Ring.C v) => [a] -> v -> [v]+approximateArgVector y x =+   scanl (+) zero (zipWith (*>) y (iterate (x*) 1))+++{- * Simple series manipulation -}++{- |+For the series of a real function @f@+compute the series for @\x -> f (-x)@+-}++alternate :: Additive.C a => [a] -> [a]+alternate = zipWith id (cycle [id, NP.negate])++{- |+For the series of a real function @f@+compute the series for @\x -> (f x + f (-x)) \/ 2@+-}++holes2 :: Additive.C a => [a] -> [a]+holes2 = zipWith id (cycle [id, P.const zero])++{- |+For the series of a real function @f@+compute the real series for @\x -> (f (i*x) + f (-i*x)) \/ 2@+-}+holes2alternate :: Additive.C a => [a] -> [a]+holes2alternate =+   zipWith id (cycle [id, P.const zero, NP.negate, P.const zero])+++{- * Series arithmetic -}++add, sub :: (Additive.C a) => [a] -> [a] -> [a]+add = Poly.add+sub = Poly.sub++negate :: (Additive.C a) => [a] -> [a]+negate = Poly.negate++scale :: Ring.C a => a -> [a] -> [a]+scale = Poly.scale++mul :: Ring.C a => [a] -> [a] -> [a]+mul = Poly.mul+++stripLeadZero :: (ZeroTestable.C a) => [a] -> [a] -> ([a],[a])+stripLeadZero (x:xs) (y:ys) =+  if isZero x && isZero y+    then stripLeadZero xs ys+    else (x:xs,y:ys)+stripLeadZero xs ys = (xs,ys)+++divMod :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> ([a],[a])+divMod xs ys =+   let (yZero,yRem) = span isZero ys+       (xMod, xRem) = Match.splitAt yZero xs+   in  (divide xRem yRem, xMod)++{- |+Divide two series where the absolute term of the divisor is non-zero.+That is, power series with leading non-zero terms are the units+in the ring of power series.++Knuth: Seminumerical algorithms+-}+divide :: (Field.C a) => [a] -> [a] -> [a]+divide (x:xs) (y:ys) =+   let zs = map (/y) (x : sub xs (mul zs ys))+   in  zs+divide [] _ = []+divide _ [] = error "PowerSeries.divide: division by empty series"++{- |+Divide two series also if the divisor has leading zeros.+-}+divideStripZero :: (ZeroTestable.C a, Field.C a) => [a] -> [a] -> [a]+divideStripZero x' y' =+   let (x0,y0) = stripLeadZero x' y'+   in  if null y0 || isZero (head y0)+         then error "PowerSeries.divideStripZero: Division by zero."+         else divide x0 y0+++progression :: Ring.C a => [a]+progression = Poly.progression++recipProgression :: (Field.C a) => [a]+recipProgression = map recip progression++differentiate :: (Ring.C a) => [a] -> [a]+differentiate = Poly.differentiate++integrate :: (Field.C a) => a -> [a] -> [a]+integrate = Poly.integrate+++{- |+We need to compute the square root only of the first term.+That is, if the first term is rational,+then all terms of the series are rational.+-}+sqrt :: Field.C a => (a -> a) -> [a] -> [a]+sqrt _ [] = []+sqrt f0 (x:xs) =+   let y  = f0 x+       ys = map (/(y+y)) (xs - (0 : mul ys ys))+   in  y:ys++{-+pow alpha t = t^alpha+(pow alpha . x)' = alpha * (pow (alpha-1) . x) * x'+alpha * (pow alpha . x) = x * x' * (pow alpha . x)'+y = pow alpha . x+alpha * y = x * x' * y'+-}++{- |+Input series must start with non-zero term.+-}+pow :: (Field.C a) => (a -> a) -> a -> [a] -> [a]+pow f0 expon x =+   let y  = integrate (f0 (head x)) y'+       y' = scale expon (divide y (mul x (differentiate x)))+   in  y+++{- |+The first term needs a transcendent computation but the others do not.+That's why we accept a function which computes the first term.++> (exp . x)' =   (exp . x) * x'+> (sin . x)' =   (cos . x) * x'+> (cos . x)' = - (sin . x) * x'+-}+exp :: Field.C a => (a -> a) -> [a] -> [a]+exp f0 x =+   let x' = differentiate x+       y  = integrate (f0 (head x)) (mul y x')+   in  y++sinCos :: Field.C a => (a -> (a,a)) -> [a] -> ([a],[a])+sinCos f0 x =+   let (y0Sin, y0Cos) = f0 (head x)+       x'   = differentiate x+       ySin = integrate y0Sin         (mul yCos x')+       yCos = integrate y0Cos (negate (mul ySin x'))+   in  (ySin, yCos)++sinCosScalar :: Transcendental.C a => a -> (a,a)+sinCosScalar x = (Transcendental.sin x, Transcendental.cos x)++sin, cos :: Field.C a => (a -> (a,a)) -> [a] -> [a]+sin f0 = fst . sinCos f0+cos f0 = snd . sinCos f0++tan :: (Field.C a) => (a -> (a,a)) -> [a] -> [a]+tan f0 = uncurry divide . sinCos f0++{-+(log x)' == x'/x+(asin x)' == (acos x) == x'/sqrt(1-x^2)+(atan x)' == x'/(1+x^2)+-}++{- |+Input series must start with non-zero term.+-}+log :: (Field.C a) => (a -> a) -> [a] -> [a]+log f0 x = integrate (f0 (head x)) (derivedLog x)++{- |+Computes @(log x)'@, that is @x'\/x@+-}+derivedLog :: (Field.C a) => [a] -> [a]+derivedLog x = divide (differentiate x) x++atan :: (Field.C a) => (a -> a) -> [a] -> [a]+atan f0 x =+   let x' = differentiate x+   in  integrate (f0 (head x)) (divide x' ([1] + mul x x))++asin, acos :: (Field.C a) =>+   (a -> a) -> (a -> a) -> [a] -> [a]+asin sqrt0 f0 x =+   let x' = differentiate x+   in  integrate (f0 (head x))+                 (divide x' (sqrt sqrt0 ([1] - mul x x)))+acos = asin++{- |+Since the inner series must start with a zero,+the first term is omitted in y.+-}+compose :: (Ring.C a) => [a] -> [a] -> [a]+compose xs y = foldr (\x acc -> x : mul y acc) [] xs+++{- |+Compose two power series where the outer series+can be developed for any expansion point.+To be more precise:+The outer series must be expanded with respect to the leading term+of the inner series.+-}+composeTaylor :: Ring.C a => (a -> [a]) -> [a] -> [a]+composeTaylor x (y:ys) = compose (x y) ys+composeTaylor x []     = x 0++++{-+(x . y) = id+(x' . y) * y' = 1+y' = 1 / (x' . y)+-}++{- |+This function returns the series of the function in the form:+(point of the expansion, power series)++This is exceptionally slow and needs cubic run-time.+-}++inv :: (Field.C a) => [a] -> (a, [a])+inv x =+   let y' = divide [1] (compose (differentiate x) (tail y))+       y  = integrate 0 y'+            -- the first term is zero, which is required for composition+   in  (head x, y)
src/MathObj/PowerSeries/DifferentialEquation.hs view
@@ -11,14 +11,14 @@  module MathObj.PowerSeries.DifferentialEquation where -import qualified MathObj.PowerSeries         as PS+import qualified MathObj.PowerSeries.Core    as PS import qualified MathObj.PowerSeries.Example as PSE  import qualified Algebra.Field        as Field import qualified Algebra.ZeroTestable as ZeroTestable -import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base   {- |
src/MathObj/PowerSeries/Example.hs view
@@ -1,7 +1,7 @@ {-# LANGUAGE NoImplicitPrelude #-} module MathObj.PowerSeries.Example where -import qualified MathObj.PowerSeries as PS+import qualified MathObj.PowerSeries.Core as PS  import qualified Algebra.Field          as Field import qualified Algebra.Ring           as Ring@@ -14,9 +14,9 @@ import Data.List (map, tail, cycle, zipWith, scanl, intersperse) import Data.List.HT (sieve) -import NumericPrelude (one, (*), (/),+import NumericPrelude.Numeric (one, (*), (/),                        fromInteger, {-fromRational,-} pi)-import PreludeBase -- (Bool, const, map, zipWith, id, (&&), (==))+import NumericPrelude.Base -- (Bool, const, map, zipWith, id, (&&), (==))   {- * Default implementations. -}
src/MathObj/PowerSeries/Mean.hs view
@@ -6,7 +6,9 @@ module MathObj.PowerSeries.Mean where  import qualified MathObj.PowerSeries2        as PS2+import qualified MathObj.PowerSeries2.Core   as PS2Core import qualified MathObj.PowerSeries         as PS+import qualified MathObj.PowerSeries.Core    as PSCore import qualified MathObj.PowerSeries.Example as PSE  import qualified Algebra.Field as Field@@ -14,8 +16,8 @@  import Data.List.HT (shearTranspose) -import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base  {- $M_f$ is a generalized $f$-mean (quasi-arithmetic) if@@ -73,8 +75,8 @@  diffComp :: (Ring.C a) => [a] -> [a] -> [a] diffComp ys x =-   map sum (shearTranspose (tail (zipWith PS.scale ys-                    (map tail (iterate (PS.mul x) [1])))))+   map sum (shearTranspose (tail (zipWith PSCore.scale ys+                    (map tail (iterate (PSCore.mul x) [1])))))  {- Now we solve@@ -87,7 +89,7 @@ logarithmic :: (Field.C a) => [a] logarithmic =    let -- series for \frac{2\cdot x}{\ln(1+2\cdot x)}-       fracLn = PS.divide [2]+       fracLn = PSCore.divide [2]                       (tail (zipWith (*) (iterate (2*) 1) PSE.log))        fDiffFracLn = diffComp f (tail fracLn)        f = 0 : 1 : zipWith (/) fDiffFracLn@@ -118,7 +120,7 @@ $M(1+t,1) = \sqrt{1+t+t^2/2}$ -} quadratic :: (Field.C a, Eq a) => [a]-quadratic = PS.sqrt (\1 -> 1) [1,1,1/2]+quadratic = PSCore.sqrt (\1 -> 1) [1,1,1/2]  quadraticMVF :: (Field.C a) => [a] quadraticMVF =@@ -126,8 +128,8 @@    -- [1,1,1,1,1/2,1/2]    [1,1,1,1,1/2,-1/14] --- map (\x -> PS.coeffs (meanValueDiff2 quadratic2 [1,1,1,1,1/2,x] !! 4) !! 2) (GNUPlot.linearScale 10 (-0.071429,-1/14::Double))--- take 20 $ Numerics.ZeroFinder.RegulaFalsi.zero (-1,0) (\x -> PS.coeffs (meanValueDiff2 quadratic2 [1::Double,1,1,1,1/2,x] !! 4) !! 2)+-- map (\x -> PSCore.coeffs (meanValueDiff2 quadratic2 [1,1,1,1,1/2,x] !! 4) !! 2) (GNUPlot.linearScale 10 (-0.071429,-1/14::Double))+-- take 20 $ Numerics.ZeroFinder.RegulaFalsi.zero (-1,0) (\x -> PSCore.coeffs (meanValueDiff2 quadratic2 [1::Double,1,1,1,1/2,x] !! 4) !! 2)  {- Result: It seems,@@ -139,8 +141,8 @@ quadraticDiff :: (Field.C a, Eq a) => [a] quadraticDiff =    let divDiffPS = tail quadraticMVF -- (f(1+t)-f(1))/((1+t)-1)-       (1, invPS) = PS.inv (PS.differentiate quadraticMVF)-       meanValuePS = PS.composeTaylor (\1 -> invPS) divDiffPS+       (1, invPS) = PSCore.inv (PSCore.differentiate quadraticMVF)+       meanValuePS = PSCore.composeTaylor (\1 -> invPS) divDiffPS        {- instead of computing an inverse series           we could also apply (compose) the derived series           to the series of the quadratic mean. -}@@ -151,11 +153,11 @@  $M(1+x,1+y) = \sqrt{1+x+y+(x^2+y^2)/2}$ -}-quadratic2 :: (Field.C a, Eq a) => PS2.Core a+quadratic2 :: (Field.C a, Eq a) => PS2Core.T a quadratic2 =-   PS2.sqrt (\1 -> 1) [[1],[1,1],[1/2,0,1/2]]+   PS2Core.sqrt (\1 -> 1) [[1],[1,1],[1/2,0,1/2]] -quadraticDiff2 :: (Field.C a, Eq a) => PS2.Core a+quadraticDiff2 :: (Field.C a, Eq a) => PS2Core.T a quadraticDiff2 =    meanValueDiff2 quadratic2 quadraticMVF @@ -176,15 +178,15 @@ {- $M(1+x,1+y) = 2/(recip (1+x) + recip (1+y))$ -}-harmonic2 :: (Field.C a, Eq a) => PS2.Core a+harmonic2 :: (Field.C a, Eq a) => PS2Core.T a harmonic2 =    let rec = PS.fromCoeffs PSE.recip-   in  PS2.divide [[2]] $+   in  PS2Core.divide [[2]] $        PS2.coeffs $           PS2.fromPowerSeries0 rec +           PS2.fromPowerSeries1 rec -harmonicDiff2 :: (Field.C a, Eq a) => PS2.Core a+harmonicDiff2 :: (Field.C a, Eq a) => PS2Core.T a harmonicDiff2 =    meanValueDiff2 harmonic2 harmonicMVF @@ -196,10 +198,10 @@ {- $M(1+x,1+y) = 1+x/2+y/2$ -}-arithmetic2 :: (Field.C a, Eq a) => PS2.Core a+arithmetic2 :: (Field.C a, Eq a) => PS2Core.T a arithmetic2 = [[1],[1/2,1/2]] -arithmeticDiff2 :: (Field.C a, Eq a) => PS2.Core a+arithmeticDiff2 :: (Field.C a, Eq a) => PS2Core.T a arithmeticDiff2 =    meanValueDiff2 arithmetic2 arithmeticMVF @@ -210,11 +212,11 @@ {- $M(1+x,1+y) = \sqrt{(1+x)·(1+y)}$ -}-geometric2 :: (Field.C a, Eq a) => PS2.Core a+geometric2 :: (Field.C a, Eq a) => PS2Core.T a geometric2 =-   PS2.sqrt (\1 -> 1) [[1],[1,1],[0,1,0]]+   PS2Core.sqrt (\1 -> 1) [[1],[1,1],[0,1,0]] -geometricDiff2 :: (Field.C a, Eq a) => PS2.Core a+geometricDiff2 :: (Field.C a, Eq a) => PS2Core.T a geometricDiff2 =    meanValueDiff2 geometric2 geometricMVF @@ -222,11 +224,11 @@   meanValueDiff2 :: (Field.C a, Eq a) =>-   PS2.Core a -> [a] -> PS2.Core a+   PS2Core.T a -> [a] -> PS2Core.T a meanValueDiff2 mean2 curve =    let -- (f(1+x)-f(1+y)) / (x-y)        divDiffPS =           zipWith replicate [1..] $ tail curve        meanValuePS =-          PS2.comp (PS.differentiate curve) (tail mean2)+          PS2Core.compose (PSCore.differentiate curve) (tail mean2)    in  meanValuePS - divDiffPS
src/MathObj/PowerSeries2.hs view
@@ -8,10 +8,10 @@  module MathObj.PowerSeries2 where -import qualified MathObj.PowerSeries    as PS-import qualified MathObj.Polynomial     as Poly+import qualified MathObj.PowerSeries2.Core as Core+import qualified MathObj.PowerSeries as PS+import qualified MathObj.Polynomial.Core as Poly -import qualified Algebra.Differential   as Differential import qualified Algebra.Vector         as Vector import qualified Algebra.Algebraic      as Algebraic import qualified Algebra.Field          as Field@@ -19,16 +19,14 @@ import qualified Algebra.Additive       as Additive import qualified Algebra.ZeroTestable   as ZeroTestable -import qualified NumericPrelude as NP-import qualified PreludeBase as P+import qualified NumericPrelude.Numeric as NP+import qualified NumericPrelude.Base as P  import Data.List (isPrefixOf, ) import qualified Data.List.Match as Match -import PreludeBase    hiding (const)-import NumericPrelude hiding (negate, stdUnit,-                              sqrt, exp, log,-                              sin, cos, tan, asin, acos, atan)+import NumericPrelude.Base    hiding (const)+import NumericPrelude.Numeric  {- | In order to handle both variables equivalently@@ -41,9 +39,8 @@ Although the sub-lists are always finite and thus are more like polynomials than power series, division and square root computation are easier to implement for power series. -}-newtype T a = Cons {coeffs :: Core a} deriving (Ord)+newtype T a = Cons {coeffs :: Core.T a} deriving (Ord) -type Core a = [[a]]  isValid :: [[a]] -> Bool isValid = flip isPrefixOf [1..] . map length@@ -73,71 +70,45 @@    map (:[]) (PS.coeffs x)  -lift0 :: Core a -> T a+lift0 :: Core.T a -> T a lift0 = Cons -lift1 :: (Core a -> Core a) -> (T a -> T a)+lift1 :: (Core.T a -> Core.T a) -> (T a -> T a) lift1 f (Cons x0) = Cons (f x0) -lift2 :: (Core a -> Core a -> Core a) -> (T a -> T a -> T a)+lift2 :: (Core.T a -> Core.T a -> Core.T a) -> (T a -> T a -> T a) lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1)  -lift0fromPowerSeries :: [PS.T a] -> Core a-lift0fromPowerSeries = map PS.coeffs--lift1fromPowerSeries :: ([PS.T a] -> [PS.T a]) -> (Core a -> Core a)-lift1fromPowerSeries f x0 = map PS.coeffs (f (map PS.fromCoeffs x0))--lift2fromPowerSeries :: ([PS.T a] -> [PS.T a] -> [PS.T a]) -> (Core a -> Core a -> Core a)-lift2fromPowerSeries f x0 x1 = map PS.coeffs (f (map PS.fromCoeffs x0) (map PS.fromCoeffs x1))-- const :: a -> T a const x = lift0 [[x]]   instance Functor T where-  fmap f (Cons xs) = Cons (map (map f) xs)+   fmap f (Cons xs) = Cons (map (map f) xs)  appPrec :: Int appPrec  = 10  instance (Show a) => Show (T a) where-  showsPrec p (Cons xs) =-    showParen (p >= appPrec) (showString "PowerSeries2.fromCoeffs " . shows xs)---{- * Series arithmetic -}--add, sub :: (Additive.C a) => Core a -> Core a -> Core a-add = PS.add-sub = PS.sub--negate :: (Additive.C a) => Core a -> Core a-negate = PS.negate+   showsPrec p (Cons xs) =+      showParen (p >= appPrec) (showString "PowerSeries2.fromCoeffs " . shows xs)   instance (Eq a, ZeroTestable.C a) => Eq (T a) where-    (Cons x) == (Cons y) = Poly.equal x y+   (Cons x) == (Cons y) = Poly.equal x y  instance (Additive.C a) => Additive.C (T a) where-    negate = lift1 PS.negate-    (+)    = lift2 PS.add-    (-)    = lift2 PS.sub-    zero   = lift0 []+   negate = lift1 Core.negate+   (+)    = lift2 Core.add+   (-)    = lift2 Core.sub+   zero   = lift0 []  -scale :: Ring.C a => a -> Core a -> Core a-scale = map . (Vector.*>)--mul :: Ring.C a => Core a -> Core a -> Core a-mul = lift2fromPowerSeries PS.mul- instance (Ring.C a) => Ring.C (T a) where-    one           = const one-    fromInteger n = const (fromInteger n)-    (*)           = lift2 mul+   one           = const one+   fromInteger n = const (fromInteger n)+   (*)           = lift2 Core.mul  instance Vector.C T where    zero  = zero@@ -145,50 +116,11 @@    (*>)  = Vector.functorScale  -divide :: (Field.C a) =>-   Core a -> Core a -> Core a-divide = lift2fromPowerSeries PS.divide-- instance (Field.C a) => Field.C (T a) where-  (/) = lift2 divide---sqrt :: (Field.C a) =>-   (a -> a) -> Core a -> Core a-sqrt fSqRt = lift1fromPowerSeries $ PS.sqrt (PS.const . (\[x] -> fSqRt x) . PS.coeffs)+   (/) = lift2 Core.divide   instance (Algebraic.C a) => Algebraic.C (T a) where-   sqrt   = lift1 (sqrt Algebraic.sqrt)---   x ^/ y = lift1 (pow (Algebraic.^/ y)+   sqrt   = lift1 (Core.sqrt Algebraic.sqrt)+--   x ^/ y = lift1 (Core.pow (Algebraic.^/ y) --                       (fromRational' y)) x---swapVariables :: Core a -> Core a-swapVariables = map reverse---differentiate0 :: (Ring.C a) => Core a -> Core a-differentiate0 =-   swapVariables . differentiate1 . swapVariables--differentiate1 :: (Ring.C a) => Core a -> Core a-differentiate1 = lift1fromPowerSeries $ map Differential.differentiate--integrate0 :: (Field.C a) => [a] -> Core a -> Core a-integrate0 cs =-   swapVariables . integrate1 cs . swapVariables--integrate1 :: (Field.C a) => [a] -> Core a -> Core a-integrate1 = zipWith PS.integrate-----{- |-Since the inner series must start with a zero,-the first term is omitted in y.--}-comp :: (Ring.C a) => [a] -> Core a -> Core a-comp = lift1fromPowerSeries . PS.comp . map PS.const
+ src/MathObj/PowerSeries2/Core.hs view
@@ -0,0 +1,89 @@+{-# LANGUAGE NoImplicitPrelude #-}+module MathObj.PowerSeries2.Core where++import qualified MathObj.PowerSeries as PS+import qualified MathObj.PowerSeries.Core as PSCore++import qualified Algebra.Differential   as Differential+import qualified Algebra.Vector         as Vector+import qualified Algebra.Field          as Field+import qualified Algebra.Ring           as Ring+import qualified Algebra.Additive       as Additive++import NumericPrelude.Base+-- import NumericPrelude.Numeric hiding (negate, sqrt, )+++type T a = [[a]]+++lift0fromPowerSeries :: [PS.T a] -> T a+lift0fromPowerSeries = map PS.coeffs++lift1fromPowerSeries ::+   ([PS.T a] -> [PS.T a]) -> (T a -> T a)+lift1fromPowerSeries f x0 =+   map PS.coeffs (f (map PS.fromCoeffs x0))++lift2fromPowerSeries ::+   ([PS.T a] -> [PS.T a] -> [PS.T a]) -> (T a -> T a -> T a)+lift2fromPowerSeries f x0 x1 =+   map PS.coeffs (f (map PS.fromCoeffs x0) (map PS.fromCoeffs x1))+++{- * Series arithmetic -}++add, sub :: (Additive.C a) => T a -> T a -> T a+add = PSCore.add+sub = PSCore.sub++negate :: (Additive.C a) => T a -> T a+negate = PSCore.negate+++scale :: Ring.C a => a -> T a -> T a+scale = map . (Vector.*>)++mul :: Ring.C a => T a -> T a -> T a+mul = lift2fromPowerSeries PSCore.mul+++divide :: (Field.C a) =>+   T a -> T a -> T a+divide = lift2fromPowerSeries PSCore.divide+++sqrt :: (Field.C a) =>+   (a -> a) -> T a -> T a+sqrt fSqRt =+   lift1fromPowerSeries $+   PSCore.sqrt (PS.const . (\[x] -> fSqRt x) . PS.coeffs)++++swapVariables :: T a -> T a+swapVariables = map reverse+++differentiate0 :: (Ring.C a) => T a -> T a+differentiate0 =+   swapVariables . differentiate1 . swapVariables++differentiate1 :: (Ring.C a) => T a -> T a+differentiate1 = lift1fromPowerSeries $ map Differential.differentiate++integrate0 :: (Field.C a) => [a] -> T a -> T a+integrate0 cs =+   swapVariables . integrate1 cs . swapVariables++integrate1 :: (Field.C a) => [a] -> T a -> T a+integrate1 = zipWith PSCore.integrate++++{- |+Since the inner series must start with a zero,+the first term is omitted in y.+-}+compose :: (Ring.C a) => [a] -> T a -> T a+compose = lift1fromPowerSeries . PSCore.compose . map PS.const
src/MathObj/PowerSum.hs view
@@ -16,8 +16,9 @@ -} module MathObj.PowerSum where -import qualified MathObj.Polynomial  as Poly-import qualified MathObj.PowerSeries as PS+import qualified MathObj.Polynomial as Poly+import qualified MathObj.Polynomial.Core as PolyCore+import qualified MathObj.PowerSeries.Core as PS  import qualified Algebra.VectorSpace  as VectorSpace import qualified Algebra.Module       as Module@@ -34,8 +35,8 @@ import qualified Data.List as List import Data.List.HT (shearTranspose, sieve) -import PreludeBase as P hiding (const)-import NumericPrelude as NP+import NumericPrelude.Base as P hiding (const)+import NumericPrelude.Numeric as NP   newtype T a = Cons {sums :: [a]}@@ -91,11 +92,11 @@ fromElemSym :: (Eq a, Ring.C a) => [a] -> [a] fromElemSym s =    fromIntegral (length s - 1) :-      Poly.alternate (divOneFlip s (Poly.differentiate s))+      PolyCore.alternate (divOneFlip s (PolyCore.differentiate s))  divOneFlip :: (Eq a, Ring.C a) => [a] -> [a] -> [a] divOneFlip (1:xs) =-   let aux (y:ys) = y : aux (ys - Poly.scale y xs)+   let aux (y:ys) = y : aux (ys - PolyCore.scale y xs)        aux [] = []    in  aux divOneFlip _ =@@ -104,19 +105,19 @@ fromElemSymDenormalized :: (Field.C a, ZeroTestable.C a) => [a] -> [a] fromElemSymDenormalized s =    fromIntegral (length s - 1) :-      Poly.alternate (PS.derivedLog s)+      PolyCore.alternate (PS.derivedLog s)   toElemSym :: (Field.C a, ZeroTestable.C a) => [a] -> [a] toElemSym p =-   let s' = Poly.mul (Poly.alternate (tail p)) s-       s  = Poly.integrate 1 s'+   let s' = PolyCore.mul (PolyCore.alternate (tail p)) s+       s  = PolyCore.integrate 1 s'    in  s  toElemSymInt :: (Integral.C a, ZeroTestable.C a) => [a] -> [a] toElemSymInt p =-   let s' = Poly.mul (Poly.alternate (tail p)) s-       s  = Poly.integrateInt 1 s'+   let s' = PolyCore.mul (PolyCore.alternate (tail p)) s+       s  = PolyCore.integrateInt 1 s'    in  s  @@ -125,11 +126,11 @@ fromPolynomial =    let aux s =           fromIntegral (length s - 1) :-             Poly.negate (PS.derivedLog s)+             PolyCore.negate (PS.derivedLog s)    in  aux . reverse . Poly.coeffs  elemSymFromPolynomial :: Additive.C a => Poly.T a -> [a]-elemSymFromPolynomial = Poly.alternate . reverse . Poly.coeffs+elemSymFromPolynomial = PolyCore.alternate . reverse . Poly.coeffs  {- toPolynomial is not possible because this had to consume the whole sum sequence. -} @@ -154,13 +155,13 @@ {- Use binomial expansion of (x+y)^n -} add :: (Ring.C a) => [a] -> [a] -> [a] add xs ys =-   let powers = shearTranspose (Poly.tensorProduct xs ys)+   let powers = shearTranspose (PolyCore.tensorProduct xs ys)    in  zipWith Ring.scalarProduct binomials powers  instance (Ring.C a) => Additive.C (T a) where    zero   = const zero    (+)    = lift2 add-   negate = lift1 Poly.alternate+   negate = lift1 PolyCore.alternate   {- * Ring -}@@ -230,4 +231,4 @@        yp = fromPolynomial (Poly.fromRoots ys)        ze = elemSymFromPolynomial (Poly.fromRoots zs)    in  zipWith (==) (toElemSym (powerOp xp yp)) ze-       -- Poly.equal (toElemSym (powerOp xp yp)) ze+       -- PolyCore.equal (toElemSym (powerOp xp yp)) ze
+ src/MathObj/RefinementMask2.hs view
@@ -0,0 +1,171 @@+{-# LANGUAGE NoImplicitPrelude #-}+module MathObj.RefinementMask2 (+   T, coeffs, fromCoeffs,+   fromPolynomial,+   toPolynomial,+   toPolynomialFast,+   refinePolynomial,+   ) where++import qualified MathObj.Polynomial as Poly+import qualified Algebra.RealField as RealField+import qualified Algebra.Field  as Field+import qualified Algebra.Ring   as Ring+import qualified Algebra.Vector as Vector++import qualified Data.List as List+import qualified Data.List.HT as ListHT+import qualified Data.List.Match as Match+import Control.Monad (liftM2, )++import qualified Test.QuickCheck as QC++import qualified NumericPrelude.List.Generic as NPList+import NumericPrelude.Base+import NumericPrelude.Numeric+++newtype T a = Cons {coeffs :: [a]}+++{-# INLINE fromCoeffs #-}+fromCoeffs :: [a] -> T a+fromCoeffs = lift0++{-# INLINE lift0 #-}+lift0 :: [a] -> T a+lift0 = Cons++{-+{-# INLINE lift1 #-}+lift1 :: ([a] -> [a]) -> (T a -> T a)+lift1 f (Cons x0) = Cons (f x0)++{-# INLINE lift2 #-}+lift2 :: ([a] -> [a] -> [a]) -> (T a -> T a -> T a)+lift2 f (Cons x0) (Cons x1) = Cons (f x0 x1)+-}++{-+Functor instance is e.g. useful for converting number types,+say 'Rational' to 'Double'.+-}++instance Functor T where+   fmap f (Cons xs) = Cons (map f xs)++{-# INLINE appPrec #-}+appPrec :: Int+appPrec  = 10++instance (Show a) => Show (T a) where+   showsPrec p (Cons xs) =+      showParen (p >= appPrec)+         (showString "RefinementMask2.fromCoeffs " . shows xs)++instance (QC.Arbitrary a, Field.C a) => QC.Arbitrary (T a) where+   arbitrary =+      liftM2+         (\degree body ->+            let s = sum body+            in  Cons $ map ((2 ^- degree - s) / NPList.lengthLeft body +) body)+         (QC.choose (-5,0)) QC.arbitrary+++{- |+Determine mask by Gauss elimination.++R - alternating binomial coefficients+L - differences of translated polynomials in columns++p2 = L * R^(-1) * m++R * L^(-1) * p2 = m+-}+fromPolynomial ::+   (Field.C a) => Poly.T a -> T a+fromPolynomial poly =+   fromCoeffs $+   foldr (\p ps ->+      ListHT.mapAdjacent (-) (p:ps++[0]))+      [] $+   foldr (\(db,dp) cs ->+      ListHT.switchR+         (error "RefinementMask2.fromPolynomial: polynomial should be non-empty")+         (\dps dpe ->+            cs ++ [(db - Ring.scalarProduct dps cs) / dpe])+         dp) [] $+   zip+      (Poly.coeffs $ Poly.dilate 2 poly)+      (List.transpose $+       Match.take (Poly.coeffs poly) $+       map Poly.coeffs $+       iterate polynomialDifference poly)++polynomialDifference ::+   (Ring.C a) => Poly.T a -> Poly.T a+polynomialDifference poly =+   Poly.fromCoeffs $ init $ Poly.coeffs $+   Poly.translate 1 poly - poly++{- |+If the mask does not sum up to a power of @1/2@+then the function returns 'Nothing'.+-}+toPolynomial ::+   (RealField.C a) => T a -> Maybe (Poly.T a)+toPolynomial (Cons []) = Just $ Poly.fromCoeffs []+toPolynomial mask =+   let s = sum $ coeffs mask+       ks = reverse $ takeWhile (<=1) $ iterate (2*) s+   in  case ks of+          1:ks0 ->+             Just $+             foldl+                (\p k ->+                   let ip = Poly.integrate zero p+                   in  ip + Poly.const (correctConstant (fmap (k/s*) mask) ip))+                (Poly.const 1) ks0+          _ -> Nothing+{-+> fmap (6 Vector.*>) $ toPolynomial (Cons [0.1, 0.02, 0.005::Rational])+Just (Polynomial.fromCoeffs [-12732 % 109375, 272 % 625, -18 % 25, 1 % 1])+-}++{-+The constant term must be zero,+higher terms must already satisfy the refinement constraint.+-}+correctConstant ::+   (Field.C a) => T a -> Poly.T a -> a+correctConstant mask poly =+   let refined = refinePolynomial mask poly+   in  head (Poly.coeffs refined) / (1 - sum (coeffs mask))++toPolynomialFast ::+   (RealField.C a) => T a -> Maybe (Poly.T a)+toPolynomialFast mask =+   let s = sum $ coeffs mask+       ks = reverse $ takeWhile (<=1) $ iterate (2*) s+   in  case ks of+          1:ks0 ->+             Just $+             foldl+                (\p k ->+                   let ip = Poly.integrate zero p+                       c = head (Poly.coeffs (refinePolynomial mask ip))+                   in  ip + Poly.const (c*k / ((1-k)*s)))+                (Poly.const 1) ks0+          _ -> Nothing++refinePolynomial ::+   (Ring.C a) => T a -> Poly.T a -> Poly.T a+refinePolynomial mask =+   Poly.shrink 2 .+   Vector.linearComb (coeffs mask) .+   iterate (Poly.translate 1)+{-+> mapM_ print $ take 50 $ iterate (refinePolynomial (Cons [0.1, 0.02, 0.005])) (Poly.fromCoeffs [0,0,0,1::Double])+...+Polynomial.fromCoeffs [-0.11640685714285712,0.4351999999999999,-0.7199999999999999,1.0]+-}
src/MathObj/RootSet.hs view
@@ -15,8 +15,9 @@ -} module MathObj.RootSet where -import qualified MathObj.Polynomial  as Poly-import qualified MathObj.PowerSum    as PowerSum+import qualified MathObj.Polynomial      as Poly+import qualified MathObj.Polynomial.Core as PolyCore+import qualified MathObj.PowerSum        as PowerSum  import qualified Algebra.Algebraic    as Algebraic import qualified Algebra.IntegralDomain as Integral@@ -28,8 +29,8 @@ import qualified Data.List.Match as Match import Control.Monad (liftM2) -import PreludeBase as P hiding (const)-import NumericPrelude as NP+import NumericPrelude.Base as P hiding (const)+import NumericPrelude.Numeric as NP   newtype T a = Cons {coeffs :: [a]}@@ -51,11 +52,11 @@ const x = Cons [1,x]  -toPolynomial :: Poly.T a -> T a-toPolynomial xs = Cons (reverse (Poly.coeffs xs))+toPolynomial :: T a -> Poly.T a+toPolynomial (Cons xs) = Poly.fromCoeffs (reverse xs) -fromPolynomial :: T a -> Poly.T a-fromPolynomial (Cons xs) = Poly.fromCoeffs (reverse xs)+fromPolynomial :: Poly.T a -> T a+fromPolynomial xs = Cons (reverse (Poly.coeffs xs))   @@ -69,7 +70,7 @@ {- | cf. 'MathObj.Polynomial.mulLinearFactor' -} addRoot :: Ring.C a => a -> [a] -> [a] addRoot x yt@(y:ys) =-   y : (ys + Poly.scale x yt)+   y : (ys + PolyCore.scale x yt) addRoot _ [] =    error "addRoot: list of elementar symmetric terms must consist at least of a 1" @@ -132,7 +133,7 @@ instance (Field.C a, ZeroTestable.C a) => Additive.C (T a) where    zero   = const zero    (+)    = lift2 add-   negate = lift1 Poly.alternate+   negate = lift1 PolyCore.alternate   {- * Ring -}
− src/MyPrelude.hs 
@@ -1,5 +0,0 @@-{-# LANGUAGE NoImplicitPrelude #-}-module MyPrelude(module NumericPrelude, module PreludeBase, max, min, abs) where-import NumericPrelude hiding (abs)-import PreludeBase hiding (max,min)-import Algebra.Lattice (max,min,abs)
src/Number/Complex.hs view
@@ -32,6 +32,7 @@         fromPolar,         cis,         signum,+        signumNorm,         toPolar,         magnitude,         magnitudeSqr,@@ -63,7 +64,8 @@ import qualified Algebra.Units              as Units import qualified Algebra.PrincipalIdealDomain as PID import qualified Algebra.IntegralDomain     as Integral-import qualified Algebra.Real               as Real+import qualified Algebra.RealRing          as RealRing+import qualified Algebra.Absolute               as Absolute import qualified Algebra.Ring               as Ring import qualified Algebra.Additive           as Additive import qualified Algebra.ZeroTestable       as ZeroTestable@@ -84,8 +86,8 @@ import Control.Monad (liftM2, )  import qualified Prelude as P-import PreludeBase-import NumericPrelude hiding (signum, exp, )+import NumericPrelude.Base+import NumericPrelude.Numeric hiding (signum, exp, ) import Text.Show.HT (showsInfixPrec, ) import Text.Read.HT (readsInfixPrec, ) @@ -181,16 +183,25 @@  {- | Scale a complex number to magnitude 1. -For a complex number @z@, @'abs' z@ is a number with the magnitude of @z@,-but oriented in the positive real direction, whereas @'signum' z@-has the phase of @z@, but unit magnitude.+For a complex number @z@,+@'abs' z@ is a number with the magnitude of @z@,+but oriented in the positive real direction,+whereas @'signum' z@ has the phase of @z@, but unit magnitude. -}+ {- SPECIALISE signum :: T Double -> T Double -}-{-# INLINE signum #-}-signum :: (Algebraic.C a, NormedEuc.C a a, ZeroTestable.C a) => T a -> T a+signum :: (Algebraic.C a, ZeroTestable.C a) => T a -> T a signum z =    if isZero z      then zero+     else scale (recip (magnitude z)) z++{- SPECIALISE signumNorm :: T Double -> T Double -}+{-# INLINE signumNorm #-}+signumNorm :: (Algebraic.C a, NormedEuc.C a a, ZeroTestable.C a) => T a -> T a+signumNorm z =+   if isZero z+     then zero      else scale (recip (NormedEuc.norm z)) z  -- | Form a complex number from polar components of magnitude and phase.@@ -305,6 +316,14 @@     {-# INLINE fromInteger #-}     fromInteger                 =  fromReal . fromInteger +instance  (Absolute.C a, Algebraic.C a) => Absolute.C (T a)  where+    {- SPECIALISE instance Absolute.C (T Float) -}+    {- SPECIALISE instance Absolute.C (T Double) -}+    {-# INLINE abs #-}+    {-# INLINE signum #-}+    abs x  = Cons (magnitude x) zero+    signum = signum+ instance Vector.C T where    {-# INLINE zero #-}    zero  = zero@@ -474,7 +493,7 @@     power  =  defltPow  -instance  (Real.C a, Algebraic.C a, Power a) =>+instance  (RealRing.C a, Algebraic.C a, Power a) =>           Algebraic.C (T a)  where     {-# INLINE sqrt #-}     sqrt z@(Cons x y)  =  if z == zero@@ -488,7 +507,7 @@     (^/) = flip power  -instance  (Real.C a, RealTrans.C a, Power a) =>+instance  (RealRing.C a, RealTrans.C a, Power a) =>           Trans.C (T a)  where     {- SPECIALISE instance Trans.C (T Float) -}     {- SPECIALISE instance Trans.C (T Double) -}
src/Number/DimensionTerm.hs view
@@ -20,7 +20,7 @@ import qualified Algebra.Module        as Module import qualified Algebra.Algebraic     as Algebraic import qualified Algebra.Field         as Field-import qualified Algebra.Real          as Real+import qualified Algebra.Absolute          as Absolute import qualified Algebra.Ring          as Ring import qualified Algebra.Additive      as Additive @@ -32,7 +32,7 @@ import System.Random (Random, randomR, random)  import Data.Tuple.HT (mapFst, )-import PreludeBase+import NumericPrelude.Base import Prelude ()  @@ -139,11 +139,11 @@ sqrt (Cons x) = Cons (Algebraic.sqrt x)  -abs :: (Dim.C u, Real.C a) => T u a -> T u a-abs (Cons x) = Cons (Real.abs x)+abs :: (Dim.C u, Absolute.C a) => T u a -> T u a+abs (Cons x) = Cons (Absolute.abs x) -absSignum :: (Dim.C u, Real.C a) => T u a -> (T u a, a)-absSignum x0@(Cons x) = (abs x0, Real.signum x)+absSignum :: (Dim.C u, Absolute.C a) => T u a -> (T u a, a)+absSignum x0@(Cons x) = (abs x0, Absolute.signum x)  scale, (*&) :: (Dim.C u, Ring.C a) =>    a -> T u a -> T u a
src/Number/DimensionTerm/SI.hs view
@@ -45,8 +45,8 @@ import qualified Number.DimensionTerm  as DN import qualified Number.SI.Unit as SI --- aimport PreludeBase hiding (length)-import NumericPrelude hiding (one)+-- aimport NumericPrelude.Base hiding (length)+import NumericPrelude.Numeric hiding (one)   second  :: Field.C a => DN.Time        a
src/Number/FixedPoint.hs view
@@ -16,7 +16,7 @@ -} module Number.FixedPoint where -import qualified Algebra.RealField    as RealField+import qualified Algebra.RealRing    as RealRing import qualified Algebra.Additive       as Additive -- import qualified Algebra.ZeroTestable   as ZeroTestable import qualified Algebra.Transcendental as Trans@@ -29,18 +29,18 @@ import Data.List (transpose, unfoldr, ) import Data.Char (intToDigit, ) -import PreludeBase-import NumericPrelude hiding (recip, sqrt, exp, sin, cos, tan,+import NumericPrelude.Base+import NumericPrelude.Numeric hiding (recip, sqrt, exp, sin, cos, tan,                               fromRational') -import qualified NumericPrelude as NP+import qualified NumericPrelude.Numeric as NP   {- ** Conversion -}  {- ** other number types -} -fromFloat :: RealField.C a => Integer -> a -> Integer+fromFloat :: RealRing.C a => Integer -> a -> Integer fromFloat den x =    round (x * NP.fromInteger den)
src/Number/FixedPoint/Check.hs view
@@ -7,18 +7,18 @@  import qualified Algebra.Transcendental as Trans import qualified Algebra.Algebraic      as Algebraic-import qualified Algebra.RealField      as RealField+import qualified Algebra.RealRing      as RealRing import qualified Algebra.Field          as Field-import qualified Algebra.Real           as Real+import qualified Algebra.Absolute           as Absolute import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive import qualified Algebra.ZeroTestable   as ZeroTestable -import PreludeBase-import NumericPrelude   hiding (fromRational')+import NumericPrelude.Base+import NumericPrelude.Numeric   hiding (fromRational')  import qualified Prelude        as P98-import qualified NumericPrelude as NP+import qualified NumericPrelude.Numeric as NP   {- * Types -}@@ -33,7 +33,7 @@  {- ** other number types -} -fromFloat :: RealField.C a => Integer -> a -> T+fromFloat :: RealRing.C a => Integer -> a -> T fromFloat den x =    cons den (FP.fromFloat den x) @@ -45,7 +45,7 @@ fromRational' den x =    cons den (round (x * NP.fromInteger den)) -fromFloatBasis :: RealField.C a => Integer -> Int -> a -> T+fromFloatBasis :: RealRing.C a => Integer -> Int -> a -> T fromFloatBasis basis numDigits =    fromFloat (ringPower numDigits basis) @@ -165,11 +165,11 @@    compare (Cons xd xn) (Cons yd yn) =       commonDenominator xd yd (compare xn yn) -instance Real.C T where+instance Absolute.C T where    abs = lift1 (const abs)-   -- use default implementation for signum+   signum = Absolute.signumOrd -instance RealField.C T where+instance RealRing.C T where    splitFraction (Cons xd xn) =       let (int, frac) = divMod xd xn       in  (fromInteger int, Cons xd frac)
src/Number/GaloisField2p32m5.hs view
@@ -25,8 +25,8 @@  import Test.QuickCheck (Arbitrary(arbitrary), ) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric   newtype T = Cons {decons :: Word32}
src/Number/NonNegative.hs view
@@ -1,4 +1,5 @@-{-# OPTIONS -XNoImplicitPrelude -fno-warn-orphans #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}  {- Rationale for -fno-warn-orphans:@@ -33,13 +34,14 @@ import qualified Algebra.NonNegative        as NonNeg import qualified Algebra.Transcendental     as Trans import qualified Algebra.Algebraic          as Algebraic-import qualified Algebra.RealField          as RealField+import qualified Algebra.RealRing          as RealRing import qualified Algebra.Field              as Field import qualified Algebra.RealIntegral       as RealIntegral import qualified Algebra.IntegralDomain     as Integral-import qualified Algebra.Real               as Real+import qualified Algebra.Absolute               as Absolute import qualified Algebra.Ring               as Ring import qualified Algebra.Additive           as Additive+import qualified Algebra.Monoid             as Monoid import qualified Algebra.ZeroTestable       as ZeroTestable  import qualified Algebra.ToInteger          as ToInteger@@ -50,9 +52,9 @@  import qualified Prelude as P -import PreludeBase+import NumericPrelude.Base import Data.Tuple.HT (mapSnd, mapPair, )-import NumericPrelude hiding (Int, Integer, Float, Double, Rational, )+import NumericPrelude.Numeric hiding (Int, Integer, Float, Double, Rational, )   {- |@@ -114,8 +116,13 @@ instance ZeroTestable.C a => ZeroTestable.C (T a) where    isZero = isZero . toNumber +instance (Additive.C a) => Monoid.C (T a) where+   idt = fromNumberUnsafe Additive.zero+   x <*> y = fromNumberUnsafe (toNumber x + toNumber y)+--   mconcat = fromNumberUnsafe . sum . map toNumber+ instance (Ord a, Additive.C a) => NonNeg.C (T a) where-   x -| y = fromNumberClip (toNumber x - toNumber y)+   split = NonNeg.splitDefault toNumber fromNumberUnsafe  instance (Ord a, Additive.C a) => Additive.C (T a) where    zero   = fromNumberUnsafe zero@@ -127,7 +134,7 @@    (*)    = lift2 (*)    fromInteger = fromNumberWrap "fromInteger" . fromInteger -instance ToRational.C a => ToRational.C (T a) where+instance (Ord a, ToRational.C a) => ToRational.C (T a) where    toRational = ToRational.toRational . toNumber  instance ToInteger.C a => ToInteger.C (T a) where@@ -167,11 +174,11 @@    (/) = lift2 (/)  -instance (ZeroTestable.C a, Real.C a) => Real.C (T a) where+instance (ZeroTestable.C a, Ord a, Absolute.C a) => Absolute.C (T a) where    abs    = lift abs    signum = lift signum -instance (RealField.C a) => RealField.C (T a) where+instance (RealRing.C a) => RealRing.C (T a) where    splitFraction = mapSnd fromNumberUnsafe . splitFraction . toNumber    truncate = truncate . toNumber    round    = round    . toNumber
src/Number/NonNegativeChunky.hs view
@@ -1,5 +1,5 @@ {- |-Copyright   :  (c) Henning Thielemann 2007+Copyright   :  (c) Henning Thielemann 2007-2010  Maintainer  :  haskell@henning-thielemann.de Stability   :  stable@@ -17,14 +17,15 @@ -} module Number.NonNegativeChunky    (T, fromChunks, toChunks, fromNumber, toNumber, fromChunky98, toChunky98,-    minMaxDiff, normalize, isNull, isPositive) where+    minMaxDiff, normalize, isNull, isPositive,+    divModLazy, divModStrict, ) where -import qualified Numeric.NonNegative.ChunkyPrivate as Chunky98+import qualified Numeric.NonNegative.Chunky as Chunky98 import qualified Numeric.NonNegative.Class as NonNeg98  import qualified Algebra.NonNegative  as NonNeg import qualified Algebra.Field        as Field-import qualified Algebra.Real         as Real+import qualified Algebra.Absolute         as Absolute import qualified Algebra.Ring         as Ring import qualified Algebra.Additive     as Additive import qualified Algebra.ToInteger    as ToInteger@@ -38,13 +39,13 @@ import qualified Data.Monoid as Mn98  import Control.Monad (liftM, liftM2, )+import Data.Tuple.HT (mapFst, mapSnd, mapPair, )  import Test.QuickCheck (Arbitrary(arbitrary)) -import NumericPrelude-import Data.Tuple.HT (mapFst, mapPair, )-import PreludeBase-import qualified Prelude     as P98+import NumericPrelude.Numeric+import NumericPrelude.Base+import qualified Prelude as P98 (Num(..), Fractional(..), )   {- |@@ -71,7 +72,7 @@ toChunks = decons  fromChunky98 :: (NonNeg.C a, NonNeg98.C a) => Chunky98.T a -> T a-fromChunky98 = fromChunks . Chunky98.toChunksUnsafe+fromChunky98 = fromChunks . Chunky98.toChunks  toChunky98 :: (NonNeg.C a, NonNeg98.C a) => T a -> Chunky98.T a toChunky98 = Chunky98.fromChunks . toChunks@@ -80,7 +81,7 @@ fromNumber = fromChunks . (:[])  toNumber :: NonNeg.C a => T a -> a-toNumber =  sum . toChunks+toNumber =  Monoid.cumulate . toChunks   @@ -92,10 +93,10 @@ Remove zero chunks. -} normalize :: NonNeg.C a => T a -> T a-normalize = fromChunks . filter (>zero) . toChunks+normalize = fromChunks . filter (> NonNeg.zero) . toChunks  isNullList :: NonNeg.C a => [a] -> Bool-isNullList = null . filter (>zero)+isNullList = null . filter (> NonNeg.zero)  isNull :: NonNeg.C a => T a -> Bool isNull = isNullList . toChunks@@ -130,48 +131,44 @@      else error ("Numeric.NonNegative.Chunky."++funcName++": negative number")  -glue :: (NonNeg.C a) => [a] -> [a] -> ([a], [a], Bool)-glue [] ys = ([], ys, True)-glue xs [] = ([], xs, False)+glue :: (NonNeg.C a) => [a] -> [a] -> ([a], (Bool, [a]))+glue [] ys = ([], (True,  ys))+glue xs [] = ([], (False, xs)) glue (x:xs) (y:ys) =-   let (z,(zs,rs,b)) =-           case compare x y of-              LT -> (x, glue xs ((y-x):ys))-              GT -> (y, glue ((x-y):xs) ys)-              EQ -> (x, glue xs ys)-   in  (z:zs,rs,b)+   let (z,~(zs,brs)) =+          flip mapSnd (NonNeg.split x y) $+          \(b,d) ->+             if b+               then glue xs $+                    if NonNeg.zero == d+                      then ys else d:ys+               else glue (d:xs) ys+   in  (z:zs,brs) -{- |-In @minMaxDiff x y == (z,r,b)@-@z@ represents @min x y@,-@r@ represents @max x y - min x y@,-and @x<y  ==>  b@ or @x>y  ==>  not b@,- for @x==y@ the value of b is arbitrary.--}-minMaxDiff :: (NonNeg.C a) => T a -> T a -> (T a, T a, Bool)+minMaxDiff :: (NonNeg.C a) => T a -> T a -> (T a, (Bool, T a)) minMaxDiff (Cons xs) (Cons ys) =-   let (zs, rs, b) = glue xs ys-   in  (Cons zs, Cons rs, b)+   let (zs, (b, rs)) = glue xs ys+   in  (Cons zs, (b, Cons rs))  equalList :: (NonNeg.C a) => [a] -> [a] -> Bool equalList x y =-   let (_,r,_) = glue x y-   in  isNullList r+   isNullList $ snd $ snd $ glue x y  compareList :: (NonNeg.C a) => [a] -> [a] -> Ordering compareList x y =-   let (_,r,b) = glue x y+   let (b,r) = snd $ glue x y    in  if isNullList r          then EQ          else if b then LT else GT  minList :: (NonNeg.C a) => [a] -> [a] -> [a] minList x y =-   let (z,_,_) = glue x y in z+   fst $ glue x y  maxList :: (NonNeg.C a) => [a] -> [a] -> [a] maxList x y =-   let (z,r,_) = glue x y in z++r+   -- matching the inner pair lazily is important+   let (z,~(_,r)) = glue x y in z++r   instance (NonNeg.C a) => Eq (T a) where@@ -184,12 +181,9 @@   instance (NonNeg.C a) => NonNeg.C (T a) where-   (-|) =-      lift2 (\x w ->-         let sub _ [] = []-             sub z (y:ys) =-                if z<y then (y-z):ys else sub (z-y) ys-         in  foldr sub x w)+   split (Cons xs) (Cons ys) =+      let (zs, ~(b, rs)) = glue xs ys+      in  (Cons zs, (b, Cons rs))  instance (ZeroTestable.C a) => ZeroTestable.C (T a) where    isZero = isNullZT@@ -197,9 +191,9 @@ instance (NonNeg.C a) => Additive.C (T a) where    zero  = Monoid.idt    (+)   = (Monoid.<*>)-   x - y =-      let (_,d,b) = glue (toChunks x) (toChunks y)-          d' = fromChunks d+   (Cons x) - (Cons y) =+      let (b,d) = snd $ glue x y+          d' = Cons d       in check "-" (not b || isNull d') d'    negate x = check "negate" (isNull x) x {-@@ -213,7 +207,7 @@    (*)   = lift2 (liftM2 (*))    fromInteger = fromNumber . fromInteger -instance (Ring.C a, ZeroTestable.C a, NonNeg.C a) => Real.C (T a) where+instance (Ring.C a, ZeroTestable.C a, NonNeg.C a) => Absolute.C (T a) where    abs    = id    signum = fromNumber . (\b -> if b then one else zero) . isPositive @@ -229,18 +223,53 @@    rem  = mod    quotRem = divMod +{- |+'divMod' is implemented in terms of 'divModStrict'.+If it is needed we could also provide a function+that accesses the divisor first in a lazy way+and then uses a strict divisor for subsequent rounds of the subtraction loop.+This way we can handle the cases \"dividend smaller than divisor\"+and \"dividend greater than divisor\" in a lazy and efficient way.+However changing the way of operation within one number is also not nice.+-} instance (Ord a, Integral.C a, NonNeg.C a) => Integral.C (T a) where-   divMod x0 y0 =-      let y = toChunks y0-          recurse x =-             let (r,d,b) = glue y x-             in  if not b-                   then ([], r)-                   else mapFst (one:) (recurse d)-      in  mapPair-             (fromChunks, fromChunks)-             (recurse (toChunks x0))+   divMod x y =+      mapSnd fromNumber $+      divModStrict x (toNumber y) +{- |+divModLazy accesses the divisor in a lazy way.+However this is only relevant if the dividend is smaller than the divisor.+For large dividends the divisor will be accessed multiple times+but since it is already fully evaluated it could also be strict.+-}+divModLazy ::+   (Ring.C a, NonNeg.C a) =>+   T a -> T a -> (T a, T a)+divModLazy x0 y0 =+   let y = toChunks y0+       recourse x =+          let (r,~(b,d)) = glue y x+          in  if not b+                then ([], r)+                else mapFst (one:) (recourse d)+   in  mapPair+          (fromChunks, fromChunks)+          (recourse (toChunks x0))++{- |+This function has a strict divisor+and maintains the chunk structure of the dividend at a smaller scale.+-}+divModStrict ::+   (Integral.C a, NonNeg.C a) =>+   T a -> a -> (T a, a)+divModStrict x0 y =+   let recourse [] r = ([], r)+       recourse (x:xs) r0 =+          case divMod (x+r0) y of+             (q,r1) -> mapFst (q:) $ recourse xs r1+   in  mapFst fromChunks $ recourse (toChunks x0) zero
src/Number/OccasionallyScalarExpression.hs view
@@ -19,7 +19,7 @@ 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.Absolute                as Absolute import qualified Algebra.Ring                as Ring import qualified Algebra.Additive            as Additive import qualified Algebra.ZeroTestable        as ZeroTestable@@ -30,8 +30,8 @@ import Data.Maybe(fromMaybe) import Data.Array(listArray,(!)) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric   {- | A value of type 'T' stores information on how to resolve unit violations.@@ -123,7 +123,7 @@ instance (ZeroTestable.C v) => ZeroTestable.C (T a v) where   isZero (Cons _ x) = isZero x -instance (Real.C v) => Real.C (T a v) where+instance (Absolute.C v) => Absolute.C (T a v) where   {- are these definitions sensible? -}   abs    = lift abs   signum = lift signum
src/Number/PartiallyTranscendental.hs view
@@ -17,8 +17,8 @@ import qualified Algebra.Additive       as Additive -- import qualified Algebra.ZeroTestable   as ZeroTestable -import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base  import qualified Prelude as P
src/Number/Peano.hs view
@@ -15,26 +15,31 @@ import qualified Algebra.Units                as Units import qualified Algebra.RealIntegral         as RealIntegral import qualified Algebra.IntegralDomain       as Integral-import qualified Algebra.Real                 as Real+import qualified Algebra.Absolute             as Absolute import qualified Algebra.Ring                 as Ring import qualified Algebra.Additive             as Additive import qualified Algebra.ZeroTestable         as ZeroTestable import qualified Algebra.Indexable            as Indexable+import qualified Algebra.Monoid               as Monoid  import qualified Algebra.ToInteger            as ToInteger import qualified Algebra.ToRational           as ToRational import qualified Algebra.NonNegative          as NonNeg +import qualified Algebra.EqualityDecision as EqDec+import qualified Algebra.OrderDecision    as OrdDec+ import Data.Maybe (catMaybes, ) import Data.Array(Ix(..))  import qualified Prelude     as P98-import qualified PreludeBase as P-import qualified NumericPrelude as NP+import qualified NumericPrelude.Base as P+import qualified NumericPrelude.Numeric as NP import Data.List.HT (mapAdjacent, shearTranspose, )+import Data.Tuple.HT (mapFst, ) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric   data T = Zero@@ -119,6 +124,26 @@    enumFromThenTo =    -} ++{- |+If all values are completely defined,+then it holds++> if b then x else y == ifLazy b x y++However if @b@ is undefined,+then it is at least known that the result is larger than @min x y@.+-}+ifLazy :: Bool -> T -> T -> T+ifLazy b (Succ x) (Succ y) = Succ (ifLazy b x y)+ifLazy b x y = if b then x else y++instance EqDec.C T where+   (==?) x y = ifLazy (x==y)++instance OrdDec.C T where+   (<=?) x y le gt = ifLazy (x<=y) le gt+ {- The default instance is good for compare, but fails for min and max.@@ -158,12 +183,12 @@ -} argMinFull :: (T,a) -> (T,a) -> (T,a) argMinFull (x0,xv) (y0,yv) =-   let recurse (Succ x) (Succ y) =-          let (z,zv) = recurse x y+   let recourse (Succ x) (Succ y) =+          let (z,zv) = recourse x y           in  (Succ z, zv)-       recurse Zero _ = (Zero,xv)-       recurse _    _ = (Zero,yv)-   in  recurse x0 y0+       recourse Zero _ = (Zero,xv)+       recourse _    _ = (Zero,yv)+   in  recourse x0 y0  {- | On equality the first operand is returned.@@ -177,12 +202,12 @@  argMaxFull :: (T,a) -> (T,a) -> (T,a) argMaxFull (x0,xv) (y0,yv) =-   let recurse (Succ x) (Succ y) =-          let (z,zv) = recurse x y+   let recourse (Succ x) (Succ y) =+          let (z,zv) = recourse x y           in  (Succ z, zv)-       recurse x Zero = (x,xv)-       recurse _ y    = (y,yv)-   in  recurse x0 y0+       recourse x Zero = (x,xv)+       recourse _ y    = (y,yv)+   in  recourse x0 y0  {- | On equality the first operand is returned.@@ -221,24 +246,23 @@  toListMaybe :: a -> T -> [Maybe a] toListMaybe a =-   let recurse Zero     = [Just a]-       recurse (Succ x) = Nothing : recurse x-   in  recurse+   let recourse Zero     = [Just a]+       recourse (Succ x) = Nothing : recourse x+   in  recourse  {- |-In @glue x y == (z,r,b)@+In @glue x y == (z,(b,r))@ @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 :: T -> T -> (T, (Bool, T))+glue Zero ys = (Zero, (True, ys))+glue xs Zero = (Zero, (False, xs)) glue (Succ xs) (Succ ys) =-   let (common, difference, sign) = glue xs ys-   in  (Succ common, difference, sign)+   mapFst Succ $ glue xs ys  {- Implementation notes:@@ -256,7 +280,7 @@    and . catMaybes . concat .    shearTranspose .    mapAdjacent (\x y ->-      let (costs0,_,le) = glue x y+      let (costs0,(le,_)) = glue x y       in  toListMaybe le costs0)  @@ -277,7 +301,7 @@ infixr 3 &&~ (&&~) :: Valuable Bool -> Valuable Bool -> Valuable Bool (&&~) (Valuable xc xb) (Valuable yc yb) =-   let (minc,difc,le) = glue xc yc+   let (minc,~(le,difc)) = glue xc yc        (bCheap,bExpensive) =           if le             then (xb,yb)@@ -296,7 +320,7 @@  leW :: T -> T -> Valuable Bool leW x y =-   let (minc,_difc,le) = glue x y+   let (minc,~(le,_difc)) = glue x y    in  Valuable minc le  isAscendingW :: [T] -> Valuable Bool@@ -313,7 +337,7 @@  -- instances -instance Real.C T where+instance Absolute.C T where    signum Zero     = zero    signum (Succ _) = one    abs             = id@@ -338,12 +362,13 @@             then (zero,x)             else let (q,r) = divMod d y in (succ q,r) +instance Monoid.C T where+   idt = zero+   (<*>) = add+   cumulate = foldr add Zero+ instance NonNeg.C T where-   (-|) x y =-      let (isNeg,d) = subNeg y x-      in  if isNeg-            then zero-            else d+   split = glue  instance Ix T where    range = uncurry enumFromTo
src/Number/Physical.hs view
@@ -23,7 +23,7 @@ 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.Absolute                as Absolute import qualified Algebra.Ring                as Ring import qualified Algebra.Additive            as Additive import qualified Algebra.ZeroTestable        as ZeroTestable@@ -39,8 +39,8 @@ import Data.Maybe.HT(toMaybe) import Data.Maybe(fromMaybe) -import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base   -- | A Physics.Quantity.Value.T combines a numeric value with a physical unit.@@ -167,7 +167,7 @@   I need absolute values for sample rates and amplitudes.   There the second interpretation is needed. -}-instance (Ord i, Real.C a) => Real.C (T i a) where+instance (Ord i, Absolute.C a) => Absolute.C (T i a) where   abs               = lift abs   signum (Cons _ x) = fromScalarSingle (signum x) @@ -179,8 +179,7 @@ instance (Ord i, Algebraic.C a) => Algebraic.C (T i a) where   sqrt (Cons xu x) = Cons (Unit.ratScale 0.5 xu) (sqrt x)   Cons xu x ^/ y =-     let y' = fromRational' (toRational y)-     in  Cons (Unit.ratScale y' xu) (x ^/ y)+     Cons (Unit.ratScale (fromRational' y) xu) (x ^/ y)  instance (Ord i, Trans.C a) => Trans.C (T i a) where   pi      = fromScalarSingle pi
src/Number/Physical/Read.hs view
@@ -22,8 +22,8 @@ import Text.ParserCombinators.Parsec import Control.Monad(liftM) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  mulPrec :: Int mulPrec = 7
src/Number/Physical/Show.hs view
@@ -24,8 +24,8 @@ import Data.List(find) import Data.Maybe(mapMaybe) -import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base   mulPrec :: Int
src/Number/Physical/Unit.hs view
@@ -21,8 +21,8 @@  import Data.Maybe.HT(toMaybe) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  {- | A Unit.T is a sparse vector with integer entries    Each map n->m means that the unit of the n-th dimension
src/Number/Physical/UnitDatabase.hs view
@@ -22,8 +22,8 @@ import Data.Maybe.HT (toMaybe) import Data.List (findIndices, partition, unfoldr, find, minimumBy) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  type T i a = [UnitSet i a]
src/Number/Positional.hs view
@@ -14,7 +14,7 @@ module Number.Positional where  import qualified MathObj.LaurentPolynomial as LPoly-import qualified MathObj.Polynomial        as Poly+import qualified MathObj.Polynomial.Core   as Poly  import qualified Algebra.IntegralDomain as Integral import qualified Algebra.Ring           as Ring@@ -22,11 +22,11 @@ import qualified Algebra.ToInteger      as ToInteger  import qualified Prelude as P98-import qualified PreludeBase as P-import qualified NumericPrelude as NP+import qualified NumericPrelude.Base as P+import qualified NumericPrelude.Numeric as NP -import PreludeBase-import NumericPrelude hiding (sqrt, tan, one, zero, )+import NumericPrelude.Base+import NumericPrelude.Numeric hiding (sqrt, tan, one, zero, )  import qualified Data.List as List import Data.Char (intToDigit)@@ -43,7 +43,7 @@   {--bugs:+FIXME:  defltBase = 10 defltExp = 4@@ -65,7 +65,7 @@  moveToZero :: Basis -> Digit -> (Digit,Digit) moveToZero b n =-   let b2 = NP.div b 2+   let b2 = div b 2        (q,r) = divMod (n+b2) b    in  (q,r-b2) @@ -492,6 +492,79 @@ equalApprox b bnd x y =    fst (trimUntil bnd (sub b x y)) == bnd ++{- |+If all values are completely defined,+then it holds++> if b then x else y == ifLazy b x y++However if @b@ is undefined,+then it is at least known that the result is between @x@ and @y@.+-}+ifLazy :: Basis -> Bool -> T -> T -> T+ifLazy b c x@(xe, _) y@(ye, _) =+   let ze = max xe ye+       xm = alignMant b ze x+       ym = alignMant b ze y+       recurse :: Mantissa -> Mantissa -> Mantissa+       recurse xs0 ys0 =+          withTwoMantissas xs0 ys0 [] $ \(x0,xs1) (y0,ys1) ->+          if abs (y0-x0) > 2+            then if c then xs0 else ys0+            else+              {-+              @x0==y0 || c@ means that in case of @x0==y0@+              we do not have to check @c@.+              -}+              withTwoMantissas xs1 ys1 ((if x0==y0 || c then x0 else y0) : []) $+                  \(x1,xs2) (y1,ys2) ->+                {-+                We can choose @z0@ only when knowing also x1 and y1.+                Because of x0x1 = 09 and y0y1 = 19+                we may always choose the larger one of x0 and y0+                in order to get admissible digit z1.+                But this would be wrong for x0x1 = 0(-9) and y0y1 = 1(-9).+                -}+                let z0  = mean2 b (x0,x1) (y0,y1)+                    x1' = x1+(x0-z0)*b+                    y1' = y1+(y0-z0)*b+                in  if abs x1' < b  &&  abs y1' < b+                      then z0 : recurse (x1':xs2) (y1':ys2)+                      else if c then xs0 else ys0+   in  (ze, recurse xm ym)++{- |+> mean2 b (x0,x1) (y0,y1)++computes @ round ((x0.x1 + y0.y1)/2) @,+where @x0.x1@ and @y0.y1@ are positional rational numbers+with respect to basis @b@+-}+{-# INLINE mean2 #-}+mean2 :: Basis -> (Digit,Digit) -> (Digit,Digit) -> Digit+mean2 b (x0,x1) (y0,y1) =+   ((x0+y0+1)*b + (x1+y1)) `div` (2*b)++{-+In a first trial I used++> zipMantissas :: Mantissa -> Mantissa -> [(Digit, Digit)]++for implementation of ifLazy.+However, this required to extract digits from the pairs+after the decision for an argument.+With withTwoMantissas we can just return a pointer to the original list.+-}+withTwoMantissas ::+   Mantissa -> Mantissa ->+   a ->+   ((Digit,Mantissa) -> (Digit,Mantissa) -> a) ->+   a+withTwoMantissas [] [] r _ = r+withTwoMantissas [] (y:ys) _ f = f (0,[]) (y,ys)+withTwoMantissas (x:xs) [] _ f = f (x,xs) (0,[])+withTwoMantissas (x:xs) (y:ys) _ f = f (x,xs) (y,ys)   align :: Basis -> Exponent -> T -> T
src/Number/Positional/Check.hs view
@@ -21,16 +21,20 @@ import qualified Algebra.Algebraic          as Algebraic import qualified Algebra.RealField          as RealField import qualified Algebra.Field              as Field-import qualified Algebra.Real               as Real+import qualified Algebra.RealRing           as RealRing+import qualified Algebra.Absolute           as Absolute import qualified Algebra.Ring               as Ring import qualified Algebra.Additive           as Additive import qualified Algebra.ZeroTestable       as ZeroTestable -import qualified PreludeBase as P+import qualified Algebra.EqualityDecision as EqDec+import qualified Algebra.OrderDecision    as OrdDec++import qualified NumericPrelude.Base as P import qualified Prelude     as P98 -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   {- |@@ -96,9 +100,16 @@  lift2 :: (Int -> Pos.T -> Pos.T -> Pos.T) -> T -> T -> T lift2 op (Cons xb xe xm) (Cons yb ye ym) =-   let zb = commonBasis xb yb-   in  uncurry (Cons zb) (op xb (xe, xm) (ye, ym))+   let b = commonBasis xb yb+   in  uncurry (Cons b) (op b (xe, xm) (ye, ym)) +{-+lift4 :: (Int -> Pos.T -> Pos.T -> Pos.T -> Pos.T -> Pos.T) -> T -> T -> T -> T -> T+lift4 op (Cons xb xe xm) (Cons yb ye ym) (Cons zb ze zm) (Cons wb we wm) =+   let b = xb `commonBasis` yb `commonBasis` zb `commonBasis` wb+   in  uncurry (Cons b) (op b (xe, xm) (ye, ym) (ze, zm) (we, wm))+-}+ commonBasis :: Pos.Basis -> Pos.Basis -> Pos.Basis commonBasis xb yb =    if xb == yb@@ -178,6 +189,19 @@    cosh   = lift1 (\b -> snd . Pos.cosSinh b)    -} +{-+The way EqDec and OrdDec are instantiated+it is possible to have different bases+for the arguments for comparison+and the arguments between we decide.+However, I would not rely on this.+-}+instance EqDec.C T where+   x==?y  =  lift2 (\b -> Pos.ifLazy b (x==y))++instance OrdDec.C T where+   x<=?y  =  lift2 (\b -> Pos.ifLazy b (x<=y))+ instance ZeroTestable.C T where    isZero (Cons xb xe xm) =       Pos.cmp xb (xe,xm) Pos.zero == EQ@@ -190,14 +214,16 @@    compare (Cons xb xe xm) (Cons yb ye ym) =       Pos.cmp (commonBasis xb yb) (xe,xm) (ye,ym) -instance Real.C T where+instance Absolute.C T where    abs = lift1 (const Pos.absolute)-   -- use default implementation for signum+   signum = Absolute.signumOrd -instance RealField.C T where+instance RealRing.C T where    splitFraction (Cons xb xe xm) =       let (int, frac) = Pos.toFixedPoint xb (xe,xm)       in  (fromInteger int, Cons xb (-1) frac)++instance RealField.C T where  instance RealTrans.C T where    atan2  = lift2 (curry . Pos.angle)
src/Number/Quaternion.hs view
@@ -61,8 +61,8 @@ import qualified Data.Array as Array  import qualified Prelude as P-import PreludeBase-import NumericPrelude hiding (signum)+import NumericPrelude.Base+import NumericPrelude.Numeric hiding (signum) import Text.Show.HT (showsInfixPrec, ) import Text.Read.HT (readsInfixPrec, )
src/Number/Ratio.hs view
@@ -25,7 +25,7 @@  import qualified Algebra.PrincipalIdealDomain as PID import qualified Algebra.Units                as Units-import qualified Algebra.Real                 as Real+import qualified Algebra.Absolute                 as Absolute import qualified Algebra.Ring                 as Ring import qualified Algebra.Additive             as Additive import qualified Algebra.ZeroTestable         as ZeroTestable@@ -50,7 +50,7 @@ import qualified Data.Ratio as Ratio98  import qualified Prelude as P-import PreludeBase+import NumericPrelude.Base   infixl 7 %@@ -72,7 +72,7 @@        (n:%d) = s%y    in  ((n*x):%d) -{- | similar to 'Algebra.RealField.splitFraction' -}+{- | similar to 'Algebra.RealRing.splitFraction' -} split :: (PID.C a) => T a -> (a, T a) split (x:%y) =    let (q,r) = divMod x y@@ -101,9 +101,9 @@     fromInteger x       =  fromValue $ fromInteger x     (x:%y) * (x':%y')   =  (x * x') % (y * y') -instance (Real.C a, PID.C a) => Real.C (T a) where-    abs (x:%y)          =  Real.abs x :% y-    signum (x:%_)       =  Real.signum x :% one+instance (Absolute.C a, PID.C a) => Absolute.C (T a) where+    abs (x:%y)          =  Absolute.abs x :% y+    signum (x:%_)       =  Absolute.signum x :% one   liftOrd :: Ring.C a => (a -> a -> b) -> (T a -> T a -> b)@@ -208,7 +208,7 @@ -- * Legacy Instances  --- | Necessary when mixing NumericPrelude Rationals with Prelude98 Rationals+-- | Necessary when mixing NumericPrelude.Numeric Rationals with Prelude98 Rationals  toRational98 :: (P.Integral a, PID.C a) => T a -> Ratio98.Ratio a toRational98 x = numerator x Ratio98.% denominator x@@ -220,16 +220,16 @@   -- instance (P.Num a, PID.C a) => P.Num (T a) where-instance (P.Num a, PID.C a, Real.C a) => P.Num (T a) where+instance (P.Num a, PID.C a, Absolute.C a) => P.Num (T a) where    fromInteger n = P.fromInteger n % 1    negate = negate -- for unary minus    (+)    = legacyInstance "(+)"    (*)    = legacyInstance "(*)"-   abs    = Real.abs -- needed for Arbitrary instance of NonNegative.Ratio+   abs    = Absolute.abs -- needed for Arbitrary instance of NonNegative.Ratio    signum = legacyInstance "signum"  -- instance (P.Num a, PID.C a) => P.Fractional (T a) where-instance (P.Num a, PID.C a, Real.C a) => P.Fractional (T a) where+instance (P.Num a, PID.C a, Absolute.C a) => P.Fractional (T a) where --   fromRational = Field.fromRational    fromRational x =       fromInteger (Ratio98.numerator x) :%
src/Number/ResidueClass.hs view
@@ -8,8 +8,8 @@  import Algebra.ZeroTestable(isZero) -import PreludeBase-import NumericPrelude hiding (recip)+import NumericPrelude.Base+import NumericPrelude.Numeric hiding (recip) import Data.Maybe.HT (toMaybe) import Data.Maybe (fromMaybe)
src/Number/ResidueClass/Check.hs view
@@ -12,8 +12,8 @@  import Algebra.ZeroTestable(isZero) -import PreludeBase-import NumericPrelude (Int, Integer, mod, )+import NumericPrelude.Base+import NumericPrelude.Numeric (Int, Integer, mod, ) import Data.Maybe.HT (toMaybe, ) import Text.Show.HT (showsInfixPrec, ) import Text.Read.HT (readsInfixPrec, )
src/Number/ResidueClass/Func.hs view
@@ -4,16 +4,18 @@ import qualified Number.ResidueClass as Res  import qualified Algebra.PrincipalIdealDomain as PID-import qualified Algebra.IntegralDomain as Integral-import qualified Algebra.Field          as Field-import qualified Algebra.Ring           as Ring-import qualified Algebra.Additive       as Additive+import qualified Algebra.IntegralDomain   as Integral+import qualified Algebra.Field            as Field+import qualified Algebra.Ring             as Ring+import qualified Algebra.Additive         as Additive+import qualified Algebra.EqualityDecision as EqDec -import PreludeBase-import NumericPrelude hiding (zero, one, )+import Algebra.EqualityDecision ((==?), )+import NumericPrelude.Base+import NumericPrelude.Numeric hiding (zero, one, )  import qualified Prelude        as P-import qualified NumericPrelude as NP+import qualified NumericPrelude.Numeric as NP  {- | Here a residue class is a representative@@ -53,6 +55,10 @@ equal :: Eq a => a -> T a -> T a -> Bool equal m (Cons x) (Cons y)  =  x m == y m ++instance  (EqDec.C a) => EqDec.C (T a)  where+    (==?) (Cons x) (Cons y) (Cons eq) (Cons noteq) =+       Cons (\m -> (x m ==? y m) (eq m) (noteq m))  instance  (Integral.C a) => Additive.C (T a)  where     zero		=  zero
src/Number/ResidueClass/Maybe.hs view
@@ -10,8 +10,8 @@  import Algebra.ZeroTestable(isZero) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  infix 7 /:, `Cons`
src/Number/ResidueClass/Reader.hs view
@@ -8,19 +8,19 @@ import qualified Algebra.Ring           as Ring import qualified Algebra.Additive       as Additive -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric  import Control.Monad (liftM2, liftM4) -- import Control.Monad.Reader (MonadReader)  import qualified Prelude        as P-import qualified NumericPrelude as NP+import qualified NumericPrelude.Numeric as NP   {- | T is a Reader monad but does not need functional dependencies-like that from the Monad Template Library.+like that from the Monad Transformer Library. -} newtype T a b = Cons {toFunc :: a -> b}
src/Number/SI.hs view
@@ -33,7 +33,7 @@ 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.Absolute                as Absolute import qualified Algebra.Ring                as Ring import qualified Algebra.Additive            as Additive import qualified Algebra.ZeroTestable        as ZeroTestable@@ -44,8 +44,8 @@  import qualified Prelude as P -import NumericPrelude-import PreludeBase+import NumericPrelude.Numeric+import NumericPrelude.Base   newtype T a v = Cons (PValue v)@@ -145,7 +145,7 @@   (<=)    = lift2Gen (<=)   (>=)    = lift2Gen (>=) -instance (Real.C v) => Real.C (T a v) where+instance (Absolute.C v) => Absolute.C (T a v) where   abs    = lift abs   signum = lift signum
src/Number/SI/Unit.hs view
@@ -20,8 +20,8 @@ import Number.Physical.UnitDatabase(initScale, initUnitSet) import Data.Maybe(catMaybes) -import PreludeBase hiding (length)-import NumericPrelude hiding (one)+import NumericPrelude.Base hiding (length)+import NumericPrelude.Numeric hiding (one)  data Dimension =    Length | Time | Mass | Charge |
src/NumericPrelude.hs view
@@ -1,44 +1,9 @@-{-# LANGUAGE NoImplicitPrelude #-}-module NumericPrelude (-    {- Additive -} (+), (-), negate, zero, subtract, sum, sum1,-    {- ZeroTestable -} isZero,-    {- Ring -} (*), one, fromInteger, (^), ringPower, sqr, product, product1,-    {- IntegralDomain -} div, mod, divMod, divides, even, odd,-    {- Field -} (/), recip, fromRational', (^-), fieldPower, fromRational,-    {- Algebraic -} (^/), sqrt,-    {- Transcendental -}-        pi, exp, log, logBase, (**), (^?), sin, cos, tan,-        asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh,-    {- Real -} abs, signum,-    {- RealIntegral -} quot, rem, quotRem,-    {- RealFrac -} splitFraction, fraction, truncate, round, ceiling, floor, approxRational,-    {- RealTrans -} atan2,-    {- ToRational -} toRational,-    {- ToInteger -} toInteger, fromIntegral,-    {- Units -} isUnit, stdAssociate, stdUnit, stdUnitInv,-    {- PID -} extendedGCD, gcd, lcm, euclid, extendedEuclid,-    {- Ratio -} Rational, (%), numerator, denominator,-    Integer, Int, Float, Double,-    {- Module -} (*>)-) where--import Number.Ratio (Rational, (%), numerator, denominator)--import Algebra.Module((*>))-import Algebra.RealTranscendental(atan2)-import Algebra.Transcendental-import Algebra.Algebraic((^/), sqrt)-import Algebra.RealField(splitFraction, fraction, truncate, round, ceiling, floor, approxRational, )-import Algebra.Field((/), (^-), recip, fromRational', fromRational, )-import Algebra.PrincipalIdealDomain (extendedGCD, gcd, lcm, euclid, extendedEuclid)-import Algebra.Units (isUnit, stdAssociate, stdUnit, stdUnitInv)-import Algebra.RealIntegral (quot, rem, quotRem, )-import Algebra.IntegralDomain (div, mod, divMod, divides, even, odd)-import Algebra.Real (abs, signum, )-import Algebra.Ring (one, fromInteger, (*), (^), sqr, product, product1)-import Algebra.Additive (zero, (+), (-), negate, subtract, sum, sum1)-import Algebra.ZeroTestable (isZero)-import Algebra.ToInteger (ringPower, fieldPower, toInteger, fromIntegral, )-import Algebra.ToRational (toRational, )+module NumericPrelude+   (module NumericPrelude.Numeric,+    module NumericPrelude.Base,+    max, min, abs, ) where -import Prelude (Int, Integer, Float, Double)+import NumericPrelude.Numeric hiding (abs, )+import NumericPrelude.Base    hiding (max, min, )+import Prelude ()+import Algebra.Lattice (max, min, abs, )
+ src/NumericPrelude/Base.hs view
@@ -0,0 +1,12 @@+{- |+The only point of this module is+to reexport items that we want from the standard Prelude.+-}++module NumericPrelude.Base (module Prelude) where+import Prelude hiding (+       Int, Integer, Float, Double, Rational, Num(..), Real(..),+       Integral(..), Fractional(..), Floating(..), RealFrac(..),+       RealFloat(..), subtract, even, odd,+       gcd, lcm, (^), (^^), sum, product,+       fromIntegral, fromRational, )
src/NumericPrelude/Elementwise.hs view
@@ -28,6 +28,19 @@ run2 :: T (x,y) a -> x -> y -> a run2 = curry . run +{-# INLINE run3 #-}+run3 :: T (x,y,z) a -> x -> y -> z -> a+run3 e x y z = run e (x,y,z)++{-# INLINE run4 #-}+run4 :: T (x,y,z,w) a -> x -> y -> z -> w -> a+run4 e x y z w = run e (x,y,z,w)++{-# INLINE run5 #-}+run5 :: T (x,y,z,u,w) a -> x -> y -> z -> u -> w -> a+run5 e x y z u w = run e (x,y,z,u,w)++ instance Functor (T v) where    {-# INLINE fmap #-}    fmap f (Cons e) =
+ src/NumericPrelude/List/Checked.hs view
@@ -0,0 +1,94 @@+{-# LANGUAGE NoImplicitPrelude #-}+{- |+Some functions that are counterparts of functions from "Data.List"+using NumericPrelude.Numeric type classes.+They are distinct in that they check for valid arguments,+e.g. the length argument of 'take' must be at most the length of the input list.+However, since many Haskell programs rely on the absence of such checks,+we did not make these the default implementations+as in "NumericPrelude.List.Generic".+-}+module NumericPrelude.List.Checked+   (take, drop, splitAt, (!!), zipWith,+   ) where++import qualified Algebra.ToInteger  as ToInteger+import qualified Algebra.Ring       as Ring+import Algebra.Ring (one, )+import Algebra.Additive (zero, (-), )++import Data.Tuple.HT (mapFst, )++import qualified NumericPrelude.List as NPList++import NumericPrelude.Base hiding (take, drop, splitAt, length, replicate, (!!), zipWith, )+++moduleError :: String -> String -> a+moduleError name msg =+   error $ "NumericPrelude.List.Left." ++ name ++ ": " ++ msg++{- |+Taken number of elements must be at most the length of the list,+otherwise the end of the list is undefined.+-}+take :: (ToInteger.C n) => n -> [a] -> [a]+take n =+   if n<=zero+     then const []+     else \xt ->+       case xt of+          [] -> moduleError "take" "index out of range"+          (x:xs) -> x : take (n-one) xs++{- |+Dropped number of elements must be at most the length of the list,+otherwise the end of the list is undefined.+-}+drop :: (ToInteger.C n) => n -> [a] -> [a]+drop n =+   if n<=zero+     then id+     else \xt ->+       case xt of+          [] -> moduleError "drop" "index out of range"+          (_:xs) -> drop (n-one) xs++{- |+Split position must be at most the length of the list,+otherwise the end of the first list and the second list are undefined.+-}+splitAt :: (ToInteger.C n) => n -> [a] -> ([a], [a])+splitAt n xt =+   if n<=zero+     then ([], xt)+     else+       case xt of+          [] -> moduleError "splitAt" "index out of range"+          (x:xs) -> mapFst (x:) $ splitAt (n-one) xs++{- |+The index must be smaller than the length of the list,+otherwise the result is undefined.+-}+(!!) :: (ToInteger.C n) => [a] -> n -> a+(!!) [] _ = moduleError "(!!)" "index out of range"+(!!) (x:xs) n =+   if n<=zero+     then x+     else (!!) xs (n-one)+++{- |+Zip two lists which must be of the same length.+This is checked only lazily, that is unequal lengths are detected only+if the list is evaluated completely.+But it is more strict than @zipWithPad undefined f@+since the latter one may succeed on unequal length list if @f@ is lazy.+-}+zipWith+   :: (a -> b -> c)   {-^ function applied to corresponding elements of the lists -}+   -> [a]+   -> [b]+   -> [c]+zipWith = NPList.zipWithMatch
+ src/NumericPrelude/List/Generic.hs view
@@ -0,0 +1,84 @@+{-# LANGUAGE NoImplicitPrelude #-}+{- |+Functions that are counterparts of the @generic@ functions in "Data.List"+using NumericPrelude.Numeric type classes.+For input arguments we use the restrictive @ToInteger@ constraint,+although in principle @RealRing@ would be enough.+However we think that @take 0.5 xs@ is rather a bug than a feature,+thus we forbid fractional types.+On the other hand fractional types as result can be quite handy,+e.g. in @average xs = sum xs / length xs@.+-}+module NumericPrelude.List.Generic+   ((!!), lengthLeft, lengthRight, replicate,+    take, drop, splitAt,+    findIndex, elemIndex, findIndices, elemIndices,+   ) where++import NumericPrelude.List.Checked ((!!), )++import qualified Algebra.ToInteger  as ToInteger+import qualified Algebra.Ring       as Ring+import Algebra.Ring (one, )+import Algebra.Additive (zero, (+), (-), )++import qualified Data.Maybe         as Maybe+import Data.Tuple.HT (mapFst, )++import NumericPrelude.Base as List+   hiding (take, drop, splitAt, length, replicate, (!!), )+++replicate :: (ToInteger.C n) => n -> a -> [a]+replicate n x = take n (List.repeat x)++take :: (ToInteger.C n) => n -> [a] -> [a]+take _ [] = []+take n (x:xs) =+   if n<=zero+     then []+     else x : take (n-one) xs++drop :: (ToInteger.C n) => n -> [a] -> [a]+drop _ [] = []+drop n xt@(_:xs) =+   if n<=zero+     then xt+     else drop (n-one) xs++splitAt :: (ToInteger.C n) => n -> [a] -> ([a], [a])+splitAt _ [] = ([], [])+splitAt n xt@(x:xs) =+   if n<=zero+     then ([], xt)+     else mapFst (x:) $ splitAt (n-one) xs+++{- |+Left associative length computation+that is appropriate for types like @Integer@.+-}+lengthLeft :: (Ring.C n) => [a] -> n+lengthLeft = List.foldl (\n _ -> n+one) zero++{- |+Right associative length computation+that is appropriate for types like @Peano@ number.+-}+lengthRight :: (Ring.C n) => [a] -> n+lengthRight = List.foldr (\_ n -> one+n) zero++elemIndex :: (Ring.C n, Eq a) => a -> [a] -> Maybe n+elemIndex e = findIndex (e==)++elemIndices :: (Ring.C n, Eq a) => a -> [a] -> [n]+elemIndices e = findIndices (e==)++findIndex :: Ring.C n => (a -> Bool) -> [a] -> Maybe n+findIndex p = Maybe.listToMaybe . findIndices p++findIndices :: Ring.C n => (a -> Bool) -> [a] -> [n]+findIndices p =+   map fst .+   filter (p . snd) .+   zip (iterate (one+) zero)
+ src/NumericPrelude/Numeric.hs view
@@ -0,0 +1,44 @@+{-# LANGUAGE NoImplicitPrelude #-}+module NumericPrelude.Numeric (+    {- Additive -} (+), (-), negate, zero, subtract, sum, sum1,+    {- ZeroTestable -} isZero,+    {- Ring -} (*), one, fromInteger, (^), ringPower, sqr, product, product1,+    {- IntegralDomain -} div, mod, divMod, divides, even, odd,+    {- Field -} (/), recip, fromRational', (^-), fieldPower, fromRational,+    {- Algebraic -} (^/), sqrt,+    {- Transcendental -}+        pi, exp, log, logBase, (**), (^?), sin, cos, tan,+        asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh,+    {- Absolute -} abs, signum,+    {- RealIntegral -} quot, rem, quotRem,+    {- RealFrac -} splitFraction, fraction, truncate, round, ceiling, floor, approxRational,+    {- RealTrans -} atan2,+    {- ToRational -} toRational,+    {- ToInteger -} toInteger, fromIntegral,+    {- Units -} isUnit, stdAssociate, stdUnit, stdUnitInv,+    {- PID -} extendedGCD, gcd, lcm, euclid, extendedEuclid,+    {- Ratio -} Rational, (%), numerator, denominator,+    Integer, Int, Float, Double,+    {- Module -} (*>)+) where++import Number.Ratio (Rational, (%), numerator, denominator)++import Algebra.Module((*>))+import Algebra.RealTranscendental(atan2)+import Algebra.Transcendental+import Algebra.Algebraic((^/), sqrt)+import Algebra.RealRing(splitFraction, fraction, truncate, round, ceiling, floor, approxRational, )+import Algebra.Field((/), (^-), recip, fromRational', fromRational, )+import Algebra.PrincipalIdealDomain (extendedGCD, gcd, lcm, euclid, extendedEuclid)+import Algebra.Units (isUnit, stdAssociate, stdUnit, stdUnitInv)+import Algebra.RealIntegral (quot, rem, quotRem, )+import Algebra.IntegralDomain (div, mod, divMod, divides, even, odd)+import Algebra.Absolute (abs, signum, )+import Algebra.Ring (one, fromInteger, (*), (^), sqr, product, product1)+import Algebra.Additive (zero, (+), (-), negate, subtract, sum, sum1)+import Algebra.ZeroTestable (isZero)+import Algebra.ToInteger (ringPower, fieldPower, toInteger, fromIntegral, )+import Algebra.ToRational (toRational, )++import Prelude (Int, Integer, Float, Double)
− src/PreludeBase.hs 
@@ -1,12 +0,0 @@-{- |-The only point of this module is-to reexport items that we want from the standard Prelude.--}--module PreludeBase (module Prelude) where-import Prelude hiding(-       Int, Integer, Float, Double, Rational, Num(..), Real(..),-       Integral(..), Fractional(..), Floating(..), RealFrac(..),-       RealFloat(..), subtract, even, odd,-       gcd, lcm, (^), (^^), sum, product,-       fromIntegral, fromRational)
test/Test.hs view
@@ -10,7 +10,7 @@            deca, hecto, kilo, mega, giga, tera, peta) import Number.OccasionallyScalarExpression as Expr -import qualified Number.Positional.Check  as Real+import qualified Number.Positional.Check  as Absolute import qualified Number.FixedPoint.Check  as FixedPoint import qualified Number.ResidueClass.Func as ResidueClass import qualified Number.Peano             as Peano@@ -28,8 +28,8 @@  import Data.List (genericTake, genericLength) -import PreludeBase-import NumericPrelude+import NumericPrelude.Base+import NumericPrelude.Numeric   {- * Physical units -}@@ -77,13 +77,13 @@ {- * Reals -}  testReal :: String-testReal = Real.defltShow (sqrt 2 + log 2 * pi)+testReal = Absolute.defltShow (sqrt 2 + log 2 * pi) -testComplexReal :: Complex.T Real.T+testComplexReal :: Complex.T Absolute.T testComplexReal = exp (0 +: pi) + exp (0 -: pi) -showReal :: Real.T -> String-showReal = Real.defltShow+showReal :: Absolute.T -> String+showReal = Absolute.defltShow   {- * Fixed point numbers -}
test/Test/MathObj/Gaussian/Bell.hs view
@@ -17,8 +17,8 @@  import Data.Function.HT (nest, ) -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   simple ::
test/Test/MathObj/Gaussian/Polynomial.hs view
@@ -24,8 +24,8 @@  -- import Debug.Trace (trace, ) -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   simple ::
test/Test/MathObj/Gaussian/Variance.hs view
@@ -15,8 +15,8 @@  import Data.Function.HT (nest, ) -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   simple ::
test/Test/MathObj/Matrix.hs view
@@ -18,8 +18,8 @@ import qualified Test.HUnit as HUnit  -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   type Seed = Int
test/Test/MathObj/PartialFraction.hs view
@@ -22,8 +22,8 @@ import qualified Test.HUnit as HUnit  -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   {- * Properties for generic types -}
test/Test/MathObj/Polynomial.hs view
@@ -1,9 +1,8 @@ {-# LANGUAGE NoImplicitPrelude #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE FlexibleInstances #-} module Test.MathObj.Polynomial where -import qualified MathObj.Polynomial as Poly+import qualified MathObj.Polynomial      as Poly+import qualified MathObj.Polynomial.Core as PolyCore  import qualified Algebra.IntegralDomain as Integral import qualified Algebra.Ring           as Ring@@ -18,18 +17,18 @@ import qualified Test.HUnit as HUnit  -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   tensorProductTranspose :: (Ring.C a, Eq a) => [a] -> [a] -> Property tensorProductTranspose xs ys =    not (null xs) && not (null ys) ==>-      Poly.tensorProduct xs ys == List.transpose (Poly.tensorProduct ys xs)+      PolyCore.tensorProduct xs ys == List.transpose (PolyCore.tensorProduct ys xs)   mul :: (Ring.C a, Eq a, ZeroTestable.C a) => [a] -> [a] -> Bool-mul xs ys  =  Poly.equal (Poly.mul xs ys) (Poly.mulShear xs ys)+mul xs ys  =  PolyCore.equal (PolyCore.mul xs ys) (PolyCore.mulShear xs ys)   test :: Testable a => (Poly.T Integer -> a) -> IO ()
test/Test/MathObj/PowerSeries.hs view
@@ -3,7 +3,7 @@ {-# LANGUAGE FlexibleInstances #-} module Test.MathObj.PowerSeries where -import qualified MathObj.PowerSeries         as PS+import qualified MathObj.PowerSeries.Core    as PS import qualified MathObj.PowerSeries.Example as PSE  import Test.NumericPrelude.Utility (equalInfLists {- , testUnit -} )@@ -11,8 +11,8 @@ import qualified Test.HUnit as HUnit  -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   identitiesExplODE, identitiesSeriesFunction, identitiesInverses ::
+ test/Test/MathObj/RefinementMask2.hs view
@@ -0,0 +1,78 @@+{-# LANGUAGE NoImplicitPrelude #-}+module Test.MathObj.RefinementMask2 where++import qualified MathObj.RefinementMask2 as Mask+import qualified Algebra.Differential    as D++import qualified MathObj.Polynomial      as Poly+import qualified MathObj.Polynomial.Core as PolyCore++import qualified Algebra.RealField      as RealField+import qualified Algebra.Ring           as Ring++import qualified Algebra.ZeroTestable   as ZeroTestable++import Data.Maybe (fromMaybe, )++import Test.NumericPrelude.Utility (testUnit)+import Test.QuickCheck (Property, quickCheck, (==>), Testable, )+import qualified Test.HUnit as HUnit+++import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP++++hasMultipleZero :: (Ring.C a, Eq a) => Int -> a -> Poly.T a -> Bool+hasMultipleZero n x poly =+   all (zero==) $ take n $+   map (flip Poly.evaluate x) $+   iterate D.differentiate poly++inverse0 :: (RealField.C a) => Mask.T a -> Property+inverse0 mask0 =+   let (b,poly) =+          case Mask.toPolynomial mask0 of+             Just p -> (True, p)+             Nothing -> (False, error "RefinementMask2.inverse0: no admissible mask")+       mask1 = Mask.fromPolynomial poly+       maskD =+          Poly.fromCoeffs (Mask.coeffs mask1) -+          Poly.fromCoeffs (Mask.coeffs mask0)+   in  b ==>+          hasMultipleZero (fromMaybe 0 $ Poly.degree poly)+             1 maskD++truncatePolynomial :: (ZeroTestable.C a) => Int -> Poly.T a -> Poly.T a+truncatePolynomial n =+   Poly.fromCoeffs . PolyCore.normalize . take n . Poly.coeffs++inverse1 :: (RealField.C a) => Poly.T a -> Bool+inverse1 poly0 =+   case Mask.toPolynomial (Mask.fromPolynomial poly0) of+      Just poly1 -> Poly.collinear poly0 poly1+      Nothing -> False++refining :: (RealField.C a) => Poly.T a -> Bool+refining poly =+   poly == Mask.refinePolynomial (Mask.fromPolynomial poly) poly++++test :: Testable a => (Poly.T Integer -> a) -> IO ()+test = quickCheck++testRat :: Testable a => (Poly.T Rational -> a) -> IO ()+testRat = quickCheck+++tests :: HUnit.Test+tests =+   HUnit.TestLabel "refinement mask" $+   HUnit.TestList $+   map testUnit $+      ("inverse0", quickCheck (inverse0 :: Mask.T Rational -> Property)) :+      ("inverse1", quickCheck (inverse1 . truncatePolynomial 5 :: Poly.T Rational -> Bool)) :+      ("refining", quickCheck (refining . truncatePolynomial 5 :: Poly.T Rational -> Bool)) :+      []
test/Test/Number/GaloisField2p32m5.hs view
@@ -10,8 +10,8 @@ import qualified Test.HUnit as HUnit  -import PreludeBase as P-import NumericPrelude as NP+import NumericPrelude.Base as P+import NumericPrelude.Numeric as NP   test :: Testable a => (GF.T -> a) -> IO ()
test/Test/Run.hs view
@@ -1,5 +1,6 @@ module Main where +import qualified Test.MathObj.RefinementMask2 as RefinementMask2 import qualified Test.Algebra.RealRing as RealRing import qualified Test.MathObj.Gaussian.Polynomial as GaussPoly import qualified Test.MathObj.Gaussian.Variance as GaussVariance@@ -15,6 +16,7 @@ main :: IO () main =    do HUnitText.runTestTT (HUnit.TestList $+         RefinementMask2.tests :          RealRing.tests :          GaussVariance.tests :          GaussBell.tests :