gcc/libgcc/config/libbid/bid64_mul.c

/* Copyright (C) 2007  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 2, or (at your option) any later
version.

In addition to the permissions in the GNU General Public License, the
Free Software Foundation gives you unlimited permission to link the
compiled version of this file into combinations with other programs,
and to distribute those combinations without any restriction coming
from the use of this file.  (The General Public License restrictions
do apply in other respects; for example, they cover modification of
the file, and distribution when not linked into a combine
executable.)

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.

You should have received a copy of the GNU General Public License
along with GCC; see the file COPYING.  If not, write to the Free
Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.  */

/*****************************************************************************
 *    BID64 multiply
 *****************************************************************************
 *
 *  Algorithm description:
 *
 *  if(number_digits(coefficient_x)+number_digits(coefficient_y) guaranteed
 *       below 16)
 *      return get_BID64(sign_x^sign_y, exponent_x + exponent_y - dec_bias,
 *                     coefficient_x*coefficient_y)
 *  else
 *      get long product: coefficient_x*coefficient_y
 *      determine number of digits to round off (extra_digits)
 *      rounding is performed as a 128x128-bit multiplication by 
 *         2^M[extra_digits]/10^extra_digits, followed by a shift
 *         M[extra_digits] is sufficiently large for required accuracy 
 *
 ****************************************************************************/

#include "bid_internal.h"

#if DECIMAL_CALL_BY_REFERENCE

void
bid64_mul (UINT64 * pres, UINT64 * px,
	   UINT64 *
	   py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	   _EXC_INFO_PARAM) {
  UINT64 x, y;
#else

UINT64
bid64_mul (UINT64 x,
	   UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
	   _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
#endif
  UINT128 P, PU, C128, Q_high, Q_low, Stemp;
  UINT64 sign_x, sign_y, coefficient_x, coefficient_y;
  UINT64 C64, remainder_h, carry, CY, res;
  UINT64 valid_x, valid_y;
  int_double tempx, tempy;
  int extra_digits, exponent_x, exponent_y, bin_expon_cx, bin_expon_cy,
    bin_expon_product;
  int rmode, digits_p, bp, amount, amount2, final_exponent, round_up;
  unsigned status, uf_status;

#if DECIMAL_CALL_BY_REFERENCE
#if !DECIMAL_GLOBAL_ROUNDING
  _IDEC_round rnd_mode = *prnd_mode;
#endif
  x = *px;
  y = *py;
#endif

  valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
  valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);

  // unpack arguments, check for NaN or Infinity
  if (!valid_x) {

#ifdef SET_STATUS_FLAGS
    if ((y & SNAN_MASK64) == SNAN_MASK64)	// y is sNaN
      __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
    // x is Inf. or NaN

    // test if x is NaN
    if ((x & NAN_MASK64) == NAN_MASK64) {
#ifdef SET_STATUS_FLAGS
      if ((x & SNAN_MASK64) == SNAN_MASK64)	// sNaN
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
      BID_RETURN (coefficient_x & QUIET_MASK64);
    }
    // x is Infinity?
    if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
      // check if y is 0
      if (((y & INFINITY_MASK64) != INFINITY_MASK64)
	  && !coefficient_y) {
#ifdef SET_STATUS_FLAGS
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
	// y==0 , return NaN
	BID_RETURN (NAN_MASK64);
      }
      // check if y is NaN
      if ((y & NAN_MASK64) == NAN_MASK64)
	// y==NaN , return NaN
	BID_RETURN (coefficient_y & QUIET_MASK64);
      // otherwise return +/-Inf
      BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64);
    }
    // x is 0
    if (((y & INFINITY_MASK64) != INFINITY_MASK64)) {
      if ((y & SPECIAL_ENCODING_MASK64) == SPECIAL_ENCODING_MASK64)
	exponent_y = ((UINT32) (y >> 51)) & 0x3ff;
      else
	exponent_y = ((UINT32) (y >> 53)) & 0x3ff;
      sign_y = y & 0x8000000000000000ull;

      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;
      BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
    }
  }
  if (!valid_y) {
    // y is Inf. or NaN

    // test if y is NaN
    if ((y & NAN_MASK64) == NAN_MASK64) {
#ifdef SET_STATUS_FLAGS
      if ((y & SNAN_MASK64) == SNAN_MASK64)	// sNaN
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
      BID_RETURN (coefficient_y & QUIET_MASK64);
    }
    // y is Infinity?
    if ((y & INFINITY_MASK64) == INFINITY_MASK64) {
      // check if x is 0
      if (!coefficient_x) {
	__set_status_flags (pfpsf, INVALID_EXCEPTION);
	// x==0, return NaN
	BID_RETURN (NAN_MASK64);
      }
      // otherwise return +/-Inf
      BID_RETURN (((x ^ y) & 0x8000000000000000ull) | INFINITY_MASK64);
    }
    // y is 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;
    BID_RETURN ((sign_x ^ sign_y) | (((UINT64) exponent_x) << 53));
  }
  //--- 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;

    res =
      get_BID64_small_mantissa (sign_x ^ sign_y,
				exponent_x + exponent_y -
				DECIMAL_EXPONENT_BIAS, C64, rnd_mode,
				pfpsf);
    BID_RETURN (res);
  } else {
    uf_status = 0;
    // 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;

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
    rmode = rnd_mode;
    if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
      rmode = 3 - rmode;
#else
    rmode = 0;
#endif
#else
    rmode = 0;
#endif

    round_up = 0;
    if (((unsigned) final_exponent) >= 3 * 256) {
      if (final_exponent < 0) {
	// underflow
	if (final_exponent + 16 < 0) {
	  res = sign_x ^ sign_y;
	  __set_status_flags (pfpsf,
			      UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
	  if (rmode == ROUNDING_UP)
	    res |= 1;
	  BID_RETURN (res);
	}

	uf_status = UNDERFLOW_EXCEPTION;
	if (final_exponent == -1) {
	  __add_128_64 (PU, P, round_const_table[rmode][extra_digits]);
	  if (__unsigned_compare_ge_128
	      (PU, power10_table_128[extra_digits + 16]))
	    uf_status = 0;
	}
	extra_digits -= final_exponent;
	final_exponent = 0;

	if (extra_digits > 17) {
	  __mul_128x128_full (Q_high, Q_low, P, reciprocals10_128[16]);

	  amount = recip_scale[16];
	  __shr_128 (P, Q_high, amount);

	  // get sticky bits
	  amount2 = 64 - amount;
	  remainder_h = 0;
	  remainder_h--;
	  remainder_h >>= amount2;
	  remainder_h = remainder_h & Q_high.w[0];

	  extra_digits -= 16;
	  if (remainder_h || (Q_low.w[1] > reciprocals10_128[16].w[1]
			      || (Q_low.w[1] ==
				  reciprocals10_128[16].w[1]
				  && Q_low.w[0] >=
				  reciprocals10_128[16].w[0]))) {
	    round_up = 1;
	    __set_status_flags (pfpsf,
				UNDERFLOW_EXCEPTION |
				INEXACT_EXCEPTION);
	    P.w[0] = (P.w[0] << 3) + (P.w[0] << 1);
	    P.w[0] |= 1;
	    extra_digits++;
	  }
	}
      } else {
	res =
	  fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent,
				   1000000000000000ull, rnd_mode,
				   pfpsf);
	BID_RETURN (res);
      }
    }


    if (extra_digits > 0) {
      // will divide by 10^(digits_p - 16)

      // add a constant to P, depending on rounding mode
      // 0.5*10^(digits_p - 16) for round-to-nearest
      __add_128_64 (P, P, round_const_table[rmode][extra_digits]);

      // get P*(2^M[extra_digits])/10^extra_digits
      __mul_128x128_full (Q_high, Q_low, P,
			  reciprocals10_128[extra_digits]);

      // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
      amount = recip_scale[extra_digits];
      __shr_128 (C128, Q_high, amount);

      C64 = __low_64 (C128);

#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
#ifndef IEEE_ROUND_NEAREST
      if (rmode == 0)	//ROUNDING_TO_NEAREST
#endif
	if ((C64 & 1) && !round_up) {
	  // check whether fractional part of initial_P/10^extra_digits 
	  // is exactly .5
	  // this is the same as fractional part of 
	  // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero

	  // get remainder
	  remainder_h = Q_high.w[0] << (64 - amount);

	  // test whether fractional part is 0
	  if (!remainder_h
	      && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		  || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		      && Q_low.w[0] <
		      reciprocals10_128[extra_digits].w[0]))) {
	    C64--;
	  }
	}
#endif

#ifdef SET_STATUS_FLAGS
      status = INEXACT_EXCEPTION | uf_status;

      // get remainder
      remainder_h = Q_high.w[0] << (64 - amount);

      switch (rmode) {
      case ROUNDING_TO_NEAREST:
      case ROUNDING_TIES_AWAY:
	// test whether fractional part is 0
	if (remainder_h == 0x8000000000000000ull
	    && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		|| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		    && Q_low.w[0] <
		    reciprocals10_128[extra_digits].w[0])))
	  status = EXACT_STATUS;
	break;
      case ROUNDING_DOWN:
      case ROUNDING_TO_ZERO:
	if (!remainder_h
	    && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
		|| (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
		    && Q_low.w[0] <
		    reciprocals10_128[extra_digits].w[0])))
	  status = EXACT_STATUS;
	break;
      default:
	// round up
	__add_carry_out (Stemp.w[0], CY, Q_low.w[0],
			 reciprocals10_128[extra_digits].w[0]);
	__add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
			    reciprocals10_128[extra_digits].w[1], CY);
	if ((remainder_h >> (64 - amount)) + carry >=
	    (((UINT64) 1) << amount))
	  status = EXACT_STATUS;
      }

      __set_status_flags (pfpsf, status);
#endif

      // convert to BID and return
      res =
	fast_get_BID64_check_OF (sign_x ^ sign_y, final_exponent, C64,
				 rmode, pfpsf);
      BID_RETURN (res);
    }
    // go to convert_format and exit
    C64 = __low_64 (P);
    res =
      get_BID64 (sign_x ^ sign_y,
		 exponent_x + exponent_y - DECIMAL_EXPONENT_BIAS, C64,
		 rmode, pfpsf);
    BID_RETURN (res);
  }
}