diff --git a/quaddtype/meson.build b/quaddtype/meson.build index 75e9e38..6c7c95c 100644 --- a/quaddtype/meson.build +++ b/quaddtype/meson.build @@ -24,6 +24,13 @@ incdir_numpy = run_command(py, check : true ).stdout().strip() +# print the version of numpy being used +numpy_version = run_command(py, + ['-c', 'import numpy; print(numpy.__version__)'], + check : true +).stdout().strip() +message('Using NumPy version: ' + numpy_version) + npymath_path = incdir_numpy / '..' / 'lib' npymath_lib = c.find_library('npymath', dirs: npymath_path) @@ -117,6 +124,8 @@ srcs = [ 'numpy_quaddtype/src/umath/promoters.hpp', 'numpy_quaddtype/src/umath/matmul.h', 'numpy_quaddtype/src/umath/matmul.cpp', + 'numpy_quaddtype/src/umath/frexp_op.h', + 'numpy_quaddtype/src/umath/frexp_op.cpp', ] py.install_sources( diff --git a/quaddtype/numpy_quaddtype/src/dtype.c b/quaddtype/numpy_quaddtype/src/dtype.c index 7d0fad3..b10d16d 100644 --- a/quaddtype/numpy_quaddtype/src/dtype.c +++ b/quaddtype/numpy_quaddtype/src/dtype.c @@ -184,7 +184,6 @@ static PyType_Slot QuadPrecDType_Slots[] = { {NPY_DT_setitem, &quadprec_setitem}, {NPY_DT_getitem, &quadprec_getitem}, {NPY_DT_default_descr, &quadprec_default_descr}, - {NPY_DT_PyArray_ArrFuncs_dotfunc, NULL}, {0, NULL}}; static PyObject * diff --git a/quaddtype/numpy_quaddtype/src/scalar.c b/quaddtype/numpy_quaddtype/src/scalar.c index 6d82d19..8c83f3f 100644 --- a/quaddtype/numpy_quaddtype/src/scalar.c +++ b/quaddtype/numpy_quaddtype/src/scalar.c @@ -253,5 +253,6 @@ PyTypeObject QuadPrecision_Type = { int init_quadprecision_scalar(void) { + QuadPrecision_Type.tp_base = &PyFloatingArrType_Type; return PyType_Ready(&QuadPrecision_Type); } \ No newline at end of file diff --git a/quaddtype/numpy_quaddtype/src/umath/frexp_op.cpp b/quaddtype/numpy_quaddtype/src/umath/frexp_op.cpp new file mode 100644 index 0000000..4e56090 --- /dev/null +++ b/quaddtype/numpy_quaddtype/src/umath/frexp_op.cpp @@ -0,0 +1,182 @@ +#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API +#define PY_UFUNC_UNIQUE_SYMBOL QuadPrecType_UFUNC_API +#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION +#define NPY_TARGET_VERSION NPY_2_0_API_VERSION +#define NO_IMPORT_ARRAY +#define NO_IMPORT_UFUNC + +extern "C" { +#include +#include + +#include "numpy/arrayobject.h" +#include "numpy/ndarraytypes.h" +#include "numpy/ufuncobject.h" +#include "numpy/dtype_api.h" +} +#include "../quad_common.h" +#include "../scalar.h" +#include "../dtype.h" +#include "../ops.hpp" + +// Forward declarations for frexp operations +static Sleef_quad quad_frexp_mantissa(const Sleef_quad *op, int *exp); +static long double ld_frexp_mantissa(const long double *op, int *exp); + +static Sleef_quad +quad_frexp_mantissa(const Sleef_quad *op, int *exp) +{ + return Sleef_frexpq1(*op, exp); +} + +static long double +ld_frexp_mantissa(const long double *op, int *exp) +{ + return frexpl(*op, exp); +} + +static NPY_CASTING +quad_frexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[], + PyArray_Descr *const given_descrs[], PyArray_Descr *loop_descrs[], + npy_intp *NPY_UNUSED(view_offset)) +{ + Py_INCREF(given_descrs[0]); + loop_descrs[0] = given_descrs[0]; + + // Output 1: QuadPrecDType (mantissa) + if (given_descrs[1] == NULL) { + Py_INCREF(given_descrs[0]); + loop_descrs[1] = given_descrs[0]; + } + else { + Py_INCREF(given_descrs[1]); + loop_descrs[1] = given_descrs[1]; + } + + // Output 2: Int32 (exponent) + if (given_descrs[2] == NULL) { + loop_descrs[2] = PyArray_DescrFromType(NPY_INT32); + } + else { + Py_INCREF(given_descrs[2]); + loop_descrs[2] = given_descrs[2]; + } + + return NPY_NO_CASTING; +} + +static int +quad_frexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + NpyAuxData *auxdata) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *mantissa_ptr = data[1]; + char *exp_ptr = data[2]; + npy_intp in_stride = strides[0]; + npy_intp mantissa_stride = strides[1]; + npy_intp exp_stride = strides[2]; + + QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; + QuadBackendType backend = descr->backend; + size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double); + + quad_value in, mantissa; + int exp; + while (N--) { + memcpy(&in, in_ptr, elem_size); + + if (backend == BACKEND_SLEEF) { + mantissa.sleef_value = quad_frexp_mantissa(&in.sleef_value, &exp); + } + else { + mantissa.longdouble_value = ld_frexp_mantissa(&in.longdouble_value, &exp); + } + + memcpy(mantissa_ptr, &mantissa, elem_size); + *(npy_int32 *)exp_ptr = (npy_int32)exp; + + in_ptr += in_stride; + mantissa_ptr += mantissa_stride; + exp_ptr += exp_stride; + } + return 0; +} + +static int +quad_frexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data[], + npy_intp const dimensions[], npy_intp const strides[], + NpyAuxData *auxdata) +{ + npy_intp N = dimensions[0]; + char *in_ptr = data[0]; + char *mantissa_ptr = data[1]; + char *exp_ptr = data[2]; + npy_intp in_stride = strides[0]; + npy_intp mantissa_stride = strides[1]; + npy_intp exp_stride = strides[2]; + + QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0]; + QuadBackendType backend = descr->backend; + + int exp; + while (N--) { + if (backend == BACKEND_SLEEF) { + *(Sleef_quad *)mantissa_ptr = quad_frexp_mantissa((Sleef_quad *)in_ptr, &exp); + } + else { + *(long double *)mantissa_ptr = ld_frexp_mantissa((long double *)in_ptr, &exp); + } + + *(npy_int32 *)exp_ptr = (npy_int32)exp; + + in_ptr += in_stride; + mantissa_ptr += mantissa_stride; + exp_ptr += exp_stride; + } + return 0; +} + +int +create_quad_frexp_ufunc(PyObject *numpy) +{ + PyObject *ufunc = PyObject_GetAttrString(numpy, "frexp"); + if (ufunc == NULL) { + return -1; + } + + PyArray_DTypeMeta *dtypes[3] = {&QuadPrecDType, &QuadPrecDType, &PyArray_Int32DType}; + + PyType_Slot slots[] = { + {NPY_METH_resolve_descriptors, (void *)&quad_frexp_resolve_descriptors}, + {NPY_METH_strided_loop, (void *)&quad_frexp_strided_loop_aligned}, + {NPY_METH_unaligned_strided_loop, (void *)&quad_frexp_strided_loop_unaligned}, + {0, NULL}}; + + PyArrayMethod_Spec Spec = { + .name = "quad_frexp", + .nin = 1, + .nout = 2, + .casting = NPY_NO_CASTING, + .flags = NPY_METH_SUPPORTS_UNALIGNED, + .dtypes = dtypes, + .slots = slots, + }; + + if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) { + return -1; + } + + Py_DECREF(ufunc); + return 0; +} + +int +init_quad_frexp(PyObject *numpy) +{ + if (create_quad_frexp_ufunc(numpy) < 0) { + return -1; + } + return 0; +} \ No newline at end of file diff --git a/quaddtype/numpy_quaddtype/src/umath/frexp_op.h b/quaddtype/numpy_quaddtype/src/umath/frexp_op.h new file mode 100644 index 0000000..87031d8 --- /dev/null +++ b/quaddtype/numpy_quaddtype/src/umath/frexp_op.h @@ -0,0 +1,9 @@ +#ifndef _QUADDTYPE_FREXP_OP_H +#define _QUADDTYPE_FREXP_OP_H + +#include + +int +init_quad_frexp(PyObject *numpy); + +#endif \ No newline at end of file diff --git a/quaddtype/numpy_quaddtype/src/umath/umath.cpp b/quaddtype/numpy_quaddtype/src/umath/umath.cpp index 0c72d32..cfd67fa 100644 --- a/quaddtype/numpy_quaddtype/src/umath/umath.cpp +++ b/quaddtype/numpy_quaddtype/src/umath/umath.cpp @@ -24,6 +24,7 @@ extern "C" { #include "binary_ops.h" #include "comparison_ops.h" #include "matmul.h" +#include "frexp_op.h" // helper debugging function static const char * @@ -113,6 +114,11 @@ init_quad_umath(void) goto err; } + if (init_quad_frexp(numpy) < 0) { + PyErr_SetString(PyExc_RuntimeError, "Failed to initialize quad frexp operation"); + goto err; + } + Py_DECREF(numpy); return 0; diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index a761df1..1ac5907 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -595,8 +595,277 @@ def test_hyperbolic_functions(op, val): rtol = 1e-13 if abs(float_result) < 1e100 else 1e-10 np.testing.assert_allclose(float(quad_result), float_result, rtol=rtol, atol=1e-15, err_msg=f"Value mismatch for {op}({val})") - + # Check sign for zero results if float_result == 0.0: assert np.signbit(float_result) == np.signbit( quad_result), f"Zero sign mismatch for {op}({val})" + + +@pytest.mark.parametrize("backend", ["sleef", "longdouble"]) +@pytest.mark.parametrize("val", [ + # Basic test values + 1.0, 2.0, 4.0, 8.0, 16.0, + # Fractional values + 0.5, 0.25, 0.125, 0.0625, + # Negative values + -1.0, -2.0, -0.5, -0.25, + # Powers of 2 near boundaries + 1024.0, 2048.0, 0.00048828125, # 2^10, 2^11, 2^-11 + # Large values + 1e10, 1e20, 1e30, + # Small values + 1e-10, 1e-20, 1e-30, + # Edge cases + np.finfo(np.float64).max, np.finfo(np.float64).min, + np.finfo(np.float64).smallest_normal, np.finfo(np.float64).smallest_subnormal, +]) +def test_frexp_finite_values(backend, val): + """Test frexp for finite values comparing against numpy's float64 frexp.""" + quad_dtype = QuadPrecDType(backend=backend) + quad_arr = np.array([val], dtype=quad_dtype) + + # Get frexp results for quad precision + quad_mantissa, quad_exponent = np.frexp(quad_arr) + + # Get frexp results for float64 for comparison + float64_mantissa, float64_exponent = np.frexp(np.array([val], dtype=np.float64)) + + # Convert results for comparison + quad_mantissa_float = float(quad_mantissa[0]) + quad_exponent_int = int(quad_exponent[0]) + + # For finite values, mantissa should be in [0.5, 1) and val = mantissa * 2^exponent + if val != 0.0: + # Convert results for comparison (do this early) + quad_mantissa_float = float(quad_mantissa[0]) if abs(val) <= 1e300 else 0.0 + + # For very large values that might overflow Python float, check using quad precision + if abs(val) > 1e300: # Too large for safe Python float conversion + # Use quad precision comparisons + half_quad = np.array([0.5], dtype=quad_dtype)[0] + one_quad = np.array([1.0], dtype=quad_dtype)[0] + abs_mantissa = np.abs(quad_mantissa[0]) + assert np.greater_equal(abs_mantissa, half_quad), f"Mantissa magnitude too small for {val}" + assert np.less(abs_mantissa, one_quad), f"Mantissa magnitude too large for {val}" + else: + # Safe to convert to Python float for smaller values + assert 0.5 <= abs(quad_mantissa_float) < 1.0, f"Mantissa {quad_mantissa_float} not in [0.5, 1) for {val}" + + # For very large values, avoid overflow in reconstruction test + # Just check that the relationship holds conceptually + if abs(quad_exponent_int) < 1000 and abs(val) <= 1e300: # Avoid overflow in 2**exponent + # Reconstruct the original value + reconstructed = quad_mantissa_float * (2.0 ** quad_exponent_int) + np.testing.assert_allclose(reconstructed, val, rtol=1e-14, atol=1e-300, + err_msg=f"frexp reconstruction failed for {val}") + + if abs(val) <= np.finfo(np.float64).max and abs(val) >= np.finfo(np.float64).smallest_normal and abs(val) <= 1e300: + np.testing.assert_allclose(quad_mantissa_float, float64_mantissa[0], rtol=1e-14, + err_msg=f"Mantissa mismatch for {val}") + assert quad_exponent_int == float64_exponent[0], f"Exponent mismatch for {val}" + + +@pytest.mark.parametrize("backend", ["sleef", "longdouble"]) +def test_frexp_special_values(backend): + """Test frexp for special values (zero, inf, nan).""" + quad_dtype = QuadPrecDType(backend=backend) + + # Test zero + zero_arr = np.array([0.0], dtype=quad_dtype) + mantissa, exponent = np.frexp(zero_arr) + assert float(mantissa[0]) == 0.0, "frexp(0.0) mantissa should be 0.0" + assert int(exponent[0]) == 0, "frexp(0.0) exponent should be 0" + + # Test negative zero + neg_zero_arr = np.array([-0.0], dtype=quad_dtype) + mantissa, exponent = np.frexp(neg_zero_arr) + mantissa_val = float(mantissa[0]) + assert mantissa_val == 0.0, "frexp(-0.0) mantissa should be 0.0" + assert int(exponent[0]) == 0, "frexp(-0.0) exponent should be 0" + + # Test positive infinity + inf_arr = np.array([np.inf], dtype=quad_dtype) + mantissa, exponent = np.frexp(inf_arr) + assert np.isinf(float(mantissa[0])), "frexp(inf) mantissa should be inf" + assert float(mantissa[0]) > 0, "frexp(inf) mantissa should be positive" + + # Test negative infinity + neg_inf_arr = np.array([-np.inf], dtype=quad_dtype) + mantissa, exponent = np.frexp(neg_inf_arr) + assert np.isinf(float(mantissa[0])), "frexp(-inf) mantissa should be inf" + assert float(mantissa[0]) < 0, "frexp(-inf) mantissa should be negative" + + # Test NaN + nan_arr = np.array([np.nan], dtype=quad_dtype) + mantissa, exponent = np.frexp(nan_arr) + assert np.isnan(float(mantissa[0])), "frexp(nan) mantissa should be nan" + + +@pytest.mark.parametrize("backend", ["sleef", "longdouble"]) +def test_frexp_array_operations(backend): + """Test frexp on arrays of different sizes and shapes.""" + quad_dtype = QuadPrecDType(backend=backend) + + # Test 1D array + values_1d = [0.5, 1.0, 2.0, 4.0, 8.0, -1.0, -0.25] + quad_arr_1d = np.array(values_1d, dtype=quad_dtype) + mantissa_1d, exponent_1d = np.frexp(quad_arr_1d) + + assert mantissa_1d.shape == quad_arr_1d.shape, "Mantissa shape should match input" + assert exponent_1d.shape == quad_arr_1d.shape, "Exponent shape should match input" + assert mantissa_1d.dtype == quad_dtype, "Mantissa should have quad precision dtype" + assert exponent_1d.dtype == np.int32, "Exponent should have int32 dtype" + + # Verify each element + for i, val in enumerate(values_1d): + if val != 0.0: + mantissa_float = float(mantissa_1d[i]) + exponent_int = int(exponent_1d[i]) + reconstructed = mantissa_float * (2.0 ** exponent_int) + np.testing.assert_allclose(reconstructed, val, rtol=1e-14) + + # Test 2D array + values_2d = np.array([[1.0, 2.0, 4.0], [0.5, 0.25, 8.0]], dtype=float) + quad_arr_2d = values_2d.astype(quad_dtype) + mantissa_2d, exponent_2d = np.frexp(quad_arr_2d) + + assert mantissa_2d.shape == quad_arr_2d.shape, "2D mantissa shape should match input" + assert exponent_2d.shape == quad_arr_2d.shape, "2D exponent shape should match input" + + +@pytest.mark.parametrize("backend", ["sleef", "longdouble"]) +def test_frexp_dtype_consistency(backend): + """Test that frexp output dtypes are consistent.""" + quad_dtype = QuadPrecDType(backend=backend) + + # Single value + arr = np.array([3.14159], dtype=quad_dtype) + mantissa, exponent = np.frexp(arr) + + # Check mantissa dtype matches input + assert mantissa.dtype == quad_dtype, f"Mantissa dtype {mantissa.dtype} should match input {quad_dtype}" + + # Check exponent dtype is int32 + assert exponent.dtype == np.int32, f"Exponent dtype should be int32, got {exponent.dtype}" + + # Test with different input shapes + for shape in [(5,), (2, 3), (2, 2, 2)]: + arr = np.ones(shape, dtype=quad_dtype) + mantissa, exponent = np.frexp(arr) + assert mantissa.dtype == quad_dtype + assert exponent.dtype == np.int32 + assert mantissa.shape == shape + assert exponent.shape == shape + + +def test_frexp_basic_functionality(): + """Test that frexp works correctly for basic cases.""" + # This is the key test that was failing before implementing frexp + quad_dtype = QuadPrecDType() + + # Test that frexp ufunc is registered and callable + quad_arr = np.array([1.0], dtype=quad_dtype) + mantissa, exponent = np.frexp(quad_arr) + + assert len(mantissa) == 1, "frexp should return mantissa array" + assert len(exponent) == 1, "frexp should return exponent array" + assert mantissa.dtype == quad_dtype, "Mantissa should have quad dtype" + assert exponent.dtype == np.int32, "Exponent should have int32 dtype" + + # Test with epsilon value directly from constants + from numpy_quaddtype import epsilon + eps_arr = np.array([float(epsilon)], dtype=quad_dtype) + eps_mantissa, eps_exponent = np.frexp(eps_arr) + + # For binary128 (quad precision), epsilon = 2^-112 + # So frexp(epsilon) should give mantissa=0.5 and exponent=-111 + eps_mantissa_float = float(eps_mantissa[0]) + eps_exponent_int = int(eps_exponent[0]) + + assert abs(eps_mantissa_float - 0.5) < 1e-15, f"Expected mantissa 0.5 for epsilon, got {eps_mantissa_float}" + assert eps_exponent_int == -111, f"Expected exponent -111 for epsilon, got {eps_exponent_int}" + + +def test_frexp_with_quad_constants(): + """Test frexp using the constants exposed by numpy_quaddtype instead of np.finfo.""" + from numpy_quaddtype import epsilon, smallest_normal, max_value, pi, e + + quad_dtype = QuadPrecDType() + + # Test various important constants + constants_to_test = [ + ("epsilon", epsilon), + ("smallest_normal", smallest_normal), + ("pi", pi), + ("e", e), + ] + + for name, constant in constants_to_test: + # Convert constant to quad array + const_arr = np.array([float(constant)], dtype=quad_dtype) + mantissa, exponent = np.frexp(const_arr) + + # Basic checks + assert mantissa.dtype == quad_dtype, f"Mantissa dtype should be quad for {name}" + assert exponent.dtype == np.int32, f"Exponent dtype should be int32 for {name}" + + # Reconstruction check (if value is in reasonable range) + mantissa_float = float(mantissa[0]) + exponent_int = int(exponent[0]) + + if abs(exponent_int) < 1000: # Avoid overflow in Python arithmetic + reconstructed = mantissa_float * (2.0 ** exponent_int) + original_float = float(constant) + np.testing.assert_allclose(reconstructed, original_float, rtol=1e-14, + err_msg=f"Reconstruction failed for {name}") + + # Check mantissa range for non-zero values + if mantissa_float != 0.0: + assert 0.5 <= abs(mantissa_float) < 1.0, f"Mantissa {mantissa_float} not in [0.5, 1) for {name}" + + +@pytest.mark.parametrize("backend", ["sleef", "longdouble"]) +def test_frexp_edge_cases_quad_precision(backend): + """Test frexp with values specific to quad precision range.""" + quad_dtype = QuadPrecDType(backend=backend) + + # Test with more reasonable values that won't overflow/underflow + test_values = [ + 1.0, # Basic case + 2.0**100, # Large but manageable + 2.0**(-100), # Small but manageable + 2.0**1000, # Very large + 2.0**(-1000), # Very small + ] + + for val in test_values: + arr = np.array([val], dtype=quad_dtype) + mantissa, exponent = np.frexp(arr) + + # Test that mantissa is in quad precision format + assert mantissa.dtype == quad_dtype, f"Mantissa should have quad dtype, got {mantissa.dtype}" + assert exponent.dtype == np.int32, f"Exponent should be int32, got {exponent.dtype}" + + # Test reconstruction using quad precision arithmetic + # mantissa * 2^exponent should equal original value + pow2_exp = np.array([2.0], dtype=quad_dtype) ** exponent[0] + reconstructed = mantissa[0] * pow2_exp + + # Convert both to the same dtype for comparison + original_quad = np.array([val], dtype=quad_dtype) + + # For reasonable values, reconstruction should be exact + if abs(val) >= 1e-300 and abs(val) <= 1e300: + np.testing.assert_array_equal(reconstructed, original_quad[0], + err_msg=f"Reconstruction failed for {val}") + + # Test that mantissa is in correct range (using quad precision comparisons) + if not np.equal(mantissa[0], 0.0): # Skip zero case + half_quad = np.array([0.5], dtype=quad_dtype)[0] + one_quad = np.array([1.0], dtype=quad_dtype)[0] + + # |mantissa| should be in [0.5, 1.0) range + abs_mantissa = np.abs(mantissa[0]) + assert np.greater_equal(abs_mantissa, half_quad), f"Mantissa magnitude too small for {val}" + assert np.less(abs_mantissa, one_quad), f"Mantissa magnitude too large for {val}"