Skip to content

Commit e8193b4

Browse files
SwayamInSyncSwayamInSync
authored
Merge pull request #125 from SwayamInSync/main
Proceeding it to merge
2 parents 0635b84 + ac7340b commit e8193b4

3 files changed

Lines changed: 166 additions & 19 deletions

File tree

‎quaddtype/numpy_quaddtype/src/ops.hpp‎

Lines changed: 83 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#define QUAD_ZERO sleef_q(+0x0000000000000LL, 0x0000000000000000ULL, -16383)
77
#define QUAD_ONE sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 0)
88
#define QUAD_POS_INF sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384)
9+
#define QUAD_NAN sleef_q(+0x1ffffffffffffLL, 0xffffffffffffffffULL, 16384)
910

1011
// Unary Quad Operations
1112
typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
@@ -174,9 +175,12 @@ ld_absolute(const long double *op)
174175
static inline long double
175176
ld_sign(const long double *op)
176177
{
177-
if (*op < 0.0) return -1.0;
178-
if (*op == 0.0) return 0.0;
179-
if (*op > 0.0) return 1.0;
178+
if (*op < 0.0)
179+
return -1.0;
180+
if (*op == 0.0)
181+
return 0.0;
182+
if (*op > 0.0)
183+
return 1.0;
180184
// sign(x=NaN) = x
181185
return *op;
182186
}
@@ -391,39 +395,80 @@ quad_pow(const Sleef_quad *a, const Sleef_quad *b)
391395
static inline Sleef_quad
392396
quad_mod(const Sleef_quad *a, const Sleef_quad *b)
393397
{
394-
return Sleef_fmodq1(*a, *b);
398+
// division by zero
399+
if (Sleef_icmpeqq1(*b, QUAD_ZERO)) {
400+
return QUAD_NAN;
401+
}
402+
403+
// NaN inputs
404+
if (Sleef_iunordq1(*a, *b)) {
405+
return Sleef_iunordq1(*a, *a) ? *a : *b; // Return the NaN
406+
}
407+
408+
// infinity dividend -> NaN
409+
if (quad_isinf(a)) {
410+
return QUAD_NAN;
411+
}
412+
413+
// finite % inf
414+
if (quad_isfinite(a) && quad_isinf(b)) {
415+
int sign_a = quad_signbit(a);
416+
int sign_b = quad_signbit(b);
417+
418+
// return a if sign_a == sign_b
419+
return (sign_a == sign_b) ? *a : *b;
420+
}
421+
422+
// NumPy mod formula: a % b = a - floor(a/b) * b
423+
Sleef_quad quotient = Sleef_divq1_u05(*a, *b);
424+
Sleef_quad floored = Sleef_floorq1(quotient);
425+
Sleef_quad product = Sleef_mulq1_u05(floored, *b);
426+
Sleef_quad result = Sleef_subq1_u05(*a, product);
427+
428+
// Handle zero result sign: when result is exactly zero,
429+
// it should have the same sign as the divisor (NumPy convention)
430+
if (Sleef_icmpeqq1(result, QUAD_ZERO)) {
431+
if (Sleef_icmpltq1(*b, QUAD_ZERO)) {
432+
return Sleef_negq1(QUAD_ZERO); // -0.0
433+
}
434+
else {
435+
return QUAD_ZERO; // +0.0
436+
}
437+
}
438+
439+
return result;
395440
}
396441

397442
static inline Sleef_quad
398443
quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2)
399444
{
400-
return Sleef_iunordq1(*in1, *in2) ? (
401-
Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
402-
) : Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
445+
return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in1, *in1) ? *in1 : *in2)
446+
: Sleef_icmpleq1(*in1, *in2) ? *in1
447+
: *in2;
403448
}
404449

405450
static inline Sleef_quad
406451
quad_maximum(const Sleef_quad *in1, const Sleef_quad *in2)
407452
{
408-
return Sleef_iunordq1(*in1, *in2) ? (
409-
Sleef_iunordq1(*in1, *in1) ? *in1 : *in2
410-
) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
453+
return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in1, *in1) ? *in1 : *in2)
454+
: Sleef_icmpgeq1(*in1, *in2) ? *in1
455+
: *in2;
411456
}
412457

413458
static inline Sleef_quad
414459
quad_fmin(const Sleef_quad *in1, const Sleef_quad *in2)
415460
{
416-
return Sleef_iunordq1(*in1, *in2) ? (
417-
Sleef_iunordq1(*in2, *in2) ? *in1 : *in2
418-
) : Sleef_icmpleq1(*in1, *in2) ? *in1 : *in2;
461+
return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in2, *in2) ? *in1 : *in2)
462+
: Sleef_icmpleq1(*in1, *in2) ? *in1
463+
: *in2;
419464
}
420465

421466
static inline Sleef_quad
422467
quad_fmax(const Sleef_quad *in1, const Sleef_quad *in2)
423468
{
424-
return Sleef_iunordq1(*in1, *in2) ? (
425-
Sleef_iunordq1(*in2, *in2) ? *in1 : *in2
426-
) : Sleef_icmpgeq1(*in1, *in2) ? *in1 : *in2;
469+
return Sleef_iunordq1(*in1, *in2) ? (Sleef_iunordq1(*in2, *in2) ? *in1 : *in2)
470+
: Sleef_icmpgeq1(*in1, *in2) ? *in1
471+
: *in2;
427472
}
428473

429474
static inline Sleef_quad
@@ -474,7 +519,28 @@ ld_pow(const long double *a, const long double *b)
474519
static inline long double
475520
ld_mod(const long double *a, const long double *b)
476521
{
477-
return fmodl(*a, *b);
522+
if (*b == 0.0L)
523+
return NAN;
524+
if (isnan(*a) || isnan(*b))
525+
return isnan(*a) ? *a : *b;
526+
if (isinf(*a))
527+
return NAN;
528+
529+
if (isfinite(*a) && isinf(*b)) {
530+
int sign_a = signbit(*a);
531+
int sign_b = signbit(*b);
532+
return (sign_a == sign_b) ? *a : *b;
533+
}
534+
535+
long double quotient = (*a) / (*b);
536+
long double floored = floorl(quotient);
537+
long double result = (*a) - floored * (*b);
538+
539+
if (result == 0.0L) {
540+
return (*b < 0.0L) ? -0.0L : 0.0L;
541+
}
542+
543+
return result;
478544
}
479545

480546
static inline long double

‎quaddtype/release_tracker.md‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
# Plan for `numpy-quaddtype` v1.0.0
2+
23
- [ ] High-Endian System support
34
- [ ] Complete Documentation
45

@@ -17,8 +18,8 @@
1718
| positive | ✅ | ✅ |
1819
| power | ✅ | ✅ |
1920
| float_power | | |
20-
| remainder | | |
21-
| mod | ✅ | ❌ _Need: basic tests + edge cases (NaN/inf/-0.0/large values)_ |
21+
| remainder | ✅ | ✅ |
22+
| mod | ✅ | ✅ |
2223
| fmod | | |
2324
| divmod | | |
2425
| absolute | ✅ | ✅ |

‎quaddtype/tests/test_quaddtype.py‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -152,3 +152,83 @@ def test_array_operations():
152152
expected = np.array(
153153
[QuadPrecision("2.0"), QuadPrecision("3.5"), QuadPrecision("5.0")])
154154
assert all(x == y for x, y in zip(result, expected))
155+
156+
157+
@pytest.mark.parametrize("backend", ["sleef", "longdouble"])
158+
@pytest.mark.parametrize("op", [np.mod, np.remainder])
159+
@pytest.mark.parametrize("a,b", [
160+
# Basic cases - positive/negative combinations
161+
(7.0, 3.0), (-7.0, 3.0), (7.0, -3.0), (-7.0, -3.0),
162+
163+
# Zero dividend cases
164+
(0.0, 3.0), (-0.0, 3.0), (0.0, -3.0), (-0.0, -3.0),
165+
166+
# Cases that result in zero (sign testing)
167+
(6.0, 3.0), (-6.0, 3.0), (6.0, -3.0), (-6.0, -3.0),
168+
(1.0, 1.0), (-1.0, 1.0), (1.0, -1.0), (-1.0, -1.0),
169+
170+
# Fractional cases
171+
(7.5, 2.5), (-7.5, 2.5), (7.5, -2.5), (-7.5, -2.5),
172+
(0.75, 0.25), (-0.1, 0.3), (0.9, -1.0), (-1.1, -1.0),
173+
174+
# Large/small numbers
175+
(1e10, 1e5), (-1e10, 1e5), (1e-10, 1e-5), (-1e-10, 1e-5),
176+
177+
# Finite % infinity cases
178+
(5.0, float('inf')), (-5.0, float('inf')),
179+
(5.0, float('-inf')), (-5.0, float('-inf')),
180+
(0.0, float('inf')), (-0.0, float('-inf')),
181+
182+
# NaN cases (should return NaN)
183+
(float('nan'), 3.0), (3.0, float('nan')), (float('nan'), float('nan')),
184+
185+
# Division by zero cases (should return NaN)
186+
(5.0, 0.0), (-5.0, 0.0), (0.0, 0.0), (-0.0, 0.0),
187+
188+
# Infinity dividend cases (should return NaN)
189+
(float('inf'), 3.0), (float('-inf'), 3.0),
190+
(float('inf'), float('inf')), (float('-inf'), float('-inf')),
191+
])
192+
def test_mod(a, b, backend, op):
193+
"""Comprehensive test for mod operation against NumPy behavior"""
194+
if backend == "sleef":
195+
quad_a = QuadPrecision(str(a))
196+
quad_b = QuadPrecision(str(b))
197+
elif backend == "longdouble":
198+
quad_a = QuadPrecision(a, backend='longdouble')
199+
quad_b = QuadPrecision(b, backend='longdouble')
200+
float_a = np.float64(a)
201+
float_b = np.float64(b)
202+
203+
quad_result = op(quad_a, quad_b)
204+
numpy_result = op(float_a, float_b)
205+
206+
# Handle NaN cases
207+
if np.isnan(numpy_result):
208+
assert np.isnan(
209+
float(quad_result)), f"Expected NaN for {a} % {b}, got {float(quad_result)}"
210+
return
211+
212+
if np.isinf(numpy_result):
213+
assert np.isinf(
214+
float(quad_result)), f"Expected inf for {a} % {b}, got {float(quad_result)}"
215+
assert np.sign(numpy_result) == np.sign(
216+
float(quad_result)), f"Infinity sign mismatch for {a} % {b}"
217+
return
218+
219+
np.testing.assert_allclose(float(quad_result), numpy_result, rtol=1e-10, atol=1e-15,
220+
err_msg=f"Value mismatch for {a} % {b}")
221+
222+
if numpy_result == 0.0:
223+
numpy_sign = np.signbit(numpy_result)
224+
quad_sign = np.signbit(quad_result)
225+
assert numpy_sign == quad_sign, f"Zero sign mismatch for {a} % {b}: numpy={numpy_sign}, quad={quad_sign}"
226+
227+
# Check that non-zero results have correct sign relative to divisor
228+
if numpy_result != 0.0 and not np.isnan(b) and not np.isinf(b) and b != 0.0:
229+
# In Python mod, non-zero result should have same sign as divisor (or be zero)
230+
result_negative = float(quad_result) < 0
231+
divisor_negative = b < 0
232+
numpy_negative = numpy_result < 0
233+
234+
assert result_negative == numpy_negative, f"Sign mismatch for {a} % {b}: quad={result_negative}, numpy={numpy_negative}"

0 commit comments

Comments
 (0)