diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,5 @@ +#!/usr/bin/env runhaskell + +> import Distribution.Simple +> main = defaultMain + diff --git a/polynomial.cabal b/polynomial.cabal new file mode 100644 --- /dev/null +++ b/polynomial.cabal @@ -0,0 +1,38 @@ +name: polynomial +version: 0.5 +stability: provisional + +cabal-version: >= 1.6 +build-type: Simple + +author: James Cook +maintainer: James Cook +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 diff --git a/src/Data/List/ZipSum.hs b/src/Data/List/ZipSum.hs new file mode 100644 --- /dev/null +++ b/src/Data/List/ZipSum.hs @@ -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 + diff --git a/src/Math/Polynomial.hs b/src/Math/Polynomial.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial.hs @@ -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) diff --git a/src/Math/Polynomial/Bernstein.hs b/src/Math/Polynomial/Bernstein.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Bernstein.hs @@ -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 diff --git a/src/Math/Polynomial/Chebyshev.hs b/src/Math/Polynomial/Chebyshev.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Chebyshev.hs @@ -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] diff --git a/src/Math/Polynomial/Interpolation.hs b/src/Math/Polynomial/Interpolation.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Interpolation.hs @@ -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 diff --git a/src/Math/Polynomial/Lagrange.hs b/src/Math/Polynomial/Lagrange.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Lagrange.hs @@ -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 + ] diff --git a/src/Math/Polynomial/Legendre.hs b/src/Math/Polynomial/Legendre.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Legendre.hs @@ -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 diff --git a/src/Math/Polynomial/Newton.hs b/src/Math/Polynomial/Newton.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Newton.hs @@ -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 + ] diff --git a/src/Math/Polynomial/NumInstance.hs b/src/Math/Polynomial/NumInstance.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/NumInstance.hs @@ -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" diff --git a/src/Math/Polynomial/Pretty.hs b/src/Math/Polynomial/Pretty.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Pretty.hs @@ -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 + diff --git a/src/Math/Polynomial/Type.hs b/src/Math/Polynomial/Type.hs new file mode 100644 --- /dev/null +++ b/src/Math/Polynomial/Type.hs @@ -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 *^)