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