packages feed

hypergeomatrix (empty) → 1.0.0.0

raw patch · 11 files changed

+580/−0 lines, 11 filesdep +arraydep +basedep +containerssetup-changed

Dependencies added: array, base, containers, cyclotomic, hypergeomatrix, tasty, tasty-hunit

Files

+ CHANGELOG.md view
@@ -0,0 +1,3 @@+1.0.0.0+-------+* initial release
+ LICENSE view
@@ -0,0 +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.
+ README.md view
@@ -0,0 +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.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ hypergeomatrix.cabal view
@@ -0,0 +1,49 @@+cabal-version:       2.2+name:                hypergeomatrix+version:             1.0.0.0+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+                     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+                     , 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+                      , tasty-hunit+                      , hypergeomatrix+  Default-Language:     Haskell2010++source-repository head+  type:     git+  location: https://github.com/stla/hypergeomatrix
+ src/Math/HypergeoMatrix.hs view
@@ -0,0 +1,3 @@+module Math.HypergeoMatrix (module X) where+import Math.HypergeoMatrix.Gaussian       as X+import Math.HypergeoMatrix.HypergeoMatrix as X
+ src/Math/HypergeoMatrix/Gaussian.hs view
@@ -0,0 +1,8 @@+module Math.HypergeoMatrix.Gaussian +  where+import Data.Complex.Cyclotomic++type GaussianRational = Cyclotomic++(+:) :: Rational -> Rational -> GaussianRational+(+:) = gaussianRat
+ src/Math/HypergeoMatrix/HypergeoMatrix.hs view
@@ -0,0 +1,158 @@+{-# 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
+ src/Math/HypergeoMatrix/Internal.hs view
@@ -0,0 +1,150 @@+{-# 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
+ tests/Approx.hs view
@@ -0,0 +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)
+ tests/Main.hs view
@@ -0,0 +1,47 @@+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)+  ]