-- |
-- Module:      Data.Poly.Orthogonal
-- Copyright:   (c) 2019 Andrew Lelechenko
-- Licence:     BSD3
-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- Classical orthogonal polynomials.
--
-- @since 0.4.0.0

{-# 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)

-- | <https://en.wikipedia.org/wiki/Legendre_polynomials Legendre polynomials>.
--
-- >>> take 3 legendre :: [Data.Poly.VPoly Double]
-- [1.0,1.0 * X + 0.0,1.5 * X^2 + 0.0 * X + (-0.5)]
--
-- @since 0.4.0.0
legendre :: (Eq a, Field a, Vector v a) => [Poly v a]
legendre :: forall a (v :: * -> *). (Eq a, Field a, Vector v a) => [Poly v a]
legendre = forall a b. (a -> b) -> [a] -> [b]
map (forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Poly v a -> Poly v a -> Poly v a
`subst'` forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [a
1 forall a. Euclidean a => a -> a -> a
`quot` a
2, a
1 forall a. Euclidean a => a -> a -> a
`quot` a
2]) forall a (v :: * -> *).
(Eq a, Euclidean a, Ring a, Vector v a) =>
[Poly v a]
legendreShifted
  where
    subst' :: (Eq a, Semiring a, Vector v a) => Poly v a -> Poly v a -> Poly v a
    subst' :: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Poly v a -> Poly v a -> Poly v a
subst' = forall a (v :: * -> *) (w :: * -> *).
(Eq a, Semiring a, Vector v a, Vector w a) =>
Poly v a -> Poly w a -> Poly w a
subst

-- | <https://en.wikipedia.org/wiki/Legendre_polynomials#Shifted_Legendre_polynomials Shifted Legendre polynomials>.
--
-- >>> take 3 legendreShifted :: [Data.Poly.VPoly Integer]
-- [1,2 * X + (-1),6 * X^2 + (-6) * X + 1]
--
-- @since 0.4.0.0
legendreShifted :: (Eq a, Euclidean a, Ring a, Vector v a) => [Poly v a]
legendreShifted :: forall a (v :: * -> *).
(Eq a, Euclidean a, Ring a, Vector v a) =>
[Poly v a]
legendreShifted = [Poly v a]
xs
  where
    xs :: [Poly v a]
xs = Poly v a
1 forall a. a -> [a] -> [a]
: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [-a
1, a
2] forall a. a -> [a] -> [a]
: forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 forall {a} {v :: * -> *}.
(Eq a, Euclidean a, Vector v a, Ring a) =>
a -> Poly v a -> Poly v a -> Poly v a
rec (forall a. (a -> a) -> a -> [a]
iterate (forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)
    rec :: a -> Poly v a -> Poly v a -> Poly v a
rec a
n Poly v a
pm1 Poly v a
p = forall a (v :: * -> *).
(Eq a, Euclidean a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
unscale' Word
0 (a
n forall a. Semiring a => a -> a -> a
+ a
1) (forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [-a
1 forall a. Ring a => a -> a -> a
- a
2 forall a. Semiring a => a -> a -> a
* a
n, a
2 forall a. Semiring a => a -> a -> a
+ a
4 forall a. Semiring a => a -> a -> a
* a
n] forall a. Semiring a => a -> a -> a
* Poly v a
p forall a. Ring a => a -> a -> a
- forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 a
n Poly v a
pm1)

-- | <https://en.wikipedia.org/wiki/Gegenbauer_polynomials Gegenbauer polynomials>.
--
-- @since 0.4.0.0
gegenbauer :: (Eq a, Field a, Vector v a) => a -> [Poly v a]
gegenbauer :: forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
a -> [Poly v a]
gegenbauer a
g = forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
a -> a -> [Poly v a]
jacobi a
a a
a
  where
    a :: a
a = a
g forall a. Ring a => a -> a -> a
- a
1 forall a. Euclidean a => a -> a -> a
`quot` a
2

-- | <https://en.wikipedia.org/wiki/Jacobi_polynomials Jacobi polynomials>.
--
-- @since 0.4.0.0
jacobi :: (Eq a, Field a, Vector v a) => a -> a -> [Poly v a]
jacobi :: forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
a -> a -> [Poly v a]
jacobi a
a a
b = [Poly v a]
xs
  where
    x0 :: Poly v a
x0 = Poly v a
1
    x1 :: Poly v a
x1 = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [(a
a forall a. Ring a => a -> a -> a
- a
b) forall a. Euclidean a => a -> a -> a
`quot` a
2, (a
a forall a. Semiring a => a -> a -> a
+ a
b forall a. Semiring a => a -> a -> a
+ a
2) forall a. Euclidean a => a -> a -> a
`quot` a
2]
    xs :: [Poly v a]
xs = Poly v a
x0 forall a. a -> [a] -> [a]
: Poly v a
x1 forall a. a -> [a] -> [a]
: forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 forall {v :: * -> *}.
Vector v a =>
a -> Poly v a -> Poly v a -> Poly v a
rec (forall a. (a -> a) -> a -> [a]
iterate (forall a. Semiring a => a -> a -> a
+ a
1) a
2) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)
    rec :: a -> Poly v a -> Poly v a -> Poly v a
rec a
n Poly v a
pm1 Poly v a
p = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [a
d, a
c] forall a. Semiring a => a -> a -> a
* Poly v a
p forall a. Ring a => a -> a -> a
- forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 a
cm1 Poly v a
pm1
      where
        cp1 :: a
cp1 = a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
* (a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Semiring a => a -> a -> a
+ a
b) forall a. Semiring a => a -> a -> a
* (a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Semiring a => a -> a -> a
+ a
b forall a. Ring a => a -> a -> a
- a
2)
        q :: a
q   = (a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Semiring a => a -> a -> a
+ a
b forall a. Ring a => a -> a -> a
- a
1) forall a. Euclidean a => a -> a -> a
`quot` a
cp1
        c :: a
c   = a
q forall a. Semiring a => a -> a -> a
* ((a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Semiring a => a -> a -> a
+ a
b) forall a. Semiring a => a -> a -> a
* (a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Semiring a => a -> a -> a
+ a
b forall a. Ring a => a -> a -> a
- a
2))
        d :: a
d   = a
q forall a. Semiring a => a -> a -> a
* (a
a forall a. Semiring a => a -> a -> a
* a
a forall a. Ring a => a -> a -> a
- a
b forall a. Semiring a => a -> a -> a
* a
b)
        cm1 :: a
cm1 = a
2 forall a. Semiring a => a -> a -> a
* (a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Ring a => a -> a -> a
- a
1) forall a. Semiring a => a -> a -> a
* (a
n forall a. Semiring a => a -> a -> a
+ a
b forall a. Ring a => a -> a -> a
- a
1) forall a. Semiring a => a -> a -> a
* (a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
+ a
a forall a. Semiring a => a -> a -> a
+ a
b) forall a. Euclidean a => a -> a -> a
`quot` a
cp1

-- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials>
-- of the first kind.
--
-- >>> take 3 chebyshev1 :: [Data.Poly.VPoly Integer]
-- [1,1 * X + 0,2 * X^2 + 0 * X + (-1)]
--
-- @since 0.4.0.0
chebyshev1 :: (Eq a, Ring a, Vector v a) => [Poly v a]
chebyshev1 :: forall a (v :: * -> *). (Eq a, Ring a, Vector v a) => [Poly v a]
chebyshev1 = [Poly v a]
xs
  where
    xs :: [Poly v a]
xs = Poly v a
1 forall a. a -> [a] -> [a]
: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
1 forall a. a -> [a] -> [a]
: forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Poly v a
pm1 Poly v a
p -> forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
1 a
2 Poly v a
p forall a. Ring a => a -> a -> a
- Poly v a
pm1) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)

-- | <https://en.wikipedia.org/wiki/Chebyshev_polynomials Chebyshev polynomials>
-- of the second kind.
--
-- >>> take 3 chebyshev2 :: [Data.Poly.VPoly Integer]
-- [1,2 * X + 0,4 * X^2 + 0 * X + (-1)]
--
-- @since 0.4.0.0
chebyshev2 :: (Eq a, Ring a, Vector v a) => [Poly v a]
chebyshev2 :: forall a (v :: * -> *). (Eq a, Ring a, Vector v a) => [Poly v a]
chebyshev2 = [Poly v a]
xs
  where
    xs :: [Poly v a]
xs = Poly v a
1 forall a. a -> [a] -> [a]
: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
2 forall a. a -> [a] -> [a]
: forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Poly v a
pm1 Poly v a
p -> forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
1 a
2 Poly v a
p forall a. Ring a => a -> a -> a
- Poly v a
pm1) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)

-- | Probabilists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>.
--
-- >>> take 3 hermiteProb :: [Data.Poly.VPoly Integer]
-- [1,1 * X + 0,1 * X^2 + 0 * X + (-1)]
--
-- @since 0.4.0.0
hermiteProb :: (Eq a, Ring a, Vector v a) => [Poly v a]
hermiteProb :: forall a (v :: * -> *). (Eq a, Ring a, Vector v a) => [Poly v a]
hermiteProb = [Poly v a]
xs
  where
    xs :: [Poly v a]
xs = Poly v a
1 forall a. a -> [a] -> [a]
: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
1 forall a. a -> [a] -> [a]
: forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 forall {a} {v :: * -> *}.
(Eq a, Ring a, Vector v a) =>
a -> Poly v a -> Poly v a -> Poly v a
rec (forall a. (a -> a) -> a -> [a]
iterate (forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)
    rec :: a -> Poly v a -> Poly v a -> Poly v a
rec a
n Poly v a
pm1 Poly v a
p = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
1 a
1 Poly v a
p forall a. Ring a => a -> a -> a
- forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 a
n Poly v a
pm1

-- | Physicists' <https://en.wikipedia.org/wiki/Hermite_polynomials Hermite polynomials>.
--
-- >>> take 3 hermitePhys :: [Data.Poly.VPoly Double]
-- [1.0,2.0 * X + 0.0,4.0 * X^2 + 0.0 * X + (-2.0)]
--
-- @since 0.4.0.0
hermitePhys :: (Eq a, Ring a, Vector v a) => [Poly v a]
hermitePhys :: forall a (v :: * -> *). (Eq a, Ring a, Vector v a) => [Poly v a]
hermitePhys = [Poly v a]
xs
  where
    xs :: [Poly v a]
xs = Poly v a
1 forall a. a -> [a] -> [a]
: forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
2 forall a. a -> [a] -> [a]
: forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 forall {a} {v :: * -> *}.
(Eq a, Ring a, Vector v a) =>
a -> Poly v a -> Poly v a -> Poly v a
rec (forall a. (a -> a) -> a -> [a]
iterate (forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)
    rec :: a -> Poly v a -> Poly v a -> Poly v a
rec a
n Poly v a
pm1 Poly v a
p = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
1 a
2 Poly v a
p forall a. Ring a => a -> a -> a
- forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 (a
2 forall a. Semiring a => a -> a -> a
* a
n) Poly v a
pm1

-- | <https://en.wikipedia.org/wiki/Laguerre_polynomials Laguerre polynomials>.
--
-- >>> take 3 laguerre :: [Data.Poly.VPoly Double]
-- [1.0,(-1.0) * X + 1.0,0.5 * X^2 + (-2.0) * X + 1.0]
--
-- @since 0.4.0.0
laguerre :: (Eq a, Field a, Vector v a) => [Poly v a]
laguerre :: forall a (v :: * -> *). (Eq a, Field a, Vector v a) => [Poly v a]
laguerre = forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
a -> [Poly v a]
laguerreGen a
0

-- | <https://en.wikipedia.org/wiki/Laguerre_polynomials#Generalized_Laguerre_polynomials Generalized Laguerre polynomials>
--
-- @since 0.4.0.0
laguerreGen :: (Eq a, Field a, Vector v a) => a -> [Poly v a]
laguerreGen :: forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
a -> [Poly v a]
laguerreGen a
a = [Poly v a]
xs
  where
    x0 :: Poly v a
x0 = Poly v a
1
    x1 :: Poly v a
x1 = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [a
1 forall a. Semiring a => a -> a -> a
+ a
a, -a
1]
    xs :: [Poly v a]
xs = Poly v a
x0 forall a. a -> [a] -> [a]
: Poly v a
x1 forall a. a -> [a] -> [a]
: forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 forall {v :: * -> *}.
Vector v a =>
a -> Poly v a -> Poly v a -> Poly v a
rec (forall a. (a -> a) -> a -> [a]
iterate (forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs (forall a. [a] -> [a]
tail [Poly v a]
xs)
    rec :: a -> Poly v a -> Poly v a -> Poly v a
rec a
n Poly v a
pm1 Poly v a
p = forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [(a
2 forall a. Semiring a => a -> a -> a
* a
n forall a. Semiring a => a -> a -> a
+ a
1 forall a. Semiring a => a -> a -> a
+ a
a) forall a. Euclidean a => a -> a -> a
`quot` (a
n forall a. Semiring a => a -> a -> a
+ a
1), -a
1 forall a. Euclidean a => a -> a -> a
`quot` (a
n forall a. Semiring a => a -> a -> a
+ a
1)] forall a. Semiring a => a -> a -> a
* Poly v a
p forall a. Ring a => a -> a -> a
- forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 ((a
n forall a. Semiring a => a -> a -> a
+ a
a) forall a. Euclidean a => a -> a -> a
`quot` (a
n forall a. Semiring a => a -> a -> a
+ a
1)) Poly v a
pm1