{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE RebindableSyntax #-}
module Data.Poly.Orthogonal
( legendre
, legendreShifted
, gegenbauer
, jacobi
, chebyshev1
, chebyshev2
, hermiteProb
, hermitePhys
, laguerre
, laguerreGen
) where
import Prelude hiding (quot, Num(..), fromIntegral)
import Data.Euclidean
import Data.Semiring
import Data.Poly.Semiring
import Data.Poly.Internal.Dense (unscale')
import Data.Vector.Generic (Vector, fromListN)
legendre :: (Eq a, Field a, Vector v a) => [Poly v a]
legendre = map (flip subst' (toPoly [1 `quot` 2, 1 `quot` 2])) legendreShifted
where
subst' :: (Eq a, Semiring a, Vector v a) => Poly v a -> Poly v a -> Poly v a
subst' = subst
legendreShifted :: (Eq a, Euclidean a, Ring a, Vector v a) => [Poly v a]
legendreShifted = xs
where
xs = 1 : toPoly [-1, 2] : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
rec n pm1 p = unscale' 0 (n + 1) (toPoly [-1 - 2 * n, 2 + 4 * n] * p - scale 0 n pm1)
gegenbauer :: (Eq a, Field a, Vector v a) => a -> [Poly v a]
gegenbauer g = jacobi a a
where
a = g - 1 `quot` 2
jacobi :: (Eq a, Field a, Vector v a) => a -> a -> [Poly v a]
jacobi a b = xs
where
x0 = 1
x1 = toPoly [(a - b) `quot` 2, (a + b + 2) `quot` 2]
xs = x0 : x1 : zipWith3 rec (iterate (+ 1) 2) xs (tail xs)
rec n pm1 p = toPoly [d, c] * p - scale 0 cm1 pm1
where
cp1 = 2 * n * (n + a + b) * (2 * n + a + b - 2)
q = (2 * n + a + b - 1) `quot` cp1
c = q * ((2 * n + a + b) * (2 * n + a + b - 2))
d = q * (a * a - b * b)
cm1 = 2 * (n + a - 1) * (n + b - 1) * (2 * n + a + b) `quot` cp1
chebyshev1 :: (Eq a, Ring a, Vector v a) => [Poly v a]
chebyshev1 = xs
where
xs = 1 : monomial 1 1 : zipWith (\pm1 p -> scale 1 2 p - pm1) xs (tail xs)
chebyshev2 :: (Eq a, Ring a, Vector v a) => [Poly v a]
chebyshev2 = xs
where
xs = 1 : monomial 1 2 : zipWith (\pm1 p -> scale 1 2 p - pm1) xs (tail xs)
hermiteProb :: (Eq a, Ring a, Vector v a) => [Poly v a]
hermiteProb = xs
where
xs = 1 : monomial 1 1 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
rec n pm1 p = scale 1 1 p - scale 0 n pm1
hermitePhys :: (Eq a, Ring a, Vector v a) => [Poly v a]
hermitePhys = xs
where
xs = 1 : monomial 1 2 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
rec n pm1 p = scale 1 2 p - scale 0 (2 * n) pm1
laguerre :: (Eq a, Field a, Vector v a) => [Poly v a]
laguerre = laguerreGen 0
laguerreGen :: (Eq a, Field a, Vector v a) => a -> [Poly v a]
laguerreGen a = xs
where
x0 = 1
x1 = toPoly [1 + a, -1]
xs = x0 : x1 : zipWith3 rec (iterate (+ 1) 1) xs (tail xs)
rec n pm1 p = toPoly [(2 * n + 1 + a) `quot` (n + 1), -1 `quot` (n + 1)] * p - scale 0 ((n + a) `quot` (n + 1)) pm1