From 4e497148b0610f86bbc0feff0da3e9627d7b691b Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Sun, 20 Jul 2025 18:49:59 +0000
Subject: [PATCH 1/7] Implement sign ufunc extend unary op tests

---
 quaddtype/numpy_quaddtype/src/ops.hpp   | 18 ++++++++++
 quaddtype/numpy_quaddtype/src/umath.cpp |  3 ++
 quaddtype/tests/test_quaddtype.py       | 48 ++++++++++++-------------
 3 files changed, 45 insertions(+), 24 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index f27f14b..62af903 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -17,6 +17,14 @@ quad_positive(Sleef_quad *op)
     return *op;
 }
 
+static inline Sleef_quad
+quad_sign(Sleef_quad *op)
+{
+    int32_t sign = Sleef_icmpq1(*op, Sleef_cast_from_doubleq1(0.0));
+    // sign(x=NaN) = x; otherwise sign(x) in { -1.0; 0.0; +1.0 }
+    return Sleef_iunordq1(*op, *op) ? *op : Sleef_cast_from_int64q1(sign);
+}
+
 static inline Sleef_quad
 quad_absolute(Sleef_quad *op)
 {
@@ -152,6 +160,16 @@ ld_absolute(long double *op)
     return fabsl(*op);
 }
 
+static inline long double
+ld_sign(long double *op)
+{
+    if (x < 0.0) return -1.0;
+    if (x == 0.0) return 0.0;
+    if (x > 0.0) return 1.0;
+    // sign(x=NaN) = x
+    return x;
+}
+
 static inline long double
 ld_rint(long double *op)
 {
diff --git a/quaddtype/numpy_quaddtype/src/umath.cpp b/quaddtype/numpy_quaddtype/src/umath.cpp
index 0058236..006e3d5 100644
--- a/quaddtype/numpy_quaddtype/src/umath.cpp
+++ b/quaddtype/numpy_quaddtype/src/umath.cpp
@@ -211,6 +211,9 @@ init_quad_unary_ops(PyObject *numpy)
     if (create_quad_unary_ufunc<quad_absolute, ld_absolute>(numpy, "absolute") < 0) {
         return -1;
     }
+    if (create_quad_unary_ufunc<quad_sign, ld_sign>(numpy, "sign") < 0) {
+        return -1;
+    }
     if (create_quad_unary_ufunc<quad_rint, ld_rint>(numpy, "rint") < 0) {
         return -1;
     }
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index 6bf4af5..c615abc 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -32,6 +32,7 @@ def test_binary_ops(op, other):
     quad_result = op_func(quad_a, quad_b)
     float_result = op_func(float_a, float_b)
 
+    # FIXME: @juntyr: replace with array_equal once isnan is supported
     with np.errstate(invalid="ignore"):
         assert (
             (np.float64(quad_result) == float_result) or
@@ -106,31 +107,30 @@ def test_array_aminmax(op, a, b):
         assert np.all((quad_res == float_res) | ((quad_res != quad_res) & (float_res != float_res)))
 
 
-@pytest.mark.parametrize("op, val, expected", [
-    ("neg", "3.0", "-3.0"),
-    ("neg", "-3.0", "3.0"),
-    ("pos", "3.0", "3.0"),
-    ("pos", "-3.0", "-3.0"),
-    ("abs", "3.0", "3.0"),
-    ("abs", "-3.0", "3.0"),
-    ("neg", "12.5", "-12.5"),
-    ("pos", "100.0", "100.0"),
-    ("abs", "-25.5", "25.5"),
-])
-def test_unary_ops(op, val, expected):
+@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign")])
+@pytest.mark.parametrize("val", ["3.0", "-3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
+def test_unary_ops(op, nop, val):
+    op_func = None if op is None else getattr(operator, op)
+    nop_func = getattr(np, nop)
+    
     quad_val = QuadPrecision(val)
-    expected_val = QuadPrecision(expected)
-
-    if op == "neg":
-        result = -quad_val
-    elif op == "pos":
-        result = +quad_val
-    elif op == "abs":
-        result = abs(quad_val)
-    else:
-        raise ValueError(f"Unsupported operation: {op}")
-
-    assert result == expected_val, f"{op}({val}) should be {expected}, but got {result}"
+    float_val = float(val)
+
+    for op_func in [op_func, nop_func]:
+        if op_func is None:
+            continue
+
+        quad_result = op_func(quad_val)
+        float_result = op_func(float_val)
+
+        # FIXME: @juntyr: replace with array_equal once isnan is supported
+        # FIXME: @juntyr: also check the signbit once that is supported
+        with np.errstate(invalid="ignore"):
+            assert (
+                (np.float64(quad_result) == float_result) or
+                ((float_result != float_result) and (quad_result != quad_result))
+            ), f"{op}({val}) should be {float_result}, but got {quad_result}"
+
 
 def test_inf():
     assert QuadPrecision("inf") > QuadPrecision("1e1000")

From ac85722bd86f7b9e037489c808e5030aefaa7f73 Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Sun, 20 Jul 2025 18:54:03 +0000
Subject: [PATCH 2/7] Add the copysign ufunc

---
 quaddtype/numpy_quaddtype/src/ops.hpp   | 12 ++++++++++++
 quaddtype/numpy_quaddtype/src/umath.cpp |  3 +++
 quaddtype/tests/test_quaddtype.py       |  4 ++--
 3 files changed, 17 insertions(+), 2 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index 62af903..0412f0a 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -339,6 +339,12 @@ quad_atan2(Sleef_quad *in1, Sleef_quad *in2)
     return Sleef_atan2q1_u10(*in1, *in2);
 }
 
+static inline Sleef_quad
+quad_copysign(Sleef_quad *in1, Sleef_quad *in2)
+{
+    return Sleef_copysignq1(*in1, *in2);
+}
+
 // Binary long double operations
 typedef long double (*binary_op_longdouble_def)(long double *, long double *);
 
@@ -396,6 +402,12 @@ ld_atan2(long double *in1, long double *in2)
     return atan2l(*in1, *in2);
 }
 
+static inline long double
+ld_copysign(long double *in1, long double *in2)
+{
+    return copysignl(*in1, *in2);
+}
+
 // comparison quad functions
 typedef npy_bool (*cmp_quad_def)(const Sleef_quad *, const Sleef_quad *);
 
diff --git a/quaddtype/numpy_quaddtype/src/umath.cpp b/quaddtype/numpy_quaddtype/src/umath.cpp
index 006e3d5..896b38d 100644
--- a/quaddtype/numpy_quaddtype/src/umath.cpp
+++ b/quaddtype/numpy_quaddtype/src/umath.cpp
@@ -556,6 +556,9 @@ init_quad_binary_ops(PyObject *numpy)
     if (create_quad_binary_ufunc<quad_atan2, ld_atan2>(numpy, "arctan2") < 0) {
         return -1;
     }
+    if (create_quad_binary_ufunc<quad_copysign, ld_copysign>(numpy, "copysign") < 0) {
+        return -1;
+    }
     return 0;
 }
 
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index c615abc..3c3b6aa 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -17,13 +17,13 @@ def test_basic_equality():
         "12.0") == QuadPrecision("12.00")
 
 
-@pytest.mark.parametrize("op", ["add", "sub", "mul", "truediv", "pow"])
+@pytest.mark.parametrize("op", ["add", "sub", "mul", "truediv", "pow", "copysign"])
 @pytest.mark.parametrize("other", ["3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
 def test_binary_ops(op, other):
     if op == "truediv" and float(other) == 0:
         pytest.xfail("float division by zero")
 
-    op_func = getattr(operator, op)
+    op_func = getattr(operator, op, None) or getattr(np, op)
     quad_a = QuadPrecision("12.5")
     quad_b = QuadPrecision(other)
     float_a = 12.5

From 2360f75fcc6804f99e0d51f23cab6b77fb709663 Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Sun, 20 Jul 2025 18:56:52 +0000
Subject: [PATCH 3/7] Fix typo

---
 quaddtype/numpy_quaddtype/src/ops.hpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index 0412f0a..f708221 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -163,11 +163,11 @@ ld_absolute(long double *op)
 static inline long double
 ld_sign(long double *op)
 {
-    if (x < 0.0) return -1.0;
-    if (x == 0.0) return 0.0;
-    if (x > 0.0) return 1.0;
+    if (*op < 0.0) return -1.0;
+    if (*op == 0.0) return 0.0;
+    if (*op > 0.0) return 1.0;
     // sign(x=NaN) = x
-    return x;
+    return *op;
 }
 
 static inline long double

From 021c770f398634426c2bb8e1ec574d0cb83efb92 Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Mon, 21 Jul 2025 05:44:38 +0000
Subject: [PATCH 4/7] Implement signbit ufunc

---
 quaddtype/numpy_quaddtype/src/ops.hpp   | 21 ++++++++++++++++++++-
 quaddtype/numpy_quaddtype/src/umath.cpp |  3 +++
 quaddtype/tests/test_quaddtype.py       |  5 +++--
 3 files changed, 26 insertions(+), 3 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index f708221..e0cd7e3 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -20,11 +20,24 @@ quad_positive(Sleef_quad *op)
 static inline Sleef_quad
 quad_sign(Sleef_quad *op)
 {
-    int32_t sign = Sleef_icmpq1(*op, Sleef_cast_from_doubleq1(0.0));
+    Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
+    int32_t sign = Sleef_icmpq1(*op, zero);
     // sign(x=NaN) = x; otherwise sign(x) in { -1.0; 0.0; +1.0 }
     return Sleef_iunordq1(*op, *op) ? *op : Sleef_cast_from_int64q1(sign);
 }
 
+static inline int32_t
+quad_signbit(Sleef_quad *op)
+{
+    // FIXME @juntyr or @SwayamInSync: replace with binary implementation
+    //  once we test big and little endian in CI
+    Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
+    Sleef_quad one = Sleef_cast_from_doubleq1(1.0);
+    Sleef_quad one_signed = Sleef_copysignq1(one, *op);
+    // signbit(x) = 1 iff copysign(1, x) == -1
+    return Sleef_icmpltq1(one_signed, zero);
+}
+
 static inline Sleef_quad
 quad_absolute(Sleef_quad *op)
 {
@@ -170,6 +183,12 @@ ld_sign(long double *op)
     return *op;
 }
 
+static inline int32_t
+ld_signbit(long double *op)
+{
+    return signbit(*op);
+}
+
 static inline long double
 ld_rint(long double *op)
 {
diff --git a/quaddtype/numpy_quaddtype/src/umath.cpp b/quaddtype/numpy_quaddtype/src/umath.cpp
index 896b38d..9a79b59 100644
--- a/quaddtype/numpy_quaddtype/src/umath.cpp
+++ b/quaddtype/numpy_quaddtype/src/umath.cpp
@@ -214,6 +214,9 @@ init_quad_unary_ops(PyObject *numpy)
     if (create_quad_unary_ufunc<quad_sign, ld_sign>(numpy, "sign") < 0) {
         return -1;
     }
+    if (create_quad_unary_ufunc<quad_signbit, ld_signbit>(numpy, "signbit") < 0) {
+        return -1;
+    }
     if (create_quad_unary_ufunc<quad_rint, ld_rint>(numpy, "rint") < 0) {
         return -1;
     }
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index 3c3b6aa..d0fdc9c 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -107,7 +107,7 @@ def test_array_aminmax(op, a, b):
         assert np.all((quad_res == float_res) | ((quad_res != quad_res) & (float_res != float_res)))
 
 
-@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign")])
+@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign"), (None, "signbit")])
 @pytest.mark.parametrize("val", ["3.0", "-3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
 def test_unary_ops(op, nop, val):
     op_func = None if op is None else getattr(operator, op)
@@ -124,11 +124,12 @@ def test_unary_ops(op, nop, val):
         float_result = op_func(float_val)
 
         # FIXME: @juntyr: replace with array_equal once isnan is supported
-        # FIXME: @juntyr: also check the signbit once that is supported
         with np.errstate(invalid="ignore"):
             assert (
                 (np.float64(quad_result) == float_result) or
                 ((float_result != float_result) and (quad_result != quad_result))
+            ) and (
+                np.signbit(float_result) == np.signbit(quad_result)
             ), f"{op}({val}) should be {float_result}, but got {quad_result}"
 
 

From 6ffcdca10ddc9a5716ca17964e7faa38423e0495 Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Mon, 21 Jul 2025 06:17:45 +0000
Subject: [PATCH 5/7] Add support for unary ufunc properties

---
 quaddtype/numpy_quaddtype/src/ops.hpp   | 180 ++++++++++++------------
 quaddtype/numpy_quaddtype/src/umath.cpp | 130 ++++++++++++++++-
 2 files changed, 220 insertions(+), 90 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index e0cd7e3..50262bf 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -3,22 +3,22 @@
 #include <cmath>
 
 // Unary Quad Operations
-typedef Sleef_quad (*unary_op_quad_def)(Sleef_quad *);
+typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
 
 static Sleef_quad
-quad_negative(Sleef_quad *op)
+quad_negative(const Sleef_quad *op)
 {
     return Sleef_negq1(*op);
 }
 
 static Sleef_quad
-quad_positive(Sleef_quad *op)
+quad_positive(const Sleef_quad *op)
 {
     return *op;
 }
 
 static inline Sleef_quad
-quad_sign(Sleef_quad *op)
+quad_sign(const Sleef_quad *op)
 {
     Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
     int32_t sign = Sleef_icmpq1(*op, zero);
@@ -26,155 +26,143 @@ quad_sign(Sleef_quad *op)
     return Sleef_iunordq1(*op, *op) ? *op : Sleef_cast_from_int64q1(sign);
 }
 
-static inline int32_t
-quad_signbit(Sleef_quad *op)
-{
-    // FIXME @juntyr or @SwayamInSync: replace with binary implementation
-    //  once we test big and little endian in CI
-    Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
-    Sleef_quad one = Sleef_cast_from_doubleq1(1.0);
-    Sleef_quad one_signed = Sleef_copysignq1(one, *op);
-    // signbit(x) = 1 iff copysign(1, x) == -1
-    return Sleef_icmpltq1(one_signed, zero);
-}
-
 static inline Sleef_quad
-quad_absolute(Sleef_quad *op)
+quad_absolute(const Sleef_quad *op)
 {
     return Sleef_fabsq1(*op);
 }
 
 static inline Sleef_quad
-quad_rint(Sleef_quad *op)
+quad_rint(const Sleef_quad *op)
 {
     return Sleef_rintq1(*op);
 }
 
 static inline Sleef_quad
-quad_trunc(Sleef_quad *op)
+quad_trunc(const Sleef_quad *op)
 {
     return Sleef_truncq1(*op);
 }
 
 static inline Sleef_quad
-quad_floor(Sleef_quad *op)
+quad_floor(const Sleef_quad *op)
 {
     return Sleef_floorq1(*op);
 }
 
 static inline Sleef_quad
-quad_ceil(Sleef_quad *op)
+quad_ceil(const Sleef_quad *op)
 {
     return Sleef_ceilq1(*op);
 }
 
 static inline Sleef_quad
-quad_sqrt(Sleef_quad *op)
+quad_sqrt(const Sleef_quad *op)
 {
     return Sleef_sqrtq1_u05(*op);
 }
 
 static inline Sleef_quad
-quad_square(Sleef_quad *op)
+quad_square(const Sleef_quad *op)
 {
     return Sleef_mulq1_u05(*op, *op);
 }
 
 static inline Sleef_quad
-quad_log(Sleef_quad *op)
+quad_log(const Sleef_quad *op)
 {
     return Sleef_logq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_log2(Sleef_quad *op)
+quad_log2(const Sleef_quad *op)
 {
     return Sleef_log2q1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_log10(Sleef_quad *op)
+quad_log10(const Sleef_quad *op)
 {
     return Sleef_log10q1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_log1p(Sleef_quad *op)
+quad_log1p(const Sleef_quad *op)
 {
     return Sleef_log1pq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_exp(Sleef_quad *op)
+quad_exp(const Sleef_quad *op)
 {
     return Sleef_expq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_exp2(Sleef_quad *op)
+quad_exp2(const Sleef_quad *op)
 {
     return Sleef_exp2q1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_sin(Sleef_quad *op)
+quad_sin(const Sleef_quad *op)
 {
     return Sleef_sinq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_cos(Sleef_quad *op)
+quad_cos(const Sleef_quad *op)
 {
     return Sleef_cosq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_tan(Sleef_quad *op)
+quad_tan(const Sleef_quad *op)
 {
     return Sleef_tanq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_asin(Sleef_quad *op)
+quad_asin(const Sleef_quad *op)
 {
     return Sleef_asinq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_acos(Sleef_quad *op)
+quad_acos(const Sleef_quad *op)
 {
     return Sleef_acosq1_u10(*op);
 }
 
 static inline Sleef_quad
-quad_atan(Sleef_quad *op)
+quad_atan(const Sleef_quad *op)
 {
     return Sleef_atanq1_u10(*op);
 }
 
 // Unary long double operations
-typedef long double (*unary_op_longdouble_def)(long double *);
+typedef long double (*unary_op_longdouble_def)(const long double *);
 
 static inline long double
-ld_negative(long double *op)
+ld_negative(const long double *op)
 {
     return -(*op);
 }
 
 static inline long double
-ld_positive(long double *op)
+ld_positive(const long double *op)
 {
     return *op;
 }
 
 static inline long double
-ld_absolute(long double *op)
+ld_absolute(const long double *op)
 {
     return fabsl(*op);
 }
 
 static inline long double
-ld_sign(long double *op)
+ld_sign(const long double *op)
 {
     if (*op < 0.0) return -1.0;
     if (*op == 0.0) return 0.0;
@@ -183,161 +171,179 @@ ld_sign(long double *op)
     return *op;
 }
 
-static inline int32_t
-ld_signbit(long double *op)
-{
-    return signbit(*op);
-}
-
 static inline long double
-ld_rint(long double *op)
+ld_rint(const long double *op)
 {
     return rintl(*op);
 }
 
 static inline long double
-ld_trunc(long double *op)
+ld_trunc(const long double *op)
 {
     return truncl(*op);
 }
 
 static inline long double
-ld_floor(long double *op)
+ld_floor(const long double *op)
 {
     return floorl(*op);
 }
 
 static inline long double
-ld_ceil(long double *op)
+ld_ceil(const long double *op)
 {
     return ceill(*op);
 }
 
 static inline long double
-ld_sqrt(long double *op)
+ld_sqrt(const long double *op)
 {
     return sqrtl(*op);
 }
 
 static inline long double
-ld_square(long double *op)
+ld_square(const long double *op)
 {
     return (*op) * (*op);
 }
 
 static inline long double
-ld_log(long double *op)
+ld_log(const long double *op)
 {
     return logl(*op);
 }
 
 static inline long double
-ld_log2(long double *op)
+ld_log2(const long double *op)
 {
     return log2l(*op);
 }
 
 static inline long double
-ld_log10(long double *op)
+ld_log10(const long double *op)
 {
     return log10l(*op);
 }
 
 static inline long double
-ld_log1p(long double *op)
+ld_log1p(const long double *op)
 {
     return log1pl(*op);
 }
 
 static inline long double
-ld_exp(long double *op)
+ld_exp(const long double *op)
 {
     return expl(*op);
 }
 
 static inline long double
-ld_exp2(long double *op)
+ld_exp2(const long double *op)
 {
     return exp2l(*op);
 }
 
 static inline long double
-ld_sin(long double *op)
+ld_sin(const long double *op)
 {
     return sinl(*op);
 }
 
 static inline long double
-ld_cos(long double *op)
+ld_cos(const long double *op)
 {
     return cosl(*op);
 }
 
 static inline long double
-ld_tan(long double *op)
+ld_tan(const long double *op)
 {
     return tanl(*op);
 }
 
 static inline long double
-ld_asin(long double *op)
+ld_asin(const long double *op)
 {
     return asinl(*op);
 }
 
 static inline long double
-ld_acos(long double *op)
+ld_acos(const long double *op)
 {
     return acosl(*op);
 }
 
 static inline long double
-ld_atan(long double *op)
+ld_atan(const long double *op)
 {
     return atanl(*op);
 }
 
+// Unary Quad properties
+typedef npy_bool (*unary_prop_quad_def)(const Sleef_quad *);
+
+static inline npy_bool
+quad_signbit(const Sleef_quad *op)
+{
+    // FIXME @juntyr or @SwayamInSync: replace with binary implementation
+    //  once we test big and little endian in CI
+    Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
+    Sleef_quad one = Sleef_cast_from_doubleq1(1.0);
+    Sleef_quad one_signed = Sleef_copysignq1(one, *op);
+    // signbit(x) = 1 iff copysign(1, x) == -1
+    return Sleef_icmpltq1(one_signed, zero);
+}
+
+// Unary Quad properties
+typedef npy_bool (*unary_prop_longdouble_def)(const long double *);
+
+static inline npy_bool
+ld_signbit(const long double *op)
+{
+    return signbit(*op);
+}
+
 // Binary Quad operations
-typedef Sleef_quad (*binary_op_quad_def)(Sleef_quad *, Sleef_quad *);
+typedef Sleef_quad (*binary_op_quad_def)(const Sleef_quad *, const Sleef_quad *);
 
 static inline Sleef_quad
-quad_add(Sleef_quad *in1, Sleef_quad *in2)
+quad_add(const Sleef_quad *in1, const Sleef_quad *in2)
 {
     return Sleef_addq1_u05(*in1, *in2);
 }
 
 static inline Sleef_quad
-quad_sub(Sleef_quad *in1, Sleef_quad *in2)
+quad_sub(const Sleef_quad *in1, const Sleef_quad *in2)
 {
     return Sleef_subq1_u05(*in1, *in2);
 }
 
 static inline Sleef_quad
-quad_mul(Sleef_quad *a, Sleef_quad *b)
+quad_mul(const Sleef_quad *a, const Sleef_quad *b)
 {
     return Sleef_mulq1_u05(*a, *b);
 }
 
 static inline Sleef_quad
-quad_div(Sleef_quad *a, Sleef_quad *b)
+quad_div(const Sleef_quad *a, const Sleef_quad *b)
 {
     return Sleef_divq1_u05(*a, *b);
 }
 
 static inline Sleef_quad
-quad_pow(Sleef_quad *a, Sleef_quad *b)
+quad_pow(const Sleef_quad *a, const Sleef_quad *b)
 {
     return Sleef_powq1_u10(*a, *b);
 }
 
 static inline Sleef_quad
-quad_mod(Sleef_quad *a, Sleef_quad *b)
+quad_mod(const Sleef_quad *a, const Sleef_quad *b)
 {
     return Sleef_fmodq1(*a, *b);
 }
 
 static inline Sleef_quad
-quad_minimum(Sleef_quad *in1, Sleef_quad *in2)
+quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2)
 {
     return Sleef_iunordq1(*in1, *in2) ? (
         Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
@@ -345,7 +351,7 @@ quad_minimum(Sleef_quad *in1, Sleef_quad *in2)
 }
 
 static inline Sleef_quad
-quad_maximum(Sleef_quad *in1, Sleef_quad *in2)
+quad_maximum(const Sleef_quad *in1, const Sleef_quad *in2)
 {
     return Sleef_iunordq1(*in1, *in2) ? (
         Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
@@ -353,76 +359,76 @@ quad_maximum(Sleef_quad *in1, Sleef_quad *in2)
 }
 
 static inline Sleef_quad
-quad_atan2(Sleef_quad *in1, Sleef_quad *in2)
+quad_atan2(const Sleef_quad *in1, const Sleef_quad *in2)
 {
     return Sleef_atan2q1_u10(*in1, *in2);
 }
 
 static inline Sleef_quad
-quad_copysign(Sleef_quad *in1, Sleef_quad *in2)
+quad_copysign(const Sleef_quad *in1, const Sleef_quad *in2)
 {
     return Sleef_copysignq1(*in1, *in2);
 }
 
 // Binary long double operations
-typedef long double (*binary_op_longdouble_def)(long double *, long double *);
+typedef long double (*binary_op_longdouble_def)(const long double *, const long double *);
 
 static inline long double
-ld_add(long double *in1, long double *in2)
+ld_add(const long double *in1, const long double *in2)
 {
     return (*in1) + (*in2);
 }
 
 static inline long double
-ld_sub(long double *in1, long double *in2)
+ld_sub(const long double *in1, const long double *in2)
 {
     return (*in1) - (*in2);
 }
 
 static inline long double
-ld_mul(long double *a, long double *b)
+ld_mul(const long double *a, const long double *b)
 {
     return (*a) * (*b);
 }
 
 static inline long double
-ld_div(long double *a, long double *b)
+ld_div(const long double *a, const long double *b)
 {
     return (*a) / (*b);
 }
 
 static inline long double
-ld_pow(long double *a, long double *b)
+ld_pow(const long double *a, const long double *b)
 {
     return powl(*a, *b);
 }
 
 static inline long double
-ld_mod(long double *a, long double *b)
+ld_mod(const long double *a, const long double *b)
 {
     return fmodl(*a, *b);
 }
 
 static inline long double
-ld_minimum(long double *in1, long double *in2)
+ld_minimum(const long double *in1, const long double *in2)
 {
     return (*in1 < *in2) ? *in1 : *in2;
 }
 
 static inline long double
-ld_maximum(long double *in1, long double *in2)
+ld_maximum(const long double *in1, const long double *in2)
 {
     return (*in1 > *in2) ? *in1 : *in2;
 }
 
 static inline long double
-ld_atan2(long double *in1, long double *in2)
+ld_atan2(const long double *in1, const long double *in2)
 {
     return atan2l(*in1, *in2);
 }
 
 static inline long double
-ld_copysign(long double *in1, long double *in2)
+ld_copysign(const long double *in1, const long double *in2)
 {
     return copysignl(*in1, *in2);
 }
@@ -503,4 +509,4 @@ static inline npy_bool
 ld_greaterequal(const long double *a, const long double *b)
 {
     return *a >= *b;
-}
\ No newline at end of file
+}
diff --git a/quaddtype/numpy_quaddtype/src/umath.cpp b/quaddtype/numpy_quaddtype/src/umath.cpp
index 9a79b59..b23f395 100644
--- a/quaddtype/numpy_quaddtype/src/umath.cpp
+++ b/quaddtype/numpy_quaddtype/src/umath.cpp
@@ -214,9 +214,6 @@ init_quad_unary_ops(PyObject *numpy)
     if (create_quad_unary_ufunc<quad_sign, ld_sign>(numpy, "sign") < 0) {
         return -1;
     }
-    if (create_quad_unary_ufunc<quad_signbit, ld_signbit>(numpy, "signbit") < 0) {
-        return -1;
-    }
     if (create_quad_unary_ufunc<quad_rint, ld_rint>(numpy, "rint") < 0) {
         return -1;
     }
@@ -274,6 +271,128 @@ init_quad_unary_ops(PyObject *numpy)
     return 0;
 }
 
+static NPY_CASTING
+quad_unary_prop_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];
+
+    loop_descrs[1] = PyArray_DescrFromType(NPY_BOOL);
+
+    return NPY_NO_CASTING;
+}
+
+template <unary_prop_quad_def sleef_op, unary_prop_longdouble_def longdouble_op>
+int
+quad_generic_unary_prop_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 *out_ptr = data[1];
+    npy_intp in_stride = strides[0];
+    npy_intp out_stride = strides[1];
+
+    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;
+    while (N--) {
+        memcpy(&in, in_ptr, elem_size);
+        npy_bool result;
+        if (backend == BACKEND_SLEEF) {
+            result = sleef_op(&in.sleef_value);
+        }
+        else {
+            result = longdouble_op(&in.longdouble_value);
+        }
+        memcpy(out_ptr, &result, sizeof(npy_bool));
+
+        in_ptr += in_stride;
+        out_ptr += out_stride;
+    }
+    return 0;
+}
+
+template <unary_prop_quad_def sleef_op, unary_prop_longdouble_def longdouble_op>
+int
+quad_generic_unary_prop_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 *out_ptr = data[1];
+    npy_intp in_stride = strides[0];
+    npy_intp out_stride = strides[1];
+
+    QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0];
+    QuadBackendType backend = descr->backend;
+
+    while (N--) {
+        npy_bool result;
+        if (backend == BACKEND_SLEEF) {
+            result = sleef_op((Sleef_quad *)in_ptr);
+        }
+        else {
+            result = longdouble_op((long double *)in_ptr);
+        }
+        *(npy_bool *)out_ptr = result;
+        in_ptr += in_stride;
+        out_ptr += out_stride;
+    }
+    return 0;
+}
+
+template <unary_prop_quad_def sleef_op, unary_prop_longdouble_def longdouble_op>
+int
+create_quad_unary_prop_ufunc(PyObject *numpy, const char *ufunc_name)
+{
+    PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name);
+    if (ufunc == NULL) {
+        return -1;
+    }
+
+    PyArray_DTypeMeta *dtypes[2] = {&QuadPrecDType, &PyArray_BoolDType};
+
+    PyType_Slot slots[] = {
+            {NPY_METH_resolve_descriptors, (void *)&quad_unary_prop_resolve_descriptors},
+            {NPY_METH_strided_loop,
+             (void *)&quad_generic_unary_prop_strided_loop_aligned<sleef_op, longdouble_op>},
+            {NPY_METH_unaligned_strided_loop,
+             (void *)&quad_generic_unary_prop_strided_loop_unaligned<sleef_op, longdouble_op>},
+            {0, NULL}};
+
+    PyArrayMethod_Spec Spec = {
+            .name = "quad_unary_prop",
+            .nin = 1,
+            .nout = 1,
+            .casting = NPY_NO_CASTING,
+            .flags = NPY_METH_SUPPORTS_UNALIGNED,
+            .dtypes = dtypes,
+            .slots = slots,
+    };
+
+    if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) {
+        return -1;
+    }
+
+    return 0;
+}
+
+int
+init_quad_unary_props(PyObject *numpy)
+{
+    if (create_quad_unary_prop_ufunc<quad_signbit, ld_signbit>(numpy, "signbit") < 0) {
+        return -1;
+    }
+    return 0;
+}
+
 // Binary ufuncs
 
 static NPY_CASTING
@@ -796,6 +915,11 @@ init_quad_umath(void)
         goto err;
     }
 
+    if (init_quad_unary_props(numpy) < 0) {
+        PyErr_SetString(PyExc_RuntimeError, "Failed to initialize quad unary properties");
+        goto err;
+    }
+
     if (init_quad_binary_ops(numpy) < 0) {
         PyErr_SetString(PyExc_RuntimeError, "Failed to initialize quad binary operations");
         goto err;

From 93f7c0ad113386670fdce945e259df958d6ad91e Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Mon, 21 Jul 2025 12:32:50 +0000
Subject: [PATCH 6/7] Implement isfinite, isinf, isnan ufuncs

---
 quaddtype/numpy_quaddtype/src/ops.hpp   | 40 +++++++++++++++++++++++++
 quaddtype/numpy_quaddtype/src/umath.cpp |  9 ++++++
 quaddtype/tests/test_quaddtype.py       | 28 ++++-------------
 3 files changed, 55 insertions(+), 22 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index 50262bf..6d76c24 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -294,6 +294,28 @@ quad_signbit(const Sleef_quad *op)
     return Sleef_icmpltq1(one_signed, zero);
 }
 
+static inline npy_bool
+quad_isfinite(const Sleef_quad *op)
+{
+    // isfinite(x) = abs(x) < inf
+    Sleef_quad inf = Sleef_cast_from_doubleq1(std::numeric_limits<double>::infinity());
+    return Sleef_icmpltq1(Sleef_fabsq1(*op), inf);
+}
+
+static inline npy_bool
+quad_isinf(const Sleef_quad *op)
+{
+    // isinf(x) = abs(x) == inf
+    Sleef_quad inf = Sleef_cast_from_doubleq1(std::numeric_limits<double>::infinity());
+    return Sleef_icmpeqq1(Sleef_fabsq1(*op), inf);
+}
+
+static inline npy_bool
+quad_isnan(const Sleef_quad *op)
+{
+    return Sleef_iunordq1(*op, *op);
+}
+
 // Unary Quad properties
 typedef npy_bool (*unary_prop_longdouble_def)(const long double *);
 
@@ -303,6 +325,24 @@ ld_signbit(const long double *op)
     return signbit(*op);
 }
 
+static inline npy_bool
+ld_isfinite(const long double *op)
+{
+    return isfinite(*op);
+}
+
+static inline npy_bool
+ld_isinf(const long double *op)
+{
+    return isinf(*op);
+}
+
+static inline npy_bool
+ld_isnan(const long double *op)
+{
+    return isnan(*op);
+}
+
 // Binary Quad operations
 typedef Sleef_quad (*binary_op_quad_def)(const Sleef_quad *, const Sleef_quad *);
 
diff --git a/quaddtype/numpy_quaddtype/src/umath.cpp b/quaddtype/numpy_quaddtype/src/umath.cpp
index b23f395..7716692 100644
--- a/quaddtype/numpy_quaddtype/src/umath.cpp
+++ b/quaddtype/numpy_quaddtype/src/umath.cpp
@@ -387,6 +387,15 @@ create_quad_unary_prop_ufunc(PyObject *numpy, const char *ufunc_name)
 int
 init_quad_unary_props(PyObject *numpy)
 {
+    if (create_quad_unary_prop_ufunc<quad_isfinite, ld_isfinite>(numpy, "isfinite") < 0) {
+        return -1;
+    }
+    if (create_quad_unary_prop_ufunc<quad_isinf, ld_isinf>(numpy, "isinf") < 0) {
+        return -1;
+    }
+    if (create_quad_unary_prop_ufunc<quad_isnan, ld_isnan>(numpy, "isnan") < 0) {
+        return -1;
+    }
     if (create_quad_unary_prop_ufunc<quad_signbit, ld_signbit>(numpy, "signbit") < 0) {
         return -1;
     }
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index d0fdc9c..98b7a94 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -32,13 +32,7 @@ def test_binary_ops(op, other):
     quad_result = op_func(quad_a, quad_b)
     float_result = op_func(float_a, float_b)
 
-    # FIXME: @juntyr: replace with array_equal once isnan is supported
-    with np.errstate(invalid="ignore"):
-        assert (
-            (np.float64(quad_result) == float_result) or
-            (np.abs(np.float64(quad_result) - float_result) < 1e-10) or
-            ((float_result != float_result) and (quad_result != quad_result))
-        )
+    np.testing.assert_allclose(np.float64(quad_result), float_result, atol=1e-10, rtol=0, equal_nan=True)
 
 
 @pytest.mark.parametrize("op", ["eq", "ne", "le", "lt", "ge", "gt"])
@@ -83,9 +77,7 @@ def test_array_minmax(op, a, b):
     quad_res = op_func(quad_a, quad_b)
     float_res = op_func(float_a, float_b)
 
-    # FIXME: @juntyr: replace with array_equal once isnan is supported
-    with np.errstate(invalid="ignore"):
-        assert np.all((quad_res == float_res) | ((quad_res != quad_res) & (float_res != float_res)))
+    np.testing.assert_array_equal(quad_res.astype(float), float_res)
 
 
 @pytest.mark.parametrize("op", ["amin", "amax", "nanmin", "nanmax"])
@@ -102,12 +94,10 @@ def test_array_aminmax(op, a, b):
     quad_res = op_func(quad_ab)
     float_res = op_func(float_ab)
 
-    # FIXME: @juntyr: replace with array_equal once isnan is supported
-    with np.errstate(invalid="ignore"):
-        assert np.all((quad_res == float_res) | ((quad_res != quad_res) & (float_res != float_res)))
+    np.testing.assert_array_equal(np.array(quad_res).astype(float), float_res)
 
 
-@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign"), (None, "signbit")])
+@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign"), (None, "signbit"), (None, "isfinite"), (None, "isinf"), (None, "isnan")])
 @pytest.mark.parametrize("val", ["3.0", "-3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
 def test_unary_ops(op, nop, val):
     op_func = None if op is None else getattr(operator, op)
@@ -123,14 +113,8 @@ def test_unary_ops(op, nop, val):
         quad_result = op_func(quad_val)
         float_result = op_func(float_val)
 
-        # FIXME: @juntyr: replace with array_equal once isnan is supported
-        with np.errstate(invalid="ignore"):
-            assert (
-                (np.float64(quad_result) == float_result) or
-                ((float_result != float_result) and (quad_result != quad_result))
-            ) and (
-                np.signbit(float_result) == np.signbit(quad_result)
-            ), f"{op}({val}) should be {float_result}, but got {quad_result}"
+        np.testing.assert_array_equal(np.array(quad_result).astype(float), float_result)
+        assert np.signbit(float_result) == np.signbit(quad_result)
 
 
 def test_inf():

From d5a15f1c0689ea1ffaee6bb33e5d13d25fc016db Mon Sep 17 00:00:00 2001
From: Juniper Tyree <50025784+juntyr@users.noreply.github.com>
Date: Tue, 22 Jul 2025 06:55:20 +0000
Subject: [PATCH 7/7] Use sleef_q constants to avoid cast warnings

---
 quaddtype/numpy_quaddtype/src/ops.hpp | 21 +++++++++++----------
 quaddtype/tests/test_quaddtype.py     |  8 ++++----
 2 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index 6d76c24..9f5767a 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -1,6 +1,12 @@
 #include <sleef.h>
 #include <sleefquad.h>
 #include <cmath>
+#include <limits>
+
+// Quad Constants, generated with qutil
+#define QUAD_ZERO sleef_q(+0x0000000000000LL, 0x0000000000000000ULL, -16383)
+#define QUAD_ONE sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 0)
+#define QUAD_POS_INF sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384)
 
 // Unary Quad Operations
 typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
@@ -20,8 +26,7 @@ quad_positive(const Sleef_quad *op)
 static inline Sleef_quad
 quad_sign(const Sleef_quad *op)
 {
-    Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
-    int32_t sign = Sleef_icmpq1(*op, zero);
+    int32_t sign = Sleef_icmpq1(*op, QUAD_ZERO);
     // sign(x=NaN) = x; otherwise sign(x) in { -1.0; 0.0; +1.0 }
     return Sleef_iunordq1(*op, *op) ? *op : Sleef_cast_from_int64q1(sign);
 }
@@ -287,27 +292,23 @@ quad_signbit(const Sleef_quad *op)
 {
     // FIXME @juntyr or @SwayamInSync: replace with binary implementation
     //  once we test big and little endian in CI
-    Sleef_quad zero = Sleef_cast_from_doubleq1(0.0);
-    Sleef_quad one = Sleef_cast_from_doubleq1(1.0);
-    Sleef_quad one_signed = Sleef_copysignq1(one, *op);
+    Sleef_quad one_signed = Sleef_copysignq1(QUAD_ONE, *op);
     // signbit(x) = 1 iff copysign(1, x) == -1
-    return Sleef_icmpltq1(one_signed, zero);
+    return Sleef_icmpltq1(one_signed, QUAD_ZERO);
 }
 
 static inline npy_bool
 quad_isfinite(const Sleef_quad *op)
 {
     // isfinite(x) = abs(x) < inf
-    Sleef_quad inf = Sleef_cast_from_doubleq1(std::numeric_limits<double>::infinity());
-    return Sleef_icmpltq1(Sleef_fabsq1(*op), inf);
+    return Sleef_icmpltq1(Sleef_fabsq1(*op), QUAD_POS_INF);
 }
 
 static inline npy_bool
 quad_isinf(const Sleef_quad *op)
 {
     // isinf(x) = abs(x) == inf
-    Sleef_quad inf = Sleef_cast_from_doubleq1(std::numeric_limits<double>::infinity());
-    return Sleef_icmpeqq1(Sleef_fabsq1(*op), inf);
+    return Sleef_icmpeqq1(Sleef_fabsq1(*op), QUAD_POS_INF);
 }
 
 static inline npy_bool
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index 98b7a94..018fda1 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -97,11 +97,11 @@ def test_array_aminmax(op, a, b):
     np.testing.assert_array_equal(np.array(quad_res).astype(float), float_res)
 
 
-@pytest.mark.parametrize("op,nop", [("neg", "negative"), ("pos", "positive"), ("abs", "absolute"), (None, "sign"), (None, "signbit"), (None, "isfinite"), (None, "isinf"), (None, "isnan")])
+@pytest.mark.parametrize("op", ["negative", "positive", "absolute", "sign", "signbit", "isfinite", "isinf", "isnan"])
 @pytest.mark.parametrize("val", ["3.0", "-3.0", "12.5", "100.0", "0.0", "-0.0", "inf", "-inf", "nan", "-nan"])
-def test_unary_ops(op, nop, val):
-    op_func = None if op is None else getattr(operator, op)
-    nop_func = getattr(np, nop)
+def test_unary_ops(op, val):
+    op_func = dict(negative=operator.neg, positive=operator.pos, absolute=operator.abs).get(op, None)
+    nop_func = getattr(np, op)
     
     quad_val = QuadPrecision(val)
     float_val = float(val)