--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libgcc/config/libbid/bid64_mul.c	Fri Jul 17 14:47:48 2009 +0900
@@ -0,0 +1,374 @@
+/* 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 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);
+  }
+}