{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE ScopedTypeVariables #-}

-- | Gaussian quadrature over [-1, 1] using the Legendre polynomials.
--
module Integration.Gaussian (
  jacobi_matrix )
where

import Algebra.Absolute ( abs )
import qualified Algebra.Algebraic as Algebraic ( C )
import qualified Algebra.Ring as Ring ( C )
import Data.Vector.Fixed ( Arity )
import NumericPrelude hiding ( abs )

import Linear.Matrix ( Mat(..), construct )
import Linear.QR ( eigenvectors_symmetric )


-- | Compute the Jacobi matrix for the Legendre polynomials over
--   [-1,1]. This matrix is derived from the three-term recurrence
--   relation which is expressable as a matrix equation. When the
--   polynomials are ortho/normal/, its Jacobi matrix should be
--   symmetric. But since the Legendre polynomials are not
--   ortho/normal/, we have to make a few modifications. The procedure
--   is described in Amparo, Gil, and Segura, \"Numerical Methods for
--   Special Functions\".
--
--   Examples:
--
--   >>> import Linear.Matrix ( Mat2, frobenius_norm, fromList )
--
--   >>> let c = 1/(sqrt 3)
--   >>> let expected = fromList [[0,c],[c, 0]] :: Mat2 Double
--   >>> let actual = jacobi_matrix :: Mat2 Double
--   >>> frobenius_norm (actual - expected) < 1e-14
--   True
--
jacobi_matrix :: forall m a. (Arity m, Algebraic.C a, Ring.C a)
              => Mat m m a
jacobi_matrix =
  construct lambda
  where
    coeff_a :: Int -> a
    coeff_a n = fromRational' $ (toRational (n + 1))/(toRational (2*n + 1))

    coeff_c :: Int -> a
    coeff_c n = fromRational' $ (toRational n)/(toRational (2*n + 1))

    max2 :: Ord b => b -> b -> b
    max2 x y = if x >= y then x else y

    lambda i j
      | j == i = fromInteger 0
      | abs (i-j) == 1 =
          let k = max2 i j
          in  sqrt $ (coeff_a (k-1)) * (coeff_c k)

      | otherwise = fromInteger 0