Merge pull request #8 from seberg/mpfdtype

ENH: Add `MPFloat.prec` and `MPFDType.prec`, also add basic scalar tests

Author: ngoldbaum

.github/workflows/ci.yml

@@ -29,6 +29,15 @@ jobs:
         working-directory: metadatadtype
         run: |
           pytest -vvv --color=yes
+      - name: install mpfdtype
+        working-directory: mpfdtype
+        run: |
+          sudo apt install libmpfr-dev -y
+          CFLAGS="-Werror" python -m pip install . --no-build-isolation
+      - name: Run mpfdtype tests
+        working-directory: mpfdtype
+        run: |
+          pytest -vvv --color=yes
       - name: Install unytdtype
         working-directory: unytdtype
         run: |

mpfdtype/mpfdtype/src/casts.cpp

@@ -182,8 +182,9 @@ numpy_to_mpf_strided_loop(PyArrayMethod_Context *context,
         T *in = (T *)in_ptr;
         mpf_load(out, out_ptr, prec_out);
 
-        int res = C_to_mpfr<conv_T, T>(*in, out);
-        // TODO: At least for ints, could flag out of bounds...
+        C_to_mpfr<conv_T, T>(*in, out);
+        // TODO: At least for ints, could flag out of bounds, the return value
+        //       of C_to_mpfr may help with that (it flags imprecisions).
         mpf_store(out_ptr, out);
 
         in_ptr += strides[0];

mpfdtype/mpfdtype/src/dtype.c

@@ -46,7 +46,7 @@ new_MPFDType_instance(mpfr_prec_t precision)
                 "storage of single float would be too large for precision.");
     }
     new->base.elsize = sizeof(mpf_storage) + size;
-    new->base.alignment = _Alignof(mpf_storage);
+    new->base.alignment = _Alignof(mpf_field);
     new->base.flags |= NPY_NEEDS_INIT;
 
     return new;
@@ -213,6 +213,22 @@ MPFDType_repr(MPFDTypeObject *self)
 }
 
 
+PyObject *
+MPFDType_get_prec(MPFDTypeObject *self)
+{
+    return PyLong_FromLong(self->precision);
+}
+
+
+NPY_NO_EXPORT PyGetSetDef mpfdtype_getsetlist[] = {
+    {"prec",
+        (getter)MPFDType_get_prec,
+        NULL,
+        NULL, NULL},
+    {NULL, NULL, NULL, NULL, NULL},  /* Sentinel */
+};
+
+
 /*
  * This is the basic things that you need to create a Python Type/Class in C.
  * However, there is a slight difference here because we create a
@@ -226,6 +242,7 @@ PyArray_DTypeMeta MPFDType = {{{
         .tp_new = MPFDType_new,
         .tp_repr = (reprfunc)MPFDType_repr,
         .tp_str = (reprfunc)MPFDType_repr,
+        .tp_getset = mpfdtype_getsetlist,
     }},
     /* rest, filled in during DTypeMeta initialization */
 };

mpfdtype/mpfdtype/src/dtype.h

@@ -23,8 +23,6 @@ typedef struct {
  * So not storeing only those (saving 16 bytes of 48 for a 128 bit number)
  * removes the need to worry about this.
  */
-static_assert(_Alignof(mpfr_t) >= _Alignof(mp_limb_t),
-              "mpfr_t storage not aligned as much as limb_t?!");
 typedef mpfr_t mpf_storage;
 
 extern PyArray_DTypeMeta MPFDType;

mpfdtype/mpfdtype/src/scalar.c

@@ -184,13 +184,22 @@ MPFloat_repr(MPFloatObject* self)
 }
 
 
-int
-init_mpf_scalar(void)
+static PyObject *
+MPFloat_get_prec(MPFloatObject *self)
 {
-    return PyType_Ready(&MPFloat_Type);
+    return PyLong_FromLong(mpfr_get_prec(self->mpf.x));
 }
 
 
+NPY_NO_EXPORT PyGetSetDef mpfloat_getsetlist[] = {
+    {"prec",
+        (getter)MPFloat_get_prec,
+        NULL,
+        NULL, NULL},
+    {NULL, NULL, NULL, NULL, NULL},  /* Sentinel */
+};
+
+
 PyTypeObject MPFloat_Type = {
     PyVarObject_HEAD_INIT(NULL, 0)
     .tp_name = "MPFDType.MPFDType",
@@ -200,5 +209,13 @@ PyTypeObject MPFloat_Type = {
     .tp_repr = (reprfunc)MPFloat_repr,
     .tp_str = (reprfunc)MPFloat_str,
     .tp_as_number = &mpf_as_number,
-    .tp_richcompare = &mpf_richcompare,
+    .tp_richcompare = (richcmpfunc)mpf_richcompare,
+    .tp_getset = mpfloat_getsetlist,
 };
+
+
+int
+init_mpf_scalar(void)
+{
+    return PyType_Ready(&MPFloat_Type);
+}

mpfdtype/mpfdtype/src/terrible_hacks.c

@@ -22,9 +22,10 @@
  * (for larger itemsizes at least).
  */
 static void
-copyswap_mpf(char *dst, char *src, int swap, PyArrayObject *arr)
+copyswap_mpf(char *dst, char *src, int swap, PyArrayObject *ap)
 {
-    PyArray_Descr *descr = PyArray_DESCR(arr);
+    /* Note that it is probably better to only get the descr from `ap` */
+    PyArray_Descr *descr = PyArray_DESCR(ap);
 
     memcpy(dst, src, descr->elsize);
 }
@@ -33,9 +34,10 @@ copyswap_mpf(char *dst, char *src, int swap, PyArrayObject *arr)
 /* Should only be used for sorting, so more complex than necessary, probably */
 int compare_mpf(char *in1_ptr, char *in2_ptr, int swap, PyArrayObject *ap)
 {
+    /* Note that it is probably better to only get the descr from `ap` */
     mpfr_prec_t precision = ((MPFDTypeObject *)PyArray_DESCR(ap))->precision;
 
-    mpfr_t in1, in2;
+    mpfr_ptr in1, in2;
 
     mpf_load(in1, in1_ptr, precision);
     mpf_load(in2, in2_ptr, precision);
@@ -63,8 +65,8 @@ init_terrible_hacks(void) {
         return -1;
     }
     /* ->f slots are the same for all instances (currently). */
-    descr->base.f->copyswap = &copyswap_mpf;
-    descr->base.f->compare = &compare_mpf;
+    descr->base.f->copyswap = (PyArray_CopySwapFunc *)&copyswap_mpf;
+    descr->base.f->compare = (PyArray_CompareFunc *)&compare_mpf;
     Py_DECREF(descr);
 
     return 0;

mpfdtype/mpfdtype/tests/conftest.py

@@ -0,0 +1,3 @@
+import os
+
+os.environ["NUMPY_EXPERIMENTAL_DTYPE_API"] = "1"

mpfdtype/mpfdtype/tests/test_scalar.py

@@ -0,0 +1,94 @@
+import pytest
+
+import sys
+import numpy as np
+import operator
+
+from mpfdtype import MPFDType, MPFloat
+
+
+def test_create_scalar_simple():
+    # currently inferring 53bit precision from float:
+    assert MPFloat(12.).prec == 53
+    # currently infers 64bit or 32bit depending on system:
+    assert MPFloat(1).prec == sys.maxsize.bit_count() + 1
+
+    assert MPFloat(MPFloat(12.)).prec == 53
+    assert MPFloat(MPFloat(1)).prec == sys.maxsize.bit_count() + 1
+
+
+def test_create_scalar_prec():
+    assert MPFloat(1, prec=100).prec == 100
+    assert MPFloat(12., prec=123).prec == 123
+    assert MPFloat("12.234", prec=1000).prec == 1000
+
+    mpf1 = MPFloat("12.4325", prec=120)
+    mpf2 = MPFloat(mpf1, prec=150)
+    assert mpf1 == mpf2
+    assert mpf2.prec == 150
+
+
+def test_basic_equality():
+    assert MPFloat(12) == MPFloat(12.) == MPFloat("12.00", prec=10)
+
+
+@pytest.mark.parametrize("val", [123532.543, 12893283.5])
+def test_scalar_repr(val):
+    # For non exponentials at least, the repr matches:
+    val_repr = f"{val:e}".upper()
+    expected = f"MPFloat('{val_repr}', prec=20)"
+    assert repr(MPFloat(val, prec=20)) == expected
+
+@pytest.mark.parametrize("op",
+        ["add", "sub", "mul", "pow"])
+@pytest.mark.parametrize("other", [3., 12.5, 100., np.nan, np.inf])
+def test_binary_ops(op, other):
+    # Generally, the math ops should behave the same as double math if they
+    # use double precision (which they currently do).
+    # (double could have errors, but not for these simple ops)
+    op = getattr(operator, op)
+    try:
+        expected = op(12.5, other)
+    except Exception as e:
+        with pytest.raises(type(e)):
+            op(MPFloat(12.5), other)
+        with pytest.raises(type(e)):
+            op(12.5, MPFloat(other))
+        with pytest.raises(type(e)):
+            op(MPFloat(12.5), MPFloat(other))
+    else:
+        if np.isnan(expected):
+            # Avoiding isnan (which was also not implemented when written)
+            res = op(MPFloat(12.5), other)
+            assert res != res
+            res = op(12.5, MPFloat(other))
+            assert res != res
+            res = op(MPFloat(12.5), MPFloat(other))
+            assert res != res
+        else:
+            assert op(MPFloat(12.5), other) == expected
+            assert op(12.5, MPFloat(other)) == expected
+            assert op(MPFloat(12.5), MPFloat(other)) == expected
+
+
+@pytest.mark.parametrize("op",
+        ["eq", "ne", "le", "lt", "ge", "gt"])
+@pytest.mark.parametrize("other", [3., 12.5, 100., np.nan, np.inf])
+def test_comparisons(op, other):
+    op = getattr(operator, op)
+    expected = op(12.5, other)
+    assert op(MPFloat(12.5), other) is expected
+    assert op(12.5, MPFloat(other)) is expected
+    assert op(MPFloat(12.5), MPFloat(other)) is expected
+
+
+@pytest.mark.parametrize("op",
+        ["neg", "pos", "abs"])
+@pytest.mark.parametrize("val", [3., 12.5, 100., np.nan, np.inf])
+def test_comparisons(op, val):
+    op = getattr(operator, op)
+    expected = op(val)
+    if np.isnan(expected):
+        assert op(MPFloat(val)) != op(MPFloat(val))
+    else:
+        assert op(MPFloat(val)) == expected