{-# 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 :: [Poly Rational]
legendres = Poly Rational
forall a. (Num a, Eq a) => Poly a
one Poly Rational -> [Poly Rational] -> [Poly Rational]
forall a. a -> [a] -> [a]
: Poly Rational
forall a. (Num a, Eq a) => Poly a
x Poly Rational -> [Poly Rational] -> [Poly Rational]
forall a. a -> [a] -> [a]
: 
    [ Poly Rational -> Poly Rational -> Poly Rational
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
multPoly
        (Endianness -> [Rational] -> Poly Rational
forall a. (Num a, Eq a) => Endianness -> [a] -> Poly a
poly Endianness
LE [Rational -> Rational
forall a. Fractional a => a -> a
recip (Rational
n' Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ 1)])
        (Poly Rational -> Poly Rational -> Poly Rational
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
addPoly (Endianness -> [Rational] -> Poly Rational
forall a. (Num a, Eq a) => Endianness -> [a] -> Poly a
poly Endianness
LE [0, 2 Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Rational
n' Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
+ 1] Poly Rational -> Poly Rational -> Poly Rational
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
`multPoly` Poly Rational
p_n)
                 (Endianness -> [Rational] -> Poly Rational
forall a. (Num a, Eq a) => Endianness -> [a] -> Poly a
poly Endianness
LE           [-Rational
n'] Poly Rational -> Poly Rational -> Poly Rational
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
`multPoly` Poly Rational
p_nm1)
        )
    | Integer
n     <- [1..], let n' :: Rational
n' = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger Integer
n
    | Poly Rational
p_n   <- [Poly Rational] -> [Poly Rational]
forall a. [a] -> [a]
tail [Poly Rational]
legendres
    | Poly Rational
p_nm1 <- [Poly Rational]
legendres
    ]

-- |Compute the coefficients of the n'th Legendre polynomial.
legendre :: (Fractional a, Eq a) => Int -> Poly a
legendre :: Int -> Poly a
legendre n :: Int
n = Endianness -> [a] -> Poly a
forall a. (Num a, Eq a) => Endianness -> [a] -> Poly a
poly Endianness
LE ([a] -> Poly a)
-> (Poly Rational -> [a]) -> Poly Rational -> Poly a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> a) -> [Rational] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Rational -> a
forall a. Fractional a => Rational -> a
fromRational ([Rational] -> [a])
-> (Poly Rational -> [Rational]) -> Poly Rational -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Endianness -> Poly Rational -> [Rational]
forall a. (Num a, Eq a) => Endianness -> Poly a -> [a]
polyCoeffs Endianness
LE (Poly Rational -> Poly a) -> Poly Rational -> Poly a
forall a b. (a -> b) -> a -> b
$ [Poly Rational]
legendres [Poly Rational] -> Int -> Poly Rational
forall a. [a] -> Int -> a
!! Int
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 :: Int -> a -> a
evalLegendre n :: Int
n t :: a
t = a -> [a]
forall a. Fractional a => a -> [a]
evalLegendres a
t [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
n

-- |Evaluate all the Legendre polynomials at a point X.
evalLegendres :: Fractional a => a -> [a]
evalLegendres :: a -> [a]
evalLegendres t :: a
t = [a]
ps
    where
       ps :: [a]
ps = 1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a
t a -> [a] -> [a]
forall a. a -> [a] -> [a]
: 
            [ ((2 a -> a -> a
forall a. Num a => a -> a -> a
* a
n a -> a -> a
forall a. Num a => a -> a -> a
+ 1) a -> a -> a
forall a. Num a => a -> a -> a
* a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
p_n a -> a -> a
forall a. Num a => a -> a -> a
- a
n a -> a -> a
forall a. Num a => a -> a -> a
* a
p_nm1) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
n a -> a -> a
forall a. Num a => a -> a -> a
+ 1)
            | a
n     <- (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (1a -> a -> a
forall a. Num a => a -> a -> a
+) 1
            | a
p_n   <- [a] -> [a]
forall a. [a] -> [a]
tail [a]
ps
            | a
p_nm1 <- [a]
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 :: Int -> a -> (a, a)
evalLegendreDeriv 0 _ = (1,0)
evalLegendreDeriv n :: Int
n t :: a
t = case Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) (a -> [a]
forall a. Fractional a => a -> [a]
evalLegendres a
t) of
    (p2 :: a
p2:p1 :: a
p1:_)   -> (a
p1, Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n a -> a -> a
forall a. Num a => a -> a -> a
* (a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
p1 a -> a -> a
forall a. Num a => a -> a -> a
- a
p2) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
ta -> a -> a
forall a. Num a => a -> a -> a
*a
t a -> a -> a
forall a. Num a => a -> a -> a
- 1))
    _ -> [Char] -> (a, a)
forall a. HasCallStack => [Char] -> a
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 :: Int -> b -> [b]
legendreRoots n :: Int
n eps :: b
eps = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
forall a. Num a => a -> a
negate [b]
mRoots [b] -> [b] -> [b]
forall a. [a] -> [a] -> [a]
++ [b] -> [b]
forall a. [a] -> [a]
reverse (Int -> [b] -> [b]
forall a. Int -> [a] -> [a]
take (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
m) [b]
mRoots)
    where
        -- the roots are symmetric in the interval so we only have to find 'm' of them.
        -- The rest are reflections.
        m :: Int
m = (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` 2
        mRoots :: [b]
mRoots = [b -> b
improveRoot (Int -> b
forall b a. (Fractional b, Integral a) => a -> b
z0 Int
i) | Int
i <- [0..Int
mInt -> Int -> Int
forall a. Num a => a -> a -> a
-1]]
        
        -- Initial guess for i'th root of the n'th Legendre polynomial
        z0 :: a -> b
z0 i :: a
i = Double -> b
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Double -> Double
forall a. Floating a => a -> a
cos (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 0.75) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 0.5)) :: Double)
        -- Improve estimate of a root by newton's method
        improveRoot :: b -> b
improveRoot z1 :: b
z1
            | b -> b
forall a. Num a => a -> a
abs (b
z2b -> b -> b
forall a. Num a => a -> a -> a
-b
z1) b -> b -> Bool
forall a. Ord a => a -> a -> Bool
<= b
eps    = b
z2
            | Bool
otherwise             = b -> b
improveRoot b
z2
            where
                (y :: b
y, dy :: b
dy) = Int -> b -> (b, b)
forall a. Fractional a => Int -> a -> (a, a)
evalLegendreDeriv Int
n b
z1
                z2 :: b
z2 = b
z1 b -> b -> b
forall a. Num a => a -> a -> a
- b
yb -> b -> b
forall a. Fractional a => a -> a -> a
/b
dy