diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,3 @@ +1.0.0.0 +------- +* initial release diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/hypergeomatrix.cabal b/hypergeomatrix.cabal new file mode 100644 --- /dev/null +++ b/hypergeomatrix.cabal @@ -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 diff --git a/src/Math/HypergeoMatrix.hs b/src/Math/HypergeoMatrix.hs new file mode 100644 --- /dev/null +++ b/src/Math/HypergeoMatrix.hs @@ -0,0 +1,3 @@ +module Math.HypergeoMatrix (module X) where +import Math.HypergeoMatrix.Gaussian as X +import Math.HypergeoMatrix.HypergeoMatrix as X diff --git a/src/Math/HypergeoMatrix/Gaussian.hs b/src/Math/HypergeoMatrix/Gaussian.hs new file mode 100644 --- /dev/null +++ b/src/Math/HypergeoMatrix/Gaussian.hs @@ -0,0 +1,8 @@ +module Math.HypergeoMatrix.Gaussian + where +import Data.Complex.Cyclotomic + +type GaussianRational = Cyclotomic + +(+:) :: Rational -> Rational -> GaussianRational +(+:) = gaussianRat diff --git a/src/Math/HypergeoMatrix/HypergeoMatrix.hs b/src/Math/HypergeoMatrix/HypergeoMatrix.hs new file mode 100644 --- /dev/null +++ b/src/Math/HypergeoMatrix/HypergeoMatrix.hs @@ -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 diff --git a/src/Math/HypergeoMatrix/Internal.hs b/src/Math/HypergeoMatrix/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Math/HypergeoMatrix/Internal.hs @@ -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 diff --git a/tests/Approx.hs b/tests/Approx.hs new file mode 100644 --- /dev/null +++ b/tests/Approx.hs @@ -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) diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -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) + ]