{-# 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 :: forall a (v :: * -> *). (Eq a, Field a, Vector v a) => [Poly v a]
legendre = (Poly v a -> Poly v a) -> [Poly v a] -> [Poly v a]
forall a b. (a -> b) -> [a] -> [b]
map (Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Poly v a -> Poly v a -> Poly v a
`subst'` v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [a
1 a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
2, a
1 a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
2]) [Poly v a]
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' = Poly v a -> Poly v a -> Poly v a
forall a (v :: * -> *) (w :: * -> *).
(Eq a, Semiring a, Vector v a, Vector w a) =>
Poly v a -> Poly w a -> Poly w a
subst
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 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [-a
1, a
2] Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (a -> Poly v a -> Poly v a -> Poly v a)
-> [a] -> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 a -> Poly v a -> Poly v a -> Poly v a
forall {a} {v :: * -> *}.
(Eq a, Euclidean a, Vector v a, Ring a) =>
a -> Poly v a -> Poly v a -> Poly v a
rec ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [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 = Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Euclidean a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
unscale' Word
0 (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1) (v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [-a
1 a -> a -> a
forall a. Ring a => a -> a -> a
- a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n, a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
4 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n] Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
* Poly v a
p Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Word -> a -> Poly v a -> Poly v 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)
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 = a -> a -> [Poly v a]
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 a -> a -> a
forall a. Ring a => a -> a -> a
- a
1 a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
2
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 = v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [(a
a a -> a -> a
forall a. Ring a => a -> a -> a
- a
b) a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
2, (a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
2) a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
2]
xs :: [Poly v a]
xs = Poly v a
x0 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: Poly v a
x1 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (a -> Poly v a -> Poly v a -> Poly v a)
-> [a] -> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 a -> Poly v a -> Poly v a -> Poly v a
forall {v :: * -> *}.
Vector v a =>
a -> Poly v a -> Poly v a -> Poly v a
rec ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1) a
2) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [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 = v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [a
d, a
c] Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
* Poly v a
p Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Word -> a -> Poly v a -> Poly v 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 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b) a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b a -> a -> a
forall a. Ring a => a -> a -> a
- a
2)
q :: a
q = (a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b a -> a -> a
forall a. Ring a => a -> a -> a
- a
1) a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
cp1
c :: a
c = a
q a -> a -> a
forall a. Semiring a => a -> a -> a
* ((a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b) a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b a -> a -> a
forall a. Ring a => a -> a -> a
- a
2))
d :: a
d = a
q a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
a a -> a -> a
forall a. Semiring a => a -> a -> a
* a
a a -> a -> a
forall a. Ring a => a -> a -> a
- a
b a -> a -> a
forall a. Semiring a => a -> a -> a
* a
b)
cm1 :: a
cm1 = a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Ring a => a -> a -> a
- a
1) a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b a -> a -> a
forall a. Ring a => a -> a -> a
- a
1) a -> a -> a
forall a. Semiring a => a -> a -> a
* (a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
b) a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` a
cp1
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 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: Word -> a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
1 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (Poly v a -> Poly v a -> Poly v a)
-> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Poly v a
pm1 Poly v a
p -> Word -> a -> Poly v a -> Poly v a
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 Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Poly v a
pm1) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [a] -> [a]
tail [Poly v a]
xs)
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 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: Word -> a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
2 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (Poly v a -> Poly v a -> Poly v a)
-> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Poly v a
pm1 Poly v a
p -> Word -> a -> Poly v a -> Poly v a
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 Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Poly v a
pm1) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [a] -> [a]
tail [Poly v a]
xs)
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 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: Word -> a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
1 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (a -> Poly v a -> Poly v a -> Poly v a)
-> [a] -> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 a -> Poly v a -> Poly v a -> Poly v a
forall {a} {v :: * -> *}.
(Eq a, Vector v a, Ring a) =>
a -> Poly v a -> Poly v a -> Poly v a
rec ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [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 = Word -> a -> Poly v a -> Poly v a
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 Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Word -> a -> Poly v a -> Poly v 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
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 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: Word -> a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a
monomial Word
1 a
2 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (a -> Poly v a -> Poly v a -> Poly v a)
-> [a] -> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 a -> Poly v a -> Poly v a -> Poly v a
forall {a} {v :: * -> *}.
(Eq a, Vector v a, Ring a) =>
a -> Poly v a -> Poly v a -> Poly v a
rec ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [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 = Word -> a -> Poly v a -> Poly v a
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 Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 (a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n) Poly v a
pm1
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 = a -> [Poly v a]
forall a (v :: * -> *).
(Eq a, Field a, Vector v a) =>
a -> [Poly v a]
laguerreGen a
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 = v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [a
1 a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a, -a
1]
xs :: [Poly v a]
xs = Poly v a
x0 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: Poly v a
x1 Poly v a -> [Poly v a] -> [Poly v a]
forall a. a -> [a] -> [a]
: (a -> Poly v a -> Poly v a -> Poly v a)
-> [a] -> [Poly v a] -> [Poly v a] -> [Poly v a]
forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3 a -> Poly v a -> Poly v a -> Poly v a
forall {v :: * -> *}.
Vector v a =>
a -> Poly v a -> Poly v a -> Poly v a
rec ((a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1) a
1) [Poly v a]
xs ([Poly v a] -> [Poly v a]
forall a. HasCallStack => [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 = v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
v a -> Poly v a
toPoly [(a
2 a -> a -> a
forall a. Semiring a => a -> a -> a
* a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1 a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a) a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1), -a
1 a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1)] Poly v a -> Poly v a -> Poly v a
forall a. Semiring a => a -> a -> a
* Poly v a
p Poly v a -> Poly v a -> Poly v a
forall a. Ring a => a -> a -> a
- Word -> a -> Poly v a -> Poly v a
forall a (v :: * -> *).
(Eq a, Semiring a, Vector v a) =>
Word -> a -> Poly v a -> Poly v a
scale Word
0 ((a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
a) a -> a -> a
forall a. Euclidean a => a -> a -> a
`quot` (a
n a -> a -> a
forall a. Semiring a => a -> a -> a
+ a
1)) Poly v a
pm1