diff --git a/quaddtype/README.md b/quaddtype/README.md index 80fa86af..0d6051e9 100644 --- a/quaddtype/README.md +++ b/quaddtype/README.md @@ -31,6 +31,7 @@ np.array([1,2,3], dtype=QuadPrecDType("longdouble")) - **CMake** (≥3.15) - **Python 3.10+** - **Git** +- **NumPy >= 2.4** (build from source) ### Linux/Unix/macOS @@ -48,11 +49,11 @@ pip install numpy pytest # export CFLAGS="-DDISABLE_QUADBLAS" # export CXXFLAGS="-DDISABLE_QUADBLAS" -python -m pip install . -v +python -m pip install . -v --no-build-isolation # Run the tests cd .. -python -m pytest +python -m pytest/quaddtype/tests/ ``` ### Windows @@ -94,14 +95,14 @@ python -m pytest ```powershell # Build and install the package - python -m pip install . -v + python -m pip install . -v --no-build-isolation ``` 5. **Test Installation** ```powershell # Run tests - pytest -s tests/ + pytest -s ..\quaddtype\tests\ ``` 6. **QBLAS Disabled**: QuadBLAS optimization is automatically disabled on Windows builds due to MSVC compatibility issues. This is handled by the `-DDISABLE_QUADBLAS` compiler flag. @@ -112,3 +113,58 @@ python -m pytest - VS 2017: `"Visual Studio 15 2017"` 8. **Architecture**: The instructions are for x64. For x86 builds, change `-A x64` to `-A Win32`. + +## Building with ThreadSanitizer (TSan) + +This is a development feature to help detect threading issues. To build `numpy-quaddtype` with TSan enabled, follow these steps: + +> Use of clang is recommended with machine NOT supporting `libquadmath` (like ARM64). Set the compiler to clang/clang++ before proceeding. +> ```bash +> export CC=clang +> export CXX=clang++ +> ``` + +1. Compile free-threaded CPython with TSan support. Follow the [Python Free-Threading Guide](https://py-free-threading.github.io/thread_sanitizer/#compile-free-threaded-cpython-with-tsan) for detailed instructions. +2. Create and activate a virtual environment using the TSan-enabled Python build. +3. Installing dependencies: + + ```bash + python -m pip install meson meson-python wheel ninja + # Need NumPy built with TSan as well + python -m pip install "numpy @ git+https://github.com/numpy/numpy" -C'setup-args=-Db_sanitize=thread' + ``` +4. Building SLEEF with TSan: + + ```bash + # clone the repository + git clone -b 3.8 https://github.com/shibatch/sleef.git + cd sleef + + # Build SLEEF with TSan + cmake \ + -DCMAKE_C_COMPILER=clang \ + -DCMAKE_CXX_COMPILER=clang++ \ + -DCMAKE_C_FLAGS="-fsanitize=thread -g -O1" \ + -DCMAKE_CXX_FLAGS="-fsanitize=thread -g -O1" \ + -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=thread" \ + -DCMAKE_SHARED_LINKER_FLAGS="-fsanitize=thread" \ + -DSLEEF_BUILD_QUAD=ON \ + -DSLEEF_BUILD_TESTS=OFF \ + -S . -B build + + cmake --build build -j + + # Install the built library and headers into the system path (/usr/local) + sudo cmake --install build --prefix=/usr/local + ``` +5. Build and install `numpy-quaddtype` with TSan: + + ```bash + # SLEEF is already installed with TSan, we need to provide proper flags to numpy-quaddtype's meson file + # So that it does not build SLEEF again and use the installed one. + + export CFLAGS="-fsanitize=thread -g -O0" + export CXXFLAGS="-fsanitize=thread -g -O0" + export LDFLAGS="-fsanitize=thread" + python -m pip install . -vv --no-build-isolation -Csetup-args=-Db_sanitize=thread + ``` \ No newline at end of file diff --git a/quaddtype/numpy_quaddtype/_quaddtype_main.pyi b/quaddtype/numpy_quaddtype/_quaddtype_main.pyi index 566ff996..831c073b 100644 --- a/quaddtype/numpy_quaddtype/_quaddtype_main.pyi +++ b/quaddtype/numpy_quaddtype/_quaddtype_main.pyi @@ -1,5 +1,5 @@ from typing import Any, Literal, TypeAlias, final, overload - +import builtins import numpy as np from numpy._typing import _128Bit # pyright: ignore[reportPrivateUsage] from typing_extensions import Never, Self, override @@ -157,9 +157,10 @@ class QuadPrecision(np.floating[_128Bit]): # NOTE: is_integer() and as_integer_ratio() are defined on numpy.floating in the # stubs, but don't exist at runtime. And because QuadPrecision does not implement # them, we use this hacky workaround to emulate their absence. - # TODO: Remove after https://github.com/numpy/numpy-user-dtypes/issues/216 - is_integer: Never # pyright: ignore[reportIncompatibleMethodOverride] - as_integer_ratio: Never # pyright: ignore[reportIncompatibleMethodOverride] + @override + def is_integer(self, /) -> builtins.bool: ... + @override + def as_integer_ratio(self, /) -> tuple[int, int]: ... # def is_longdouble_128() -> bool: ... diff --git a/quaddtype/numpy_quaddtype/src/scalar.c b/quaddtype/numpy_quaddtype/src/scalar.c index 2be50789..11c5529c 100644 --- a/quaddtype/numpy_quaddtype/src/scalar.c +++ b/quaddtype/numpy_quaddtype/src/scalar.c @@ -22,6 +22,18 @@ // src: https://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format #define SLEEF_QUAD_DECIMAL_DIG 36 +#if PY_VERSION_HEX < 0x30d00b3 +static PyThread_type_lock sleef_lock; +#define LOCK_SLEEF PyThread_acquire_lock(sleef_lock, WAIT_LOCK) +#define UNLOCK_SLEEF PyThread_release_lock(sleef_lock) +#else +static PyMutex sleef_lock = {0}; +#define LOCK_SLEEF PyMutex_Lock(&sleef_lock) +#define UNLOCK_SLEEF PyMutex_Unlock(&sleef_lock) +#endif + + + QuadPrecisionObject * QuadPrecision_raw_new(QuadBackendType backend) @@ -383,6 +395,187 @@ QuadPrecision_get_imag(QuadPrecisionObject *self, void *closure) return (PyObject *)QuadPrecision_raw_new(self->backend); } +// Method implementations for float compatibility +static PyObject * +QuadPrecision_is_integer(QuadPrecisionObject *self, PyObject *Py_UNUSED(ignored)) +{ + Sleef_quad value; + + if (self->backend == BACKEND_SLEEF) { + value = self->value.sleef_value; + } + else { + // lets also tackle ld from sleef functions as well + value = Sleef_cast_from_doubleq1((double)self->value.longdouble_value); + } + + if(Sleef_iunordq1(value, value)) { + Py_RETURN_FALSE; + } + + // Check if value is finite (not inf or nan) + Sleef_quad abs_value = Sleef_fabsq1(value); + Sleef_quad pos_inf = sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384); + int32_t is_finite = Sleef_icmpltq1(abs_value, pos_inf); + + if (!is_finite) { + Py_RETURN_FALSE; + } + + // Check if value equals its truncated version + Sleef_quad truncated = Sleef_truncq1(value); + int32_t is_equal = Sleef_icmpeqq1(value, truncated); + + if (is_equal) { + Py_RETURN_TRUE; + } + else { + Py_RETURN_FALSE; + } +} + +PyObject* quad_to_pylong(Sleef_quad value) +{ + char buffer[128]; + + // Sleef_snprintf call is thread-unsafe + LOCK_SLEEF; + // Format as integer (%.0Qf gives integer with no decimal places) + // Q modifier means pass Sleef_quad by value + int written = Sleef_snprintf(buffer, sizeof(buffer), "%.0Qf", value); + UNLOCK_SLEEF; + if (written < 0 || written >= sizeof(buffer)) { + PyErr_SetString(PyExc_RuntimeError, "Failed to convert quad to string"); + return NULL; + } + + PyObject *result = PyLong_FromString(buffer, NULL, 10); + + if (result == NULL) { + PyErr_SetString(PyExc_RuntimeError, "Failed to parse integer string"); + return NULL; + } + + return result; +} + +// inspired by the CPython implementation +// https://github.com/python/cpython/blob/ac1ffd77858b62d169a08040c08aa5de26e145ac/Objects/floatobject.c#L1503C1-L1572C2 +static PyObject * +QuadPrecision_as_integer_ratio(QuadPrecisionObject *self, PyObject *Py_UNUSED(ignored)) +{ + + Sleef_quad value; + Sleef_quad pos_inf = sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384); + const int FLOAT128_PRECISION = 113; + + if (self->backend == BACKEND_SLEEF) { + value = self->value.sleef_value; + } + else { + // lets also tackle ld from sleef functions as well + value = Sleef_cast_from_doubleq1((double)self->value.longdouble_value); + } + + if(Sleef_iunordq1(value, value)) { + PyErr_SetString(PyExc_ValueError, "Cannot convert NaN to integer ratio"); + return NULL; + } + if(Sleef_icmpgeq1(Sleef_fabsq1(value), pos_inf)) { + PyErr_SetString(PyExc_OverflowError, "Cannot convert infinite value to integer ratio"); + return NULL; + } + + // Sleef_value == float_part * 2**exponent exactly + int exponent; + Sleef_quad mantissa = Sleef_frexpq1(value, &exponent); // within [0.5, 1.0) + + /* + CPython loops for 300 (some huge number) to make sure + float_part gets converted to the floor(float_part) i.e. near integer as + + for (i=0; i<300 && float_part != floor(float_part) ; i++) { + float_part *= 2.0; + exponent--; + } + + It seems highly inefficient from performance perspective, maybe they pick 300 for future-proof + or If FLT_RADIX != 2, the 300 steps may leave a tiny fractional part + + Another way can be doing as: + ``` + mantissa = ldexpq(mantissa, FLOAT128_PRECISION); + exponent -= FLOAT128_PRECISION; + ``` + This should work but give non-simplified, huge integers (although they also come down to same representation) + We can also do gcd to find simplified values, but it'll add more O(log(N)) + For the sake of simplicity and fixed 128-bit nature, we will loop till 113 only + */ + + for (int i = 0; i < FLOAT128_PRECISION && !Sleef_icmpeqq1(mantissa, Sleef_floorq1(mantissa)); i++) { + mantissa = Sleef_mulq1_u05(mantissa, Sleef_cast_from_doubleq1(2.0)); + exponent--; + } + + // numerator and denominators can't fit in int + // convert items to PyLongObject from string instead + PyObject *py_exp = PyLong_FromLongLong(Py_ABS(exponent)); + if(py_exp == NULL) + { + return NULL; + } + + PyObject *numerator = quad_to_pylong(mantissa); + if(numerator == NULL) + { + Py_DECREF(py_exp); + return NULL; + } + PyObject *denominator = PyLong_FromLong(1); + if (denominator == NULL) { + Py_DECREF(py_exp); + Py_DECREF(numerator); + return NULL; + } + + // fold in 2**exponent + if(exponent > 0) + { + PyObject *new_num = PyNumber_Lshift(numerator, py_exp); + Py_DECREF(numerator); + if(new_num == NULL) + { + Py_DECREF(denominator); + Py_DECREF(py_exp); + return NULL; + } + numerator = new_num; + } + else + { + PyObject *new_denom = PyNumber_Lshift(denominator, py_exp); + Py_DECREF(denominator); + if(new_denom == NULL) + { + Py_DECREF(numerator); + Py_DECREF(py_exp); + return NULL; + } + denominator = new_denom; + } + + Py_DECREF(py_exp); + return PyTuple_Pack(2, numerator, denominator); +} + +static PyMethodDef QuadPrecision_methods[] = { + {"is_integer", (PyCFunction)QuadPrecision_is_integer, METH_NOARGS, + "Return True if the value is an integer."}, + {"as_integer_ratio", (PyCFunction)QuadPrecision_as_integer_ratio, METH_NOARGS, + "Return a pair of integers whose ratio is exactly equal to the original value."}, + {NULL, NULL, 0, NULL} /* Sentinel */ +}; + static PyGetSetDef QuadPrecision_getset[] = { {"real", (getter)QuadPrecision_get_real, NULL, "Real part of the scalar", NULL}, {"imag", (getter)QuadPrecision_get_imag, NULL, "Imaginary part of the scalar (always 0 for real types)", NULL}, @@ -400,12 +593,20 @@ PyTypeObject QuadPrecision_Type = { .tp_as_number = &quad_as_scalar, .tp_as_buffer = &QuadPrecision_as_buffer, .tp_richcompare = (richcmpfunc)quad_richcompare, + .tp_methods = QuadPrecision_methods, .tp_getset = QuadPrecision_getset, }; int init_quadprecision_scalar(void) { +#if PY_VERSION_HEX < 0x30d00b3 + sleef_lock = PyThread_allocate_lock(); + if (sleef_lock == NULL) { + PyErr_NoMemory(); + return -1; + } +#endif QuadPrecision_Type.tp_base = &PyFloatingArrType_Type; return PyType_Ready(&QuadPrecision_Type); } \ No newline at end of file diff --git a/quaddtype/reinstall.sh b/quaddtype/reinstall.sh index 9144f05b..b6f5e043 100755 --- a/quaddtype/reinstall.sh +++ b/quaddtype/reinstall.sh @@ -8,7 +8,11 @@ if [ -d "build/" ]; then rm -rf subprojects/sleef fi -# export CFLAGS="-g -O0" -# export CXXFLAGS="-g -O0" python -m pip uninstall -y numpy_quaddtype python -m pip install . -vv --no-build-isolation 2>&1 | tee build_log.txt + +# for debugging and TSAN builds, comment the above line and uncomment all below: +# export CFLAGS="-fsanitize=thread -g -O0" +# export CXXFLAGS="-fsanitize=thread -g -O0" +# export LDFLAGS="-fsanitize=thread" +# CC=clang CXX=clang++ python -m pip install . -vv --no-build-isolation -Csetup-args=-Db_sanitize=thread 2>&1 | tee build_log.txt \ No newline at end of file diff --git a/quaddtype/subprojects/packagefiles/sleef/meson.build b/quaddtype/subprojects/packagefiles/sleef/meson.build index 20faeff4..68d9384f 100644 --- a/quaddtype/subprojects/packagefiles/sleef/meson.build +++ b/quaddtype/subprojects/packagefiles/sleef/meson.build @@ -12,6 +12,7 @@ if host_machine.system() == 'windows' parallel_flag = [] endif +# For building sleef with TSan, delete the sleef subproject and follow the README instructions to build sleef externally. sleef_configure = run_command([ cmake, '-S', meson.current_source_dir(), diff --git a/quaddtype/tests/test_multithreading.py b/quaddtype/tests/test_multithreading.py new file mode 100644 index 00000000..90eb257b --- /dev/null +++ b/quaddtype/tests/test_multithreading.py @@ -0,0 +1,32 @@ +import concurrent.futures +import threading + +import pytest + +import numpy as np +from numpy._core import _rational_tests +from numpy._core.tests.test_stringdtype import random_unicode_string_list +from numpy.testing import IS_64BIT, IS_WASM +from numpy.testing._private.utils import run_threaded + +if IS_WASM: + pytest.skip(allow_module_level=True, reason="no threading support in wasm") + +from numpy_quaddtype import * + + +def test_as_integer_ratio_reconstruction(): + """Multi-threaded test that as_integer_ratio() can reconstruct the original value.""" + + def test(barrier): + barrier.wait() # All threads start simultaneously + values = ["3.14", "0.1", "1.414213562373095", "2.718281828459045", + "-1.23456789", "1000.001", "0.0001", "1e20", "1.23e15", "1e-30", pi] + for val in values: + quad_val = QuadPrecision(val) + num, denom = quad_val.as_integer_ratio() + # todo: can remove str converstion after merging PR #213 + reconstructed = QuadPrecision(str(num)) / QuadPrecision(str(denom)) + assert reconstructed == quad_val + + run_threaded(test, pass_barrier=True, max_workers=64, outer_iterations=100) \ No newline at end of file diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py index d657d8df..ae8d5b5b 100644 --- a/quaddtype/tests/test_quaddtype.py +++ b/quaddtype/tests/test_quaddtype.py @@ -3507,4 +3507,119 @@ def test_fromfile_single_value(self): expected = np.array([42.0], dtype=QuadPrecDType(backend='sleef')) np.testing.assert_array_equal(result, expected) finally: - os.unlink(fname) \ No newline at end of file + os.unlink(fname) + + +class Test_Is_Integer_Methods: + """Test suite for float compatibility methods: is_integer() and as_integer_ratio().""" + + @pytest.mark.parametrize("value,expected", [ + # Positive integers + ("1.0", True), + ("42.0", True), + ("1000.0", True), + # Negative integers + ("-1.0", True), + ("-42.0", True), + # Zero + ("0.0", True), + ("-0.0", True), + # Large integers + ("1e20", True), + ("123456789012345678901234567890", True), + # Fractional values + ("1.5", False), + ("3.14", False), + ("-2.5", False), + ("0.1", False), + ("0.0001", False), + # Values close to integers but not exact + ("1.0000000000001", False), + ("0.9999999999999", False), + # Special values + ("inf", False), + ("-inf", False), + ("nan", False), + ]) + def test_is_integer(self, value, expected): + """Test is_integer() returns correct result for various values.""" + assert QuadPrecision(value).is_integer() == expected + + @pytest.mark.parametrize("value", ["0.0", "1.0", "1.5", "-3.0", "-3.7", "42.0"]) + def test_is_integer_compatibility_with_float(self, value): + """Test is_integer() matches behavior of Python's float.""" + quad_val = QuadPrecision(value) + float_val = float(value) + assert quad_val.is_integer() == float_val.is_integer() + + @pytest.mark.parametrize("value,expected_num,expected_denom", [ + ("1.0", 1, 1), + ("42.0", 42, 1), + ("-5.0", -5, 1), + ("0.0", 0, 1), + ("-0.0", 0, 1), + ]) + def test_as_integer_ratio_integers(self, value, expected_num, expected_denom): + """Test as_integer_ratio() for integer values.""" + num, denom = QuadPrecision(value).as_integer_ratio() + assert num == expected_num and denom == expected_denom + + @pytest.mark.parametrize("value,expected_ratio", [ + ("0.5", 0.5), + ("0.25", 0.25), + ("1.5", 1.5), + ("-2.5", -2.5), + ]) + def test_as_integer_ratio_fractional(self, value, expected_ratio): + """Test as_integer_ratio() for fractional values.""" + num, denom = QuadPrecision(value).as_integer_ratio() + assert QuadPrecision(str(num)) / QuadPrecision(str(denom)) == QuadPrecision(str(expected_ratio)) + assert denom > 0 # Denominator should always be positive + + @pytest.mark.parametrize("value", [ + "3.14", "0.1", "1.414213562373095", "2.718281828459045", + "-1.23456789", "1000.001", "0.0001", "1e20", "1.23e15", "1e-30", quad_pi + ]) + def test_as_integer_ratio_reconstruction(self, value): + """Test that as_integer_ratio() can reconstruct the original value.""" + quad_val = QuadPrecision(value) + num, denom = quad_val.as_integer_ratio() + # todo: can remove str converstion after merging PR #213 + reconstructed = QuadPrecision(str(num)) / QuadPrecision(str(denom)) + assert reconstructed == quad_val + + def test_as_integer_ratio_return_types(self): + """Test that as_integer_ratio() returns Python ints.""" + num, denom = QuadPrecision("3.14").as_integer_ratio() + assert isinstance(num, int) + assert isinstance(denom, int) + + @pytest.mark.parametrize("value", ["-1.0", "-3.14", "-0.5", "1.0", "3.14", "0.5"]) + def test_as_integer_ratio_denominator_positive(self, value): + """Test that denominator is always positive.""" + num, denom = QuadPrecision(value).as_integer_ratio() + assert denom > 0 + + @pytest.mark.parametrize("value,exception,match", [ + ("inf", OverflowError, "Cannot convert infinite value to integer ratio"), + ("-inf", OverflowError, "Cannot convert infinite value to integer ratio"), + ("nan", ValueError, "Cannot convert NaN to integer ratio"), + ]) + def test_as_integer_ratio_special_values_raise(self, value, exception, match): + """Test that as_integer_ratio() raises appropriate errors for special values.""" + with pytest.raises(exception, match=match): + QuadPrecision(value).as_integer_ratio() + + @pytest.mark.parametrize("value", ["1.0", "0.5", "3.14", "-2.5", "0.0"]) + def test_as_integer_ratio_compatibility_with_float(self, value): + """Test as_integer_ratio() matches behavior of Python's float where possible.""" + quad_val = QuadPrecision(value) + float_val = float(value) + + quad_num, quad_denom = quad_val.as_integer_ratio() + float_num, float_denom = float_val.as_integer_ratio() + + # The ratios should be equal + quad_ratio = quad_num / quad_denom + float_ratio = float_num / float_denom + assert abs(quad_ratio - float_ratio) < 1e-15 \ No newline at end of file