diff --git a/docs/NOTES b/docs/NOTES
--- a/docs/NOTES
+++ b/docs/NOTES
@@ -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
 
diff --git a/numeric-prelude.cabal b/numeric-prelude.cabal
--- a/numeric-prelude.cabal
+++ b/numeric-prelude.cabal
@@ -1,5 +1,5 @@
 Name:           numeric-prelude
-Version:        0.1.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
diff --git a/src/Algebra/Absolute.hs b/src/Algebra/Absolute.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/Absolute.hs
@@ -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
+
diff --git a/src/Algebra/Additive.hs b/src/Algebra/Additive.hs
--- a/src/Algebra/Additive.hs
+++ b/src/Algebra/Additive.hs
@@ -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
diff --git a/src/Algebra/AffineSpace.hs b/src/Algebra/AffineSpace.hs
--- a/src/Algebra/AffineSpace.hs
+++ b/src/Algebra/AffineSpace.hs
@@ -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 ()
 
 {- |
diff --git a/src/Algebra/Algebraic.hs b/src/Algebra/Algebraic.hs
--- a/src/Algebra/Algebraic.hs
+++ b/src/Algebra/Algebraic.hs
@@ -12,7 +12,7 @@
 import Algebra.Ring ((*), (^), fromInteger)
 import Algebra.Additive((+))
 
-import PreludeBase
+import NumericPrelude.Base
 import qualified Prelude as P
 
 
diff --git a/src/Algebra/Differential.hs b/src/Algebra/Differential.hs
--- a/src/Algebra/Differential.hs
+++ b/src/Algebra/Differential.hs
@@ -3,7 +3,7 @@
 
 import qualified Algebra.Ring as Ring
 
--- import NumericPrelude
+-- import NumericPrelude.Numeric
 import qualified Prelude
 
 {- |
diff --git a/src/Algebra/EqualityDecision.hs b/src/Algebra/EqualityDecision.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/EqualityDecision.hs
@@ -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)
diff --git a/src/Algebra/Field.hs b/src/Algebra/Field.hs
--- a/src/Algebra/Field.hs
+++ b/src/Algebra/Field.hs
@@ -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)
diff --git a/src/Algebra/GenerateRules.hs b/src/Algebra/GenerateRules.hs
--- a/src/Algebra/GenerateRules.hs
+++ b/src/Algebra/GenerateRules.hs
@@ -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 >>
 
diff --git a/src/Algebra/IntegralDomain.hs b/src/Algebra/IntegralDomain.hs
--- a/src/Algebra/IntegralDomain.hs
+++ b/src/Algebra/IntegralDomain.hs
@@ -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
 
diff --git a/src/Algebra/Lattice.hs b/src/Algebra/Lattice.hs
--- a/src/Algebra/Lattice.hs
+++ b/src/Algebra/Lattice.hs
@@ -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`
diff --git a/src/Algebra/Module.hs b/src/Algebra/Module.hs
--- a/src/Algebra/Module.hs
+++ b/src/Algebra/Module.hs
@@ -60,6 +60,7 @@
     (*>) :: a -> v -> v
 
 
+{-# INLINE (<*>.*>) #-}
 (<*>.*>) ::
    (C a x) =>
    Elem.T (a,v) (x -> c) -> (v -> x) -> Elem.T (a,v) c
diff --git a/src/Algebra/Monoid.hs b/src/Algebra/Monoid.hs
--- a/src/Algebra/Monoid.hs
+++ b/src/Algebra/Monoid.hs
@@ -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
diff --git a/src/Algebra/NonNegative.hs b/src/Algebra/NonNegative.hs
--- a/src/Algebra/NonNegative.hs
+++ b/src/Algebra/NonNegative.hs
@@ -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
+-}
diff --git a/src/Algebra/NormedSpace/Euclidean.hs b/src/Algebra/NormedSpace/Euclidean.hs
--- a/src/Algebra/NormedSpace/Euclidean.hs
+++ b/src/Algebra/NormedSpace/Euclidean.hs
@@ -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
diff --git a/src/Algebra/NormedSpace/Maximum.hs b/src/Algebra/NormedSpace/Maximum.hs
--- a/src/Algebra/NormedSpace/Maximum.hs
+++ b/src/Algebra/NormedSpace/Maximum.hs
@@ -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
+-}
diff --git a/src/Algebra/NormedSpace/Sum.hs b/src/Algebra/NormedSpace/Sum.hs
--- a/src/Algebra/NormedSpace/Sum.hs
+++ b/src/Algebra/NormedSpace/Sum.hs
@@ -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
diff --git a/src/Algebra/OccasionallyScalar.hs b/src/Algebra/OccasionallyScalar.hs
--- a/src/Algebra/OccasionallyScalar.hs
+++ b/src/Algebra/OccasionallyScalar.hs
@@ -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)
diff --git a/src/Algebra/OrderDecision.hs b/src/Algebra/OrderDecision.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/OrderDecision.hs
@@ -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)
diff --git a/src/Algebra/PrincipalIdealDomain.hs b/src/Algebra/PrincipalIdealDomain.hs
--- a/src/Algebra/PrincipalIdealDomain.hs
+++ b/src/Algebra/PrincipalIdealDomain.hs
@@ -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)
diff --git a/src/Algebra/Real.hs b/src/Algebra/Real.hs
deleted file mode 100644
--- a/src/Algebra/Real.hs
+++ /dev/null
@@ -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
-
diff --git a/src/Algebra/RealField.hs b/src/Algebra/RealField.hs
--- a/src/Algebra/RealField.hs
+++ b/src/Algebra/RealField.hs
@@ -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
diff --git a/src/Algebra/RealIntegral.hs b/src/Algebra/RealIntegral.hs
--- a/src/Algebra/RealIntegral.hs
+++ b/src/Algebra/RealIntegral.hs
@@ -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)
 
diff --git a/src/Algebra/RealRing.hs b/src/Algebra/RealRing.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/RealRing.hs
@@ -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
diff --git a/src/Algebra/RealTranscendental.hs b/src/Algebra/RealTranscendental.hs
--- a/src/Algebra/RealTranscendental.hs
+++ b/src/Algebra/RealTranscendental.hs
@@ -12,7 +12,7 @@
 import Data.Bool.HT (select, )
 
 import qualified Prelude as P
-import PreludeBase
+import NumericPrelude.Base
 
 
 
diff --git a/src/Algebra/RightModule.hs b/src/Algebra/RightModule.hs
--- a/src/Algebra/RightModule.hs
+++ b/src/Algebra/RightModule.hs
@@ -6,7 +6,7 @@
 import qualified Algebra.Ring     as Ring
 import qualified Algebra.Additive as Additive
 
--- import NumericPrelude
+-- import NumericPrelude.Numeric
 import qualified Prelude
 
 
diff --git a/src/Algebra/Ring.hs b/src/Algebra/Ring.hs
--- a/src/Algebra/Ring.hs
+++ b/src/Algebra/Ring.hs
@@ -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
diff --git a/src/Algebra/ToInteger.hs b/src/Algebra/ToInteger.hs
--- a/src/Algebra/ToInteger.hs
+++ b/src/Algebra/ToInteger.hs
@@ -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.
diff --git a/src/Algebra/ToRational.hs b/src/Algebra/ToRational.hs
--- a/src/Algebra/ToRational.hs
+++ b/src/Algebra/ToRational.hs
@@ -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;
+  #-}
diff --git a/src/Algebra/Transcendental.hs b/src/Algebra/Transcendental.hs
--- a/src/Algebra/Transcendental.hs
+++ b/src/Algebra/Transcendental.hs
@@ -13,7 +13,7 @@
 import Algebra.Additive  ((+), (-), negate)
 
 import qualified Prelude as P
-import PreludeBase
+import NumericPrelude.Base
 
 
 infixr 8  **, ^?
diff --git a/src/Algebra/Units.hs b/src/Algebra/Units.hs
--- a/src/Algebra/Units.hs
+++ b/src/Algebra/Units.hs
@@ -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
diff --git a/src/Algebra/VectorSpace.hs b/src/Algebra/VectorSpace.hs
--- a/src/Algebra/VectorSpace.hs
+++ b/src/Algebra/VectorSpace.hs
@@ -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
 
 
diff --git a/src/Algebra/ZeroTestable.hs b/src/Algebra/ZeroTestable.hs
--- a/src/Algebra/ZeroTestable.hs
+++ b/src/Algebra/ZeroTestable.hs
@@ -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:
diff --git a/src/MathObj/Algebra.hs b/src/MathObj/Algebra.hs
--- a/src/MathObj/Algebra.hs
+++ b/src/MathObj/Algebra.hs
@@ -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)
 
 
diff --git a/src/MathObj/DiscreteMap.hs b/src/MathObj/DiscreteMap.hs
--- a/src/MathObj/DiscreteMap.hs
+++ b/src/MathObj/DiscreteMap.hs
@@ -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.
diff --git a/src/MathObj/Gaussian/Bell.hs b/src/MathObj/Gaussian/Bell.hs
--- a/src/MathObj/Gaussian/Bell.hs
+++ b/src/MathObj/Gaussian/Bell.hs
@@ -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))
diff --git a/src/MathObj/Gaussian/Example.hs b/src/MathObj/Gaussian/Example.hs
--- a/src/MathObj/Gaussian/Example.hs
+++ b/src/MathObj/Gaussian/Example.hs
@@ -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
 
 
diff --git a/src/MathObj/Gaussian/Polynomial.hs b/src/MathObj/Gaussian/Polynomial.hs
--- a/src/MathObj/Gaussian/Polynomial.hs
+++ b/src/MathObj/Gaussian/Polynomial.hs
@@ -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,
diff --git a/src/MathObj/Gaussian/Variance.hs b/src/MathObj/Gaussian/Variance.hs
--- a/src/MathObj/Gaussian/Variance.hs
+++ b/src/MathObj/Gaussian/Variance.hs
@@ -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))
 
diff --git a/src/MathObj/LaurentPolynomial.hs b/src/MathObj/LaurentPolynomial.hs
--- a/src/MathObj/LaurentPolynomial.hs
+++ b/src/MathObj/LaurentPolynomial.hs
@@ -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
diff --git a/src/MathObj/Matrix.hs b/src/MathObj/Matrix.hs
--- a/src/MathObj/Matrix.hs
+++ b/src/MathObj/Matrix.hs
@@ -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, )
 
 
 {- |
diff --git a/src/MathObj/Monoid.hs b/src/MathObj/Monoid.hs
--- a/src/MathObj/Monoid.hs
+++ b/src/MathObj/Monoid.hs
@@ -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.
diff --git a/src/MathObj/PartialFraction.hs b/src/MathObj/PartialFraction.hs
--- a/src/MathObj/PartialFraction.hs
+++ b/src/MathObj/PartialFraction.hs
@@ -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)
 
 
 
diff --git a/src/MathObj/Permutation.hs b/src/MathObj/Permutation.hs
--- a/src/MathObj/Permutation.hs
+++ b/src/MathObj/Permutation.hs
@@ -16,8 +16,8 @@
 
 import Data.Array(Ix)
 
--- import NumericPrelude (Integer)
-import PreludeBase
+-- import NumericPrelude.Numeric (Integer)
+import NumericPrelude.Base
 
 
 {- |
diff --git a/src/MathObj/Permutation/CycleList.hs b/src/MathObj/Permutation/CycleList.hs
--- a/src/MathObj/Permutation/CycleList.hs
+++ b/src/MathObj/Permutation/CycleList.hs
@@ -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]
diff --git a/src/MathObj/Permutation/CycleList/Check.hs b/src/MathObj/Permutation/CycleList/Check.hs
--- a/src/MathObj/Permutation/CycleList/Check.hs
+++ b/src/MathObj/Permutation/CycleList/Check.hs
@@ -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
diff --git a/src/MathObj/Permutation/Table.hs b/src/MathObj/Permutation/Table.hs
--- a/src/MathObj/Permutation/Table.hs
+++ b/src/MathObj/Permutation/Table.hs
@@ -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
diff --git a/src/MathObj/Polynomial.hs b/src/MathObj/Polynomial.hs
--- a/src/MathObj/Polynomial.hs
+++ b/src/MathObj/Polynomial.hs
@@ -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 -}
diff --git a/src/MathObj/Polynomial/Core.hs b/src/MathObj/Polynomial/Core.hs
new file mode 100644
--- /dev/null
+++ b/src/MathObj/Polynomial/Core.hs
@@ -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
+-}
+
diff --git a/src/MathObj/PowerSeries.hs b/src/MathObj/PowerSeries.hs
--- a/src/MathObj/PowerSeries.hs
+++ b/src/MathObj/PowerSeries.hs
@@ -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)
diff --git a/src/MathObj/PowerSeries/Core.hs b/src/MathObj/PowerSeries/Core.hs
new file mode 100644
--- /dev/null
+++ b/src/MathObj/PowerSeries/Core.hs
@@ -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)
diff --git a/src/MathObj/PowerSeries/DifferentialEquation.hs b/src/MathObj/PowerSeries/DifferentialEquation.hs
--- a/src/MathObj/PowerSeries/DifferentialEquation.hs
+++ b/src/MathObj/PowerSeries/DifferentialEquation.hs
@@ -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
 
 
 {- |
diff --git a/src/MathObj/PowerSeries/Example.hs b/src/MathObj/PowerSeries/Example.hs
--- a/src/MathObj/PowerSeries/Example.hs
+++ b/src/MathObj/PowerSeries/Example.hs
@@ -1,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. -}
diff --git a/src/MathObj/PowerSeries/Mean.hs b/src/MathObj/PowerSeries/Mean.hs
--- a/src/MathObj/PowerSeries/Mean.hs
+++ b/src/MathObj/PowerSeries/Mean.hs
@@ -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
diff --git a/src/MathObj/PowerSeries2.hs b/src/MathObj/PowerSeries2.hs
--- a/src/MathObj/PowerSeries2.hs
+++ b/src/MathObj/PowerSeries2.hs
@@ -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
diff --git a/src/MathObj/PowerSeries2/Core.hs b/src/MathObj/PowerSeries2/Core.hs
new file mode 100644
--- /dev/null
+++ b/src/MathObj/PowerSeries2/Core.hs
@@ -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
diff --git a/src/MathObj/PowerSum.hs b/src/MathObj/PowerSum.hs
--- a/src/MathObj/PowerSum.hs
+++ b/src/MathObj/PowerSum.hs
@@ -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
diff --git a/src/MathObj/RefinementMask2.hs b/src/MathObj/RefinementMask2.hs
new file mode 100644
--- /dev/null
+++ b/src/MathObj/RefinementMask2.hs
@@ -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]
+-}
diff --git a/src/MathObj/RootSet.hs b/src/MathObj/RootSet.hs
--- a/src/MathObj/RootSet.hs
+++ b/src/MathObj/RootSet.hs
@@ -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 -}
diff --git a/src/MyPrelude.hs b/src/MyPrelude.hs
deleted file mode 100644
--- a/src/MyPrelude.hs
+++ /dev/null
@@ -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)
diff --git a/src/Number/Complex.hs b/src/Number/Complex.hs
--- a/src/Number/Complex.hs
+++ b/src/Number/Complex.hs
@@ -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) -}
diff --git a/src/Number/DimensionTerm.hs b/src/Number/DimensionTerm.hs
--- a/src/Number/DimensionTerm.hs
+++ b/src/Number/DimensionTerm.hs
@@ -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
diff --git a/src/Number/DimensionTerm/SI.hs b/src/Number/DimensionTerm/SI.hs
--- a/src/Number/DimensionTerm/SI.hs
+++ b/src/Number/DimensionTerm/SI.hs
@@ -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
diff --git a/src/Number/FixedPoint.hs b/src/Number/FixedPoint.hs
--- a/src/Number/FixedPoint.hs
+++ b/src/Number/FixedPoint.hs
@@ -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)
 
diff --git a/src/Number/FixedPoint/Check.hs b/src/Number/FixedPoint/Check.hs
--- a/src/Number/FixedPoint/Check.hs
+++ b/src/Number/FixedPoint/Check.hs
@@ -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)
diff --git a/src/Number/GaloisField2p32m5.hs b/src/Number/GaloisField2p32m5.hs
--- a/src/Number/GaloisField2p32m5.hs
+++ b/src/Number/GaloisField2p32m5.hs
@@ -25,8 +25,8 @@
 
 import Test.QuickCheck (Arbitrary(arbitrary), )
 
-import PreludeBase
-import NumericPrelude
+import NumericPrelude.Base
+import NumericPrelude.Numeric
 
 
 newtype T = Cons {decons :: Word32}
diff --git a/src/Number/NonNegative.hs b/src/Number/NonNegative.hs
--- a/src/Number/NonNegative.hs
+++ b/src/Number/NonNegative.hs
@@ -1,4 +1,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
diff --git a/src/Number/NonNegativeChunky.hs b/src/Number/NonNegativeChunky.hs
--- a/src/Number/NonNegativeChunky.hs
+++ b/src/Number/NonNegativeChunky.hs
@@ -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
 
 
 
diff --git a/src/Number/OccasionallyScalarExpression.hs b/src/Number/OccasionallyScalarExpression.hs
--- a/src/Number/OccasionallyScalarExpression.hs
+++ b/src/Number/OccasionallyScalarExpression.hs
@@ -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
diff --git a/src/Number/PartiallyTranscendental.hs b/src/Number/PartiallyTranscendental.hs
--- a/src/Number/PartiallyTranscendental.hs
+++ b/src/Number/PartiallyTranscendental.hs
@@ -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
 
diff --git a/src/Number/Peano.hs b/src/Number/Peano.hs
--- a/src/Number/Peano.hs
+++ b/src/Number/Peano.hs
@@ -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
diff --git a/src/Number/Physical.hs b/src/Number/Physical.hs
--- a/src/Number/Physical.hs
+++ b/src/Number/Physical.hs
@@ -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
diff --git a/src/Number/Physical/Read.hs b/src/Number/Physical/Read.hs
--- a/src/Number/Physical/Read.hs
+++ b/src/Number/Physical/Read.hs
@@ -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
diff --git a/src/Number/Physical/Show.hs b/src/Number/Physical/Show.hs
--- a/src/Number/Physical/Show.hs
+++ b/src/Number/Physical/Show.hs
@@ -24,8 +24,8 @@
 import Data.List(find)
 import Data.Maybe(mapMaybe)
 
-import NumericPrelude
-import PreludeBase
+import NumericPrelude.Numeric
+import NumericPrelude.Base
 
 
 mulPrec :: Int
diff --git a/src/Number/Physical/Unit.hs b/src/Number/Physical/Unit.hs
--- a/src/Number/Physical/Unit.hs
+++ b/src/Number/Physical/Unit.hs
@@ -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
diff --git a/src/Number/Physical/UnitDatabase.hs b/src/Number/Physical/UnitDatabase.hs
--- a/src/Number/Physical/UnitDatabase.hs
+++ b/src/Number/Physical/UnitDatabase.hs
@@ -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]
 
diff --git a/src/Number/Positional.hs b/src/Number/Positional.hs
--- a/src/Number/Positional.hs
+++ b/src/Number/Positional.hs
@@ -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
diff --git a/src/Number/Positional/Check.hs b/src/Number/Positional/Check.hs
--- a/src/Number/Positional/Check.hs
+++ b/src/Number/Positional/Check.hs
@@ -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)
diff --git a/src/Number/Quaternion.hs b/src/Number/Quaternion.hs
--- a/src/Number/Quaternion.hs
+++ b/src/Number/Quaternion.hs
@@ -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, )
 
diff --git a/src/Number/Ratio.hs b/src/Number/Ratio.hs
--- a/src/Number/Ratio.hs
+++ b/src/Number/Ratio.hs
@@ -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) :%
diff --git a/src/Number/ResidueClass.hs b/src/Number/ResidueClass.hs
--- a/src/Number/ResidueClass.hs
+++ b/src/Number/ResidueClass.hs
@@ -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)
 
diff --git a/src/Number/ResidueClass/Check.hs b/src/Number/ResidueClass/Check.hs
--- a/src/Number/ResidueClass/Check.hs
+++ b/src/Number/ResidueClass/Check.hs
@@ -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, )
diff --git a/src/Number/ResidueClass/Func.hs b/src/Number/ResidueClass/Func.hs
--- a/src/Number/ResidueClass/Func.hs
+++ b/src/Number/ResidueClass/Func.hs
@@ -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
diff --git a/src/Number/ResidueClass/Maybe.hs b/src/Number/ResidueClass/Maybe.hs
--- a/src/Number/ResidueClass/Maybe.hs
+++ b/src/Number/ResidueClass/Maybe.hs
@@ -10,8 +10,8 @@
 
 import Algebra.ZeroTestable(isZero)
 
-import PreludeBase
-import NumericPrelude
+import NumericPrelude.Base
+import NumericPrelude.Numeric
 
 infix 7 /:, `Cons`
 
diff --git a/src/Number/ResidueClass/Reader.hs b/src/Number/ResidueClass/Reader.hs
--- a/src/Number/ResidueClass/Reader.hs
+++ b/src/Number/ResidueClass/Reader.hs
@@ -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}
 
diff --git a/src/Number/SI.hs b/src/Number/SI.hs
--- a/src/Number/SI.hs
+++ b/src/Number/SI.hs
@@ -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
 
diff --git a/src/Number/SI/Unit.hs b/src/Number/SI/Unit.hs
--- a/src/Number/SI/Unit.hs
+++ b/src/Number/SI/Unit.hs
@@ -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 |
diff --git a/src/NumericPrelude.hs b/src/NumericPrelude.hs
--- a/src/NumericPrelude.hs
+++ b/src/NumericPrelude.hs
@@ -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, )
diff --git a/src/NumericPrelude/Base.hs b/src/NumericPrelude/Base.hs
new file mode 100644
--- /dev/null
+++ b/src/NumericPrelude/Base.hs
@@ -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, )
diff --git a/src/NumericPrelude/Elementwise.hs b/src/NumericPrelude/Elementwise.hs
--- a/src/NumericPrelude/Elementwise.hs
+++ b/src/NumericPrelude/Elementwise.hs
@@ -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) =
diff --git a/src/NumericPrelude/List/Checked.hs b/src/NumericPrelude/List/Checked.hs
new file mode 100644
--- /dev/null
+++ b/src/NumericPrelude/List/Checked.hs
@@ -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
diff --git a/src/NumericPrelude/List/Generic.hs b/src/NumericPrelude/List/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/NumericPrelude/List/Generic.hs
@@ -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)
diff --git a/src/NumericPrelude/Numeric.hs b/src/NumericPrelude/Numeric.hs
new file mode 100644
--- /dev/null
+++ b/src/NumericPrelude/Numeric.hs
@@ -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)
diff --git a/src/PreludeBase.hs b/src/PreludeBase.hs
deleted file mode 100644
--- a/src/PreludeBase.hs
+++ /dev/null
@@ -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)
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -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 -}
diff --git a/test/Test/MathObj/Gaussian/Bell.hs b/test/Test/MathObj/Gaussian/Bell.hs
--- a/test/Test/MathObj/Gaussian/Bell.hs
+++ b/test/Test/MathObj/Gaussian/Bell.hs
@@ -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 ::
diff --git a/test/Test/MathObj/Gaussian/Polynomial.hs b/test/Test/MathObj/Gaussian/Polynomial.hs
--- a/test/Test/MathObj/Gaussian/Polynomial.hs
+++ b/test/Test/MathObj/Gaussian/Polynomial.hs
@@ -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 ::
diff --git a/test/Test/MathObj/Gaussian/Variance.hs b/test/Test/MathObj/Gaussian/Variance.hs
--- a/test/Test/MathObj/Gaussian/Variance.hs
+++ b/test/Test/MathObj/Gaussian/Variance.hs
@@ -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 ::
diff --git a/test/Test/MathObj/Matrix.hs b/test/Test/MathObj/Matrix.hs
--- a/test/Test/MathObj/Matrix.hs
+++ b/test/Test/MathObj/Matrix.hs
@@ -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
diff --git a/test/Test/MathObj/PartialFraction.hs b/test/Test/MathObj/PartialFraction.hs
--- a/test/Test/MathObj/PartialFraction.hs
+++ b/test/Test/MathObj/PartialFraction.hs
@@ -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 -}
diff --git a/test/Test/MathObj/Polynomial.hs b/test/Test/MathObj/Polynomial.hs
--- a/test/Test/MathObj/Polynomial.hs
+++ b/test/Test/MathObj/Polynomial.hs
@@ -1,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 ()
diff --git a/test/Test/MathObj/PowerSeries.hs b/test/Test/MathObj/PowerSeries.hs
--- a/test/Test/MathObj/PowerSeries.hs
+++ b/test/Test/MathObj/PowerSeries.hs
@@ -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 ::
diff --git a/test/Test/MathObj/RefinementMask2.hs b/test/Test/MathObj/RefinementMask2.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/MathObj/RefinementMask2.hs
@@ -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)) :
+      []
diff --git a/test/Test/Number/GaloisField2p32m5.hs b/test/Test/Number/GaloisField2p32m5.hs
--- a/test/Test/Number/GaloisField2p32m5.hs
+++ b/test/Test/Number/GaloisField2p32m5.hs
@@ -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 ()
diff --git a/test/Test/Run.hs b/test/Test/Run.hs
--- a/test/Test/Run.hs
+++ b/test/Test/Run.hs
@@ -1,5 +1,6 @@
 module Main where
 
+import qualified Test.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 :