lib/totalorderl.c

   1 /* Total order 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.  */
  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 totalorderl (long double const *x, long double const *y)
  28 {
  29   return totalorder ((double const *) x, (double const *) y);
  30 }
  31
  32 #else
  33
  34 # include <float.h>
  35
  36 # ifndef LDBL_SIGNBIT_WORD
  37 #  define LDBL_SIGNBIT_WORD (-1)
  38 # endif
  39
  40 int
  41 totalorderl (long double const *x, long double const *y)
  42 {
  43   /* If the sign bits of *X and *Y differ, the one with non-zero sign bit
  44      is "smaller" than the one with sign bit == 0.  */
  45   int xs = signbit (*x);
  46   int ys = signbit (*y);
  47   if (!xs != !ys)
  48     return xs;
  49
  50   /* If one of *X, *Y is a NaN and the other isn't, the answer is easy
  51      as well: the negative NaN is "smaller", the positive NaN is "greater"
  52      than the other argument.  */
  53   int xn = isnanl (*x);
  54   int yn = isnanl (*y);
  55   if (!xn != !yn)
  56     return !xn == !xs;
  57   /* If none of *X, *Y is a NaN, the '<=' operator does the job, including
  58      for -Infinity and +Infinity.  */
  59   if (!xn)
  60     return *x <= *y;
  61
  62   /* At this point, *X and *Y are NaNs with the same sign bit.  */
  63
  64   unsigned long long extended_sign = -!!xs;
  65   unsigned long long extended_sign_hi = extended_sign;
  66 # if defined __hppa || (defined __mips__ && !MIPS_NAN2008_LONG_DOUBLE) || defined __sh__
  67   /* Invert the most significant bit of the mantissa field.  Cf. snan.h.  */
  68   extended_sign_hi ^=
  69     (1ULL << (LDBL_MANT_DIG == 106
  70               ? 51                          /* double-double representation */
  71               : (LDBL_MANT_DIG - 2) - 64)); /* quad precision representation */
  72 # endif
  73   union u { unsigned long long i[2]; long double f; } volatile xu, yu;
  74   /* Although it is tempting to initialize with {0}, Solaris cc (Sun C 5.8)
  75      on x86_64 miscompiles {0}: it initializes only the lower 80 bits,
  76      not the entire 128 bits.  */
  77   xu.i[0] = 0; xu.i[1] = 0;
  78   yu.i[0] = 0; yu.i[1] = 0;
  79   xu.f = *x;
  80   yu.f = *y;
  81
  82   /* Set BIGENDIAN to true if and only if the most significant bits of
  83      xu.f's fraction are in xu.i[0].  Use the sign bit's location to
  84      infer BIGENDIAN.  */
  85   bool bigendian;
  86   if (LDBL_SIGNBIT_WORD < 0)
  87     {
  88       union u volatile zu;
  89       zu.i[0] = 0; zu.i[1] = 0;
  90       zu.f = -zu.f;
  91       bigendian = !!zu.i[0];
  92     }
  93   else
  94     bigendian = LDBL_SIGNBIT_WORD < sizeof xu.i[0] / sizeof (unsigned);
  95
  96   unsigned long long
  97     xhi = xu.i[!bigendian] ^ extended_sign_hi,
  98     yhi = yu.i[!bigendian] ^ extended_sign_hi,
  99     xlo = xu.i[ bigendian] ^ extended_sign,
 100     ylo = yu.i[ bigendian] ^ extended_sign;
 101   return (xhi < yhi) | ((xhi == yhi) & (xlo <= ylo));
 102 }
 103
 104 #endif