lib/totalordermagl.c

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