#What is a Polynomial?
Here, we will implement standard mathematical operations on polynomials.
If you do not know what a polynomial is, here is an example:
``` haskell
f x = x**2 + 3*x + 2
```
In the example, `f` is a _polynomial function_, and the right-hand-side is an expression of a polynomial.
##How will we represent a Polynomial?
More specifically, we will represent polynomials as a list of coefficients:
``` haskell
(x**2 + 3*x + 2) {- is represented by: -} [2,3,1]
(20) {- is represented by: -} [20]
(x) {- is represented by: -} [0,1]
(c0 + c1*x + c2*x**2 + ...) {- is represented by: -} [c0,c1,c2,..]
```
##Proper and Improper Representations
Mathematically, something like (`0*x**2 + x + 0`) is equal to (`x`), but one of them contains more symbols than the other.
For simplicity, we will say:
``` haskell
(0*x**2 + x + 0) {- is represented by: -} [0,1,0]
(x) {- is represented by: -} [0,1]
```
so that by looking at the expression of a polynomial, we can immediately get a list of coefficients.
This decision implies that a polynomial like `0` has an infinity of representations:
``` haskell
listOfAllRepresentationsOfZero = [[],[0],[0,0],[0,0,0],..]
```
To avoid unnecessary memory usage, we would like to use `[]` to represent `0`; it is the minimal representation. So, we say that the minimal representation is the proper representation, and any other representation is improper.
It is possible to define a function which takes any finite representation, and turns it into a proper representation:
``` active haskell
import Data.List
-- show This gets rid of all the zeroes at the end of the representation.
toProper [] = []
toProper p = if (last p /= 0) -- this means there is nothing to do.
then p
else toProper $ init p -- this means the last zero can be discarded.
main = do print $ toProper []
print $ toProper [0,0,0]
print $ toProper [0,1]
print $ toProper [0,1,0]
-- /show
```
#Operations
First, we implement operations like `(+)` for representations which are not necessarily proper. Afterwards, we will use `toProper` to "close" these operations for proper polynomials.
##You can add two Polynomials
``` active haskell
addPoly p1 p2 = if (length p1 >= length p2)
then zipWith (+) p1 (p2 ++ repeat)
else addPoly p2 p1
{-
if you have
(2*x**2 + 3*x + 2) + (x)
you rewrite
(2*x**2 + 3*x + 2) + (0*x**2 + x + 0) = (2+0)*x**2 + (3+1)*x + (2+0)
-}
main = print $ [0,1] `addPoly` [2,3,2]
```
##You can multiply a polynomial by a constant
``` active haskell
multiplyBy a p1 = map (a*) p1
{-
if you have
4 * (2*x**2 + 3*x + 2) = (4*2)*x**2 + (4*3)*x + (4*2)
-}
main = do print $ 4 `multiplyBy` [2,3,2]
```
We see that this operation is already closed over proper representations.
##You can multiply by `x`
``` active haskell
multiplyByX p = 0:p
{-
if you have
(x) * (x**2 + 3*x + 4) = (x**3 + 3*x**2 + 4*x + 0)
-}
main = do print $ multiplyByX [4,3,1]
```
We see that this is not closed over proper representations. However, it is almost closed.
Exercise: Find the only case where this transforms a proper representation into an improper one.
##We can multiply two polynomials
``` active haskell
-- show
multPoly [] p2 = []
multPoly (p:p1) p2 = let pTimesP2 = multiplyBy p p2
xTimesP1Timesp2 = multiplyByX $ multPoly p1 p2
in addPoly pTimesP2 xTimesP1Timesp2
{-
if we have
(2*x + 1) * (x**2 + 3*x + 2) = 1 * (x**2 + 3*x + 2) + (x) * (2) * (x**2 + 3*x + 2)
-}
main = print $ multPoly [1,2] [2,3,1]
-- /show
addPoly p1 p2 = if (length p1 >= length p2)
then zipWith (+) p1 (p2 ++ repeat 0)
else addPoly p2 p1
multiplyBy a p1 = map (a*) p1
multiplyByX p = 0:p
```
##We can negate a polynomial (coefficient-wise)
``` active haskell
negatePoly = map negate
main = print $ negatePoly [0,1]
```
##We can derive a polynomial
``` active haskell
derive [] = []
derive (_:ps) = zipWith (*) ps [1..]
main = do print $ derive [0,1]
print $ derive [0,1,4,2]
```
#Instances
Here are three important instances: `Num`, `Show`, `Eq`.
``` active haskell
import Data.List
multPoly [] p2 = []
multPoly (p:p1) p2 = let pTimesP2 = multiplyBy p p2
xTimesP1Timesp2 = multiplyByX $ multPoly p1 p2
in addPoly pTimesP2 xTimesP1Timesp2
addPoly p1 p2 = if (length p1 >= length p2)
then zipWith (+) p1 (p2 ++ repeat 0)
else addPoly p2 p1
multiplyBy a p1 = map (a*) p1
multiplyByX p = 0:p
negatePoly :: Num a => [a] -> [a]
negatePoly = map negate
toProper [] = []
toProper p = if (last p /= 0) -- this means there is nothing to do.
then p
else toProper $ init p -- this means the last zero can be discarded.
-- show This allows us to instantiate typeclasses to which regular lists belong.
newtype Poly a = P [a]
-- /show
-- show A polynomial is sort of a number. We enforce closure by properP, which only creates proper Ps.
--By using modules, we could only export properP, not P. This would simply get rid
--of all those improper representations. However, some things are more efficient
--with improper representations, so we do not want to lose that power. Depending
--on the application, this constraint might be removed.
properP :: (Num a, Eq a) => [a] -> Poly a
properP = P . toProper
instance (Num a, Eq a) => Num (Poly a) where
(P a) + (P b) = properP $ addPoly a b
(P a) * (P b) = properP $ multPoly a b
negate (P a) = properP $ negatePoly a
abs = undefined
signum = undefined
fromInteger i = properP [fromIntegral i]
-- /show
-- show Shows polynomials how we write them mathematically.
-- It would be harder to parse.
showPoly [] = show 0
showPoly p = let cOs = zip p [0..]
nonZeroCOs = filter (\(c,_) -> c /= 0) cOs
cShow c = if c == 1
then ""
else show c
nShow n = case n of
0 -> ""
1 -> "x"
m -> "x^" ++ show m
cnShow c n = if c == 1 && n == 0
then show 1
else intercalate " " $ filter (/="") [cShow c, nShow n]
terms = map (\(c,n) -> cnShow c n) nonZeroCOs
in intercalate " + " (reverse terms)
instance (Show a, Eq a, Num a) => Show (Poly a) where
show (P a) = showPoly $ toProper a
-- /show
-- show Eqs two polynomials
instance (Num a, Eq a) => Eq (Poly a) where
(P a) == (P b) = toProper a == toProper b
-- /show
-- show Try it!
main = do print (P [1,3,4] - P [0,1,0,0])
print $ P []
print $ P [0,0]
print (P [3,2,5] * P [4,1,1])
print (P [3,2,5] * P [4,1,1] == properP [4,1,1] * P [3,2,5,0,0])
-- /show
```
Exercise: 1) modify the instance of Show so that multiplication is denoted by `*` (as in `2*x + 1`). 2) (Bonus) modify the instance of Show so that things like `2 x^2 + -3` become `2 x^2 - 3`.
#Polynomials as functions
``` haskell
-- show
listOfPowers x = map (\n -> x**n) [0..]
makeFunction p = \x -> sum $ zipWith (*) p (listOfPowers x)
-- /show
```
#Important Families of Polynomials
``` active haskell
import Data.List
multPoly [] p2 = []
multPoly (p:p1) p2 = let pTimesP2 = multiplyBy p p2
xTimesP1Timesp2 = multiplyByX $ multPoly p1 p2
in addPoly pTimesP2 xTimesP1Timesp2
addPoly p1 p2 = if (length p1 >= length p2)
then zipWith (+) p1 (p2 ++ repeat 0)
else addPoly p2 p1
multiplyBy a p1 = map (a*) p1
multiplyByX p = 0:p
negatePoly :: Num a => [a] -> [a]
negatePoly = map negate
toProper [] = []
toProper p = if (last p /= 0) -- this means there is nothing to do.
then p
else toProper $ init p -- this means the last zero can be discarded.
newtype Poly a = P [a]
properP :: (Num a, Eq a) => [a] -> Poly a
properP = P . toProper
instance (Num a, Eq a) => Num (Poly a) where
(P a) + (P b) = properP $ addPoly a b
(P a) * (P b) = properP $ multPoly a b
negate (P a) = properP $ negatePoly a
abs = undefined
signum = undefined
fromInteger i = properP [fromIntegral i]
showPoly [] = show 0
showPoly p = let cOs = zip p [0..]
nonZeroCOs = filter (\(c,_) -> c /= 0) cOs
cShow c = if c == 1
then ""
else show c
nShow n = case n of
0 -> ""
1 -> "x"
m -> "x^" ++ show m
cnShow c n = if c == 1 && n == 0
then show 1
else intercalate " " $ filter (/="") [cShow c, nShow n]
terms = map (\(c,n) -> cnShow c n) nonZeroCOs
in intercalate " + " (reverse terms)
instance (Show a, Eq a, Num a) => Show (Poly a) where
show (P a) = showPoly $ toProper a
instance (Num a, Eq a) => Eq (Poly a) where
(P a) == (P b) = toProper a == toProper b
-- show Legendre Polynomials
x = P [0,1]
coef a = properP [a]
coeff `mono` order = properP (replicate order 0 ++ [coeff])
--this is the zeroth legendre polynomial
leg0 = coef 1
--this is the 1st legendre poly.
leg1 = x
--this is the relationship between the current legendre
--polynomial and the two previous ones.
legN n legP legPP = coef (1/n) * (coef (2*n-1) * x * legP + coef (1-n) * legPP)
--the infinite list of legendre polynomials
legs = leg0 : leg1 : zipWith3 legN [2..] (tail legs) legs
getLegendre n = legs !! n
derive (P []) = (P [])
derive (P (_:ps)) = P $ zipWith (*) ps [1..]
getLegendre' = derive . getLegendre
-- /show
--show Chebyshev Polynomials
ch0 = coef 1
ch1 = x
chN chPrev chPrevPrev = 2 * x * chPrev - chPrevPrev
chebys = ch0:ch1: zipWith chN (tail chebys) chebys
getCheby n = chebys !! n
-- /show
-- show Try it!
main = do putStrLn "Legendre polynomials:"
putStrLn $ intercalate "\n" $ map show $ take 10 legs
putStrLn "Associated Legendre polynomials:"
putStrLn $ intercalate "\n" $ map show $ map getLegendre' [0..9]
putStrLn "Chebyshev Polynomials:"
putStrLn $ intercalate "\n" $ map show $ take 10 chebys
-- /show
```
Exercise: Define the set of all legendre polynomials using 1) a recursion, 2) an iteration. 3) (Bonus) Why is it a bad idea to do this recursively for n = 1000? 4) Is it a good idea to use an iteration for n = 1000 ?
#Conclusion
Now that we have the power of polynomials, there are many important applications to mention. See my tutorial on Numerical Methods for the actual applications.
We can do interpolation with polynomials. The idea is to construct a polynomial which interpolates a set of points.
We can find the roots of simple polynomials exactly:
``` haskell
realquadroots (P [c,b,a]) = let disc = b**2 - 4 * a * c
in if disc < 0
then []
else if disc == 0
then [-b/(2*a)]
else [(-b + sqrt disc)/(2*a),(-b - sqrt disc)/(2*a)]
```
We can also use `makeFunction` to find the roots of simple polynomials approximately (with Newton's method, for instance).
You can also use the legendre polynomials in more complex techniques such as the Gaussian Quadrature.
#Bonus Round: Getting rid of Improper Representations.
Here is how you ensure the utter annilation of Improper Polynomials. This technique can be used for your own applications with a similar structure.
The idea is that we will make `P` dissapear, leaving only `properP`. You see, if all of our operations are closed, and there is no way to introduce Improper representations, then _it becomes impossible to get an Improper polynomial_.
First, we have to introduce explicitely a function we used without knowing it.
``` haskell
toList (P p) = p
```
Because we want to still be able to work directly on the list (this is necessary if we want to define additional functions).
So now, we put what we have done so far in a module, and only export the functions in the parantesis.
The user of our module will work in another file (Main.hs).
``` active haskell
{-# START_FILE Poly.hs #-}
-- show
module Poly (toList, makeFunction, properP, Poly(), x, coef, mono, legs, getLegendre, getLegendre', derive, chebys, getCheby) where
--Note: the instances for Show, Num, and Eq are exported.
--There is actually more stuff here, but it is all a repeat of what we saw previously.
-- /show
import Data.List
toList (P p) = p
listOfPowers x = map (\n -> x**n) [0..]
makeFunction (P p) = \x -> sum $ zipWith (*) p (listOfPowers x)
multPoly [] p2 = []
multPoly (p:p1) p2 = let pTimesP2 = multiplyBy p p2
xTimesP1Timesp2 = multiplyByX $ multPoly p1 p2
in addPoly pTimesP2 xTimesP1Timesp2
addPoly p1 p2 = if (length p1 >= length p2)
then zipWith (+) p1 (p2 ++ repeat 0)
else addPoly p2 p1
multiplyBy a p1 = map (a*) p1
multiplyByX p = 0:p
negatePoly :: Num a => [a] -> [a]
negatePoly = map negate
toProper [] = []
toProper p = if (last p /= 0) -- this means there is nothing to do.
then p
else toProper $ init p -- this means the last zero can be discarded.
newtype Poly a = P [a]
properP :: (Num a, Eq a) => [a] -> Poly a
properP = P . toProper
instance (Num a, Eq a) => Num (Poly a) where
(P a) + (P b) = properP $ addPoly a b
(P a) * (P b) = properP $ multPoly a b
negate (P a) = properP $ negatePoly a
abs = undefined
signum = undefined
fromInteger i = properP [fromIntegral i]
showPoly [] = show 0
showPoly p = let cOs = zip p [0..]
nonZeroCOs = filter (\(c,_) -> c /= 0) cOs
cShow c = if c == 1
then ""
else show c
nShow n = case n of
0 -> ""
1 -> "x"
m -> "x^" ++ show m
cnShow c n = if c == 1 && n == 0
then show 1
else intercalate " " $ filter (/="") [cShow c, nShow n]
terms = map (\(c,n) -> cnShow c n) nonZeroCOs
in intercalate " + " (reverse terms)
instance (Show a, Eq a, Num a) => Show (Poly a) where
show (P a) = showPoly $ toProper a
instance (Num a, Eq a) => Eq (Poly a) where
(P a) == (P b) = toProper a == toProper b
x = P [0,1]
coef a = properP [a]
coeff `mono` order = properP (replicate order 0 ++ [coeff])
--this is the zeroth legendre polynomial
leg0 = coef 1
--this is the 1st legendre poly.
leg1 = x
--this is the relationship between the current legendre
--polynomial and the two previous ones.
legN n legP legPP = coef (1/n) * (coef (2*n-1) * x * legP + coef (1-n) * legPP)
--the infinite list of legendre polynomials
legs = leg0 : leg1 : zipWith3 legN [2..] (tail legs) legs
getLegendre n = legs !! n
derive (P []) = (P [])
derive (P (_:ps)) = P $ zipWith (*) ps [1..]
getLegendre' = derive . getLegendre
ch0 = coef 1
ch1 = x
chN chPrev chPrevPrev = 2 * x * chPrev - chPrevPrev
chebys = ch0:ch1: zipWith chN (tail chebys) chebys
getCheby n = chebys !! n
{-# START_FILE Main.hs #-}
import Poly
import Data.List
main = do putStrLn "Legendre polynomials:"
putStrLn $ intercalate "\n" $ map show $ take 10 legs
putStrLn "Testing the arithmetic operations"
print (properP [1,3,4] - properP [0,1,0,0])
print $ properP []
print $ properP [0,0]
print (properP [3,2,5] * properP [4,1,1])
putStrLn "Testing the Equality test"
print (properP [3,2,5] * properP [4,1,1] == properP [4,1,1] * properP [3,2,5,0,0])
putStrLn "Testing derivation"
print (derive (properP [1,0,2]))
```
Exercise: 1) write a test for `makeFunction`, 2) try to replace `properP` by `P` in Main.hs; what happens?
3) (Bonus \*) Try to implement Polynomials as a Vector (from Data.Vector.Unboxed) of coefficient instead of a list (do all this in Poly.hs), but you must keep Main.hs the way it is, and have all the tests give you the same thing as they give you now.
Have fun.
# Bonus Round 2: Integration
In the same way that polynomials can be derived exactly according to calculus, so too can they be integrated.
The integral of a polynomial is a polynomial, which is rather nice. (as I will discuss later, this kind of property is often very good if you want your code to be powerful and relatively small)
Of course, integration (indefinite integration) does not fully specify what the result should be. We know that the following property must hold:
``` haskell
derive (integrate p) == p
```
This only specifies integration up to a constant. In other words,
``` haskell
integrate (derive p) == p + c
```
where `c` is a constant polynomial.
For simplicity, we use the convention that integration always gives a polynomial going through the origin. For instance,
``` haskell
integrate (P [1]) == P [0,1]
```
## Why not just do definite integration?
The ultimate goal is definite integration, but suppose we want to visualise the integral of some polynomial. In that case, it is more efficient (and more interesting) to first calculate the indefinite integral, and then evaluate it between various points plot the curve.
Exercise: Try to implement polynomial integration.
Solution:
@@@
``` haskell
integrate (P p) = P (0 : zipWith (/) p (map fromIntegral [1..length p]))
```
@@@
# Solution to using vectors and optimising a bit.
As you can see below, programs using vectors are a little less descriptive of the meaning, and more descriptive of the procedures. However, after understanding the list version, you can use the vector version without difficulty. Notice also how I didn't convert the `show` function to a vector style. The `show` function doesn't need to be super-fast, and the tools to work with `String`s are using lists (because `String` is just `[Char]`). This is the kind of compromise in favor of convenience you should always be making.
A good rule of thumb is this: functions that take some `a` and makes a new `a` are potential bottlenecks, because some bigger function could call this smaller function a billion times in some program. However, a function that takes some `a` and returns something that only *you* will understand ... well, this will never be the bottleneck; you and your understanding will (you can't look at a billion strings representing a polynomial faster than your computer can print them).
@@@
``` haskell
module PolyV (toVec, makeFunction, properP, Poly(), x, coef, mono , legs, getLegendre, getLegendre', derive, integrate, chebys, getCheby) where
import Data.List (intercalate)
import Numeric.Polynomial
import qualified Data.Vector.Generic as G
import Data.Vector.Generic (Vector)
import qualified Data.Vector.Unboxed as U
toVec :: Poly t -> U.Vector t
toVec (P p) = p
{-# INLINE toVec #-}
makeFunction :: (Num a, U.Unbox a) => Poly a -> (a -> a)
makeFunction (P p) = \x -> evaluatePolynomial x p
{-# INLINE makeFunction #-}
multPoly :: (Eq a, Num a, U.Unbox a) => U.Vector a -> U.Vector a -> U.Vector a
multPoly xs myYs = if U.null myYs
then U.empty
else let
y = U.head myYs
ys = U.tail myYs
mul 0 bs = 0 `U.cons` bs
mul x bs = (x*y) `U.cons` addPoly (U.map (x*) ys) bs
in U.foldr mul U.empty xs
{-# INLINE multPoly #-}
addPoly :: (Num a, U.Unbox a) => U.Vector a -> U.Vector a -> U.Vector a
addPoly p1 p2 = let
n1 = U.length p1
n2 = U.length p2
in if n1 >= n2
then let
(fittingPart, rest) = U.splitAt n2 p1
in (U.zipWith (+) fittingPart p2) U.++ rest
else addPoly p2 p1
{-# INLINE addPoly #-}
multiplyBy :: (Num b, U.Unbox b) => b -> U.Vector b -> U.Vector b
multiplyBy a p1 = U.map (a*) p1
{-# INLINE multiplyBy #-}
multiplyByX :: (Num a, U.Unbox a) => U.Vector a -> U.Vector a
multiplyByX p = 0 `U.cons` p
{-# INLINE multiplyByX #-}
negatePoly ::(U.Unbox b0, Num b0) => U.Vector b0 -> U.Vector b0
negatePoly = U.map negate
{-# INLINE negatePoly #-}
toProper :: (Eq a, Num a, U.Unbox a) => U.Vector a -> U.Vector a
toProper p = if (U.null p || U.last p /= 0)
then p -- this means there is nothing to do.
else toProper $ U.init p -- this means the last zero can be discarded.
{-# INLINE toProper #-}
newtype Poly a = P (U.Vector a)
properP :: (Eq a, Num a, U.Unbox a) => U.Vector a -> Poly a
properP = P . toProper
{-# INLINE properP #-}
showPoly mp = if U.null mp
then show 0
else let
p = U.toList mp
cOs = zip p [0..]
nonZeroCOs = filter (\(c,_) -> c /= 0) cOs
cShow c = if c == 1
then ""
else show c
nShow n = case n of
0 -> ""
1 -> "x"
m -> "x^" ++ show m
cnShow c n = if c == 1 && n == 0
then show 1
else intercalate " " $ filter (/="") [cShow c, nShow n]
terms = map (\(c,n) -> cnShow c n) nonZeroCOs
in intercalate " + " (reverse terms)
instance (Num a, Eq a, U.Unbox a) => Num (Poly a) where
(P a) + (P b) = properP $ addPoly a b
(P a) * (P b) = properP $ multPoly a b
negate (P a) = properP $ negatePoly a
abs = undefined
signum = undefined
fromInteger i = properP $ U.singleton $ fromIntegral i
{-
instance (Show a, Eq a, Num a) => Show (Poly a) where
show (P a) = showPoly $ toProper a
-}
instance (Num a, Eq a, U.Unbox a) => Eq (Poly a) where
(P a) == (P b) = toProper a == toProper b
x :: (Num a0, U.Unbox a0) => Poly a0
x = P $ 0 `U.cons` U.singleton 1
coef :: (Eq a, Num a, U.Unbox a) => a -> Poly a
coef a = properP $ U.singleton a
mono :: (Eq a, Num a, U.Unbox a) => a -> Int -> Poly a
coeff `mono` order = properP (U.replicate order 0 `U.snoc` coeff)
--this is the zeroth legendre polynomial
leg0 :: (Eq a, Num a, U.Unbox a) => Poly a
leg0 = coef 1
--this is the 1st legendre poly.
leg1 :: (Num a, U.Unbox a) => Poly a
leg1 = x
--this is the relationship between the current legendre
--polynomial and the two previous ones.
legN n legP legPP = coef (1/n) * (coef (2*n-1) * x * legP + coef (1-n) * legPP)
--the infinite list of legendre polynomials
legs :: (Eq a, Num a, U.Unbox a, Fractional a) => [Poly a]
legs = leg0 : leg1 : zipWith3 legN (map fromIntegral [2..]) (tail legs) legs
getLegendre :: (Eq a, Num a, U.Unbox a, Fractional a) => Int -> Poly a
getLegendre n = legs !! n
derivePoly p = if U.null p
then p
else U.zipWith (*) (U.tail p) (U.enumFromN 1 (U.length p - 1))
derive :: (Fractional a, U.Unbox a) => Poly a -> Poly a
derive = P . derivePoly . toVec
integratePoly p = 0 `U.cons` U.zipWith (/) p (U.enumFromN 1 (U.length p))
integrate :: (Fractional a, U.Unbox a) => Poly a -> Poly a
integrate = P . integratePoly . toVec
getLegendre' :: (Eq a, Num a, U.Unbox a, Fractional a) => Int -> Poly a
getLegendre' = derive . getLegendre
ch0 :: (Eq a, Num a, U.Unbox a) => Poly a
ch0 = coef 1
ch1 :: (Num a, U.Unbox a) => Poly a
ch1 = x
chN chPrev chPrevPrev = 2 * x * chPrev - chPrevPrev
chebys :: (Eq a, Num a, U.Unbox a) => [Poly a]
chebys = ch0:ch1: zipWith chN (tail chebys) chebys
getCheby n = chebys !! n
```
@@@