Seeking More Advice on Library

[Matthew Donadio]

--Boundary_(ID_uwdC8KuP1Ecf2hP7+d2r1w)
Content-type: text/plain
Content-transfer-encoding: 7BIT

Hi all,

One of the needs for a few of my DSP routines is support for functions
over polynomials.  The attached file is what I threw together.

Since this definetly has a broader scope than just me, I wanted to get
some opionions on it.

1.  The module is pretty type-free right now.  I wanted to get some
opionions on the best type for polynomials and names for the type and
constructor(s).

The big question here (in my mind) is how to represent polynomials.  The
most obvious way is as a list of the coefficients.  They can also be
represented as a list of the roots.  Should there be a single type with
two constructors or two types each with one constructor.

2.  Are there any functions that are missing?  I have a root finder in a
different module (the one from Numeric Quest; Laguerre's method), but I
am going to rewrite it in terms of the module that I attached. 
Functions for inflation and deflation by a single root will need to be
added.

3.  Should I add infix operators?  I'm not sure what to do about this
because of namespace clutter.

Any other thoughts?

Once again, thanks.

-- 
Matthew Donadio <m.p.donadio@ieee.org>

--Boundary_(ID_uwdC8KuP1Ecf2hP7+d2r1w)
Content-type: text/x-haskell; name=Poly.hs; charset=ISO-8859-1
Content-transfer-encoding: 7BIT
Content-disposition: attachment; filename=Poly.hs

-----------------------------------------------------------------------------
-- |
-- Module      :  Poly.Poly
-- Copyright   :  (c) Matthew Donadio 2002
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Simple module for handling polynomials.
--
-----------------------------------------------------------------------------

-- TODO: We should really create a datatype for polynomials...

-- TODO: Should polydiv return the quotient and the remainder as a tuple?

module Poly.Poly where

-- * Types

-- | Polynomials are lists of numbers:
-- [ a0, a1, ... , an ] == an*x^n + ... + a1*x + a0
-- and negative exponents are currently verboten.

data Polynomial a 
    = P [a]  -- ^ builds a polynomial from a list
    deriving (Read,Show)

-- * Functions

-- | Evaluate a polynomial using Horner's method.

polyeval :: Num a => [a] -> a -> a
polyeval []     x = 0
polyeval (p:ps) x = p + x * polyeval ps x

-- | Add two polynomials

polyadd :: Num a => [a] -> [a] -> [a]
polyadd [] []          = []
polyadd [] ys          = ys
polyadd xs []          = xs
polyadd (x:xs) (y:ys)  = (x+y) : polyadd xs ys

-- | Subtract two polynomials

polysub :: Num a => [a] -> [a] -> [a]
polysub [] []          = []
polysub [] ys          = map negate ys
polysub xs []          = xs
polysub (x:xs) (y:ys)  = (x-y) : polysub xs ys

-- | Scale a polynomial

polyscale :: Num a => a -> [a] -> [a]
polyscale a x = map (a*) x

-- | Multiply two polynomials

polymult :: Num a => [a] -> [a] -> [a]
polymult (x:[]) ys = map (x*) ys
polymult (x:xs) ys = polyadd (map (x*) ys) (polymult xs (0:ys))

-- | Divide two polynomials

polydiv :: Fractional a => [a] -> [a] -> [a]
polydiv x y = reverse $ polydiv' (reverse x) (reverse y)
    where polydiv' (x:xs) y | length (x:xs) < length y = []
			    | otherwise = z : (polydiv' (tail (polysub (x:xs) (polymult [z] y))) y)
	      where z = x / head y

-- | Modulus of two polynomials (remainder of division)

polymod :: Fractional a => [a] -> [a] -> [a]
polymod x y = reverse $ polymod' (reverse x) (reverse y)
    where polymod' (x:xs) y | length (x:xs) < length y = (x:xs)
	                    | otherwise = polymod' (tail (polysub (x:xs) (polymult [z] y))) y
	      where z = x / head y

-- | Raise a polynomial to a non-negative integer power

polypow :: (Num a, Integral b) => [a] -> b -> [a]
polypow x 0 = [ 1 ]
polypow x 1 = x
polypow x 2 = polymult x x
polypow x n | even n = polymult x2 x2
	    | odd n  = polymult x (polymult x2 x2)
    where x2 = polypow x (n `div` 2)

-- | Polynomial substitution y(n) = x(w(n))

polysubst :: Num a => [a] -> [a] -> [a]
polysubst w x = foldr polyadd [0] (polysubst' 0 w x )
    where polysubst' _ _ []     = []
          polysubst' n w (x:xs) = map (x*) (polypow w n) : polysubst' (n+1) w xs

-- | Polynomial derivative

polyderiv :: Num a => [a] -> [a]
polyderiv (x:xs) = polyderiv' 1 xs
    where polyderiv' _ []     = []
          polyderiv' n (x:xs) = n * x : polyderiv' (n+1) xs

-- | Polynomial integration (C is set to 0)

polyinteg :: Fractional a => [a] -> [a]
polyinteg x = 0 : polyinteg' 1 x
    where polyinteg' _ []     = []
          polyinteg' n (x:xs) = x / fromIntegral n : polyinteg' (n+1) xs

-- | Convert roots to a polynomial

roots2poly :: Num a => [a] -> [a]
roots2poly (r:[]) = [-r, 1]
roots2poly (r:rs) = polymult [-r, 1] (roots2poly rs)

--Boundary_(ID_uwdC8KuP1Ecf2hP7+d2r1w)--