Skip to content

Commit d90519b

Browse files
SwayamInSyncSwayamInSync
authored
Merge pull request #200 from SwayamInSync/frexp
2 parents 6cfc296 + 6727239 commit d90519b

File tree

5 files changed

+417
-3
lines changed

5 files changed

+417
-3
lines changed

‎quaddtype/numpy_quaddtype/src/ops.hpp‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1051,6 +1051,10 @@ quad_spacing(const Sleef_quad *x)
10511051
typedef Sleef_quad (*ldexp_op_quad_def)(const Sleef_quad *, const int *);
10521052
typedef long double (*ldexp_op_longdouble_def)(const long double *, const int *);
10531053

1054+
// Frexp operations: quad -> (quad mantissa, int exponent)
1055+
typedef Sleef_quad (*frexp_op_quad_def)(const Sleef_quad *, int *);
1056+
typedef long double (*frexp_op_longdouble_def)(const long double *, int *);
1057+
10541058
static inline Sleef_quad
10551059
quad_ldexp(const Sleef_quad *x, const int *exp)
10561060
{
@@ -1103,6 +1107,67 @@ ld_ldexp(const long double *x, const int *exp)
11031107
return result;
11041108
}
11051109

1110+
static inline Sleef_quad
1111+
quad_frexp(const Sleef_quad *x, int *exp)
1112+
{
1113+
// frexp(x) returns mantissa m and exponent e such that x = m * 2^e
1114+
// where 0.5 <= |m| < 1.0
1115+
// NumPy's documentation says "between -1 and 1" but actual behavior is:
1116+
// - Positive x: mantissa in [0.5, 1.0)
1117+
// - Negative x: mantissa in (-1.0, -0.5]
1118+
// This matches SLEEF's Sleef_frexpq1 behavior exactly.
1119+
1120+
// NaN input -> NaN output with exponent 0
1121+
if (Sleef_iunordq1(*x, *x)) {
1122+
*exp = 0;
1123+
return *x;
1124+
}
1125+
1126+
// ±0 -> mantissa ±0 with exponent 0 (preserves sign of zero)
1127+
if (Sleef_icmpeqq1(*x, QUAD_ZERO)) {
1128+
*exp = 0;
1129+
return *x;
1130+
}
1131+
1132+
// ±inf -> mantissa ±inf with exponent 0 (preserves sign of infinity)
1133+
if (quad_isinf(x)) {
1134+
*exp = 0;
1135+
return *x;
1136+
}
1137+
1138+
Sleef_quad mantissa = Sleef_frexpq1(*x, exp);
1139+
1140+
return mantissa;
1141+
}
1142+
1143+
static inline long double
1144+
ld_frexp(const long double *x, int *exp)
1145+
{
1146+
// frexp(x) returns mantissa m and exponent e such that x = m * 2^e
1147+
1148+
// NaN input -> NaN output with exponent 0
1149+
if (isnan(*x)) {
1150+
*exp = 0;
1151+
return *x;
1152+
}
1153+
1154+
// ±0 -> mantissa ±0 with exponent 0 (preserves sign of zero)
1155+
if (*x == 0.0L) {
1156+
*exp = 0;
1157+
return *x;
1158+
}
1159+
1160+
// ±inf -> mantissa ±inf with exponent 0 (preserves sign of infinity)
1161+
if (isinf(*x)) {
1162+
*exp = 0;
1163+
return *x;
1164+
}
1165+
1166+
long double mantissa = frexpl(*x, exp);
1167+
1168+
return mantissa;
1169+
}
1170+
11061171
// Binary long double operations
11071172
typedef long double (*binary_op_longdouble_def)(const long double *, const long double *);
11081173
// Binary long double operations with 2 outputs (for divmod, modf, frexp)

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -454,7 +454,7 @@ create_quad_ldexp_ufunc(PyObject *numpy, const char *ufunc_name)
454454
return 0;
455455
}
456456

457-
// Create binary ufunc with 2 outputs (generic for divmod, modf, frexp, etc.)
457+
// Create binary ufunc with 2 outputs
458458
template <binary_op_2out_quad_def sleef_op, binary_op_2out_longdouble_def longdouble_op>
459459
int
460460
create_quad_binary_2out_ufunc(PyObject *numpy, const char *ufunc_name)

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

Lines changed: 164 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -408,6 +408,165 @@ create_quad_unary_2out_ufunc(PyObject *numpy, const char *ufunc_name)
408408
return 0;
409409
}
410410

411+
// Frexp-specific resolver: QuadPrecDType -> (QuadPrecDType, int32)
412+
static NPY_CASTING
413+
quad_frexp_resolve_descriptors(PyObject *self, PyArray_DTypeMeta *const dtypes[],
414+
PyArray_Descr *const given_descrs[], PyArray_Descr *loop_descrs[],
415+
npy_intp *NPY_UNUSED(view_offset))
416+
{
417+
// Input descriptor
418+
Py_INCREF(given_descrs[0]);
419+
loop_descrs[0] = given_descrs[0];
420+
421+
// Output 1: mantissa (same dtype as input)
422+
if (given_descrs[1] == NULL) {
423+
Py_INCREF(given_descrs[0]);
424+
loop_descrs[1] = given_descrs[0];
425+
}
426+
else {
427+
Py_INCREF(given_descrs[1]);
428+
loop_descrs[1] = given_descrs[1];
429+
}
430+
431+
// Output 2: exponent (int32)
432+
if (given_descrs[2] == NULL) {
433+
loop_descrs[2] = PyArray_DescrFromType(NPY_INT32);
434+
}
435+
else {
436+
Py_INCREF(given_descrs[2]);
437+
loop_descrs[2] = given_descrs[2];
438+
}
439+
440+
QuadPrecDTypeObject *descr_in = (QuadPrecDTypeObject *)given_descrs[0];
441+
QuadPrecDTypeObject *descr_out1 = (QuadPrecDTypeObject *)loop_descrs[1];
442+
443+
if (descr_in->backend != descr_out1->backend) {
444+
return NPY_UNSAFE_CASTING;
445+
}
446+
447+
return NPY_NO_CASTING;
448+
}
449+
450+
// Strided loop for frexp (unaligned)
451+
template <frexp_op_quad_def sleef_op, frexp_op_longdouble_def longdouble_op>
452+
int
453+
quad_frexp_strided_loop_unaligned(PyArrayMethod_Context *context, char *const data[],
454+
npy_intp const dimensions[], npy_intp const strides[],
455+
NpyAuxData *auxdata)
456+
{
457+
npy_intp N = dimensions[0];
458+
char *in_ptr = data[0];
459+
char *out_mantissa_ptr = data[1];
460+
char *out_exp_ptr = data[2];
461+
npy_intp in_stride = strides[0];
462+
npy_intp out_mantissa_stride = strides[1];
463+
npy_intp out_exp_stride = strides[2];
464+
465+
QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0];
466+
QuadBackendType backend = descr->backend;
467+
size_t elem_size = (backend == BACKEND_SLEEF) ? sizeof(Sleef_quad) : sizeof(long double);
468+
469+
quad_value in, out_mantissa;
470+
int out_exp;
471+
472+
while (N--) {
473+
memcpy(&in, in_ptr, elem_size);
474+
if (backend == BACKEND_SLEEF) {
475+
out_mantissa.sleef_value = sleef_op(&in.sleef_value, &out_exp);
476+
}
477+
else {
478+
out_mantissa.longdouble_value = longdouble_op(&in.longdouble_value, &out_exp);
479+
}
480+
memcpy(out_mantissa_ptr, &out_mantissa, elem_size);
481+
memcpy(out_exp_ptr, &out_exp, sizeof(int));
482+
483+
in_ptr += in_stride;
484+
out_mantissa_ptr += out_mantissa_stride;
485+
out_exp_ptr += out_exp_stride;
486+
}
487+
return 0;
488+
}
489+
490+
// Strided loop for frexp (aligned)
491+
template <frexp_op_quad_def sleef_op, frexp_op_longdouble_def longdouble_op>
492+
int
493+
quad_frexp_strided_loop_aligned(PyArrayMethod_Context *context, char *const data[],
494+
npy_intp const dimensions[], npy_intp const strides[],
495+
NpyAuxData *auxdata)
496+
{
497+
npy_intp N = dimensions[0];
498+
char *in_ptr = data[0];
499+
char *out_mantissa_ptr = data[1];
500+
char *out_exp_ptr = data[2];
501+
npy_intp in_stride = strides[0];
502+
npy_intp out_mantissa_stride = strides[1];
503+
npy_intp out_exp_stride = strides[2];
504+
505+
QuadPrecDTypeObject *descr = (QuadPrecDTypeObject *)context->descriptors[0];
506+
QuadBackendType backend = descr->backend;
507+
508+
int out_exp;
509+
510+
while (N--) {
511+
if (backend == BACKEND_SLEEF) {
512+
Sleef_quad mantissa = sleef_op((Sleef_quad *)in_ptr, &out_exp);
513+
memcpy(out_mantissa_ptr, &mantissa, sizeof(Sleef_quad));
514+
}
515+
else {
516+
long double mantissa = longdouble_op((long double *)in_ptr, &out_exp);
517+
memcpy(out_mantissa_ptr, &mantissa, sizeof(long double));
518+
}
519+
memcpy(out_exp_ptr, &out_exp, sizeof(int));
520+
521+
in_ptr += in_stride;
522+
out_mantissa_ptr += out_mantissa_stride;
523+
out_exp_ptr += out_exp_stride;
524+
}
525+
return 0;
526+
}
527+
528+
529+
template <frexp_op_quad_def sleef_op, frexp_op_longdouble_def longdouble_op>
530+
int
531+
create_quad_frexp_ufunc(PyObject *numpy, const char *ufunc_name)
532+
{
533+
PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name);
534+
if (ufunc == NULL) {
535+
return -1;
536+
}
537+
538+
// 1 input (QuadPrecDType), 2 outputs (QuadPrecDType, int32)
539+
PyArray_DTypeMeta *dtypes[3] = {
540+
&QuadPrecDType,
541+
&QuadPrecDType,
542+
&PyArray_Int32DType
543+
};
544+
545+
PyType_Slot slots[] = {
546+
{NPY_METH_resolve_descriptors, (void *)&quad_frexp_resolve_descriptors},
547+
{NPY_METH_strided_loop,
548+
(void *)&quad_frexp_strided_loop_aligned<sleef_op, longdouble_op>},
549+
{NPY_METH_unaligned_strided_loop,
550+
(void *)&quad_frexp_strided_loop_unaligned<sleef_op, longdouble_op>},
551+
{0, NULL}};
552+
553+
PyArrayMethod_Spec Spec = {
554+
.name = "quad_frexp",
555+
.nin = 1,
556+
.nout = 2,
557+
.casting = NPY_NO_CASTING,
558+
.flags = NPY_METH_SUPPORTS_UNALIGNED,
559+
.dtypes = dtypes,
560+
.slots = slots,
561+
};
562+
563+
if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) {
564+
return -1;
565+
}
566+
567+
return 0;
568+
}
569+
411570
int
412571
init_quad_unary_ops(PyObject *numpy)
413572
{
@@ -540,5 +699,10 @@ init_quad_unary_ops(PyObject *numpy)
540699
return -1;
541700
}
542701

702+
// Frexp: special case with (QuadPrecDType, int32) outputs
703+
if (create_quad_frexp_ufunc<quad_frexp, ld_frexp>(numpy, "frexp") < 0) {
704+
return -1;
705+
}
706+
543707
return 0;
544708
}

‎quaddtype/release_tracker.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -81,7 +81,7 @@
8181
| spacing | ✅ | ✅ |
8282
| modf | ✅ | ✅ |
8383
| ldexp | ✅ | ✅ |
84-
| frexp | | |
84+
| frexp | ✅ | ✅ |
8585
| floor | ✅ | ✅ |
8686
| ceil | ✅ | ✅ |
8787
| trunc | ✅ | ✅ |

0 commit comments

Comments
 (0)