aboutsummaryrefslogtreecommitdiffstats
path: root/gcc-4.4.3/libgcc/config/libbid/bid64_sqrt.c
diff options
context:
space:
mode:
Diffstat (limited to 'gcc-4.4.3/libgcc/config/libbid/bid64_sqrt.c')
-rw-r--r--gcc-4.4.3/libgcc/config/libbid/bid64_sqrt.c552
1 files changed, 552 insertions, 0 deletions
diff --git a/gcc-4.4.3/libgcc/config/libbid/bid64_sqrt.c b/gcc-4.4.3/libgcc/config/libbid/bid64_sqrt.c
new file mode 100644
index 000000000..310ab9459
--- /dev/null
+++ b/gcc-4.4.3/libgcc/config/libbid/bid64_sqrt.c
@@ -0,0 +1,552 @@
+/* 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 square root
+ *****************************************************************************
+ *
+ * Algorithm description:
+ *
+ * if(exponent_x is odd)
+ * scale coefficient_x by 10, adjust exponent
+ * - get lower estimate for number of digits in coefficient_x
+ * - scale coefficient x to between 31 and 33 decimal digits
+ * - in parallel, check for exact case and return if true
+ * - get high part of result coefficient using double precision sqrt
+ * - compute remainder and refine coefficient in one iteration (which
+ * modifies it by at most 1)
+ * - result exponent is easy to compute from the adjusted arg. exponent
+ *
+ ****************************************************************************/
+
+#include "bid_internal.h"
+#include "bid_sqrt_macros.h"
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+#include <fenv.h>
+
+#define FE_ALL_FLAGS FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW|FE_UNDERFLOW|FE_INEXACT
+#endif
+
+extern double sqrt (double);
+
+#if DECIMAL_CALL_BY_REFERENCE
+
+void
+bid64_sqrt (UINT64 * pres,
+ UINT64 *
+ px _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
+ _EXC_INFO_PARAM) {
+ UINT64 x;
+#else
+
+UINT64
+bid64_sqrt (UINT64 x _RND_MODE_PARAM _EXC_FLAGS_PARAM
+ _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
+#endif
+ UINT128 CA, CT;
+ UINT64 sign_x, coefficient_x;
+ UINT64 Q, Q2, A10, C4, R, R2, QE, res;
+ SINT64 D;
+ int_double t_scale;
+ int_float tempx;
+ double da, dq, da_h, da_l, dqe;
+ int exponent_x, exponent_q, bin_expon_cx;
+ int digits_x;
+ int scale;
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+ fexcept_t binaryflags = 0;
+#endif
+
+#if DECIMAL_CALL_BY_REFERENCE
+#if !DECIMAL_GLOBAL_ROUNDING
+ _IDEC_round rnd_mode = *prnd_mode;
+#endif
+ x = *px;
+#endif
+
+ // unpack arguments, check for NaN or Infinity
+ if (!unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x)) {
+ // x is Inf. or NaN or 0
+ if ((x & INFINITY_MASK64) == INFINITY_MASK64) {
+ res = coefficient_x;
+ if ((coefficient_x & SSNAN_MASK64) == SINFINITY_MASK64) // -Infinity
+ {
+ res = NAN_MASK64;
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
+#endif
+ }
+#ifdef SET_STATUS_FLAGS
+ if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
+#endif
+ BID_RETURN (res & QUIET_MASK64);
+ }
+ // x is 0
+ exponent_x = (exponent_x + DECIMAL_EXPONENT_BIAS) >> 1;
+ res = sign_x | (((UINT64) exponent_x) << 53);
+ BID_RETURN (res);
+ }
+ // x<0?
+ if (sign_x && coefficient_x) {
+ res = NAN_MASK64;
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
+#endif
+ BID_RETURN (res);
+ }
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+ (void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
+#endif
+ //--- get number of bits in the coefficient of x ---
+ tempx.d = (float) coefficient_x;
+ bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
+ digits_x = estimate_decimal_digits[bin_expon_cx];
+ // add test for range
+ if (coefficient_x >= power10_index_binexp[bin_expon_cx])
+ digits_x++;
+
+ A10 = coefficient_x;
+ if (exponent_x & 1) {
+ A10 = (A10 << 2) + A10;
+ A10 += A10;
+ }
+
+ dqe = sqrt ((double) A10);
+ QE = (UINT32) dqe;
+ if (QE * QE == A10) {
+ res =
+ very_fast_get_BID64 (0, (exponent_x + DECIMAL_EXPONENT_BIAS) >> 1,
+ QE);
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+ (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
+#endif
+ BID_RETURN (res);
+ }
+ // if exponent is odd, scale coefficient by 10
+ scale = 31 - digits_x;
+ exponent_q = exponent_x - scale;
+ scale += (exponent_q & 1); // exp. bias is even
+
+ CT = power10_table_128[scale];
+ __mul_64x128_short (CA, coefficient_x, CT);
+
+ // 2^64
+ t_scale.i = 0x43f0000000000000ull;
+ // convert CA to DP
+ da_h = CA.w[1];
+ da_l = CA.w[0];
+ da = da_h * t_scale.d + da_l;
+
+ dq = sqrt (da);
+
+ Q = (UINT64) dq;
+
+ // get sign(sqrt(CA)-Q)
+ R = CA.w[0] - Q * Q;
+ R = ((SINT64) R) >> 63;
+ D = R + R + 1;
+
+ exponent_q = (exponent_q + DECIMAL_EXPONENT_BIAS) >> 1;
+
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INEXACT_EXCEPTION);
+#endif
+
+#ifndef IEEE_ROUND_NEAREST
+#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
+ if (!((rnd_mode) & 3)) {
+#endif
+#endif
+
+ // midpoint to check
+ Q2 = Q + Q + D;
+ C4 = CA.w[0] << 2;
+
+ // get sign(-sqrt(CA)+Midpoint)
+ R2 = Q2 * Q2 - C4;
+ R2 = ((SINT64) R2) >> 63;
+
+ // adjust Q if R!=R2
+ Q += (D & (R ^ R2));
+#ifndef IEEE_ROUND_NEAREST
+#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
+ } else {
+ C4 = CA.w[0];
+ Q += D;
+ if ((SINT64) (Q * Q - C4) > 0)
+ Q--;
+ if (rnd_mode == ROUNDING_UP)
+ Q++;
+ }
+#endif
+#endif
+
+ res = fast_get_BID64 (0, exponent_q, Q);
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+ (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
+#endif
+ BID_RETURN (res);
+}
+
+
+TYPE0_FUNCTION_ARG1 (UINT64, bid64q_sqrt, x)
+
+ UINT256 M256, C4, C8;
+ UINT128 CX, CX2, A10, S2, T128, CS, CSM, CS2, C256, CS1,
+ mul_factor2_long = { {0x0ull, 0x0ull} }, QH, Tmp, TP128, Qh, Ql;
+UINT64 sign_x, Carry, B10, res, mul_factor, mul_factor2 = 0x0ull, CS0;
+SINT64 D;
+int_float fx, f64;
+int exponent_x, bin_expon_cx, done = 0;
+int digits, scale, exponent_q = 0, exact = 1, amount, extra_digits;
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+fexcept_t binaryflags = 0;
+#endif
+
+ // unpack arguments, check for NaN or Infinity
+if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
+ res = CX.w[1];
+ // NaN ?
+ if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
+#ifdef SET_STATUS_FLAGS
+ if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
+#endif
+ Tmp.w[1] = (CX.w[1] & 0x00003fffffffffffull);
+ Tmp.w[0] = CX.w[0];
+ TP128 = reciprocals10_128[18];
+ __mul_128x128_full (Qh, Ql, Tmp, TP128);
+ amount = recip_scale[18];
+ __shr_128 (Tmp, Qh, amount);
+ res = (CX.w[1] & 0xfc00000000000000ull) | Tmp.w[0];
+ BID_RETURN (res);
+ }
+ // x is Infinity?
+ if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
+ if (sign_x) {
+ // -Inf, return NaN
+ res = 0x7c00000000000000ull;
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
+#endif
+ }
+ BID_RETURN (res);
+ }
+ // x is 0 otherwise
+
+ exponent_x =
+ ((exponent_x - DECIMAL_EXPONENT_BIAS_128) >> 1) +
+ DECIMAL_EXPONENT_BIAS;
+ if (exponent_x < 0)
+ exponent_x = 0;
+ if (exponent_x > DECIMAL_MAX_EXPON_64)
+ exponent_x = DECIMAL_MAX_EXPON_64;
+ //res= sign_x | (((UINT64)exponent_x)<<53);
+ res = get_BID64 (sign_x, exponent_x, 0, rnd_mode, pfpsf);
+ BID_RETURN (res);
+}
+if (sign_x) {
+ res = 0x7c00000000000000ull;
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INVALID_EXCEPTION);
+#endif
+ BID_RETURN (res);
+}
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);
+#endif
+
+ // 2^64
+f64.i = 0x5f800000;
+
+ // fx ~ CX
+fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
+bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
+digits = estimate_decimal_digits[bin_expon_cx];
+
+A10 = CX;
+if (exponent_x & 1) {
+ A10.w[1] = (CX.w[1] << 3) | (CX.w[0] >> 61);
+ A10.w[0] = CX.w[0] << 3;
+ CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
+ CX2.w[0] = CX.w[0] << 1;
+ __add_128_128 (A10, A10, CX2);
+}
+
+C256.w[1] = A10.w[1];
+C256.w[0] = A10.w[0];
+CS.w[0] = short_sqrt128 (A10);
+CS.w[1] = 0;
+mul_factor = 0;
+ // check for exact result
+if (CS.w[0] < 10000000000000000ull) {
+ if (CS.w[0] * CS.w[0] == A10.w[0]) {
+ __sqr64_fast (S2, CS.w[0]);
+ if (S2.w[1] == A10.w[1]) // && S2.w[0]==A10.w[0])
+ {
+ res =
+ get_BID64 (0,
+ ((exponent_x - DECIMAL_EXPONENT_BIAS_128) >> 1) +
+ DECIMAL_EXPONENT_BIAS, CS.w[0], rnd_mode, pfpsf);
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+ (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
+#endif
+ BID_RETURN (res);
+ }
+ }
+ if (CS.w[0] >= 1000000000000000ull) {
+ done = 1;
+ exponent_q = exponent_x;
+ C256.w[1] = A10.w[1];
+ C256.w[0] = A10.w[0];
+ }
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INEXACT_EXCEPTION);
+#endif
+ exact = 0;
+} else {
+ B10 = 0x3333333333333334ull;
+ __mul_64x64_to_128_full (CS2, CS.w[0], B10);
+ CS0 = CS2.w[1] >> 1;
+ if (CS.w[0] != ((CS0 << 3) + (CS0 << 1))) {
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INEXACT_EXCEPTION);
+#endif
+ exact = 0;
+ }
+ done = 1;
+ CS.w[0] = CS0;
+ exponent_q = exponent_x + 2;
+ mul_factor = 10;
+ mul_factor2 = 100;
+ if (CS.w[0] >= 10000000000000000ull) {
+ __mul_64x64_to_128_full (CS2, CS.w[0], B10);
+ CS0 = CS2.w[1] >> 1;
+ if (CS.w[0] != ((CS0 << 3) + (CS0 << 1))) {
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INEXACT_EXCEPTION);
+#endif
+ exact = 0;
+ }
+ exponent_q += 2;
+ CS.w[0] = CS0;
+ mul_factor = 100;
+ mul_factor2 = 10000;
+ }
+ if (exact) {
+ CS0 = CS.w[0] * mul_factor;
+ __sqr64_fast (CS1, CS0)
+ if ((CS1.w[0] != A10.w[0]) || (CS1.w[1] != A10.w[1])) {
+#ifdef SET_STATUS_FLAGS
+ __set_status_flags (pfpsf, INEXACT_EXCEPTION);
+#endif
+ exact = 0;
+ }
+ }
+}
+
+if (!done) {
+ // get number of digits in 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]))
+ digits++;
+
+ // if exponent is odd, scale coefficient by 10
+ scale = 31 - digits;
+ exponent_q = exponent_x - scale;
+ scale += (exponent_q & 1); // exp. bias is even
+
+ T128 = power10_table_128[scale];
+ __mul_128x128_low (C256, CX, T128);
+
+
+ CS.w[0] = short_sqrt128 (C256);
+}
+ //printf("CS=%016I64x\n",CS.w[0]);
+
+exponent_q =
+ ((exponent_q - DECIMAL_EXPONENT_BIAS_128) >> 1) +
+ DECIMAL_EXPONENT_BIAS;
+if ((exponent_q < 0) && (exponent_q + MAX_FORMAT_DIGITS >= 0)) {
+ extra_digits = -exponent_q;
+ exponent_q = 0;
+
+ // get coeff*(2^M[extra_digits])/10^extra_digits
+ __mul_64x64_to_128 (QH, CS.w[0], reciprocals10_64[extra_digits]);
+
+ // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
+ amount = short_recip_scale[extra_digits];
+
+ CS0 = QH.w[1] >> amount;
+
+#ifdef SET_STATUS_FLAGS
+ if (exact) {
+ if (CS.w[0] != CS0 * power10_table_128[extra_digits].w[0])
+ exact = 0;
+ }
+ if (!exact)
+ __set_status_flags (pfpsf, UNDERFLOW_EXCEPTION | INEXACT_EXCEPTION);
+#endif
+
+ CS.w[0] = CS0;
+ if (!mul_factor)
+ mul_factor = 1;
+ mul_factor *= power10_table_128[extra_digits].w[0];
+ __mul_64x64_to_128 (mul_factor2_long, mul_factor, mul_factor);
+ if (mul_factor2_long.w[1])
+ mul_factor2 = 0;
+ else
+ mul_factor2 = mul_factor2_long.w[1];
+}
+ // 4*C256
+C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62);
+C4.w[0] = C256.w[0] << 2;
+
+#ifndef IEEE_ROUND_NEAREST
+#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
+if (!((rnd_mode) & 3)) {
+#endif
+#endif
+ // compare to midpoints
+ CSM.w[0] = (CS.w[0] + CS.w[0]) | 1;
+ //printf("C256=%016I64x %016I64x, CSM=%016I64x %016I64x %016I64x\n",C4.w[1],C4.w[0],CSM.w[1],CSM.w[0], CS.w[0]);
+ if (mul_factor)
+ CSM.w[0] *= mul_factor;
+ // CSM^2
+ __mul_64x64_to_128 (M256, CSM.w[0], CSM.w[0]);
+ //__mul_128x128_to_256(M256, CSM, CSM);
+
+ if (C4.w[1] > M256.w[1] ||
+ (C4.w[1] == M256.w[1] && C4.w[0] > M256.w[0])) {
+ // round up
+ CS.w[0]++;
+ } else {
+ C8.w[0] = CS.w[0] << 3;
+ C8.w[1] = 0;
+ if (mul_factor) {
+ if (mul_factor2) {
+ __mul_64x64_to_128 (C8, C8.w[0], mul_factor2);
+ } else {
+ __mul_64x128_low (C8, C8.w[0], mul_factor2_long);
+ }
+ }
+ // M256 - 8*CSM
+ __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
+ M256.w[1] = M256.w[1] - C8.w[1] - Carry;
+
+ // if CSM' > C256, round up
+ if (M256.w[1] > C4.w[1] ||
+ (M256.w[1] == C4.w[1] && M256.w[0] > C4.w[0])) {
+ // round down
+ if (CS.w[0])
+ CS.w[0]--;
+ }
+ }
+#ifndef IEEE_ROUND_NEAREST
+#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
+} else {
+ CS.w[0]++;
+ CSM.w[0] = CS.w[0];
+ C8.w[0] = CSM.w[0] << 1;
+ if (mul_factor)
+ CSM.w[0] *= mul_factor;
+ __mul_64x64_to_128 (M256, CSM.w[0], CSM.w[0]);
+ C8.w[1] = 0;
+ if (mul_factor) {
+ if (mul_factor2) {
+ __mul_64x64_to_128 (C8, C8.w[0], mul_factor2);
+ } else {
+ __mul_64x128_low (C8, C8.w[0], mul_factor2_long);
+ }
+ }
+ //printf("C256=%016I64x %016I64x, CSM=%016I64x %016I64x %016I64x\n",C256.w[1],C256.w[0],M256.w[1],M256.w[0], CS.w[0]);
+
+ if (M256.w[1] > C256.w[1] ||
+ (M256.w[1] == C256.w[1] && M256.w[0] > C256.w[0])) {
+ __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
+ M256.w[1] = M256.w[1] - Carry - C8.w[1];
+ M256.w[0]++;
+ if (!M256.w[0]) {
+ M256.w[1]++;
+
+ }
+
+ if ((M256.w[1] > C256.w[1] ||
+ (M256.w[1] == C256.w[1] && M256.w[0] > C256.w[0]))
+ && (CS.w[0] > 1)) {
+
+ CS.w[0]--;
+
+ if (CS.w[0] > 1) {
+ __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
+ M256.w[1] = M256.w[1] - Carry - C8.w[1];
+ M256.w[0]++;
+ if (!M256.w[0]) {
+ M256.w[1]++;
+ }
+
+ if (M256.w[1] > C256.w[1] ||
+ (M256.w[1] == C256.w[1] && M256.w[0] > C256.w[0]))
+ CS.w[0]--;
+ }
+ }
+ }
+
+ else {
+ /*__add_carry_out(M256.w[0], Carry, M256.w[0], C8.w[0]);
+ M256.w[1] = M256.w[1] + Carry + C8.w[1];
+ M256.w[0]++;
+ if(!M256.w[0])
+ {
+ M256.w[1]++;
+ }
+ CS.w[0]++;
+ if(M256.w[1]<C256.w[1] ||
+ (M256.w[1]==C256.w[1] && M256.w[0]<=C256.w[0]))
+ {
+ CS.w[0]++;
+ }*/
+ CS.w[0]++;
+ }
+ //printf("C256=%016I64x %016I64x, CSM=%016I64x %016I64x %016I64x %d\n",C4.w[1],C4.w[0],M256.w[1],M256.w[0], CS.w[0], exact);
+ // RU?
+ if (((rnd_mode) != ROUNDING_UP) || exact) {
+ if (CS.w[0])
+ CS.w[0]--;
+ }
+
+}
+#endif
+#endif
+ //printf("C256=%016I64x %016I64x, CSM=%016I64x %016I64x %016I64x %d\n",C4.w[1],C4.w[0],M256.w[1],M256.w[0], CS.w[0], exact);
+
+res = get_BID64 (0, exponent_q, CS.w[0], rnd_mode, pfpsf);
+#ifdef UNCHANGED_BINARY_STATUS_FLAGS
+(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
+#endif
+BID_RETURN (res);
+
+
+}