module Math.Polynomial.Hermite where

import Math.Polynomial
import Data.VectorSpace

probHermite :: [Poly Integer]
probHermite :: [Poly Integer]
probHermite 
    = Poly Integer
forall a. (Num a, Eq a) => Poly a
one
    Poly Integer -> [Poly Integer] -> [Poly Integer]
forall a. a -> [a] -> [a]
: [ Poly Integer -> Poly Integer -> Poly Integer
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
multPoly Poly Integer
forall a. (Num a, Eq a) => Poly a
x Poly Integer
h_n Poly Integer -> Poly Integer -> Poly Integer
forall v. AdditiveGroup v => v -> v -> v
^-^ Poly Integer -> Poly Integer
forall a. (Num a, Eq a) => Poly a -> Poly a
polyDeriv Poly Integer
h_n
      | Poly Integer
h_n <- [Poly Integer]
probHermite
      ]

physHermite :: [Poly Integer]
physHermite :: [Poly Integer]
physHermite
    = Poly Integer
forall a. (Num a, Eq a) => Poly a
one
    Poly Integer -> [Poly Integer] -> [Poly Integer]
forall a. a -> [a] -> [a]
: [ Integer -> Poly Integer -> Poly Integer
forall a. (Num a, Eq a) => a -> Poly a -> Poly a
scalePoly 2 (Poly Integer -> Poly Integer -> Poly Integer
forall a. (Num a, Eq a) => Poly a -> Poly a -> Poly a
multPoly Poly Integer
forall a. (Num a, Eq a) => Poly a
x Poly Integer
h_n) Poly Integer -> Poly Integer -> Poly Integer
forall v. AdditiveGroup v => v -> v -> v
^-^ Poly Integer -> Poly Integer
forall a. (Num a, Eq a) => Poly a -> Poly a
polyDeriv Poly Integer
h_n
      | Poly Integer
h_n <- [Poly Integer]
physHermite
      ]

evalProbHermite :: (Integral a, Num b) => a -> b -> b
evalProbHermite :: a -> b -> b
evalProbHermite n :: a
n = (b, b) -> b
forall a b. (a, b) -> a
fst ((b, b) -> b) -> (b -> (b, b)) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> b -> (b, b)
forall a b. (Integral a, Num b) => a -> b -> (b, b)
evalProbHermiteDeriv a
n

evalProbHermiteDeriv :: (Integral a, Num b) => a -> b -> (b, b)
evalProbHermiteDeriv :: a -> b -> (b, b)
evalProbHermiteDeriv 0 _ = (1, 0)
evalProbHermiteDeriv 1 x :: b
x = (b
x, 1)
evalProbHermiteDeriv n :: a
n x :: b
x
    | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0     = [Char] -> (b, b)
forall a. HasCallStack => [Char] -> a
error "evalProbHermite: n < 0"
    | Bool
otherwise = a -> b -> b -> (b, b)
loop 1 b
x 1
    where
        loop :: a -> b -> b -> (b, b)
loop k :: a
k h_k :: b
h_k h_km1 :: b
h_km1
            | a
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
n    = (b
h_k, b
k' b -> b -> b
forall a. Num a => a -> a -> a
* b
h_km1)
            | Bool
otherwise = a -> b -> b -> (b, b)
loop (a
ka -> a -> a
forall a. Num a => a -> a -> a
+1) (b
x b -> b -> b
forall a. Num a => a -> a -> a
* b
h_k b -> b -> b
forall a. Num a => a -> a -> a
- b
k' b -> b -> b
forall a. Num a => a -> a -> a
* b
h_km1) b
h_k
            where k' :: b
k' = a -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k

evalPhysHermite :: (Integral a, Num b) => a -> b -> b
evalPhysHermite :: a -> b -> b
evalPhysHermite n :: a
n = (b, b) -> b
forall a b. (a, b) -> a
fst ((b, b) -> b) -> (b -> (b, b)) -> b -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> b -> (b, b)
forall a b. (Integral a, Num b) => a -> b -> (b, b)
evalPhysHermiteDeriv a
n

evalPhysHermiteDeriv :: (Integral a, Num b) => a -> b -> (b,b)
evalPhysHermiteDeriv :: a -> b -> (b, b)
evalPhysHermiteDeriv 0 _ = (1,   0)
evalPhysHermiteDeriv 1 x :: b
x = (2b -> b -> b
forall a. Num a => a -> a -> a
*b
x, 2)
evalPhysHermiteDeriv n :: a
n x :: b
x
    | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0     = [Char] -> (b, b)
forall a. HasCallStack => [Char] -> a
error "evalProbHermite: n < 0"
    | Bool
otherwise = a -> b -> b -> (b, b)
loop 1 (2b -> b -> b
forall a. Num a => a -> a -> a
*b
x) 1
    where
        loop :: a -> b -> b -> (b, b)
loop k :: a
k h_k :: b
h_k h_km1 :: b
h_km1
            | a
k a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
n    = (b
h_k, 2 b -> b -> b
forall a. Num a => a -> a -> a
* b
k' b -> b -> b
forall a. Num a => a -> a -> a
* b
h_km1)
            | Bool
otherwise = a -> b -> b -> (b, b)
loop (a
ka -> a -> a
forall a. Num a => a -> a -> a
+1) (2 b -> b -> b
forall a. Num a => a -> a -> a
* (b
x b -> b -> b
forall a. Num a => a -> a -> a
* b
h_k b -> b -> b
forall a. Num a => a -> a -> a
- b
k' b -> b -> b
forall a. Num a => a -> a -> a
* b
h_km1)) b
h_k
            where k' :: b
k' = a -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
k