Skip to content

Commit 4256e12

Browse files
SwayamInSyncSwayamInSync
committed
fixing to_quad
1 parent 4199b57 commit 4256e12

2 files changed

Lines changed: 84 additions & 57 deletions

File tree

‎quaddtype/numpy_quaddtype/src/casts.cpp‎

Lines changed: 36 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1079,8 +1079,18 @@ inline quad_value
10791079
to_quad<float>(float x, QuadBackendType backend)
10801080
{
10811081
quad_value result;
1082-
if (backend == BACKEND_SLEEF) {
1083-
result.sleef_value = Sleef_cast_from_doubleq1(x);
1082+
if (backend == BACKEND_SLEEF)
1083+
{
1084+
if (std::isnan(x)) {
1085+
result.sleef_value = std::signbit(x) ? QUAD_PRECISION_NEG_NAN : QUAD_PRECISION_NAN;
1086+
}
1087+
else if (std::isinf(x)) {
1088+
result.sleef_value = (x > 0) ? QUAD_PRECISION_INF : QUAD_PRECISION_NINF;
1089+
}
1090+
else {
1091+
Sleef_quad temp = Sleef_cast_from_doubleq1(static_cast<double>(x));
1092+
std::memcpy(&result.sleef_value, &temp, sizeof(Sleef_quad));
1093+
}
10841094
}
10851095
else {
10861096
result.longdouble_value = (long double)x;
@@ -1093,8 +1103,18 @@ inline quad_value
10931103
to_quad<double>(double x, QuadBackendType backend)
10941104
{
10951105
quad_value result;
1096-
if (backend == BACKEND_SLEEF) {
1097-
result.sleef_value = Sleef_cast_from_doubleq1(x);
1106+
if (backend == BACKEND_SLEEF)
1107+
{
1108+
if (std::isnan(x)) {
1109+
result.sleef_value = std::signbit(x) ? QUAD_PRECISION_NEG_NAN : QUAD_PRECISION_NAN;
1110+
}
1111+
else if (std::isinf(x)) {
1112+
result.sleef_value = (x > 0) ? QUAD_PRECISION_INF : QUAD_PRECISION_NINF;
1113+
}
1114+
else {
1115+
Sleef_quad temp = Sleef_cast_from_doubleq1(x);
1116+
std::memcpy(&result.sleef_value, &temp, sizeof(Sleef_quad));
1117+
}
10981118
}
10991119
else {
11001120
result.longdouble_value = (long double)x;
@@ -1107,8 +1127,18 @@ inline quad_value
11071127
to_quad<long double>(long double x, QuadBackendType backend)
11081128
{
11091129
quad_value result;
1110-
if (backend == BACKEND_SLEEF) {
1111-
result.sleef_value = Sleef_cast_from_doubleq1(x);
1130+
if (backend == BACKEND_SLEEF)
1131+
{
1132+
if (std::isnan(x)) {
1133+
result.sleef_value = std::signbit(x) ? QUAD_PRECISION_NEG_NAN : QUAD_PRECISION_NAN;
1134+
}
1135+
else if (std::isinf(x)) {
1136+
result.sleef_value = (x > 0) ? QUAD_PRECISION_INF : QUAD_PRECISION_NINF;
1137+
}
1138+
else {
1139+
Sleef_quad temp = Sleef_cast_from_doubleq1(static_cast<double>(x));
1140+
std::memcpy(&result.sleef_value, &temp, sizeof(Sleef_quad));
1141+
}
11121142
}
11131143
else {
11141144
result.longdouble_value = x;

‎quaddtype/numpy_quaddtype/src/ops.hpp‎

Lines changed: 48 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,7 @@
44
#include "constants.hpp"
55

66
// Quad Constants, generated with qutil
7-
#define QUAD_ZERO sleef_q(+0x0000000000000LL, 0x0000000000000000ULL, -16383)
87
#define QUAD_ONE sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 0)
9-
#define QUAD_POS_INF sleef_q(+0x1000000000000LL, 0x0000000000000000ULL, 16384)
10-
#define QUAD_NAN sleef_q(+0x1ffffffffffffLL, 0xffffffffffffffffULL, 16384)
118

129
// Unary Quad Operations
1310
typedef Sleef_quad (*unary_op_quad_def)(const Sleef_quad *);
@@ -29,7 +26,7 @@ quad_positive(const Sleef_quad *op)
2926
static inline Sleef_quad
3027
quad_sign(const Sleef_quad *op)
3128
{
32-
int sign = Sleef_icmpq1(*op, QUAD_ZERO);
29+
int sign = Sleef_icmpq1(*op, QUAD_PRECISION_ZERO);
3330
// sign(x=NaN) = x; otherwise sign(x) in { -1.0; 0.0; +1.0 }
3431
return Sleef_iunordq1(*op, *op) ? *op : Sleef_cast_from_int64q1(sign);
3532
}
@@ -97,11 +94,11 @@ quad_cbrt(const Sleef_quad *op)
9794
if (Sleef_iunordq1(*op, *op)) {
9895
return *op; // NaN
9996
}
100-
if (Sleef_icmpeqq1(*op, QUAD_ZERO)) {
97+
if (Sleef_icmpeqq1(*op, QUAD_PRECISION_ZERO)) {
10198
return *op; // ±0
10299
}
103100
// Check if op is ±inf: isinf(x) = abs(x) == inf
104-
if (Sleef_icmpeqq1(Sleef_fabsq1(*op), QUAD_POS_INF)) {
101+
if (Sleef_icmpeqq1(Sleef_fabsq1(*op), QUAD_PRECISION_INF)) {
105102
return *op; // ±inf
106103
}
107104

@@ -110,7 +107,7 @@ quad_cbrt(const Sleef_quad *op)
110107
Sleef_quad one_third = Sleef_divq1_u05(QUAD_ONE, three);
111108

112109
// Handle negative values: cbrt(-x) = -cbrt(x)
113-
if (Sleef_icmpltq1(*op, QUAD_ZERO)) {
110+
if (Sleef_icmpltq1(*op, QUAD_PRECISION_ZERO)) {
114111
Sleef_quad abs_val = Sleef_fabsq1(*op);
115112
Sleef_quad result = Sleef_powq1_u10(abs_val, one_third);
116113
return Sleef_negq1(result);
@@ -497,21 +494,21 @@ quad_signbit(const Sleef_quad *op)
497494
// once we test big and little endian in CI
498495
Sleef_quad one_signed = Sleef_copysignq1(QUAD_ONE, *op);
499496
// signbit(x) = 1 iff copysign(1, x) == -1
500-
return Sleef_icmpltq1(one_signed, QUAD_ZERO);
497+
return Sleef_icmpltq1(one_signed, QUAD_PRECISION_ZERO);
501498
}
502499

503500
static inline npy_bool
504501
quad_isfinite(const Sleef_quad *op)
505502
{
506503
// isfinite(x) = abs(x) < inf
507-
return Sleef_icmpltq1(Sleef_fabsq1(*op), QUAD_POS_INF);
504+
return Sleef_icmpltq1(Sleef_fabsq1(*op), QUAD_PRECISION_INF);
508505
}
509506

510507
static inline npy_bool
511508
quad_isinf(const Sleef_quad *op)
512509
{
513510
// isinf(x) = abs(x) == inf
514-
return Sleef_icmpeqq1(Sleef_fabsq1(*op), QUAD_POS_INF);
511+
return Sleef_icmpeqq1(Sleef_fabsq1(*op), QUAD_PRECISION_INF);
515512
}
516513

517514
static inline npy_bool
@@ -586,13 +583,13 @@ quad_floor_divide(const Sleef_quad *a, const Sleef_quad *b)
586583

587584
// inf / finite_nonzero or -inf / finite_nonzero -> NaN
588585
// But inf / 0 -> inf
589-
if (quad_isinf(a) && quad_isfinite(b) && !Sleef_icmpeqq1(*b, QUAD_ZERO)) {
590-
return QUAD_NAN;
586+
if (quad_isinf(a) && quad_isfinite(b) && !Sleef_icmpeqq1(*b, QUAD_PRECISION_ZERO)) {
587+
return QUAD_PRECISION_NAN;
591588
}
592589

593590
// 0 / 0 (including -0.0 / 0.0, 0.0 / -0.0, -0.0 / -0.0) -> NaN
594-
if (Sleef_icmpeqq1(*a, QUAD_ZERO) && Sleef_icmpeqq1(*b, QUAD_ZERO)) {
595-
return QUAD_NAN;
591+
if (Sleef_icmpeqq1(*a, QUAD_PRECISION_ZERO) && Sleef_icmpeqq1(*b, QUAD_PRECISION_ZERO)) {
592+
return QUAD_PRECISION_NAN;
596593
}
597594

598595
Sleef_quad quotient = Sleef_divq1_u05(*a, *b);
@@ -601,8 +598,8 @@ quad_floor_divide(const Sleef_quad *a, const Sleef_quad *b)
601598
// floor_divide semantics: when result is -0.0 from non-zero numerator, convert to -1.0
602599
// This happens when: (negative & non-zero)/+inf, (positive & non-zero)/-inf
603600
// But NOT when numerator is ±0.0 (then result stays as ±0.0)
604-
if (Sleef_icmpeqq1(result, QUAD_ZERO) && quad_signbit(&result) &&
605-
!Sleef_icmpeqq1(*a, QUAD_ZERO)) {
601+
if (Sleef_icmpeqq1(result, QUAD_PRECISION_ZERO) && quad_signbit(&result) &&
602+
!Sleef_icmpeqq1(*a, QUAD_PRECISION_ZERO)) {
606603
return Sleef_negq1(QUAD_ONE); // -1.0
607604
}
608605

@@ -619,8 +616,8 @@ static inline Sleef_quad
619616
quad_mod(const Sleef_quad *a, const Sleef_quad *b)
620617
{
621618
// division by zero
622-
if (Sleef_icmpeqq1(*b, QUAD_ZERO)) {
623-
return QUAD_NAN;
619+
if (Sleef_icmpeqq1(*b, QUAD_PRECISION_ZERO)) {
620+
return QUAD_PRECISION_NAN;
624621
}
625622

626623
// NaN inputs
@@ -630,7 +627,7 @@ quad_mod(const Sleef_quad *a, const Sleef_quad *b)
630627

631628
// infinity dividend -> NaN
632629
if (quad_isinf(a)) {
633-
return QUAD_NAN;
630+
return QUAD_PRECISION_NAN;
634631
}
635632

636633
// finite % inf
@@ -650,12 +647,12 @@ quad_mod(const Sleef_quad *a, const Sleef_quad *b)
650647

651648
// Handle zero result sign: when result is exactly zero,
652649
// it should have the same sign as the divisor (NumPy convention)
653-
if (Sleef_icmpeqq1(result, QUAD_ZERO)) {
654-
if (Sleef_icmpltq1(*b, QUAD_ZERO)) {
655-
return Sleef_negq1(QUAD_ZERO); // -0.0
650+
if (Sleef_icmpeqq1(result, QUAD_PRECISION_ZERO)) {
651+
if (Sleef_icmpltq1(*b, QUAD_PRECISION_ZERO)) {
652+
return Sleef_negq1(QUAD_PRECISION_ZERO); // -0.0
656653
}
657654
else {
658-
return QUAD_ZERO; // +0.0
655+
return QUAD_PRECISION_ZERO; // +0.0
659656
}
660657
}
661658

@@ -671,13 +668,13 @@ quad_fmod(const Sleef_quad *a, const Sleef_quad *b)
671668
}
672669

673670
// Division by zero -> NaN
674-
if (Sleef_icmpeqq1(*b, QUAD_ZERO)) {
675-
return QUAD_NAN;
671+
if (Sleef_icmpeqq1(*b, QUAD_PRECISION_ZERO)) {
672+
return QUAD_PRECISION_NAN;
676673
}
677674

678675
// Infinity dividend -> NaN
679676
if (quad_isinf(a)) {
680-
return QUAD_NAN;
677+
return QUAD_PRECISION_NAN;
681678
}
682679

683680
// Finite % infinity -> return dividend (same as a)
@@ -688,14 +685,14 @@ quad_fmod(const Sleef_quad *a, const Sleef_quad *b)
688685
// x - trunc(x/y) * y
689686
Sleef_quad result = Sleef_fmodq1(*a, *b);
690687

691-
if (Sleef_icmpeqq1(result, QUAD_ZERO)) {
688+
if (Sleef_icmpeqq1(result, QUAD_PRECISION_ZERO)) {
692689
// Preserve sign of dividend (first argument)
693690
Sleef_quad sign_test = Sleef_copysignq1(QUAD_ONE, *a);
694-
if (Sleef_icmpltq1(sign_test, QUAD_ZERO)) {
695-
return Sleef_negq1(QUAD_ZERO); // -0.0
691+
if (Sleef_icmpltq1(sign_test, QUAD_PRECISION_ZERO)) {
692+
return Sleef_negq1(QUAD_PRECISION_ZERO); // -0.0
696693
}
697694
else {
698-
return QUAD_ZERO; // +0.0
695+
return QUAD_PRECISION_ZERO; // +0.0
699696
}
700697
}
701698

@@ -717,7 +714,7 @@ quad_minimum(const Sleef_quad *in1, const Sleef_quad *in2)
717714
return Sleef_iunordq1(*in1, *in1) ? *in1 : *in2;
718715
}
719716
// minimum(-0.0, +0.0) = -0.0
720-
if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) {
717+
if (Sleef_icmpeqq1(*in1, QUAD_PRECISION_ZERO) && Sleef_icmpeqq1(*in2, QUAD_PRECISION_ZERO)) {
721718
return Sleef_icmpleq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2;
722719
}
723720
return Sleef_fminq1(*in1, *in2);
@@ -730,7 +727,7 @@ quad_maximum(const Sleef_quad *in1, const Sleef_quad *in2)
730727
return Sleef_iunordq1(*in1, *in1) ? *in1 : *in2;
731728
}
732729
// maximum(-0.0, +0.0) = +0.0
733-
if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) {
730+
if (Sleef_icmpeqq1(*in1, QUAD_PRECISION_ZERO) && Sleef_icmpeqq1(*in2, QUAD_PRECISION_ZERO)) {
734731
return Sleef_icmpgeq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2;
735732
}
736733
return Sleef_fmaxq1(*in1, *in2);
@@ -743,7 +740,7 @@ quad_fmin(const Sleef_quad *in1, const Sleef_quad *in2)
743740
return Sleef_iunordq1(*in2, *in2) ? *in1 : *in2;
744741
}
745742
// fmin(-0.0, +0.0) = -0.0
746-
if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) {
743+
if (Sleef_icmpeqq1(*in1, QUAD_PRECISION_ZERO) && Sleef_icmpeqq1(*in2, QUAD_PRECISION_ZERO)) {
747744
return Sleef_icmpleq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2;
748745
}
749746
return Sleef_fminq1(*in1, *in2);
@@ -756,7 +753,7 @@ quad_fmax(const Sleef_quad *in1, const Sleef_quad *in2)
756753
return Sleef_iunordq1(*in2, *in2) ? *in1 : *in2;
757754
}
758755
// maximum(-0.0, +0.0) = +0.0
759-
if (Sleef_icmpeqq1(*in1, QUAD_ZERO) && Sleef_icmpeqq1(*in2, QUAD_ZERO)) {
756+
if (Sleef_icmpeqq1(*in1, QUAD_PRECISION_ZERO) && Sleef_icmpeqq1(*in2, QUAD_PRECISION_ZERO)) {
760757
return Sleef_icmpgeq1(Sleef_copysignq1(QUAD_ONE, *in1), Sleef_copysignq1(QUAD_ONE, *in2)) ? *in1 : *in2;
761758
}
762759
return Sleef_fmaxq1(*in1, *in2);
@@ -787,14 +784,14 @@ quad_logaddexp(const Sleef_quad *x, const Sleef_quad *y)
787784

788785
// Handle infinities
789786
// If both are -inf, result is -inf
790-
Sleef_quad neg_inf = Sleef_negq1(QUAD_POS_INF);
787+
Sleef_quad neg_inf = Sleef_negq1(QUAD_PRECISION_INF);
791788
if (Sleef_icmpeqq1(*x, neg_inf) && Sleef_icmpeqq1(*y, neg_inf)) {
792789
return neg_inf;
793790
}
794791

795792
// If either is +inf, result is +inf
796-
if (Sleef_icmpeqq1(*x, QUAD_POS_INF) || Sleef_icmpeqq1(*y, QUAD_POS_INF)) {
797-
return QUAD_POS_INF;
793+
if (Sleef_icmpeqq1(*x, QUAD_PRECISION_INF) || Sleef_icmpeqq1(*y, QUAD_PRECISION_INF)) {
794+
return QUAD_PRECISION_INF;
798795
}
799796

800797
// If one is -inf, result is the other value
@@ -829,14 +826,14 @@ quad_logaddexp2(const Sleef_quad *x, const Sleef_quad *y)
829826

830827
// Handle infinities
831828
// If both are -inf, result is -inf
832-
Sleef_quad neg_inf = Sleef_negq1(QUAD_POS_INF);
829+
Sleef_quad neg_inf = Sleef_negq1(QUAD_PRECISION_INF);
833830
if (Sleef_icmpeqq1(*x, neg_inf) && Sleef_icmpeqq1(*y, neg_inf)) {
834831
return neg_inf;
835832
}
836833

837834
// If either is +inf, result is +inf
838-
if (Sleef_icmpeqq1(*x, QUAD_POS_INF) || Sleef_icmpeqq1(*y, QUAD_POS_INF)) {
839-
return QUAD_POS_INF;
835+
if (Sleef_icmpeqq1(*x, QUAD_PRECISION_INF) || Sleef_icmpeqq1(*y, QUAD_PRECISION_INF)) {
836+
return QUAD_PRECISION_INF;
840837
}
841838

842839
// If one is -inf, result is the other value
@@ -868,10 +865,10 @@ quad_heaviside(const Sleef_quad *x1, const Sleef_quad *x2)
868865
return *x1; // x1 is NaN, return NaN
869866
}
870867

871-
if (Sleef_icmpltq1(*x1, QUAD_ZERO)) {
872-
return QUAD_ZERO;
868+
if (Sleef_icmpltq1(*x1, QUAD_PRECISION_ZERO)) {
869+
return QUAD_PRECISION_ZERO;
873870
}
874-
else if (Sleef_icmpeqq1(*x1, QUAD_ZERO)) {
871+
else if (Sleef_icmpeqq1(*x1, QUAD_PRECISION_ZERO)) {
875872
return *x2; // When x1 == 0, return x2 (even if x2 is NaN)
876873
}
877874
else {
@@ -1026,17 +1023,17 @@ quad_spacing(const Sleef_quad *x)
10261023

10271024
// Handle infinity -> NaN (numpy convention)
10281025
if (quad_isinf(x)) {
1029-
return QUAD_NAN;
1026+
return QUAD_PRECISION_NAN;
10301027
}
10311028

10321029
// Determine direction based on sign of x
10331030
Sleef_quad direction;
1034-
if (Sleef_icmpltq1(*x, QUAD_ZERO)) {
1031+
if (Sleef_icmpltq1(*x, QUAD_PRECISION_ZERO)) {
10351032
// Negative: move toward -inf
1036-
direction = Sleef_negq1(QUAD_POS_INF);
1033+
direction = Sleef_negq1(QUAD_PRECISION_INF);
10371034
} else {
10381035
// Positive or zero: move toward +inf
1039-
direction = QUAD_POS_INF;
1036+
direction = QUAD_PRECISION_INF;
10401037
}
10411038

10421039
// Compute nextafter(x, direction)
@@ -1068,7 +1065,7 @@ quad_ldexp(const Sleef_quad *x, const int *exp)
10681065
}
10691066

10701067
// ±0 * 2^exp = ±0 (preserves sign of zero)
1071-
if (Sleef_icmpeqq1(*x, QUAD_ZERO)) {
1068+
if (Sleef_icmpeqq1(*x, QUAD_PRECISION_ZERO)) {
10721069
return *x;
10731070
}
10741071

@@ -1125,7 +1122,7 @@ quad_frexp(const Sleef_quad *x, int *exp)
11251122
}
11261123

11271124
// ±0 -> mantissa ±0 with exponent 0 (preserves sign of zero)
1128-
if (Sleef_icmpeqq1(*x, QUAD_ZERO)) {
1125+
if (Sleef_icmpeqq1(*x, QUAD_PRECISION_ZERO)) {
11291126
*exp = 0;
11301127
return *x;
11311128
}
@@ -1588,7 +1585,7 @@ quad_is_nonzero(const Sleef_quad *a)
15881585
{
15891586
// A value is falsy if it's exactly zero (positive or negative)
15901587
// NaN and inf are truthy
1591-
npy_bool is_zero = Sleef_icmpeqq1(*a, QUAD_ZERO);
1588+
npy_bool is_zero = Sleef_icmpeqq1(*a, QUAD_PRECISION_ZERO);
15921589
return !is_zero;
15931590
}
15941591

@@ -1665,7 +1662,7 @@ ld_logical_not(const long double *a)
16651662
static inline double
16661663
cast_sleef_to_double(const Sleef_quad in)
16671664
{
1668-
if (Sleef_icmpeqq1(in, QUAD_ZERO))
1665+
if (Sleef_icmpeqq1(in, QUAD_PRECISION_ZERO))
16691666
{
16701667
return quad_signbit(&in) ? -0.0 : 0.0;
16711668
}

0 commit comments

Comments
 (0)