Branch
Hash :
7b089321
Author :
Date :
2025-01-01T09:24:36
maint: run 'make update-copyright'
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
/* Total order of absolute value for 'long double'.
Copyright 2023-2025 Free Software Foundation, Inc.
This file 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 3 of the
License, or (at your option) any later version.
This file 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 this program. If not, see <https://www.gnu.org/licenses/>. */
/* Written by Paul Eggert and Bruno Haible. */
#include <config.h>
/* Specification. */
#include <math.h>
#if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
int
totalordermagl (long double const *x, long double const *y)
{
return totalordermag ((double const *) x, (double const *) y);
}
#else
# include <float.h>
# include <stdint.h>
# include <string.h>
# include "verify.h"
int
totalordermagl (long double const *x, long double const *y)
{
/* If one of *X, *Y is a NaN and the other isn't, the answer is easy:
the NaN is "greater" than the other argument. */
int xn = isnanl (*x);
int yn = isnanl (*y);
if (!xn != !yn)
return yn;
/* If none of *X, *Y is a NaN, the '<=' operator on the absolute values
does the job, including for -Infinity and +Infinity. */
if (!xn)
return (signbit (*x) ? - *x : *x) <= (signbit (*y) ? - *y : *y);
/* At this point, *X and *Y are NaNs. */
# if defined LDBL_SIGNBIT_WORD && defined LDBL_SIGNBIT_BIT
/* The use of a union to extract the bits of the representation of a
'double' is safe in practice, despite of the "aliasing rules" of
C99, because the GCC docs say
"Even with '-fstrict-aliasing', type-punning is allowed, provided the
memory is accessed through the union type."
and similarly for other compilers. */
# define NWORDS \
((sizeof (long double) + sizeof (unsigned int) - 1) / sizeof (unsigned int))
# if LDBL_MANT_DIG == 64 /* on i386, x86_64, ia64, m68k */
/* A single uint64_t is enough to hold the payload, since we ignore the sign
bit and the exponent has the maximum possible value. */
verify (NWORDS >= 3);
# if LDBL_SIGNBIT_WORD == 2 && LDBL_SIGNBIT_BIT == 15 /* on i386, x86_64, ia64 */
/* Also, LDBL_EXPBIT0_WORD == 2 && LDBL_EXPBIT0_BIT == 0. */
union { unsigned int word[NWORDS]; uint64_t i; long double value; } xu, yu;
xu.value = *x;
yu.value = *y;
return xu.i <= yu.i;
# elif LDBL_SIGNBIT_WORD == 0 && LDBL_SIGNBIT_BIT == 31 /* on m68k */
/* Also, LDBL_EXPBIT0_WORD == 0 && LDBL_EXPBIT0_BIT == 16. */
union { unsigned int word[NWORDS]; long double value; } xu, yu;
xu.value = *x;
yu.value = *y;
/* The payload is in word[1] (high part) and word[2] (low part). */
return (xu.word[1] < yu.word[1])
| ((xu.word[1] == yu.word[1]) & (xu.word[2] <= yu.word[2]));
# else
# error "Please port gnulib totalordermagl.c to your platform!"
# endif
# elif LDBL_MANT_DIG == 106 /* on powerpc, powerpc64, powerpc64le */
/* Two uint64_t are needed to hold the payload.
In the double-double representation, each of the two uint64_t holds
a sign bit. */
verify (NWORDS == 4);
union { unsigned int word[NWORDS]; uint64_t i[2]; long double value; }
xu = {0}, yu = {0};
xu.value = *x;
yu.value = *y;
uint64_t xhi = xu.i[0] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32));
uint64_t xlo = xu.i[1] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32));
uint64_t yhi = yu.i[0] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32));
uint64_t ylo = yu.i[1] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32));
return (xhi < yhi) | ((xhi == yhi) & (xlo <= ylo));
# else
/* Here, LDBL_MANT_DIG == 113 (alpha, arm64, loongarch64, mips64, riscv64,
s390x, sparc64). */
/* Two uint64_t are needed to hold the payload. */
verify (NWORDS == 4);
union { unsigned int word[NWORDS]; uint64_t i[2]; long double value; }
xu = {0}, yu = {0};
# if 0
xu.value = *x;
yu.value = *y;
# else
/* On Linux/SPARC in 64-bit mode, gcc-6.4.0 -O2 miscompiles the simple
assignments above. Use memcpy as a workaround. */
memcpy (&xu.value, x, sizeof (long double));
memcpy (&yu.value, y, sizeof (long double));
# endif
uint64_t xhi;
uint64_t xlo;
uint64_t yhi;
uint64_t ylo;
# if LDBL_SIGNBIT_WORD < 2 /* big-endian */
xhi = xu.i[0] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32 * (1 - LDBL_SIGNBIT_WORD)));
xlo = xu.i[1];
yhi = yu.i[0] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32 * (1 - LDBL_SIGNBIT_WORD)));
ylo = yu.i[1];
# else /* little-endian */
xhi = xu.i[1] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32 * (LDBL_SIGNBIT_WORD - 2)));
xlo = xu.i[0];
yhi = yu.i[1] & ~((uint64_t) 1 << (LDBL_SIGNBIT_BIT + 32 * (LDBL_SIGNBIT_WORD - 2)));
ylo = yu.i[0];
# endif
# if defined __hppa || (defined __mips__ && !MIPS_NAN2008_LONG_DOUBLE) || defined __sh__
/* Invert the most significant bit of the mantissa field. Cf. snan.h. */
xhi ^=
(1ULL << (LDBL_MANT_DIG == 106
? 51 /* double-double representation */
: (LDBL_MANT_DIG - 2) - 64)); /* quad precision representation */
yhi ^=
(1ULL << (LDBL_MANT_DIG == 106
? 51 /* double-double representation */
: (LDBL_MANT_DIG - 2) - 64)); /* quad precision representation */
# endif
return (xhi < yhi) | ((xhi == yhi) & (xlo <= ylo));
# endif
# else
# error "Please port gnulib totalordermagl.c to your platform!"
# endif
}
#endif