lib/totalorder.c

   1 /* Total order for 'double'
   2    Copyright 2023-2025 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.  */
  18
  19 #include <config.h>
  20
  21 /* Specification.  */
  22 #include <math.h>
  23
  24 #include <string.h>
  25
  26 int
  27 totalorder (double const *x, double const *y)
  28 {
  29   /* If the sign bits of *X and *Y differ, the one with non-zero sign bit
  30      is "smaller" than the one with sign bit == 0.  */
  31   int xs = signbit (*x);
  32   int ys = signbit (*y);
  33   if (!xs != !ys)
  34     return xs;
  35
  36   /* If one of *X, *Y is a NaN and the other isn't, the answer is easy
  37      as well: the negative NaN is "smaller", the positive NaN is "greater"
  38      than the other argument.  */
  39   int xn = isnand (*x);
  40   int yn = isnand (*y);
  41   if (!xn != !yn)
  42     return !xn == !xs;
  43   /* If none of *X, *Y is a NaN, the '<=' operator does the job, including
  44      for -Infinity and +Infinity.  */
  45   if (!xn)
  46     return *x <= *y;
  47
  48   /* At this point, *X and *Y are NaNs with the same sign bit.  */
  49
  50   unsigned long long extended_sign = -!!xs;
  51 #if defined __hppa || (defined __mips__ && !MIPS_NAN2008_DOUBLE) || defined __sh__
  52   /* Invert the most significant bit of the mantissa field.  Cf. snan.h.  */
  53   extended_sign ^= (1ULL << 51);
  54 #endif
  55   union { unsigned long long i; double f; } xu = {0}, yu = {0};
  56 #if 0
  57   xu.f = *x;
  58   yu.f = *y;
  59 #else
  60 # if defined __GNUC__ || defined __clang__
  61   /* Prevent gcc and clang from reusing the values of *x and *y (fetched above)
  62      in optimized inlined memcpy expansions.
  63      Seen with gcc <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=114659>
  64      and with clang 16.0.6 on OpenBSD 7.5.  */
  65   __asm__ __volatile__ ("" : : : "memory");
  66 # endif
  67   /* On 32-bit x86 processors, as well as on x86_64 processors with
  68      CC="gcc -mfpmath=387", the evaluation of *x and *y above is done through
  69      an 'fldl' instruction, which converts a signalling NaN to a quiet NaN. See
  70      <https://lists.gnu.org/archive/html/bug-gnulib/2023-10/msg00060.html>
  71      for details.  Use memcpy to avoid this.  */
  72   memcpy (&xu.f, x, sizeof (double));
  73   memcpy (&yu.f, y, sizeof (double));
  74 #endif
  75   return (xu.i ^ extended_sign) <= (yu.i ^ extended_sign);
  76 }