02110-1301, USA. */
#define BID_128RES
-#include "div_macros.h"
+#include "bid_div_macros.h"
-BID128_FUNCTION_ARG2(__bid128_rem, x, y)
+BID128_FUNCTION_ARG2_NORND_CUSTOMRESTYPE (UINT128, bid128_rem, x, y)
- UINT256 P256;
- UINT128 CX, CY, CX2, CQ, CR, T, CXS, P128, res;
- UINT64 sign_x, sign_y;
- SINT64 D;
- int_float f64, fx;
- int exponent_x = 0, exponent_y, diff_expon, bin_expon_cx, scale, scale0;
+ UINT256 P256;
+ UINT128 CX, CY, CX2, CQ, CR, T, CXS, P128, res;
+ UINT64 sign_x, sign_y, valid_y;
+ SINT64 D;
+ int_float f64, fx;
+ int exponent_x, exponent_y, diff_expon, bin_expon_cx, scale,
+ scale0;
// unpack arguments, check for NaN or Infinity
- if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
+
+valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
+
+if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
#ifdef SET_STATUS_FLAGS
- if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
- __set_status_flags (pfpsf, INVALID_EXCEPTION);
+if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
// test if x is NaN
- if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
+if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
#ifdef SET_STATUS_FLAGS
- if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
- __set_status_flags (pfpsf, INVALID_EXCEPTION);
+ if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
- res.w[1] = x.w[1] & QUIET_MASK64;
- res.w[0] = x.w[0];
- BID_RETURN (res);
- }
+ res.w[1] = CX.w[1] & QUIET_MASK64;
+ res.w[0] = CX.w[0];
+ BID_RETURN (res);
+}
// x is Infinity?
- if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
- // check if y is Inf.
- if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
- // return NaN
- {
+if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
+ // check if y is Inf.
+ if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
+ // return NaN
+ {
#ifdef SET_STATUS_FLAGS
- // set status flags
- __set_status_flags (pfpsf, INVALID_EXCEPTION);
+ // set status flags
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
- res.w[1] = 0x7c00000000000000ull;
- res.w[0] = 0;
- BID_RETURN (res);
- }
+ res.w[1] = 0x7c00000000000000ull;
+ res.w[0] = 0;
+ BID_RETURN (res);
+ }
- }
+}
// x is 0
- if ((y.w[1] & INFINITY_MASK64) != INFINITY_MASK64) { // not Inf/NaN
- if ((!y.w[0]) && !(y.w[1] & 0x0001ffffffffffffull)) {
+if ((!CY.w[1]) && (!CY.w[0])) {
#ifdef SET_STATUS_FLAGS
- // set status flags
- __set_status_flags (pfpsf, INVALID_EXCEPTION);
+ // set status flags
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
- // x=y=0, return NaN
- res.w[1] = 0x7c00000000000000ull;
- res.w[0] = 0;
- BID_RETURN (res);
- }
- // return 0
- exponent_y = (y.w[1] >> 49) & 0x3fff;
- if (exponent_x > exponent_y)
- exponent_x = exponent_y;
-
- res.w[1] = sign_x | (((UINT64) exponent_x) << 49);
- res.w[0] = 0;
- BID_RETURN (res);
- }
- }
- if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
- // y is Inf. or NaN
+ // x=y=0, return NaN
+ res.w[1] = 0x7c00000000000000ull;
+ res.w[0] = 0;
+ BID_RETURN (res);
+}
+if (valid_y || ((y.w[1] & NAN_MASK64) == INFINITY_MASK64)) {
+ // return 0
+ if ((exponent_x > exponent_y)
+ && ((y.w[1] & NAN_MASK64) != INFINITY_MASK64))
+ exponent_x = exponent_y;
+
+ res.w[1] = sign_x | (((UINT64) exponent_x) << 49);
+ res.w[0] = 0;
+ BID_RETURN (res);
+}
+}
+if (!valid_y) {
+ // y is Inf. or NaN
- // test if y is NaN
- if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
+ // test if y is NaN
+ if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
#ifdef SET_STATUS_FLAGS
- if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
- __set_status_flags (pfpsf, INVALID_EXCEPTION);
+ if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64) // y is sNaN
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
- res.w[1] = y.w[1] & QUIET_MASK64;
- res.w[0] = y.w[0];
- BID_RETURN (res);
- }
- // y is Infinity?
- if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
- // return x
- res.w[1] = x.w[1];
- res.w[0] = x.w[0];
- BID_RETURN (res);
- }
- // y is 0
+ res.w[1] = CY.w[1] & QUIET_MASK64;
+ res.w[0] = CY.w[0];
+ BID_RETURN (res);
+ }
+ // y is Infinity?
+ if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
+ // return x
+ res.w[1] = x.w[1];
+ res.w[0] = x.w[0];
+ BID_RETURN (res);
+ }
+ // y is 0
#ifdef SET_STATUS_FLAGS
- // set status flags
- __set_status_flags (pfpsf, INVALID_EXCEPTION);
+ // set status flags
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
- res.w[1] = 0x7c00000000000000ull;
- res.w[0] = 0;
+ res.w[1] = 0x7c00000000000000ull;
+ res.w[0] = 0;
+ BID_RETURN (res);
+}
+
+diff_expon = exponent_x - exponent_y;
+
+if (diff_expon <= 0) {
+ diff_expon = -diff_expon;
+
+ if (diff_expon > 34) {
+ // |x|<|y| in this case
+ res = x;
BID_RETURN (res);
}
+ // set exponent of y to exponent_x, scale coefficient_y
+ T = power10_table_128[diff_expon];
+ __mul_128x128_to_256 (P256, CY, T);
- diff_expon = exponent_x - exponent_y;
-
- if (diff_expon <= 0) {
- diff_expon = -diff_expon;
-
- if (diff_expon > 34) {
- // |x|<|y| in this case
- res = x;
- BID_RETURN (res);
- }
- // set exponent of y to exponent_x, scale coefficient_y
- T = __bid_power10_table_128[diff_expon];
- __mul_128x128_to_256 (P256, CY, T);
-
- if (P256.w[2] || P256.w[3]) {
- // |x|<|y| in this case
- res = x;
- BID_RETURN (res);
- }
-
- CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
- CX2.w[0] = CX.w[0] << 1;
- if (__unsigned_compare_ge_128 (P256, CX2)) {
- // |x|<|y| in this case
- res = x;
- BID_RETURN (res);
- }
-
- P128.w[0] = P256.w[0];
- P128.w[1] = P256.w[1];
- __div_128_by_128 (&CQ, &CR, CX, P128);
-
- CX2.w[1] = (CR.w[1] << 1) | (CR.w[0] >> 63);
- CX2.w[0] = CR.w[0] << 1;
- if ((__unsigned_compare_gt_128 (CX2, P256))
- || (CX2.w[1] == P256.w[1] && CX2.w[0] == P256.w[0]
- && (CQ.w[0] & 1))) {
- __sub_128_128 (CR, P256, CR);
- sign_x ^= 0x8000000000000000ull;
- }
-
- get_BID128_very_fast (&res, sign_x, exponent_x, CR);
+ if (P256.w[2] || P256.w[3]) {
+ // |x|<|y| in this case
+ res = x;
BID_RETURN (res);
}
- // 2^64
- f64.i = 0x5f800000;
-
- scale0 = 38;
- if (!CY.w[1])
- scale0 = 34;
-
- while (diff_expon > 0) {
- // get number of digits in CX and scale=38-digits
- // fx ~ CX
- fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
- bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
- scale = scale0 - __bid_estimate_decimal_digits[bin_expon_cx];
- // scale = 38-__bid_estimate_decimal_digits[bin_expon_cx];
- D = CX.w[1] - __bid_power10_index_binexp_128[bin_expon_cx].w[1];
- if (D > 0
- || (!D
- && CX.w[0] >= __bid_power10_index_binexp_128[bin_expon_cx].w[0]))
- scale--;
-
- if (diff_expon >= scale)
- diff_expon -= scale;
- else {
- scale = diff_expon;
- diff_expon = 0;
- }
-
- T = __bid_power10_table_128[scale];
- __mul_128x128_low (CXS, CX, T);
-
- __div_128_by_128 (&CQ, &CX, CXS, CY);
-
- // check for remainder == 0
- if (!CX.w[1] && !CX.w[0]) {
- get_BID128_very_fast (&res, sign_x, exponent_y, CX);
- BID_RETURN (res);
- }
- }
CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
CX2.w[0] = CX.w[0] << 1;
- if ((__unsigned_compare_gt_128 (CX2, CY))
- || (CX2.w[1] == CY.w[1] && CX2.w[0] == CY.w[0] && (CQ.w[0] & 1))) {
- __sub_128_128 (CX, CY, CX);
+ if (__unsigned_compare_ge_128 (P256, CX2)) {
+ // |x|<|y| in this case
+ res = x;
+ BID_RETURN (res);
+ }
+
+ P128.w[0] = P256.w[0];
+ P128.w[1] = P256.w[1];
+ __div_128_by_128 (&CQ, &CR, CX, P128);
+
+ CX2.w[1] = (CR.w[1] << 1) | (CR.w[0] >> 63);
+ CX2.w[0] = CR.w[0] << 1;
+ if ((__unsigned_compare_gt_128 (CX2, P256))
+ || (CX2.w[1] == P256.w[1] && CX2.w[0] == P256.w[0]
+ && (CQ.w[0] & 1))) {
+ __sub_128_128 (CR, P256, CR);
sign_x ^= 0x8000000000000000ull;
}
- get_BID128_very_fast (&res, sign_x, exponent_y, CX);
+ get_BID128_very_fast (&res, sign_x, exponent_x, CR);
BID_RETURN (res);
}
+ // 2^64
+f64.i = 0x5f800000;
+
+scale0 = 38;
+if (!CY.w[1])
+ scale0 = 34;
+
+while (diff_expon > 0) {
+ // get number of digits in CX and scale=38-digits
+ // fx ~ CX
+ fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
+ bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
+ scale = scale0 - estimate_decimal_digits[bin_expon_cx];
+ // scale = 38-estimate_decimal_digits[bin_expon_cx];
+ D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];
+ if (D > 0
+ || (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))
+ scale--;
+
+ if (diff_expon >= scale)
+ diff_expon -= scale;
+ else {
+ scale = diff_expon;
+ diff_expon = 0;
+ }
+
+ T = power10_table_128[scale];
+ __mul_128x128_low (CXS, CX, T);
+
+ __div_128_by_128 (&CQ, &CX, CXS, CY);
+
+ // check for remainder == 0
+ if (!CX.w[1] && !CX.w[0]) {
+ get_BID128_very_fast (&res, sign_x, exponent_y, CX);
+ BID_RETURN (res);
+ }
+}
+
+CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
+CX2.w[0] = CX.w[0] << 1;
+if ((__unsigned_compare_gt_128 (CX2, CY))
+ || (CX2.w[1] == CY.w[1] && CX2.w[0] == CY.w[0] && (CQ.w[0] & 1))) {
+ __sub_128_128 (CX, CY, CX);
+ sign_x ^= 0x8000000000000000ull;
+}
+
+get_BID128_very_fast (&res, sign_x, exponent_y, CX);
+BID_RETURN (res);
+}
OSDN Git Service
