From c0b5b7d931046729d743e2c0a27e34fa8b67b6d2 Mon Sep 17 00:00:00 2001 From: Anton Soloviev Date: Thu, 15 Oct 2026 10:04:55 -0700 Subject: [PATCH] Promote integer samples before finite-difference arithmetic --- pynumdiff/finite_difference.py | 3 ++- pynumdiff/tests/test_integer_differences.py | 24 +++++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) create mode 100644 pynumdiff/tests/test_integer_differences.py diff --git a/pynumdiff/finite_difference.py b/pynumdiff/finite_difference.py index 7dc170a..e683e9b 100644 --- a/pynumdiff/finite_difference.py +++ b/pynumdiff/finite_difference.py @@ -23,7 +23,8 @@ def finitediff(x, dt, num_iterations=1, order=2, axis=0): if order not in [1, 2, 4]: raise ValueError("order must be 1, 2, or 4") x = np.moveaxis(x, axis, 0) # move the axis of interest to the front to simplify differencing indexing - x_hat = np.asarray(x) # allows for array-like. Preserve reference to x, for finding the final constant of integration + # Promote integer samples before subtraction and weighted stencil sums can overflow. + x_hat = x.astype(float) if np.issubdtype(x.dtype, np.integer) else np.asarray(x) dxdt_hat = np.zeros(x.shape) # preallocate reusable memory # For all but the last iteration, do the differentate->integrate smoothing loop, being careful with endpoints diff --git a/pynumdiff/tests/test_integer_differences.py b/pynumdiff/tests/test_integer_differences.py new file mode 100644 index 0000000..428eefd --- /dev/null +++ b/pynumdiff/tests/test_integer_differences.py @@ -0,0 +1,24 @@ +"""Integer samples must not overflow before floating derivative scaling.""" +import numpy as np +import pytest + +from pynumdiff import finitediff + + +@pytest.mark.parametrize('dtype', [np.int8, np.uint8, np.int16, np.uint16, np.int32, np.uint32, np.int64, np.uint64]) +@pytest.mark.parametrize('order', [1, 2, 4]) +@pytest.mark.parametrize('num_iterations', [1, 3]) +def test_integer_finite_differences(dtype, order, num_iterations): + """Check analytic quadratic derivatives and equivalent floating samples.""" + # The exact quadratic and its derivative fit in the input dtype, while + # weighted stencil intermediates can overflow even on unsigned inputs. + t = np.arange(11) + signal = np.stack([t**2, (10-t)**2]).astype(dtype) + original = signal.copy() + _, derivative = finitediff(signal, .25, order=order, num_iterations=num_iterations, axis=-1) + expected = finitediff(signal.astype(float), .25, order=order, + num_iterations=num_iterations, axis=-1) + np.testing.assert_array_equal(derivative, expected[1]) + if num_iterations == 1 and order > 1: + np.testing.assert_array_equal(derivative, np.stack([8*t, -8*(10-t)])) + np.testing.assert_array_equal(signal, original)