packages feed

polynomial (empty) → 0.5

raw patch · 13 files changed

+848/−0 lines, 13 filesdep +basedep +prettydep +prettyclasssetup-changed

Dependencies added: base, pretty, prettyclass, vector-space

Files

+ Setup.lhs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain+
+ polynomial.cabal view
@@ -0,0 +1,38 @@+name:                   polynomial+version:                0.5+stability:              provisional++cabal-version:          >= 1.6+build-type:             Simple++author:                 James Cook <mokus@deepbondi.net>+maintainer:             James Cook <mokus@deepbondi.net>+license:                PublicDomain+homepage:               /dev/null++category:               Math, Numerical+synopsis:               Polynomials+description:            A type for representing polynomials, several functions+                        for manipulating and evaluating them, and several+                        interesting polynomial sequences.++source-repository head+  type: darcs+  location: http://code.haskell.org/~mokus/polynomial++Library+  ghc-options:          -Wall -fno-warn-name-shadowing+  hs-source-dirs:       src+  exposed-modules:      Math.Polynomial+                        Math.Polynomial.Bernstein+                        Math.Polynomial.Chebyshev+                        Math.Polynomial.Interpolation+                        Math.Polynomial.Lagrange+                        Math.Polynomial.Legendre+                        Math.Polynomial.Newton+                        Math.Polynomial.NumInstance+  other-modules:        Data.List.ZipSum+                        Math.Polynomial.Type+                        Math.Polynomial.Pretty+                        +  build-depends:        base >= 3 && <5, pretty, prettyclass, vector-space
+ src/Data/List/ZipSum.hs view
@@ -0,0 +1,18 @@+module Data.List.ZipSum where++import Data.AdditiveGroup++-- like @zipWith (+)@ except that when the end of either list is+-- reached, the rest of the output is the rest of the longer input list.+zipSum :: Num t => [t] -> [t] -> [t]+zipSum xs [] = xs+zipSum [] ys = ys+zipSum (x:xs) (y:ys) = (x+y) : zipSum xs ys++-- like @zipWith (^+^)@ except that when the end of either list is+-- reached, the rest of the output is the rest of the longer input list.+zipSumV :: AdditiveGroup t => [t] -> [t] -> [t]+zipSumV xs [] = xs+zipSumV [] ys = ys+zipSumV (x:xs) (y:ys) = (x^+^y) : zipSumV xs ys+
+ src/Math/Polynomial.hs view
@@ -0,0 +1,159 @@+{-# LANGUAGE ParallelListComp, ViewPatterns, FlexibleContexts #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}+module Math.Polynomial+    ( Endianness(..)+    , Poly, poly, polyCoeffs, polyIsZero, polyIsOne+    , zero, one, x+    , scalePoly, negatePoly+    , addPoly, sumPolys, multPoly, powPoly+    , quotRemPoly, quotPoly, remPoly+    , evalPoly, evalPolyDeriv, evalPolyDerivs+    , contractPoly+    , gcdPoly, separateRoots+    , polyDeriv, polyIntegral+    ) where++import Math.Polynomial.Type+import Math.Polynomial.Pretty ({- instance -})++import Data.List+import Data.List.ZipSum++zero :: Num a => Poly a+zero = poly LE []++one :: Num a => Poly a+one = poly LE [1]++x :: Num a => Poly a+x = poly LE [0,1]++scalePoly :: Num a => a -> Poly a -> Poly a+scalePoly s p = fmap (s*) p++negatePoly :: Num a => Poly a -> Poly a+negatePoly = fmap negate++addPoly :: Num a => Poly a -> Poly a -> Poly a+addPoly (polyCoeffs LE ->  a) (polyCoeffs LE ->  b) = poly LE (zipSum a b)++{-# RULES+  "sum Poly"    forall ps. foldl addPoly zero ps = sumPolys ps+  #-}+sumPolys :: Num a => [Poly a] -> Poly a+sumPolys [] = zero+sumPolys ps = poly LE (foldl1 zipSum (map (polyCoeffs LE) ps))++multPoly :: Num a => Poly a -> Poly a -> Poly a+multPoly (polyCoeffs LE -> xs) (polyCoeffs LE -> ys) = poly LE $ multX ys+    where+        multX (0:ys) = 0:multX ys+        multX ys = foldl zipSum []+            [ shift ++ map (x *) ys+            | (x, shift) <- zip xs (inits (repeat 0))+            , x /= 0+            ]++powPoly :: (Num a, Integral b) => Poly a -> b -> Poly a+powPoly _ 0 = poly LE [1]+powPoly p 1 = p+powPoly p n+    | odd n     = p `multPoly` powPoly p (n-1)+    | otherwise = (\x -> multPoly x x) (powPoly p (n`div`2))++quotRemPoly :: Fractional a => Poly a -> Poly a -> (Poly a, Poly a)+quotRemPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)+    = go [] u (length u - length v)+    where+        v0  | null v    = 0+            | otherwise = head v+        go q u n+            | null u || n < 0   = (poly LE q, poly BE u)+            | otherwise         = go (q0:q) u' (n-1)+            where+                q0 = head u / v0+                u' = tail (zipSum u (map (negate q0 *) v))++quotPoly :: Fractional a => Poly a -> Poly a -> Poly a+quotPoly u v = fst (quotRemPoly u v)+remPoly :: Fractional a => Poly a -> Poly a -> Poly a+remPoly (polyCoeffs BE -> u) (polyCoeffs BE -> v)+    = go u (length u - length v)+    where+        v0  | null v    = 0+            | otherwise = head v+        go u n+            | null u || n < 0   = poly BE u+            | otherwise         = go u' (n-1)+            where+                q0 = head u / v0+                u' = tail (zipSum u (map (negate q0 *) v))+++evalPoly :: Num a => Poly a -> a -> a+evalPoly (polyCoeffs LE -> cs) x = foldr mul 0 cs+    where+        mul c acc = c + acc * x++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)++evalPolyDerivs :: Num a => Poly a -> a -> [a]+evalPolyDerivs (polyCoeffs LE -> cs) x = trunc . zipWith (*) factorials $ foldr mul [] cs+    where+        trunc list = zipWith const list cs+        factorials = scanl (*) 1 (iterate (+1) 1)+        mul c pds@(p:pd) = (p * x + c) : map (x *) pd `zipSum` pds+        mul c [] = [c]++-- |\"Contract\" a polynomial by attempting to divide out a root.+--+-- @contractPoly p a@ returns @(q,r)@ such that @q*(x-a) + r == p@+contractPoly :: Num a => Poly a -> a -> (Poly a, a)+contractPoly (polyCoeffs LE -> cs) a = (poly LE q, r)+    where+        cut remainder swap = (swap + remainder * a, remainder)+        (r,q) = mapAccumR cut 0 cs++gcdPoly :: Fractional a => Poly a -> Poly a -> Poly a+gcdPoly a b +    | polyIsZero b  = if polyIsZero a+        then error "gcdPoly: gcdPoly zero zero is undefined"+        else monic a+    | otherwise     = gcdPoly b (a `remPoly` b)++-- |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)+++polyDeriv :: Num a => Poly a -> Poly a+polyDeriv (polyCoeffs LE -> cs) = poly LE+    [ c * n+    | c <- drop 1 cs+    | n <- iterate (1+) 1+    ]++polyIntegral :: Fractional a => Poly a -> Poly a+polyIntegral (polyCoeffs LE -> cs) = poly LE $ 0 :+    [ c / n+    | c <- cs+    | n <- iterate (1+) 1+    ]++-- |Separate a 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.+--+-- Useful when applicable as a way to simplify root-finding problems.+separateRoots :: Fractional a => Poly a -> [Poly a]+separateRoots p+    | polyIsOne q   = [p]+    | otherwise     = p `quotPoly` q : separateRoots q+    where+        q = gcdPoly p (polyDeriv p)
+ src/Math/Polynomial/Bernstein.hs view
@@ -0,0 +1,80 @@+{-# LANGUAGE ParallelListComp #-}+module Math.Polynomial.Bernstein+    ( bernstein+    , evalBernstein+    , bernsteinFit+    , evalBernsteinSeries+    , deCasteljau+    , splitBernsteinSeries+    ) where++import Math.Polynomial+import Data.List++-- |The Bernstein basis polynomials.  The @n@th inner list is a basis for +-- the polynomials of order @n@ or lower.  The @n@th basis consists of @n@+-- polynomials of order @n@ which sum to @1@, and have roots of varying +-- multiplicities at @0@ and @1@.+bernstein :: [[Poly Integer]]+bernstein = +    [ [ scalePoly nCv p `multPoly` q+      | q <- reverse qs+      | p <- ps+      | nCv  <- bico+      ]+    | ps <- tail $ inits [poly BE (1 : zs) | zs <- inits (repeat 0)]+    | qs <- tail $ inits (iterate (multPoly (poly LE [1,-1])) one)+    | bico <- ptri+    ]+    where+        -- pascal's triangle+        ptri = [1] : [ 1 : zipWith (+) row (tail row) ++ [1] | row <- ptri]++-- |@evalBernstein n v x@ evaluates the @v@'th Bernstein polynomial of order @n@+-- at the point @x@.+evalBernstein :: (Integral a, Num b) => a -> a -> b -> b+evalBernstein n v t+    | n < 0 || v > n    = 0+    | otherwise         = fromInteger nCv * t^v * (1-t)^(n-v)+    where+        n' = toInteger n+        v' = toInteger v+        nCv = product [1..n'] `div` (product [1..v'] * product [1..n'-v'])++-- |@bernsteinFit n f@: Approximate a function @f@ as a linear combination of+-- Bernstein polynomials of order @n@.  This approximation converges slowly+-- but uniformly to @f@ on the interval [0,1].+bernsteinFit :: (Fractional b, Integral a) => a -> (b -> b) -> [b]+bernsteinFit n f = [f (fromIntegral v / fromIntegral n) | v <- [0..n]]++-- |Evaluate a polynomial given as a list of @n@ coefficients for the @n@th+-- Bernstein basis.  Roughly:+-- +-- > evalBernsteinSeries cs = sum (zipWith scalePoly cs (bernstein !! (length cs - 1)))+evalBernsteinSeries :: Num a => [a] -> a -> a+evalBernsteinSeries [] = const 0+evalBernsteinSeries cs = head . last . deCasteljau cs++-- |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+            ]++-- |Given a polynomial in Bernstein form (that is, a list of coefficients+-- for a basis set from 'bernstein', such as is returned by 'bernsteinFit')+-- and a parameter value @x@, split the polynomial into two halves, mapping+-- @[0,x]@ and @[x,1]@ respectively onto @[0,1]@.+--+-- A typical use for this operation would be to split a Bezier curve +-- (inserting a new knot at @x@).+splitBernsteinSeries :: Num a => [a] -> a -> ([a], [a])+splitBernsteinSeries cs t = (map head betas, map last (reverse betas))+    where+        betas = deCasteljau cs t
+ src/Math/Polynomial/Chebyshev.hs view
@@ -0,0 +1,109 @@+{-# LANGUAGE ParallelListComp #-}+module Math.Polynomial.Chebyshev where++import Math.Polynomial+import Data.List++-- |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)+              (poly LE [-1,0,1] `multPoly` u_n)+    | t_n <- ts+    | u_n <- poly LE [0] : us+    ]++-- The Chebyshev polynomials of the second kind with 'Integer' coefficients.+us :: [Poly Integer]+us = +    [ addPoly t_n (multPoly u_n (poly LE [0,1]))+    | 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++-- |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++-- |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++-- |Evaluate all the Chebyshev polynomials of the first kind at a point X.+evalTs :: Num a => a -> [a]+evalTs = fst . evalTsUs++-- |Evaluate the n'th Chebyshev polynomial of the second kind at a point X.  +-- 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++-- |Evaluate all the Chebyshev polynomials of the second kind at a point X.+evalUs :: Num a => a -> [a]+evalUs = snd . evalTsUs++-- |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++-- |Evaluate all the Chebyshev polynomials of the both kinds at a point X.+evalTsUs :: Num a => a -> ([a], [a])+evalTsUs x = (ts, tail us)+    where+        ts = 1 : [x * t_n - (1-x*x)*u_n  | t_n <- ts | u_n <- us]+        us = 0 : [x * u_n + t_n          | t_n <- ts | u_n <- us]++-- |Compute the roots of the n'th Chebyshev polynomial of the first kind.+tRoots :: Floating a => Int -> [a]+tRoots   n = [cos (pi / fromIntegral n * (fromIntegral k + 0.5)) | k <- [0..n-1]]++-- |Compute the extreme points of the n'th Chebyshev polynomial of the first kind.+tExtrema :: Floating a => Int -> [a]+tExtrema n = [cos (pi / fromIntegral n *  fromIntegral k       ) | k <- [0..n]]++-- |@chebyshevFit n f@ returns a list of N coefficients @cs@ such that +-- @f x@ ~= @sum (zipWith (*) cs (evalTs x))@ on the interval -1 < x < 1.+-- +-- The N roots of the N'th Chebyshev polynomial are the fitting points at +-- which the function will be evaluated and at which the approximation will be+-- exact.  These points always lie within the interval -1 < x < 1.  Outside+-- this interval, the approximation will diverge quickly.+--+-- This function deviates from most chebyshev-fit implementations in that it+-- returns the first coefficient pre-scaled so that the series evaluation +-- operation is a simple inner product, since in most other algorithms+-- operating on chebyshev series, that factor is almost always a nuissance.+chebyshevFit :: Floating a => Int -> (a -> a) -> [a]+chebyshevFit n f = +    [ oneOrTwo / fromIntegral n +    * sum (zipWith (*) ts fxs)+    | ts <- transpose txs+    | oneOrTwo <- 1 : repeat 2+    ]+    where+        txs = map (take n . evalTs) xs+        fxs = map f xs+        xs = tRoots n++-- |Evaluate a Chebyshev series expansion with a finite number of terms.+-- +-- Note that this function expects the first coefficient to be pre-scaled+-- by 1/2, which is what is produced by 'chebyshevFit'.  Thus, this computes+-- a simple inner product of the given list with a matching-length sequence of+-- chebyshev polynomials.+evalChebyshevSeries :: Num a => [a] -> a -> a+evalChebyshevSeries     []  _ = 0+evalChebyshevSeries (c0:cs) x = +        let b1:b2:_ = reverse bs+         in x*b1 - b2 + c0+    where+        -- Clenshaw's recurrence formula+        bs = 0 : 0 : [2*x*b1 - b2 + c | b2:b1:_ <- tails bs | c <- reverse cs]
+ src/Math/Polynomial/Interpolation.hs view
@@ -0,0 +1,97 @@+{-# LANGUAGE ParallelListComp #-}+module Math.Polynomial.Interpolation where++import Math.Polynomial+import Math.Polynomial.Lagrange+import Data.List++-- |Evaluate a polynomial passing through the specified set of points.  The+-- order of the interpolating polynomial will (at most) be one less than +-- the number of points given.+polyInterp :: Fractional a => [(a,a)] -> a -> a+polyInterp xys = head . last . neville xys++-- |Computes the tableau generated by Neville's algorithm.  Each successive+-- row of the table is a list of interpolants one order higher than the previous,+-- using a range of input points starting at the same position in the input+-- list as the interpolant's position in the output list.+neville :: Fractional a => [(a,a)] -> a -> [[a]]+neville xys x = table+    where+        (xs,ys) = unzip xys+        table = ys :+            [ [ ((x - x_j) * p1 + (x_i - x) * p0) / (x_i - x_j)+              | p0:p1:_ <- tails row+              | x_j     <- xs+              | x_i     <- x_is+              ]+            | row  <- table+            | x_is <- tail (tails xs)+            , not (null x_is)+            ]++-- |Computes the tableau generated by a modified form of Neville's algorithm+-- described in Numerical Recipes, Ch. 3, Sec. 2, which records the differences+-- between interpolants at each level.  Each pair (c,d) is the amount to add+-- to the previous level's interpolant at either the same or the subsequent+-- position (respectively) in order to obtain the new level's interpolant.+-- Mathematically, either sum yields the same value, but due to numerical+-- errors they may differ slightly, and some \"paths\" through the table+-- may yield more accurate final results than others.+nevilleDiffs :: Fractional a => [(a,a)] -> a -> [[(a,a)]]+nevilleDiffs xys x = table+    where+        (xs,ys) = unzip xys+        table = zip ys ys :+            [ [ ( {-c-} (x_j - x) * (c1 - d0) / (x_j - x_i)+                , {-d-} (x_i - x) * (c1 - d0) / (x_j - x_i)+                )+              | (_c0,d0):(c1,_d1):_ <- tails row+              | x_j     <- xs+              | x_i     <- x_is+              ]+            | row  <- table+            | x_is <- tail (tails xs)+            , not (null x_is)+            ]++-- |Fit a polynomial to a set of points by iteratively evaluating the +-- interpolated polynomial (using 'polyInterp') at 0 to establish the+-- constant coefficient and reducing the polynomial by subtracting that+-- coefficient from all y's and dividing by their corresponding x's.+-- +-- Slower than 'lagrangePolyFit' but stable under different sets of +-- conditions.+-- +-- Note that computing the coefficients of a fitting polynomial is an +-- inherently ill-conditioned problem.  In most cases it is both faster and +-- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating+-- a fitted polynomial.+iterativePolyFit :: Fractional a => [(a,a)] -> Poly a+iterativePolyFit = poly LE . loop+    where+        loop  [] = []+        loop xys = c0 : loop (drop 1 xys')+            where+                c0   = polyInterp xys 0+                xys' = +                    [ (x,(y - c0) / x)+                    | (x,y) <- xys+                    ]++-- |Fit a polynomial to a set of points using barycentric Lagrange polynomials.+-- +-- Note that computing the coefficients of a fitting polynomial is an +-- inherently ill-conditioned problem.  In most cases it is both faster and +-- more accurate to use 'polyInterp' or 'nevilleDiffs' instead of evaluating+-- a fitted polynomial.+lagrangePolyFit :: Fractional a => [(a,a)] -> Poly a+lagrangePolyFit xys = sumPolys+    [ scalePoly f (fst (contractPoly p x))+    | f <- zipWith (/) ys phis+    | x <- xs+    ]+    where+        (xs,ys) = unzip xys+        p = lagrange xs+        phis = map (snd . evalPolyDeriv p) xs
+ src/Math/Polynomial/Lagrange.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE ParallelListComp #-}+module Math.Polynomial.Lagrange+    ( lagrangeBasis+    , lagrange+    , lagrangeWeights+    ) where++import Math.Polynomial++-- given a list, return one list containing each element of the original list+-- paired with all the other elements of the list.+select :: [a] -> [(a,[a])]+select [] = []+select (x:xs) = (x, xs) : [(y, x:ys) | (y, ys) <- select xs]++-- |Returns the Lagrange basis set of polynomials associated with a set of +-- points. This is the set of polynomials each of which is @1@ at its +-- corresponding point in the input list and @0@ at all others.+--+-- These polynomials are especially convenient, mathematically, for +-- interpolation.  The interpolating polynomial for a set of points  @(x,y)@ +-- is given by using the @y@s as coefficients for the basis given by +-- @lagrangeBasis xs@.  Computationally, this is not an especially stable +-- procedure though.  'Math.Polynomial.Interpolation.lagrangePolyFit'+-- implements a slightly better algorithm based on the same idea.  +-- +-- Generally it is better to not compute the coefficients at all.  +-- 'Math.Polynomial.Interpolation.polyInterp' evaluates the interpolating+-- polynomial directly, and is both quicker and more stable than any method+-- I know of that computes the coefficients.+lagrangeBasis :: Fractional a => [a] -> [Poly a]+lagrangeBasis xs =+    [ foldl1 multPoly+        [ if q /= 0+            then poly LE [negate x_j/q, 1/q]+            else error ("lagrangeBasis: duplicate root: " ++ show x_i)+        | x_j <- otherXs+        , let q = x_i - x_j+        ]+    | (x_i, otherXs) <- select xs+    ]++-- |Construct the Lagrange "master polynomial" for the Lagrange barycentric form:+-- That is, the monic polynomial with a root at each point in the input list.+lagrange :: Num a => [a] -> Poly a+lagrange [] = one+lagrange xs = foldl1 multPoly+    [ poly LE [negate x_i, 1]+    | x_i <- xs+    ]++-- |Compute the weights associated with each abscissa in the Lagrange+-- barycentric form.+lagrangeWeights :: Fractional a => [a] -> [a]+lagrangeWeights xs = +    [ recip $ product+        [ x_i - x_j+        | x_j <- otherXs+        ]+    | (x_i, otherXs) <- select xs+    ]
+ src/Math/Polynomial/Legendre.hs view
@@ -0,0 +1,69 @@+{-# LANGUAGE ParallelListComp #-}+module Math.Polynomial.Legendre where++import Math.Polynomial++-- |The Legendre polynomials with 'Rational' coefficients.  These polynomials +-- form an orthogonal basis of the space of all polynomials, relative to the +-- L2 inner product on [-1,1] (which is given by integrating the product of+-- 2 polynomials over that range).+legendres :: [Poly Rational]+legendres = one : x : +    [ multPoly+        (poly LE [recip (n' + 1)])+        (addPoly (poly LE [0, 2 * n' + 1] `multPoly` p_n)+                 (poly LE           [-n'] `multPoly` p_nm1)+        )+    | n     <- [1..], let n' = fromInteger n+    | p_n   <- tail legendres+    | p_nm1 <- legendres+    ]++-- |Compute the coefficients of the n'th Legendre polynomial.+legendre :: Fractional a => Int -> Poly a+legendre n = poly LE . map fromRational . polyCoeffs LE $ legendres !! n++-- |Evaluate the n'th Legendre polynomial at a point X.  Both more efficient+-- and more numerically stable than computing the coefficients and evaluating+-- the polynomial.+evalLegendre :: Fractional a => Int -> a -> a+evalLegendre n t = evalLegendres t !! n++-- |Evaluate all the Legendre polynomials at a point X.+evalLegendres :: Fractional a => a -> [a]+evalLegendres t = ps+    where+       ps = 1 : t : +            [ ((2 * n + 1) * t * p_n - n * p_nm1) / (n + 1)+            | n     <- iterate (1+) 1+            | p_n   <- tail ps+            | p_nm1 <- ps+            ]++-- |Evaluate the n'th Legendre polynomial and its derivative at a point X.  +-- Both more efficient and more numerically stable than computing the+-- coefficients and evaluating the polynomial.+evalLegendreDeriv :: Fractional a => Int -> a -> (a,a)+evalLegendreDeriv 0 _ = (1,0)+evalLegendreDeriv n t = case drop (n-1) (evalLegendres t) of+    (p2:p1:_)   -> (p1, fromIntegral n * (t * p1 - p2) / (t*t - 1))+    _ -> error "evalLegendreDeriv: evalLegendres didn't return a long enough list" {- should be infinite -}++-- |Zeroes of the n'th Legendre polynomial.+legendreRoots :: (Fractional b, Ord b) => Int -> b -> [b]+legendreRoots n eps = map negate mRoots ++ reverse (take (n-m) mRoots)+    where+        -- the roots are symmetric in the interval so we only have to find 'm' of them.+        -- The rest are reflections.+        m = (n + 1) `div` 2+        mRoots = [improveRoot (z0 i) | i <- [0..m-1]]+        +        -- Initial guess for i'th root of the n'th Legendre polynomial+        z0 i = realToFrac (cos (pi * (fromIntegral i + 0.75) / (fromIntegral n + 0.5)) :: Double)+        -- Improve estimate of a root by newton's method+        improveRoot z1+            | abs (z2-z1) <= eps    = z2+            | otherwise             = improveRoot z2+            where+                (y, dy) = evalLegendreDeriv n z1+                z2 = z1 - y/dy
+ src/Math/Polynomial/Newton.hs view
@@ -0,0 +1,16 @@+module Math.Polynomial.Newton where++import Math.Polynomial+import Data.List++-- |Returns the Newton basis set of polynomials associated with a set of +-- abscissas.  This is the set of monic polynomials each of which is @0@ +-- at all previous points in the input list.+newtonBasis :: Num a => [a] -> [Poly a]+newtonBasis xs = +    [ foldl multPoly (poly LE [1]) +        [ poly LE [-x_i, 1]+        | x_i <- xs'+        ]+    | xs' <- inits xs+    ]
+ src/Math/Polynomial/NumInstance.hs view
@@ -0,0 +1,19 @@+{-# OPTIONS_GHC -fno-warn-orphans #-}+-- |This module exports a 'Num' instance for the 'Poly' type.+-- This instance does not implement all operations, because 'abs' and 'signum'+-- are simply not definable, so I have placed it into a separate module so+-- that I can make people read this caveat ;).+--+-- Use at your own risk.+module Math.Polynomial.NumInstance where++import Math.Polynomial++instance Num a => Num (Poly a) where+    fromInteger i = poly LE [fromInteger i]+    (+) = addPoly+    negate = negatePoly+    (*) = multPoly++    abs     = error    "abs cannot be defined for the Poly type"+    signum  = error "signum cannot be defined for the Poly type"
+ src/Math/Polynomial/Pretty.hs view
@@ -0,0 +1,72 @@+{-# LANGUAGE +        ParallelListComp, ViewPatterns,+        FlexibleInstances, FlexibleContexts, IncoherentInstances+  #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+-- This code is a big ugly mess, but it more or less works.  Someday I might+-- get around to cleaning it up.++-- |This module exports a 'Pretty' instance for the 'Poly' type.+module Math.Polynomial.Pretty () where++import Math.Polynomial.Type++import Data.Complex++import Text.PrettyPrint+import Text.PrettyPrint.HughesPJClass++instance (Pretty a, Num a, Ord a) => Pretty (Poly a) where+    pPrintPrec l p x = ppr+        where+            ppr    = pPrintPolyWith p BE (pPrintOrdTerm pPrNum 'x') x+            pPrNum = pPrintPrec l 11++instance (RealFloat a, Pretty a) => Pretty (Complex a) where+    pPrintPrec l p (a :+ b) = ppr+        where+            x = poly LE [a,b]+            ppr = pPrintPolyWith p LE (pPrintOrdTerm pPrNum 'i') x+            pPrNum = pPrintPrec l 11++instance (RealFloat a, Pretty (Complex a)) => Pretty (Poly (Complex a)) where+    pPrintPrec l p x = ppr+        where+            ppr    = pPrintPolyWith p BE (pPrintUnOrdTerm pPrNum 'x') x+            pPrNum = pPrintPrec l 11++pPrintPolyWith prec end v p = parenSep (prec > 5) $ filter (not . isEmpty)+    [ v first coeff exp+    | (coeff, exp) <- +        (if end == BE then reverse else dropWhile ((0==).fst))+        (zip (polyCoeffs LE p) [0..])+    | first <- True : repeat False+    ]++parenSep p xs = +    prettyParen (p && not (null (drop 1 xs)))   +        (hsep xs)++pPrintOrdTerm   _ _ _ 0 _ = empty+pPrintOrdTerm num _ f c 0 = sign f c <> num (abs c)+pPrintOrdTerm   _ v f c 1   | abs c == 1    = sign f c <> char v+pPrintOrdTerm num v f c 1 = sign f c <> num (abs c) <> char v+pPrintOrdTerm   _ v f c e   | abs c == 1    = sign f c <> char v <> text "^" <> int e+pPrintOrdTerm num v f c e = sign f c <> num (abs c) <> char v <> text "^" <> int e++sign True x+    | x < 0     = char '-'+    | otherwise = empty+sign False x+    | x < 0     = text "- "+    | otherwise = text "+ "++pPrintUnOrdTerm   _ _ _ 0 _ = empty+pPrintUnOrdTerm num _ f c 0 = sign f 1 <> num c+pPrintUnOrdTerm   _ v f 1 1 = sign f 1 <> char v+pPrintUnOrdTerm num v f c 1 = sign f 1 <> num c <> char v+pPrintUnOrdTerm   _ v f 1 e = sign f 1 <> char v <> text "^" <> int e+pPrintUnOrdTerm num v f c e = sign f 1 <> num c <> char v <> text "^" <> int e+
+ src/Math/Polynomial/Type.hs view
@@ -0,0 +1,105 @@+{-# LANGUAGE ViewPatterns, TypeFamilies #-}+module Math.Polynomial.Type +    ( Endianness(..)+    , Poly, poly, polyCoeffs+    , polyIsZero, polyIsOne+    ) where++-- import Data.List.Extras.LazyLength+import Data.AdditiveGroup+import Data.VectorSpace+import Data.List.ZipSum++dropEnd :: (a -> Bool) -> [a] -> [a]+-- dropEnd p = reverse . dropWhile p . reverse+dropEnd p = go id+    where+        go t (x:xs)+            -- if p x, stash x (will only be used if 'not (any p xs)')+            | p x       =        go (t.(x:))  xs+            -- otherwise insert x and all stashed values in output and reset the stash+            | otherwise = t (x : go  id       xs)+        -- 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)++-- |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)++-- |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++polyIsZero :: Num a => Poly a -> Bool+polyIsZero = null . coeffs . trim++polyIsOne :: Num a => Poly a -> Bool+polyIsOne = ([1]==) . coeffs . trim++data Endianness +    = BE +    -- ^ Big-Endian (head is highest-order term)+    | LE+    -- ^ Little-Endian (head is const term)+    deriving (Eq, Ord, Enum, Bounded, Show)++data Poly a = Poly +    { endianness :: !Endianness+    , _trimmed   :: !Bool+    , coeffs     :: ![a]+    }++instance Num a => Show (Poly a) where+    showsPrec p (trim -> Poly end _ cs) +        = showParen (p > 10) +            ( showString "poly "+            . showsPrec 11 end+            . showChar ' '+            . showsPrec 11 cs+            )++instance (Num a, Eq a) => Eq (Poly a) where+    p == q  +        | endianness p == endianness q+        = coeffs (trim p) == coeffs (trim q)+        | otherwise +        = polyCoeffs BE p == polyCoeffs BE q+        ++-- -- Ord would be nice for some purposes, but it really just doesn't+-- -- make sense (there is no natural order that is much better than any+-- -- other, AFAIK), so I'm leaving it out.+-- instance (Num a, Ord a) => Ord (Poly a) where+--     compare p q = mconcat+--             [ lengthCompare pCoeffs qCoeffs+--             , compare       pCoeffs qCoeffs+--             ]+--         where+--             pCoeffs = polyCoeffs BE p+--             qCoeffs = polyCoeffs BE q++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)+    negateV = fmap negateV++instance VectorSpace a => VectorSpace (Poly a) where+    type Scalar (Poly a) = Scalar a+    (*^) s = fmap (s *^)