From 74ae891766d5bd015ce881c0d305ecbf5328a438 Mon Sep 17 00:00:00 2001
From: skymanbp <57272723+skymanbp@users.noreply.github.com>
Date: Wed, 19 Aug 2026 03:00:05 -0400
Subject: [PATCH 1/5] fix: statistics ignored the null bitmap

Stats read the raw payload of nullable columns, so the sentinel stored
at null slots entered every result (a zero for unboxed columns, an
error thunk for boxed ones):

- mean/median/percentile/genericPercentile/stddev/skewness/variance/
  IQR/sum now drop null slots via a shared dropNulls view.
- correlation does listwise deletion over both columns.
- valueCounts/valueProportions (and so frequencies) no longer count
  the sentinel as a category.
- Maybe-typed views are untouched: their nulls are real values.
---
 .../src-internal/DataFrame/Internal/Column.hs |  18 +-
 .../src/DataFrame/Operations/Core.hs          |  20 +-
 .../src/DataFrame/Operations/Statistics.hs    |  64 +++---
 tests/Operations/Statistics.hs                | 183 ++++++++++++++++++
 4 files changed, 254 insertions(+), 31 deletions(-)

diff --git a/dataframe-core/src-internal/DataFrame/Internal/Column.hs b/dataframe-core/src-internal/DataFrame/Internal/Column.hs
index 1395b9b4..6f76f1df 100644
--- a/dataframe-core/src-internal/DataFrame/Internal/Column.hs
+++ b/dataframe-core/src-internal/DataFrame/Internal/Column.hs
@@ -39,7 +39,7 @@ import Data.Bits (
  )
 import Data.Kind (Type)
 import Data.Maybe
-import Data.Type.Equality (TestEquality (..))
+import Data.Type.Equality (TestEquality (..), type (:~~:) (HRefl))
 import Data.Word (Word8)
 import DataFrame.Errors
 import DataFrame.Internal.PackedText (
@@ -209,6 +209,22 @@ columnBitmap (UnboxedColumn bm _) = bm
 columnBitmap (PackedText bm _) = bm
 columnBitmap (MergedColumn _ _) = Nothing
 
+{- | Drop the null slots of a payload-typed view of a nullable column: those
+slots hold a sentinel, not a value. A @Maybe@-typed view already encodes the
+nulls, so it is returned untouched. Backpermute never forces the kept-out
+slots, so boxed error thunks at null slots are safe.
+-}
+dropNulls ::
+    forall v a.
+    (Typeable a, VG.Vector v a, VG.Vector v Int) => Maybe Bitmap -> v a -> v a
+dropNulls Nothing xs = xs
+dropNulls (Just bm) xs = case typeRep @a of
+    App m _ | Just HRefl <- eqTypeRep m (typeRep @Maybe) -> xs
+    _ -> VG.backpermute xs keep
+  where
+    keep = VG.fromList [i | i <- [0 .. VG.length xs - 1], bitmapTestBit bm i]
+{-# INLINE dropNulls #-}
+
 {- | Decode a 'PackedText' into a @BoxedColumn Text@ (bit-identical to
 materializing at freeze). Identity on every other column.
 -}
diff --git a/dataframe-operations/src/DataFrame/Operations/Core.hs b/dataframe-operations/src/DataFrame/Operations/Core.hs
index 77df934f..ef7e66c9 100644
--- a/dataframe-operations/src/DataFrame/Operations/Core.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Core.hs
@@ -67,8 +67,10 @@ import DataFrame.Internal.Column (
     Column (..),
     Columnable,
     TypedColumn (..),
+    columnBitmap,
     columnLength,
     columnTypeString,
+    dropNulls,
     fromList,
     fromVector,
     materializeMerged,
@@ -688,7 +690,7 @@ valueCounts ::
     forall a. (Ord a, Columnable a) => Expr a -> DataFrame -> [(a, Int)]
 valueCounts expr df
     | null df = throw (EmptyDataSetException "valueCounts")
-    | otherwise = case columnAsVector expr df of
+    | otherwise = case columnAsVectorNonNull expr df of
         Left e -> throw e
         Right column' ->
             let
@@ -696,6 +698,20 @@ valueCounts expr df
              in
                 M.toAscList column
 
+-- | As 'columnAsVector', minus null slots: a sentinel is not a category.
+columnAsVectorNonNull ::
+    forall a.
+    (Columnable a) => Expr a -> DataFrame -> Either DataFrameException (V.Vector a)
+columnAsVectorNonNull expr df = case expr of
+    Col name -> case getColumn name df of
+        Just col -> withColumnName name (dropNulls (columnBitmap col) <$> toVector col)
+        Nothing ->
+            Left $
+                ColumnsNotFoundException [name] "valueCounts" (M.keys $ columnIndices df)
+    _ -> case interpret df expr of
+        Left e -> throw e
+        Right (TColumn col) -> dropNulls (columnBitmap col) <$> toVector col
+
 {- | O (k * n) Shows the proportions of each value in a given column.
 
 ==== __Example__
@@ -712,7 +728,7 @@ valueProportions ::
     forall a. (Ord a, Columnable a) => Expr a -> DataFrame -> [(a, Double)]
 valueProportions expr df
     | null df = throw (EmptyDataSetException "valueCounts")
-    | otherwise = case columnAsVector expr df of
+    | otherwise = case columnAsVectorNonNull expr df of
         Left e -> throw e
         Right column' ->
             let
diff --git a/dataframe-operations/src/DataFrame/Operations/Statistics.hs b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
index 7322c0dc..c67d6d96 100644
--- a/dataframe-operations/src/DataFrame/Operations/Statistics.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
@@ -120,7 +120,7 @@ mean expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> mean' xs
+        Right xs -> mean' (dropNulls (columnBitmap col) xs)
 
 meanMaybe ::
     forall a. (Columnable a, Real a) => Expr (Maybe a) -> DataFrame -> Double
@@ -137,13 +137,13 @@ meanMaybe expr df = case interpret @(Maybe a) df expr of
 median ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
 median (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> median' xs
+    Right xs -> median' (dropNulls (colBitmap name df) xs)
     Left e -> throw e
 median expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> median' xs
+        Right xs -> median' (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the median of a given column (containing optional values) as a standalone value.
 medianMaybe ::
@@ -162,50 +162,50 @@ percentile ::
     forall a.
     (Columnable a, Real a, VU.Unbox a) => Int -> Expr a -> DataFrame -> Double
 percentile n (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> percentile' n xs
+    Right xs -> percentile' n (dropNulls (colBitmap name df) xs)
     Left e -> throw e
 percentile n expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> percentile' n xs
+        Right xs -> percentile' n (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the nth percentile of a given column as a standalone value.
 genericPercentile ::
     forall a.
     (Columnable a, Ord a) => Int -> Expr a -> DataFrame -> a
 genericPercentile n (Col name) df = case columnAsVector (Col @a name) df of
-    Right xs -> percentileOrd' n xs
+    Right xs -> percentileOrd' n (dropNulls (colBitmap name df) xs)
     Left e -> throw e
 genericPercentile n expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toVector @a col of
         Left e -> throw e
-        Right xs -> percentileOrd' n xs
+        Right xs -> percentileOrd' n (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the standard deviation of a given column as a standalone value.
 standardDeviation ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
 standardDeviation (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> (sqrt . variance') xs
+    Right xs -> (sqrt . variance') (dropNulls (colBitmap name df) xs)
     Left e -> throw e
 standardDeviation expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> (sqrt . variance') xs
+        Right xs -> (sqrt . variance') (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the skewness of a given column as a standalone value.
 skewness ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
 skewness (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> skewness' xs
+    Right xs -> skewness' (dropNulls (colBitmap name df) xs)
     Left e -> throw e
 skewness expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> skewness' xs
+        Right xs -> skewness' (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the variance of a given column as a standalone value.
 variance ::
@@ -217,36 +217,44 @@ variance expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> variance' xs
+        Right xs -> variance' (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the inter-quartile range of a given column as a standalone value.
 interQuartileRange ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
 interQuartileRange (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> interQuartileRange' xs
+    Right xs -> interQuartileRange' (dropNulls (colBitmap name df) xs)
     Left e -> throw e
 interQuartileRange expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn col) -> case toUnboxedVector @a col of
         Left e -> throw e
-        Right xs -> interQuartileRange' xs
+        Right xs -> interQuartileRange' (dropNulls (columnBitmap col) xs)
 
 -- | Calculates the Pearson's correlation coefficient between two given columns as a standalone value.
 correlation :: T.Text -> T.Text -> DataFrame -> Maybe Double
 correlation first second df = do
-    f <- _getColumnAsDouble first df
-    s <- _getColumnAsDouble second df
+    -- Listwise deletion: a null in either column drops the pair.
+    let df' = filterJust first (filterJust second df)
+    f <- _getColumnAsDouble first df'
+    s <- _getColumnAsDouble second df'
     correlation' f s
 
+-- | Bitmap of a named column, if it has one.
+colBitmap :: T.Text -> DataFrame -> Maybe Bitmap
+colBitmap name df = columnBitmap =<< getColumn name df
+
 _getColumnAsDouble :: T.Text -> DataFrame -> Maybe (VU.Vector Double)
 _getColumnAsDouble name df = case getColumn name df of
-    Just (UnboxedColumn _ (f :: VU.Vector a)) -> case testEquality (typeRep @a) (typeRep @Double) of
-        Just Refl -> Just f
-        Nothing -> case sIntegral @a of
-            STrue -> Just (VU.map fromIntegral f)
-            SFalse -> case sFloating @a of
-                STrue -> Just (VU.map realToFrac f)
-                SFalse -> Nothing
+    Just (UnboxedColumn bm (f' :: VU.Vector a)) ->
+        let f = dropNulls bm f'
+         in case testEquality (typeRep @a) (typeRep @Double) of
+                Just Refl -> Just f
+                Nothing -> case sIntegral @a of
+                    STrue -> Just (VU.map fromIntegral f)
+                    SFalse -> case sFloating @a of
+                        STrue -> Just (VU.map realToFrac f)
+                        SFalse -> Nothing
     Nothing ->
         throw $
             ColumnsNotFoundException [name] "_getColumnAsDouble" (M.keys $ columnIndices df)
@@ -265,11 +273,11 @@ sum ::
     forall a. (Columnable a, Num a) => Expr a -> DataFrame -> a
 sum (Col name) df = case getColumn name df of
     Nothing -> throw $ ColumnsNotFoundException [name] "sum" (M.keys $ columnIndices df)
-    Just ((UnboxedColumn _ (column :: VU.Vector a'))) -> case testEquality (typeRep @a') (typeRep @a) of
-        Just Refl -> VG.sum column
+    Just ((UnboxedColumn bm (column :: VU.Vector a'))) -> case testEquality (typeRep @a') (typeRep @a) of
+        Just Refl -> VG.sum (dropNulls bm column)
         Nothing -> 0
-    Just ((BoxedColumn _ (column :: V.Vector a'))) -> case testEquality (typeRep @a') (typeRep @a) of
-        Just Refl -> VG.sum column
+    Just ((BoxedColumn bm (column :: V.Vector a'))) -> case testEquality (typeRep @a') (typeRep @a) of
+        Just Refl -> VG.sum (dropNulls bm column)
         Nothing -> 0
     Just (PackedText _ _) -> 0
     Just (MergedColumn _ _) -> 0 -- matches the old eager These column (type never Num)
@@ -277,7 +285,7 @@ sum expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn xs) -> case toVector @a @V.Vector xs of
         Left e -> throw e
-        Right xs' -> VG.sum xs'
+        Right xs' -> VG.sum (dropNulls (columnBitmap xs) xs')
 
 {- | /O(n)/ Impute missing values in a column using a derived scalar.
 
diff --git a/tests/Operations/Statistics.hs b/tests/Operations/Statistics.hs
index db907654..dc67124f 100644
--- a/tests/Operations/Statistics.hs
+++ b/tests/Operations/Statistics.hs
@@ -6,6 +6,7 @@ module Operations.Statistics where
 
 import qualified Data.Vector.Unboxed as VU
 import qualified DataFrame as D
+import qualified DataFrame.Functions as F
 import qualified DataFrame.Internal.Column as DI
 import qualified DataFrame.Internal.Statistics as D
 
@@ -254,6 +255,168 @@ correlationMissingColumn =
             (print $ D.correlation "x" "missingcol" correlationDf)
         )
 
+-- Nullable columns: statistics must skip null slots, not read the sentinel
+-- stored there.
+
+nullableDf :: D.DataFrame
+nullableDf =
+    D.fromNamedColumns [("x", DI.fromList [Just (10 :: Double), Nothing, Just 20])]
+
+nullableIntDf :: D.DataFrame
+nullableIntDf =
+    D.fromNamedColumns [("n", DI.fromList [Just (10 :: Int), Nothing, Just 20])]
+
+nullableBoxedDf :: D.DataFrame
+nullableBoxedDf =
+    D.fromNamedColumns [("b", DI.fromList [Just (10 :: Integer), Nothing, Just 20])]
+
+meanIgnoresNulls :: Test
+meanIgnoresNulls =
+    TestCase
+        (assertEqual "mean skips nulls" 15.0 (D.mean (F.col @Double "x") nullableDf))
+
+meanExprIgnoresNulls :: Test
+meanExprIgnoresNulls =
+    TestCase
+        ( assertEqual
+            "mean over a derived nullable expression skips nulls"
+            30.0
+            (D.mean (F.lift (* 2) (F.col @Double "x")) nullableDf)
+        )
+
+medianIgnoresNulls :: Test
+medianIgnoresNulls =
+    TestCase
+        (assertEqual "median skips nulls" 15.0 (D.median (F.col @Double "x") nullableDf))
+
+percentileIgnoresNulls :: Test
+percentileIgnoresNulls =
+    TestCase
+        ( assertEqual
+            "percentile skips nulls"
+            15.0
+            (D.percentile 50 (F.col @Double "x") nullableDf)
+        )
+
+stdDevIgnoresNulls :: Test
+stdDevIgnoresNulls =
+    TestCase
+        ( assertBool
+            "standard deviation skips nulls"
+            ( abs (D.standardDeviation (F.col @Double "x") nullableDf - 7.0710678118654755)
+                < 1e-12
+            )
+        )
+
+varianceIgnoresNulls :: Test
+varianceIgnoresNulls =
+    TestCase
+        ( assertEqual
+            "variance skips nulls"
+            50.0
+            (D.variance (F.col @Double "x") nullableDf)
+        )
+
+varianceExprIgnoresNulls :: Test
+varianceExprIgnoresNulls =
+    TestCase
+        ( assertEqual
+            "variance over a derived nullable expression skips nulls"
+            200.0
+            (D.variance (F.lift (* 2) (F.col @Double "x")) nullableDf)
+        )
+
+iqrIgnoresNulls :: Test
+iqrIgnoresNulls =
+    TestCase
+        ( assertEqual
+            "inter-quartile range skips nulls"
+            5.0
+            (D.interQuartileRange (F.col @Double "x") nullableDf)
+        )
+
+skewnessIgnoresNulls :: Test
+skewnessIgnoresNulls =
+    TestCase
+        ( let skewDf =
+                D.fromNamedColumns
+                    [("s", DI.fromList [Just (10 :: Double), Nothing, Just 20, Just 100, Just 11])]
+           in assertBool
+                "skewness skips nulls"
+                ( abs
+                    ( D.skewness (F.col @Double "s") skewDf
+                        - D.skewness' (VU.fromList [10 :: Double, 20, 100, 11])
+                    )
+                    < 1e-12
+                )
+        )
+
+genericPercentileIgnoresNulls :: Test
+genericPercentileIgnoresNulls =
+    TestCase
+        ( assertEqual
+            "genericPercentile skips the sentinel"
+            10
+            (D.genericPercentile 10 (F.col @Int "n") nullableIntDf)
+        )
+
+genericPercentileBoxedNullableDoesNotThrow :: Test
+genericPercentileBoxedNullableDoesNotThrow =
+    TestCase
+        ( assertEqual
+            "genericPercentile on a boxed nullable column skips the error thunk"
+            20
+            (D.genericPercentile 100 (F.col @Integer "b") nullableBoxedDf)
+        )
+
+genericPercentileMaybeViewKeepsNothing :: Test
+genericPercentileMaybeViewKeepsNothing =
+    TestCase
+        ( assertEqual
+            "a Maybe-typed view still sees its Nothings"
+            (Nothing :: Maybe Int)
+            (D.genericPercentile 0 (F.col @(Maybe Int) "n") nullableIntDf)
+        )
+
+sumUnboxedNullable :: Test
+sumUnboxedNullable =
+    TestCase
+        (assertEqual "sum skips null slots" 30.0 (D.sum (F.col @Double "x") nullableDf))
+
+sumBoxedNullableDoesNotThrow :: Test
+sumBoxedNullableDoesNotThrow =
+    TestCase
+        ( assertEqual
+            "sum on a boxed nullable column skips the error thunk"
+            (30 :: Integer)
+            (D.sum (F.col @Integer "b") nullableBoxedDf)
+        )
+
+correlationIgnoresNullRows :: Test
+correlationIgnoresNullRows =
+    TestCase
+        ( let dfc =
+                D.fromNamedColumns
+                    [ ("a", DI.fromList [Just (1 :: Double), Nothing, Just 3])
+                    , ("c", DI.fromList [1 :: Double, 2, 3])
+                    ]
+           in case D.correlation "a" "c" dfc of
+                Nothing -> assertFailure "Expected Just 1.0, got Nothing"
+                Just r ->
+                    assertBool
+                        "null rows are dropped pairwise"
+                        (abs (r - 1.0) < 1e-10)
+        )
+
+frequenciesSkipsNulls :: Test
+frequenciesSkipsNulls =
+    TestCase
+        ( assertEqual
+            "frequencies has no sentinel category"
+            3 -- Statistic, 10 and 20
+            (D.nColumns (D.frequencies (F.col @Int "n") nullableIntDf))
+        )
+
 tests :: [Test]
 tests =
     [ TestLabel "medianOfOddLengthDataSet" medianOfOddLengthDataSet
@@ -281,4 +444,24 @@ tests =
     , TestLabel "correlationPerfectNegative" correlationPerfectNegative
     , TestLabel "correlationSelfIdentity" correlationSelfIdentity
     , TestLabel "correlationMissingColumn" correlationMissingColumn
+    , TestLabel "meanIgnoresNulls" meanIgnoresNulls
+    , TestLabel "meanExprIgnoresNulls" meanExprIgnoresNulls
+    , TestLabel "medianIgnoresNulls" medianIgnoresNulls
+    , TestLabel "percentileIgnoresNulls" percentileIgnoresNulls
+    , TestLabel "stdDevIgnoresNulls" stdDevIgnoresNulls
+    , TestLabel "varianceIgnoresNulls" varianceIgnoresNulls
+    , TestLabel "varianceExprIgnoresNulls" varianceExprIgnoresNulls
+    , TestLabel "iqrIgnoresNulls" iqrIgnoresNulls
+    , TestLabel "skewnessIgnoresNulls" skewnessIgnoresNulls
+    , TestLabel "genericPercentileIgnoresNulls" genericPercentileIgnoresNulls
+    , TestLabel
+        "genericPercentileBoxedNullableDoesNotThrow"
+        genericPercentileBoxedNullableDoesNotThrow
+    , TestLabel
+        "genericPercentileMaybeViewKeepsNothing"
+        genericPercentileMaybeViewKeepsNothing
+    , TestLabel "sumUnboxedNullable" sumUnboxedNullable
+    , TestLabel "sumBoxedNullableDoesNotThrow" sumBoxedNullableDoesNotThrow
+    , TestLabel "correlationIgnoresNullRows" correlationIgnoresNullRows
+    , TestLabel "frequenciesSkipsNulls" frequenciesSkipsNulls
     ]

From cbefe71687e21ca9b0a623e96d4769c8228af55a Mon Sep 17 00:00:00 2001
From: skymanbp <57272723+skymanbp@users.noreply.github.com>
Date: Wed, 19 Aug 2026 03:34:03 -0400
Subject: [PATCH 2/5] fix: numeric correctness in statistics, metrics and
 parsing

- Int mean/median widen per element instead of wrapping at 2^63.
- correlation uses the centered two-pass form, clamped: the one-pass
  form returned |r| > 1, NaN, or a flipped sign on offset data.
- skewness computes population g1, the formula the docs define; the
  old factor was off by sqrt((n-1)/n).
- variance of n < 2 is NaN (scatter kernels included), not a fake 0
  that made singleton groups look constant.
- mse/mae/r2/accuracy/logLoss average over compared pairs and refuse
  to score zero predictions; meanSquaredError guards length mismatch.
- readInt rejects overflow instead of wrapping.
- Double parsing is correctly rounded: exact byte-level reference
  (clamped exponents), Clinger fast path, Infinity round-trips.
---
 .../DataFrame/Internal/AggKernel.hs           |   4 +-
 .../DataFrame/Internal/AggKernelPar.hs        |   4 +-
 dataframe-learn/src/DataFrame/Metrics.hs      |  39 ++++--
 .../tests-internal/Learn/EdgeCases.hs         |  35 ++++-
 .../DataFrame/Internal/Statistics.hs          |  63 +++++----
 .../src/DataFrame/Operations/Inference.hs     |  10 +-
 .../src/DataFrame/Operations/Statistics.hs    |   5 +-
 .../src/DataFrame/Operations/Typing.hs        |   4 +-
 dataframe-parsing/dataframe-parsing.cabal     |   1 -
 .../DataFrame/Internal/Parsing.hs             | 120 ++++++++++++++----
 .../DataFrame/Internal/Parsing/Fast/Double.hs | 108 ++++++++++------
 .../tests/Properties/FastParsing.hs           |   6 +-
 dataframe-parsing/tests/Unit/FastParsing.hs   |  41 +++++-
 docs/exploratory_data_analysis_primer.md      |   2 +-
 tests/Internal/Parsing.hs                     |  98 ++++++++++++++
 tests/Learn/MetricsTests.hs                   |  10 ++
 tests/Operations/ParallelGroupBy.hs           |   7 +-
 tests/Operations/Statistics.hs                |  58 ++++++++-
 tests/Operations/VectorKernel.hs              |   7 +-
 19 files changed, 487 insertions(+), 135 deletions(-)

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 =

From 1aa33291d55d1af6648f4cb494a9cb393f8a761b Mon Sep 17 00:00:00 2001
From: skymanbp <57272723+skymanbp@users.noreply.github.com>
Date: Wed, 19 Aug 2026 12:26:40 -0400
Subject: [PATCH 3/5] Address review: dropNulls on Column, metrics throw,
 revert hot-path changes

---
 .../src-internal/DataFrame/Internal/Column.hs |  34 +-
 dataframe-learn/src/DataFrame/Metrics.hs      |  19 +-
 .../DataFrame/Internal/Statistics.hs          |  43 ++-
 .../src/DataFrame/Operations/Core.hs          |  19 +-
 .../src/DataFrame/Operations/Inference.hs     |  10 +-
 .../src/DataFrame/Operations/Statistics.hs    | 122 ++++---
 .../src/DataFrame/Operations/Typing.hs        |   4 +-
 dataframe-parsing/dataframe-parsing.cabal     |   1 +
 .../DataFrame/Internal/Parsing.hs             | 120 ++-----
 .../DataFrame/Internal/Parsing/Fast/Double.hs | 108 ++----
 .../tests/Properties/FastParsing.hs           |   6 +-
 dataframe-parsing/tests/Unit/FastParsing.hs   |  41 +--
 dataframe.cabal                               |   1 -
 docs/exploratory_data_analysis_primer.md      |   2 +-
 tests/Internal/Parsing.hs                     | 335 ------------------
 tests/Learn/MetricsTests.hs                   |  11 +-
 tests/Main.hs                                 |   2 -
 tests/Operations/Statistics.hs                |  55 ---
 18 files changed, 221 insertions(+), 712 deletions(-)
 delete mode 100644 tests/Internal/Parsing.hs

diff --git a/dataframe-core/src-internal/DataFrame/Internal/Column.hs b/dataframe-core/src-internal/DataFrame/Internal/Column.hs
index 6f76f1df..19784282 100644
--- a/dataframe-core/src-internal/DataFrame/Internal/Column.hs
+++ b/dataframe-core/src-internal/DataFrame/Internal/Column.hs
@@ -39,7 +39,7 @@ import Data.Bits (
  )
 import Data.Kind (Type)
 import Data.Maybe
-import Data.Type.Equality (TestEquality (..), type (:~~:) (HRefl))
+import Data.Type.Equality (TestEquality (..))
 import Data.Word (Word8)
 import DataFrame.Errors
 import DataFrame.Internal.PackedText (
@@ -209,22 +209,28 @@ columnBitmap (UnboxedColumn bm _) = bm
 columnBitmap (PackedText bm _) = bm
 columnBitmap (MergedColumn _ _) = Nothing
 
-{- | Drop the null slots of a payload-typed view of a nullable column: those
-slots hold a sentinel, not a value. A @Maybe@-typed view already encodes the
-nulls, so it is returned untouched. Backpermute never forces the kept-out
-slots, so boxed error thunks at null slots are safe.
+{- | Drops the null values in a nullable column: those slots hold a sentinel,
+not a value. Identity on columns without a bitmap. 'VG.ifilter' inspects only
+the index, so boxed error thunks at null slots are never forced.
 -}
-dropNulls ::
-    forall v a.
-    (Typeable a, VG.Vector v a, VG.Vector v Int) => Maybe Bitmap -> v a -> v a
-dropNulls Nothing xs = xs
-dropNulls (Just bm) xs = case typeRep @a of
-    App m _ | Just HRefl <- eqTypeRep m (typeRep @Maybe) -> xs
-    _ -> VG.backpermute xs keep
-  where
-    keep = VG.fromList [i | i <- [0 .. VG.length xs - 1], bitmapTestBit bm i]
+dropNulls :: Column -> Column
+dropNulls (BoxedColumn (Just bm) xs) =
+    BoxedColumn Nothing (VG.ifilter (\i _ -> bitmapTestBit bm i) xs)
+dropNulls (UnboxedColumn (Just bm) xs) =
+    UnboxedColumn Nothing (VG.ifilter (\i _ -> bitmapTestBit bm i) xs)
+dropNulls c@(PackedText (Just _) _) = dropNulls (materializePacked c)
+dropNulls c = c
 {-# INLINE dropNulls #-}
 
+{- | 'dropNulls', unless the view type @a@ is @Maybe@-headed: a @Maybe@-typed
+view encodes the nulls as values, so the column passes through untouched.
+-}
+dropNullsExceptMaybe :: forall a. (Typeable a) => Column -> Column
+dropNullsExceptMaybe c = case typeRep @a of
+    App m _ | Just HRefl <- eqTypeRep m (typeRep @Maybe) -> c
+    _ -> dropNulls c
+{-# INLINE dropNullsExceptMaybe #-}
+
 {- | Decode a 'PackedText' into a @BoxedColumn Text@ (bit-identical to
 materializing at freeze). Identity on every other column.
 -}
diff --git a/dataframe-learn/src/DataFrame/Metrics.hs b/dataframe-learn/src/DataFrame/Metrics.hs
index e7be7844..eb60725e 100644
--- a/dataframe-learn/src/DataFrame/Metrics.hs
+++ b/dataframe-learn/src/DataFrame/Metrics.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE OverloadedStrings #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 
@@ -45,6 +46,7 @@ import Data.Ord (comparing)
 import qualified Data.Text as T
 import qualified Data.Vector.Unboxed as VU
 
+import DataFrame.Errors (DataFrameException (..))
 import DataFrame.Internal.Column (TypedColumn (..), toVector)
 import DataFrame.Internal.DataFrame (DataFrame)
 import DataFrame.Internal.Expression (Expr)
@@ -82,9 +84,8 @@ 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
+    -- No compared pairs is not a score of any kind.
+    | n == 0 = throw (EmptyDataSetException "mse")
     | otherwise =
         VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n
   where
@@ -97,8 +98,7 @@ rmse preds truth = sqrt (mse preds truth)
 -- | Mean absolute error.
 mae :: Metric
 mae preds truth
-    | VU.null truth = 0
-    | n == 0 = 0 / 0
+    | n == 0 = throw (EmptyDataSetException "mae")
     | otherwise = VU.sum (VU.zipWith (\p t -> abs (p - t)) preds truth) / n
   where
     n = nCompared preds truth
@@ -106,8 +106,7 @@ mae preds truth
 -- | Coefficient of determination @R²@.
 r2 :: Metric
 r2 preds truth
-    | VU.null truth = 0
-    | n == 0 = 0 / 0
+    | n == 0 = throw (EmptyDataSetException "r2")
     | ssTot == 0 = 0
     | otherwise = 1 - ssRes / ssTot
   where
@@ -120,8 +119,7 @@ r2 preds truth
 -- | Fraction of exact matches.
 accuracy :: Metric
 accuracy preds truth
-    | VU.null truth = 0
-    | n == 0 = 0 / 0
+    | n == 0 = throw (EmptyDataSetException "accuracy")
     | otherwise =
         fromIntegral (VU.length (VU.filter id (VU.zipWith (==) preds truth))) / n
   where
@@ -130,8 +128,7 @@ accuracy preds truth
 -- | Binary log loss; probabilities clamped away from @0@/@1@.
 logLoss :: Metric
 logLoss probs truth
-    | VU.null truth = 0
-    | n == 0 = 0 / 0
+    | n == 0 = throw (EmptyDataSetException "logLoss")
     | otherwise =
         negate
             ( VU.sum
diff --git a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
index 190e4ec2..d7403bf9 100644
--- a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
+++ b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
@@ -17,9 +17,7 @@ import DataFrame.Errors (DataFrameException (..))
 mean' :: (Real a, VU.Unbox a) => VU.Vector a -> Double
 mean' samp
     | VU.null samp = throw $ EmptyDataSetException "mean"
-    -- Widen per element: summing at 'a' wraps for Int columns.
-    | otherwise =
-        VU.foldl' (\acc x -> acc + rtf x) 0 samp / fromIntegral (VU.length samp)
+    | otherwise = rtf (VU.sum samp) / fromIntegral (VU.length samp)
 {-# INLINE [0] mean' #-}
 
 meanDouble' :: VU.Vector Double -> Double
@@ -31,9 +29,7 @@ meanDouble' samp
 meanInt' :: VU.Vector Int -> Double
 meanInt' samp
     | VU.null samp = throw $ EmptyDataSetException "mean"
-    | otherwise =
-        VU.foldl' (\acc x -> acc + fromIntegral x) 0 samp
-            / fromIntegral (VU.length samp)
+    | otherwise = fromIntegral (VU.sum samp) / fromIntegral (VU.length samp)
 {-# INLINE meanInt' #-}
 
 {-# RULES
@@ -58,8 +54,7 @@ median' samp
             then pure (rtf middleElement)
             else do
                 prev <- VUM.read mutableSamp (middleIndex - 1)
-                -- Widen before adding: 'a' addition wraps for Int.
-                pure ((rtf middleElement + rtf prev) / 2)
+                pure (rtf (middleElement + prev) / 2)
 {-# INLINE median' #-}
 
 -- accumulator: count, mean, m2
@@ -121,29 +116,31 @@ skewness' :: (VU.Unbox a, Real a, Num a) => VU.Vector a -> Double
 skewness' = computeSkewness . VU.foldl' skewnessStep (SkewAcc 0 0 0 0)
 {-# INLINE skewness' #-}
 
-{- | Centered two-pass form: the one-pass @n*Sxy - Sx*Sy@ form cancels
-catastrophically on low-variance columns and can report |r| > 1.
--}
+data CorrelationStats
+    = CorrelationStats
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+        {-# UNPACK #-} !Double
+
 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
-            !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)))
+            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)
   where
     n = VU.length xs
-    -- 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
+    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)
 {-# INLINE correlation' #-}
 
 quantiles' ::
diff --git a/dataframe-operations/src/DataFrame/Operations/Core.hs b/dataframe-operations/src/DataFrame/Operations/Core.hs
index ef7e66c9..2bbeee73 100644
--- a/dataframe-operations/src/DataFrame/Operations/Core.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Core.hs
@@ -67,10 +67,9 @@ import DataFrame.Internal.Column (
     Column (..),
     Columnable,
     TypedColumn (..),
-    columnBitmap,
     columnLength,
     columnTypeString,
-    dropNulls,
+    dropNullsExceptMaybe,
     fromList,
     fromVector,
     materializeMerged,
@@ -690,7 +689,7 @@ valueCounts ::
     forall a. (Ord a, Columnable a) => Expr a -> DataFrame -> [(a, Int)]
 valueCounts expr df
     | null df = throw (EmptyDataSetException "valueCounts")
-    | otherwise = case columnAsVectorNonNull expr df of
+    | otherwise = case nonNullElements expr df of
         Left e -> throw e
         Right column' ->
             let
@@ -698,19 +697,21 @@ valueCounts expr df
              in
                 M.toAscList column
 
--- | As 'columnAsVector', minus null slots: a sentinel is not a category.
-columnAsVectorNonNull ::
+{- | The column's present values: 'columnAsVector' after 'dropNulls'. Shorter
+than the column when it has nulls; a sentinel is not a category.
+-}
+nonNullElements ::
     forall a.
     (Columnable a) => Expr a -> DataFrame -> Either DataFrameException (V.Vector a)
-columnAsVectorNonNull expr df = case expr of
+nonNullElements expr df = case expr of
     Col name -> case getColumn name df of
-        Just col -> withColumnName name (dropNulls (columnBitmap col) <$> toVector col)
+        Just col -> withColumnName name (toVector (dropNullsExceptMaybe @a col))
         Nothing ->
             Left $
                 ColumnsNotFoundException [name] "valueCounts" (M.keys $ columnIndices df)
     _ -> case interpret df expr of
         Left e -> throw e
-        Right (TColumn col) -> dropNulls (columnBitmap col) <$> toVector col
+        Right (TColumn col) -> toVector (dropNullsExceptMaybe @a col)
 
 {- | O (k * n) Shows the proportions of each value in a given column.
 
@@ -728,7 +729,7 @@ valueProportions ::
     forall a. (Ord a, Columnable a) => Expr a -> DataFrame -> [(a, Double)]
 valueProportions expr df
     | null df = throw (EmptyDataSetException "valueCounts")
-    | otherwise = case columnAsVectorNonNull expr df of
+    | otherwise = case nonNullElements expr df of
         Left e -> throw e
         Right column' ->
             let
diff --git a/dataframe-operations/src/DataFrame/Operations/Inference.hs b/dataframe-operations/src/DataFrame/Operations/Inference.hs
index aec40f5f..6464d4a2 100644
--- a/dataframe-operations/src/DataFrame/Operations/Inference.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Inference.hs
@@ -65,9 +65,15 @@ byteStringDateParser "%Y-%m-%d" = parseDateField
 byteStringDateParser fmt = readByteStringDate fmt
 {-# INLINE byteStringDateParser #-}
 
--- | Alias for 'readInt', which now rejects overflow itself.
+{- | '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.
+-}
 readIntStrict :: T.Text -> Maybe Int
-readIntStrict = readInt
+readIntStrict t
+    | T.length t <= 18 = readInt t
+    | otherwise = parseIntField (TE.encodeUtf8 t)
 {-# 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 6ce9475c..1eacc8c9 100644
--- a/dataframe-operations/src/DataFrame/Operations/Statistics.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
@@ -118,9 +118,9 @@ mean (Col name) df = case _getColumnAsDouble name df of
     Nothing -> error "[INTERNAL ERROR] Column is non-numeric"
 mean expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> mean' (dropNulls (columnBitmap col) xs)
+        Right xs -> mean' xs
 
 meanMaybe ::
     forall a. (Columnable a, Real a) => Expr (Maybe a) -> DataFrame -> Double
@@ -136,14 +136,17 @@ meanMaybe expr df = case interpret @(Maybe a) df expr of
 -- | Calculates the median of a given column as a standalone value.
 median ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
-median (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> median' (dropNulls (colBitmap name df) xs)
-    Left e -> throw e
+median (Col name) df = case getColumn name df of
+    Just col -> case toUnboxedVector @a (dropNulls col) of
+        Right xs -> median' xs
+        Left e -> throw e
+    Nothing ->
+        throw $ ColumnsNotFoundException [name] "median" (M.keys $ columnIndices df)
 median expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> median' (dropNulls (columnBitmap col) xs)
+        Right xs -> median' xs
 
 -- | Calculates the median of a given column (containing optional values) as a standalone value.
 medianMaybe ::
@@ -161,51 +164,65 @@ medianMaybe expr df = case interpret @(Maybe a) df expr of
 percentile ::
     forall a.
     (Columnable a, Real a, VU.Unbox a) => Int -> Expr a -> DataFrame -> Double
-percentile n (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> percentile' n (dropNulls (colBitmap name df) xs)
-    Left e -> throw e
+percentile n (Col name) df = case getColumn name df of
+    Just col -> case toUnboxedVector @a (dropNulls col) of
+        Right xs -> percentile' n xs
+        Left e -> throw e
+    Nothing ->
+        throw $ ColumnsNotFoundException [name] "percentile" (M.keys $ columnIndices df)
 percentile n expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> percentile' n (dropNulls (columnBitmap col) xs)
+        Right xs -> percentile' n xs
 
 -- | Calculates the nth percentile of a given column as a standalone value.
 genericPercentile ::
     forall a.
     (Columnable a, Ord a) => Int -> Expr a -> DataFrame -> a
-genericPercentile n (Col name) df = case columnAsVector (Col @a name) df of
-    Right xs -> percentileOrd' n (dropNulls (colBitmap name df) xs)
-    Left e -> throw e
+genericPercentile n (Col name) df = case getColumn name df of
+    Just col -> case toVector @a (dropNullsExceptMaybe @a col) of
+        Right xs -> percentileOrd' n xs
+        Left e -> throw e
+    Nothing ->
+        throw $
+            ColumnsNotFoundException [name] "genericPercentile" (M.keys $ columnIndices df)
 genericPercentile n expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toVector @a col of
+    Right (TColumn col) -> case toVector @a (dropNullsExceptMaybe @a col) of
         Left e -> throw e
-        Right xs -> percentileOrd' n (dropNulls (columnBitmap col) xs)
+        Right xs -> percentileOrd' n xs
 
 -- | Calculates the standard deviation of a given column as a standalone value.
 standardDeviation ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
-standardDeviation (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> (sqrt . variance') (dropNulls (colBitmap name df) xs)
-    Left e -> throw e
+standardDeviation (Col name) df = case getColumn name df of
+    Just col -> case toUnboxedVector @a (dropNulls col) of
+        Right xs -> (sqrt . variance') xs
+        Left e -> throw e
+    Nothing ->
+        throw $
+            ColumnsNotFoundException [name] "standardDeviation" (M.keys $ columnIndices df)
 standardDeviation expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> (sqrt . variance') (dropNulls (columnBitmap col) xs)
+        Right xs -> (sqrt . variance') xs
 
 -- | Calculates the skewness of a given column as a standalone value.
 skewness ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
-skewness (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> skewness' (dropNulls (colBitmap name df) xs)
-    Left e -> throw e
+skewness (Col name) df = case getColumn name df of
+    Just col -> case toUnboxedVector @a (dropNulls col) of
+        Right xs -> skewness' xs
+        Left e -> throw e
+    Nothing ->
+        throw $ ColumnsNotFoundException [name] "skewness" (M.keys $ columnIndices df)
 skewness expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> skewness' (dropNulls (columnBitmap col) xs)
+        Right xs -> skewness' xs
 
 -- | Calculates the variance of a given column as a standalone value.
 variance ::
@@ -215,21 +232,25 @@ variance (Col name) df = case _getColumnAsDouble name df of
     Nothing -> error "[INTERNAL ERROR] Column is non-numeric"
 variance expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> variance' (dropNulls (columnBitmap col) xs)
+        Right xs -> variance' xs
 
 -- | Calculates the inter-quartile range of a given column as a standalone value.
 interQuartileRange ::
     forall a. (Columnable a, Real a, VU.Unbox a) => Expr a -> DataFrame -> Double
-interQuartileRange (Col name) df = case columnAsUnboxedVector (Col @a name) df of
-    Right xs -> interQuartileRange' (dropNulls (colBitmap name df) xs)
-    Left e -> throw e
+interQuartileRange (Col name) df = case getColumn name df of
+    Just col -> case toUnboxedVector @a (dropNulls col) of
+        Right xs -> interQuartileRange' xs
+        Left e -> throw e
+    Nothing ->
+        throw $
+            ColumnsNotFoundException [name] "interQuartileRange" (M.keys $ columnIndices df)
 interQuartileRange expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn col) -> case toUnboxedVector @a col of
+    Right (TColumn col) -> case toUnboxedVector @a (dropNulls col) of
         Left e -> throw e
-        Right xs -> interQuartileRange' (dropNulls (columnBitmap col) xs)
+        Right xs -> interQuartileRange' xs
 
 -- | Calculates the Pearson's correlation coefficient between two given columns as a standalone value.
 correlation :: T.Text -> T.Text -> DataFrame -> Maybe Double
@@ -240,25 +261,21 @@ correlation first second df = do
     s <- _getColumnAsDouble second df'
     correlation' f s
 
--- | Bitmap of a named column, if it has one.
-colBitmap :: T.Text -> DataFrame -> Maybe Bitmap
-colBitmap name df = columnBitmap =<< getColumn name df
-
 _getColumnAsDouble :: T.Text -> DataFrame -> Maybe (VU.Vector Double)
 _getColumnAsDouble name df = case getColumn name df of
-    Just (UnboxedColumn bm (f' :: VU.Vector a)) ->
-        let f = dropNulls bm f'
-         in case testEquality (typeRep @a) (typeRep @Double) of
+    Just col -> case dropNulls col of
+        UnboxedColumn _ (f :: VU.Vector a) ->
+            case testEquality (typeRep @a) (typeRep @Double) of
                 Just Refl -> Just f
                 Nothing -> case sIntegral @a of
                     STrue -> Just (VU.map fromIntegral f)
                     SFalse -> case sFloating @a of
                         STrue -> Just (VU.map realToFrac f)
                         SFalse -> Nothing
+        _ -> Nothing
     Nothing ->
         throw $
             ColumnsNotFoundException [name] "_getColumnAsDouble" (M.keys $ columnIndices df)
-    _ -> Nothing
 {-# INLINE _getColumnAsDouble #-}
 
 optionalToDoubleVector :: (Real a) => V.Vector (Maybe a) -> VU.Vector Double
@@ -273,19 +290,20 @@ sum ::
     forall a. (Columnable a, Num a) => Expr a -> DataFrame -> a
 sum (Col name) df = case getColumn name df of
     Nothing -> throw $ ColumnsNotFoundException [name] "sum" (M.keys $ columnIndices df)
-    Just ((UnboxedColumn bm (column :: VU.Vector a'))) -> case testEquality (typeRep @a') (typeRep @a) of
-        Just Refl -> VG.sum (dropNulls bm column)
-        Nothing -> 0
-    Just ((BoxedColumn bm (column :: V.Vector a'))) -> case testEquality (typeRep @a') (typeRep @a) of
-        Just Refl -> VG.sum (dropNulls bm column)
-        Nothing -> 0
-    Just (PackedText _ _) -> 0
-    Just (MergedColumn _ _) -> 0 -- matches the old eager These column (type never Num)
+    Just c -> case dropNulls c of
+        UnboxedColumn _ (column :: VU.Vector a') -> case testEquality (typeRep @a') (typeRep @a) of
+            Just Refl -> VG.sum column
+            Nothing -> 0
+        BoxedColumn _ (column :: V.Vector a') -> case testEquality (typeRep @a') (typeRep @a) of
+            Just Refl -> VG.sum column
+            Nothing -> 0
+        PackedText _ _ -> 0
+        MergedColumn _ _ -> 0 -- matches the old eager These column (type never Num)
 sum expr df = case interpret df expr of
     Left e -> throw e
-    Right (TColumn xs) -> case toVector @a @V.Vector xs of
+    Right (TColumn xs) -> case toVector @a @V.Vector (dropNulls xs) of
         Left e -> throw e
-        Right xs' -> VG.sum (dropNulls (columnBitmap xs) xs')
+        Right xs' -> VG.sum xs'
 
 {- | /O(n)/ Impute missing values in a column using a derived scalar.
 
diff --git a/dataframe-operations/src/DataFrame/Operations/Typing.hs b/dataframe-operations/src/DataFrame/Operations/Typing.hs
index e314f296..4cc90f96 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. 'readInt' rejects overflow
+parsing as neither demotes the column to Text. 'readIntStrict' 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) readInt readDouble cols of
+    case promoteIntColumn (\_ t -> isNull t) readIntStrict 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 476cc287..6e8928a8 100644
--- a/dataframe-parsing/dataframe-parsing.cabal
+++ b/dataframe-parsing/dataframe-parsing.cabal
@@ -43,6 +43,7 @@ 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 b2d50558..a7126247 100644
--- a/dataframe-parsing/src-internal/DataFrame/Internal/Parsing.hs
+++ b/dataframe-parsing/src-internal/DataFrame/Internal/Parsing.hs
@@ -8,15 +8,13 @@ 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.Char (isDigit)
+import Data.ByteString.Lex.Fractional
 import Data.Foldable (fold)
-import Data.Text.Read (decimal, signed)
+import Data.Text.Read (decimal, double, signed)
 import Data.Time (Day, defaultTimeLocale, parseTimeM)
 import GHC.Stack (HasCallStack)
 import System.IO (Handle, IOMode (..), hIsEOF, hTell, withFile)
@@ -63,23 +61,11 @@ 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
-    | 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
+readInt s = case signed decimal (T.strip s) of
+    Left _ -> Nothing
+    Right (value, "") -> Just value
+    Right (_value, _) -> Nothing
 {-# INLINE readInt #-}
 
 readByteStringInt :: (HasCallStack) => C.ByteString -> Maybe Int
@@ -107,68 +93,23 @@ 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
-    | t == "Infinity" = Just (1 / 0)
-    | t == "-Infinity" = Just (-1 / 0)
-    | otherwise = parseDoubleExact t
-  where
-    t = C.strip s
+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
 {-# 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 = readByteStringDouble . TE.encodeUtf8
+readDouble s =
+    case signed double s of
+        Left _ -> Nothing
+        Right (value, "") -> Just value
+        Right (_value, _) -> Nothing
 {-# INLINE readDouble #-}
 
 readIntegerEither :: (HasCallStack) => T.Text -> Either T.Text Integer
@@ -178,24 +119,19 @@ 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
-    | 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
+readIntEither s = case signed decimal (T.strip s) of
+    Left _ -> Left s
+    Right (value, "") -> Right value
+    Right (_value, _) -> Left s
 {-# INLINE readIntEither #-}
 
 readDoubleEither :: (HasCallStack) => T.Text -> Either T.Text Double
-readDoubleEither s = maybe (Left s) Right (readDouble s)
+readDoubleEither s =
+    case signed double s of
+        Left _ -> Left s
+        Right (value, "") -> Right value
+        Right (_value, _) -> Left 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 1fe3cde9..d5dce287 100644
--- a/dataframe-parsing/src-internal/DataFrame/Internal/Parsing/Fast/Double.hs
+++ b/dataframe-parsing/src-internal/DataFrame/Internal/Parsing/Fast/Double.hs
@@ -3,9 +3,8 @@
 {-# LANGUAGE UnboxedTuples #-}
 
 {- | Fast @Double@ slice parser, bit-exact with @readByteStringDouble@.
-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.
+Replays the reference parser's exact floating-point operations via 'Word64'
+digit accumulation and 10^k tables, falling back when exactness is in doubt.
 -}
 module DataFrame.Internal.Parsing.Fast.Double (parseDoubleField#) where
 
@@ -19,31 +18,33 @@ import GHC.Exts (Double (..), Double#, Int#)
 import DataFrame.Internal.Parsing (readByteStringDouble)
 import DataFrame.Internal.Parsing.Fast.Common
 
-{- | @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.
+{- | @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.
 -}
 pow10Table :: VU.Vector Double
-pow10Table = VU.generate (exactPow10 + 1) (10 ^)
+pow10Table = VU.generate (tableMax + 1) (10 ^)
 {-# NOINLINE pow10Table #-}
 
--- | @10 ^ k@ for @k <= 19@ in 'Word64'; @10^19 < 2^64@.
-pow10w :: VU.Vector Word64
-pow10w = VU.generate 20 (10 ^)
-{-# NOINLINE pow10w #-}
+-- | @recip (10 ^ k)@, replaying @10 ^^ negate k@ bit-for-bit.
+recipPow10Table :: VU.Vector Double
+recipPow10Table = VU.map recip pow10Table
+{-# NOINLINE recipPow10Table #-}
 
--- | Largest @k@ with @10^k@ exactly representable as a 'Double'.
-exactPow10 :: Int
-exactPow10 = 22
+tableMax :: Int
+tableMax = 1024
 
-{- | 'Word64' to 'Double'; callers only pass values @<= 2^53@, where the
-conversion is exact.
+{- | '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.
 -}
 w2d :: Word64 -> Double
-w2d w = fromIntegral (fromIntegral w :: Int)
+w2d w
+    | w <= 9007199254740991 = fromIntegral (fromIntegral w :: Int)
+    | otherwise = fromInteger (toInteger w)
 {-# INLINE w2d #-}
 
--- | Outside the exact window: hand the raw slice to the reference parser.
+-- | Exactness in doubt: 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
@@ -57,9 +58,6 @@ 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
@@ -71,27 +69,14 @@ 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 - iz) wEnd
+                                else afterWhole neg w wEnd
   where
     !i0 = skipStrip bs start end0
     !end = skipStripEnd bs i0 end0
 
     none = (# 0#, 0.0## #)
 
-    -- "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
+    afterWhole !neg !w !i
         | i < end && BSU.unsafeIndex bs i == 0x2E =
             let !f0 = i + 1
                 !fz = skipZeroes bs f0 end
@@ -99,19 +84,13 @@ parseDoubleField# bs start end0
                     if fEnd == f0
                         then none
                         else
-                            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)
+                            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
         | BSU.unsafeIndex bs i == 0x65 || BSU.unsafeIndex bs i == 0x45 =
             let !i1 = i + 1
                 !eneg = i1 < end && BSU.unsafeIndex bs i1 == 0x2D
@@ -126,29 +105,18 @@ parseDoubleField# bs start end0
                                     else
                                         if eEnd - ez > 18
                                             then referenceSlice bs start end0
-                                            else
-                                                let !ev = fromIntegral e
-                                                 in finish
-                                                        neg
-                                                        sig
-                                                        ( (if eneg then negate ev else ev)
-                                                            - fracLen
-                                                        )
+                                            else done neg (val * scale eneg (fromIntegral e))
         | otherwise = none
 
-    -- 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 #-}
+    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 #-}
 
     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 7d207bd2..eb75dff2 100644
--- a/dataframe-parsing/tests/Properties/FastParsing.hs
+++ b/dataframe-parsing/tests/Properties/FastParsing.hs
@@ -185,12 +185,10 @@ 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 =
-    not (isNaN d) ==>
-        let f = C.pack (show d)
-         in bitsOf (parseDoubleField f) === bitsOf (Just d)
+    let f = C.pack (show d)
+     in bitsOf (parseDoubleField f) === bitsOf (readByteStringDouble f)
 
 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 0ce3771a..b09c9676 100644
--- a/dataframe-parsing/tests/Unit/FastParsing.hs
+++ b/dataframe-parsing/tests/Unit/FastParsing.hs
@@ -4,13 +4,11 @@
 
 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. Doubles are additionally pinned
-against 'read', which is correctly rounded.
+pins a concrete value, the value itself.
 -}
 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)
@@ -96,9 +94,6 @@ doubleCases =
     , "1e-1024"
     , "1e-1025"
     , "Infinity"
-    , "-Infinity"
-    , "+Infinity"
-    , " Infinity "
     , "NaN"
     , "1e"
     , "1e+"
@@ -119,28 +114,11 @@ pinnedDoubles =
     , ("5.", Nothing)
     , ("1e3.5", Nothing)
     , ("1_000", Nothing)
-    , ("Infinity", Just (1 / 0))
-    , ("-Infinity", Just (-1 / 0))
+    , ("Infinity", Nothing)
     , ("NaN", Nothing)
-    , ("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"
+    , -- 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)
     ]
 
 boolCases :: [BS.ByteString]
@@ -219,14 +197,6 @@ 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
@@ -234,7 +204,6 @@ 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/dataframe.cabal b/dataframe.cabal
index b7310c6f..5089bdf9 100644
--- a/dataframe.cabal
+++ b/dataframe.cabal
@@ -319,7 +319,6 @@ test-suite tests
                    Internal.DictEncode,
                    Internal.Markdown,
                    Internal.PackedText,
-                   Internal.Parsing,
                    PackedTextMigration,
                    PrettyPrint,
                    Learn.Denotation,
diff --git a/docs/exploratory_data_analysis_primer.md b/docs/exploratory_data_analysis_primer.md
index f6c5b4eb..f97944eb 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.9776922140978362
+> 0.977668529406543
 
 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
deleted file mode 100644
index 512e31f1..00000000
--- a/tests/Internal/Parsing.hs
+++ /dev/null
@@ -1,335 +0,0 @@
-{-# LANGUAGE OverloadedStrings #-}
-
-module Internal.Parsing where
-
-import DataFrame.Internal.Parsing
-import Test.HUnit
-
--- isNullish: recognized null strings
-
-isNullishEmptyString :: Test
-isNullishEmptyString =
-    TestCase (assertBool "empty string is nullish" (isNullish ""))
-
-isNullishNA :: Test
-isNullishNA = TestCase (assertBool "NA is nullish" (isNullish "NA"))
-
-isNullishNULL :: Test
-isNullishNULL = TestCase (assertBool "NULL is nullish" (isNullish "NULL"))
-
-isNullishNull :: Test
-isNullishNull = TestCase (assertBool "null is nullish" (isNullish "null"))
-
-isNullishNaN :: Test
-isNullishNaN = TestCase (assertBool "nan is nullish" (isNullish "nan"))
-
-isNullishNaNMixed :: Test
-isNullishNaNMixed = TestCase (assertBool "NaN is nullish" (isNullish "NaN"))
-
-isNullishNANUpper :: Test
-isNullishNANUpper = TestCase (assertBool "NAN is nullish" (isNullish "NAN"))
-
-isNullishNothing :: Test
-isNullishNothing =
-    TestCase (assertBool "Nothing is nullish" (isNullish "Nothing"))
-
-isNullishSpace :: Test
-isNullishSpace =
-    TestCase (assertBool "single space is nullish" (isNullish " "))
-
-isNullishNSlashA :: Test
-isNullishNSlashA = TestCase (assertBool "N/A is nullish" (isNullish "N/A"))
-
--- isNullish: values that are NOT null
-
-notNullishNumber :: Test
-notNullishNumber =
-    TestCase (assertBool "\"42\" is not nullish" (not (isNullish "42")))
-
-notNullishText :: Test
-notNullishText =
-    TestCase (assertBool "\"hello\" is not nullish" (not (isNullish "hello")))
-
-notNullishTrue :: Test
-notNullishTrue =
-    TestCase (assertBool "\"True\" is not nullish" (not (isNullish "True")))
-
-notNullishDouble :: Test
-notNullishDouble =
-    TestCase (assertBool "\"3.14\" is not nullish" (not (isNullish "3.14")))
-
-notNullishZero :: Test
-notNullishZero =
-    TestCase (assertBool "\"0\" is not nullish" (not (isNullish "0")))
-
--- readBool: positive cases
-
-readBoolTrue :: Test
-readBoolTrue =
-    TestCase (assertEqual "readBool \"True\"" (Just True) (readBool "True"))
-
-readBoolTrueLower :: Test
-readBoolTrueLower =
-    TestCase (assertEqual "readBool \"true\"" (Just True) (readBool "true"))
-
-readBoolTrueUpper :: Test
-readBoolTrueUpper =
-    TestCase (assertEqual "readBool \"TRUE\"" (Just True) (readBool "TRUE"))
-
-readBoolFalse :: Test
-readBoolFalse =
-    TestCase (assertEqual "readBool \"False\"" (Just False) (readBool "False"))
-
-readBoolFalseLower :: Test
-readBoolFalseLower =
-    TestCase (assertEqual "readBool \"false\"" (Just False) (readBool "false"))
-
-readBoolFalseUpper :: Test
-readBoolFalseUpper =
-    TestCase (assertEqual "readBool \"FALSE\"" (Just False) (readBool "FALSE"))
-
--- readBool: values that are not booleans
-
-readBoolDigit :: Test
-readBoolDigit =
-    TestCase (assertEqual "readBool \"1\" is Nothing" Nothing (readBool "1"))
-
-readBoolYes :: Test
-readBoolYes =
-    TestCase (assertEqual "readBool \"yes\" is Nothing" Nothing (readBool "yes"))
-
-readBoolEmpty :: Test
-readBoolEmpty =
-    TestCase (assertEqual "readBool \"\" is Nothing" Nothing (readBool ""))
-
-readBoolPartialTrue :: Test
-readBoolPartialTrue =
-    TestCase (assertEqual "readBool \"Tru\" is Nothing" Nothing (readBool "Tru"))
-
--- readInt
-
-readIntPositive :: Test
-readIntPositive =
-    TestCase (assertEqual "readInt \"42\"" (Just 42) (readInt "42"))
-
-readIntNegative :: Test
-readIntNegative =
-    TestCase (assertEqual "readInt \"-17\"" (Just (-17)) (readInt "-17"))
-
-readIntZero :: Test
-readIntZero =
-    TestCase (assertEqual "readInt \"0\"" (Just 0) (readInt "0"))
-
--- readInt strips whitespace before parsing
-readIntLeadingSpace :: Test
-readIntLeadingSpace =
-    TestCase
-        ( assertEqual
-            "readInt \" 5 \" (strips whitespace)"
-            (Just 5)
-            (readInt " 5 ")
-        )
-
-readIntFloat :: Test
-readIntFloat =
-    TestCase
-        (assertEqual "readInt \"3.14\" is Nothing" Nothing (readInt "3.14"))
-
-readIntText :: Test
-readIntText =
-    TestCase (assertEqual "readInt \"abc\" is Nothing" Nothing (readInt "abc"))
-
-readIntEmpty :: Test
-readIntEmpty =
-    TestCase (assertEqual "readInt \"\" is Nothing" Nothing (readInt ""))
-
--- trailing non-digits must make the parse fail
-readIntPartialSuffix :: Test
-readIntPartialSuffix =
-    TestCase
-        ( assertEqual
-            "readInt \"42abc\" is Nothing"
-            Nothing
-            (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
-readDoublePositive =
-    TestCase
-        (assertEqual "readDouble \"3.14\"" (Just 3.14) (readDouble "3.14"))
-
-readDoubleNegative :: Test
-readDoubleNegative =
-    TestCase
-        (assertEqual "readDouble \"-1.5\"" (Just (-1.5)) (readDouble "-1.5"))
-
-readDoubleWholeNumber :: Test
-readDoubleWholeNumber =
-    TestCase
-        ( assertEqual
-            "readDouble \"42\" parses as 42.0"
-            (Just 42.0)
-            (readDouble "42")
-        )
-
-readDoubleText :: Test
-readDoubleText =
-    TestCase
-        (assertEqual "readDouble \"abc\" is Nothing" Nothing (readDouble "abc"))
-
-readDoubleEmpty :: Test
-readDoubleEmpty =
-    TestCase
-        (assertEqual "readDouble \"\" is Nothing" Nothing (readDouble ""))
-
-readDoublePartialSuffix :: Test
-readDoublePartialSuffix =
-    TestCase
-        ( assertEqual
-            "readDouble \"3.14abc\" is Nothing"
-            Nothing
-            (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
-    , TestLabel "isNullishNA" isNullishNA
-    , TestLabel "isNullishNULL" isNullishNULL
-    , TestLabel "isNullishNull" isNullishNull
-    , TestLabel "isNullishNaN" isNullishNaN
-    , TestLabel "isNullishNaNMixed" isNullishNaNMixed
-    , TestLabel "isNullishNANUpper" isNullishNANUpper
-    , TestLabel "isNullishNothing" isNullishNothing
-    , TestLabel "isNullishSpace" isNullishSpace
-    , TestLabel "isNullishNSlashA" isNullishNSlashA
-    , TestLabel "notNullishNumber" notNullishNumber
-    , TestLabel "notNullishText" notNullishText
-    , TestLabel "notNullishTrue" notNullishTrue
-    , TestLabel "notNullishDouble" notNullishDouble
-    , TestLabel "notNullishZero" notNullishZero
-    , TestLabel "readBoolTrue" readBoolTrue
-    , TestLabel "readBoolTrueLower" readBoolTrueLower
-    , TestLabel "readBoolTrueUpper" readBoolTrueUpper
-    , TestLabel "readBoolFalse" readBoolFalse
-    , TestLabel "readBoolFalseLower" readBoolFalseLower
-    , TestLabel "readBoolFalseUpper" readBoolFalseUpper
-    , TestLabel "readBoolDigit" readBoolDigit
-    , TestLabel "readBoolYes" readBoolYes
-    , TestLabel "readBoolEmpty" readBoolEmpty
-    , TestLabel "readBoolPartialTrue" readBoolPartialTrue
-    , TestLabel "readIntPositive" readIntPositive
-    , TestLabel "readIntNegative" readIntNegative
-    , TestLabel "readIntZero" readIntZero
-    , TestLabel "readIntLeadingSpace" readIntLeadingSpace
-    , TestLabel "readIntFloat" readIntFloat
-    , 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 2c018a17..fec81802 100644
--- a/tests/Learn/MetricsTests.hs
+++ b/tests/Learn/MetricsTests.hs
@@ -1,8 +1,11 @@
 {-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 
 module Learn.MetricsTests (tests) where
 
+import qualified Control.Exception as E
+
 import qualified DataFrame as D
 import qualified DataFrame.Functions as F
 import qualified DataFrame.Internal.Column as DI
@@ -50,9 +53,11 @@ testRegressionMetrics = TestCase $ do
     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])))
+    -- No compared pairs throws instead of scoring.
+    r <- E.try (E.evaluate (mse VU.empty (VU.fromList [5, 5, 5])))
+    case r of
+        Left (_ :: E.SomeException) -> pure ()
+        Right v -> assertFailure ("mse with no pairs returned " ++ show v)
 
 testMulticlassMetrics :: Test
 testMulticlassMetrics = TestCase $ do
diff --git a/tests/Main.hs b/tests/Main.hs
index 7689e2f7..3263a53d 100644
--- a/tests/Main.hs
+++ b/tests/Main.hs
@@ -17,7 +17,6 @@ import qualified Internal.ColumnBuilder
 import qualified Internal.DictEncode
 import qualified Internal.Markdown
 import qualified Internal.PackedText
-import qualified Internal.Parsing
 import qualified LazyParity
 import qualified LazyParquet
 import qualified LazyProjection
@@ -75,7 +74,6 @@ tests =
             ++ Internal.DictEncode.tests
             ++ Internal.Markdown.tests
             ++ Internal.PackedText.tests
-            ++ Internal.Parsing.tests
             ++ Learn.Denotation.tests
             ++ Learn.Models.tests
             ++ Learn.TypedModel.tests
diff --git a/tests/Operations/Statistics.hs b/tests/Operations/Statistics.hs
index c27cad86..636d6e94 100644
--- a/tests/Operations/Statistics.hs
+++ b/tests/Operations/Statistics.hs
@@ -183,57 +183,6 @@ 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
@@ -491,10 +440,6 @@ 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

From bd873ccf6c131f6d55a1cad3e449b7cbfd9807ef Mon Sep 17 00:00:00 2001
From: skymanbp <57272723+skymanbp@users.noreply.github.com>
Date: Thu, 20 Aug 2026 00:03:14 -0400
Subject: [PATCH 4/5] Shorten or drop inline comments

---
 .../src-internal/DataFrame/Internal/AggKernel.hs         | 3 +--
 .../src-internal/DataFrame/Internal/AggKernelPar.hs      | 3 +--
 dataframe-learn/src/DataFrame/Metrics.hs                 | 1 -
 .../src-internal/DataFrame/Internal/Statistics.hs        | 4 ++--
 .../src/DataFrame/Operations/Statistics.hs               | 9 +++++----
 tests/Learn/MetricsTests.hs                              | 2 --
 tests/Operations/ParallelGroupBy.hs                      | 3 +--
 tests/Operations/Statistics.hs                           | 5 ++---
 tests/Operations/VectorKernel.hs                         | 3 +--
 9 files changed, 13 insertions(+), 20 deletions(-)

diff --git a/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs b/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs
index 2b773c06..a01ef9f0 100644
--- a/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs
+++ b/dataframe-core/src-internal/DataFrame/Internal/AggKernel.hs
@@ -217,8 +217,7 @@ varScatter takeSqrt g nGroups v = runST $ do
             | otherwise = do
                 c <- VUM.unsafeRead cnt k
                 mm <- VUM.unsafeRead m2 k
-                -- Sample variance is undefined at n = 1: NaN, matching
-                -- 'computeVariance'.
+                -- NaN at n = 1, 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)
diff --git a/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs b/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs
index 05308f7b..cd10c924 100644
--- a/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs
+++ b/dataframe-core/src-internal/DataFrame/Internal/AggKernelPar.hs
@@ -290,8 +290,7 @@ varPar takeSqrt vis offs nGroups v caps bounds = do
             | otherwise = do
                 c <- VUM.unsafeRead cnt k
                 mm <- VUM.unsafeRead m2 k
-                -- Sample variance is undefined at n = 1: NaN, matching
-                -- 'computeVariance'.
+                -- NaN at n = 1, 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)
diff --git a/dataframe-learn/src/DataFrame/Metrics.hs b/dataframe-learn/src/DataFrame/Metrics.hs
index eb60725e..c23a6d21 100644
--- a/dataframe-learn/src/DataFrame/Metrics.hs
+++ b/dataframe-learn/src/DataFrame/Metrics.hs
@@ -84,7 +84,6 @@ nCompared preds truth = fromIntegral (min (VU.length preds) (VU.length truth))
 -- | Mean squared error.
 mse :: Metric
 mse preds truth
-    -- No compared pairs is not a score of any kind.
     | n == 0 = throw (EmptyDataSetException "mse")
     | otherwise =
         VU.sum (VU.zipWith (\p t -> (p - t) ^ (2 :: Int)) preds truth) / n
diff --git a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
index d7403bf9..e79e8e84 100644
--- a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
+++ b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
@@ -74,7 +74,7 @@ varianceStep (VarAcc !n !meanVal !m2) !x =
 computeVariance :: VarAcc -> Double
 computeVariance (VarAcc !n _ !m2)
     | n == 0 = throw $ EmptyDataSetException "variance"
-    -- Sample variance is undefined at n = 1: NaN, not a spurious 0.
+    -- undefined at n = 1: NaN, not 0
     | n < 2 = 0 / 0
     | otherwise = m2 / fromIntegral (n - 1)
 {-# INLINE computeVariance #-}
@@ -108,7 +108,7 @@ skewnessStep (SkewAcc !n !meanVal !m2 !m3) !x' =
 computeSkewness :: SkewAcc -> Double
 computeSkewness (SkewAcc n _ m2 m3)
     | n < 3 = 0 -- or error "skewness of <3 samples"
-    -- m2, m3 are raw sums, so population g1 = sqrt n * m3 / m2^(3/2).
+    -- raw sums: g1 = sqrt n * m3 / m2^(3/2)
     | otherwise = (sqrt (fromIntegral n) * m3) / sqrt (m2 ^ (3 :: Int))
 {-# INLINE computeSkewness #-}
 
diff --git a/dataframe-operations/src/DataFrame/Operations/Statistics.hs b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
index 1eacc8c9..1526d500 100644
--- a/dataframe-operations/src/DataFrame/Operations/Statistics.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Statistics.hs
@@ -252,10 +252,11 @@ interQuartileRange expr df = case interpret df expr of
         Left e -> throw e
         Right xs -> interQuartileRange' xs
 
--- | Calculates the Pearson's correlation coefficient between two given columns as a standalone value.
+{- | Calculates the Pearson's correlation coefficient between two given columns as a standalone value.
+Pairs with a null in either column are dropped.
+-}
 correlation :: T.Text -> T.Text -> DataFrame -> Maybe Double
 correlation first second df = do
-    -- Listwise deletion: a null in either column drops the pair.
     let df' = filterJust first (filterJust second df)
     f <- _getColumnAsDouble first df'
     s <- _getColumnAsDouble second df'
@@ -298,7 +299,7 @@ sum (Col name) df = case getColumn name df of
             Just Refl -> VG.sum column
             Nothing -> 0
         PackedText _ _ -> 0
-        MergedColumn _ _ -> 0 -- matches the old eager These column (type never Num)
+        MergedColumn _ _ -> 0 -- never numeric
 sum expr df = case interpret df expr of
     Left e -> throw e
     Right (TColumn xs) -> case toVector @a @V.Vector (dropNulls xs) of
@@ -443,7 +444,7 @@ summarize df =
 -- | Round a @Double@ to Specified Precision
 roundTo :: Int -> Double -> Double
 roundTo n x
-    -- 'round' on NaN yields garbage; keep NaN visible in summaries.
+    -- round on NaN is garbage
     | isNaN x = x
     | otherwise = fromInteger (round $ x * 10 ^ n) / 10.0 ^^ n
 
diff --git a/tests/Learn/MetricsTests.hs b/tests/Learn/MetricsTests.hs
index fec81802..87089b05 100644
--- a/tests/Learn/MetricsTests.hs
+++ b/tests/Learn/MetricsTests.hs
@@ -46,14 +46,12 @@ 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)
-    -- No compared pairs throws instead of scoring.
     r <- E.try (E.evaluate (mse VU.empty (VU.fromList [5, 5, 5])))
     case r of
         Left (_ :: E.SomeException) -> pure ()
diff --git a/tests/Operations/ParallelGroupBy.hs b/tests/Operations/ParallelGroupBy.hs
index 5cff0b02..c808387b 100644
--- a/tests/Operations/ParallelGroupBy.hs
+++ b/tests/Operations/ParallelGroupBy.hs
@@ -109,8 +109,7 @@ aggParityFor n =
                 ]
             seqDf = D.aggregate aggs (groupBySeq ["ki", "kt"] df)
             parDf = D.aggregate aggs (groupByPar ["ki", "kt"] df)
-         in -- Rendered comparison: singleton-group stddev is NaN, which the
-            -- Eq instance treats as unequal.
+         in -- render: NaN /= NaN under Eq
             assertEqual
                 ("aggregate parity n=" ++ show n)
                 (D.toMarkdown seqDf)
diff --git a/tests/Operations/Statistics.hs b/tests/Operations/Statistics.hs
index 636d6e94..e2185f5a 100644
--- a/tests/Operations/Statistics.hs
+++ b/tests/Operations/Statistics.hs
@@ -57,7 +57,7 @@ skewnessOfSymmetricDataSet =
             0
         )
 
--- Population skewness g1, the form the docs define (matches scipy.stats.skew).
+-- g1, matching scipy.stats.skew
 skewnessOfSimpleDataSet :: Test
 skewnessOfSimpleDataSet =
     TestCase
@@ -256,8 +256,7 @@ correlationMissingColumn =
             (print $ D.correlation "x" "missingcol" correlationDf)
         )
 
--- Nullable columns: statistics must skip null slots, not read the sentinel
--- stored there.
+-- stats must skip null slots
 
 nullableDf :: D.DataFrame
 nullableDf =
diff --git a/tests/Operations/VectorKernel.hs b/tests/Operations/VectorKernel.hs
index cc8f5247..802e015d 100644
--- a/tests/Operations/VectorKernel.hs
+++ b/tests/Operations/VectorKernel.hs
@@ -123,8 +123,7 @@ 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 -- Full render: value-exact via 'show', and NaN cells (singleton
-            -- group stddev/variance) compare equal, which 'Eq' refuses.
+         in -- render: NaN /= NaN under Eq
             assertEqual
                 ("kernel==interpreter " ++ label)
                 (D.toMarkdown ref)

From 321b3c6af970d964f3a10fa78de2b2bb388f7e40 Mon Sep 17 00:00:00 2001
From: skymanbp <57272723+skymanbp@users.noreply.github.com>
Date: Thu, 20 Aug 2026 03:33:27 -0400
Subject: [PATCH 5/5] Remove nonNullElements; nullable columns read as Maybe

---
 .../DataFrame/Internal/Statistics.hs          |  1 -
 .../src/DataFrame/Operations/Core.hs          | 21 ++-----------------
 tests/Operations/Statistics.hs                | 12 +++++------
 3 files changed, 8 insertions(+), 26 deletions(-)

diff --git a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
index e79e8e84..ef929d83 100644
--- a/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
+++ b/dataframe-operations/src-internal/DataFrame/Internal/Statistics.hs
@@ -108,7 +108,6 @@ skewnessStep (SkewAcc !n !meanVal !m2 !m3) !x' =
 computeSkewness :: SkewAcc -> Double
 computeSkewness (SkewAcc n _ m2 m3)
     | n < 3 = 0 -- or error "skewness of <3 samples"
-    -- raw sums: g1 = sqrt n * m3 / m2^(3/2)
     | otherwise = (sqrt (fromIntegral n) * m3) / sqrt (m2 ^ (3 :: Int))
 {-# INLINE computeSkewness #-}
 
diff --git a/dataframe-operations/src/DataFrame/Operations/Core.hs b/dataframe-operations/src/DataFrame/Operations/Core.hs
index 2bbeee73..77df934f 100644
--- a/dataframe-operations/src/DataFrame/Operations/Core.hs
+++ b/dataframe-operations/src/DataFrame/Operations/Core.hs
@@ -69,7 +69,6 @@ import DataFrame.Internal.Column (
     TypedColumn (..),
     columnLength,
     columnTypeString,
-    dropNullsExceptMaybe,
     fromList,
     fromVector,
     materializeMerged,
@@ -689,7 +688,7 @@ valueCounts ::
     forall a. (Ord a, Columnable a) => Expr a -> DataFrame -> [(a, Int)]
 valueCounts expr df
     | null df = throw (EmptyDataSetException "valueCounts")
-    | otherwise = case nonNullElements expr df of
+    | otherwise = case columnAsVector expr df of
         Left e -> throw e
         Right column' ->
             let
@@ -697,22 +696,6 @@ valueCounts expr df
              in
                 M.toAscList column
 
-{- | The column's present values: 'columnAsVector' after 'dropNulls'. Shorter
-than the column when it has nulls; a sentinel is not a category.
--}
-nonNullElements ::
-    forall a.
-    (Columnable a) => Expr a -> DataFrame -> Either DataFrameException (V.Vector a)
-nonNullElements expr df = case expr of
-    Col name -> case getColumn name df of
-        Just col -> withColumnName name (toVector (dropNullsExceptMaybe @a col))
-        Nothing ->
-            Left $
-                ColumnsNotFoundException [name] "valueCounts" (M.keys $ columnIndices df)
-    _ -> case interpret df expr of
-        Left e -> throw e
-        Right (TColumn col) -> toVector (dropNullsExceptMaybe @a col)
-
 {- | O (k * n) Shows the proportions of each value in a given column.
 
 ==== __Example__
@@ -729,7 +712,7 @@ valueProportions ::
     forall a. (Ord a, Columnable a) => Expr a -> DataFrame -> [(a, Double)]
 valueProportions expr df
     | null df = throw (EmptyDataSetException "valueCounts")
-    | otherwise = case nonNullElements expr df of
+    | otherwise = case columnAsVector expr df of
         Left e -> throw e
         Right column' ->
             let
diff --git a/tests/Operations/Statistics.hs b/tests/Operations/Statistics.hs
index e2185f5a..ae253d8f 100644
--- a/tests/Operations/Statistics.hs
+++ b/tests/Operations/Statistics.hs
@@ -408,13 +408,13 @@ correlationIgnoresNullRows =
                         (abs (r - 1.0) < 1e-10)
         )
 
-frequenciesSkipsNulls :: Test
-frequenciesSkipsNulls =
+frequenciesNullableAsMaybe :: Test
+frequenciesNullableAsMaybe =
     TestCase
         ( assertEqual
-            "frequencies has no sentinel category"
-            3 -- Statistic, 10 and 20
-            (D.nColumns (D.frequencies (F.col @Int "n") nullableIntDf))
+            "nulls are a Nothing category"
+            4 -- Statistic, Nothing, Just 10, Just 20
+            (D.nColumns (D.frequencies (F.col @(Maybe Int) "n") nullableIntDf))
         )
 
 tests :: [Test]
@@ -463,5 +463,5 @@ tests =
     , TestLabel "sumUnboxedNullable" sumUnboxedNullable
     , TestLabel "sumBoxedNullableDoesNotThrow" sumBoxedNullableDoesNotThrow
     , TestLabel "correlationIgnoresNullRows" correlationIgnoresNullRows
-    , TestLabel "frequenciesSkipsNulls" frequenciesSkipsNulls
+    , TestLabel "frequenciesNullableAsMaybe" frequenciesNullableAsMaybe
     ]