Fitting polynomials on epoch timestamps: subtract t0 first, even with Polynomial.fit

NumPy 2.5.4. x = 1.7e9 + 60*i (1000 minute timestamps), with y an exact polynomial of degree 1-3 in the scaled time; max relative residual over 10 random coefficient sets. np.polyfit(x, y, deg): deg 1 1e-13 to 3e-11, deg 2 9e-9 to 1e-5 with no warning, deg 3 6e-5 to 0.1 (RankWarning only here). np.polyfit(x - t0, y, deg): <= 2e-15 for every degree. Polynomial.fit(x, y, deg) maps the domain to [-1, 1] and so avoids the disaster, but not fully: it reached 3e-12 to 9e-12. That is consistent with cancellation in its map u = off + scl*x, where |off| ~ 5.7e4 and |off|*eps ~ 6e-12. Polynomial.fit(x - t0, ...) gave <= 2e-15. So subtract a reference time in float64 before any fit, and keep the shifted x for evaluation too. Never store such timestamps as float32: the spacing at 1.7e9 is 128 s, and only 469 of the 1000 minute stamps stayed distinct.