libgcc/config/libbid/bid128_quantize.c

   1 /* Copyright (C) 2007  Free Software Foundation, Inc.
   2
   3 This file is part of GCC.
   4
   5 GCC is free software; you can redistribute it and/or modify it under
   6 the terms of the GNU General Public License as published by the Free
   7 Software Foundation; either version 2, or (at your option) any later
   8 version.
   9
  10 In addition to the permissions in the GNU General Public License, the
  11 Free Software Foundation gives you unlimited permission to link the
  12 compiled version of this file into combinations with other programs,
  13 and to distribute those combinations without any restriction coming
  14 from the use of this file.  (The General Public License restrictions
  15 do apply in other respects; for example, they cover modification of
  16 the file, and distribution when not linked into a combine
  17 executable.)
  18
  19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  21 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  22 for more details.
  23
  24 You should have received a copy of the GNU General Public License
  25 along with GCC; see the file COPYING.  If not, write to the Free
  26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  27 02110-1301, USA.  */
  28
  29 #define BID_128RES
  30 #include "bid_internal.h"
  31
  32 BID128_FUNCTION_ARG2(__bid128_quantize, x, y)
  33
  34   UINT256 CT;
  35   UINT128 CX, CY, T, CX2, CR, Stemp, res, REM_H, C2N;
  36   UINT64 sign_x, sign_y, remainder_h, carry, CY64;
  37   int_float tempx;
  38   int exponent_x = 0, exponent_y = 0, digits_x, extra_digits, amount;
  39   int expon_diff, total_digits, bin_expon_cx, rmode, status;
  40
  41   // unpack arguments, check for NaN or Infinity
  42   if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
  43     // y is Inf. or NaN
  44 #ifdef SET_STATUS_FLAGS
  45     if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
  46       __set_status_flags (pfpsf, INVALID_EXCEPTION);
  47 #endif
  48
  49     // test if y is NaN
  50     if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
  51 #ifdef SET_STATUS_FLAGS
  52       if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
  53         // set status flags
  54         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  55       }
  56 #endif
  57       res.w[1] = y.w[1] & QUIET_MASK64;
  58       res.w[0] = y.w[0];
  59       BID_RETURN (res);
  60     }
  61     // y is Infinity?
  62     if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
  63       // check if x is not Inf.
  64       if (((x.w[1] & 0x7c00000000000000ull) < 0x7800000000000000ull)) {
  65         // return NaN
  66 #ifdef SET_STATUS_FLAGS
  67         // set status flags
  68         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  69 #endif
  70         res.w[1] = 0x7c00000000000000ull;
  71         res.w[0] = 0;
  72         BID_RETURN (res);
  73       } else if (((x.w[1] & 0x7c00000000000000ull) <= 0x7800000000000000ull)) {
  74         res.w[1] = x.w[1];
  75         res.w[0] = x.w[0];
  76         BID_RETURN (res);
  77       }
  78     }
  79   }
  80
  81   if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
  82     // test if x is NaN or Inf
  83     if ((x.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull) {
  84 #ifdef SET_STATUS_FLAGS
  85       // set status flags
  86       __set_status_flags (pfpsf, INVALID_EXCEPTION);
  87 #endif
  88       res.w[1] = 0x7c00000000000000ull;
  89       res.w[0] = x.w[0];
  90       BID_RETURN (res);
  91     } else if ((x.w[1] & 0x7c00000000000000ull) ==
  92                0x7c00000000000000ull) {
  93       if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
  94 #ifdef SET_STATUS_FLAGS
  95         // set status flags
  96         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  97 #endif
  98       }
  99       res.w[1] = x.w[1] & QUIET_MASK64;
 100       res.w[0] = x.w[0];
 101       BID_RETURN (res);
 102     }
 103     if (!CX.w[1] && !CX.w[0]) {
 104       get_BID128_very_fast (&res, sign_x, exponent_y, CX);
 105       BID_RETURN (res);
 106     }
 107   }
 108   // get number of decimal digits in coefficient_x
 109   if (CX.w[1]) {
 110     tempx.d = (float) CX.w[1];
 111     bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f + 64;
 112   } else {
 113     tempx.d = (float) CX.w[0];
 114     bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
 115   }
 116   digits_x = __bid_estimate_decimal_digits[bin_expon_cx];
 117   if (CX.w[1] > __bid_power10_table_128[digits_x].w[1]
 118       || (CX.w[1] == __bid_power10_table_128[digits_x].w[1]
 119           && CX.w[0] >= __bid_power10_table_128[digits_x].w[0]))
 120     digits_x++;
 121
 122   expon_diff = exponent_x - exponent_y;
 123   total_digits = digits_x + expon_diff;
 124
 125   if ((UINT32) total_digits <= 34) {
 126     if (expon_diff >= 0) {
 127       T = __bid_power10_table_128[expon_diff];
 128       __mul_128x128_low (CX2, T, CX);
 129       get_BID128_very_fast (&res, sign_x, exponent_y, CX2);
 130       BID_RETURN (res);
 131     }
 132 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 133 #ifndef IEEE_ROUND_NEAREST
 134     rmode = rnd_mode;
 135     if (sign_x && (unsigned) (rmode - 1) < 2)
 136       rmode = 3 - rmode;
 137 #else
 138     rmode = 0;
 139 #endif
 140 #else
 141     rmode = 0;
 142 #endif
 143     // must round off -expon_diff digits
 144     extra_digits = -expon_diff;
 145     __add_128_128 (CX, CX, __bid_round_const_table_128[rmode][extra_digits]);
 146
 147     // get P*(2^M[extra_digits])/10^extra_digits
 148     __mul_128x128_to_256 (CT, CX, __bid_reciprocals10_128[extra_digits]);
 149
 150     // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 151     amount = __bid_recip_scale[extra_digits];
 152     CX2.w[0] = CT.w[2];
 153     CX2.w[1] = CT.w[3];
 154     if (amount >= 64) {
 155       CR.w[1] = 0;
 156       CR.w[0] = CX2.w[1] >> (amount - 64);
 157     } else {
 158       __shr_128 (CR, CX2, amount);
 159     }
 160
 161 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 162 #ifndef IEEE_ROUND_NEAREST
 163     if (rnd_mode == 0)
 164 #endif
 165       if (CR.w[0] & 1) {
 166         // check whether fractional part of initial_P/10^extra_digits is
 167         // exactly .5 this is the same as fractional part of
 168         // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
 169
 170         // get remainder
 171         if (amount >= 64) {
 172           remainder_h = CX2.w[0] | (CX2.w[1] << (128 - amount));
 173         } else
 174           remainder_h = CX2.w[0] << (64 - amount);
 175
 176         // test whether fractional part is 0
 177         if (!remainder_h
 178             && (CT.w[1] < __bid_reciprocals10_128[extra_digits].w[1]
 179                 || (CT.w[1] == __bid_reciprocals10_128[extra_digits].w[1]
 180                     && CT.w[0] < __bid_reciprocals10_128[extra_digits].w[0]))) {
 181           CR.w[0]--;
 182         }
 183       }
 184 #endif
 185
 186 #ifdef SET_STATUS_FLAGS
 187     status = INEXACT_EXCEPTION;
 188
 189     // get remainder
 190     if (amount >= 64) {
 191       REM_H.w[1] = (CX2.w[1] << (128 - amount));
 192       REM_H.w[0] = CX2.w[0];
 193     } else {
 194       REM_H.w[1] = CX2.w[0] << (64 - amount);
 195       REM_H.w[0] = 0;
 196     }
 197
 198     switch (rmode) {
 199     case ROUNDING_TO_NEAREST:
 200     case ROUNDING_TIES_AWAY:
 201       // test whether fractional part is 0
 202       if (REM_H.w[1] == 0x8000000000000000ull && !REM_H.w[0]
 203           && (CT.w[1] < __bid_reciprocals10_128[extra_digits].w[1]
 204               || (CT.w[1] == __bid_reciprocals10_128[extra_digits].w[1]
 205                   && CT.w[0] < __bid_reciprocals10_128[extra_digits].w[0])))
 206         status = EXACT_STATUS;
 207       break;
 208     case ROUNDING_DOWN:
 209     case ROUNDING_TO_ZERO:
 210       if (!(REM_H.w[1] | REM_H.w[0])
 211           && (CT.w[1] < __bid_reciprocals10_128[extra_digits].w[1]
 212               || (CT.w[1] == __bid_reciprocals10_128[extra_digits].w[1]
 213                   && CT.w[0] < __bid_reciprocals10_128[extra_digits].w[0])))
 214         status = EXACT_STATUS;
 215       break;
 216     default:
 217       // round up
 218       __add_carry_out (Stemp.w[0], CY64, CT.w[0],
 219                        __bid_reciprocals10_128[extra_digits].w[0]);
 220       __add_carry_in_out (Stemp.w[1], carry, CT.w[1],
 221                           __bid_reciprocals10_128[extra_digits].w[1], CY64);
 222       if (amount < 64) {
 223         C2N.w[1] = 0;
 224         C2N.w[0] = ((UINT64) 1) << amount;
 225         REM_H.w[0] = REM_H.w[1] >> (64 - amount);
 226         REM_H.w[1] = 0;
 227       } else {
 228         C2N.w[1] = ((UINT64) 1) << (amount - 64);
 229         C2N.w[0] = 0;
 230         REM_H.w[1] >>= (128 - amount);
 231       }
 232       REM_H.w[0] += carry;
 233       if (REM_H.w[0] < carry)
 234         REM_H.w[1]++;
 235       if (__unsigned_compare_ge_128 (REM_H, C2N))
 236         status = EXACT_STATUS;
 237     }
 238
 239     __set_status_flags (pfpsf, status);
 240
 241 #endif
 242     get_BID128_very_fast (&res, sign_x, exponent_y, CR);
 243     BID_RETURN (res);
 244   }
 245   if (total_digits < 0) {
 246     CR.w[1] = CR.w[0] = 0;
 247 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 248 #ifndef IEEE_ROUND_NEAREST
 249     rmode = rnd_mode;
 250     if (sign_x && (unsigned) (rmode - 1) < 2)
 251       rmode = 3 - rmode;
 252     if (rmode == ROUNDING_UP)
 253       CR.w[0] = 1;
 254 #endif
 255 #endif
 256 #ifdef SET_STATUS_FLAGS
 257     __set_status_flags (pfpsf, INEXACT_EXCEPTION);
 258 #endif
 259     get_BID128_very_fast (&res, sign_x, exponent_y, CR);
 260     BID_RETURN (res);
 261   }
 262   // else  more than 34 digits in coefficient
 263 #ifdef SET_STATUS_FLAGS
 264   __set_status_flags (pfpsf, INVALID_EXCEPTION);
 265 #endif
 266   res.w[1] = 0x7c00000000000000ull;
 267   res.w[0] = 0;
 268   BID_RETURN (res);
 269
 270 }