math/fromfp.h

   1 /* Round to integer type.  Common helper functions.
   2    Copyright (C) 2016-2024 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <https://www.gnu.org/licenses/>.  */
  18
  19 #include <errno.h>
  20 #include <fenv.h>
  21 #include <float.h>
  22 #include <math.h>
  23 #include <math-barriers.h>
  24 #include <stdbool.h>
  25 #include <stdint.h>
  26
  27 /* The including file should have defined UNSIGNED to 0 (signed return
  28    type) or 1 (unsigned return type), INEXACT to 0 (no inexact
  29    exceptions) or 1 (raise inexact exceptions) and RET_TYPE to the
  30    return type (intmax_t or uintmax_t).  */
  31
  32 /* Return the maximum unbiased exponent for an argument (negative if
  33    NEGATIVE is set) that might be in range for a call to a fromfp
  34    function with width WIDTH (greater than 0, and not exceeding that
  35    of intmax_t).  The truncated argument may still be out of range in
  36    the case of negative arguments, and if not out of range it may
  37    become out of range as a result of rounding.  */
  38
  39 static int
  40 fromfp_max_exponent (bool negative, int width)
  41 {
  42   if (UNSIGNED)
  43     return negative ? -1 : width - 1;
  44   else
  45     return negative ? width - 1 : width - 2;
  46 }
  47
  48 /* Return the result of rounding an integer value X (passed as the
  49    absolute value; NEGATIVE is true if the value is negative), where
  50    HALF_BIT is true if the bit with value 0.5 is set and MORE_BITS is
  51    true if any lower bits are set, in the rounding direction
  52    ROUND.  */
  53
  54 static uintmax_t
  55 fromfp_round (bool negative, uintmax_t x, bool half_bit, bool more_bits,
  56               int round)
  57 {
  58   switch (round)
  59     {
  60     case FP_INT_UPWARD:
  61       return x + (!negative && (half_bit || more_bits));
  62
  63     case FP_INT_DOWNWARD:
  64       return x + (negative && (half_bit || more_bits));
  65
  66     case FP_INT_TOWARDZERO:
  67     default:
  68       /* Unknown rounding directions are defined to mean unspecified
  69          rounding; treat this as truncation.  */
  70       return x;
  71
  72     case FP_INT_TONEARESTFROMZERO:
  73       return x + half_bit;
  74
  75     case FP_INT_TONEAREST:
  76       return x + (half_bit && ((x & 1) || more_bits));
  77     }
  78 }
  79
  80 /* Integer rounding, of a value whose exponent EXPONENT did not exceed
  81    the maximum exponent MAX_EXPONENT and so did not necessarily
  82    overflow, has produced X (possibly wrapping to 0); the sign is
  83    negative if NEGATIVE is true.  Return whether this overflowed the
  84    allowed width.  */
  85
  86 static bool
  87 fromfp_overflowed (bool negative, uintmax_t x, int exponent, int max_exponent)
  88 {
  89   if (UNSIGNED)
  90     {
  91       if (negative)
  92         return x != 0;
  93       else if (max_exponent == INTMAX_WIDTH - 1)
  94         return exponent == INTMAX_WIDTH - 1 && x == 0;
  95       else
  96         return x == (1ULL << (max_exponent + 1));
  97     }
  98   else
  99     {
 100       if (negative)
 101         return exponent == max_exponent && x != (1ULL << max_exponent);
 102       else
 103         return x == (1ULL << (max_exponent + 1));
 104     }
 105 }
 106
 107 /* Handle a domain error for a call to a fromfp function with an
 108    argument which is negative if NEGATIVE is set, and specified width
 109    (not exceeding that of intmax_t) WIDTH.  The return value is
 110    unspecified (with it being unclear if the result needs to fit
 111    within WIDTH bits in this case); we choose to saturate to the given
 112    number of bits (treating NaNs like any other value).  */
 113
 114 static RET_TYPE
 115 fromfp_domain_error (bool negative, unsigned int width)
 116 {
 117   feraiseexcept (FE_INVALID);
 118   __set_errno (EDOM);
 119   /* The return value is unspecified; we choose to saturate to the
 120      given number of bits (treating NaNs like any other value).  */
 121   if (UNSIGNED)
 122     {
 123       if (negative)
 124         return 0;
 125       else if (width == INTMAX_WIDTH)
 126         return -1;
 127       else
 128         return (1ULL << width) - 1;
 129     }
 130   else
 131     {
 132       if (width == 0)
 133         return 0;
 134       else if (negative)
 135         return -(1ULL << (width - 1));
 136       else
 137         return (1ULL << (width - 1)) - 1;
 138     }
 139 }
 140
 141 /* Given X, the absolute value of a floating-point number (negative if
 142    NEGATIVE is set) truncated towards zero, where HALF_BIT is true if
 143    the bit with value 0.5 is set and MORE_BITS is true if any lower
 144    bits are set, round it in the rounding direction ROUND, handle
 145    errors and exceptions and return the appropriate return value for a
 146    fromfp function.  X originally had floating-point exponent
 147    EXPONENT, which does not exceed MAX_EXPONENT, the return value from
 148    fromfp_max_exponent with width WIDTH.  */
 149
 150 static RET_TYPE
 151 fromfp_round_and_return (bool negative, uintmax_t x, bool half_bit,
 152                          bool more_bits, int round, int exponent,
 153                          int max_exponent, unsigned int width)
 154 {
 155   uintmax_t uret = fromfp_round (negative, x, half_bit, more_bits, round);
 156   if (fromfp_overflowed (negative, uret, exponent, max_exponent))
 157     return fromfp_domain_error (negative, width);
 158
 159   if (INEXACT && (half_bit || more_bits))
 160     {
 161       /* There is no need for this to use the specific floating-point
 162          type for which this header is included, and there is no need
 163          for this header to know that type at all, so just use float
 164          here.  */
 165       float force_inexact = 1.0f + FLT_MIN;
 166       math_force_eval (force_inexact);
 167     }
 168   if (UNSIGNED)
 169     /* A negative argument not rounding to zero will already have
 170        produced a domain error.  */
 171     return uret;
 172   else
 173     return negative ? -uret : uret;
 174 }