content: > {-# LANGUAGE TypeSynonymInstances #-}
  > module Lecture10 where
  
  Lecture 10: Polynomials and Power Series
  ========================================
  
  0. Preliminaries
  ----------------
  
  Last time, we defined a Num structure on pairs (Double, Double) by
  requiring the operations to be compatible with the interpretation (f
  a, f' a).  For example
  
  < (x, x') * (y, y') = (x * y, x' * y + x * y')
  
  There is nothing in the "nature" of pairs of Double that forces this
  definition upon us.  We chose it, because of the intended
  interpretation.
  
  This multiplication is obviously not the one we need for *complex
  numbers*:
  
  < (x, x') * (y, y') = (x * y - x' * y', x * y' + x' * y)
  
  Again, there is nothing in the nature of pairs that foists this
  operation on us.  In particular, it is, strictly speaking, incorrect
  to say that a complex number *is* a pair of real numbers.  The correct
  interpretation is that a complex number can be *represented* by a pair
  of real numbers, provided we define the operations on these pairs in a
  suitable way.
  
  The distinction between definition and representation is similar to
  the one between specification and implementation, and, in a certain
  sense, to the one between syntax and semantics.  All these
  distinctions are frequently obscured, for example, because of
  prototyping (working with representations/implementations/concrete
  objects in order to find out what definition/specification/syntax is
  most adequate).  They can also be context-dependent (one man's
  specification is another man's implementation).  Insisting on the
  difference between definition and representation can also appear quite
  pedantic (as in the discussion of complex numbers above).  In general
  though, it is a good idea to be aware of these distinctions, even if
  they are suppressed for reasons of brevity or style.
  
  1. Polynomials
  --------------
  
  From @adams2010calculus, page 55:
  
    > A **polynomial** is a function $P$ whose value at $x$ is
  
     $P x = a_n x^n + a_{n-1} x^{n - 1} + ... a_1 x + a_0$
  
    > where $a_n$, $a_{n-1}$, ..., $a_1$, and $a_0$, called the
    > **coefficients** of the polymonial [original spelling], are
    > constants and, if $n > 0$, then $a_n ≠ 0$.  The number $n$, the
    > degree of the highest power of $x$ in the polynomial, is called
    > the **degree** of the polynomial. (The degree of the zero
    > polynomial is not defined.)
  
  This definition raises a number of questions, for example "what is the
  zero polynomial?".
  
  The types of the elements involved in the definition appear to be
  
    > $P : ℝ → ℝ$, $x ∈ ℝ$, $a_0$, ... $a_n ∈ ℝ$ with $a_n ≠ 0$ if $n > 0$
  
  The phrasing should be "whose value at *any* $x$ is".  The remark
  that the $a_i$ are constants is probably meant to indicate that they
  do not depend on $x$, otherwise every function would be a polynomial.
  The zero polynomial is, according to this definition, the `const 0`
  function.  Thus, what is meant is
  
    > A **polynomial** is a function $P : ℝ → ℝ$ which is either constant zero,
    > or there exist $a_0$, ..., $a_n$ ∈ ℝ with $a_n ≠ 0$ such that, for
    > any $x ∈ ℝ$
  
    > $P x = a_n x^n + a_{n-1} x^{n - 1} + ... a_1 x + a_0$
  
  
  Obviously, given the coefficients $a_i$ we can evaluate $P$ at any
  given $x$.  Assuming the coefficients are given as
  
  <      as = [a0, a1, ..., an]
  
  (we prefer counting up), then the evaluation function is written
  
  <      eval   :: [Real] -> Real -> Real
  <      eval []  x       =  0
  <      eval (a : as) x  =  a + x * eval as x
  
  Not every list is valid according to the definition.  In particular,
  the empty list is not a valid list of coefficients, so we have a
  conceptual, if not empirical, type error in our evaluator.
  
  The valid lists are those *finite* lists in the set
  
  <      {[0]} ∪ {(a : as) | last (a : as) ≠ 0}
  
  We cannot express the "last (a : as) ≠ 0" in Haskell, but we can
  express the condition that the list should not be empty:
  
  > data Poly a  =  Single a  |  Cons a (Poly a)
  >                 deriving (Eq, Ord)
  
  The relationship between `Poly a` and `[a]` is given by the following
  functions:
  
  > toList :: Poly a -> [a]
  > toList (Single a)   =  [a]
  > toList (Cons a as)  =  a : toList as
  
  > fromList :: [a] -> Poly a
  > fromList [a] = Single a
  > fromList (a0 : a1 : as) = Cons a0 (fromList (a1 : as))
  
  > instance Show a => Show (Poly a) where
  >   show = show . toList
  
  Since we only use the arithmetical operations, we can generalise our
  evaluator:
  
  > evalPoly ::  Num a => Poly a -> a -> a
  > evalPoly (Single a)     x   =  a
  > evalPoly (Cons a as)    x   =  a + x * evalPoly as x
  
  Since we have Num a, there is a Num structure on `a -> a`, and
  `evalPoly` looks like a homomorphism.  Question: is there a Num
  structure on `Poly a`, such that `evalPoly` is a homomorphism?
  
  For example, the homomorphism condition gives for `(+)`
  
  <  evalPoly as + evalPoly bs = evalPoly (as + bs)
  
  Both sides are functions, they are equal iff they are equal for every
  argument.  For an arbitrary `x`
  
  <  (evalPoly as + evalPoly bs) x = evalPoly (as + bs) x
  
  <<=> {+ on functions is point-wise}
  
  <  evalPoly as x + evalPoly bs x = evalPoly (as + bs) x
  
  To proceed further, we need to consider the various cases in the
  definition of `evalPoly`.  We give here the computation for the last
  case, using the traditional list notation (:) for brevity.
  
  <  evalPoly (a : as) x  +  evalPoly (b : bs) x =
  <  evalPoly ((a : as)  +  (b : bs))
  
  For the left-hand side, we have:
  
  <  evalPoly (a : as) x  +  evalPoly (b : bs) x
  
  <=  {def evalPoly}
  
  <  a + x * evalPoly as x + b + x * eval bs x
  
  <=  {properties of +, valid in any ring}
  
  <  a + b + x * (evalPoly as x + evalPoly bs x)
  
  <=  {homomorphism condition}
  
  <  a + b + x * evalPoly (as + bs)
  
  <=  {def evalPoly}
  
  <  evalPoly ((a + b) : as + bs) x
  
  The homomorphism condition will hold for every `x` if we define
  
  <  (a : as) + (b : bs) = (a + b) : (as + bs)
  
  We leave the derivation of the other cases and operations as an
  exercise.  Here, we just give the corresponding definitions.
  
  > instance Num a => Num (Poly a) where
  >   Single a   +  Single b   =  Single (a + b)
  >   Single a   +  Cons b bs  =  Cons (a + b) bs
  >   Cons a as  +  Single b   =  Cons (a + b) as
  >   Cons a as  +  Cons b bs  =  Cons (a + b) (as + bs)
  >
  >   Single a   *  Single b   =  Single (a * b)
  >   Single a   *  Cons b bs  =  Cons (a * b) (Single a * bs)
  >   Cons a as  *  Single b   =  Cons (a * b) (as * Single b)
  >   Cons a as  *  Cons b bs  =  Cons (a * b) (as * Cons b bs + Single a * bs)
  >
  >   negate (Single a)        =  Single (negate a)
  >   negate (Cons a as)       =  Cons (negate a) (negate as)
  >
  >   fromInteger              =  Single . fromInteger
  
  
  Therefore, we *can* define a ring structure (the mathematical
  counterpart of Num) on `Poly a`, and we have arrived at the canonical
  definition of polynomials, as found in any algebra book (see, for
  example, @rotman2006first for a very readable text):
  
   > Given a commutative ring A, the commutative ring given by the set
   > `Poly A` together with the operations defined above is the ring of
   > **polynomials** with coefficients in A.
  
  The functions `evalPoly as` are known as *polynomial functions*.
  
  **Caveat:** The canonical representation of polynomials in algebra
    does not use finite lists, but the equivalent
  
    > `Poly' A` = { a : ℕ → A | a has only a finite number of non-zero values }
  
  Exercise: what are the ring operations on `Poly' A`?  For example,
  here is addition:
  
  <          a + b = c  <=>  a n + b n = c n    ∀ n : ℕ
  
  **Observations:**
  
    a) Polynomials are not, in general, isomorphic (in one-to-one
    correspondence) with polynomial functions.  For any finite ring A,
    there is a finite number of functions |A → A|, but there is a
    countable number of polynomials.  That means that the same
    polynomial function on A will be the evaluation of many different
    polynomials.
  
    > For example, consider the ring ℤ₂ ({0, 1} with addition and
      multiplication modulo 2).  In this ring, we have
  
  <    evalPoly [0, -1, 1] = const 0 = evalPoly [0]
  
    > but
  
  <    [0, -1, 1] ≠ [0]  in Poly ℤ₂
  
    > Therefore, it is not generally a good idea to confuse polynomials
      with polynomial functions.
  
    b) In keeping with the DSL terminology, we can say that the
    polynomial functions are the semantics of the language of
    polynomials.  We started with polynomial functions, we wrote the
    evaluation function and realised that we have the makings of a
    homomorphism.  That suggested that we could create an adequate
    language for polynomial functions.  Indeed, this turns out to be the
    case; in so doing, we have recreated an important mathematical
    achievement: the algebraic definition of polynomials.
  
  Let
  
  > x = Cons 0 (Single 1)
  
  Then (again, using the list notation for brevity) for any polynomial
  `as = [a0, a1, ..., an]` we have
  
  <  as = a0 + a1 * x + a2 * x^2 + ... + an * x^n
  
  Exercise: check this.
  
  This justifies the standard notation
  
  <  as = Σ_{i = 0}^n a_i * x^i
  
  2. Power series
  ----------------------
  
  Power series are obtained from polynomials by removing in `Poly'` the restriction
  that there should be a finite number of non-zero coefficients; or, in,
  the case of `Poly`, by going from lists to streams.
  
    PowerSeries' a = { f : ℕ → a }
  
  > type PowerSeries a = Poly a -- finite and infinite non-empty lists
  
  The operations are still defined as before.  If we consider only
  infinite lists, then only the equations which do not contain the
  patterns for singleton lists will apply.
  
  Power series are usually denoted
  
      Σ_{n = 0}^{\infty} a_n * x^n
  
  the interpretation of `x` being the same as before.
  
  The evaluation of a power series represented by a : ℕ → A is defined,
  in case the necessary operations make sense on A, as a function
  
  <   eval a : A -> A
  <   eval a x  =  lim s where s n = Σ_{i = 0}^n (a i) * x^i
  
  eval a is, in general a partial function (the limit might not exist).
  
  We will consider, as is usual, only the case in which A = ℝ or A = ℂ.
  
  The term *formal* refers to the independence of the definition of
  power series from the ideas of convergence and evaluation.  In
  particular, two power series represented by a and b, respectively, are
  equal only if a = b (as functions).  If a ≠ b, then the power series
  are different, even if eval a = eval b.
  
  Since we cannot compute limits, we can use an "approximative" `eval`,
  by evaluating the polynomial resulting from an initial segment of the
  power series.
  
  > eval n as x = evalPoly (takePoly n as) x
  
  > takePoly :: Integer -> Poly a -> Poly a
  > takePoly n (Single a)   =  Single a
  > takePoly n (Cons a as)  =  if n <= 1
  >                               then  Single a
  >                               else  Cons a (takePoly (n-1) as)
  
  3. Operations on power series
  -----------------------------
  
  Power series have a richer structure than polynomials.  For example,
  we also have division (this is similar to the move from ℤ to ℚ).
  Assume that a * b ≠ 0.  Then (again, using list notation for brevity),
  we want to find, for any given (a : as) and (b : bs), the series (c :
  cs) satisfying
  
  <  (a : as) / (b : bs) = c : cs
  
  <<=>
  
  <  (a : as) = (c : cs) * (b : bs)
  
  <<=>
  
  <  (a : as) = c * b : (cs * (b : bs) +
  <                              [c]*bs)
  
  <<=>
  
  <  c = a / b and
  <  as = cs * (b : bs) + [a/b] * bs
  
  <<=>
  
  <  c = a / b and
  <  cs =  (as - [a/b] * bs) / (b : bs)
  
  This leads to the implementation:
  
  > instance (Eq a, Fractional a) => Fractional (PowerSeries a) where
  >   as / Single b           =  as * Single (1 / b)
  >   Single a / Cons b bs    =  if a == 0 then Single 0 else Cons a (Single 0) / Cons b bs
  >   Cons a as / Cons b bs   =  let  q = a / b
  >                              in   Cons q  ((as - Single q * bs) / Cons b bs)
  >   fromRational            =  Single . fromRational
  
  The first two equations allow us to also use division on polynomials,
  but the result will, in general, be a power series, not a polynomial.
  The first one should be self-explanatory.  The second one extends a
  constant polynomial, in a process similar to that of long division.
  
  For example:
  
  > test0, test1, test2 :: PowerSeries Double
  > test0 = 1 / (1 - x)
  > test1 = 1 / (1 - x)^2
  > test2 = (x^2 - 2 * x + 1) / (x - 1)
  
  Every test is the result of a division of polynomials: the first two
  return power series, the third is a polynomial (almost: it has a
  trailing 0.0).
  
  > x :: Num a => Poly a -- defined earlier
  
  4. Formal derivative
  --------------------
  
  Considering the analogy between power series and polynomial functions
  (via polynomials), we can define a formal derivative for power series
  according to the formula
  
  <    (Σ a_n * X^n)'  =  Σ (n * a_n * X^(n-1))
  
  We can implement this, for example, as
  
  > deriv (Single a)   =  Single 0
  > deriv (Cons a as)  =  deriv' as 1
  >                       where deriv' (Single a)  n  =  Single (n * a)
  >                             deriv' (Cons a as) n  =  Cons (n * a) (deriv' as (n+1))
  
  Side note: we cannot in general implement a Boolean equality test for
  PowerSeries. For example, we know that `deriv test0` equals `test1`
  but we cannot compute `True` in finite time by comparing the
  coefficients of the two power series.
  