diff --git a/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs b/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs
index 185df392..2b773c06 100644
--- a/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs
+++ b/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs
@@ -217,7 +217,9 @@ varScatter takeSqrt g nGroups v = runST $ do
             | otherwise = do
                 c <- VUM.unsafeRead cnt k
                 mm <- VUM.unsafeRead m2 k
-                let var = if c < 2 then 0 else mm / fromIntegral (c - 1)
+                -- Sample variance is undefined at n = 1: NaN, matching
+                -- 'computeVariance'.
+                let var = if c < 2 then 0 / 0 else mm / fromIntegral (c - 1)
                 VUM.unsafeWrite out k (if takeSqrt then sqrt var else var)
                 fin (k + 1)
     fin 0
diff --git a/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs b/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs
index db8ae7a7..05308f7b 100644
--- a/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs
+++ b/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs
@@ -290,7 +290,9 @@ varPar takeSqrt vis offs nGroups v caps bounds = do
             | otherwise = do
                 c <- VUM.unsafeRead cnt k
                 mm <- VUM.unsafeRead m2 k
-                let var = if c < 2 then 0 else mm / fromIntegral (c - 1)
+                -- Sample variance is undefined at n = 1: NaN, matching
+                -- 'computeVariance'.
+                let var = if c < 2 then 0 / 0 else mm / fromIntegral (c - 1)
                 VUM.unsafeWrite out k (if takeSqrt then sqrt var else var)
                 fin (k + 1)
     fin 0
diff --git a/dataframe-learn/src/DataFrame/Metrics.hs b/dataframe-learn/src/DataFrame/Metrics.hs
index 3982efba..e7be7844 100644
--- a/dataframe-learn/src/DataFrame/Metrics.hs
+++ b/dataframe-learn/src/DataFrame/Metrics.hs
@@ -73,15 +73,22 @@ columnOf df e = case interpret @Double df e of
     Right (TColumn c) -> fromRight VU.empty (toVector @Double @VU.Vector c)
     Left err -> throw err
 
-n2 :: VU.Vector Double -> Double
-n2 = fromIntegral . VU.length
+{- | Compared pairs: 'VU.zipWith' truncates to the shorter vector, so every
+mean below divides by this, never by the length of 'truth' alone.
+-}
+nCompared :: VU.Vector Double -> VU.Vector Double -> Double
+nCompared preds truth = fromIntegral (min (VU.length preds) (VU.length truth))
 
 -- | Mean squared error.
 mse :: Metric
 mse preds truth
     | VU.null truth = 0
+    -- No predictions is not a perfect score.
+    | n == 0 = 0 / 0
     | otherwise =
-        VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n2 truth
+        VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n
+  where
+    n = nCompared preds truth
 
 -- | Root mean squared error.
 rmse :: Metric
@@ -91,29 +98,40 @@ rmse preds truth = sqrt (mse preds truth)
 mae :: Metric
 mae preds truth
     | VU.null truth = 0
-    | otherwise = VU.sum (VU.zipWith (\p t -> abs (p - t)) preds truth) / n2 truth
+    | n == 0 = 0 / 0
+    | otherwise = VU.sum (VU.zipWith (\p t -> abs (p - t)) preds truth) / n
+  where
+    n = nCompared preds truth
 
 -- | Coefficient of determination @R²@.
 r2 :: Metric
 r2 preds truth
-    | VU.null truth || ssTot == 0 = 0
+    | VU.null truth = 0
+    | n == 0 = 0 / 0
+    | ssTot == 0 = 0
     | otherwise = 1 - ssRes / ssTot
   where
-    mean = VU.sum truth / n2 truth
-    ssRes = VU.sum (VU.zipWith (\p t -> (t - p) ^ (2 :: Int)) preds truth)
-    ssTot = VU.sum (VU.map (\t -> (t - mean) ^ (2 :: Int)) truth)
+    n = nCompared preds truth
+    truth' = VU.take (min (VU.length preds) (VU.length truth)) truth
+    mean = VU.sum truth' / n
+    ssRes = VU.sum (VU.zipWith (\p t -> (t - p) ^ (2 :: Int)) preds truth')
+    ssTot = VU.sum (VU.map (\t -> (t - mean) ^ (2 :: Int)) truth')
 
 -- | Fraction of exact matches.
 accuracy :: Metric
 accuracy preds truth
     | VU.null truth = 0
+    | n == 0 = 0 / 0
     | otherwise =
-        fromIntegral (VU.length (VU.filter id (VU.zipWith (==) preds truth))) / n2 truth
+        fromIntegral (VU.length (VU.filter id (VU.zipWith (==) preds truth))) / n
+  where
+    n = nCompared preds truth
 
 -- | Binary log loss; probabilities clamped away from @0@/@1@.
 logLoss :: Metric
 logLoss probs truth
     | VU.null truth = 0
+    | n == 0 = 0 / 0
     | otherwise =
         negate
             ( VU.sum
@@ -123,8 +141,9 @@ logLoss probs truth
                     truth
                 )
             )
-            / n2 truth
+            / n
   where
+    n = nCompared probs truth
     clampP p = max 1e-15 (min (1 - 1e-15) p)
 
 -- | Averaging strategy for multiclass precision/recall/F1.
diff --git a/dataframe-learn/tests-internal/Learn/EdgeCases.hs b/dataframe-learn/tests-internal/Learn/EdgeCases.hs
index 6a991570..a3012c0e 100644
--- a/dataframe-learn/tests-internal/Learn/EdgeCases.hs
+++ b/dataframe-learn/tests-internal/Learn/EdgeCases.hs
@@ -28,7 +28,7 @@ import DataFrame.LinearModel
 import DataFrame.LinearSolver (sigmoid)
 import DataFrame.PCA
 
-import DataFrame.Internal.Statistics (correlation', variance')
+import DataFrame.Internal.Statistics (correlation', meanSquaredError, variance')
 
 import Test.HUnit
 
@@ -124,14 +124,13 @@ testVarianceConstant = TestCase $ do
     let v = variance' (VU.replicate 100 (7.0 :: Double))
     assertEqual "variance of constant column is 0" 0 v
 
-{- Variance of fewer than two samples is defined to be 0 (computeVariance guard),
-   not NaN from a /0. -}
+{- Sample variance of one observation is undefined: NaN, so a singleton
+   group can never look as tight as a genuinely constant column. -}
 testVarianceSingleton :: Test
 testVarianceSingleton = TestCase $ do
-    assertEqual
-        "variance of one sample is 0"
-        0
-        (variance' (VU.fromList [3.5 :: Double]))
+    assertBool
+        "variance of one sample is NaN"
+        (isNaN (variance' (VU.fromList [3.5 :: Double])))
 
 {- Correlation of a perfectly linear pair is exactly +1 (and -1 reversed),
    computed stably. y = 2x+1 over a spread of x. -}
@@ -169,6 +168,27 @@ testCorrelationTooFew = TestCase $ do
         Nothing
         (correlation' (VU.fromList [1]) (VU.fromList [2]))
 
+{- meanSquaredError refuses length mismatches and empty inputs rather than
+   averaging over terms it never summed (or indexing out of bounds). -}
+testMeanSquaredErrorGuards :: Test
+testMeanSquaredErrorGuards = TestCase $ do
+    assertEqual
+        "mse of mismatched lengths is Nothing"
+        Nothing
+        (meanSquaredError (VU.fromList [0, 0, 0, 0]) (VU.fromList [2, 2]))
+    assertEqual
+        "mse with the longer prediction does not index out of bounds"
+        Nothing
+        (meanSquaredError (VU.fromList [1]) (VU.fromList [1, 2, 3]))
+    assertEqual
+        "mse of empty inputs is Nothing"
+        Nothing
+        (meanSquaredError VU.empty VU.empty)
+    assertEqual
+        "mse of equal-length inputs is the plain mean"
+        (Just 4.0)
+        (meanSquaredError (VU.fromList [0, 0]) (VU.fromList [2, 2]))
+
 -- ===========================================================================
 -- Category 8: stability inside the model expr layer
 -- ===========================================================================
@@ -425,6 +445,7 @@ tests =
     , testCorrelationPerfect
     , testCorrelationConstantColumnIsNaN
     , testCorrelationTooFew
+    , testMeanSquaredErrorGuards
     , testLogisticProbsExtremeFeatures
     , testOLSOneRow
     , testLogisticSingleClass
diff --git a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
index 0c80ad4a..190e4ec2 100644
--- a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
+++ b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
@@ -17,7 +17,9 @@ import DataFrame.Errors (DataFrameException (..))
 mean' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
 mean' samp
     | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = rtf (VU.sum samp) / fromIntegral (VU.length samp)
+    -- Widen per element: summing at 'a' wraps for Int columns.
+    | otherwise =
+        VU.foldl' (\acc x -> acc + rtf x) 0 samp / fromIntegral (VU.length samp)
 {-# INLINE [0] mean' #-}
 
 meanDouble' :: VU.Vector Double -> Double
@@ -29,7 +31,9 @@ meanDouble' samp
 meanInt' :: VU.Vector Int -> Double
 meanInt' samp
     | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise = fromIntegral (VU.sum samp) / fromIntegral (VU.length samp)
+    | otherwise =
+        VU.foldl' (\acc x -> acc + fromIntegral x) 0 samp
+            / fromIntegral (VU.length samp)
 {-# INLINE meanInt' #-}
 
 {-# RULES
@@ -54,7 +58,8 @@ median' samp
             then pure (rtf middleElement)
             else do
                 prev <- VUM.read mutableSamp (middleIndex - 1)
-                pure (rtf (middleElement + prev) / 2)
+                -- Widen before adding: 'a' addition wraps for Int.
+                pure ((rtf middleElement + rtf prev) / 2)
 {-# INLINE median' #-}
 
 -- accumulator: count, mean, m2
@@ -73,7 +78,9 @@ varianceStep (VarAcc !n !meanVal !m2) !x =
 
 computeVariance :: VarAcc -> Double
 computeVariance (VarAcc !n _ !m2)
-    | n < 2 = 0 -- or error "variance of <2 samples"
+    | n == 0 = throw $ EmptyDataSetException "variance"
+    -- Sample variance is undefined at n = 1: NaN, not a spurious 0.
+    | n < 2 = 0 / 0
     | otherwise = m2 / fromIntegral (n - 1)
 {-# INLINE computeVariance #-}
 
@@ -106,38 +113,37 @@ skewnessStep (SkewAcc !n !meanVal !m2 !m3) !x' =
 computeSkewness :: SkewAcc -> Double
 computeSkewness (SkewAcc n _ m2 m3)
     | n < 3 = 0 -- or error "skewness of <3 samples"
-    | otherwise = (sqrt (fromIntegral n - 1) * m3) / sqrt (m2 ^ (3 :: Int))
+    -- m2, m3 are raw sums, so population g1 = sqrt n * m3 / m2^(3/2).
+    | otherwise = (sqrt (fromIntegral n) * m3) / sqrt (m2 ^ (3 :: Int))
 {-# INLINE computeSkewness #-}
 
 skewness' :: (VU.Unbox a, Real a, Num a) => VU.Vector a -> Double
 skewness' = computeSkewness . VU.foldl' skewnessStep (SkewAcc 0 0 0 0)
 {-# INLINE skewness' #-}
 
-data CorrelationStats
-    = CorrelationStats
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-        {-# UNPACK #-} !Double
-
+{- | Centered two-pass form: the one-pass @n*Sxy - Sx*Sy@ form cancels
+catastrophically on low-variance columns and can report |r| > 1.
+-}
 correlation' :: VU.Vector Double -> VU.Vector Double -> Maybe Double
 correlation' xs ys
     | n < 2 = Nothing
     | VU.length xs /= VU.length ys = Nothing
     | otherwise =
         let nf = fromIntegral n
-            initial = CorrelationStats 0 0 0 0 0
-            (CorrelationStats sumX sumY sumXX sumYY sumXY) = VU.ifoldl' step initial xs
-
-            !num = nf * sumXY - sumX * sumY
-            !den = sqrt ((nf * sumXX - sumX * sumX) * (nf * sumYY - sumY * sumY))
-         in Just (num / den)
+            !mx = VU.sum xs / nf
+            !my = VU.sum ys / nf
+            !sxy = VU.sum (VU.zipWith (\x y -> (x - mx) * (y - my)) xs ys)
+            !sxx = VU.sum (VU.map (\x -> (x - mx) * (x - mx)) xs)
+            !syy = VU.sum (VU.map (\y -> (y - my) * (y - my)) ys)
+         in Just (clamp (sxy / sqrt (sxx * syy)))
   where
     n = VU.length xs
-    step (CorrelationStats sx sy sxx syy sxy) i x =
-        let !y = VU.unsafeIndex ys i
-         in CorrelationStats (sx + x) (sy + y) (sxx + x * x) (syy + y * y) (sxy + x * y)
+    -- Cauchy-Schwarz: clamp the last rounding, never the formula. Explicit
+    -- guards so the zero-variance NaN passes through ('max' would eat it).
+    clamp r
+        | r > 1 = 1
+        | r < -1 = -1
+        | otherwise = r
 {-# INLINE correlation' #-}
 
 quantiles' ::
@@ -202,11 +208,14 @@ interQuartileRange' samp =
 {-# INLINE interQuartileRange' #-}
 
 meanSquaredError :: VU.Vector Double -> VU.Vector Double -> Maybe Double
-meanSquaredError target prediction =
-    let
-        squareDiff = VU.ifoldl' (\sq i e -> (e - target VU.! i) ^ (2 :: Int) + sq) 0 prediction
-     in
-        Just $ squareDiff / fromIntegral (max (VU.length target) (VU.length prediction))
+meanSquaredError target prediction
+    | VU.length target /= VU.length prediction = Nothing
+    | VU.null target = Nothing
+    | otherwise =
+        Just
+            ( VU.sum (VU.zipWith (\t p -> (p - t) ^ (2 :: Int)) target prediction)
+                / fromIntegral (VU.length target)
+            )
 {-# INLINE meanSquaredError #-}
 
 mutualInformationBinned ::
diff --git a/dataframe-operations/src/DataFrame/Operations/Inference.hs b/dataframe-operations/src/DataFrame/Operations/Inference.hs
index 6464d4a2..aec40f5f 100644
--- a/dataframe-operations/src/DataFrame/Operations/Inference.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Inference.hs
@@ -65,15 +65,9 @@ byteStringDateParser "%Y-%m-%d" = parseDateField
 byteStringDateParser fmt = readByteStringDate fmt
 {-# INLINE byteStringDateParser #-}
 
-{- | 'DataFrame.Internal.Parsing.readInt' that rejects overflow instead of
-wrapping. Fields of <= 18 chars cannot overflow and keep the Text-level
-parse; longer (rare) fields take the exact byte-level parser, so an
-overflowing cell demotes\/promotes instead of silently wrapping.
--}
+-- | Alias for 'readInt', which now rejects overflow itself.
 readIntStrict :: T.Text -> Maybe Int
-readIntStrict t
-    | T.length t <= 18 = readInt t
-    | otherwise = parseIntField (TE.encodeUtf8 t)
+readIntStrict = readInt
 {-# INLINE readIntStrict #-}
 
 {- | Candidate-mask priority, reproducing the documented fallback order:
diff --git a/dataframe-operations/src/DataFrame/Operations/Statistics.hs b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
index 7322c0dc..10b6cbeb 100644
--- a/dataframe-operations/src/DataFrame/Operations/Statistics.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
@@ -416,7 +416,10 @@ summarize df =
 
 -- | Round a @Double@ to Specified Precision
 roundTo :: Int -> Double -> Double
-roundTo n x = fromInteger (round $ x * 10 ^ n) / 10.0 ^^ n
+roundTo n x
+    -- 'round' on NaN yields garbage; keep NaN visible in summaries.
+    | isNaN x = x
+    | otherwise = fromInteger (round $ x * 10 ^ n) / 10.0 ^^ n
 
 toPct2dp :: Double -> String
 toPct2dp x
diff --git a/dataframe-operations/src/DataFrame/Operations/Typing.hs b/dataframe-operations/src/DataFrame/Operations/Typing.hs
index 4cc90f96..e314f296 100644
--- a/dataframe-operations/src/DataFrame/Operations/Typing.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Typing.hs
@@ -217,12 +217,12 @@ handleBoolAssumption isNull cols =
         (handleTextAssumption isNull cols)
 
 {- | Int columns: one fused pass with in-place Int -> Double promotion; a cell
-parsing as neither demotes the column to Text. 'readIntStrict' rejects overflow
+parsing as neither demotes the column to Text. 'readInt' rejects overflow
 so a huge integer promotes to 'Double' rather than wrapping.
 -}
 handleIntAssumption :: (T.Text -> Bool) -> V.Vector T.Text -> Column
 handleIntAssumption isNull cols =
-    case promoteIntColumn (\_ t -> isNull t) readIntStrict readDouble cols of
+    case promoteIntColumn (\_ t -> isNull t) readInt readDouble cols of
         Just col -> col
         Nothing -> handleTextAssumption isNull cols
 
diff --git a/dataframe-parsing/dataframe-parsing.cabal b/dataframe-parsing/dataframe-parsing.cabal
index 6e8928a8..476cc287 100644
--- a/dataframe-parsing/dataframe-parsing.cabal
+++ b/dataframe-parsing/dataframe-parsing.cabal
@@ -43,7 +43,6 @@ library
     build-depends:      base >= 4 && < 5,
                         attoparsec >= 0.12 && < 0.16,
                         bytestring >= 0.11 && < 0.14,
-                        bytestring-lexing >= 0.5 && < 0.7,
                         containers >= 0.6.7 && < 0.10,
                         dataframe-core >= 2.4 && < 2.5,
                         text >= 2.1 && < 3,
diff --git a/dataframe-parsing/src-internal/DataFrame/Internal/Parsing.hs b/dataframe-parsing/src-internal/DataFrame/Internal/Parsing.hs
index a7126247..b2d50558 100644
--- a/dataframe-parsing/src-internal/DataFrame/Internal/Parsing.hs
+++ b/dataframe-parsing/src-internal/DataFrame/Internal/Parsing.hs
@@ -8,13 +8,15 @@ module DataFrame.Internal.Parsing where
 import qualified Data.ByteString.Char8 as C
 import qualified Data.Set as S
 import qualified Data.Text as T
+import qualified Data.Text.Encoding as TE
 import qualified Data.Text.IO as TIO
 
 import Control.Applicative (many, (<|>))
+import Control.Monad (guard)
 import Data.Attoparsec.Text hiding (decimal, double, signed)
-import Data.ByteString.Lex.Fractional
+import Data.Char (isDigit)
 import Data.Foldable (fold)
-import Data.Text.Read (decimal, double, signed)
+import Data.Text.Read (decimal, signed)
 import Data.Time (Day, defaultTimeLocale, parseTimeM)
 import GHC.Stack (HasCallStack)
 import System.IO (Handle, IOMode (..), hIsEOF, hTell, withFile)
@@ -61,11 +63,23 @@ readInteger s = case signed decimal (T.strip s) of
     Right (value, "") -> Just value
     Right (_value, _) -> Nothing
 
+{- | 'Data.Text.Read.decimal' at 'Int' wraps on overflow. Fields of <= 18
+chars cannot overflow and keep the fast path; longer (rare) ones go
+through 'Integer' with a range check, mirroring 'readByteStringInt'.
+-}
 readInt :: (HasCallStack) => T.Text -> Maybe Int
-readInt s = case signed decimal (T.strip s) of
-    Left _ -> Nothing
-    Right (value, "") -> Just value
-    Right (_value, _) -> Nothing
+readInt s
+    | T.length t <= 18 = case signed decimal t of
+        Right (value, "") -> Just value
+        _ -> Nothing
+    | otherwise = case signed decimal t :: Either String (Integer, T.Text) of
+        Right (value, "")
+            | value >= toInteger (minBound :: Int)
+            , value <= toInteger (maxBound :: Int) ->
+                Just (fromInteger value)
+        _ -> Nothing
+  where
+    t = T.strip s
 {-# INLINE readInt #-}
 
 readByteStringInt :: (HasCallStack) => C.ByteString -> Maybe Int
@@ -93,23 +107,68 @@ readByteStringInt s
 #endif
 {-# INLINE readByteStringInt #-}
 
+{- | Exact decimal -> 'Double': one correctly-rounded conversion, matching
+strtod\/'read' on every input, subnormals and overflow included. Also
+takes back the Infinity\/-Infinity tokens 'show' emits, so written CSV
+round-trips.
+-}
 readByteStringDouble :: (HasCallStack) => C.ByteString -> Maybe Double
-readByteStringDouble s =
-    let
-        readFunc = if C.any (\c -> c == 'e' || c == 'E') s then readExponential else readDecimal
-     in
-        case readSigned readFunc (C.strip s) of
-            Nothing -> Nothing
-            Just (value, "") -> Just value
-            Just (_value, _) -> Nothing
+readByteStringDouble s
+    | t == "Infinity" = Just (1 / 0)
+    | t == "-Infinity" = Just (-1 / 0)
+    | otherwise = parseDoubleExact t
+  where
+    t = C.strip s
 {-# INLINE readByteStringDouble #-}
 
+parseDoubleExact :: C.ByteString -> Maybe Double
+parseDoubleExact t0 = do
+    let (neg, t1) = case C.uncons t0 of
+            Just ('-', r) -> (True, r)
+            Just ('+', r) -> (False, r)
+            _ -> (False, t0)
+        (ws, t2) = C.span isDigit t1
+    guard (not (C.null ws))
+    (fs, t3) <- case C.uncons t2 of
+        Just ('.', r) ->
+            let (ds, rest) = C.span isDigit r
+             in if C.null ds then Nothing else Just (ds, rest)
+        _ -> Just (C.empty, t2)
+    e <- case C.uncons t3 of
+        Nothing -> Just 0
+        Just (c, r)
+            | c == 'e' || c == 'E' -> do
+                let (eneg, r1) = case C.uncons r of
+                        Just ('-', r') -> (True, r')
+                        Just ('+', r') -> (False, r')
+                        _ -> (False, r)
+                guard (maybe False (isDigit . fst) (C.uncons r1))
+                (ev, rest) <- C.readInteger r1
+                guard (C.null rest)
+                Just (if eneg then negate ev else ev)
+            | otherwise -> Nothing
+    let sigDigits = ws <> fs
+    (sig, _) <- C.readInteger sigDigits
+    let dexp = e - toInteger (C.length fs)
+        nd = toInteger (C.length (C.dropWhile (== '0') sigDigits))
+    Just (mkDouble neg sig nd dexp)
+
+{- | Nearest 'Double' for @sig * 10^dexp@ in ONE rounding step. The clamps
+bound the 'Rational' so an absurd exponent cannot allocate its 10^e.
+-}
+mkDouble :: Bool -> Integer -> Integer -> Integer -> Double
+mkDouble neg sig nd dexp
+    | sig == 0 = sgn 0
+    | dexp > 350 = sgn (1 / 0)
+    -- Even nd digits scaled this low sit under half the smallest denormal.
+    | nd + dexp < -350 = sgn 0
+    | otherwise =
+        sgn (fromRational (fromInteger sig * 10 ^^ (fromInteger dexp :: Int)))
+  where
+    sgn = if neg then negate else id
+
 readDouble :: (HasCallStack) => T.Text -> Maybe Double
-readDouble s =
-    case signed double s of
-        Left _ -> Nothing
-        Right (value, "") -> Just value
-        Right (_value, _) -> Nothing
+readDouble = readByteStringDouble . TE.encodeUtf8
 {-# INLINE readDouble #-}
 
 readIntegerEither :: (HasCallStack) => T.Text -> Either T.Text Integer
@@ -119,19 +178,24 @@ readIntegerEither s = case signed decimal (T.strip s) of
     Right (_value, _) -> Left s
 {-# INLINE readIntegerEither #-}
 
+-- | As 'readInt': overflow is a failed parse, not a wrapped value.
 readIntEither :: (HasCallStack) => T.Text -> Either T.Text Int
-readIntEither s = case signed decimal (T.strip s) of
-    Left _ -> Left s
-    Right (value, "") -> Right value
-    Right (_value, _) -> Left s
+readIntEither s
+    | T.length t <= 18 = case signed decimal t of
+        Right (value, "") -> Right value
+        _ -> Left s
+    | otherwise = case signed decimal t :: Either String (Integer, T.Text) of
+        Right (value, "")
+            | value >= toInteger (minBound :: Int)
+            , value <= toInteger (maxBound :: Int) ->
+                Right (fromInteger value)
+        _ -> Left s
+  where
+    t = T.strip s
 {-# INLINE readIntEither #-}
 
 readDoubleEither :: (HasCallStack) => T.Text -> Either T.Text Double
-readDoubleEither s =
-    case signed double s of
-        Left _ -> Left s
-        Right (value, "") -> Right value
-        Right (_value, _) -> Left s
+readDoubleEither s = maybe (Left s) Right (readDouble s)
 {-# INLINE readDoubleEither #-}
 
 -- ---------------------------------------------------------------------------
diff --git a/dataframe-parsing/src-internal/DataFrame/Internal/Parsing/Fast/Double.hs b/dataframe-parsing/src-internal/DataFrame/Internal/Parsing/Fast/Double.hs
index d5dce287..1fe3cde9 100644
--- a/dataframe-parsing/src-internal/DataFrame/Internal/Parsing/Fast/Double.hs
+++ b/dataframe-parsing/src-internal/DataFrame/Internal/Parsing/Fast/Double.hs
@@ -3,8 +3,9 @@
 {-# LANGUAGE UnboxedTuples #-}
 
 {- | Fast @Double@ slice parser, bit-exact with @readByteStringDouble@.
-Replays the reference parser's exact floating-point operations via 'Word64'
-digit accumulation and 10^k tables, falling back when exactness is in doubt.
+Fields whose significand and scale fit Clinger's exact window are one
+correctly-rounded machine op; everything else falls back to the exact
+reference parser.
 -}
 module DataFrame.Internal.Parsing.Fast.Double (parseDoubleField#) where
 
@@ -18,33 +19,31 @@ import GHC.Exts (Double (..), Double#, Int#)
 import DataFrame.Internal.Parsing (readByteStringDouble)
 import DataFrame.Internal.Parsing.Fast.Common
 
-{- | @10 ^ k@ for @k <= tableMax@, computed with the same @(^)@ the reference
-parser uses, so every entry is bit-identical. Entries from @10^309@ up are
-@Infinity@, so clamping larger exponents to 'tableMax' is exact.
+{- | @10 ^ k@ for @k <= 22@: every entry is exactly representable (5^22
+fits in 53 bits), so multiplying or dividing by one is a single
+correctly-rounded op.
 -}
 pow10Table :: VU.Vector Double
-pow10Table = VU.generate (tableMax + 1) (10 ^)
+pow10Table = VU.generate (exactPow10 + 1) (10 ^)
 {-# NOINLINE pow10Table #-}
 
--- | @recip (10 ^ k)@, replaying @10 ^^ negate k@ bit-for-bit.
-recipPow10Table :: VU.Vector Double
-recipPow10Table = VU.map recip pow10Table
-{-# NOINLINE recipPow10Table #-}
+-- | @10 ^ k@ for @k <= 19@ in 'Word64'; @10^19 < 2^64@.
+pow10w :: VU.Vector Word64
+pow10w = VU.generate 20 (10 ^)
+{-# NOINLINE pow10w #-}
 
-tableMax :: Int
-tableMax = 1024
+-- | Largest @k@ with @10^k@ exactly representable as a 'Double'.
+exactPow10 :: Int
+exactPow10 = 22
 
-{- | 'Word64' to 'Double' exactly as the reference parser's
-'fromInteger' rounds it; values up to @2^53@ take the exact
-'Int' conversion (int2Double#), larger ones the 'Integer' route.
+{- | 'Word64' to 'Double'; callers only pass values @<= 2^53@, where the
+conversion is exact.
 -}
 w2d :: Word64 -> Double
-w2d w
-    | w <= 9007199254740991 = fromIntegral (fromIntegral w :: Int)
-    | otherwise = fromInteger (toInteger w)
+w2d w = fromIntegral (fromIntegral w :: Int)
 {-# INLINE w2d #-}
 
--- | Exactness in doubt: hand the raw slice to the reference parser.
+-- | Outside the exact window: hand the raw slice to the reference parser.
 referenceSlice :: BS.ByteString -> Int -> Int -> (# Int#, Double# #)
 referenceSlice bs start end =
     case readByteStringDouble (BSU.unsafeTake (end - start) (BSU.unsafeDrop start bs)) of
@@ -58,6 +57,9 @@ referenceSlice bs start end =
 parseDoubleField# :: BS.ByteString -> Int -> Int -> (# Int#, Double# #)
 parseDoubleField# bs start end0
     | i0 >= end = none
+    -- The one non-numeric token the reference accepts ('show' writes it).
+    | isInfinityAt i0 = done False (1 / 0)
+    | BSU.unsafeIndex bs i0 == 0x2D && isInfinityAt (i0 + 1) = done True (1 / 0)
     | otherwise =
         let !c0 = BSU.unsafeIndex bs i0
             !neg = c0 == 0x2D
@@ -69,14 +71,27 @@ parseDoubleField# bs start end0
                      in takeDigits64 bs iz end $ \wEnd w ->
                             if wEnd - iz > 19
                                 then referenceSlice bs start end0
-                                else afterWhole neg w wEnd
+                                else afterWhole neg w (wEnd - iz) wEnd
   where
     !i0 = skipStrip bs start end0
     !end = skipStripEnd bs i0 end0
 
     none = (# 0#, 0.0## #)
 
-    afterWhole !neg !w !i
+    -- "Infinity", exactly, to the end of the slice.
+    isInfinityAt i =
+        end - i == 8
+            && BSU.unsafeIndex bs i == 0x49
+            && BSU.unsafeIndex bs (i + 1) == 0x6E
+            && BSU.unsafeIndex bs (i + 2) == 0x66
+            && BSU.unsafeIndex bs (i + 3) == 0x69
+            && BSU.unsafeIndex bs (i + 4) == 0x6E
+            && BSU.unsafeIndex bs (i + 5) == 0x69
+            && BSU.unsafeIndex bs (i + 6) == 0x74
+            && BSU.unsafeIndex bs (i + 7) == 0x79
+
+    -- Fold the fraction into the significand so there is one scale step.
+    afterWhole !neg !w !wd !i
         | i < end && BSU.unsafeIndex bs i == 0x2E =
             let !f0 = i + 1
                 !fz = skipZeroes bs f0 end
@@ -84,13 +99,19 @@ parseDoubleField# bs start end0
                     if fEnd == f0
                         then none
                         else
-                            if fEnd - fz > 19
-                                then referenceSlice bs start end0
-                                else afterExponent neg (w2d w + (w2d p / pow10 (fEnd - f0))) fEnd
-        | otherwise = afterExponent neg (w2d w) i
-
-    afterExponent !neg !val !i
-        | i >= end = done neg val
+                            let !fracLen = fEnd - f0
+                             in if wd + fracLen > 19
+                                    then referenceSlice bs start end0
+                                    else
+                                        afterExponent
+                                            neg
+                                            (w * VU.unsafeIndex pow10w fracLen + p)
+                                            fracLen
+                                            fEnd
+        | otherwise = afterExponent neg w 0 i
+
+    afterExponent !neg !sig !fracLen !i
+        | i >= end = finish neg sig (negate fracLen)
         | BSU.unsafeIndex bs i == 0x65 || BSU.unsafeIndex bs i == 0x45 =
             let !i1 = i + 1
                 !eneg = i1 < end && BSU.unsafeIndex bs i1 == 0x2D
@@ -105,18 +126,29 @@ parseDoubleField# bs start end0
                                     else
                                         if eEnd - ez > 18
                                             then referenceSlice bs start end0
-                                            else done neg (val * scale eneg (fromIntegral e))
+                                            else
+                                                let !ev = fromIntegral e
+                                                 in finish
+                                                        neg
+                                                        sig
+                                                        ( (if eneg then negate ev else ev)
+                                                            - fracLen
+                                                        )
         | otherwise = none
 
-    scale !eneg !ex
-        | eneg = VU.unsafeIndex recipPow10Table k
-        | otherwise = VU.unsafeIndex pow10Table k
-      where
-        !k = min ex tableMax
-    {-# INLINE scale #-}
-
-    pow10 !k = VU.unsafeIndex pow10Table (min k tableMax)
-    {-# INLINE pow10 #-}
+    -- Clinger: exact significand times an exact power of ten is one
+    -- correctly-rounded op; anything wider goes to the reference.
+    finish !neg !sig !dexp
+        | sig == 0 = done neg 0
+        | sig <= 9007199254740991 && dexp >= negate exactPow10 && dexp <= exactPow10 =
+            done
+                neg
+                ( if dexp >= 0
+                    then w2d sig * VU.unsafeIndex pow10Table dexp
+                    else w2d sig / VU.unsafeIndex pow10Table (negate dexp)
+                )
+        | otherwise = referenceSlice bs start end0
+    {-# INLINE finish #-}
 
     done !neg !v = case if neg then negate v else v of
         D# d -> (# 1#, d #)
diff --git a/dataframe-parsing/tests/Properties/FastParsing.hs b/dataframe-parsing/tests/Properties/FastParsing.hs
index eb75dff2..7d207bd2 100644
--- a/dataframe-parsing/tests/Properties/FastParsing.hs
+++ b/dataframe-parsing/tests/Properties/FastParsing.hs
@@ -185,10 +185,12 @@ prop_doubleParity :: Field -> Property
 prop_doubleParity (Field f) =
     bitsOf (parseDoubleField f) === bitsOf (readByteStringDouble f)
 
+-- | 'show' output must parse back to the exact same double (NaN is nullish).
 prop_doubleShowRoundTrip :: Double -> Property
 prop_doubleShowRoundTrip d =
-    let f = C.pack (show d)
-     in bitsOf (parseDoubleField f) === bitsOf (readByteStringDouble f)
+    not (isNaN d) ==>
+        let f = C.pack (show d)
+         in bitsOf (parseDoubleField f) === bitsOf (Just d)
 
 prop_doubleSlice :: Field -> Field -> Field -> Property
 prop_doubleSlice (Field pre) (Field f) (Field post) =
diff --git a/dataframe-parsing/tests/Unit/FastParsing.hs b/dataframe-parsing/tests/Unit/FastParsing.hs
index b09c9676..0ce3771a 100644
--- a/dataframe-parsing/tests/Unit/FastParsing.hs
+++ b/dataframe-parsing/tests/Unit/FastParsing.hs
@@ -4,11 +4,13 @@
 
 Each case asserts parity with the reference parsers plus, where the
 parser-semantics contract (benchmarks/csv/results/code-audit.md section 2)
-pins a concrete value, the value itself.
+pins a concrete value, the value itself. Doubles are additionally pinned
+against 'read', which is correctly rounded.
 -}
 module Unit.FastParsing (tests) where
 
 import qualified Data.ByteString as BS
+import qualified Data.ByteString.Char8 as C
 
 import Data.Word (Word64)
 import GHC.Float (castDoubleToWord64)
@@ -94,6 +96,9 @@ doubleCases =
     , "1e-1024"
     , "1e-1025"
     , "Infinity"
+    , "-Infinity"
+    , "+Infinity"
+    , " Infinity "
     , "NaN"
     , "1e"
     , "1e+"
@@ -114,11 +119,28 @@ pinnedDoubles =
     , ("5.", Nothing)
     , ("1e3.5", Nothing)
     , ("1_000", Nothing)
-    , ("Infinity", Nothing)
+    , ("Infinity", Just (1 / 0))
+    , ("-Infinity", Just (-1 / 0))
     , ("NaN", Nothing)
-    , -- The reference parser is NOT correctly rounded here; bit-exact
-      -- parity means we must reproduce its 2.2250738585072e-308.
-      ("2.2250738585072011e-308", Just 2.2250738585072e-308)
+    , ("2.2250738585072011e-308", Just 2.225073858507201e-308)
+    ]
+
+-- | Correct rounding: every case must equal 'read' bit for bit.
+strtodDoubles :: [BS.ByteString]
+strtodDoubles =
+    [ "1.7976931348623157e308"
+    , "4.9406564584124654e-324"
+    , "1e-310"
+    , "1e-309"
+    , "1e-308"
+    , "2.2250871628565451e-249"
+    , "1.602176634e-19"
+    , "0.30000000000000004"
+    , "9007199254740993.9007199254740993"
+    , "123456789012345678901234567890e-25"
+    , "1e999"
+    , "1e-999"
+    , "-0"
     ]
 
 boolCases :: [BS.ByteString]
@@ -197,6 +219,14 @@ pinnedCase input expected =
             (bitsOf expected)
             (bitsOf (parseDoubleField input))
 
+readParityCase :: BS.ByteString -> Test
+readParityCase input =
+    TestLabel ("read parity: " ++ show input) . TestCase $
+        assertEqual
+            ("read parity on " ++ show input)
+            (bitsOf (Just (read (C.unpack input))))
+            (bitsOf (parseDoubleField input))
+
 tests :: [Test]
 tests =
     map (parityCase "int" parseIntField readByteStringInt) intCases
@@ -204,6 +234,7 @@ tests =
             (parityCase "double" (bitsOf . parseDoubleField) (bitsOf . readByteStringDouble))
             doubleCases
         ++ map (uncurry pinnedCase) pinnedDoubles
+        ++ map readParityCase strtodDoubles
         ++ map (parityCase "bool" parseBoolField readByteStringBool) boolCases
         ++ map (parityCase "missing" isMissingField isNullishBS) missingCases
         ++ map (parityCase "date" parseDateField (readByteStringDate "%Y-%m-%d")) dateCases
diff --git a/docs/exploratory_data_analysis_primer.md b/docs/exploratory_data_analysis_primer.md
index f97944eb..f6c5b4eb 100644
--- a/docs/exploratory_data_analysis_primer.md
+++ b/docs/exploratory_data_analysis_primer.md
@@ -249,7 +249,7 @@ D.skewness (F.col @Double "median_house_value") df
 
 ```
 
-> 0.977668529406543
+> 0.9776922140978362
 
 So the median house value is moderately skewed to the left. That is, there are more houses that are cheaper than the mean values and a tail of expensive outliers. Having lived in California, I can confirm that this data reflects reality.
 
diff --git a/tests/Internal/Parsing.hs b/tests/Internal/Parsing.hs
index 4ad13e53..512e31f1 100644
--- a/tests/Internal/Parsing.hs
+++ b/tests/Internal/Parsing.hs
@@ -153,6 +153,52 @@ readIntPartialSuffix =
             (readInt "42abc")
         )
 
+-- overflow must be a failed parse, not a wrapped value
+readIntOverflow :: Test
+readIntOverflow =
+    TestCase
+        ( assertEqual
+            "readInt \"9223372036854775808\" is Nothing"
+            Nothing
+            (readInt "9223372036854775808")
+        )
+
+readIntUnderflow :: Test
+readIntUnderflow =
+    TestCase
+        ( assertEqual
+            "readInt \"-9223372036854775809\" is Nothing"
+            Nothing
+            (readInt "-9223372036854775809")
+        )
+
+readIntWordWrap :: Test
+readIntWordWrap =
+    TestCase
+        ( assertEqual
+            "readInt \"18446744073709551616\" is Nothing"
+            Nothing
+            (readInt "18446744073709551616")
+        )
+
+readIntMaxBound :: Test
+readIntMaxBound =
+    TestCase
+        ( assertEqual
+            "readInt maxBound"
+            (Just (maxBound :: Int))
+            (readInt "9223372036854775807")
+        )
+
+readIntMinBound :: Test
+readIntMinBound =
+    TestCase
+        ( assertEqual
+            "readInt minBound"
+            (Just (minBound :: Int))
+            (readInt "-9223372036854775808")
+        )
+
 -- readDouble
 
 readDoublePositive :: Test
@@ -193,6 +239,48 @@ readDoublePartialSuffix =
             (readDouble "3.14abc")
         )
 
+-- correct rounding: parses must equal 'read' bit for bit
+readDoubleSubnormal :: Test
+readDoubleSubnormal =
+    TestCase
+        ( assertEqual
+            "readDouble on the smallest denormal"
+            (Just (read "4.9406564584124654e-324"))
+            (readDouble "4.9406564584124654e-324")
+        )
+
+readDoubleMaxFinite :: Test
+readDoubleMaxFinite =
+    TestCase
+        ( assertEqual
+            "readDouble on the largest finite double stays finite"
+            (Just (read "1.7976931348623157e308"))
+            (readDouble "1.7976931348623157e308")
+        )
+
+readDoubleOverflowIsInfinity :: Test
+readDoubleOverflowIsInfinity =
+    TestCase
+        (assertEqual "readDouble \"1e999\"" (Just (1 / 0)) (readDouble "1e999"))
+
+readDoubleInfinityToken :: Test
+readDoubleInfinityToken =
+    TestCase
+        ( assertEqual
+            "readDouble takes back what 'show' wrote"
+            (Just (-1 / 0))
+            (readDouble "-Infinity")
+        )
+
+readDoubleAbsurdExponent :: Test
+readDoubleAbsurdExponent =
+    TestCase
+        ( assertEqual
+            "readDouble \"1e18446744073709551617\" clamps, not allocates"
+            (Just (1 / 0))
+            (readDouble "1e18446744073709551617")
+        )
+
 tests :: [Test]
 tests =
     [ TestLabel "isNullishEmptyString" isNullishEmptyString
@@ -228,10 +316,20 @@ tests =
     , TestLabel "readIntText" readIntText
     , TestLabel "readIntEmpty" readIntEmpty
     , TestLabel "readIntPartialSuffix" readIntPartialSuffix
+    , TestLabel "readIntOverflow" readIntOverflow
+    , TestLabel "readIntUnderflow" readIntUnderflow
+    , TestLabel "readIntWordWrap" readIntWordWrap
+    , TestLabel "readIntMaxBound" readIntMaxBound
+    , TestLabel "readIntMinBound" readIntMinBound
     , TestLabel "readDoublePositive" readDoublePositive
     , TestLabel "readDoubleNegative" readDoubleNegative
     , TestLabel "readDoubleWholeNumber" readDoubleWholeNumber
     , TestLabel "readDoubleText" readDoubleText
     , TestLabel "readDoubleEmpty" readDoubleEmpty
     , TestLabel "readDoublePartialSuffix" readDoublePartialSuffix
+    , TestLabel "readDoubleSubnormal" readDoubleSubnormal
+    , TestLabel "readDoubleMaxFinite" readDoubleMaxFinite
+    , TestLabel "readDoubleOverflowIsInfinity" readDoubleOverflowIsInfinity
+    , TestLabel "readDoubleInfinityToken" readDoubleInfinityToken
+    , TestLabel "readDoubleAbsurdExponent" readDoubleAbsurdExponent
     ]
diff --git a/tests/Learn/MetricsTests.hs b/tests/Learn/MetricsTests.hs
index d43ccea6..2c018a17 100644
--- a/tests/Learn/MetricsTests.hs
+++ b/tests/Learn/MetricsTests.hs
@@ -43,6 +43,16 @@ testRegressionMetrics = TestCase $ do
     assertBool "rmse" (close 1e-9 (rmse p t) 0.5)
     assertBool "mae" (close 1e-9 (mae p t) 0.25)
     assertBool "r2 in range" (r2 p t <= 1)
+    -- zipWith truncates: the mean is over compared pairs, not all of truth.
+    assertBool
+        "mse averages over compared pairs"
+        (close 1e-9 (mse (VU.fromList [2, 2]) (VU.fromList [0, 0, 0, 0])) 4)
+    assertBool
+        "mae averages over compared pairs"
+        (close 1e-9 (mae (VU.fromList [2, 2]) (VU.fromList [0, 0, 0, 0])) 2)
+    assertBool
+        "no predictions is not a perfect score"
+        (isNaN (mse VU.empty (VU.fromList [5, 5, 5])))
 
 testMulticlassMetrics :: Test
 testMulticlassMetrics = TestCase $ do
diff --git a/tests/Operations/ParallelGroupBy.hs b/tests/Operations/ParallelGroupBy.hs
index 4b5b9310..5cff0b02 100644
--- a/tests/Operations/ParallelGroupBy.hs
+++ b/tests/Operations/ParallelGroupBy.hs
@@ -109,7 +109,12 @@ aggParityFor n =
                 ]
             seqDf = D.aggregate aggs (groupBySeq ["ki", "kt"] df)
             parDf = D.aggregate aggs (groupByPar ["ki", "kt"] df)
-         in assertEqual ("aggregate parity n=" ++ show n) seqDf parDf
+         in -- Rendered comparison: singleton-group stddev is NaN, which the
+            -- Eq instance treats as unequal.
+            assertEqual
+                ("aggregate parity n=" ++ show n)
+                (D.toMarkdown seqDf)
+                (D.toMarkdown parDf)
 
 collisionParity :: Test
 collisionParity =
diff --git a/tests/Operations/Statistics.hs b/tests/Operations/Statistics.hs
index db907654..9cbb509e 100644
--- a/tests/Operations/Statistics.hs
+++ b/tests/Operations/Statistics.hs
@@ -56,6 +56,7 @@ skewnessOfSymmetricDataSet =
             0
         )
 
+-- Population skewness g1, the form the docs define (matches scipy.stats.skew).
 skewnessOfSimpleDataSet :: Test
 skewnessOfSimpleDataSet =
     TestCase
@@ -63,7 +64,7 @@ skewnessOfSimpleDataSet =
             "Skewness of a simple data set"
             ( abs
                 ( D.skewness' (VU.fromList [25 :: Int, 28, 26, 30, 40, 50, 40])
-                    - 0.566_731_633_676
+                    - 0.612_140_127_240_396_6
                 )
                 < 1e-12
             )
@@ -181,6 +182,57 @@ wrongQuantileIndex =
             (print $ D.quantiles' (VU.fromList [5]) 4 (VU.fromList [1 :: Int, 2, 3, 4, 5]))
         )
 
+-- Int aggregation must widen before summing, not wrap at 2^63.
+medianOfLargeIntDataSet :: Test
+medianOfLargeIntDataSet =
+    TestCase
+        ( assertEqual
+            "Median of an even length Int data set near maxBound"
+            9.223372036854776e18
+            (D.median' (VU.fromList [maxBound - 1, maxBound :: Int]))
+        )
+
+meanOfLargeIntDataSet :: Test
+meanOfLargeIntDataSet =
+    TestCase
+        ( assertEqual
+            "Mean of an Int data set summing past 2^63"
+            9.223372036854776e18
+            (D.meanInt' (VU.fromList [maxBound, maxBound :: Int]))
+        )
+
+-- The one-pass correlation form cancelled catastrophically here (gave 2.0).
+correlationLowVarianceBounded :: Test
+correlationLowVarianceBounded =
+    TestCase
+        ( let df =
+                D.fromNamedColumns
+                    [ ("a", DI.fromList [1e8 :: Double, 1e8, 1.00000002e8])
+                    , ("b", DI.fromList [1e8 :: Double, 1e8, 1.00000003e8])
+                    ]
+           in case D.correlation "a" "b" df of
+                Nothing -> assertFailure "Expected Just 1.0, got Nothing"
+                Just r ->
+                    assertBool
+                        "collinear large-offset columns give r = 1"
+                        (abs (r - 1.0) < 1e-9)
+        )
+
+-- Self-correlation with a large offset flipped sign (gave -1.0).
+correlationSelfLargeOffset :: Test
+correlationSelfLargeOffset =
+    TestCase
+        ( let df =
+                D.fromNamedColumns
+                    [("a", DI.fromList [1e11 :: Double, 1e11 + 1, 1e11 + 2])]
+           in case D.correlation "a" "a" df of
+                Nothing -> assertFailure "Expected Just 1.0, got Nothing"
+                Just r ->
+                    assertBool
+                        "self correlation at a large offset is 1"
+                        (abs (r - 1.0) < 1e-9)
+        )
+
 summarizeOptional :: Test
 summarizeOptional =
     TestCase
@@ -276,6 +328,10 @@ tests =
         interQuartileRangeOfEvenLengthDataSet
     , TestLabel "wrongQuantileNumber" wrongQuantileNumber
     , TestLabel "wrongQuantileIndex" wrongQuantileIndex
+    , TestLabel "medianOfLargeIntDataSet" medianOfLargeIntDataSet
+    , TestLabel "meanOfLargeIntDataSet" meanOfLargeIntDataSet
+    , TestLabel "correlationLowVarianceBounded" correlationLowVarianceBounded
+    , TestLabel "correlationSelfLargeOffset" correlationSelfLargeOffset
     , TestLabel "summarizeOptional" summarizeOptional
     , TestLabel "correlationPerfectPositive" correlationPerfectPositive
     , TestLabel "correlationPerfectNegative" correlationPerfectNegative
diff --git a/tests/Operations/VectorKernel.hs b/tests/Operations/VectorKernel.hs
index b9c898cb..cc8f5247 100644
--- a/tests/Operations/VectorKernel.hs
+++ b/tests/Operations/VectorKernel.hs
@@ -123,7 +123,12 @@ parityCase n keys aggs =
             fast = D.aggregate aggs gdf
             ref = interpretOnly aggs gdf
             label = "n=" ++ show n ++ " keys=" ++ show keys ++ " #aggs=" ++ show (length aggs)
-         in assertEqual ("kernel==interpreter " ++ label) ref fast
+         in -- Full render: value-exact via 'show', and NaN cells (singleton
+            -- group stddev/variance) compare equal, which 'Eq' refuses.
+            assertEqual
+                ("kernel==interpreter " ++ label)
+                (D.toMarkdown ref)
+                (D.toMarkdown fast)
 
 tests :: [Test]
 tests =