Skip to content

Commit 4c2863d

Browse files
SwayamInSyncSwayamInSync
authored
Merge pull request #178 from SwayamInSync/fmod
2 parents 4e5d7d3 + a83573a commit 4c2863d

File tree

4 files changed

+210
-1
lines changed

4 files changed

+210
-1
lines changed

‎quaddtype/numpy_quaddtype/src/ops.hpp‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -553,6 +553,46 @@ quad_mod(const Sleef_quad *a, const Sleef_quad *b)
553553
return result;
554554
}
555555

556+
static inline Sleef_quad
557+
quad_fmod(const Sleef_quad *a, const Sleef_quad *b)
558+
{
559+
// Handle NaN inputs
560+
if (Sleef_iunordq1(*a, *b)) {
561+
return Sleef_iunordq1(*a, *a) ? *a : *b;
562+
}
563+
564+
// Division by zero -> NaN
565+
if (Sleef_icmpeqq1(*b, QUAD_ZERO)) {
566+
return QUAD_NAN;
567+
}
568+
569+
// Infinity dividend -> NaN
570+
if (quad_isinf(a)) {
571+
return QUAD_NAN;
572+
}
573+
574+
// Finite % infinity -> return dividend (same as a)
575+
if (quad_isfinite(a) && quad_isinf(b)) {
576+
return *a;
577+
}
578+
579+
// x - trunc(x/y) * y
580+
Sleef_quad result = Sleef_fmodq1(*a, *b);
581+
582+
if (Sleef_icmpeqq1(result, QUAD_ZERO)) {
583+
// Preserve sign of dividend (first argument)
584+
Sleef_quad sign_test = Sleef_copysignq1(QUAD_ONE, *a);
585+
if (Sleef_icmpltq1(sign_test, QUAD_ZERO)) {
586+
return Sleef_negq1(QUAD_ZERO); // -0.0
587+
}
588+
else {
589+
return QUAD_ZERO; // +0.0
590+
}
591+
}
592+
593+
return result;
594+
}
595+
556596
static inline Sleef_quad
557597
quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2)
558598
{
@@ -794,6 +834,43 @@ ld_mod(const long double *a, const long double *b)
794834
return result;
795835
}
796836

837+
static inline long double
838+
ld_fmod(const long double *a, const long double *b)
839+
{
840+
// Handle NaN inputs
841+
if (isnan(*a) || isnan(*b)) {
842+
return isnan(*a) ? *a : *b;
843+
}
844+
845+
// Division by zero -> NaN
846+
if (*b == 0.0L) {
847+
return NAN;
848+
}
849+
850+
// Infinity dividend -> NaN
851+
if (isinf(*a)) {
852+
return NAN;
853+
}
854+
855+
// Finite % infinity -> return dividend
856+
if (isfinite(*a) && isinf(*b)) {
857+
return *a;
858+
}
859+
860+
long double result = fmodl(*a, *b);
861+
862+
if (result == 0.0L) {
863+
// Preserve sign of dividend
864+
if (signbit(*a)) {
865+
return -0.0L;
866+
} else {
867+
return 0.0L;
868+
}
869+
}
870+
871+
return result;
872+
}
873+
797874
static inline long double
798875
ld_minimum(const long double *in1, const long double *in2)
799876
{

‎quaddtype/numpy_quaddtype/src/umath/binary_ops.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -232,6 +232,9 @@ init_quad_binary_ops(PyObject *numpy)
232232
if (create_quad_binary_ufunc<quad_mod, ld_mod>(numpy, "mod") < 0) {
233233
return -1;
234234
}
235+
if (create_quad_binary_ufunc<quad_fmod, ld_fmod>(numpy, "fmod") < 0) {
236+
return -1;
237+
}
235238
if (create_quad_binary_ufunc<quad_minimum, ld_minimum>(numpy, "minimum") < 0) {
236239
return -1;
237240
}

‎quaddtype/release_tracker.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@
2020
| float_power | ✅ | ✅ |
2121
| remainder | ✅ | ✅ |
2222
| mod | ✅ | ✅ |
23-
| fmod | | |
23+
| fmod | ✅ | ✅ |
2424
| divmod | | |
2525
| absolute | ✅ | ✅ |
2626
| fabs | | |

‎quaddtype/tests/test_quaddtype.py‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -837,6 +837,135 @@ def test_floor_divide_special_properties():
837837
np.testing.assert_allclose(float(result_floor_divide), float(result_floor_true_divide), rtol=1e-30)
838838

839839

840+
@pytest.mark.parametrize("x_val,y_val", [
841+
(x, y) for x in [-1e10, -100.0, -7.0, -1.0, -0.5, -0.0, 0.0, 0.5, 1.0, 7.0, 100.0, 1e10,
842+
float('inf'), float('-inf'), float('nan'),
843+
-6.0, 6.0, -0.1, 0.1, -3.14159, 3.14159]
844+
for y in [-1e10, -100.0, -3.0, -1.0, -0.5, -0.0, 0.0, 0.5, 1.0, 3.0, 100.0, 1e10,
845+
float('inf'), float('-inf'), float('nan'),
846+
-2.0, 2.0, -0.25, 0.25, -1.5, 1.5]
847+
])
848+
def test_fmod(x_val, y_val):
849+
"""Test fmod ufunc with comprehensive edge cases"""
850+
x_quad = QuadPrecision(str(x_val))
851+
y_quad = QuadPrecision(str(y_val))
852+
853+
# Compute using QuadPrecision
854+
result_quad = np.fmod(x_quad, y_quad)
855+
856+
# Compute using float64 for comparison
857+
result_float64 = np.fmod(np.float64(x_val), np.float64(y_val))
858+
859+
# Compare results
860+
if np.isnan(result_float64):
861+
assert np.isnan(float(result_quad)), f"Expected NaN for fmod({x_val}, {y_val})"
862+
elif np.isinf(result_float64):
863+
assert np.isinf(float(result_quad)), f"Expected inf for fmod({x_val}, {y_val})"
864+
assert np.sign(float(result_quad)) == np.sign(result_float64), f"Sign mismatch for fmod({x_val}, {y_val})"
865+
else:
866+
# For finite results, check relative tolerance
867+
atol = max(1e-10, abs(result_float64) * 1e-9) if abs(result_float64) > 1e6 else 1e-10
868+
np.testing.assert_allclose(
869+
float(result_quad), result_float64, rtol=1e-12, atol=atol,
870+
err_msg=f"Mismatch for fmod({x_val}, {y_val})"
871+
)
872+
873+
# Critical: Check sign preservation for zero results
874+
if result_float64 == 0.0:
875+
assert np.signbit(result_quad) == np.signbit(result_float64), \
876+
f"Sign mismatch for zero result: fmod({x_val}, {y_val}), " \
877+
f"expected signbit={np.signbit(result_float64)}, got signbit={np.signbit(result_quad)}"
878+
879+
880+
def test_fmod_special_properties():
881+
"""Test special mathematical properties of fmod"""
882+
# fmod(x, 1) gives fractional part of x (with sign preserved)
883+
x = QuadPrecision("42.7")
884+
result = np.fmod(x, QuadPrecision("1.0"))
885+
np.testing.assert_allclose(float(result), 0.7, rtol=1e-15, atol=1e-15)
886+
887+
# fmod(0, non-zero) = 0 with correct sign
888+
result = np.fmod(QuadPrecision("0.0"), QuadPrecision("5.0"))
889+
assert float(result) == 0.0 and not np.signbit(result)
890+
891+
result = np.fmod(QuadPrecision("-0.0"), QuadPrecision("5.0"))
892+
assert float(result) == 0.0 and np.signbit(result)
893+
894+
# fmod by 0 gives NaN
895+
result = np.fmod(QuadPrecision("1.0"), QuadPrecision("0.0"))
896+
assert np.isnan(float(result))
897+
898+
result = np.fmod(QuadPrecision("-1.0"), QuadPrecision("0.0"))
899+
assert np.isnan(float(result))
900+
901+
# 0 fmod 0 = NaN
902+
result = np.fmod(QuadPrecision("0.0"), QuadPrecision("0.0"))
903+
assert np.isnan(float(result))
904+
905+
# inf fmod x = NaN
906+
result = np.fmod(QuadPrecision("inf"), QuadPrecision("100.0"))
907+
assert np.isnan(float(result))
908+
909+
result = np.fmod(QuadPrecision("-inf"), QuadPrecision("100.0"))
910+
assert np.isnan(float(result))
911+
912+
# x fmod inf = x (for finite x)
913+
result = np.fmod(QuadPrecision("100.0"), QuadPrecision("inf"))
914+
np.testing.assert_allclose(float(result), 100.0, rtol=1e-30)
915+
916+
result = np.fmod(QuadPrecision("-100.0"), QuadPrecision("inf"))
917+
np.testing.assert_allclose(float(result), -100.0, rtol=1e-30)
918+
919+
# inf fmod inf = NaN
920+
result = np.fmod(QuadPrecision("inf"), QuadPrecision("inf"))
921+
assert np.isnan(float(result))
922+
923+
# fmod uses truncated division (rounds toward zero)
924+
# Result has same sign as dividend (first argument)
925+
result = np.fmod(QuadPrecision("7.0"), QuadPrecision("3.0"))
926+
assert float(result) == 1.0 # 7 - trunc(7/3)*3 = 7 - 2*3 = 1
927+
928+
result = np.fmod(QuadPrecision("-7.0"), QuadPrecision("3.0"))
929+
assert float(result) == -1.0 # -7 - trunc(-7/3)*3 = -7 - (-2)*3 = -1
930+
931+
result = np.fmod(QuadPrecision("7.0"), QuadPrecision("-3.0"))
932+
assert float(result) == 1.0 # 7 - trunc(7/-3)*(-3) = 7 - (-2)*(-3) = 1
933+
934+
result = np.fmod(QuadPrecision("-7.0"), QuadPrecision("-3.0"))
935+
assert float(result) == -1.0 # -7 - trunc(-7/-3)*(-3) = -7 - 2*(-3) = -1
936+
937+
# Sign preservation when result is exactly zero
938+
result = np.fmod(QuadPrecision("6.0"), QuadPrecision("3.0"))
939+
assert float(result) == 0.0 and not np.signbit(result)
940+
941+
result = np.fmod(QuadPrecision("-6.0"), QuadPrecision("3.0"))
942+
assert float(result) == 0.0 and np.signbit(result)
943+
944+
result = np.fmod(QuadPrecision("6.0"), QuadPrecision("-3.0"))
945+
assert float(result) == 0.0 and not np.signbit(result)
946+
947+
result = np.fmod(QuadPrecision("-6.0"), QuadPrecision("-3.0"))
948+
assert float(result) == 0.0 and np.signbit(result)
949+
950+
# Difference from mod/remainder (which uses floor division)
951+
# fmod result has sign of dividend, mod result has sign of divisor
952+
x = QuadPrecision("-7.0")
953+
y = QuadPrecision("3.0")
954+
fmod_result = np.fmod(x, y)
955+
mod_result = np.remainder(x, y)
956+
957+
assert float(fmod_result) == -1.0 # sign of dividend (negative)
958+
assert float(mod_result) == 2.0 # sign of divisor (positive)
959+
960+
# Relationship: x = trunc(x/y) * y + fmod(x, y)
961+
x = QuadPrecision("10.5")
962+
y = QuadPrecision("3.2")
963+
quotient = np.trunc(np.true_divide(x, y))
964+
remainder = np.fmod(x, y)
965+
reconstructed = np.add(np.multiply(quotient, y), remainder)
966+
np.testing.assert_allclose(float(reconstructed), float(x), rtol=1e-30)
967+
968+
840969
def test_inf():
841970
assert QuadPrecision("inf") > QuadPrecision("1e1000")
842971
assert np.signbit(QuadPrecision("inf")) == 0

0 commit comments

Comments
 (0)