Subversion Repositories openrisc

/* Copyright (C) 2007, 2009  Free Software Foundation, Inc.

This file is part of GCC.

GCC is free software; you can redistribute it and/or modify it under

the terms of the GNU General Public License as published by the Free

Software Foundation; either version 3, or (at your option) any later

version.

GCC is distributed in the hope that it will be useful, but WITHOUT ANY

WARRANTY; without even the implied warranty of MERCHANTABILITY or

FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

for more details.

Under Section 7 of GPL version 3, you are granted additional

permissions described in the GCC Runtime Library Exception, version

3.1, as published by the Free Software Foundation.

You should have received a copy of the GNU General Public License and

a copy of the GCC Runtime Library Exception along with this program;

see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see

<http://www.gnu.org/licenses/>.  */

/*****************************************************************************

 *    BID64 fma

 *****************************************************************************

 *

 *  Algorithm description:

 *

 *  if multiplication is guranteed exact (short coefficients)

 *     call the unpacked arg. equivalent of bid64_add(x*y, z)

 *  else

 *     get full coefficient_x*coefficient_y product

 *     call subroutine to perform addition of 64-bit argument

 *                                         to 128-bit product

 *

 ****************************************************************************/

#include "bid_inline_add.h"

#if DECIMAL_CALL_BY_REFERENCE

extern void bid64_mul (UINT64 * pres, UINT64 * px,

                       UINT64 *

                       py _RND_MODE_PARAM _EXC_FLAGS_PARAM

                       _EXC_MASKS_PARAM _EXC_INFO_PARAM);

#else

extern UINT64 bid64_mul (UINT64 x,

                         UINT64 y _RND_MODE_PARAM

                         _EXC_FLAGS_PARAM _EXC_MASKS_PARAM

                         _EXC_INFO_PARAM);

#endif

#if DECIMAL_CALL_BY_REFERENCE

void

bid64_fma (UINT64 * pres, UINT64 * px, UINT64 * py,

           UINT64 *

           pz _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM

           _EXC_INFO_PARAM) {

  UINT64 x, y, z;

#else

UINT64

bid64_fma (UINT64 x, UINT64 y,

           UINT64 z _RND_MODE_PARAM _EXC_FLAGS_PARAM

           _EXC_MASKS_PARAM _EXC_INFO_PARAM) {

#endif

  UINT128 P, PU, CT, CZ;

  UINT64 sign_x, sign_y, coefficient_x, coefficient_y, sign_z,

    coefficient_z;

  UINT64 C64, remainder_y, res;

  UINT64 CYh, CY0L, T, valid_x, valid_y, valid_z;

  int_double tempx, tempy;

  int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,

    bin_expon_product, rmode;

  int digits_p, bp, final_exponent, exponent_z, digits_z, ez, ey,

    scale_z, uf_status;

#if DECIMAL_CALL_BY_REFERENCE

#if !DECIMAL_GLOBAL_ROUNDING

  _IDEC_round rnd_mode = *prnd_mode;

#endif

  x = *px;

  y = *py;

  z = *pz;

#endif

  valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);

  valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);

  valid_z = unpack_BID64 (&sign_z, &exponent_z, &coefficient_z, z);

  // unpack arguments, check for NaN, Infinity, or 0

  if (!valid_x || !valid_y || !valid_z) {

    if ((y & MASK_NAN) == MASK_NAN) {   // y is NAN

      // if x = {0, f, inf, NaN}, y = NaN, z = {0, f, inf, NaN} then res = Q (y)

      // check first for non-canonical NaN payload

      y = y & 0xfe03ffffffffffffull;    // clear G6-G12

      if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {

        y = y & 0xfe00000000000000ull;  // clear G6-G12 and the payload bits

      }

      if ((y & MASK_SNAN) == MASK_SNAN) {       // y is SNAN

        // set invalid flag

        *pfpsf |= INVALID_EXCEPTION;

        // return quiet (y)

        res = y & 0xfdffffffffffffffull;

      } else {  // y is QNaN

        // return y

        res = y;

        // if z = SNaN or x = SNaN signal invalid exception

        if ((z & MASK_SNAN) == MASK_SNAN

            || (x & MASK_SNAN) == MASK_SNAN) {

          // set invalid flag

          *pfpsf |= INVALID_EXCEPTION;

        }

      }

      BID_RETURN (res)

    } else if ((z & MASK_NAN) == MASK_NAN) {    // z is NAN

      // if x = {0, f, inf, NaN}, y = {0, f, inf}, z = NaN then res = Q (z)

      // check first for non-canonical NaN payload

      z = z & 0xfe03ffffffffffffull;    // clear G6-G12

      if ((z & 0x0003ffffffffffffull) > 999999999999999ull) {

        z = z & 0xfe00000000000000ull;  // clear G6-G12 and the payload bits

      }

      if ((z & MASK_SNAN) == MASK_SNAN) {       // z is SNAN

        // set invalid flag

        *pfpsf |= INVALID_EXCEPTION;

        // return quiet (z)

        res = z & 0xfdffffffffffffffull;

      } else {  // z is QNaN

        // return z

        res = z;

        // if x = SNaN signal invalid exception

        if ((x & MASK_SNAN) == MASK_SNAN) {

          // set invalid flag

          *pfpsf |= INVALID_EXCEPTION;

        }

      }

      BID_RETURN (res)

    } else if ((x & MASK_NAN) == MASK_NAN) {    // x is NAN

      // if x = NaN, y = {0, f, inf}, z = {0, f, inf} then res = Q (x)

      // check first for non-canonical NaN payload

      x = x & 0xfe03ffffffffffffull;    // clear G6-G12

      if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {

        x = x & 0xfe00000000000000ull;  // clear G6-G12 and the payload bits

      }

      if ((x & MASK_SNAN) == MASK_SNAN) {       // x is SNAN

        // set invalid flag

        *pfpsf |= INVALID_EXCEPTION;

        // return quiet (x)

        res = x & 0xfdffffffffffffffull;

      } else {  // x is QNaN

        // return x

        res = x;        // clear out G[6]-G[16]

      }

      BID_RETURN (res)

    }

    if (!valid_x) {

      // x is Inf. or 0

      // x is Infinity?

      if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {

        // check if y is 0

        if (!coefficient_y) {

          // y==0, return NaN

#ifdef SET_STATUS_FLAGS

          if ((z & 0x7e00000000000000ull) != 0x7c00000000000000ull)

            __set_status_flags (pfpsf, INVALID_EXCEPTION);

#endif

          BID_RETURN (0x7c00000000000000ull);

        }

        // test if z is Inf of oposite sign

        if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull)

            && (((x ^ y) ^ z) & 0x8000000000000000ull)) {

          // return NaN 

#ifdef SET_STATUS_FLAGS

          __set_status_flags (pfpsf, INVALID_EXCEPTION);

#endif

          BID_RETURN (0x7c00000000000000ull);

        }

        // otherwise return +/-Inf

        BID_RETURN (((x ^ y) & 0x8000000000000000ull) |

                    0x7800000000000000ull);

      }

      // x is 0

      if (((y & 0x7800000000000000ull) != 0x7800000000000000ull)

          && ((z & 0x7800000000000000ull) != 0x7800000000000000ull)) {

        if (coefficient_z) {

          exponent_y = exponent_x - DECIMAL_EXPONENT_BIAS + exponent_y;

          sign_z = z & 0x8000000000000000ull;

          if (exponent_y >= exponent_z)

            BID_RETURN (z);

          res =

            add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z,

                        &rnd_mode, pfpsf);

          BID_RETURN (res);

        }

      }

    }

    if (!valid_y) {

      // y is Inf. or 0

      // y is Infinity?

      if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {

        // check if x is 0

        if (!coefficient_x) {

          // y==0, return NaN

#ifdef SET_STATUS_FLAGS

          __set_status_flags (pfpsf, INVALID_EXCEPTION);

#endif

          BID_RETURN (0x7c00000000000000ull);

        }

        // test if z is Inf of oposite sign

        if (((z & 0x7c00000000000000ull) == 0x7800000000000000ull)

            && (((x ^ y) ^ z) & 0x8000000000000000ull)) {

#ifdef SET_STATUS_FLAGS

          __set_status_flags (pfpsf, INVALID_EXCEPTION);

#endif

          // return NaN

          BID_RETURN (0x7c00000000000000ull);

        }

        // otherwise return +/-Inf

        BID_RETURN (((x ^ y) & 0x8000000000000000ull) |

                    0x7800000000000000ull);

      }

      // y is 0 

      if (((z & 0x7800000000000000ull) != 0x7800000000000000ull)) {

        if (coefficient_z) {

          exponent_y += exponent_x - DECIMAL_EXPONENT_BIAS;

          sign_z = z & 0x8000000000000000ull;

          if (exponent_y >= exponent_z)

            BID_RETURN (z);

          res =

            add_zero64 (exponent_y, sign_z, exponent_z, coefficient_z,

                        &rnd_mode, pfpsf);

          BID_RETURN (res);

        }

      }

    }

    if (!valid_z) {

      // y is Inf. or 0

      // test if y is NaN/Inf

      if ((z & 0x7800000000000000ull) == 0x7800000000000000ull) {

        BID_RETURN (coefficient_z & QUIET_MASK64);

      }

      // z is 0, return x*y

      if ((!coefficient_x) || (!coefficient_y)) {

        //0+/-0

        exponent_x += exponent_y - DECIMAL_EXPONENT_BIAS;

        if (exponent_x > DECIMAL_MAX_EXPON_64)

          exponent_x = DECIMAL_MAX_EXPON_64;

        else if (exponent_x < 0)

          exponent_x = 0;

        if (exponent_x <= exponent_z)

          res = ((UINT64) exponent_x) << 53;

        else

          res = ((UINT64) exponent_z) << 53;

        if ((sign_x ^ sign_y) == sign_z)

          res |= sign_z;

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY

#ifndef IEEE_ROUND_NEAREST

        else if (rnd_mode == ROUNDING_DOWN)

          res |= 0x8000000000000000ull;

#endif

#endif

        BID_RETURN (res);

      }

    }

  }

  /* get binary coefficients of x and y */

  //--- get number of bits in the coefficients of x and y ---

  // version 2 (original)

  tempx.d = (double) coefficient_x;

  bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52);

  tempy.d = (double) coefficient_y;

  bin_expon_cy = ((tempy.i & MASK_BINARY_EXPONENT) >> 52);

  // magnitude estimate for coefficient_x*coefficient_y is 

  //        2^(unbiased_bin_expon_cx + unbiased_bin_expon_cx)

  bin_expon_product = bin_expon_cx + bin_expon_cy;

  // check if coefficient_x*coefficient_y<2^(10*k+3)

  // equivalent to unbiased_bin_expon_cx + unbiased_bin_expon_cx < 10*k+1

  if (bin_expon_product < UPPER_EXPON_LIMIT + 2 * BINARY_EXPONENT_BIAS) {

    //  easy multiply

    C64 = coefficient_x * coefficient_y;

    final_exponent = exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS;

    if ((final_exponent > 0) || (!coefficient_z)) {

      res =

        get_add64 (sign_x ^ sign_y,

                   final_exponent, C64, sign_z, exponent_z, coefficient_z, rnd_mode, pfpsf);

      BID_RETURN (res);

    } else {

      P.w[0] = C64;

      P.w[1] = 0;

      extra_digits = 0;

    }

  } else {

    if (!coefficient_z) {

#if DECIMAL_CALL_BY_REFERENCE

      bid64_mul (&res, px,

                 py _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG

                 _EXC_INFO_ARG);

#else

      res =

        bid64_mul (x,

                   y _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG

                   _EXC_INFO_ARG);

#endif

      BID_RETURN (res);

    }

    // get 128-bit product: coefficient_x*coefficient_y

    __mul_64x64_to_128 (P, coefficient_x, coefficient_y);

    // tighten binary range of P:  leading bit is 2^bp

    // unbiased_bin_expon_product <= bp <= unbiased_bin_expon_product+1

    bin_expon_product -= 2 * BINARY_EXPONENT_BIAS;

    __tight_bin_range_128 (bp, P, bin_expon_product);

    // get number of decimal digits in the product

    digits_p = estimate_decimal_digits[bp];

    if (!(__unsigned_compare_gt_128 (power10_table_128[digits_p], P)))

      digits_p++;       // if power10_table_128[digits_p] <= P

    // determine number of decimal digits to be rounded out

    extra_digits = digits_p - MAX_FORMAT_DIGITS;

    final_exponent =

      exponent_x + exponent_y + extra_digits - DECIMAL_EXPONENT_BIAS;

  }

  if (((unsigned) final_exponent) >= 3 * 256) {

    if (final_exponent < 0) {

      //--- get number of bits in the coefficients of z  ---

      tempx.d = (double) coefficient_z;

      bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;

      // get number of decimal digits in the coeff_x

      digits_z = estimate_decimal_digits[bin_expon_cx];

      if (coefficient_z >= power10_table_128[digits_z].w[0])

        digits_z++;

      // underflow

      if ((final_exponent + 16 < 0)

          || (exponent_z + digits_z > 33 + final_exponent)) {

        res =

          BID_normalize (sign_z, exponent_z, coefficient_z,

                         sign_x ^ sign_y, 1, rnd_mode, pfpsf);

        BID_RETURN (res);

      }

      ez = exponent_z + digits_z - 16;

      if (ez < 0)

        ez = 0;

      scale_z = exponent_z - ez;

      coefficient_z *= power10_table_128[scale_z].w[0];

      ey = final_exponent - extra_digits;

      extra_digits = ez - ey;

      if (extra_digits > 33) {

        res =

          BID_normalize (sign_z, exponent_z, coefficient_z,

                         sign_x ^ sign_y, 1, rnd_mode, pfpsf);

        BID_RETURN (res);

      }

      //else  // extra_digits<=32

      if (extra_digits > 17) {

        CYh = __truncate (P, 16);

        // get remainder

        T = power10_table_128[16].w[0];

        __mul_64x64_to_64 (CY0L, CYh, T);

        remainder_y = P.w[0] - CY0L;

        extra_digits -= 16;

        P.w[0] = CYh;

        P.w[1] = 0;

      } else

        remainder_y = 0;

      // align coeff_x, CYh

      __mul_64x64_to_128 (CZ, coefficient_z,

                          power10_table_128[extra_digits].w[0]);

      if (sign_z == (sign_y ^ sign_x)) {

        __add_128_128 (CT, CZ, P);

        if (__unsigned_compare_ge_128

            (CT, power10_table_128[16 + extra_digits])) {

          extra_digits++;

          ez++;

        }

      } else {

        if (remainder_y && (__unsigned_compare_ge_128 (CZ, P))) {

          P.w[0]++;

          if (!P.w[0])

            P.w[1]++;

        }

        __sub_128_128 (CT, CZ, P);

        if (((SINT64) CT.w[1]) < 0) {

          sign_z = sign_y ^ sign_x;

          CT.w[0] = 0 - CT.w[0];

          CT.w[1] = 0 - CT.w[1];

          if (CT.w[0])

            CT.w[1]--;

        } else if(!(CT.w[1]|CT.w[0]))

                sign_z = (rnd_mode!=ROUNDING_DOWN)? 0: 0x8000000000000000ull;

        if (ez

            &&

            (__unsigned_compare_gt_128

             (power10_table_128[15 + extra_digits], CT))) {

          extra_digits--;

          ez--;

        }

      }

#ifdef SET_STATUS_FLAGS

      uf_status = 0;

      if ((!ez)

          &&

          __unsigned_compare_gt_128 (power10_table_128

                                     [extra_digits + 15], CT)) {

        rmode = rnd_mode;

        if (sign_z && (unsigned) (rmode - 1) < 2)

          rmode = 3 - rmode;

        //__add_128_64(PU, CT, round_const_table[rmode][extra_digits]);

        PU = power10_table_128[extra_digits + 15];

        PU.w[0]--;

        if (__unsigned_compare_gt_128 (PU, CT)

            || (rmode == ROUNDING_DOWN)

            || (rmode == ROUNDING_TO_ZERO))

          uf_status = UNDERFLOW_EXCEPTION;

        else if (extra_digits < 2) {

          if ((rmode == ROUNDING_UP)) {

            if (!extra_digits)

              uf_status = UNDERFLOW_EXCEPTION;

            else {

              if (remainder_y && (sign_z != (sign_y ^ sign_x)))

                remainder_y = power10_table_128[16].w[0] - remainder_y;

              if (power10_table_128[15].w[0] > remainder_y)

                uf_status = UNDERFLOW_EXCEPTION;

            }

          } else        // RN or RN_away

          {

            if (remainder_y && (sign_z != (sign_y ^ sign_x)))

              remainder_y = power10_table_128[16].w[0] - remainder_y;

            if (!extra_digits) {

              remainder_y += round_const_table[rmode][15];

              if (remainder_y < power10_table_128[16].w[0])

                uf_status = UNDERFLOW_EXCEPTION;

            } else {

              if (remainder_y < round_const_table[rmode][16])

                uf_status = UNDERFLOW_EXCEPTION;

            }

          }

          //__set_status_flags (pfpsf, uf_status);

        }

      }

#endif

      res =

        __bid_full_round64_remainder (sign_z, ez - extra_digits, CT,

                                      extra_digits, remainder_y,

                                      rnd_mode, pfpsf, uf_status);

      BID_RETURN (res);

    } else {

      if ((sign_z == (sign_x ^ sign_y))

          || (final_exponent > 3 * 256 + 15)) {

        res =

          fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,

                                   1000000000000000ull, rnd_mode,

                                   pfpsf);

        BID_RETURN (res);

      }

    }

  }

  if (extra_digits > 0) {

    res =

      get_add128 (sign_z, exponent_z, coefficient_z, sign_x ^ sign_y,

                  final_exponent, P, extra_digits, rnd_mode, pfpsf);

    BID_RETURN (res);

  }

  // go to convert_format and exit

  else {

    C64 = __low_64 (P);

    res =

      get_add64 (sign_x ^ sign_y,

                 exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,

                 sign_z, exponent_z, coefficient_z,

                 rnd_mode, pfpsf);

    BID_RETURN (res);

  }

}