sysdeps/ieee754/ldbl-128ibm/s_fromfpl_main.c

   1 /* Round to integer type.  ldbl-128ibm version.
   2    Copyright (C) 2016-2023 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 <math.h>
  22 #include <math_private.h>
  23 #include <stdbool.h>
  24 #include <stdint.h>
  25
  26 #define BIAS 0x3ff
  27 #define MANT_DIG 53
  28
  29 #if UNSIGNED
  30 # define RET_TYPE uintmax_t
  31 #else
  32 # define RET_TYPE intmax_t
  33 #endif
  34
  35 #include <fromfp.h>
  36
  37 RET_TYPE
  38 FUNC (long double x, int round, unsigned int width)
  39 {
  40   double hi, lo;
  41   if (width > INTMAX_WIDTH)
  42     width = INTMAX_WIDTH;
  43   uint64_t hx, lx;
  44   ldbl_unpack (x, &hi, &lo);
  45   EXTRACT_WORDS64 (hx, hi);
  46   EXTRACT_WORDS64 (lx, lo);
  47   bool negative = (hx & 0x8000000000000000ULL) != 0;
  48   bool lo_negative = (lx & 0x8000000000000000ULL) != 0;
  49   if (width == 0)
  50     return fromfp_domain_error (negative, width);
  51   hx &= 0x7fffffffffffffffULL;
  52   lx &= 0x7fffffffffffffffULL;
  53   if ((hx | lx) == 0)
  54     return 0;
  55   int hi_exponent = hx >> (MANT_DIG - 1);
  56   hi_exponent -= BIAS;
  57   int exponent = hi_exponent;
  58   hx &= ((1ULL << (MANT_DIG - 1)) - 1);
  59   if (hx == 0 && lx != 0 && lo_negative != negative)
  60     exponent--;
  61   int max_exponent = fromfp_max_exponent (negative, width);
  62   if (exponent > max_exponent)
  63     return fromfp_domain_error (negative, width);
  64   int lo_exponent = lx >> (MANT_DIG - 1);
  65   lo_exponent -= BIAS;
  66
  67   /* Convert the high part to integer.  */
  68   hx |= 1ULL << (MANT_DIG - 1);
  69   uintmax_t uret;
  70   bool half_bit, more_bits;
  71   if (hi_exponent >= MANT_DIG - 1)
  72     {
  73       uret = hx;
  74       uret <<= hi_exponent - (MANT_DIG - 1);
  75       half_bit = false;
  76       more_bits = false;
  77     }
  78   else if (hi_exponent >= -1)
  79     {
  80       uint64_t h = 1ULL << (MANT_DIG - 2 - hi_exponent);
  81       half_bit = (hx & h) != 0;
  82       more_bits = (hx & (h - 1)) != 0;
  83       uret = hx >> (MANT_DIG - 1 - hi_exponent);
  84     }
  85   else
  86     {
  87       uret = 0;
  88       half_bit = false;
  89       more_bits = true;
  90     }
  91
  92   /* Likewise, the low part.  */
  93   if (lx != 0)
  94     {
  95       uintmax_t lo_uret;
  96       bool lo_half_bit, lo_more_bits;
  97       lx &= ((1ULL << (MANT_DIG - 1)) - 1);
  98       lx |= 1ULL << (MANT_DIG - 1);
  99       /* The high part exponent is at most 64, so the low part
 100          exponent is at most 11.  */
 101       if (lo_exponent >= -1)
 102         {
 103           uint64_t h = 1ULL << (MANT_DIG - 2 - lo_exponent);
 104           lo_half_bit = (lx & h) != 0;
 105           lo_more_bits = (lx & (h - 1)) != 0;
 106           lo_uret = lx >> (MANT_DIG - 1 - lo_exponent);
 107         }
 108       else
 109         {
 110           lo_uret = 0;
 111           lo_half_bit = false;
 112           lo_more_bits = true;
 113         }
 114       if (lo_negative == negative)
 115         {
 116           uret += lo_uret;
 117           half_bit |= lo_half_bit;
 118           more_bits |= lo_more_bits;
 119         }
 120       else
 121         {
 122           uret -= lo_uret;
 123           if (lo_half_bit)
 124             {
 125               uret--;
 126               half_bit = true;
 127             }
 128           if (lo_more_bits && !more_bits)
 129             {
 130               if (half_bit)
 131                 {
 132                   half_bit = false;
 133                   more_bits = true;
 134                 }
 135               else
 136                 {
 137                   uret--;
 138                   half_bit = true;
 139                   more_bits = true;
 140                 }
 141             }
 142         }
 143     }
 144
 145   return fromfp_round_and_return (negative, uret, half_bit, more_bits, round,
 146                                   exponent, max_exponent, width);
 147 }