hypergeomatrix 1.1.0.1 → 1.1.0.2
raw patch · 11 files changed
+596/−593 lines, 11 filesdep ~arraydep ~containersdep ~cyclotomicsetup-changedPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: array, containers, cyclotomic, tasty, tasty-hunit
API changes (from Hackage documentation)
Files
- CHANGELOG.md +18/−12
- LICENSE +29/−29
- README.md +123/−123
- Setup.hs +2/−2
- hypergeomatrix.cabal +50/−55
- src/Math/HypergeoMatrix.hs +2/−2
- src/Math/HypergeoMatrix/Gaussian.hs +7/−7
- src/Math/HypergeoMatrix/HypergeoMatrix.hs +160/−158
- src/Math/HypergeoMatrix/Internal.hs +151/−150
- tests/Approx.hs +8/−8
- tests/Main.hs +46/−47
CHANGELOG.md view
@@ -1,12 +1,18 @@-1.0.0.0---------* initial release--1.1.0.0---------* upgrade version bounds of the 'containers' dependency-* fixed LaTeX code in README--1.1.0.1---------* removed the upper bounds+1.0.0.0 +------- +* initial release + +1.1.0.0 +------- +* upgraded version bounds of the 'containers' dependency +* fixed LaTeX code in README + +1.1.0.1 +------- +* removed the upper bounds + +1.1.0.2 +------- +* added the upper bounds +* added the 'TypeOperators' extension +* use of pattern matching to avoid warning regarding the use of `head`
LICENSE view
@@ -1,29 +1,29 @@-BSD 3-Clause License--Copyright (c) 2022, Stéphane Laurent-All rights reserved.--Redistribution and use in source and binary forms, with or without-modification, are permitted provided that the following conditions are met:--1. Redistributions of source code must retain the above copyright notice, this- list of conditions and the following disclaimer.--2. Redistributions in binary form must reproduce the above copyright notice,- this list of conditions and the following disclaimer in the documentation- and/or other materials provided with the distribution.--3. Neither the name of the copyright holder nor the names of its- contributors may be used to endorse or promote products derived from- this software without specific prior written permission.--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"-AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE-IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE-DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE-FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL-DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR-SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER-CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,-OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+BSD 3-Clause License + +Copyright (c) 2022, Stéphane Laurent +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
README.md view
@@ -1,123 +1,123 @@-# hypergeomatrix--## Evaluation of the hypergeometric function of a matrix argument (Koev & Edelman's algorithm)--Let $(a_1, \ldots, a_p)$ and $(b_1, \ldots, b_q)$ be two vectors of real or -complex numbers, possibly empty, $\alpha > 0$ and $X$ a real symmetric or a -complex Hermitian matrix. -The corresponding *hypergeometric function of a matrix argument* is defined by --$${}_pF_q^{(\alpha)} \left(\begin{matrix} a_1, \ldots, a_p \\\\ b_1, \ldots, b_q\end{matrix}; X\right) = \sum_{k=0}^{\infty}\sum_{\kappa \vdash k} \frac{{(a_1)}_{\kappa}^{(\alpha)} \cdots {(a_p)}_{\kappa}^{(\alpha)}} {{(b_1)}_{\kappa}^{(\alpha)} \cdots {(b_q)}_{\kappa}^{(\alpha)}} \frac{C_{\kappa}^{(\alpha)}(X)}{k!}.$$--The inner sum is over the integer partitions $\kappa$ of $k$ (which we also -denote by $|\kappa| = k$). The symbol ${(\cdot)}_{\kappa}^{(\alpha)}$ is the -*generalized Pochhammer symbol*, defined by--$${(c)}^{(\alpha)}_{\kappa} = \prod_{i=1}^{\ell}\prod_{j=1}^{\kappa_i} \left(c - \frac{i-1}{\alpha} + j-1\right)$$--when $\kappa = (\kappa_1, \ldots, \kappa_\ell)$. -Finally, $C_{\kappa}^{(\alpha)}$ is a *Jack function*. -Given an integer partition $\kappa$ and $\alpha > 0$, and a -real symmetric or complex Hermitian matrix $X$ of order $n$, -the Jack function --$$C_{\kappa}^{(\alpha)}(X) = C_{\kappa}^{(\alpha)}(x_1, \ldots, x_n)$$--is a symmetric homogeneous polynomial of degree $|\kappa|$ in the -eigen values $x_1$, $\ldots$, $x_n$ of $X$. --The series defining the hypergeometric function does not always converge. -See the references for a discussion about the convergence. --The inner sum in the definition of the hypergeometric function is over -all partitions $\kappa \vdash k$ but actually -$C_{\kappa}^{(\alpha)}(X) = 0$ when $\ell(\kappa)$, the number of non-zero -entries of $\kappa$, is strictly greater than $n$.--For $\alpha=1$, $C_{\kappa}^{(\alpha)}$ is a *Schur polynomial* and it is -a *zonal polynomial* for $\alpha = 2$. -In random matrix theory, the hypergeometric function appears for $\alpha=2$ -and $\alpha$ is omitted from the notation, implicitely assumed to be $2$. --Koev and Edelman (2006) provided an efficient algorithm for the evaluation -of the truncated series --$$\sideset{_p^m}{_q^{(\alpha)}}F \left(\begin{matrix} a_1, \ldots, a_p \\\\ b_1, \ldots, b_q\end{matrix}; X\right) = \sum_{k=0}^{m}\sum_{\kappa \vdash k} \frac{{(a_1)}_{\kappa}^{(\alpha)} \cdots {(a_p)}_{\kappa}^{(\alpha)}} {{(b_1)}_{\kappa}^{(\alpha)} \cdots {(b_q)}_{\kappa}^{(\alpha)}} -\frac{C_{\kappa}^{(\alpha)}(X)}{k!}.$$--Hereafter, $m$ is called the *truncation weight of the summation* -(because $|\kappa|$ is called the weight of $\kappa$), the vector -$(a_1, \ldots, a_p)$ is called the vector of *upper parameters* while -the vector $(b_1, \ldots, b_q)$ is called the vector of *lower parameters*. -The user has to supply the vector $(x_1, \ldots, x_n)$ of the eigenvalues -of $X$. --For example, to compute--$$\sideset{_2^{15}}{_3^{(2)}}F \left(\begin{matrix} 3, 4 \\\\ 5, 6, 7\end{matrix}; 0.1, 0.4\right)$$--you have to enter --```haskell-hypergeomat 15 2 [3.0, 4.0], [5.0, 6.0, 7.0] [0.1, 0.4]-```--We said that the hypergeometric function is defined for a real symmetric -matrix or a complex Hermitian matrix $X$. Thus the eigenvalues of $X$ -are real. However we do not impose this restriction in `hypergeomatrix`. -The user can enter any list of real or complex numbers for the eigenvalues. --### Gaussian rational numbers--The library allows to use **Gaussian rational numbers**, i.e. complex numbers -with a rational real part and a rational imaginary part. The Gaussian rational -number $a + ib$ is obtained with `a +: b`, e.g. `(2%3) +: (5%2)`. The imaginary -unit usually denoted by $i$ is represented by `e(4)`:--```haskell-ghci> import Math.HypergeoMatrix-ghci> import Data.Ratio-ghci> alpha = 2%1-ghci> a = (2%7) +: (1%2)-ghci> b = (1%2) +: (0%1)-ghci> c = (2%1) +: (3%1)-ghci> x1 = (1%3) +: (1%4)-ghci> x2 = (1%5) +: (1%6)-ghci> hypergeomat 3 alpha [a, b] [c] [x1, x2]-26266543409/25159680000 + 155806638989/3698472960000*e(4)-```--### Univariate case--For $n = 1$, the hypergeometric function of a matrix argument is known as the -[generalized hypergeometric function](https://mathworld.wolfram.com/HypergeometricFunction.html). -It does not depend on $\alpha$. The case of $\sideset{_{2\thinspace}^{}}{_1^{}}F$ is the most known, -this is the Gauss hypergeometric function. Let's check a value. It is known that--$$\sideset{_{2\thinspace}^{}}{_1^{}}F \left(\begin{matrix} 1/4, 1/2 \\\\ 3/4\end{matrix}; 80/81\right) = 1.8.$$--Since $80/81$ is close to $1$, the convergence is slow. We compute the truncated series below -for $m = 300$.--```haskell-ghci> h <- hypergeomat 300 2 [1/4, 1/2] [3/4] [80/81]-ghci> h-1.7990026528192298-```---## References--- Plamen Koev and Alan Edelman. -*The efficient evaluation of the hypergeometric function of a matrix argument*.-Mathematics of computation, vol. 75, n. 254, 833-846, 2006.--- Robb Muirhead. -*Aspects of multivariate statistical theory*. -Wiley series in probability and mathematical statistics. -Probability and mathematical statistics. -John Wiley & Sons, New York, 1982.--- A. K. Gupta and D. K. Nagar. -*Matrix variate distributions*. -Chapman and Hall, 1999.+# hypergeomatrix + +## Evaluation of the hypergeometric function of a matrix argument (Koev & Edelman's algorithm) + +Let $(a_1, \ldots, a_p)$ and $(b_1, \ldots, b_q)$ be two vectors of real or +complex numbers, possibly empty, $\alpha > 0$ and $X$ a real symmetric or a +complex Hermitian matrix. +The corresponding *hypergeometric function of a matrix argument* is defined by + +$${}_pF_q^{(\alpha)} \left(\begin{matrix} a_1, \ldots, a_p \\\\ b_1, \ldots, b_q\end{matrix}; X\right) = \sum_{k=0}^{\infty}\sum_{\kappa \vdash k} \frac{{(a_1)}_{\kappa}^{(\alpha)} \cdots {(a_p)}_{\kappa}^{(\alpha)}} {{(b_1)}_{\kappa}^{(\alpha)} \cdots {(b_q)}_{\kappa}^{(\alpha)}} \frac{C_{\kappa}^{(\alpha)}(X)}{k!}.$$ + +The inner sum is over the integer partitions $\kappa$ of $k$ (which we also +denote by $|\kappa| = k$). The symbol ${(\cdot)}_{\kappa}^{(\alpha)}$ is the +*generalized Pochhammer symbol*, defined by + +$${(c)}^{(\alpha)}_{\kappa} = \prod_{i=1}^{\ell}\prod_{j=1}^{\kappa_i} \left(c - \frac{i-1}{\alpha} + j-1\right)$$ + +when $\kappa = (\kappa_1, \ldots, \kappa_\ell)$. +Finally, $C_{\kappa}^{(\alpha)}$ is a *Jack function*. +Given an integer partition $\kappa$ and $\alpha > 0$, and a +real symmetric or complex Hermitian matrix $X$ of order $n$, +the Jack function + +$$C_{\kappa}^{(\alpha)}(X) = C_{\kappa}^{(\alpha)}(x_1, \ldots, x_n)$$ + +is a symmetric homogeneous polynomial of degree $|\kappa|$ in the +eigen values $x_1$, $\ldots$, $x_n$ of $X$. + +The series defining the hypergeometric function does not always converge. +See the references for a discussion about the convergence. + +The inner sum in the definition of the hypergeometric function is over +all partitions $\kappa \vdash k$ but actually +$C_{\kappa}^{(\alpha)}(X) = 0$ when $\ell(\kappa)$, the number of non-zero +entries of $\kappa$, is strictly greater than $n$. + +For $\alpha=1$, $C_{\kappa}^{(\alpha)}$ is a *Schur polynomial* and it is +a *zonal polynomial* for $\alpha = 2$. +In random matrix theory, the hypergeometric function appears for $\alpha=2$ +and $\alpha$ is omitted from the notation, implicitely assumed to be $2$. + +Koev and Edelman (2006) provided an efficient algorithm for the evaluation +of the truncated series + +$$\sideset{_p^m}{_q^{(\alpha)}}F \left(\begin{matrix} a_1, \ldots, a_p \\\\ b_1, \ldots, b_q\end{matrix}; X\right) = \sum_{k=0}^{m}\sum_{\kappa \vdash k} \frac{{(a_1)}_{\kappa}^{(\alpha)} \cdots {(a_p)}_{\kappa}^{(\alpha)}} {{(b_1)}_{\kappa}^{(\alpha)} \cdots {(b_q)}_{\kappa}^{(\alpha)}} +\frac{C_{\kappa}^{(\alpha)}(X)}{k!}.$$ + +Hereafter, $m$ is called the *truncation weight of the summation* +(because $|\kappa|$ is called the weight of $\kappa$), the vector +$(a_1, \ldots, a_p)$ is called the vector of *upper parameters* while +the vector $(b_1, \ldots, b_q)$ is called the vector of *lower parameters*. +The user has to supply the vector $(x_1, \ldots, x_n)$ of the eigenvalues +of $X$. + +For example, to compute + +$$\sideset{_2^{15}}{_3^{(2)}}F \left(\begin{matrix} 3, 4 \\\\ 5, 6, 7\end{matrix}; 0.1, 0.4\right)$$ + +you have to enter + +```haskell +hypergeomat 15 2 [3.0, 4.0], [5.0, 6.0, 7.0] [0.1, 0.4] +``` + +We said that the hypergeometric function is defined for a real symmetric +matrix or a complex Hermitian matrix $X$. Thus the eigenvalues of $X$ +are real. However we do not impose this restriction in `hypergeomatrix`. +The user can enter any list of real or complex numbers for the eigenvalues. + +### Gaussian rational numbers + +The library allows to use **Gaussian rational numbers**, i.e. complex numbers +with a rational real part and a rational imaginary part. The Gaussian rational +number $a + ib$ is obtained with `a +: b`, e.g. `(2%3) +: (5%2)`. The imaginary +unit usually denoted by $i$ is represented by `e(4)`: + +```haskell +ghci> import Math.HypergeoMatrix +ghci> import Data.Ratio +ghci> alpha = 2%1 +ghci> a = (2%7) +: (1%2) +ghci> b = (1%2) +: (0%1) +ghci> c = (2%1) +: (3%1) +ghci> x1 = (1%3) +: (1%4) +ghci> x2 = (1%5) +: (1%6) +ghci> hypergeomat 3 alpha [a, b] [c] [x1, x2] +26266543409/25159680000 + 155806638989/3698472960000*e(4) +``` + +### Univariate case + +For $n = 1$, the hypergeometric function of a matrix argument is known as the +[generalized hypergeometric function](https://mathworld.wolfram.com/HypergeometricFunction.html). +It does not depend on $\alpha$. The case of $\sideset{_{2\thinspace}^{}}{_1^{}}F$ is the most known, +this is the Gauss hypergeometric function. Let's check a value. It is known that + +$$\sideset{_{2\thinspace}^{}}{_1^{}}F \left(\begin{matrix} 1/4, 1/2 \\\\ 3/4\end{matrix}; 80/81\right) = 1.8.$$ + +Since $80/81$ is close to $1$, the convergence is slow. We compute the truncated series below +for $m = 300$. + +```haskell +ghci> h <- hypergeomat 300 2 [1/4, 1/2] [3/4] [80/81] +ghci> h +1.7990026528192298 +``` + + +## References + +- Plamen Koev and Alan Edelman. +*The efficient evaluation of the hypergeometric function of a matrix argument*. +Mathematics of computation, vol. 75, n. 254, 833-846, 2006. + +- Robb Muirhead. +*Aspects of multivariate statistical theory*. +Wiley series in probability and mathematical statistics. +Probability and mathematical statistics. +John Wiley & Sons, New York, 1982. + +- A. K. Gupta and D. K. Nagar. +*Matrix variate distributions*. +Chapman and Hall, 1999.
Setup.hs view
@@ -1,2 +1,2 @@-import Distribution.Simple-main = defaultMain+import Distribution.Simple +main = defaultMain
hypergeomatrix.cabal view
@@ -1,55 +1,50 @@-cabal-version: 2.2-name: hypergeomatrix-version: 1.1.0.1-license: BSD-3-Clause-license-file: LICENSE-copyright: 2022 Stéphane Laurent-maintainer: laurent_step@outlook.fr-author: Stéphane Laurent-homepage: https://github.com/stla/hypergeomatrix#readme-synopsis: Hypergeometric function of a matrix argument-description:- Evaluation of hypergeometric functions of a matrix argument,- following Koev & Edelman's algorithm.--category: Math, Numeric-build-type: Simple-extra-source-files:- README.md- CHANGELOG.md--source-repository head- type: git- location: https://github.com/stla/hypergeomatrix--library- exposed-modules: Math.HypergeoMatrix- hs-source-dirs: src- other-modules:- Math.HypergeoMatrix.HypergeoMatrix- Math.HypergeoMatrix.Internal- Math.HypergeoMatrix.Gaussian-- default-language: Haskell2010- other-extensions:- BangPatterns DefaultSignatures ScopedTypeVariables TypeFamilies- TypeSynonymInstances-- ghc-options: -Wall- build-depends:- base >=4.7 && <5,- array >=0.5.4.0,- containers >=0.6.5.1,- cyclotomic >=1.1.1--test-suite unit-tests- type: exitcode-stdio-1.0- main-is: Main.hs- hs-source-dirs: tests/- other-modules: Approx- default-language: Haskell2010- build-depends:- base >=4.7 && <5,- tasty >=1.4.2.3,- tasty-hunit >=0.10.0.3,- hypergeomatrix -any+cabal-version: 2.2 +name: hypergeomatrix +version: 1.1.0.2 +synopsis: Hypergeometric function of a matrix argument +description: Evaluation of hypergeometric functions of a matrix argument, + following Koev & Edelman's algorithm. +homepage: https://github.com/stla/hypergeomatrix#readme +license: BSD-3-Clause +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2022 Stéphane Laurent +category: Math, Numeric +build-type: Simple +extra-source-files: README.md +extra-doc-files: CHANGELOG.md + +library + hs-source-dirs: src + exposed-modules: Math.HypergeoMatrix + other-modules: Math.HypergeoMatrix.HypergeoMatrix + , Math.HypergeoMatrix.Internal + , Math.HypergeoMatrix.Gaussian + build-depends: base >= 4.7 && < 5 + , array >= 0.5.4.0 && < 0.6 + , containers >= 0.6.4.1 && < 0.7 + , cyclotomic >= 1.1.1 && < 1.2 + other-extensions: BangPatterns + , DefaultSignatures + , ScopedTypeVariables + , TypeFamilies + , TypeOperators + , TypeSynonymInstances + default-language: Haskell2010 + ghc-options: -Wall + +test-suite unit-tests + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: tests/ + other-modules: Approx + Build-Depends: base >= 4.7 && < 5 + , tasty >= 1.4 && < 1.6 + , tasty-hunit >= 0.10 && < 0.11 + , hypergeomatrix + Default-Language: Haskell2010 + +source-repository head + type: git + location: https://github.com/stla/hypergeomatrix
src/Math/HypergeoMatrix.hs view
@@ -1,3 +1,3 @@-module Math.HypergeoMatrix (module X) where-import Math.HypergeoMatrix.Gaussian as X+module Math.HypergeoMatrix (module X) where +import Math.HypergeoMatrix.Gaussian as X import Math.HypergeoMatrix.HypergeoMatrix as X
src/Math/HypergeoMatrix/Gaussian.hs view
@@ -1,8 +1,8 @@-module Math.HypergeoMatrix.Gaussian - where-import Data.Complex.Cyclotomic--type GaussianRational = Cyclotomic--(+:) :: Rational -> Rational -> GaussianRational+module Math.HypergeoMatrix.Gaussian + where +import Data.Complex.Cyclotomic + +type GaussianRational = Cyclotomic + +(+:) :: Rational -> Rational -> GaussianRational (+:) = gaussianRat
src/Math/HypergeoMatrix/HypergeoMatrix.hs view
@@ -1,158 +1,160 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE ScopedTypeVariables #-}--module Math.HypergeoMatrix.HypergeoMatrix (hypergeomat) where-import Control.Monad (when)-import Data.Array hiding (index)-import Data.Array.IO hiding (index)-import Data.Sequence (Seq, index, update, (!?), (|>))-import qualified Data.Sequence as S-import Math.HypergeoMatrix.Internal --hypergeoI :: forall a. (Eq a, Fractional a, BaseFrac a)- => Int -> BaseFracType a -> [a] -> [a] -> Int -> a -> a-hypergeoI m alpha a b n x =- 1 + summation' 0 1 m []- where- summation' :: Fractional a => Int -> a -> Int -> [Int] -> a- summation' i z j kappa = go 1 z 0- where- go :: Int -> a -> a -> a- go kappai zz s- | i == 0 && kappai > j || i>0 && kappai > min (kappa!!(i-1)) j = s- | otherwise = go (kappai + 1) z' s''- where- kappa' = kappa ++ [kappai]- t = _T alpha a b (S.fromList $ filter (> 0) kappa') -- inutile de filtrer- z' = zz * x *- (fromIntegral (n-i) + inject alpha * (fromIntegral kappai-1)) * t- s' = if j > kappai && i <= n- then s + summation' (i+1) z' (j-kappai) kappa'- else s- s'' = s' + z'---summation :: forall a. (Fractional a, Eq a, BaseFrac a)- => [a] -> [a] -> [a] -> Seq (Maybe Int) -> Int -> BaseFracType a -> Int- -> a -> Int -> Seq Int -> IOArray (Int, Int) a -> IO a-summation a b x dico n alpha i z j kappa jarray- = if i == n- then- return 0- else do- let lkappa = kappa `index` (S.length kappa - 1)- let go :: Int -> a -> a -> IO a- go kappai !z' !s- | i == 0 && kappai > j || i > 0 && kappai > min lkappa j =- return s- | otherwise = do- let kappa' = kappa |> kappai- nkappa = _nkappa dico kappa'- z'' = z' * _T alpha a b kappa'- lkappa' = S.length kappa'- when (nkappa > 1 && (lkappa' == 1 || kappa' !? 1 == Just 0)) $ do- entry <- readArray jarray (nkappa - 1, 1)- let kap0m1' = fromIntegral (kappa' `index` 0 - 1)- newval = head x * (1 + inject alpha * kap0m1') * entry- writeArray jarray (nkappa, 1) newval- let go' :: Int -> IO ()- go' t- | t == n + 1 = return ()- | otherwise = do- _ <- jack alpha x dico 0 1 0 t kappa' jarray kappa' nkappa- go' (t + 1)- _ <- go' 2- entry' <- readArray jarray (nkappa, n)- let s' = s + z'' * entry'- if j > kappai && i <= n- then do- s'' <-- summation- a- b- x- dico- n- alpha- (i + 1)- z''- (j - kappai)- kappa'- jarray- go (kappai + 1) z'' (s' + s'')- else go (kappai + 1) z'' s'- go 1 z 0--jack :: (Fractional a, BaseFrac a)- => BaseFracType a -> [a] -> Seq (Maybe Int) -> Int -> a -> Int -> Int- -> Seq Int -> IOArray (Int, Int) a -> Seq Int -> Int -> IO ()-jack alpha x dico k beta c t mu jarray kappa nkappa = do- let i0 = max k 1- i1 = S.length (cleanPart mu) + 1- go :: Int -> IO ()- go i- | i == i1 = return ()- | otherwise- = do- let u = mu `index` (i - 1)- when (S.length mu == i || u > mu `index` i) $ do- let gamma = beta * _betaratio kappa mu i alpha- mu' = cleanPart $ update (i-1) (u - 1) mu- nmu = _nkappa dico mu'- if S.length mu' >= i && u > 1 -- "not (S.null mu')" useless because i>=1- then- jack alpha x dico i gamma (c + 1) t mu' jarray kappa nkappa- else- when (nkappa > 1) $ do- entry' <- readArray jarray (nkappa, t)- if not (S.null mu') -- any (> 0) mu'- then do- entry <- readArray jarray (nmu, t - 1)- writeArray- jarray- (nkappa, t)- (entry' + gamma * entry * x !! (t - 1) ^ (c + 1))- else writeArray- jarray- (nkappa, t)- (entry' + gamma * x !! (t - 1) ^ (c + 1))- go (i + 1)- _ <- go i0- entry1 <- readArray jarray (nkappa, t)- if k == 0- then- when (nkappa > 1) $ do- entry2 <- readArray jarray (nkappa, t - 1)- writeArray jarray (nkappa, t) (entry1 + entry2)- else do- entry2 <- readArray jarray (_nkappa dico mu, t - 1)- writeArray jarray (nkappa, t) (entry1 + beta * x !! (t - 1) ^ c * entry2)---- | Hypergeometric function of a matrix argument.--- Actually the matrix argument is given by the eigenvalues of the matrix.--- For a type `a` of real numbers, `BaseFracType a = a`. If `a = Complex b` --- is a type of complex numbers, then `BaseFracType a = b`. Thus `alpha` --- parameter cannot be a complex number.-hypergeomat :: forall a. (Eq a, Fractional a, BaseFrac a)- => Int -- ^ truncation weight- -> BaseFracType a -- ^ alpha parameter (usually 2)- -> [a] -- ^ upper parameters- -> [a] -- ^ lower parameters- -> [a] -- ^ variables (the eigenvalues)- -> IO a-hypergeomat m alpha a b x = do- let n = length x- if all (== head x) x- then- return $ hypergeoI m alpha a b n (head x)- else do- let pmn = _P m n- dico = _dico pmn m- xrange = [1 .. n]- line1 = zipWith (\i u -> ((1, i), u)) xrange (scanl1 (+) x)- otherlines = concatMap (\j -> [((j, i), 0) | i <- xrange]) [2 .. pmn]- arr0 =- array ((1, 1), (pmn, n)) (line1 ++ otherlines)- jarray <- thaw arr0- s <- summation a b x dico n alpha 0 1 m S.empty jarray- return $ s + 1+{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.HypergeoMatrix.HypergeoMatrix (hypergeomat) where +import Control.Monad (when) +import Data.Array hiding (index) +import Data.Array.IO hiding (index) +import Data.Sequence (Seq, index, update, (!?), (|>)) +import qualified Data.Sequence as S +import Math.HypergeoMatrix.Internal + +hypergeoI :: forall a. (Eq a, Fractional a, BaseFrac a) + => Int -> BaseFracType a -> [a] -> [a] -> Int -> a -> a +hypergeoI m alpha a b n x = + 1 + summation' 0 1 m [] + where + summation' :: Fractional a => Int -> a -> Int -> [Int] -> a + summation' i z j kappa = go 1 z 0 + where + go :: Int -> a -> a -> a + go kappai zz s + | i == 0 && kappai > j || i>0 && kappai > min (kappa!!(i-1)) j = s + | otherwise = go (kappai + 1) z' s'' + where + kappa' = kappa ++ [kappai] + t = _T alpha a b (S.fromList $ filter (> 0) kappa') -- inutile de filtrer + z' = zz * x * + (fromIntegral (n-i) + inject alpha * (fromIntegral kappai-1)) * t + s' = if j > kappai && i <= n + then s + summation' (i+1) z' (j-kappai) kappa' + else s + s'' = s' + z' + + +summation :: forall a. (Fractional a, Eq a, BaseFrac a) + => [a] -> [a] -> [a] -> Seq (Maybe Int) -> Int -> BaseFracType a -> Int + -> a -> Int -> Seq Int -> IOArray (Int, Int) a -> IO a +summation _ _ [] _ _ _ _ _ _ _ _ = error "summation: empty list of variables." +summation a b x@(x0:_) dico n alpha i z j kappa jarray + = if i == n + then + return 0 + else do + let lkappa = kappa `index` (S.length kappa - 1) + let go :: Int -> a -> a -> IO a + go kappai !z' !s + | i == 0 && kappai > j || i > 0 && kappai > min lkappa j = + return s + | otherwise = do + let kappa' = kappa |> kappai + nkappa = _nkappa dico kappa' + z'' = z' * _T alpha a b kappa' + lkappa' = S.length kappa' + when (nkappa > 1 && (lkappa' == 1 || kappa' !? 1 == Just 0)) $ do + entry <- readArray jarray (nkappa - 1, 1) + let kap0m1' = fromIntegral (kappa' `index` 0 - 1) + newval = x0 * (1 + inject alpha * kap0m1') * entry + writeArray jarray (nkappa, 1) newval + let go' :: Int -> IO () + go' t + | t == n + 1 = return () + | otherwise = do + _ <- jack alpha x dico 0 1 0 t kappa' jarray kappa' nkappa + go' (t + 1) + _ <- go' 2 + entry' <- readArray jarray (nkappa, n) + let s' = s + z'' * entry' + if j > kappai && i <= n + then do + s'' <- + summation + a + b + x + dico + n + alpha + (i + 1) + z'' + (j - kappai) + kappa' + jarray + go (kappai + 1) z'' (s' + s'') + else go (kappai + 1) z'' s' + go 1 z 0 + +jack :: (Fractional a, BaseFrac a) + => BaseFracType a -> [a] -> Seq (Maybe Int) -> Int -> a -> Int -> Int + -> Seq Int -> IOArray (Int, Int) a -> Seq Int -> Int -> IO () +jack alpha x dico k beta c t mu jarray kappa nkappa = do + let i0 = max k 1 + i1 = S.length (cleanPart mu) + 1 + go :: Int -> IO () + go i + | i == i1 = return () + | otherwise + = do + let u = mu `index` (i - 1) + when (S.length mu == i || u > mu `index` i) $ do + let gamma = beta * _betaratio kappa mu i alpha + mu' = cleanPart $ update (i-1) (u - 1) mu + nmu = _nkappa dico mu' + if S.length mu' >= i && u > 1 -- "not (S.null mu')" useless because i>=1 + then + jack alpha x dico i gamma (c + 1) t mu' jarray kappa nkappa + else + when (nkappa > 1) $ do + entry' <- readArray jarray (nkappa, t) + if not (S.null mu') -- any (> 0) mu' + then do + entry <- readArray jarray (nmu, t - 1) + writeArray + jarray + (nkappa, t) + (entry' + gamma * entry * x !! (t - 1) ^ (c + 1)) + else writeArray + jarray + (nkappa, t) + (entry' + gamma * x !! (t - 1) ^ (c + 1)) + go (i + 1) + _ <- go i0 + entry1 <- readArray jarray (nkappa, t) + if k == 0 + then + when (nkappa > 1) $ do + entry2 <- readArray jarray (nkappa, t - 1) + writeArray jarray (nkappa, t) (entry1 + entry2) + else do + entry2 <- readArray jarray (_nkappa dico mu, t - 1) + writeArray jarray (nkappa, t) (entry1 + beta * x !! (t - 1) ^ c * entry2) + +-- | Hypergeometric function of a matrix argument. +-- Actually the matrix argument is given by the eigenvalues of the matrix. +-- For a type \`a\` of real numbers, \`BaseFracType a = a\`. If \`a = Complex b\` +-- is a type of complex numbers, then \`BaseFracType a = b\`. Thus the \`alpha\` +-- parameter cannot be a complex number. +hypergeomat :: forall a. (Eq a, Fractional a, BaseFrac a) + => Int -- ^ truncation weight + -> BaseFracType a -- ^ alpha parameter (usually 2) + -> [a] -- ^ upper parameters + -> [a] -- ^ lower parameters + -> [a] -- ^ variables (the eigenvalues) + -> IO a +hypergeomat _ _ _ _ [] = error "hypergeomat: empty list of variables." +hypergeomat m alpha a b x@(x0:_) = do + let n = length x + if all (== x0) x + then + return $ hypergeoI m alpha a b n x0 + else do + let pmn = _P m n + dico = _dico pmn m + xrange = [1 .. n] + line1 = zipWith (\i u -> ((1, i), u)) xrange (scanl1 (+) x) + otherlines = concatMap (\j -> [((j, i), 0) | i <- xrange]) [2 .. pmn] + arr0 = + array ((1, 1), (pmn, n)) (line1 ++ otherlines) + jarray <- thaw arr0 + s <- summation a b x dico n alpha 0 1 m S.empty jarray + return $ s + 1
src/Math/HypergeoMatrix/Internal.hs view
@@ -1,150 +1,151 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE DefaultSignatures #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeSynonymInstances #-}--module Math.HypergeoMatrix.Internal where-import Data.Complex-import Data.Ratio-import Data.Maybe-import Data.Sequence (Seq ((:<|), (:|>), Empty), elemIndexL,- index, (!?), (><), (|>))-import qualified Data.Sequence as S-import Math.HypergeoMatrix.Gaussian--class BaseFrac a where- type family BaseFracType a- type BaseFracType a = a -- Default type family instance (unless overridden)- inject :: BaseFracType a -> a- default inject :: BaseFracType a ~ a => BaseFracType a -> a- inject = id--instance Integral a => BaseFrac (Ratio a)-instance BaseFrac Float-instance BaseFrac Double-instance BaseFrac GaussianRational where- type BaseFracType GaussianRational = Rational- inject x = x +: 0-instance Num a => BaseFrac (Complex a) where- type BaseFracType (Complex a) = a- inject x = x :+ 0--_diffSequence :: Seq Int -> Seq Int-_diffSequence (x :<| ys@(y :<| _)) = (x - y) :<| _diffSequence ys-_diffSequence x = x--_dualPartition :: Seq Int -> Seq Int-_dualPartition Empty = S.empty-_dualPartition xs = go 0 (_diffSequence xs) S.empty- where- go !i (d :<| ds) acc = go (i + 1) ds (d :<| acc)- go n Empty acc = finish n acc- finish !j (k :<| ks) = S.replicate k j >< finish (j - 1) ks- finish _ Empty = S.empty--_betaratio :: (Fractional a, BaseFrac a)- => Seq Int -> Seq Int -> Int -> BaseFracType a -> a-_betaratio kappa mu k alpha = alpha' * prod1 * prod2 * prod3- where- alpha' = inject alpha- t = fromIntegral k - alpha' * fromIntegral (mu `index` (k - 1))- ss = S.fromList [1 .. k - 1]- sss = ss |> k- u =- S.zipWith- (\s kap -> t + 1 - fromIntegral s + alpha' * fromIntegral kap)- sss (S.take k kappa)- v =- S.zipWith- (\s m -> t - fromIntegral s + alpha' * fromIntegral m)- ss (S.take (k - 1) mu)- l = mu `index` (k - 1) - 1- mu' = S.take l (_dualPartition mu)- w =- S.zipWith- (\s m -> fromIntegral m - t - alpha' * fromIntegral s)- (S.fromList [1 .. l]) mu'- prod1 = product $ fmap (\x -> x / (x + alpha' - 1)) u- prod2 = product $ fmap (\x -> (x + alpha') / x) v- prod3 = product $ fmap (\x -> (x + alpha') / x) w---_T :: (Fractional a, Eq a, BaseFrac a)- => BaseFracType a -> [a] -> [a] -> Seq Int -> a-_T alpha a b kappa- | S.null kappa || kappa !? 0 == Just 0 = 1- | prod1_den == 0 = 0- | otherwise = prod1_num/prod1_den * prod2 * prod3- where- alpha' = inject alpha- lkappa = S.length kappa - 1- kappai = kappa `index` lkappa- kappai' = fromIntegral kappai- i = fromIntegral lkappa- c = kappai' - 1 - i / alpha'- d = kappai' * alpha' - i - 1- s = fmap fromIntegral (S.fromList [1 .. kappai - 1])- kappa' = fromIntegral <$> S.take kappai (_dualPartition kappa)- e = S.zipWith (\x y -> d - x * alpha' + y) s kappa'- g = fmap (+ 1) e- s' = fmap fromIntegral (S.fromList [1 .. lkappa])- f = S.zipWith (\x y -> y * alpha' - x - d) s' (fmap fromIntegral kappa)- h = fmap (+ alpha') f- l = S.zipWith (*) h f- prod1_num = product (fmap (+ c) a)- prod1_den = product (fmap (+ c) b)- prod2 =- product $ S.zipWith (\x y -> (y - alpha') * x / y / (x + alpha')) e g- prod3 = product $ S.zipWith3 (\x y z -> (z - x) / (z + y)) f h l--a008284 :: [[Int]]-a008284 = [1] : f [[1]]- where- f xss = ys : f (ys : xss)- where- ys = map sum (zipWith take [1 ..] xss) ++ [1]--_P :: Int -> Int -> Int-_P m n = sum (concatMap (take (min m n)) (take m a008284))--_dico :: Int -> Int -> Seq (Maybe Int)-_dico pmn m = go False S.empty- where- go :: Bool -> Seq (Maybe Int) -> Seq (Maybe Int)- go k !d'- | k = d'- | otherwise = inner 0 [0] [m] [m] 0 d' Nothing- where- inner :: Int -> [Int] -> [Int] -> [Int] -> Int- -> Seq (Maybe Int) -> Maybe Int -> Seq (Maybe Int)- inner i !a !b !c !end !d !dlast- | dlast == Just pmn = go True d- | otherwise =- let bi = b !! i- in if bi > 0- then let l = min bi (c !! i)- in let ddlast = Just $ end + 1- in let dd = d |> ddlast- in let range1l = [1 .. l]- in inner- (i + 1)- (a ++ [end + 1 .. end + l])- (b ++ map (bi -) range1l)- (c ++ range1l)- (end + l)- dd- ddlast- else inner (i + 1) a b c end (d |> Nothing) Nothing--_nkappa :: Seq (Maybe Int) -> Seq Int -> Int-_nkappa dico (kappa0 :|> kappan) =- fromJust (dico `S.index` _nkappa dico kappa0) + kappan - 1-_nkappa _ Empty = 0--cleanPart :: Seq Int -> Seq Int-cleanPart kappa =- let i = elemIndexL 0 kappa- in if isJust i- then S.take (fromJust i) kappa- else kappa+{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE TypeOperators #-} + +module Math.HypergeoMatrix.Internal where +import Data.Complex +import Data.Ratio +import Data.Maybe +import Data.Sequence (Seq ((:<|), (:|>), Empty), elemIndexL, + index, (!?), (><), (|>)) +import qualified Data.Sequence as S +import Math.HypergeoMatrix.Gaussian + +class BaseFrac a where + type family BaseFracType a + type BaseFracType a = a -- Default type family instance (unless overridden) + inject :: BaseFracType a -> a + default inject :: BaseFracType a ~ a => BaseFracType a -> a + inject = id + +instance Integral a => BaseFrac (Ratio a) +instance BaseFrac Float +instance BaseFrac Double +instance BaseFrac GaussianRational where + type BaseFracType GaussianRational = Rational + inject x = x +: 0 +instance Num a => BaseFrac (Complex a) where + type BaseFracType (Complex a) = a + inject x = x :+ 0 + +_diffSequence :: Seq Int -> Seq Int +_diffSequence (x :<| ys@(y :<| _)) = (x - y) :<| _diffSequence ys +_diffSequence x = x + +_dualPartition :: Seq Int -> Seq Int +_dualPartition Empty = S.empty +_dualPartition xs = go 0 (_diffSequence xs) S.empty + where + go !i (d :<| ds) acc = go (i + 1) ds (d :<| acc) + go n Empty acc = finish n acc + finish !j (k :<| ks) = S.replicate k j >< finish (j - 1) ks + finish _ Empty = S.empty + +_betaratio :: (Fractional a, BaseFrac a) + => Seq Int -> Seq Int -> Int -> BaseFracType a -> a +_betaratio kappa mu k alpha = alpha' * prod1 * prod2 * prod3 + where + alpha' = inject alpha + t = fromIntegral k - alpha' * fromIntegral (mu `index` (k - 1)) + ss = S.fromList [1 .. k - 1] + sss = ss |> k + u = + S.zipWith + (\s kap -> t + 1 - fromIntegral s + alpha' * fromIntegral kap) + sss (S.take k kappa) + v = + S.zipWith + (\s m -> t - fromIntegral s + alpha' * fromIntegral m) + ss (S.take (k - 1) mu) + l = mu `index` (k - 1) - 1 + mu' = S.take l (_dualPartition mu) + w = + S.zipWith + (\s m -> fromIntegral m - t - alpha' * fromIntegral s) + (S.fromList [1 .. l]) mu' + prod1 = product $ fmap (\x -> x / (x + alpha' - 1)) u + prod2 = product $ fmap (\x -> (x + alpha') / x) v + prod3 = product $ fmap (\x -> (x + alpha') / x) w + + +_T :: (Fractional a, Eq a, BaseFrac a) + => BaseFracType a -> [a] -> [a] -> Seq Int -> a +_T alpha a b kappa + | S.null kappa || kappa !? 0 == Just 0 = 1 + | prod1_den == 0 = 0 + | otherwise = prod1_num/prod1_den * prod2 * prod3 + where + alpha' = inject alpha + lkappa = S.length kappa - 1 + kappai = kappa `index` lkappa + kappai' = fromIntegral kappai + i = fromIntegral lkappa + c = kappai' - 1 - i / alpha' + d = kappai' * alpha' - i - 1 + s = fmap fromIntegral (S.fromList [1 .. kappai - 1]) + kappa' = fromIntegral <$> S.take kappai (_dualPartition kappa) + e = S.zipWith (\x y -> d - x * alpha' + y) s kappa' + g = fmap (+ 1) e + s' = fmap fromIntegral (S.fromList [1 .. lkappa]) + f = S.zipWith (\x y -> y * alpha' - x - d) s' (fmap fromIntegral kappa) + h = fmap (+ alpha') f + l = S.zipWith (*) h f + prod1_num = product (fmap (+ c) a) + prod1_den = product (fmap (+ c) b) + prod2 = + product $ S.zipWith (\x y -> (y - alpha') * x / y / (x + alpha')) e g + prod3 = product $ S.zipWith3 (\x y z -> (z - x) / (z + y)) f h l + +a008284 :: [[Int]] +a008284 = [1] : f [[1]] + where + f xss = ys : f (ys : xss) + where + ys = map sum (zipWith take [1 ..] xss) ++ [1] + +_P :: Int -> Int -> Int +_P m n = sum (concatMap (take (min m n)) (take m a008284)) + +_dico :: Int -> Int -> Seq (Maybe Int) +_dico pmn m = go False S.empty + where + go :: Bool -> Seq (Maybe Int) -> Seq (Maybe Int) + go k !d' + | k = d' + | otherwise = inner 0 [0] [m] [m] 0 d' Nothing + where + inner :: Int -> [Int] -> [Int] -> [Int] -> Int + -> Seq (Maybe Int) -> Maybe Int -> Seq (Maybe Int) + inner i !a !b !c !end !d !dlast + | dlast == Just pmn = go True d + | otherwise = + let bi = b !! i + in if bi > 0 + then let l = min bi (c !! i) + in let ddlast = Just $ end + 1 + in let dd = d |> ddlast + in let range1l = [1 .. l] + in inner + (i + 1) + (a ++ [end + 1 .. end + l]) + (b ++ map (bi -) range1l) + (c ++ range1l) + (end + l) + dd + ddlast + else inner (i + 1) a b c end (d |> Nothing) Nothing + +_nkappa :: Seq (Maybe Int) -> Seq Int -> Int +_nkappa dico (kappa0 :|> kappan) = + fromJust (dico `S.index` _nkappa dico kappa0) + kappan - 1 +_nkappa _ Empty = 0 + +cleanPart :: Seq Int -> Seq Int +cleanPart kappa = + let i = elemIndexL 0 kappa + in if isJust i + then S.take (fromJust i) kappa + else kappa
tests/Approx.hs view
@@ -1,8 +1,8 @@-module Approx where-import Data.Complex--approx :: Int -> Double -> Double-approx n x = fromInteger (round $ x * (10^n)) / (10.0^^n)--approx' :: Int -> Complex Double -> Complex Double-approx' n z = approx n (realPart z) :+ approx n (imagPart z)+module Approx where +import Data.Complex + +approx :: Int -> Double -> Double +approx n x = fromInteger (round $ x * (10^n)) / (10.0^^n) + +approx' :: Int -> Complex Double -> Complex Double +approx' n z = approx n (realPart z) :+ approx n (imagPart z)
tests/Main.hs view
@@ -1,47 +1,46 @@-module Main where-import Approx-import Data.Complex-import Data.Ratio-import Math.HypergeoMatrix-import Test.Tasty (defaultMain, testGroup)-import Test.Tasty.HUnit (assertEqual, testCase)-main :: IO ()-main = defaultMain $- testGroup "Tests"- [ testCase "a 2F1 value" $ do- let alpha = 2 :: Double- h <- hypergeomat 10 2 [1,2] [3] [0.2, 0.5]- assertEqual ""- (approx 8 1.79412894456143)- (approx 8 h),-- testCase "a complex 2F1 value" $ do- let c = 2 :+ 3 :: Complex Double- h <- hypergeomat 10 2 [1,2] [c] [0.2 :+ 1, 0.5]- assertEqual ""- (approx' 6 (1.887753 :+ 0.566665))- (approx' 6 h),-- testCase "compare with rational" $ do- h1 <- hypergeomat 10 2 [1%2, 3] [3%2, 1%3, 2] [1%5, 1%4, 1%8]- let h1' = fromRational h1- h2 <- hypergeomat 10 (2::Double) [1/2, 3] [3/2, 1/3, 2] [1/5, 1/4, 1/8]- assertEqual ""- (approx 15 h1')- (approx 15 h2),-- testCase "0F0 = exponential of trace" $ do- let x = [0.1, 0.2, 0.1 :+ 0.3] :: [Complex Double]- h <- hypergeomat 20 2 [] [] x- assertEqual ""- (approx' 10 (exp(sum x)))- (approx' 10 h),-- testCase "1F0 is det(I-X)^(-a)" $ do- let x = [0.4, 0.45, 0.5] :: [Double]- a = 2 :: Double- h <- hypergeomat 35 2 [a] [] x- assertEqual ""- (approx 4 (product(map (1 -) x)**(-a)))- (approx 4 h)- ]+module Main where +import Approx +import Data.Complex +import Data.Ratio +import Math.HypergeoMatrix +import Test.Tasty (defaultMain, testGroup) +import Test.Tasty.HUnit (assertEqual, testCase) +main :: IO () +main = defaultMain $ + testGroup "Tests" + [ testCase "a 2F1 value" $ do + h <- hypergeomat 10 2 [1,2] [3] [0.2, 0.5] + assertEqual "" + (approx 8 1.79412894456143) + (approx 8 h), + + testCase "a complex 2F1 value" $ do + let c = 2 :+ 3 :: Complex Double + h <- hypergeomat 10 2 [1,2] [c] [0.2 :+ 1, 0.5] + assertEqual "" + (approx' 6 (1.887753 :+ 0.566665)) + (approx' 6 h), + + testCase "compare with rational" $ do + h1 <- hypergeomat 10 2 [1%2, 3] [3%2, 1%3, 2] [1%5, 1%4, 1%8] + let h1' = fromRational h1 + h2 <- hypergeomat 10 (2::Double) [1/2, 3] [3/2, 1/3, 2] [1/5, 1/4, 1/8] + assertEqual "" + (approx 15 h1') + (approx 15 h2), + + testCase "0F0 = exponential of trace" $ do + let x = [0.1, 0.2, 0.1 :+ 0.3] :: [Complex Double] + h <- hypergeomat 20 2 [] [] x + assertEqual "" + (approx' 10 (exp(sum x))) + (approx' 10 h), + + testCase "1F0 is det(I-X)^(-a)" $ do + let x = [0.4, 0.45, 0.5] :: [Double] + a = 2 :: Double + h <- hypergeomat 35 2 [a] [] x + assertEqual "" + (approx 4 (product(map (1 -) x)**(-a))) + (approx 4 h) + ]