diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -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`
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -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.
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -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.
diff --git a/Setup.hs b/Setup.hs
--- a/Setup.hs
+++ b/Setup.hs
@@ -1,2 +1,2 @@
-import Distribution.Simple
-main = defaultMain
+import Distribution.Simple
+main = defaultMain
diff --git a/hypergeomatrix.cabal b/hypergeomatrix.cabal
--- a/hypergeomatrix.cabal
+++ b/hypergeomatrix.cabal
@@ -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
diff --git a/src/Math/HypergeoMatrix.hs b/src/Math/HypergeoMatrix.hs
--- a/src/Math/HypergeoMatrix.hs
+++ b/src/Math/HypergeoMatrix.hs
@@ -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
diff --git a/src/Math/HypergeoMatrix/Gaussian.hs b/src/Math/HypergeoMatrix/Gaussian.hs
--- a/src/Math/HypergeoMatrix/Gaussian.hs
+++ b/src/Math/HypergeoMatrix/Gaussian.hs
@@ -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
diff --git a/src/Math/HypergeoMatrix/HypergeoMatrix.hs b/src/Math/HypergeoMatrix/HypergeoMatrix.hs
--- a/src/Math/HypergeoMatrix/HypergeoMatrix.hs
+++ b/src/Math/HypergeoMatrix/HypergeoMatrix.hs
@@ -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
diff --git a/src/Math/HypergeoMatrix/Internal.hs b/src/Math/HypergeoMatrix/Internal.hs
--- a/src/Math/HypergeoMatrix/Internal.hs
+++ b/src/Math/HypergeoMatrix/Internal.hs
@@ -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
diff --git a/tests/Approx.hs b/tests/Approx.hs
--- a/tests/Approx.hs
+++ b/tests/Approx.hs
@@ -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)
diff --git a/tests/Main.hs b/tests/Main.hs
--- a/tests/Main.hs
+++ b/tests/Main.hs
@@ -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)
+  ]