gcc/libgcc/config/libbid/bid128_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.  */

#include "bid_internal.h"

#if DECIMAL_CALL_BY_REFERENCE
void
bid64dq_mul (UINT64 * pres, UINT64 * px, UINT128 * py
	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	     _EXC_INFO_PARAM) {
  UINT64 x = *px;
#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;
#endif
#else
UINT64
bid64dq_mul (UINT64 x, UINT128 y
	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	     _EXC_INFO_PARAM) {
#endif
  UINT64 res = 0xbaddbaddbaddbaddull;
  UINT128 x1;

#if DECIMAL_CALL_BY_REFERENCE
  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid64qq_mul (&res, &x1, py
	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	       _EXC_INFO_ARG);
#else
  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  res = bid64qq_mul (x1, y
		     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		     _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}


#if DECIMAL_CALL_BY_REFERENCE
void
bid64qd_mul (UINT64 * pres, UINT128 * px, UINT64 * py
	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	     _EXC_INFO_PARAM) {
  UINT64 y = *py;
#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;
#endif
#else
UINT64
bid64qd_mul (UINT128 x, UINT64 y
	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	     _EXC_INFO_PARAM) {
#endif
  UINT64 res = 0xbaddbaddbaddbaddull;
  UINT128 y1;

#if DECIMAL_CALL_BY_REFERENCE
  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid64qq_mul (&res, px, &y1
	       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	       _EXC_INFO_ARG);
#else
  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  res = bid64qq_mul (x, y1
		     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		     _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}


#if DECIMAL_CALL_BY_REFERENCE
void
bid64qq_mul (UINT64 * pres, UINT128 * px, UINT128 * py
	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	     _EXC_INFO_PARAM) {
  UINT128 x = *px, y = *py;
#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;
#endif
#else
UINT64
bid64qq_mul (UINT128 x, UINT128 y
	     _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	     _EXC_INFO_PARAM) {
#endif

  UINT128 z = { {0x0000000000000000ull, 0x5ffe000000000000ull}
  };
  UINT64 res = 0xbaddbaddbaddbaddull;
  UINT64 x_sign, y_sign, p_sign;
  UINT64 x_exp, y_exp, p_exp;
  int true_p_exp;
  UINT128 C1, C2;

  BID_SWAP128 (z);
  // skip cases where at least one operand is NaN or infinity
  if (!(((x.w[HIGH_128W] & MASK_NAN) == MASK_NAN) ||
	((y.w[HIGH_128W] & MASK_NAN) == MASK_NAN) ||
	((x.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF) ||
	((y.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF))) {
    // x, y are 0 or f but not inf or NaN => unpack the arguments and check
    // for non-canonical values

    x_sign = x.w[HIGH_128W] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
    C1.w[1] = x.w[HIGH_128W] & MASK_COEFF;
    C1.w[0] = x.w[LOW_128W];
    // check for non-canonical values - treated as zero
    if ((x.w[HIGH_128W] & 0x6000000000000000ull) ==
	0x6000000000000000ull) {
      // G0_G1=11 => non-canonical
      x_exp = (x.w[HIGH_128W] << 2) & MASK_EXP;	// biased and shifted left 49 bits
      C1.w[1] = 0;	// significand high
      C1.w[0] = 0;	// significand low
    } else {	// G0_G1 != 11
      x_exp = x.w[HIGH_128W] & MASK_EXP;	// biased and shifted left 49 bits
      if (C1.w[1] > 0x0001ed09bead87c0ull ||
	  (C1.w[1] == 0x0001ed09bead87c0ull &&
	   C1.w[0] > 0x378d8e63ffffffffull)) {
	// x is non-canonical if coefficient is larger than 10^34 -1
	C1.w[1] = 0;
	C1.w[0] = 0;
      } else {	// canonical          
	;
      }
    }
    y_sign = y.w[HIGH_128W] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
    C2.w[1] = y.w[HIGH_128W] & MASK_COEFF;
    C2.w[0] = y.w[LOW_128W];
    // check for non-canonical values - treated as zero
    if ((y.w[HIGH_128W] & 0x6000000000000000ull) ==
	0x6000000000000000ull) {
      // G0_G1=11 => non-canonical
      y_exp = (y.w[HIGH_128W] << 2) & MASK_EXP;	// biased and shifted left 49 bits
      C2.w[1] = 0;	// significand high
      C2.w[0] = 0;	// significand low 
    } else {	// G0_G1 != 11
      y_exp = y.w[HIGH_128W] & MASK_EXP;	// biased and shifted left 49 bits
      if (C2.w[1] > 0x0001ed09bead87c0ull ||
	  (C2.w[1] == 0x0001ed09bead87c0ull &&
	   C2.w[0] > 0x378d8e63ffffffffull)) {
	// y is non-canonical if coefficient is larger than 10^34 -1
	C2.w[1] = 0;
	C2.w[0] = 0;
      } else {	// canonical
	;
      }
    }
    p_sign = x_sign ^ y_sign;	// sign of the product

    true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
    // true_p_exp, p_exp are used only for 0 * 0, 0 * f, or f * 0 
    if (true_p_exp < -398)
      p_exp = 0;	// cannot be less than EXP_MIN
    else if (true_p_exp > 369)
      p_exp = (UINT64) (369 + 398) << 53;	// cannot be more than EXP_MAX
    else
      p_exp = (UINT64) (true_p_exp + 398) << 53;

    if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
	(C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
      // x = 0 or y = 0
      // the result is 0
      res = p_sign | p_exp;	// preferred exponent in [EXP_MIN, EXP_MAX]
      BID_RETURN (res)
    }	// else continue
  }
  // swap x and y - ensure that a NaN in x has 'higher precedence' than one in y
#if DECIMAL_CALL_BY_REFERENCE
  bid64qqq_fma (&res, &y, &x, &z
		_RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		_EXC_INFO_ARG);
#else
  res = bid64qqq_fma (y, x, z
		      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		      _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}


#if DECIMAL_CALL_BY_REFERENCE
void
bid128dd_mul (UINT128 * pres, UINT64 * px, UINT64 * py
	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
  UINT64 x = *px, y = *py;
#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;
#endif
#else
UINT128
bid128dd_mul (UINT64 x, UINT64 y
	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
#endif
  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
  };
  UINT128 x1, y1;

#if DECIMAL_CALL_BY_REFERENCE
  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid128_mul (&res, &x1, &y1
	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	      _EXC_INFO_ARG);
#else
  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  res = bid128_mul (x1, y1
		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		    _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}


#if DECIMAL_CALL_BY_REFERENCE
void
bid128dq_mul (UINT128 * pres, UINT64 * px, UINT128 * py
	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
  UINT64 x = *px;
#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;
#endif
#else
UINT128
bid128dq_mul (UINT64 x, UINT128 y
	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
#endif
  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
  };
  UINT128 x1;

#if DECIMAL_CALL_BY_REFERENCE
  bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid128_mul (&res, &x1, py
	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	      _EXC_INFO_ARG);
#else
  x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  res = bid128_mul (x1, y
		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		    _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}


#if DECIMAL_CALL_BY_REFERENCE
void
bid128qd_mul (UINT128 * pres, UINT128 * px, UINT64 * py
	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
  UINT64 y = *py;
#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;
#endif
#else
UINT128
bid128qd_mul (UINT128 x, UINT64 y
	      _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	      _EXC_INFO_PARAM) {
#endif
  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
  };
  UINT128 y1;

#if DECIMAL_CALL_BY_REFERENCE
  bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  bid128_mul (&res, px, &y1
	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	      _EXC_INFO_ARG);
#else
  y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  res = bid128_mul (x, y1
		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		    _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}


// bid128_mul stands for bid128qq_mul
#if DECIMAL_CALL_BY_REFERENCE
void
bid128_mul (UINT128 * pres, UINT128 * px,
	    UINT128 *
	    py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	    _EXC_INFO_PARAM) {
  UINT128 x = *px, y = *py;

#if !DECIMAL_GLOBAL_ROUNDING
  unsigned int rnd_mode = *prnd_mode;

#endif
#else
UINT128
bid128_mul (UINT128 x,
	    UINT128 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
	    _EXC_INFO_PARAM) {

#endif
  UINT128 z = { {0x0000000000000000ull, 0x5ffe000000000000ull}
  };
  UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
  };
  UINT64 x_sign, y_sign, p_sign;
  UINT64 x_exp, y_exp, p_exp;
  int true_p_exp;
  UINT128 C1, C2;

  BID_SWAP128 (x);
  BID_SWAP128 (y);
  // skip cases where at least one operand is NaN or infinity
  if (!(((x.w[1] & MASK_NAN) == MASK_NAN) ||
	((y.w[1] & MASK_NAN) == MASK_NAN) ||
	((x.w[1] & MASK_ANY_INF) == MASK_INF) ||
	((y.w[1] & MASK_ANY_INF) == MASK_INF))) {
    // x, y are 0 or f but not inf or NaN => unpack the arguments and check
    // for non-canonical values

    x_sign = x.w[1] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
    C1.w[1] = x.w[1] & MASK_COEFF;
    C1.w[0] = x.w[0];
    // check for non-canonical values - treated as zero
    if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
      // G0_G1=11 => non-canonical
      x_exp = (x.w[1] << 2) & MASK_EXP;	// biased and shifted left 49 bits
      C1.w[1] = 0;	// significand high
      C1.w[0] = 0;	// significand low
    } else {	// G0_G1 != 11
      x_exp = x.w[1] & MASK_EXP;	// biased and shifted left 49 bits
      if (C1.w[1] > 0x0001ed09bead87c0ull ||
	  (C1.w[1] == 0x0001ed09bead87c0ull &&
	   C1.w[0] > 0x378d8e63ffffffffull)) {
	// x is non-canonical if coefficient is larger than 10^34 -1
	C1.w[1] = 0;
	C1.w[0] = 0;
      } else {	// canonical          
	;
      }
    }
    y_sign = y.w[1] & MASK_SIGN;	// 0 for positive, MASK_SIGN for negative
    C2.w[1] = y.w[1] & MASK_COEFF;
    C2.w[0] = y.w[0];
    // check for non-canonical values - treated as zero
    if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
      // G0_G1=11 => non-canonical
      y_exp = (y.w[1] << 2) & MASK_EXP;	// biased and shifted left 49 bits
      C2.w[1] = 0;	// significand high
      C2.w[0] = 0;	// significand low 
    } else {	// G0_G1 != 11
      y_exp = y.w[1] & MASK_EXP;	// biased and shifted left 49 bits
      if (C2.w[1] > 0x0001ed09bead87c0ull ||
	  (C2.w[1] == 0x0001ed09bead87c0ull &&
	   C2.w[0] > 0x378d8e63ffffffffull)) {
	// y is non-canonical if coefficient is larger than 10^34 -1
	C2.w[1] = 0;
	C2.w[0] = 0;
      } else {	// canonical
	;
      }
    }
    p_sign = x_sign ^ y_sign;	// sign of the product

    true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
    // true_p_exp, p_exp are used only for 0 * 0, 0 * f, or f * 0 
    if (true_p_exp < -6176)
      p_exp = 0;	// cannot be less than EXP_MIN
    else if (true_p_exp > 6111)
      p_exp = (UINT64) (6111 + 6176) << 49;	// cannot be more than EXP_MAX
    else
      p_exp = (UINT64) (true_p_exp + 6176) << 49;

    if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
	(C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
      // x = 0 or y = 0
      // the result is 0
      res.w[1] = p_sign | p_exp;	// preferred exponent in [EXP_MIN, EXP_MAX]
      res.w[0] = 0x0;
      BID_SWAP128 (res);
      BID_RETURN (res)
    }	// else continue
  }

  BID_SWAP128 (x);
  BID_SWAP128 (y);
  BID_SWAP128 (z);
  // swap x and y - ensure that a NaN in x has 'higher precedence' than one in y
#if DECIMAL_CALL_BY_REFERENCE
  bid128_fma (&res, &y, &x, &z
	      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
	      _EXC_INFO_ARG);
#else
  res = bid128_fma (y, x, z
		    _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
		    _EXC_INFO_ARG);
#endif
  BID_RETURN (res);
}