module Math.Polynomial.Legendre where
import Math.Polynomial
legendres :: [Poly Rational]
legendres = one : x :
[ multPoly
(poly LE [recip (n' + 1)])
(addPoly (poly LE [0, 2 * n' + 1] `multPoly` p_n)
(poly LE [n'] `multPoly` p_nm1)
)
| n <- [1..], let n' = fromInteger n
| p_n <- tail legendres
| p_nm1 <- legendres
]
legendre :: (Fractional a, Eq a) => Int -> Poly a
legendre n = poly LE . map fromRational . polyCoeffs LE $ legendres !! n
evalLegendre :: Fractional a => Int -> a -> a
evalLegendre n t = evalLegendres t !! n
evalLegendres :: Fractional a => a -> [a]
evalLegendres t = ps
where
ps = 1 : t :
[ ((2 * n + 1) * t * p_n n * p_nm1) / (n + 1)
| n <- iterate (1+) 1
| p_n <- tail ps
| p_nm1 <- ps
]
evalLegendreDeriv :: Fractional a => Int -> a -> (a,a)
evalLegendreDeriv 0 _ = (1,0)
evalLegendreDeriv n t = case drop (n1) (evalLegendres t) of
(p2:p1:_) -> (p1, fromIntegral n * (t * p1 p2) / (t*t 1))
_ -> error "evalLegendreDeriv: evalLegendres didn't return a long enough list"
legendreRoots :: (Fractional b, Ord b) => Int -> b -> [b]
legendreRoots n eps = map negate mRoots ++ reverse (take (nm) mRoots)
where
m = (n + 1) `div` 2
mRoots = [improveRoot (z0 i) | i <- [0..m1]]
z0 i = realToFrac (cos (pi * (fromIntegral i + 0.75) / (fromIntegral n + 0.5)) :: Double)
improveRoot z1
| abs (z2z1) <= eps = z2
| otherwise = improveRoot z2
where
(y, dy) = evalLegendreDeriv n z1
z2 = z1 y/dy