Skip to content
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
181 changes: 181 additions & 0 deletions quaddtype/numpy_quaddtype/src/scalar.c
Original file line number Diff line number Diff line change
Expand Up @@ -383,6 +383,186 @@ 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;
}
}

// this is thread-unsafe
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ngoldbaum precisely, Sleef_snprintf here is thread-unsafe, I believe just GIL won't help here as this is C routine and GIL only protects the Python objects.
Not 100% sure so need your opinion that whether GIL would be enough or we can lock this region with pthread_mutex

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess that means that it's not safe to concurrently call Sleef_snprintf simultaneously in two threads (e.g. it's not re-entrant)?

In that case, yeah, I think you need a global lock.

I'd avoid using pthreads directly because then you'd need to do something else on Windows. Instead, I'd use PyMutex on Python 3.13 and newer and PyThread_type_lock on 3.12 and older. See e.g. the use of lapack_lite_lock in numpy/linalg/umath_linalg.cpp in NumPy, which solves a similar problem with the non-reentrant lapack_lite library.

You could also switch to C++ and use a C++ standard library but then you need to be careful you don't deadlock with the GIL by making sure you release it before doing any possibly blocking calls. The built-in lock types have deadlock protection against the GIL so you don't need to go through that trouble.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also it'd be nice to add a multithreaded test for this as well. You can look at test_multithreading.py in NumPy for some patterns to use for multithreaded tests.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pushed a multithreaded test for now only for testing this. Will push the lock after the tests get done

PyObject* quad_to_pylong(Sleef_quad value)
{
char buffer[128];
// 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);
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
// NOTE: a 128-bit
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)) {which in theory seem better}
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(numerator);
return NULL;
}
PyObject *denominator = PyLong_FromLong(1);
if (denominator == NULL) {
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},
Expand All @@ -400,6 +580,7 @@ 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,
};

Expand Down
117 changes: 116 additions & 1 deletion quaddtype/tests/test_quaddtype.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
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