{-# LANGUAGE ScopedTypeVariables #-}

-- | The polynomials over [a,b] form a basis for L_2[a,b]. But often
--   the \"obvious\" choice of a basis {1, x, x^2,...} is not
--   convenient, because its elements are not orthogonal.
--
--   There are several sets of polynomials that are known to be
--   orthogonal over various intervals. Here we provide formulae for a
--   few popular sets (bases).
--
--   First we present the simple formulae that generates polynomials
--   orthogonal over a specific interval; say, [-1, 1] (these are what
--   you'll find in a textbook). Then, we present a \"prime\" version
--   that produces polynomials orthogonal over an arbitrary interval
--   [a, b]. This is accomplished via an affine transform which shifts
--   and scales the variable but preserves orthogonality.
--
--   The choice of basis depends not only on the interval, but on the
--   weight function. The inner product used is not necessarily the
--   standard one in L_2, rather,
--
--     \<f, g\> = \int_{a}^{b} w(x)*f(x)*g(x) dx
--
--   where w(x) is some non-negative (or non-positive) weight function.
--
module Polynomials.Orthogonal (
  legendre )
where

import NumericPrelude
import qualified Algebra.RealField as RealField ( C )
import Prelude ()


-- | The @n@th Legendre polynomial in @x@ over [-1,1]. These are
--   orthogonal over [-1,1] with the weight function w(x)=1.
--
--   Examples:
--
--   >>> let points = [0, 1, 2, 3, 4, 5] :: [Double]
--
--   The zeroth polynomial is the constant function p(x)=1.
--
--   >>> map (legendre 0) points
--   [1.0,1.0,1.0,1.0,1.0,1.0]
--
--   The first polynomial is the identity:
--
--   >>> map (legendre 1) points
--   [0.0,1.0,2.0,3.0,4.0,5.0]
--
--   The second polynomial should be (3x^2 - 1)/2:
--
--   >>> let actual = map (legendre 2) points
--   >>> let f x = (3*x^2 - 1)/2
--   >>> let expected = map f points
--   >>> actual == expected
--   True
--
--   The third polynomial should be (5x^3 - 3x)/2:
--
--   >>> let actual = map (legendre 3) points
--   >>> let f x = (5*x^3 - 3*x)/2
--   >>> let expected = map f points
--   >>> actual == expected
--   True
--
legendre :: forall a. (RealField.C a)
         => Integer -- ^ The degree (i.e. the number of) the polynomial you want.
         -> a -- ^ The dependent variable @x@.
         -> a
legendre 0 _ = fromInteger 1
legendre 1 x = x
legendre n x =
  (c1*x * (legendre (n-1) x) - c2*(legendre (n-2) x)) / (fromInteger n)
  where
    c1 = fromInteger $ 2*n - 1 :: a
    c2 = fromInteger $ n - 1 :: a