diff --git a/README.markdown b/README.markdown
--- a/README.markdown
+++ b/README.markdown
@@ -18,18 +18,11 @@
 # Get involved!
 
 Please report bugs via the
-[github issue tracker](https://github.com/bos/statistics/issues).
-
-Master [git mirror](https://github.com/bos/statistics):
-
-* `git clone git://github.com/bos/statistics.git`
-
-There's also a [Mercurial mirror](https://bitbucket.org/bos/statistics):
-
-* `hg clone https://bitbucket.org/bos/statistics`
+[github issue tracker](https://github.com/haskell/statistics/issues).
 
-(You can create and contribute changes using either Mercurial or git.)
+Master [git mirror](https://github.com/haskell/statistics):
 
+* `git clone git://github.com/haskell/statistics.git`
 
 # Authors
 
diff --git a/Statistics/Distribution/NegativeBinomial.hs b/Statistics/Distribution/NegativeBinomial.hs
--- a/Statistics/Distribution/NegativeBinomial.hs
+++ b/Statistics/Distribution/NegativeBinomial.hs
@@ -47,7 +47,7 @@
 gChoose :: Double -> Int -> Double
 gChoose n k
     | k < 0             = 0
-    | k' >= 50          = exp $ logChooseFast n k' 
+    | k' >= 50          = exp $ logChooseFast n k'
     | otherwise         = foldl' (*) 1 factors
     where   factors = [ (n - k' + j) / j | j <- [1..k'] ]
             k' = fromIntegral k
@@ -151,7 +151,7 @@
   | k < 0        = 1
   | otherwise    = incompleteBeta (fromIntegral (k+1)) r (1 - p)
   where
-    k = (floor x)::Integer
+    k = floor x :: Integer
 
 mean :: NegativeBinomialDistribution -> Double
 mean (NBD r p) = r * (1 - p)/p
@@ -166,14 +166,14 @@
   dropWhile (>= -m_epsilon) $
   [ let x = probability d k in x * log x | k <- [0..]]
 
--- | Construct negative binomial distribution. Number of failures /r/
+-- | Construct negative binomial distribution. Number of successes /r/
 --   must be positive and probability must be in (0,1] range
 negativeBinomial :: Double              -- ^ Number of successes.
                  -> Double              -- ^ Success probability.
                  -> NegativeBinomialDistribution
 negativeBinomial r p = maybe (error $ errMsg r p) id $ negativeBinomialE r p
 
--- | Construct negative binomial distribution. Number of failures /r/
+-- | Construct negative binomial distribution. Number of successes /r/
 --   must be positive and probability must be in (0,1] range
 negativeBinomialE :: Double              -- ^ Number of successes.
                   -> Double              -- ^ Success probability.
diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs
--- a/Statistics/Distribution/StudentT.hs
+++ b/Statistics/Distribution/StudentT.hs
@@ -26,7 +26,7 @@
 import Data.Data           (Data, Typeable)
 import GHC.Generics        (Generic)
 import Numeric.SpecFunctions (
-  logBeta, incompleteBeta, invIncompleteBeta, digamma)
+  logBeta, incompleteBeta, invIncompleteBeta, digamma, log1p)
 
 import qualified Statistics.Distribution as D
 import Statistics.Distribution.Transform (LinearTransform (..))
@@ -94,8 +94,9 @@
 
 
 logDensityUnscaled :: StudentT -> Double -> Double
-logDensityUnscaled (StudentT ndf) x =
-    log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)
+logDensityUnscaled (StudentT ndf) x
+  = log1p (x*x/ndf) * (-(0.5 * (1 + ndf)))
+  - logBeta 0.5 (0.5 * ndf)
 
 quantile :: StudentT -> Double -> Double
 quantile (StudentT ndf) p
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE BangPatterns #-}
 -- |
 -- Module    : Statistics.Sample
 -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan
@@ -452,8 +453,9 @@
 
 -- (^) operator from Prelude is just slow.
 (^) :: Double -> Int -> Double
-x ^ 1 = x
-x ^ n = x * (x ^ (n-1))
+x0 ^ n0 = go (n0-1) x0 where
+    go 0 !acc = acc
+    go n  acc = go (n-1) (acc*x0)
 {-# INLINE (^) #-}
 
 -- don't support polymorphism, as we can't get unboxed returns if we use it.
diff --git a/Statistics/Sample/Internal.hs b/Statistics/Sample/Internal.hs
--- a/Statistics/Sample/Internal.hs
+++ b/Statistics/Sample/Internal.hs
@@ -14,9 +14,10 @@
     (
       robustSumVar
     , sum
+    , sumF
     ) where
 
-import Numeric.Sum (kbn, sumVector)
+import qualified Numeric.Sum as Sum
 import Prelude hiding (sum)
 import Statistics.Function (square)
 import qualified Data.Vector.Generic as G
@@ -26,5 +27,9 @@
 {-# INLINE robustSumVar #-}
 
 sum :: (G.Vector v Double) => v Double -> Double
-sum = sumVector kbn
+sum = Sum.sumVector Sum.kbn
 {-# INLINE sum #-}
+
+sumF :: Foldable f => f Double -> Double
+sumF = Sum.sum Sum.kbn
+{-# INLINE sumF #-}
diff --git a/Statistics/Test/Bartlett.hs b/Statistics/Test/Bartlett.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Test/Bartlett.hs
@@ -0,0 +1,99 @@
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-|
+Module      : Statistics.Test.Bartlett
+Description : Bartlett's test for homogeneity of variances.
+Copyright   : (c) Praneya Kumar, Alexey Khudyakov, 2025
+License     : BSD-3-Clause
+
+Bartlett's test is used to check that multiple groups of observations
+come from distributions with equal variances. This test assumes that
+samples come from normal distribution. If this is not the case it may
+simple test for non-normality and Levene's ("Statistics.Test.Levene")
+is preferred
+
+>>> import qualified Data.Vector.Unboxed as VU
+>>> import Statistics.Test.Bartlett
+>>> :{
+let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
+    c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
+in bartlettTest [a,b,c]
+:}
+Right (Test {testSignificance = mkPValue 1.1254782518843598e-5, testStatistics = 22.789434813726768, testDistribution = chiSquared 2})
+
+-}
+module Statistics.Test.Bartlett (
+    bartlettTest,
+    module Statistics.Distribution.ChiSquared
+) where
+
+import qualified Data.Vector           as V
+import qualified Data.Vector.Unboxed   as VU
+import qualified Data.Vector.Generic   as VG
+import qualified Data.Vector.Storable  as VS
+import qualified Data.Vector.Primitive as VP
+#if MIN_VERSION_vector(0,13,2)
+import qualified Data.Vector.Strict    as VV
+#endif
+
+import Statistics.Distribution (complCumulative)
+import Statistics.Distribution.ChiSquared (chiSquared, ChiSquared(..))
+import Statistics.Sample (varianceUnbiased)
+import Statistics.Types (mkPValue)
+import Statistics.Test.Types (Test(..))
+
+-- | Perform Bartlett's test for equal variances. The input is a list
+--   of vectors, where each vector represents a group of observations.
+bartlettTest :: VG.Vector v Double => [v Double] -> Either String (Test ChiSquared)
+bartlettTest groups
+  | length groups < 2                 = Left "At least two groups are required for Bartlett's test."
+  | any ((< 2) . VG.length) groups    = Left "Each group must have at least two observations."
+  | any ((<= 0) . var) groupVariances = Left "All groups must have positive variance."
+  | otherwise = Right Test
+      { testSignificance = pValue
+      , testStatistics   = tStatistic
+      , testDistribution = chiDist
+      }
+  where
+    -- Number of groups
+    k = length groups
+    -- Sample sizes for each group
+    ni  = map (fromIntegral . VG.length) groups
+    -- Total number of observations across all groups
+    n_tot = sum $ fromIntegral . VG.length <$> groups
+    -- Variance estimates
+    groupVariances = toVar <$> groups
+    sumWeightedVars = sum [ (n - 1) * v | Var{sampleN=n, var=v} <- groupVariances ]
+    pooledVariance  = sumWeightedVars / fromIntegral (n_tot - k)
+    -- Numerator of Bartlett's statistic
+    numerator =
+      fromIntegral (n_tot - k) * log pooledVariance -
+      sum [ (n - 1) * log v | Var{sampleN=n, var=v} <- groupVariances ]
+    -- Denominator correction term
+    sumReciprocals = sum [1 / (n - 1) | n <- ni]
+    denomCorrection =
+      1 + (sumReciprocals - 1 / fromIntegral (n_tot - k)) / (3 * (fromIntegral k - 1))
+
+    -- Test statistic and test distrubution
+    tStatistic = max 0 $ numerator / denomCorrection
+    chiDist    = chiSquared (k - 1)
+    pValue     = mkPValue $ complCumulative chiDist tStatistic
+{-# SPECIALIZE bartlettTest :: [V.Vector  Double] -> Either String (Test ChiSquared) #-}
+{-# SPECIALIZE bartlettTest :: [VU.Vector Double] -> Either String (Test ChiSquared) #-}
+{-# SPECIALIZE bartlettTest :: [VS.Vector Double] -> Either String (Test ChiSquared) #-}
+{-# SPECIALIZE bartlettTest :: [VP.Vector Double] -> Either String (Test ChiSquared) #-}
+#if MIN_VERSION_vector(0,13,2)
+{-# SPECIALIZE bartlettTest :: [VV.Vector Double] -> Either String (Test ChiSquared) #-}
+#endif
+
+-- Estimate of variance
+data Var = Var
+  { sampleN :: !Double -- ^ N of elements
+  , var     :: !Double -- ^ Sample variance
+  }
+
+toVar :: VG.Vector v Double => v Double -> Var
+toVar xs = Var { sampleN = fromIntegral $ VG.length xs
+               , var     = varianceUnbiased xs
+               }
diff --git a/Statistics/Test/ChiSquared.hs b/Statistics/Test/ChiSquared.hs
--- a/Statistics/Test/ChiSquared.hs
+++ b/Statistics/Test/ChiSquared.hs
@@ -17,8 +17,8 @@
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
-
-
+import qualified Data.Vector.Fusion.Bundle as F
+import qualified Numeric.Sum as Sum
 
 -- | Generic form of Pearson chi squared tests for binned data. Data
 --   sample is supplied in form of tuples (observed quantity,
@@ -26,7 +26,7 @@
 --
 --   This test should be used only if all bins have expected values of
 --   at least 5.
-chi2test :: (G.Vector v (Int,Double), G.Vector v Double)
+chi2test :: (G.Vector v (Int,Double))
          => Int                 -- ^ Number of additional degrees of
                                 --   freedom. One degree of freedom
                                 --   is due to the fact that the are
@@ -44,7 +44,10 @@
   | otherwise = Nothing
   where
     n     = G.length vec - ndf - 1
-    chi2  = sum $ G.map (\(o,e) -> square (fromIntegral o - e) / e) vec
+    chi2  = Sum.kbn
+          $ F.foldl' Sum.add Sum.zero
+          $ F.map (\(o,e) -> square (fromIntegral o - e) / e)
+          $ G.stream vec
     d     = chiSquared n
 {-# INLINABLE  chi2test #-}
 {-# SPECIALIZE
@@ -56,7 +59,7 @@
 -- | Chi squared test for data with normal errors. Data is supplied in
 --   form of pair (observation with error, and expectation).
 chi2testCont
-  :: (G.Vector v (Estimate NormalErr Double, Double), G.Vector v Double)
+  :: (G.Vector v (Estimate NormalErr Double, Double))
   => Int                                   -- ^ Number of additional
                                            --   degrees of freedom.
   -> v (Estimate NormalErr Double, Double) -- ^ Observation and expectation.
@@ -71,5 +74,8 @@
   | otherwise = Nothing
   where
     n     = G.length vec - ndf - 1
-    chi2  = sum $ G.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s) vec
+    chi2  = Sum.kbn
+          $ F.foldl' Sum.add Sum.zero
+          $ F.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s)
+          $ G.stream vec
     d     = chiSquared n
diff --git a/Statistics/Test/Levene.hs b/Statistics/Test/Levene.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Test/Levene.hs
@@ -0,0 +1,153 @@
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-|
+Module      : Statistics.Test.Levene
+Description : Levene's test for homogeneity of variances.
+Copyright   : (c) Praneya Kumar, Alexey Khudyakov, 2025
+License     : BSD-3-Clause
+
+Levene's test used to check whether samples have equal variance. Null
+hypothesis is all samples are from distributions with same variance
+(homoscedacity). Test is robust to non-normality, and versatile with
+mean or median centering.
+
+>>> import qualified Data.Vector.Unboxed as VU
+>>> import Statistics.Test.Levene
+>>> :{
+let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
+    c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
+in levenesTest Median [a, b, c]
+:}
+Right (Test {testSignificance = mkPValue 2.4315059672496814e-3, testStatistics = 7.584952754501659, testDistribution = fDistributionReal 2.0 27.0})
+-}
+module Statistics.Test.Levene (
+    Center(..),
+    levenesTest
+) where
+
+import Control.Monad
+import qualified Data.Vector           as V
+import qualified Data.Vector.Unboxed   as VU
+import qualified Data.Vector.Generic   as VG
+import qualified Data.Vector.Storable  as VS
+import qualified Data.Vector.Primitive as VP
+#if MIN_VERSION_vector(0,13,2)
+import qualified Data.Vector.Strict    as VV
+#endif
+import Statistics.Distribution (complCumulative)
+import Statistics.Distribution.FDistribution (fDistribution, FDistribution)
+import Statistics.Types      (mkPValue)
+import Statistics.Test.Types (Test(..))
+import Statistics.Function   (gsort)
+import Statistics.Sample     (mean)
+
+import qualified Statistics.Sample.Internal as IS
+import Statistics.Quantile
+
+
+-- | Center calculation method
+data Center
+  = Mean             -- ^ Use arithmetic mean
+  | Median           -- ^ Use median
+  | Trimmed !Double  -- ^ Trimmed mean with given proportion to cut from each end
+  deriving (Eq, Show)
+
+-- | Main Levene's test function with full error handling
+levenesTest
+  :: (VG.Vector v Double)
+  => Center      -- ^ Centering method
+  -> [v Double]  -- ^ Input samples
+  -> Either String (Test FDistribution)
+{-# INLINABLE levenesTest #-}
+levenesTest center samples
+  | length samples < 2 = Left "At least two samples required"
+  -- NOTE: We don't have nice way of computing mean of a list!
+  | otherwise = do
+      let residuals = computeResiduals center <$> samples
+      -- Average of all Z
+      let n_tot = sum $ VG.length . vecZ <$> residuals -- Total number of samples
+      let zbar = IS.sumF [ meanZ z * sampleN z
+                         | z <- residuals]
+               / fromIntegral n_tot
+      -- Numerator: Sum over (ni * (Z[i] - Z)^2)
+      let numerator = IS.sumF [ sampleN z * sqr (meanZ z - zbar)
+                              | z <- residuals]
+      -- Denominator: Sum over Σ((dev_ij - zbari)^2)
+      let denominator = IS.sumF
+            [ IS.sum $ VU.map (sqr . subtract (meanZ z)) (vecZ z)
+            | z <- residuals
+            ]
+      -- Handle division by zero and invalid values
+      when (denominator <= 0 || isNaN denominator || isInfinite denominator)
+        $ Left "Invalid denominator in W-statistic calculation"
+      let wStat = (fromIntegral (n_tot - k) / fromIntegral (k - 1)) * (numerator / denominator)
+          fDist = fDistribution (k - 1) (n_tot - k)
+      Right Test { testStatistics   = wStat
+                 , testSignificance = mkPValue $ complCumulative fDist wStat
+                 , testDistribution = fDist
+                 }
+  where
+    k = length samples -- Number of groups
+{-# SPECIALIZE levenesTest :: Center -> [V.Vector  Double] -> Either String (Test FDistribution) #-}
+{-# SPECIALIZE levenesTest :: Center -> [VU.Vector Double] -> Either String (Test FDistribution) #-}
+{-# SPECIALIZE levenesTest :: Center -> [VS.Vector Double] -> Either String (Test FDistribution) #-}
+{-# SPECIALIZE levenesTest :: Center -> [VP.Vector Double] -> Either String (Test FDistribution) #-}
+#if MIN_VERSION_vector(0,13,2)
+{-# SPECIALIZE levenesTest :: Center -> [VV.Vector Double] -> Either String (Test FDistribution) #-}
+#endif
+
+----------------------------------------------------------------
+-- Implementation
+----------------------------------------------------------------
+
+-- | Trim data from both ends with error handling and performance optimization
+trimboth :: (Ord a, Fractional a, VG.Vector v a)
+         => v a
+         -> Double
+         -> v a
+{-# INLINE trimboth #-}
+trimboth vec p
+  | p < 0 || p >= 0.5 = error "Statistics.Test.Levene: trimming: proportion must be between 0 and 0.5"
+  | VG.null vec       = vec
+  | otherwise         = VG.slice lowerCut (upperCut - lowerCut) sorted
+  where
+    n        = VG.length vec
+    sorted   = gsort vec
+    lowerCut = ceiling $ p * fromIntegral n
+    upperCut = n - lowerCut
+
+data Residuals = Residuals
+  { sampleN :: !Double
+  , meanZ   :: !Double
+  , vecZ    :: !(VU.Vector Double)
+  }
+
+computeResiduals
+  :: VG.Vector v Double
+  => Center
+  -> v Double
+  -> Residuals
+{-# INLINE computeResiduals #-}
+computeResiduals method xs = case method of
+  Mean   ->
+    let c  = mean xs
+        zs = VU.map (\x -> abs (x - c)) $ VU.convert xs
+    in makeR zs
+  Median ->
+    let c  = median medianUnbiased xs
+        zs = VU.map (\x -> abs (x - c)) $ VU.convert xs
+    in makeR zs
+  Trimmed p ->
+    let trimmed = trimboth xs p
+        c       = mean trimmed
+        zs      = VU.map (\x -> abs (x - c)) $ VU.convert trimmed
+    in makeR zs
+  where
+    makeR zs = Residuals { sampleN = fromIntegral $ VU.length zs
+                         , meanZ   = mean zs
+                         , vecZ    = zs
+                         }
+
+sqr :: Double -> Double
+sqr x = x * x
diff --git a/benchmark/Main.hs b/benchmark/Main.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/Main.hs
@@ -0,0 +1,77 @@
+module Main where
+
+import Data.Complex
+import Statistics.Sample
+import Statistics.Transform
+import Statistics.Correlation
+import System.Random.MWC
+import qualified Data.Vector.Unboxed as VU
+import qualified Data.Vector.Unboxed.Mutable as MVU
+
+import Bench
+
+
+-- Test sample
+sample :: VU.Vector Double
+sample = VU.create $ do g <- create
+                        MVU.replicateM 10000 (uniform g)
+
+-- Weighted test sample
+sampleW :: VU.Vector (Double,Double)
+sampleW = VU.zip sample (VU.reverse sample)
+
+-- Complex vector for FFT tests
+sampleC :: VU.Vector (Complex Double)
+sampleC = VU.zipWith (:+) sample (VU.reverse sample)
+
+
+-- Simple benchmark for functions from Statistics.Sample
+main :: IO ()
+main =
+  defaultMain
+  [ bgroup "sample"
+    [ bench "range"            $ nf (\x -> range x)            sample
+      -- Mean
+    , bench "mean"             $ nf (\x -> mean x)             sample
+    , bench "meanWeighted"     $ nf (\x -> meanWeighted x)     sampleW
+    , bench "harmonicMean"     $ nf (\x -> harmonicMean x)     sample
+    , bench "geometricMean"    $ nf (\x -> geometricMean x)    sample
+      -- Variance
+    , bench "variance"         $ nf (\x -> variance x)         sample
+    , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample
+    , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW
+      -- Correlation
+    , bench "pearson"          $ nf pearson     sampleW
+    , bench "covariance"       $ nf covariance  sampleW
+    , bench "correlation"      $ nf correlation sampleW
+    , bench "covariance2"      $ nf (covariance2  sample) sample
+    , bench "correlation2"     $ nf (correlation2 sample) sample
+      -- Other
+    , bench "stdDev"           $ nf (\x -> stdDev x)           sample
+    , bench "skewness"         $ nf (\x -> skewness x)         sample
+    , bench "kurtosis"         $ nf (\x -> kurtosis x)         sample
+      -- Central moments
+    , bench "C.M. 2"           $ nf (\x -> centralMoment 2 x)  sample
+    , bench "C.M. 3"           $ nf (\x -> centralMoment 3 x)  sample
+    , bench "C.M. 4"           $ nf (\x -> centralMoment 4 x)  sample
+    , bench "C.M. 5"           $ nf (\x -> centralMoment 5 x)  sample
+    ]
+  , bgroup "FFT"
+    [ bgroup "fft"
+      [ bench  (show n) $ whnf fft   (VU.take n sampleC) | n <- fftSizes ]
+    , bgroup "ifft"
+      [ bench  (show n) $ whnf ifft  (VU.take n sampleC) | n <- fftSizes ]
+    , bgroup "dct"
+      [ bench  (show n) $ whnf dct   (VU.take n sample)  | n <- fftSizes ]
+    , bgroup "dct_"
+      [ bench  (show n) $ whnf dct_  (VU.take n sampleC) | n <- fftSizes ]
+    , bgroup "idct"
+      [ bench  (show n) $ whnf idct  (VU.take n sample)  | n <- fftSizes ]
+    , bgroup "idct_"
+      [ bench  (show n) $ whnf idct_ (VU.take n sampleC) | n <- fftSizes ]
+    ]
+  ]
+
+
+fftSizes :: [Int]
+fftSizes = [32,128,512,2048]
diff --git a/benchmark/bench.hs b/benchmark/bench.hs
deleted file mode 100644
--- a/benchmark/bench.hs
+++ /dev/null
@@ -1,75 +0,0 @@
-import Data.Complex
-import Statistics.Sample
-import Statistics.Transform
-import Statistics.Correlation
-import System.Random.MWC
-import qualified Data.Vector.Unboxed as VU
-import qualified Data.Vector.Unboxed.Mutable as MVU
-
-import Bench
-
-
--- Test sample
-sample :: VU.Vector Double
-sample = VU.create $ do g <- create
-                        MVU.replicateM 10000 (uniform g)
-
--- Weighted test sample
-sampleW :: VU.Vector (Double,Double)
-sampleW = VU.zip sample (VU.reverse sample)
-
--- Complex vector for FFT tests
-sampleC :: VU.Vector (Complex Double)
-sampleC = VU.zipWith (:+) sample (VU.reverse sample)
-
-
--- Simple benchmark for functions from Statistics.Sample
-main :: IO ()
-main =
-  defaultMain
-  [ bgroup "sample"
-    [ bench "range"            $ nf (\x -> range x)            sample
-      -- Mean
-    , bench "mean"             $ nf (\x -> mean x)             sample
-    , bench "meanWeighted"     $ nf (\x -> meanWeighted x)     sampleW
-    , bench "harmonicMean"     $ nf (\x -> harmonicMean x)     sample
-    , bench "geometricMean"    $ nf (\x -> geometricMean x)    sample
-      -- Variance
-    , bench "variance"         $ nf (\x -> variance x)         sample
-    , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample
-    , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW
-      -- Correlation
-    , bench "pearson"          $ nf pearson     sampleW
-    , bench "covariance"       $ nf covariance  sampleW
-    , bench "correlation"      $ nf correlation sampleW
-    , bench "covariance2"      $ nf (covariance2  sample) sample
-    , bench "correlation2"     $ nf (correlation2 sample) sample
-      -- Other
-    , bench "stdDev"           $ nf (\x -> stdDev x)           sample
-    , bench "skewness"         $ nf (\x -> skewness x)         sample
-    , bench "kurtosis"         $ nf (\x -> kurtosis x)         sample
-      -- Central moments
-    , bench "C.M. 2"           $ nf (\x -> centralMoment 2 x)  sample
-    , bench "C.M. 3"           $ nf (\x -> centralMoment 3 x)  sample
-    , bench "C.M. 4"           $ nf (\x -> centralMoment 4 x)  sample
-    , bench "C.M. 5"           $ nf (\x -> centralMoment 5 x)  sample
-    ]
-  , bgroup "FFT"
-    [ bgroup "fft"
-      [ bench  (show n) $ whnf fft   (VU.take n sampleC) | n <- fftSizes ]
-    , bgroup "ifft"
-      [ bench  (show n) $ whnf ifft  (VU.take n sampleC) | n <- fftSizes ]
-    , bgroup "dct"
-      [ bench  (show n) $ whnf dct   (VU.take n sample)  | n <- fftSizes ]
-    , bgroup "dct_"
-      [ bench  (show n) $ whnf dct_  (VU.take n sampleC) | n <- fftSizes ]
-    , bgroup "idct"
-      [ bench  (show n) $ whnf idct  (VU.take n sample)  | n <- fftSizes ]
-    , bgroup "idct_"
-      [ bench  (show n) $ whnf idct_ (VU.take n sampleC) | n <- fftSizes ]
-    ]
-  ]
-
-
-fftSizes :: [Int]
-fftSizes = [32,128,512,2048]
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,13 @@
+## Changes in 0.16.4.0 [2025.10.23]
+
+ * Bartlett's test (`Statistics.Test.Bartlett`) and Levene's test
+   (`Statistics.Test.Levene`) for homogeneity of variances is added.
+
+ * Improved performance in calculation of moments.
+
+ * Improved precision in calculation of `logDensity` of Student T distribution.
+
+
 ## Changes in 0.16.3.0
 
  * `S.Sample.correlation`, `S.Sample.covariance`,
@@ -61,11 +71,11 @@
 
  * Computation of CDF and quantiles of Cauchy distribution is now numerically
    stable.
- 
+
  * Fix loss of precision in computing of CDF of gamma distribution
 
  * Log-normal and Weibull distributions added.
- 
+
  * `DiscreteGen` instance added for `DiscreteUniform`
 
 
@@ -130,7 +140,7 @@
 ## Changes in 0.14.0.0
 
 Breaking update. It seriously changes parts of API. It adds new data types for
-dealing with with estimates, confidence intervals, confidence levels and
+dealing with estimates, confidence intervals, confidence levels and
 p-value. Also API for statistical tests is changed.
 
  * Module `Statistis.Types` now contains new data types for estimates,
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -2,7 +2,7 @@
 build-type:     Simple
 
 name:           statistics
-version:        0.16.3.0
+version:        0.16.4.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -36,7 +36,6 @@
 
 extra-source-files:
   README.markdown
-  changelog.md
   examples/kde/KDE.hs
   examples/kde/data/faithful.csv
   examples/kde/kde.html
@@ -44,6 +43,9 @@
   tests/utils/Makefile
   tests/utils/fftw.c
 
+extra-doc-files:
+  changelog.md
+
 tested-with:
   GHC ==8.4.4
    || ==8.6.5
@@ -52,9 +54,10 @@
    || ==9.0.2
    || ==9.2.8
    || ==9.4.8
-   || ==9.6.6
+   || ==9.6.7
    || ==9.8.4
-   || ==9.10.1
+   || ==9.10.2
+   || ==9.12.2
 
 source-repository head
   type:     git
@@ -105,6 +108,8 @@
     Statistics.Sample.KernelDensity.Simple
     Statistics.Sample.Normalize
     Statistics.Sample.Powers
+    Statistics.Test.Bartlett
+    Statistics.Test.Levene
     Statistics.Test.ChiSquared
     Statistics.Test.KolmogorovSmirnov
     Statistics.Test.KruskalWallis
@@ -132,7 +137,7 @@
                , binary                  >= 0.5.1.0
                , primitive               >= 0.3
                , dense-linear-algebra    >= 0.1 && <0.2
-               , parallel                >= 3.2.2.0 && <3.3
+               , parallel                >= 3.2.2.0 && <3.4
                , vector                  >= 0.10
                , vector-algorithms       >= 0.4
                , vector-th-unbox
@@ -169,6 +174,8 @@
     Tests.Quantile
   ghc-options:
     -Wall -threaded -rtsopts -fsimpl-tick-factor=500
+  if impl(ghc >= 9.8)
+    ghc-options: -Wno-x-partial
   build-depends: base
                , statistics
                , dense-linear-algebra
@@ -201,7 +208,7 @@
   build-depends:
             base       -any
           , statistics -any
-          , doctest    >=0.15 && <0.24
+          , doctest    >=0.15 && <0.25
 
 -- We want to be able to build benchmarks using both tasty-bench and tasty-papi.
 -- They have similar API so we just create two shim modules which reexport
@@ -219,7 +226,7 @@
   import:         bench-stanza
   type:           exitcode-stdio-1.0
   hs-source-dirs: benchmark bench-time
-  main-is:        bench.hs
+  main-is:        Main.hs
   Other-modules:  Bench
   build-depends:  tasty-bench >= 0.3
 
@@ -229,6 +236,6 @@
   if impl(ghcjs) || !flag(BenchPAPI)
      buildable: False
   hs-source-dirs: benchmark bench-papi
-  main-is:        bench.hs
+  main-is:        Main.hs
   Other-modules:  Bench
   build-depends:  tasty-papi >= 0.1.2
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
--- a/tests/Tests/Distribution.hs
+++ b/tests/Tests/Distribution.hs
@@ -340,7 +340,7 @@
   quantileIsInvCDF_enabled _ = False
 
 instance Param BetaDistribution where
-  -- FIXME: See https://github.com/bos/statistics/issues/161 for details
+  -- FIXME: See https://github.com/haskell/statistics/issues/161 for details
   quantileIsInvCDF_enabled _ = False
 
 instance Param FDistribution where
diff --git a/tests/Tests/Parametric.hs b/tests/Tests/Parametric.hs
--- a/tests/Tests/Parametric.hs
+++ b/tests/Tests/Parametric.hs
@@ -4,12 +4,18 @@
 import Statistics.Test.StudentT
 import Statistics.Types
 import qualified Data.Vector.Unboxed as U
-import Test.Tasty (testGroup)
-import Tests.Helpers  (testEquality)
+import qualified Data.Vector as V
+import Test.Tasty (testGroup, TestTree)
+import Test.Tasty.HUnit (testCase, assertBool)
+import Tests.Helpers (testEquality)
 import qualified Test.Tasty as Tst
 
+import Statistics.Test.Levene
+import Statistics.Test.Bartlett
+
+
 tests :: Tst.TestTree
-tests = testGroup "Parametric tests" studentTTests
+tests = testGroup "Parametric tests" [studentTTests, bartlettTests, leveneTests]
 
 -- 2 samples x 20 obs data
 --
@@ -77,9 +83,9 @@
   , testEquality name (isSignificant (mkPValue $ pValue pVal + 1e-5) test)
     Significant
   ]
-  
-studentTTests :: [Tst.TestTree]
-studentTTests = concat
+
+studentTTests :: Tst.TestTree
+studentTTests = testGroup "StudentT test" $ concat
   [ -- R: t.test(sample1, sample2, alt="two.sided", var.equal=T)
     testTTest "two-sample t-test SamplesDiffer Student"
       (mkPValue 0.03410) (fromJust $ studentTTest SamplesDiffer sample1 sample2)
@@ -100,3 +106,119 @@
       (mkPValue 0.01705) (fromJust $ pairedTTest BGreater sample12)
   ]
   where sample12 = U.zip sample1 sample2
+
+
+------------------------------------------------------------
+-- Bartlett's Test
+------------------------------------------------------------
+
+bartlettTests :: TestTree
+bartlettTests = testGroup "Bartlett's test"
+  [ testCase "a,b,c" $ testBartlettTest [a,b,c] 1.8027132567760222   0.40601846976301237
+  , testCase "a,b"   $ testBartlettTest [a,b]   0.005221063776321886 0.9423974408021293
+  , testCase "a,c"   $ testBartlettTest [a,c]   1.1531619271845452   0.2828882244527482
+  , testCase "a,a"   $ testBartlettTest [a,a]   0.0                  1.0
+  ]
+  where
+    a = U.fromList [9.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = U.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 9.36, 9.18, 8.67, 9.05]
+    c = U.fromList [8.95, 8.12, 8.95, 8.85, 8.03, 8.84, 8.07, 8.98, 8.86, 8.98]
+
+testBartlettTest
+  :: [U.Vector Double]
+  -> Double
+  -> Double
+  -> IO ()
+testBartlettTest samples w p = do
+  r <- case bartlettTest samples of
+    Left  _ -> error "Bartlett's test failed"
+    Right r -> pure r
+  approxEqual "W" 1e-9 (testStatistics r)            w
+  approxEqual "p" 1e-9 (pValue $ testSignificance r) p
+
+------------------------------------------------------------
+-- Levene's Test (Trimmed Mean)
+------------------------------------------------------------
+
+leveneTests :: TestTree
+leveneTests = testGroup "Levene test"
+  -- Statistics' value and p-values are computed using 
+  [ testCase "a,b,c Mean"    $ testLeveneTest [a,b,c] Mean   7.905194483442054 0.001983795817472731
+  , testCase "a,b   Mean"    $ testLeveneTest [a,b]   Mean   8.83873787256358  0.008149720958328811
+  , testCase "a,a   Mean"    $ testLeveneTest [a,a]   Mean   0.0               1.0
+  , testCase "a,b,c Median"  $ testLeveneTest [a,b,c] Median 7.584952754501659 0.002431505967249681
+  , testCase "a,b   Median"  $ testLeveneTest [a,b]   Median 8.461374333228711 0.009364737715584399
+  , testCase "aL,bL Mean"    $ testLeveneTest [aL,bL] Mean   5.84424549939465  0.01653410652558999
+  , testCase "aL,bL Trimmed" $ testLeveneTest [aL,bL] (Trimmed 0.05) 8.368311226366314 0.004294953946529551
+  ]
+  where
+    a = V.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = V.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
+    c = V.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
+    -- Large samples for testing trimmed
+    aL = V.fromList [
+      -0.18919252, -1.62837673,  5.21332355, -0.00962043, -0.28417847,
+      -0.88128233,  1.49698436,  6.1780359 , -1.22301348,  3.34598245,
+       5.33227264, -0.88732069,  0.14487346,  2.61060215,  4.22033907,
+       2.53139215, -0.72131061,  0.53063607, -0.60510374, -0.73230842,
+       1.54037043, -2.81103963,  3.40763063,  0.49005324,  2.13085513,
+       5.68650547,  4.16397279, -0.17325097,  1.12664972,  4.23297516,
+       4.15943436, -1.01452078,  2.40391646,  0.83019962,  0.29665879,
+      -3.83031046, -1.98576933,  1.5356527 ,  1.30773365,  0.292818  ,
+       2.45877828,  1.06482289, -0.63241873,  1.58465379,  1.96577614,
+       2.25791943,  4.13769848, -2.38595767, -0.65801423, -2.54007791,
+       3.17428087,  4.32096964,  0.92240335, -2.38101319,  1.35692587,
+       1.48279101, -0.04438309,  0.50296642,  2.08261495,  1.33181215,
+      -1.95427198,  4.95406809,  1.51294898, -2.68536129, -0.2441218 ,
+       2.41142613,  4.71051493,  2.66618697,  1.12668301, -0.25732583,
+       1.25021838, -1.27523641,  5.01638744,  3.38864442,  0.17979744,
+      -0.88481645,  3.89346357, -0.51512217, -1.60542888,  0.88378679,
+      -2.12962732, -1.35989539,  5.09215112, -1.37442481,  0.83578405,
+       0.13829571,  1.25171481,  3.60552158, -3.24051591, -0.44301834,
+       0.78253445,  1.76098254,  1.79677434, -0.19010505,  3.07640466,
+       3.02853882,  1.24849063,  4.84505382,  6.82274999,  2.24063474]
+    bL = V.fromList [
+        2.15584101, -2.74876744, -0.82231894,  1.97518087,  2.59280595,
+        1.28703417,  2.40450278,  1.9761031 ,  2.35186598,  1.15611047,
+        2.26709318,  1.2832138 , -2.1486074 ,  0.27563011, -0.51816861,
+        0.89658424,  3.27069545,  1.72846646,  3.84454277,  5.58301459,
+       -0.40878188,  3.41602853,  1.1281526 ,  0.9665913 ,  0.76567084,
+        1.69522855,  1.69133014,  0.70529264,  2.65243202, -1.0088019 ,
+       -0.62431026,  3.76667396,  3.66225181,  0.73217579,  0.04478736,
+        0.4169833 ,  0.77065631, -1.31484093,  1.23858618, -0.08339456,
+        3.14154286,  1.84358218, -0.53511423, -3.4919477 ,  0.24076997,
+        3.59381684,  1.99497806,  2.95499775,  1.67157731,  0.0214764 ,
+        3.32161612, -2.64762427,  0.06486472,  0.19653897,  1.34954235,
+        1.18568747, -0.54434597, -3.35544223,  1.41933109,  0.95100195,
+        2.7182116 ,  1.1334068 , -0.95297806, -0.05421818,  1.42248799,
+       -3.96201277, -3.21309254, -0.21209211,  0.9689551 ,  0.13526401,
+       -0.88656198,  0.41331783, -3.18766064,  4.34948246,  1.35656384,
+        0.41920101, -0.46578994,  1.55181583,  2.43937014,  2.49040644,
+        4.10505494,  1.68856296,  1.31503895,  0.41123368,  0.73242999,
+        0.2804349 , -1.83494592, -0.31073195,  2.61185513,  2.91645094,
+        1.26097638,  2.64197134,  3.88931972,  0.03783002,  2.55209729,
+        3.46869549,  0.96348003,  2.27658242,  2.7613171 , -0.1372434 ]
+
+    
+testLeveneTest
+  :: [V.Vector Double]
+  -> Center
+  -> Double
+  -> Double
+  -> IO ()
+testLeveneTest samples center w p = do
+  r <- case levenesTest center samples of
+    Left  _ -> error "Levene's test failed"
+    Right r -> pure r
+  approxEqual "W" 1e-9 (testStatistics r)            w
+  approxEqual "p" 1e-9 (pValue $ testSignificance r) p
+
+
+----------------------------------------------------------------
+
+approxEqual :: String -> Double -> Double -> Double -> IO ()
+approxEqual name epsilon actual expected =
+  assertBool (name ++ ": expected ≈ " ++ show expected ++ ", got " ++ show actual)
+             (diff < epsilon)
+  where
+    diff = abs (actual - expected)