libgcc/config/libbid/bid128_mul.c

   1 /* Copyright (C) 2007-2020 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 3, or (at your option) any later
   8 version.
   9
  10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13 for more details.
  14
  15 Under Section 7 of GPL version 3, you are granted additional
  16 permissions described in the GCC Runtime Library Exception, version
  17 3.1, as published by the Free Software Foundation.
  18
  19 You should have received a copy of the GNU General Public License and
  20 a copy of the GCC Runtime Library Exception along with this program;
  21 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  22 <http://www.gnu.org/licenses/>.  */
  23
  24 #include "bid_internal.h"
  25
  26 #if DECIMAL_CALL_BY_REFERENCE
  27 void
  28 bid64dq_mul (UINT64 * pres, UINT64 * px, UINT128 * py
  29              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  30              _EXC_INFO_PARAM) {
  31   UINT64 x = *px;
  32 #if !DECIMAL_GLOBAL_ROUNDING
  33   unsigned int rnd_mode = *prnd_mode;
  34 #endif
  35 #else
  36 UINT64
  37 bid64dq_mul (UINT64 x, UINT128 y
  38              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  39              _EXC_INFO_PARAM) {
  40 #endif
  41   UINT64 res = 0xbaddbaddbaddbaddull;
  42   UINT128 x1;
  43
  44 #if DECIMAL_CALL_BY_REFERENCE
  45   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  46   bid64qq_mul (&res, &x1, py
  47                _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
  48                _EXC_INFO_ARG);
  49 #else
  50   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  51   res = bid64qq_mul (x1, y
  52                      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
  53                      _EXC_INFO_ARG);
  54 #endif
  55   BID_RETURN (res);
  56 }
  57
  58
  59 #if DECIMAL_CALL_BY_REFERENCE
  60 void
  61 bid64qd_mul (UINT64 * pres, UINT128 * px, UINT64 * py
  62              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  63              _EXC_INFO_PARAM) {
  64   UINT64 y = *py;
  65 #if !DECIMAL_GLOBAL_ROUNDING
  66   unsigned int rnd_mode = *prnd_mode;
  67 #endif
  68 #else
  69 UINT64
  70 bid64qd_mul (UINT128 x, UINT64 y
  71              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  72              _EXC_INFO_PARAM) {
  73 #endif
  74   UINT64 res = 0xbaddbaddbaddbaddull;
  75   UINT128 y1;
  76
  77 #if DECIMAL_CALL_BY_REFERENCE
  78   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  79   bid64qq_mul (&res, px, &y1
  80                _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
  81                _EXC_INFO_ARG);
  82 #else
  83   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
  84   res = bid64qq_mul (x, y1
  85                      _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
  86                      _EXC_INFO_ARG);
  87 #endif
  88   BID_RETURN (res);
  89 }
  90
  91
  92 #if DECIMAL_CALL_BY_REFERENCE
  93 void
  94 bid64qq_mul (UINT64 * pres, UINT128 * px, UINT128 * py
  95              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  96              _EXC_INFO_PARAM) {
  97   UINT128 x = *px, y = *py;
  98 #if !DECIMAL_GLOBAL_ROUNDING
  99   unsigned int rnd_mode = *prnd_mode;
 100 #endif
 101 #else
 102 UINT64
 103 bid64qq_mul (UINT128 x, UINT128 y
 104              _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 105              _EXC_INFO_PARAM) {
 106 #endif
 107
 108   UINT128 z = { {0x0000000000000000ull, 0x5ffe000000000000ull}
 109   };
 110   UINT64 res = 0xbaddbaddbaddbaddull;
 111   UINT64 x_sign, y_sign, p_sign;
 112   UINT64 x_exp, y_exp, p_exp;
 113   int true_p_exp;
 114   UINT128 C1, C2;
 115
 116   BID_SWAP128 (z);
 117   // skip cases where at least one operand is NaN or infinity
 118   if (!(((x.w[HIGH_128W] & MASK_NAN) == MASK_NAN) ||
 119         ((y.w[HIGH_128W] & MASK_NAN) == MASK_NAN) ||
 120         ((x.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF) ||
 121         ((y.w[HIGH_128W] & MASK_ANY_INF) == MASK_INF))) {
 122     // x, y are 0 or f but not inf or NaN => unpack the arguments and check
 123     // for non-canonical values
 124
 125     x_sign = x.w[HIGH_128W] & MASK_SIGN;        // 0 for positive, MASK_SIGN for negative
 126     C1.w[1] = x.w[HIGH_128W] & MASK_COEFF;
 127     C1.w[0] = x.w[LOW_128W];
 128     // check for non-canonical values - treated as zero
 129     if ((x.w[HIGH_128W] & 0x6000000000000000ull) ==
 130         0x6000000000000000ull) {
 131       // G0_G1=11 => non-canonical
 132       x_exp = (x.w[HIGH_128W] << 2) & MASK_EXP; // biased and shifted left 49 bits
 133       C1.w[1] = 0;      // significand high
 134       C1.w[0] = 0;      // significand low
 135     } else {    // G0_G1 != 11
 136       x_exp = x.w[HIGH_128W] & MASK_EXP;        // biased and shifted left 49 bits
 137       if (C1.w[1] > 0x0001ed09bead87c0ull ||
 138           (C1.w[1] == 0x0001ed09bead87c0ull &&
 139            C1.w[0] > 0x378d8e63ffffffffull)) {
 140         // x is non-canonical if coefficient is larger than 10^34 -1
 141         C1.w[1] = 0;
 142         C1.w[0] = 0;
 143       } else {  // canonical
 144         ;
 145       }
 146     }
 147     y_sign = y.w[HIGH_128W] & MASK_SIGN;        // 0 for positive, MASK_SIGN for negative
 148     C2.w[1] = y.w[HIGH_128W] & MASK_COEFF;
 149     C2.w[0] = y.w[LOW_128W];
 150     // check for non-canonical values - treated as zero
 151     if ((y.w[HIGH_128W] & 0x6000000000000000ull) ==
 152         0x6000000000000000ull) {
 153       // G0_G1=11 => non-canonical
 154       y_exp = (y.w[HIGH_128W] << 2) & MASK_EXP; // biased and shifted left 49 bits
 155       C2.w[1] = 0;      // significand high
 156       C2.w[0] = 0;      // significand low
 157     } else {    // G0_G1 != 11
 158       y_exp = y.w[HIGH_128W] & MASK_EXP;        // biased and shifted left 49 bits
 159       if (C2.w[1] > 0x0001ed09bead87c0ull ||
 160           (C2.w[1] == 0x0001ed09bead87c0ull &&
 161            C2.w[0] > 0x378d8e63ffffffffull)) {
 162         // y is non-canonical if coefficient is larger than 10^34 -1
 163         C2.w[1] = 0;
 164         C2.w[0] = 0;
 165       } else {  // canonical
 166         ;
 167       }
 168     }
 169     p_sign = x_sign ^ y_sign;   // sign of the product
 170
 171     true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
 172     // true_p_exp, p_exp are used only for 0 * 0, 0 * f, or f * 0
 173     if (true_p_exp < -398)
 174       p_exp = 0;        // cannot be less than EXP_MIN
 175     else if (true_p_exp > 369)
 176       p_exp = (UINT64) (369 + 398) << 53;       // cannot be more than EXP_MAX
 177     else
 178       p_exp = (UINT64) (true_p_exp + 398) << 53;
 179
 180     if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
 181         (C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
 182       // x = 0 or y = 0
 183       // the result is 0
 184       res = p_sign | p_exp;     // preferred exponent in [EXP_MIN, EXP_MAX]
 185       BID_RETURN (res)
 186     }   // else continue
 187   }
 188   // swap x and y - ensure that a NaN in x has 'higher precedence' than one in y
 189 #if DECIMAL_CALL_BY_REFERENCE
 190   bid64qqq_fma (&res, &y, &x, &z
 191                 _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 192                 _EXC_INFO_ARG);
 193 #else
 194   res = bid64qqq_fma (y, x, z
 195                       _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 196                       _EXC_INFO_ARG);
 197 #endif
 198   BID_RETURN (res);
 199 }
 200
 201
 202 #if DECIMAL_CALL_BY_REFERENCE
 203 void
 204 bid128dd_mul (UINT128 * pres, UINT64 * px, UINT64 * py
 205               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 206               _EXC_INFO_PARAM) {
 207   UINT64 x = *px, y = *py;
 208 #if !DECIMAL_GLOBAL_ROUNDING
 209   unsigned int rnd_mode = *prnd_mode;
 210 #endif
 211 #else
 212 UINT128
 213 bid128dd_mul (UINT64 x, UINT64 y
 214               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 215               _EXC_INFO_PARAM) {
 216 #endif
 217   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
 218   };
 219   UINT128 x1, y1;
 220
 221 #if DECIMAL_CALL_BY_REFERENCE
 222   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 223   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 224   bid128_mul (&res, &x1, &y1
 225               _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 226               _EXC_INFO_ARG);
 227 #else
 228   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 229   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 230   res = bid128_mul (x1, y1
 231                     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 232                     _EXC_INFO_ARG);
 233 #endif
 234   BID_RETURN (res);
 235 }
 236
 237
 238 #if DECIMAL_CALL_BY_REFERENCE
 239 void
 240 bid128dq_mul (UINT128 * pres, UINT64 * px, UINT128 * py
 241               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 242               _EXC_INFO_PARAM) {
 243   UINT64 x = *px;
 244 #if !DECIMAL_GLOBAL_ROUNDING
 245   unsigned int rnd_mode = *prnd_mode;
 246 #endif
 247 #else
 248 UINT128
 249 bid128dq_mul (UINT64 x, UINT128 y
 250               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 251               _EXC_INFO_PARAM) {
 252 #endif
 253   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
 254   };
 255   UINT128 x1;
 256
 257 #if DECIMAL_CALL_BY_REFERENCE
 258   bid64_to_bid128 (&x1, &x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 259   bid128_mul (&res, &x1, py
 260               _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 261               _EXC_INFO_ARG);
 262 #else
 263   x1 = bid64_to_bid128 (x _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 264   res = bid128_mul (x1, y
 265                     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 266                     _EXC_INFO_ARG);
 267 #endif
 268   BID_RETURN (res);
 269 }
 270
 271
 272 #if DECIMAL_CALL_BY_REFERENCE
 273 void
 274 bid128qd_mul (UINT128 * pres, UINT128 * px, UINT64 * py
 275               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 276               _EXC_INFO_PARAM) {
 277   UINT64 y = *py;
 278 #if !DECIMAL_GLOBAL_ROUNDING
 279   unsigned int rnd_mode = *prnd_mode;
 280 #endif
 281 #else
 282 UINT128
 283 bid128qd_mul (UINT128 x, UINT64 y
 284               _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 285               _EXC_INFO_PARAM) {
 286 #endif
 287   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
 288   };
 289   UINT128 y1;
 290
 291 #if DECIMAL_CALL_BY_REFERENCE
 292   bid64_to_bid128 (&y1, &y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 293   bid128_mul (&res, px, &y1
 294               _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 295               _EXC_INFO_ARG);
 296 #else
 297   y1 = bid64_to_bid128 (y _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 298   res = bid128_mul (x, y1
 299                     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 300                     _EXC_INFO_ARG);
 301 #endif
 302   BID_RETURN (res);
 303 }
 304
 305
 306 // bid128_mul stands for bid128qq_mul
 307 #if DECIMAL_CALL_BY_REFERENCE
 308 void
 309 bid128_mul (UINT128 * pres, UINT128 * px,
 310             UINT128 *
 311             py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 312             _EXC_INFO_PARAM) {
 313   UINT128 x = *px, y = *py;
 314
 315 #if !DECIMAL_GLOBAL_ROUNDING
 316   unsigned int rnd_mode = *prnd_mode;
 317
 318 #endif
 319 #else
 320 UINT128
 321 bid128_mul (UINT128 x,
 322             UINT128 y _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 323             _EXC_INFO_PARAM) {
 324
 325 #endif
 326   UINT128 z = { {0x0000000000000000ull, 0x5ffe000000000000ull}
 327   };
 328   UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
 329   };
 330   UINT64 x_sign, y_sign, p_sign;
 331   UINT64 x_exp, y_exp, p_exp;
 332   int true_p_exp;
 333   UINT128 C1, C2;
 334
 335   BID_SWAP128 (x);
 336   BID_SWAP128 (y);
 337   // skip cases where at least one operand is NaN or infinity
 338   if (!(((x.w[1] & MASK_NAN) == MASK_NAN) ||
 339         ((y.w[1] & MASK_NAN) == MASK_NAN) ||
 340         ((x.w[1] & MASK_ANY_INF) == MASK_INF) ||
 341         ((y.w[1] & MASK_ANY_INF) == MASK_INF))) {
 342     // x, y are 0 or f but not inf or NaN => unpack the arguments and check
 343     // for non-canonical values
 344
 345     x_sign = x.w[1] & MASK_SIGN;        // 0 for positive, MASK_SIGN for negative
 346     C1.w[1] = x.w[1] & MASK_COEFF;
 347     C1.w[0] = x.w[0];
 348     // check for non-canonical values - treated as zero
 349     if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
 350       // G0_G1=11 => non-canonical
 351       x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
 352       C1.w[1] = 0;      // significand high
 353       C1.w[0] = 0;      // significand low
 354     } else {    // G0_G1 != 11
 355       x_exp = x.w[1] & MASK_EXP;        // biased and shifted left 49 bits
 356       if (C1.w[1] > 0x0001ed09bead87c0ull ||
 357           (C1.w[1] == 0x0001ed09bead87c0ull &&
 358            C1.w[0] > 0x378d8e63ffffffffull)) {
 359         // x is non-canonical if coefficient is larger than 10^34 -1
 360         C1.w[1] = 0;
 361         C1.w[0] = 0;
 362       } else {  // canonical
 363         ;
 364       }
 365     }
 366     y_sign = y.w[1] & MASK_SIGN;        // 0 for positive, MASK_SIGN for negative
 367     C2.w[1] = y.w[1] & MASK_COEFF;
 368     C2.w[0] = y.w[0];
 369     // check for non-canonical values - treated as zero
 370     if ((y.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {
 371       // G0_G1=11 => non-canonical
 372       y_exp = (y.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
 373       C2.w[1] = 0;      // significand high
 374       C2.w[0] = 0;      // significand low
 375     } else {    // G0_G1 != 11
 376       y_exp = y.w[1] & MASK_EXP;        // biased and shifted left 49 bits
 377       if (C2.w[1] > 0x0001ed09bead87c0ull ||
 378           (C2.w[1] == 0x0001ed09bead87c0ull &&
 379            C2.w[0] > 0x378d8e63ffffffffull)) {
 380         // y is non-canonical if coefficient is larger than 10^34 -1
 381         C2.w[1] = 0;
 382         C2.w[0] = 0;
 383       } else {  // canonical
 384         ;
 385       }
 386     }
 387     p_sign = x_sign ^ y_sign;   // sign of the product
 388
 389     true_p_exp = (x_exp >> 49) - 6176 + (y_exp >> 49) - 6176;
 390     // true_p_exp, p_exp are used only for 0 * 0, 0 * f, or f * 0
 391     if (true_p_exp < -6176)
 392       p_exp = 0;        // cannot be less than EXP_MIN
 393     else if (true_p_exp > 6111)
 394       p_exp = (UINT64) (6111 + 6176) << 49;     // cannot be more than EXP_MAX
 395     else
 396       p_exp = (UINT64) (true_p_exp + 6176) << 49;
 397
 398     if ((C1.w[1] == 0x0 && C1.w[0] == 0x0) ||
 399         (C2.w[1] == 0x0 && C2.w[0] == 0x0)) {
 400       // x = 0 or y = 0
 401       // the result is 0
 402       res.w[1] = p_sign | p_exp;        // preferred exponent in [EXP_MIN, EXP_MAX]
 403       res.w[0] = 0x0;
 404       BID_SWAP128 (res);
 405       BID_RETURN (res)
 406     }   // else continue
 407   }
 408
 409   BID_SWAP128 (x);
 410   BID_SWAP128 (y);
 411   BID_SWAP128 (z);
 412   // swap x and y - ensure that a NaN in x has 'higher precedence' than one in y
 413 #if DECIMAL_CALL_BY_REFERENCE
 414   bid128_fma (&res, &y, &x, &z
 415               _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 416               _EXC_INFO_ARG);
 417 #else
 418   res = bid128_fma (y, x, z
 419                     _RND_MODE_ARG _EXC_FLAGS_ARG _EXC_MASKS_ARG
 420                     _EXC_INFO_ARG);
 421 #endif
 422   BID_RETURN (res);
 423 }