diff --git a/quaddtype/numpy_quaddtype/src/ops.hpp b/quaddtype/numpy_quaddtype/src/ops.hpp
index 1a45bac..ae38dbe 100644
--- a/quaddtype/numpy_quaddtype/src/ops.hpp
+++ b/quaddtype/numpy_quaddtype/src/ops.hpp
@@ -6,6 +6,7 @@
 #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)
+#define QUAD_NAN sleef_q(+0x1ffffffffffffLL, 0xffffffffffffffffULL, 16384)
 
 // Unary Quad Operations
 typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
@@ -174,9 +175,12 @@ ld_absolute(const long double *op)
 static inline long double
 ld_sign(const long double *op)
 {
-    if (*op < 0.0) return -1.0;
-    if (*op == 0.0) return 0.0;
-    if (*op > 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 *op;
 }
@@ -391,39 +395,80 @@ quad_pow(const Sleef_quad *a, const Sleef_quad *b)
 static inline Sleef_quad
 quad_mod(const Sleef_quad *a, const Sleef_quad *b)
 {
-    return Sleef_fmodq1(*a, *b);
+    // division by zero
+    if (Sleef_icmpeqq1(*b, QUAD_ZERO)) {
+        return QUAD_NAN;
+    }
+
+    // NaN inputs
+    if (Sleef_iunordq1(*a, *b)) {
+        return Sleef_iunordq1(*a, *a) ? *a : *b;  // Return the NaN
+    }
+
+    // infinity dividend -> NaN
+    if (quad_isinf(a)) {
+        return QUAD_NAN;
+    }
+
+    // finite % inf
+    if (quad_isfinite(a) && quad_isinf(b)) {
+        int sign_a = quad_signbit(a);
+        int sign_b = quad_signbit(b);
+
+        // return a if sign_a == sign_b
+        return (sign_a == sign_b) ? *a : *b;
+    }
+
+    // NumPy mod formula: a % b = a - floor(a/b) * b
+    Sleef_quad quotient = Sleef_divq1_u05(*a, *b);
+    Sleef_quad floored = Sleef_floorq1(quotient);
+    Sleef_quad product = Sleef_mulq1_u05(floored, *b);
+    Sleef_quad result = Sleef_subq1_u05(*a, product);
+
+    // Handle zero result sign: when result is exactly zero,
+    // it should have the same sign as the divisor (NumPy convention)
+    if (Sleef_icmpeqq1(result, QUAD_ZERO)) {
+        if (Sleef_icmpltq1(*b, QUAD_ZERO)) {
+            return Sleef_negq1(QUAD_ZERO);  // -0.0
+        }
+        else {
+            return QUAD_ZERO;  // +0.0
+        }
+    }
+
+    return result;
 }
 
 static inline Sleef_quad
 quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2)
 {
-    return Sleef_iunordq1(*in1, *in2) ? (
-        Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
-    ) : Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
+    return Sleef_iunordq1(*in1, *in2)   ? (Sleef_iunordq1(*in1, *in1) ? *in1 : *in2)
+           : Sleef_icmpleq1(*in1, *in2) ? *in1
+                                        : *in2;
 }
 
 static inline Sleef_quad
 quad_maximum(const Sleef_quad *in1, const Sleef_quad *in2)
 {
-    return Sleef_iunordq1(*in1, *in2) ? (
-        Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
-    ) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
+    return Sleef_iunordq1(*in1, *in2)   ? (Sleef_iunordq1(*in1, *in1) ? *in1 : *in2)
+           : Sleef_icmpgeq1(*in1, *in2) ? *in1
+                                        : *in2;
 }
 
 static inline Sleef_quad
 quad_fmin(const Sleef_quad *in1, const Sleef_quad *in2)
 {
-    return Sleef_iunordq1(*in1, *in2) ? (
-        Sleef_iunordq1(*in2, *in2) ? *in1 : *in2
-    ) : Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
+    return Sleef_iunordq1(*in1, *in2)   ? (Sleef_iunordq1(*in2, *in2) ? *in1 : *in2)
+           : Sleef_icmpleq1(*in1, *in2) ? *in1
+                                        : *in2;
 }
 
 static inline Sleef_quad
 quad_fmax(const Sleef_quad *in1, const Sleef_quad *in2)
 {
-    return Sleef_iunordq1(*in1, *in2) ? (
-        Sleef_iunordq1(*in2, *in2) ? *in1 : *in2
-    ) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
+    return Sleef_iunordq1(*in1, *in2)   ? (Sleef_iunordq1(*in2, *in2) ? *in1 : *in2)
+           : Sleef_icmpgeq1(*in1, *in2) ? *in1
+                                        : *in2;
 }
 
 static inline Sleef_quad
@@ -474,7 +519,28 @@ ld_pow(const long double *a, const long double *b)
 static inline long double
 ld_mod(const long double *a, const long double *b)
 {
-    return fmodl(*a, *b);
+    if (*b == 0.0L)
+        return NAN;
+    if (isnan(*a) || isnan(*b))
+        return isnan(*a) ? *a : *b;
+    if (isinf(*a))
+        return NAN;
+
+    if (isfinite(*a) && isinf(*b)) {
+        int sign_a = signbit(*a);
+        int sign_b = signbit(*b);
+        return (sign_a == sign_b) ? *a : *b;
+    }
+
+    long double quotient = (*a) / (*b);
+    long double floored = floorl(quotient);
+    long double result = (*a) - floored * (*b);
+
+    if (result == 0.0L) {
+        return (*b < 0.0L) ? -0.0L : 0.0L;
+    }
+
+    return result;
 }
 
 static inline long double
diff --git a/quaddtype/release_tracker.md b/quaddtype/release_tracker.md
index b01a1c7..6f7f776 100644
--- a/quaddtype/release_tracker.md
+++ b/quaddtype/release_tracker.md
@@ -1,4 +1,5 @@
 # Plan for `numpy-quaddtype` v1.0.0
+
 - [ ] High-Endian System support
 - [ ] Complete Documentation
 
@@ -17,8 +18,8 @@
 | positive      | ✅    | ✅                                                                      |
 | power         | ✅    | ✅                                                                      |
 | float_power   |       |                                                                         |
-| remainder     |       |                                                                         |
-| mod           | ✅    | ❌ _Need: basic tests + edge cases (NaN/inf/-0.0/large values)_         |
+| remainder     | ✅    | ✅                                                                      |
+| mod           | ✅    | ✅                                                                      |
 | fmod          |       |                                                                         |
 | divmod        |       |                                                                         |
 | absolute      | ✅    | ✅                                                                      |
diff --git a/quaddtype/tests/test_quaddtype.py b/quaddtype/tests/test_quaddtype.py
index 7a755a3..77a9149 100644
--- a/quaddtype/tests/test_quaddtype.py
+++ b/quaddtype/tests/test_quaddtype.py
@@ -152,3 +152,83 @@ def test_array_operations():
     expected = np.array(
         [QuadPrecision("2.0"), QuadPrecision("3.5"), QuadPrecision("5.0")])
     assert all(x == y for x, y in zip(result, expected))
+
+
+@pytest.mark.parametrize("backend", ["sleef", "longdouble"])
+@pytest.mark.parametrize("op", [np.mod, np.remainder])
+@pytest.mark.parametrize("a,b", [
+    # Basic cases - positive/negative combinations
+    (7.0, 3.0), (-7.0, 3.0), (7.0, -3.0), (-7.0, -3.0),
+
+    # Zero dividend cases
+    (0.0, 3.0), (-0.0, 3.0), (0.0, -3.0), (-0.0, -3.0),
+
+    # Cases that result in zero (sign testing)
+    (6.0, 3.0), (-6.0, 3.0), (6.0, -3.0), (-6.0, -3.0),
+    (1.0, 1.0), (-1.0, 1.0), (1.0, -1.0), (-1.0, -1.0),
+
+    # Fractional cases
+    (7.5, 2.5), (-7.5, 2.5), (7.5, -2.5), (-7.5, -2.5),
+    (0.75, 0.25), (-0.1, 0.3), (0.9, -1.0), (-1.1, -1.0),
+
+    # Large/small numbers
+    (1e10, 1e5), (-1e10, 1e5), (1e-10, 1e-5), (-1e-10, 1e-5),
+
+    # Finite % infinity cases
+    (5.0, float('inf')), (-5.0, float('inf')),
+    (5.0, float('-inf')), (-5.0, float('-inf')),
+    (0.0, float('inf')), (-0.0, float('-inf')),
+
+    # NaN cases (should return NaN)
+    (float('nan'), 3.0), (3.0, float('nan')), (float('nan'), float('nan')),
+
+    # Division by zero cases (should return NaN)
+    (5.0, 0.0), (-5.0, 0.0), (0.0, 0.0), (-0.0, 0.0),
+
+    # Infinity dividend cases (should return NaN)
+    (float('inf'), 3.0), (float('-inf'), 3.0),
+    (float('inf'), float('inf')), (float('-inf'), float('-inf')),
+])
+def test_mod(a, b, backend, op):
+    """Comprehensive test for mod operation against NumPy behavior"""
+    if backend == "sleef":
+        quad_a = QuadPrecision(str(a))
+        quad_b = QuadPrecision(str(b))
+    elif backend == "longdouble":
+        quad_a = QuadPrecision(a, backend='longdouble')
+        quad_b = QuadPrecision(b, backend='longdouble')
+    float_a = np.float64(a)
+    float_b = np.float64(b)
+
+    quad_result = op(quad_a, quad_b)
+    numpy_result = op(float_a, float_b)
+
+    # Handle NaN cases
+    if np.isnan(numpy_result):
+        assert np.isnan(
+            float(quad_result)), f"Expected NaN for {a} % {b}, got {float(quad_result)}"
+        return
+
+    if np.isinf(numpy_result):
+        assert np.isinf(
+            float(quad_result)), f"Expected inf for {a} % {b}, got {float(quad_result)}"
+        assert np.sign(numpy_result) == np.sign(
+            float(quad_result)), f"Infinity sign mismatch for {a} % {b}"
+        return
+
+    np.testing.assert_allclose(float(quad_result), numpy_result, rtol=1e-10, atol=1e-15,
+                               err_msg=f"Value mismatch for {a} % {b}")
+
+    if numpy_result == 0.0:
+        numpy_sign = np.signbit(numpy_result)
+        quad_sign = np.signbit(quad_result)
+        assert numpy_sign == quad_sign, f"Zero sign mismatch for {a} % {b}: numpy={numpy_sign}, quad={quad_sign}"
+
+    # Check that non-zero results have correct sign relative to divisor
+    if numpy_result != 0.0 and not np.isnan(b) and not np.isinf(b) and b != 0.0:
+        # In Python mod, non-zero result should have same sign as divisor (or be zero)
+        result_negative = float(quad_result) < 0
+        divisor_negative = b < 0
+        numpy_negative = numpy_result < 0
+
+        assert result_negative == numpy_negative, f"Sign mismatch for {a} % {b}: quad={result_negative}, numpy={numpy_negative}"