BUG: datetimelike arithmetic landing on the NaT sentinel returns NaT, plus float64 precision loss in timedelta64 arithmetic/reductions

[jbrockmendel]

Follow-up to GH-66510 / GH-66520, which fixed the parsing and construction half of the "value lands on the NaT sentinel" problem (parsers.pyx, conversion.pyx, strptime.pyx, Timestamp.replace). The arithmetic half is still open, and while investigating it I found an adjacent family of float64 precision bugs in timedelta64 arithmetic and reductions.

All reproducers below are against 082e77392bd (main as of 2026-07-29).

Two distinct root causes, filed together because they were found together and overlap in the code they touch — happy to split into separate issues if that's easier to triage.

---

1. NPY_NAT is int64.min, so arithmetic landing on it silently becomes NaT

A computed value equal to NPY_NAT is perfectly representable in int64, but is indistinguishable from NaT once stored. Today it is silently read back as NaT rather than raising. This affects scalars and arrays across all four datetimelike types.

The scalar and vectorized paths also disagree about which error to raise one step further out, so you can't paper over it with a try/except.

import pandas as pd
from pandas import Timedelta, Timestamp
from pandas.core.dtypes.dtypes import PeriodDtype

# one ns past the bound -> NaT; two ns past -> raises
Timedelta.min - Timedelta(1, "ns")                  # NaT          <- should raise
Timedelta.min - Timedelta(2, "ns")                  # OverflowError

pd.array([Timedelta.min]) - Timedelta(1, "ns")      # <TimedeltaArray> [NaT]   <- should raise

Timestamp.min.as_unit("ns") - Timedelta(1, "ns")    # NaT          <- should raise
pd.to_datetime([Timestamp.min]).as_unit("ns") - Timedelta(1, "ns")
                                                    # DatetimeIndex(['NaT'])   <- should raise

per = pd.Period._from_ordinal(-(2**63) + 1, PeriodDtype("ns"))
per - 1                                             # NaT          <- should raise
pd.PeriodIndex([per]) - 1                           # PeriodIndex(['NaT'])     <- should raise

Notes on the vectorized side: mul_overflowsafe got an explicit sentinel check in GH-43178 / GH-65515, but add_overflowsafe never did — which is why the timedelta64, datetime64 and Period-ordinal array paths above all still return NaT.

1b. The same condition surfaces as a message-less AssertionError

_timedelta_from_value_and_reso guards with a bare assert value != NPY_NAT, so a user hitting this gets an AssertionError with no message (and nothing at all under python -O):

Timestamp.min - Timestamp(1)      # AssertionError:
Timedelta(2, "ns") * (-(2**62))   # AssertionError:

---

2. float64 precision loss in timedelta64 arithmetic and reductions

Separate root cause: several paths round-trip int64 values through float64, which only has a 53-bit mantissa. This produces wrong answers well away from the bounds, not just at them.

2a. Scalar Timedelta.__mul__ / __truediv__

Integral operands are applied in float64, so an exactly-representable result comes back wrong, or wrongly raises:

Timedelta.min / 1      # AssertionError    <- should be Timedelta.min
Timedelta.max / 1      # OverflowError     <- should be Timedelta.max
Timedelta.min * 1.0    # AssertionError    <- should be Timedelta.min

# off-by-one above 2**53; numpy and Series/2 both give ...983
Timedelta(36028797018963967, "ns") / 2                          # 18014398509481984 ns
np.timedelta64(36028797018963967, "ns") / 2                     # 18014398509481983 ns

__mul__ also never normalizes numpy integer scalars, so the product stays an int64 multiply and wraps silently — while the Python-int and vectorized forms both raise:

Timedelta(2**62, "ns") * np.int64(4)   # Timedelta('0 days')   <- silent wrap (RuntimeWarning only)
Timedelta(2**62, "ns") * 4             # OverflowError
pd.Series([Timedelta(2**62, "ns")]) * np.int64(4)
                                       # OutOfBoundsTimedelta: Overflow in int64 multiplication

__truediv__ already does if isinstance(other, cnp.integer): other = int(other); __mul__ does not.

2b. Series.sum / DataFrame.sum accumulate in float64

nansum sets dtype_sum = np.float64 for kind == "m", so every timedelta64 sum loses precision above 2**53. A single representable value does not survive a round trip:

Series([Timedelta.min]).sum()
# NaT                       <- should be Timedelta.min

Series([Timedelta(2**62 + 1, "ns"), Timedelta(-(2**62), "ns")]).sum()
# Timedelta('0 days')       <- should be 1 ns

Series([Timedelta(2**53 + 1, "ns"), Timedelta(0, "ns")]).sum()._value
# 9007199254740992          <- should be 9007199254740993

The existing np.fabs(result) > lib.i8max guard in _wrap_results cannot catch a sum of exactly int64.min, because fabs(-2**63) and float(i8max) are the same float64.

2c. Series.cumsum / DataFrame.cumsum have no overflow check at all

_cum_func runs np.cumsum straight on the i8 view, so a running total past the bounds wraps silently to an unrelated value:

Series([Timedelta.max] * 2).cumsum()
# [Timedelta('106751 days 23:47:16.854775807'), Timedelta('-1 days +23:59:59.999999998')]

Series([Timedelta(2**62, "ns")] * 3).cumsum()
# [Timedelta('53375 days ...'), NaT, Timedelta('-53376 days ...')]

---

Expected behaviour

For (1): reject the sentinel rather than returning NaT, in both the scalar and vectorized paths, and replace the bare assert with a real exception carrying a message.

Worth deciding explicitly: the scalar Timedelta add/sub path should probably raise OverflowError rather than OutOfBoundsTimedelta, since OutOfBoundsTimedelta is a ValueError subclass and so would not actually make the scalar and array paths agree — which is the point of the change.

For (2): compute exactly in int64 and raise on genuine overflow, instead of round-tripping through float64.

Not covered above

groupby(...).sum() has the same sentinel hole, via a separate accumulator in libgroupby's group_sum rather than nanops:

half = Timedelta(-(2**62), "ns")
Series([half, half]).groupby([0, 0]).sum()   # [NaT]