module Numeric.Integration.IntegratePolynomialOnSimplex where
import Algebra.Ring ( C )
import Data.List ( transpose )
import Data.Matrix ( detLU
, fromLists
)
import Math.Algebra.Hspray ( (*^)
, Spray
, (^+^)
, bombieriSpray
, composeSpray
, constantSpray
, lone
, toList
)
integratePolynomialOnSimplex
:: (C a, Fractional a, Ord a)
=> Spray a
-> [[a]]
-> a
integratePolynomialOnSimplex :: forall a. (C a, Fractional a, Ord a) => Spray a -> [[a]] -> a
integratePolynomialOnSimplex Spray a
p [[a]]
simplex =
a
s forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
abs (forall a. (Ord a, Fractional a) => Matrix a -> a
detLU forall a b. (a -> b) -> a -> b
$ forall a. [[a]] -> Matrix a
fromLists [[a]]
b) forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int
2 .. Int
n])
where
v :: [a]
v = forall a. [a] -> a
last [[a]]
simplex
n :: Int
n = forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
v
b :: [[a]]
b = forall a b. (a -> b) -> [a] -> [b]
map (\[a]
column -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [a]
column [a]
v) (forall a. Int -> [a] -> [a]
take Int
n [[a]]
simplex)
vb :: [(a, [a])]
vb = forall a b. [a] -> [b] -> [(a, b)]
zip [a]
v (forall a. [[a]] -> [[a]]
transpose [[a]]
b)
variables :: [Spray a]
variables = forall a b. (a -> b) -> [a] -> [b]
map forall a. C a => Int -> Spray a
lone [Int
1 .. Int
n]
newvariables :: [Spray a]
newvariables = forall a b. (a -> b) -> [a] -> [b]
map
(\(a
vi, [a]
bi) ->
(forall a. (C a, Eq a) => a -> Spray a
constantSpray a
vi) forall a. (C a, Eq a) => Spray a -> Spray a -> Spray a
^+^ forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall a. (C a, Eq a) => Spray a -> Spray a -> Spray a
(^+^) (forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. (C a, Eq a) => a -> Spray a -> Spray a
(*^) [a]
bi [Spray a]
variables)
)
[(a, [a])]
vb
q :: Spray a
q = forall a. (C a, Eq a) => Spray a -> [Spray a] -> Spray a
composeSpray Spray a
p [Spray a]
newvariables
qterms :: [([Int], a)]
qterms = forall a. Spray a -> [([Int], a)]
toList forall a b. (a -> b) -> a -> b
$ forall a. C a => Spray a -> Spray a
bombieriSpray Spray a
q
s :: a
s = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall {b} {t :: * -> *}.
(Fractional b, Foldable t) =>
(t Int, b) -> b
f [([Int], a)]
qterms
where
f :: (t Int, b) -> b
f (t Int
exponents, b
coef) = if Int
d forall a. Eq a => a -> a -> Bool
== Int
0
then b
coef
else b
coef forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Int
n forall a. Num a => a -> a -> a
+ Int
1 .. Int
n forall a. Num a => a -> a -> a
+ Int
d])
where d :: Int
d = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum t Int
exponents