diff --git a/polynomial.cabal b/polynomial.cabal
--- a/polynomial.cabal
+++ b/polynomial.cabal
@@ -1,5 +1,5 @@
 name:                   polynomial
-version:                0.5
+version:                0.6
 stability:              provisional
 
 cabal-version:          >= 1.6
@@ -35,4 +35,4 @@
                         Math.Polynomial.Type
                         Math.Polynomial.Pretty
                         
-  build-depends:        base >= 3 && <5, pretty, prettyclass, vector-space
+  build-depends:        base >= 3 && <5, deepseq, pretty, prettyclass, vector-space
diff --git a/src/Math/Polynomial.hs b/src/Math/Polynomial.hs
--- a/src/Math/Polynomial.hs
+++ b/src/Math/Polynomial.hs
@@ -3,8 +3,9 @@
 module Math.Polynomial
     ( Endianness(..)
     , Poly, poly, polyCoeffs, polyIsZero, polyIsOne
-    , zero, one, x
+    , zero, one, constPoly, x
     , scalePoly, negatePoly
+    , composePoly
     , addPoly, sumPolys, multPoly, powPoly
     , quotRemPoly, quotPoly, remPoly
     , evalPoly, evalPolyDeriv, evalPolyDerivs
@@ -19,21 +20,40 @@
 import Data.List
 import Data.List.ZipSum
 
+-- |The polynomial \"0\"
 zero :: Num a => Poly a
 zero = poly LE []
 
+-- |The polynomial \"1\"
 one :: Num a => Poly a
-one = poly LE [1]
+one = constPoly 1
 
+-- |The polynomial (in x) \"x\"
 x :: Num a => Poly a
 x = poly LE [0,1]
 
+-- |Given some constant 'k', construct the polynomial whose value is 
+-- constantly 'k'.
+constPoly :: Num a => a -> Poly a
+constPoly x = poly LE [x]
+
+-- |Given some scalar 's' and a polynomial 'f', computes the polynomial 'g'
+-- such that:
+-- 
+-- > evalPoly g x = s * evalPoly f x
 scalePoly :: Num a => a -> Poly a -> Poly a
+scalePoly 0 _ = zero
 scalePoly s p = fmap (s*) p
 
+-- |Given some polynomial 'f', computes the polynomial 'g' such that:
+-- 
+-- > evalPoly g x = negate (evalPoly f x)
 negatePoly :: Num a => Poly a -> Poly a
 negatePoly = fmap negate
 
+-- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that:
+-- 
+-- > evalPoly h x = evalPoly f x + evalPoly g x
 addPoly :: Num a => Poly a -> Poly a -> Poly a
 addPoly (polyCoeffs LE ->  a) (polyCoeffs LE ->  b) = poly LE (zipSum a b)
 
@@ -44,24 +64,38 @@
 sumPolys [] = zero
 sumPolys ps = poly LE (foldl1 zipSum (map (polyCoeffs LE) ps))
 
+-- |Given polynomials 'f' and 'g', computes the polynomial 'h' such that:
+-- 
+-- > evalPoly h x = evalPoly f x * evalPoly g x
 multPoly :: Num a => Poly a -> Poly a -> Poly a
-multPoly (polyCoeffs LE -> xs) (polyCoeffs LE -> ys) = poly LE $ multX ys
+multPoly (polyCoeffs LE -> xs) (polyCoeffs LE -> ys) = poly LE (multPolyLE xs ys)
+
+-- |(Internal): multiply polynomials in LE order.  O(length xs * length ys).
+multPolyLE :: Num a => [a] -> [a] -> [a]
+multPolyLE _  []     = []
+multPolyLE xs (y:ys) = foldr mul [] xs
     where
-        multX (0:ys) = 0:multX ys
-        multX ys = foldl zipSum []
-            [ shift ++ map (x *) ys
-            | (x, shift) <- zip xs (inits (repeat 0))
-            , x /= 0
-            ]
+        mul 0 bs = 0 : bs
+        mul x bs = (x*y) : zipSum (map (x*) ys) bs
 
+-- |Given a polynomial 'f' and exponent 'n', computes the polynomial 'g'
+-- such that:
+-- 
+-- > evalPoly g x = evalPoly f x ^ n
 powPoly :: (Num a, Integral b) => Poly a -> b -> Poly a
 powPoly _ 0 = poly LE [1]
 powPoly p 1 = p
 powPoly p n
+    | n < 0     = error "powPoly: negative exponent"
     | odd n     = p `multPoly` powPoly p (n-1)
     | otherwise = (\x -> multPoly x x) (powPoly p (n`div`2))
 
+-- |Given polynomials @a@ and @b@, with @b@ not 'zero', computes polynomials
+-- @q@ and @r@ such that:
+-- 
+-- > addPoly (multPoly q b) r == a
 quotRemPoly :: Fractional a => Poly a -> Poly a -> (Poly a, Poly a)
+quotRemPoly _ b | polyIsZero b = error "quotRemPoly: divide by zero"
 quotRemPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)
     = go [] u (length u - length v)
     where
@@ -75,8 +109,11 @@
                 u' = tail (zipSum u (map (negate q0 *) v))
 
 quotPoly :: Fractional a => Poly a -> Poly a -> Poly a
-quotPoly u v = fst (quotRemPoly u v)
+quotPoly u v
+    | polyIsZero v  = error "quotPoly: divide by zero"
+    | otherwise     = fst (quotRemPoly u v)
 remPoly :: Fractional a => Poly a -> Poly a -> Poly a
+remPoly _ b | polyIsZero b = error "remPoly: divide by zero"
 remPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)
     = go u (length u - length v)
     where
@@ -89,17 +126,59 @@
                 q0 = head u / v0
                 u' = tail (zipSum u (map (negate q0 *) v))
 
+-- |@composePoly f g@ constructs the polynomial 'h' such that:
+-- 
+-- > evalPoly h = evalPoly f . evalPoly g
+-- 
+-- This is a very expensive operation and, in general, returns a polynomial 
+-- that is quite a bit more expensive to evaluate than @f@ and @g@ together
+-- (because it is of a much higher order than either).  Unless your 
+-- polynomials are quite small or you are quite certain you need the
+-- coefficients of the composed polynomial, it is recommended that you 
+-- simply evaluate @f@ and @g@ and explicitly compose the resulting 
+-- functions.  This will usually be much more efficient.
+composePoly :: Num a => Poly a -> Poly a -> Poly a
+composePoly (polyCoeffs LE -> cs) (polyCoeffs LE -> ds) = poly LE (foldr mul [] cs)
+    where
+        -- Implementation note: this is a hand-inlining of the following
+        -- (with the 'Num' instance in "Math.Polynomial.NumInstance"):
+        -- > composePoly f g = evalPoly (fmap constPoly f) g
+        -- 
+        -- This is a very expensive operation, something like
+        -- O(length cs ^ 2 * length ds) I believe. There may be some more 
+        -- tricks to improve that, but I suspect there isn't much room for 
+        -- improvement. The number of terms in the resulting polynomial is 
+        -- O(length cs * length ds) already, and each one is the sum of 
+        -- quite a few terms.
+        mul c acc = addScalarLE c (multPolyLE acc ds)
 
+-- |(internal) add a scalar to a list of polynomial coefficients in LE order
+addScalarLE :: Num a => a -> [a] -> [a]
+addScalarLE 0 bs = bs
+addScalarLE a [] = [a]
+addScalarLE a (b:bs) = (a + b) : bs
+
+-- |Evaluate a polynomial at a point or, equivalently, convert a polynomial
+-- to the function it represents.  For example, @evalPoly 'x' = 'id'@ and 
+-- @evalPoly ('constPoly' k) = 'const' k.@
 evalPoly :: Num a => Poly a -> a -> a
+evalPoly (polyCoeffs LE -> cs) 0
+    | null cs   = 0
+    | otherwise = head cs
 evalPoly (polyCoeffs LE -> cs) x = foldr mul 0 cs
     where
         mul c acc = c + acc * x
 
+-- |Evaluate a polynomial and its derivative (respectively) at a point.
 evalPolyDeriv :: Num a => Poly a -> a -> (a,a)
 evalPolyDeriv (polyCoeffs LE -> cs) x = foldr mul (0,0) cs
     where
         mul c (p, dp) = (p * x + c, dp * x + p)
 
+-- |Evaluate a polynomial and all of its nonzero derivatives at a point.  
+-- This is roughly equivalent to:
+-- 
+-- > evalPolyDerivs p x = map (`evalPoly` x) (takeWhile (not . polyIsZero) (iterate polyDeriv p))
 evalPolyDerivs :: Num a => Poly a -> a -> [a]
 evalPolyDerivs (polyCoeffs LE -> cs) x = trunc . zipWith (*) factorials $ foldr mul [] cs
     where
@@ -117,6 +196,9 @@
         cut remainder swap = (swap + remainder * a, remainder)
         (r,q) = mapAccumR cut 0 cs
 
+-- |@gcdPoly a b@ computes the highest order monic polynomial that is a 
+-- divisor of both @a@ and @b@.  If both @a@ and @b@ are 'zero', the 
+-- result is undefined.
 gcdPoly :: Fractional a => Poly a -> Poly a -> Poly a
 gcdPoly a b 
     | polyIsZero b  = if polyIsZero a
@@ -124,13 +206,13 @@
         else monic a
     | otherwise     = gcdPoly b (a `remPoly` b)
 
--- |Normalize a polynomial so that its highest-order coefficient is 1
+-- |(internal) Normalize a polynomial so that its highest-order coefficient is 1
 monic :: Fractional a => Poly a -> Poly a
 monic p = case polyCoeffs BE p of
     []      -> poly BE []
     (c:cs)  -> poly BE (1:map (/c) cs)
 
-
+-- |Compute the derivative of a polynomial.
 polyDeriv :: Num a => Poly a -> Poly a
 polyDeriv (polyCoeffs LE -> cs) = poly LE
     [ c * n
@@ -138,6 +220,7 @@
     | n <- iterate (1+) 1
     ]
 
+-- |Compute the definite integral (from 0 to x) of a polynomial.
 polyIntegral :: Fractional a => Poly a -> Poly a
 polyIntegral (polyCoeffs LE -> cs) = poly LE $ 0 :
     [ c / n
@@ -145,7 +228,7 @@
     | n <- iterate (1+) 1
     ]
 
--- |Separate a polynomial into a set of factors none of which have
+-- |Separate a nonzero polynomial into a set of factors none of which have
 -- multiple roots, and the product of which is the original polynomial.
 -- Note that if division is not exact, it may fail to separate roots.  
 -- Rational coefficients is a good idea.
@@ -153,6 +236,7 @@
 -- Useful when applicable as a way to simplify root-finding problems.
 separateRoots :: Fractional a => Poly a -> [Poly a]
 separateRoots p
+    | polyIsZero q  = error "separateRoots: zero polynomial"
     | polyIsOne q   = [p]
     | otherwise     = p `quotPoly` q : separateRoots q
     where
diff --git a/src/Math/Polynomial/Bernstein.hs b/src/Math/Polynomial/Bernstein.hs
--- a/src/Math/Polynomial/Bernstein.hs
+++ b/src/Math/Polynomial/Bernstein.hs
@@ -58,14 +58,9 @@
 -- |de Casteljau's algorithm, returning the whole tableau.  Used both for
 -- evaluating and splitting polynomials in Bernstein form.
 deCasteljau :: Num a => [a] -> a -> [[a]]
-deCasteljau cs t = takeWhile (not.null) table
-    where
-        table = cs : 
-            [ [ b_i * (1-t) + b_ip1 * t
-              | b_i:b_ip1:_ <- tails row
-              ]
-            | row <- table
-            ]
+deCasteljau [] _ = []
+deCasteljau cs t = cs : deCasteljau (zipWith (interp t) cs (tail cs)) t
+    where interp t x0 x1 = (1-t)*x0 + t*x1
 
 -- |Given a polynomial in Bernstein form (that is, a list of coefficients
 -- for a basis set from 'bernstein', such as is returned by 'bernsteinFit')
diff --git a/src/Math/Polynomial/Chebyshev.hs b/src/Math/Polynomial/Chebyshev.hs
--- a/src/Math/Polynomial/Chebyshev.hs
+++ b/src/Math/Polynomial/Chebyshev.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE ParallelListComp #-}
+{-# LANGUAGE ParallelListComp, BangPatterns #-}
 module Math.Polynomial.Chebyshev where
 
 import Math.Polynomial
@@ -7,7 +7,7 @@
 -- |The Chebyshev polynomials of the first kind with 'Integer' coefficients.
 ts :: [Poly Integer]
 ts = poly LE [1] : 
-    [ addPoly (poly LE [0, 1]    `multPoly` t_n)
+    [ addPoly (x                `multPoly` t_n)
               (poly LE [-1,0,1] `multPoly` u_n)
     | t_n <- ts
     | u_n <- poly LE [0] : us
@@ -16,24 +16,26 @@
 -- The Chebyshev polynomials of the second kind with 'Integer' coefficients.
 us :: [Poly Integer]
 us = 
-    [ addPoly t_n (multPoly u_n (poly LE [0,1]))
+    [ addPoly t_n (multPoly x u_n)
     | t_n <- ts
     | u_n <- poly LE [0] : us
     ]
 
 -- |Compute the coefficients of the n'th Chebyshev polynomial of the first kind.
 t :: Num a => Int -> Poly a
-t n = poly LE . map fromInteger . polyCoeffs LE $ ts !! n
+t n | n >= 0    = poly LE . map fromInteger . polyCoeffs LE $ ts !! n
+    | otherwise = error "t: negative index"
 
 -- |Compute the coefficients of the n'th Chebyshev polynomial of the second kind.
 u :: Num a => Int -> Poly a
-u n = poly LE . map fromInteger . polyCoeffs LE $ us !! n
+u n | n >= 0    = poly LE . map fromInteger . polyCoeffs LE $ us !! n
+    | otherwise = error "u: negative index"
 
 -- |Evaluate the n'th Chebyshev polynomial of the first kind at a point X.  
 -- Both more efficient and more numerically stable than computing the 
 -- coefficients and evaluating the polynomial.
 evalT :: Num a => Int -> a -> a
-evalT n x = evalTs x !! n
+evalT n x = fst (evalTU n x)
 
 -- |Evaluate all the Chebyshev polynomials of the first kind at a point X.
 evalTs :: Num a => a -> [a]
@@ -43,7 +45,7 @@
 -- Both more efficient and more numerically stable than computing the 
 -- coefficients and evaluating the polynomial.
 evalU :: Num a => Int -> a -> a
-evalU n x = evalUs x !! n
+evalU n x = snd (evalTU n x)
 
 -- |Evaluate all the Chebyshev polynomials of the second kind at a point X.
 evalUs :: Num a => a -> [a]
@@ -51,10 +53,15 @@
 
 -- |Evaluate the n'th Chebyshev polynomials of both kinds at a point X.
 evalTU :: Num a => Int -> a -> (a,a)
-evalTU n x = (ts!!n, us!!n)
-    where (ts,us) = evalTsUs x
+evalTU n x = go n 1 0
+    where
+        go !0 !t_n !u_n = (t_n, u_n)
+        go !n !t_n !u_n = go (n-1) t_np1 u_np1
+            where
+                t_np1 = x * t_n - (1-x*x)*u_n
+                u_np1 = x * u_n + t_n
 
--- |Evaluate all the Chebyshev polynomials of the both kinds at a point X.
+-- |Evaluate all the Chebyshev polynomials of both kinds at a point X.
 evalTsUs :: Num a => a -> ([a], [a])
 evalTsUs x = (ts, tail us)
     where
diff --git a/src/Math/Polynomial/Type.hs b/src/Math/Polynomial/Type.hs
--- a/src/Math/Polynomial/Type.hs
+++ b/src/Math/Polynomial/Type.hs
@@ -1,10 +1,13 @@
 {-# LANGUAGE ViewPatterns, TypeFamilies #-}
+-- |Low-level interface for the 'Poly' type.
 module Math.Polynomial.Type 
     ( Endianness(..)
     , Poly, poly, polyCoeffs
+    , trim, rawPoly, rawPolyCoeffs
     , polyIsZero, polyIsOne
     ) where
 
+import Control.DeepSeq
 -- import Data.List.Extras.LazyLength
 import Data.AdditiveGroup
 import Data.VectorSpace
@@ -22,26 +25,47 @@
         -- at end of string discard the stash
         go _ [] = []
 
-trim :: Num a => Poly a -> Poly a
-trim p@(Poly _ True _) = p
-trim   (Poly LE _ cs) = Poly LE True (dropEnd   (==0) cs)
-trim   (Poly BE _ cs) = Poly BE True (dropWhile (==0) cs)
+-- |Trim zeroes from a polynomial (given a predicate for identifying zero).
+-- In particular, drops zeroes from the highest-order coefficients, so that
+-- @0x^n + 0x^(n-1) + 0x^(n-2) + ... + ax^k + ...@, @a /= 0@
+-- is normalized to @ax^k + ...@.  
+-- 
+-- The 'Eq' instance for 'Poly' and all the standard constructors / destructors
+-- are defined using @trim (0==)@.
+trim :: (a -> Bool) -> Poly a -> Poly a
+trim      _ p@(Poly _ True _) = p
+trim isZero   (Poly LE _ cs) = Poly LE True (dropEnd   isZero cs)
+trim isZero   (Poly BE _ cs) = Poly BE True (dropWhile isZero cs)
 
 -- |Make a 'Poly' from a list of coefficients using the specified coefficient order.
 poly :: Num a => Endianness -> [a] -> Poly a
-poly end cs = trim (Poly end False cs)
+poly end cs = trim (0==) (rawPoly end cs)
 
+-- |Make a 'Poly' from a list of coefficients using the specified coefficient order,
+-- without the 'Num' context (and therefore without trimming zeroes from the 
+-- coefficient list)
+rawPoly :: Endianness -> [a] -> Poly a
+rawPoly end cs = Poly end False cs 
+
 -- |Get the coefficients of a a 'Poly' in the specified order.
 polyCoeffs :: Num a => Endianness -> Poly a -> [a]
-polyCoeffs end p = case trim p of
-    Poly e _ cs | e == end  -> cs
-                | otherwise -> reverse cs
+polyCoeffs end p = rawPolyCoeffs end (trim (0==) p)
 
+-- |Get the coefficients of a a 'Poly' in the specified order, without the 'Num'
+-- constraint (and therefore without trimming zeroes).
+-- 
+-- This function does not respect the 'Eq' instance:
+--   @x == y@ =/=> @rawPolyCoeffs e x == rawPolyCoeffs e y@.
+rawPolyCoeffs :: Endianness -> Poly a -> [a]
+rawPolyCoeffs end (Poly e _ cs)
+    | e == end  = cs
+    | otherwise = reverse cs
+
 polyIsZero :: Num a => Poly a -> Bool
-polyIsZero = null . coeffs . trim
+polyIsZero = null . coeffs . trim (0==)
 
 polyIsOne :: Num a => Poly a -> Bool
-polyIsOne = ([1]==) . coeffs . trim
+polyIsOne = ([1]==) . coeffs . trim (0==)
 
 data Endianness 
     = BE 
@@ -50,14 +74,20 @@
     -- ^ Little-Endian (head is const term)
     deriving (Eq, Ord, Enum, Bounded, Show)
 
+instance NFData Endianness where
+    rnf x = seq x ()
+
 data Poly a = Poly 
     { endianness :: !Endianness
     , _trimmed   :: !Bool
     , coeffs     :: ![a]
     }
 
+instance NFData a => NFData (Poly a) where
+    rnf (Poly e t c) = rnf e `seq` rnf t `seq` rnf c
+
 instance Num a => Show (Poly a) where
-    showsPrec p (trim -> Poly end _ cs) 
+    showsPrec p (trim (0==) -> Poly end _ cs) 
         = showParen (p > 10) 
             ( showString "poly "
             . showsPrec 11 end
@@ -68,7 +98,7 @@
 instance (Num a, Eq a) => Eq (Poly a) where
     p == q  
         | endianness p == endianness q
-        = coeffs (trim p) == coeffs (trim q)
+        = coeffs (trim (0==) p) == coeffs (trim (0==) q)
         | otherwise 
         = polyCoeffs BE p == polyCoeffs BE q
         
@@ -88,16 +118,10 @@
 instance Functor Poly where
     fmap f (Poly end _ cs) = Poly end False (map f cs)
 
-
--- Local-use-only: extract coefficients in LE order, without Num constraint
--- (and therefore without trimming)
-le :: Poly a -> [a]
-le p@(endianness -> LE) = coeffs p
-le p                    = reverse (coeffs p)
-
 instance AdditiveGroup a => AdditiveGroup (Poly a) where
     zeroV = Poly LE True []
-    (le ->  a) ^+^ (le ->  b) = Poly LE False (zipSumV a b)
+    (rawPolyCoeffs LE ->  a) ^+^ (rawPolyCoeffs LE ->  b) 
+        = Poly LE False (zipSumV a b)
     negateV = fmap negateV
 
 instance VectorSpace a => VectorSpace (Poly a) where