summaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/ldbl-128ibm/math_ldbl.h')
-rw-r--r--sysdeps/ieee754/ldbl-128ibm/math_ldbl.h64
1 files changed, 61 insertions, 3 deletions
diff --git a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
index 051352f9f7..ccb646620e 100644
--- a/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
+++ b/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h
@@ -1,6 +1,23 @@
-#ifndef _MATH_PRIVATE_H_
-#error "Never use <math_ldbl.h> directly; include <math_private.h> instead."
-#endif
+/* Manipulation of the bit representation of 'long double' quantities.
+ Copyright (C) 2006-2018 Free Software Foundation, Inc.
+ This file is part of the GNU C Library.
+
+ The GNU C Library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser General Public
+ License as published by the Free Software Foundation; either
+ version 2.1 of the License, or (at your option) any later version.
+
+ The GNU C Library is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public
+ License along with the GNU C Library; if not, see
+ <http://www.gnu.org/licenses/>. */
+
+#ifndef _MATH_LDBL_H_
+#define _MATH_LDBL_H_ 1
#include <ieee754.h>
#include <stdint.h>
@@ -230,3 +247,44 @@ ldbl_nearbyint (double a)
}
return a;
}
+
+/* Canonicalize a result from an integer rounding function, in any
+ rounding mode. *A and *AA are finite and integers, with *A being
+ nonzero; if the result is not already canonical, *AA is plus or
+ minus a power of 2 that does not exceed the least set bit in
+ *A. */
+static inline void
+ldbl_canonicalize_int (double *a, double *aa)
+{
+ /* Previously we used EXTRACT_WORDS64 from math_private.h, but in order
+ to avoid including internal headers we duplicate that code here. */
+ uint64_t ax, aax;
+ union { double value; uint64_t word; } extractor;
+ extractor.value = *a;
+ ax = extractor.word;
+ extractor.value = *aa;
+ aax = extractor.word;
+
+ int expdiff = ((ax >> 52) & 0x7ff) - ((aax >> 52) & 0x7ff);
+ if (expdiff <= 53)
+ {
+ if (expdiff == 53)
+ {
+ /* Half way between two double values; noncanonical iff the
+ low bit of A's mantissa is 1. */
+ if ((ax & 1) != 0)
+ {
+ *a += 2 * *aa;
+ *aa = -*aa;
+ }
+ }
+ else
+ {
+ /* The sum can be represented in a single double. */
+ *a += *aa;
+ *aa = 0;
+ }
+ }
+}
+
+#endif /* math_ldbl.h */