URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gcc-4.5.1/] [libgcc/] [config/] [libbid/] [bid128_sqrt.c] - Blame information for rev 826

Details | Compare with Previous | View Log


/* 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/>.  */
 
#define BID_128RES
#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
 
BID128_FUNCTION_ARG1 (bid128_sqrt, x)
 
     UINT256 M256, C256, C4, C8;
     UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res;
     UINT64 sign_x, Carry;
     SINT64 D;
     int_float fx, f64;
     int exponent_x, bin_expon_cx;
     int digits, scale, exponent_q;
#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.w[1] = CX.w[1];
res.w[0] = CX.w[0];
    // 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
  res.w[1] = CX.w[1] & QUIET_MASK64;
  BID_RETURN (res);
}
    // x is Infinity?
if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
  res.w[1] = CX.w[1];
  if (sign_x) {
    // -Inf, return NaN
    res.w[1] = 0x7c00000000000000ull;
#ifdef SET_STATUS_FLAGS
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
  }
  BID_RETURN (res);
}
    // x is 0 otherwise
 
res.w[1] =
  sign_x |
  ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1) << 49);
res.w[0] = 0;
BID_RETURN (res);
}
if (sign_x) {
  res.w[1] = 0x7c00000000000000ull;
  res.w[0] = 0;
#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);
}
 
CS.w[0] = short_sqrt128 (A10);
CS.w[1] = 0;
  // check for exact result
if (CS.w[0] * CS.w[0] == A10.w[0]) {
  __mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]);
  if (S2.w[1] == A10.w[1])      // && S2.w[0]==A10.w[0])
  {
    get_BID128_very_fast (&res, 0,
                          (exponent_x +
                           DECIMAL_EXPONENT_BIAS_128) >> 1, CS);
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
#endif
    BID_RETURN (res);
  }
}
  // 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 = 67 - digits;
exponent_q = exponent_x - scale;
scale += (exponent_q & 1);      // exp. bias is even
 
if (scale > 38) {
  T128 = power10_table_128[scale - 37];
  __mul_128x128_low (CX1, CX, T128);
 
  TP128 = power10_table_128[37];
  __mul_128x128_to_256 (C256, CX1, TP128);
} else {
  T128 = power10_table_128[scale];
  __mul_128x128_to_256 (C256, CX, T128);
}
 
 
  // 4*C256
C4.w[3] = (C256.w[3] << 2) | (C256.w[2] >> 62);
C4.w[2] = (C256.w[2] << 2) | (C256.w[1] >> 62);
C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62);
C4.w[0] = C256.w[0] << 2;
 
long_sqrt128 (&CS, C256);
 
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
if (!((rnd_mode) & 3)) {
#endif
#endif
  // compare to midpoints
  CSM.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
  CSM.w[0] = (CS.w[0] + CS.w[0]) | 1;
  // CSM^2
  //__mul_128x128_to_256(M256, CSM, CSM);
  __sqr128_to_256 (M256, CSM);
 
  if (C4.w[3] > M256.w[3]
      || (C4.w[3] == M256.w[3]
          && (C4.w[2] > M256.w[2]
              || (C4.w[2] == M256.w[2]
                  && (C4.w[1] > M256.w[1]
                      || (C4.w[1] == M256.w[1]
                          && C4.w[0] > M256.w[0])))))) {
    // round up
    CS.w[0]++;
    if (!CS.w[0])
      CS.w[1]++;
  } else {
    C8.w[1] = (CS.w[1] << 3) | (CS.w[0] >> 61);
    C8.w[0] = CS.w[0] << 3;
    // M256 - 8*CSM
    __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
    __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
    __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
    M256.w[3] = M256.w[3] - Carry;
 
    // if CSM' > C256, round up
    if (M256.w[3] > C4.w[3]
        || (M256.w[3] == C4.w[3]
            && (M256.w[2] > C4.w[2]
                || (M256.w[2] == C4.w[2]
                    && (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[1]--;
      CS.w[0]--;
    }
  }
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
} else {
  __sqr128_to_256 (M256, CS);
  C8.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
  C8.w[0] = CS.w[0] << 1;
  if (M256.w[3] > C256.w[3]
      || (M256.w[3] == C256.w[3]
          && (M256.w[2] > C256.w[2]
              || (M256.w[2] == C256.w[2]
                  && (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]);
    __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
    __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
    M256.w[3] = M256.w[3] - Carry;
    M256.w[0]++;
    if (!M256.w[0]) {
      M256.w[1]++;
      if (!M256.w[1]) {
        M256.w[2]++;
        if (!M256.w[2])
          M256.w[3]++;
      }
    }
 
    if (!CS.w[0])
      CS.w[1]--;
    CS.w[0]--;
 
    if (M256.w[3] > C256.w[3]
        || (M256.w[3] == C256.w[3]
            && (M256.w[2] > C256.w[2]
                || (M256.w[2] == C256.w[2]
                    && (M256.w[1] > C256.w[1]
                        || (M256.w[1] == C256.w[1]
                            && M256.w[0] > C256.w[0])))))) {
 
      if (!CS.w[0])
        CS.w[1]--;
      CS.w[0]--;
    }
  }
 
  else {
    __add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
    __add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
    __add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
    M256.w[3] = M256.w[3] + Carry;
    M256.w[0]++;
    if (!M256.w[0]) {
      M256.w[1]++;
      if (!M256.w[1]) {
        M256.w[2]++;
        if (!M256.w[2])
          M256.w[3]++;
      }
    }
    if (M256.w[3] < C256.w[3]
        || (M256.w[3] == C256.w[3]
            && (M256.w[2] < C256.w[2]
                || (M256.w[2] == C256.w[2]
                    && (M256.w[1] < C256.w[1]
                        || (M256.w[1] == C256.w[1]
                            && M256.w[0] <= C256.w[0])))))) {
 
      CS.w[0]++;
      if (!CS.w[0])
        CS.w[1]++;
    }
  }
  // RU?
  if ((rnd_mode) == ROUNDING_UP) {
    CS.w[0]++;
    if (!CS.w[0])
      CS.w[1]++;
  }
 
}
#endif
#endif
 
#ifdef SET_STATUS_FLAGS
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
#endif
get_BID128_fast (&res, 0,
                 (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1, CS);
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
#endif
BID_RETURN (res);
}
 
 
 
BID128_FUNCTION_ARGTYPE1 (bid128d_sqrt, UINT64, x)
 
     UINT256 M256, C256, C4, C8;
     UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res;
     UINT64 sign_x, Carry;
     SINT64 D;
     int_float fx, f64;
     int exponent_x, bin_expon_cx;
     int digits, scale, exponent_q;
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
     fexcept_t binaryflags = 0;
#endif
 
        // unpack arguments, check for NaN or Infinity
   // unpack arguments, check for NaN or Infinity
CX.w[1] = 0;
if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {
res.w[1] = CX.w[0];
res.w[0] = 0;
           // NaN ?
if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
#ifdef SET_STATUS_FLAGS
  if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN
    __set_status_flags (pfpsf, INVALID_EXCEPTION);
#endif
  res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);
  __mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);
  res.w[1] |= ((CX.w[0]) & 0xfc00000000000000ull);
  BID_RETURN (res);
}
           // x is Infinity?
if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
  if (sign_x) {
    // -Inf, return NaN
    res.w[1] = 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 + DECIMAL_EXPONENT_BIAS_128;
res.w[1] =
  sign_x | ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1)
            << 49);
res.w[0] = 0;
BID_RETURN (res);
}
if (sign_x) {
  res.w[1] = 0x7c00000000000000ull;
  res.w[0] = 0;
#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
exponent_x =
  exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128;
 
           // 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);
}
 
CS.w[0] = short_sqrt128 (A10);
CS.w[1] = 0;
           // check for exact result
if (CS.w[0] * CS.w[0] == A10.w[0]) {
  __mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]);
  if (S2.w[1] == A10.w[1]) {
    get_BID128_very_fast (&res, 0,
                          (exponent_x + DECIMAL_EXPONENT_BIAS_128) >> 1,
                          CS);
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
    (void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
#endif
    BID_RETURN (res);
  }
}
           // 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 = 67 - digits;
exponent_q = exponent_x - scale;
scale += (exponent_q & 1);      // exp. bias is even
 
if (scale > 38) {
  T128 = power10_table_128[scale - 37];
  __mul_128x128_low (CX1, CX, T128);
 
  TP128 = power10_table_128[37];
  __mul_128x128_to_256 (C256, CX1, TP128);
} else {
  T128 = power10_table_128[scale];
  __mul_128x128_to_256 (C256, CX, T128);
}
 
 
           // 4*C256
C4.w[3] = (C256.w[3] << 2) | (C256.w[2] >> 62);
C4.w[2] = (C256.w[2] << 2) | (C256.w[1] >> 62);
C4.w[1] = (C256.w[1] << 2) | (C256.w[0] >> 62);
C4.w[0] = C256.w[0] << 2;
 
long_sqrt128 (&CS, C256);
 
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
if (!((rnd_mode) & 3)) {
#endif
#endif
  // compare to midpoints
  CSM.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
  CSM.w[0] = (CS.w[0] + CS.w[0]) | 1;
  // CSM^2
  //__mul_128x128_to_256(M256, CSM, CSM);
  __sqr128_to_256 (M256, CSM);
 
  if (C4.w[3] > M256.w[3]
      || (C4.w[3] == M256.w[3]
          && (C4.w[2] > M256.w[2]
              || (C4.w[2] == M256.w[2]
                  && (C4.w[1] > M256.w[1]
                      || (C4.w[1] == M256.w[1]
                          && C4.w[0] > M256.w[0])))))) {
    // round up
    CS.w[0]++;
    if (!CS.w[0])
      CS.w[1]++;
  } else {
    C8.w[1] = (CS.w[1] << 3) | (CS.w[0] >> 61);
    C8.w[0] = CS.w[0] << 3;
    // M256 - 8*CSM
    __sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
    __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
    __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
    M256.w[3] = M256.w[3] - Carry;
 
    // if CSM' > C256, round up
    if (M256.w[3] > C4.w[3]
        || (M256.w[3] == C4.w[3]
            && (M256.w[2] > C4.w[2]
                || (M256.w[2] == C4.w[2]
                    && (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[1]--;
      CS.w[0]--;
    }
  }
#ifndef IEEE_ROUND_NEAREST
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
} else {
  __sqr128_to_256 (M256, CS);
  C8.w[1] = (CS.w[1] << 1) | (CS.w[0] >> 63);
  C8.w[0] = CS.w[0] << 1;
  if (M256.w[3] > C256.w[3]
      || (M256.w[3] == C256.w[3]
          && (M256.w[2] > C256.w[2]
              || (M256.w[2] == C256.w[2]
                  && (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]);
    __sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
    __sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
    M256.w[3] = M256.w[3] - Carry;
    M256.w[0]++;
    if (!M256.w[0]) {
      M256.w[1]++;
      if (!M256.w[1]) {
        M256.w[2]++;
        if (!M256.w[2])
          M256.w[3]++;
      }
    }
 
    if (!CS.w[0])
      CS.w[1]--;
    CS.w[0]--;
 
    if (M256.w[3] > C256.w[3]
        || (M256.w[3] == C256.w[3]
            && (M256.w[2] > C256.w[2]
                || (M256.w[2] == C256.w[2]
                    && (M256.w[1] > C256.w[1]
                        || (M256.w[1] == C256.w[1]
                            && M256.w[0] > C256.w[0])))))) {
 
      if (!CS.w[0])
        CS.w[1]--;
      CS.w[0]--;
    }
  }
 
  else {
    __add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]);
    __add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);
    __add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);
    M256.w[3] = M256.w[3] + Carry;
    M256.w[0]++;
    if (!M256.w[0]) {
      M256.w[1]++;
      if (!M256.w[1]) {
        M256.w[2]++;
        if (!M256.w[2])
          M256.w[3]++;
      }
    }
    if (M256.w[3] < C256.w[3]
        || (M256.w[3] == C256.w[3]
            && (M256.w[2] < C256.w[2]
                || (M256.w[2] == C256.w[2]
                    && (M256.w[1] < C256.w[1]
                        || (M256.w[1] == C256.w[1]
                            && M256.w[0] <= C256.w[0])))))) {
 
      CS.w[0]++;
      if (!CS.w[0])
        CS.w[1]++;
    }
  }
  // RU?
  if ((rnd_mode) == ROUNDING_UP) {
    CS.w[0]++;
    if (!CS.w[0])
      CS.w[1]++;
  }
 
}
#endif
#endif
 
#ifdef SET_STATUS_FLAGS
__set_status_flags (pfpsf, INEXACT_EXCEPTION);
#endif
get_BID128_fast (&res, 0, (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1,
                 CS);
#ifdef UNCHANGED_BINARY_STATUS_FLAGS
(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);
#endif
BID_RETURN (res);
 
 
}

Browse

Tools

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-stable/] [gcc-4.5.1/] [libgcc/] [config/] [libbid/] [bid128_sqrt.c] - Blame information for rev 826

Line No.	Rev	Author	Line
1	272	jeremybenn	`/* Copyright (C) 2007, 2009 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			`#define BID_128RES`
25			`#include "bid_internal.h"`
26			`#include "bid_sqrt_macros.h"`
27			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
28			`#include <fenv.h>`
29
30			`#define FE_ALL_FLAGS FE_INVALID\|FE_DIVBYZERO\|FE_OVERFLOW\|FE_UNDERFLOW\|FE_INEXACT`
31			`#endif`
32
33			`BID128_FUNCTION_ARG1 (bid128_sqrt, x)`
34
35			`UINT256 M256, C256, C4, C8;`
36			`UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res;`
37			`UINT64 sign_x, Carry;`
38			`SINT64 D;`
39			`int_float fx, f64;`
40			`int exponent_x, bin_expon_cx;`
41			`int digits, scale, exponent_q;`
42			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
43			`fexcept_t binaryflags = 0;`
44			`#endif`
45
46			`// unpack arguments, check for NaN or Infinity`
47			`if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {`
48			`res.w[1] = CX.w[1];`
49			`res.w[0] = CX.w[0];`
50			`// NaN ?`
51			`if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {`
52			`#ifdef SET_STATUS_FLAGS`
53			`if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN`
54			`__set_status_flags (pfpsf, INVALID_EXCEPTION);`
55			`#endif`
56			`res.w[1] = CX.w[1] & QUIET_MASK64;`
57			`BID_RETURN (res);`
58			`}`
59			`// x is Infinity?`
60			`if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {`
61			`res.w[1] = CX.w[1];`
62			`if (sign_x) {`
63			`// -Inf, return NaN`
64			`res.w[1] = 0x7c00000000000000ull;`
65			`#ifdef SET_STATUS_FLAGS`
66			`__set_status_flags (pfpsf, INVALID_EXCEPTION);`
67			`#endif`
68			`}`
69			`BID_RETURN (res);`
70			`}`
71			`// x is 0 otherwise`
72
73			`res.w[1] =`
74			`sign_x \|`
75			`((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1) << 49);`
76			`res.w[0] = 0;`
77			`BID_RETURN (res);`
78			`}`
79			`if (sign_x) {`
80			`res.w[1] = 0x7c00000000000000ull;`
81			`res.w[0] = 0;`
82			`#ifdef SET_STATUS_FLAGS`
83			`__set_status_flags (pfpsf, INVALID_EXCEPTION);`
84			`#endif`
85			`BID_RETURN (res);`
86			`}`
87			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
88			`(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);`
89			`#endif`
90			`// 2^64`
91			`f64.i = 0x5f800000;`
92
93			`// fx ~ CX`
94			`fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];`
95			`bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;`
96			`digits = estimate_decimal_digits[bin_expon_cx];`
97
98			`A10 = CX;`
99			`if (exponent_x & 1) {`
100			`A10.w[1] = (CX.w[1] << 3) \| (CX.w[0] >> 61);`
101			`A10.w[0] = CX.w[0] << 3;`
102			`CX2.w[1] = (CX.w[1] << 1) \| (CX.w[0] >> 63);`
103			`CX2.w[0] = CX.w[0] << 1;`
104			`__add_128_128 (A10, A10, CX2);`
105			`}`
106
107			`CS.w[0] = short_sqrt128 (A10);`
108			`CS.w[1] = 0;`
109			`// check for exact result`
110			`if (CS.w[0] * CS.w[0] == A10.w[0]) {`
111			`__mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]);`
112			`if (S2.w[1] == A10.w[1]) // && S2.w[0]==A10.w[0])`
113			`{`
114			`get_BID128_very_fast (&res, 0,`
115			`(exponent_x +`
116			`DECIMAL_EXPONENT_BIAS_128) >> 1, CS);`
117			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
118			`(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);`
119			`#endif`
120			`BID_RETURN (res);`
121			`}`
122			`}`
123			`// get number of digits in CX`
124			`D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];`
125			`if (D > 0`
126			`\|\| (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))`
127			`digits++;`
128
129			`// if exponent is odd, scale coefficient by 10`
130			`scale = 67 - digits;`
131			`exponent_q = exponent_x - scale;`
132			`scale += (exponent_q & 1); // exp. bias is even`
133
134			`if (scale > 38) {`
135			`T128 = power10_table_128[scale - 37];`
136			`__mul_128x128_low (CX1, CX, T128);`
137
138			`TP128 = power10_table_128[37];`
139			`__mul_128x128_to_256 (C256, CX1, TP128);`
140			`} else {`
141			`T128 = power10_table_128[scale];`
142			`__mul_128x128_to_256 (C256, CX, T128);`
143			`}`
144
145
146			`// 4*C256`
147			`C4.w[3] = (C256.w[3] << 2) \| (C256.w[2] >> 62);`
148			`C4.w[2] = (C256.w[2] << 2) \| (C256.w[1] >> 62);`
149			`C4.w[1] = (C256.w[1] << 2) \| (C256.w[0] >> 62);`
150			`C4.w[0] = C256.w[0] << 2;`
151
152			`long_sqrt128 (&CS, C256);`
153
154			`#ifndef IEEE_ROUND_NEAREST`
155			`#ifndef IEEE_ROUND_NEAREST_TIES_AWAY`
156			`if (!((rnd_mode) & 3)) {`
157			`#endif`
158			`#endif`
159			`// compare to midpoints`
160			`CSM.w[1] = (CS.w[1] << 1) \| (CS.w[0] >> 63);`
161			`CSM.w[0] = (CS.w[0] + CS.w[0]) \| 1;`
162			`// CSM^2`
163			`//__mul_128x128_to_256(M256, CSM, CSM);`
164			`__sqr128_to_256 (M256, CSM);`
165
166			`if (C4.w[3] > M256.w[3]`
167			`\|\| (C4.w[3] == M256.w[3]`
168			`&& (C4.w[2] > M256.w[2]`
169			`\|\| (C4.w[2] == M256.w[2]`
170			`&& (C4.w[1] > M256.w[1]`
171			`\|\| (C4.w[1] == M256.w[1]`
172			`&& C4.w[0] > M256.w[0])))))) {`
173			`// round up`
174			`CS.w[0]++;`
175			`if (!CS.w[0])`
176			`CS.w[1]++;`
177			`} else {`
178			`C8.w[1] = (CS.w[1] << 3) \| (CS.w[0] >> 61);`
179			`C8.w[0] = CS.w[0] << 3;`
180			`// M256 - 8*CSM`
181			`__sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);`
182			`__sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);`
183			`__sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);`
184			`M256.w[3] = M256.w[3] - Carry;`
185
186			`// if CSM' > C256, round up`
187			`if (M256.w[3] > C4.w[3]`
188			`\|\| (M256.w[3] == C4.w[3]`
189			`&& (M256.w[2] > C4.w[2]`
190			`\|\| (M256.w[2] == C4.w[2]`
191			`&& (M256.w[1] > C4.w[1]`
192			`\|\| (M256.w[1] == C4.w[1]`
193			`&& M256.w[0] > C4.w[0])))))) {`
194			`// round down`
195			`if (!CS.w[0])`
196			`CS.w[1]--;`
197			`CS.w[0]--;`
198			`}`
199			`}`
200			`#ifndef IEEE_ROUND_NEAREST`
201			`#ifndef IEEE_ROUND_NEAREST_TIES_AWAY`
202			`} else {`
203			`__sqr128_to_256 (M256, CS);`
204			`C8.w[1] = (CS.w[1] << 1) \| (CS.w[0] >> 63);`
205			`C8.w[0] = CS.w[0] << 1;`
206			`if (M256.w[3] > C256.w[3]`
207			`\|\| (M256.w[3] == C256.w[3]`
208			`&& (M256.w[2] > C256.w[2]`
209			`\|\| (M256.w[2] == C256.w[2]`
210			`&& (M256.w[1] > C256.w[1]`
211			`\|\| (M256.w[1] == C256.w[1]`
212			`&& M256.w[0] > C256.w[0])))))) {`
213			`__sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);`
214			`__sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);`
215			`__sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);`
216			`M256.w[3] = M256.w[3] - Carry;`
217			`M256.w[0]++;`
218			`if (!M256.w[0]) {`
219			`M256.w[1]++;`
220			`if (!M256.w[1]) {`
221			`M256.w[2]++;`
222			`if (!M256.w[2])`
223			`M256.w[3]++;`
224			`}`
225			`}`
226
227			`if (!CS.w[0])`
228			`CS.w[1]--;`
229			`CS.w[0]--;`
230
231			`if (M256.w[3] > C256.w[3]`
232			`\|\| (M256.w[3] == C256.w[3]`
233			`&& (M256.w[2] > C256.w[2]`
234			`\|\| (M256.w[2] == C256.w[2]`
235			`&& (M256.w[1] > C256.w[1]`
236			`\|\| (M256.w[1] == C256.w[1]`
237			`&& M256.w[0] > C256.w[0])))))) {`
238
239			`if (!CS.w[0])`
240			`CS.w[1]--;`
241			`CS.w[0]--;`
242			`}`
243			`}`
244
245			`else {`
246			`__add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]);`
247			`__add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);`
248			`__add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);`
249			`M256.w[3] = M256.w[3] + Carry;`
250			`M256.w[0]++;`
251			`if (!M256.w[0]) {`
252			`M256.w[1]++;`
253			`if (!M256.w[1]) {`
254			`M256.w[2]++;`
255			`if (!M256.w[2])`
256			`M256.w[3]++;`
257			`}`
258			`}`
259			`if (M256.w[3] < C256.w[3]`
260			`\|\| (M256.w[3] == C256.w[3]`
261			`&& (M256.w[2] < C256.w[2]`
262			`\|\| (M256.w[2] == C256.w[2]`
263			`&& (M256.w[1] < C256.w[1]`
264			`\|\| (M256.w[1] == C256.w[1]`
265			`&& M256.w[0] <= C256.w[0])))))) {`
266
267			`CS.w[0]++;`
268			`if (!CS.w[0])`
269			`CS.w[1]++;`
270			`}`
271			`}`
272			`// RU?`
273			`if ((rnd_mode) == ROUNDING_UP) {`
274			`CS.w[0]++;`
275			`if (!CS.w[0])`
276			`CS.w[1]++;`
277			`}`
278
279			`}`
280			`#endif`
281			`#endif`
282
283			`#ifdef SET_STATUS_FLAGS`
284			`__set_status_flags (pfpsf, INEXACT_EXCEPTION);`
285			`#endif`
286			`get_BID128_fast (&res, 0,`
287			`(exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1, CS);`
288			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
289			`(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);`
290			`#endif`
291			`BID_RETURN (res);`
292			`}`
293
294
295
296			`BID128_FUNCTION_ARGTYPE1 (bid128d_sqrt, UINT64, x)`
297
298			`UINT256 M256, C256, C4, C8;`
299			`UINT128 CX, CX1, CX2, A10, S2, T128, TP128, CS, CSM, res;`
300			`UINT64 sign_x, Carry;`
301			`SINT64 D;`
302			`int_float fx, f64;`
303			`int exponent_x, bin_expon_cx;`
304			`int digits, scale, exponent_q;`
305			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
306			`fexcept_t binaryflags = 0;`
307			`#endif`
308
309			`// unpack arguments, check for NaN or Infinity`
310			`// unpack arguments, check for NaN or Infinity`
311			`CX.w[1] = 0;`
312			`if (!unpack_BID64 (&sign_x, &exponent_x, &CX.w[0], x)) {`
313			`res.w[1] = CX.w[0];`
314			`res.w[0] = 0;`
315			`// NaN ?`
316			`if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) {`
317			`#ifdef SET_STATUS_FLAGS`
318			`if ((x & SNAN_MASK64) == SNAN_MASK64) // sNaN`
319			`__set_status_flags (pfpsf, INVALID_EXCEPTION);`
320			`#endif`
321			`res.w[0] = (CX.w[0] & 0x0003ffffffffffffull);`
322			`__mul_64x64_to_128 (res, res.w[0], power10_table_128[18].w[0]);`
323			`res.w[1] \|= ((CX.w[0]) & 0xfc00000000000000ull);`
324			`BID_RETURN (res);`
325			`}`
326			`// x is Infinity?`
327			`if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {`
328			`if (sign_x) {`
329			`// -Inf, return NaN`
330			`res.w[1] = 0x7c00000000000000ull;`
331			`#ifdef SET_STATUS_FLAGS`
332			`__set_status_flags (pfpsf, INVALID_EXCEPTION);`
333			`#endif`
334			`}`
335			`BID_RETURN (res);`
336			`}`
337			`// x is 0 otherwise`
338
339			`exponent_x =`
340			`exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128;`
341			`res.w[1] =`
342			`sign_x \| ((((UINT64) (exponent_x + DECIMAL_EXPONENT_BIAS_128)) >> 1)`
343			`<< 49);`
344			`res.w[0] = 0;`
345			`BID_RETURN (res);`
346			`}`
347			`if (sign_x) {`
348			`res.w[1] = 0x7c00000000000000ull;`
349			`res.w[0] = 0;`
350			`#ifdef SET_STATUS_FLAGS`
351			`__set_status_flags (pfpsf, INVALID_EXCEPTION);`
352			`#endif`
353			`BID_RETURN (res);`
354			`}`
355			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
356			`(void) fegetexceptflag (&binaryflags, FE_ALL_FLAGS);`
357			`#endif`
358			`exponent_x =`
359			`exponent_x - DECIMAL_EXPONENT_BIAS + DECIMAL_EXPONENT_BIAS_128;`
360
361			`// 2^64`
362			`f64.i = 0x5f800000;`
363
364			`// fx ~ CX`
365			`fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];`
366			`bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;`
367			`digits = estimate_decimal_digits[bin_expon_cx];`
368
369			`A10 = CX;`
370			`if (exponent_x & 1) {`
371			`A10.w[1] = (CX.w[1] << 3) \| (CX.w[0] >> 61);`
372			`A10.w[0] = CX.w[0] << 3;`
373			`CX2.w[1] = (CX.w[1] << 1) \| (CX.w[0] >> 63);`
374			`CX2.w[0] = CX.w[0] << 1;`
375			`__add_128_128 (A10, A10, CX2);`
376			`}`
377
378			`CS.w[0] = short_sqrt128 (A10);`
379			`CS.w[1] = 0;`
380			`// check for exact result`
381			`if (CS.w[0] * CS.w[0] == A10.w[0]) {`
382			`__mul_64x64_to_128_fast (S2, CS.w[0], CS.w[0]);`
383			`if (S2.w[1] == A10.w[1]) {`
384			`get_BID128_very_fast (&res, 0,`
385			`(exponent_x + DECIMAL_EXPONENT_BIAS_128) >> 1,`
386			`CS);`
387			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
388			`(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);`
389			`#endif`
390			`BID_RETURN (res);`
391			`}`
392			`}`
393			`// get number of digits in CX`
394			`D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];`
395			`if (D > 0`
396			`\|\| (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))`
397			`digits++;`
398
399			`// if exponent is odd, scale coefficient by 10`
400			`scale = 67 - digits;`
401			`exponent_q = exponent_x - scale;`
402			`scale += (exponent_q & 1); // exp. bias is even`
403
404			`if (scale > 38) {`
405			`T128 = power10_table_128[scale - 37];`
406			`__mul_128x128_low (CX1, CX, T128);`
407
408			`TP128 = power10_table_128[37];`
409			`__mul_128x128_to_256 (C256, CX1, TP128);`
410			`} else {`
411			`T128 = power10_table_128[scale];`
412			`__mul_128x128_to_256 (C256, CX, T128);`
413			`}`
414
415
416			`// 4*C256`
417			`C4.w[3] = (C256.w[3] << 2) \| (C256.w[2] >> 62);`
418			`C4.w[2] = (C256.w[2] << 2) \| (C256.w[1] >> 62);`
419			`C4.w[1] = (C256.w[1] << 2) \| (C256.w[0] >> 62);`
420			`C4.w[0] = C256.w[0] << 2;`
421
422			`long_sqrt128 (&CS, C256);`
423
424			`#ifndef IEEE_ROUND_NEAREST`
425			`#ifndef IEEE_ROUND_NEAREST_TIES_AWAY`
426			`if (!((rnd_mode) & 3)) {`
427			`#endif`
428			`#endif`
429			`// compare to midpoints`
430			`CSM.w[1] = (CS.w[1] << 1) \| (CS.w[0] >> 63);`
431			`CSM.w[0] = (CS.w[0] + CS.w[0]) \| 1;`
432			`// CSM^2`
433			`//__mul_128x128_to_256(M256, CSM, CSM);`
434			`__sqr128_to_256 (M256, CSM);`
435
436			`if (C4.w[3] > M256.w[3]`
437			`\|\| (C4.w[3] == M256.w[3]`
438			`&& (C4.w[2] > M256.w[2]`
439			`\|\| (C4.w[2] == M256.w[2]`
440			`&& (C4.w[1] > M256.w[1]`
441			`\|\| (C4.w[1] == M256.w[1]`
442			`&& C4.w[0] > M256.w[0])))))) {`
443			`// round up`
444			`CS.w[0]++;`
445			`if (!CS.w[0])`
446			`CS.w[1]++;`
447			`} else {`
448			`C8.w[1] = (CS.w[1] << 3) \| (CS.w[0] >> 61);`
449			`C8.w[0] = CS.w[0] << 3;`
450			`// M256 - 8*CSM`
451			`__sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);`
452			`__sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);`
453			`__sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);`
454			`M256.w[3] = M256.w[3] - Carry;`
455
456			`// if CSM' > C256, round up`
457			`if (M256.w[3] > C4.w[3]`
458			`\|\| (M256.w[3] == C4.w[3]`
459			`&& (M256.w[2] > C4.w[2]`
460			`\|\| (M256.w[2] == C4.w[2]`
461			`&& (M256.w[1] > C4.w[1]`
462			`\|\| (M256.w[1] == C4.w[1]`
463			`&& M256.w[0] > C4.w[0])))))) {`
464			`// round down`
465			`if (!CS.w[0])`
466			`CS.w[1]--;`
467			`CS.w[0]--;`
468			`}`
469			`}`
470			`#ifndef IEEE_ROUND_NEAREST`
471			`#ifndef IEEE_ROUND_NEAREST_TIES_AWAY`
472			`} else {`
473			`__sqr128_to_256 (M256, CS);`
474			`C8.w[1] = (CS.w[1] << 1) \| (CS.w[0] >> 63);`
475			`C8.w[0] = CS.w[0] << 1;`
476			`if (M256.w[3] > C256.w[3]`
477			`\|\| (M256.w[3] == C256.w[3]`
478			`&& (M256.w[2] > C256.w[2]`
479			`\|\| (M256.w[2] == C256.w[2]`
480			`&& (M256.w[1] > C256.w[1]`
481			`\|\| (M256.w[1] == C256.w[1]`
482			`&& M256.w[0] > C256.w[0])))))) {`
483			`__sub_borrow_out (M256.w[0], Carry, M256.w[0], C8.w[0]);`
484			`__sub_borrow_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);`
485			`__sub_borrow_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);`
486			`M256.w[3] = M256.w[3] - Carry;`
487			`M256.w[0]++;`
488			`if (!M256.w[0]) {`
489			`M256.w[1]++;`
490			`if (!M256.w[1]) {`
491			`M256.w[2]++;`
492			`if (!M256.w[2])`
493			`M256.w[3]++;`
494			`}`
495			`}`
496
497			`if (!CS.w[0])`
498			`CS.w[1]--;`
499			`CS.w[0]--;`
500
501			`if (M256.w[3] > C256.w[3]`
502			`\|\| (M256.w[3] == C256.w[3]`
503			`&& (M256.w[2] > C256.w[2]`
504			`\|\| (M256.w[2] == C256.w[2]`
505			`&& (M256.w[1] > C256.w[1]`
506			`\|\| (M256.w[1] == C256.w[1]`
507			`&& M256.w[0] > C256.w[0])))))) {`
508
509			`if (!CS.w[0])`
510			`CS.w[1]--;`
511			`CS.w[0]--;`
512			`}`
513			`}`
514
515			`else {`
516			`__add_carry_out (M256.w[0], Carry, M256.w[0], C8.w[0]);`
517			`__add_carry_in_out (M256.w[1], Carry, M256.w[1], C8.w[1], Carry);`
518			`__add_carry_in_out (M256.w[2], Carry, M256.w[2], 0, Carry);`
519			`M256.w[3] = M256.w[3] + Carry;`
520			`M256.w[0]++;`
521			`if (!M256.w[0]) {`
522			`M256.w[1]++;`
523			`if (!M256.w[1]) {`
524			`M256.w[2]++;`
525			`if (!M256.w[2])`
526			`M256.w[3]++;`
527			`}`
528			`}`
529			`if (M256.w[3] < C256.w[3]`
530			`\|\| (M256.w[3] == C256.w[3]`
531			`&& (M256.w[2] < C256.w[2]`
532			`\|\| (M256.w[2] == C256.w[2]`
533			`&& (M256.w[1] < C256.w[1]`
534			`\|\| (M256.w[1] == C256.w[1]`
535			`&& M256.w[0] <= C256.w[0])))))) {`
536
537			`CS.w[0]++;`
538			`if (!CS.w[0])`
539			`CS.w[1]++;`
540			`}`
541			`}`
542			`// RU?`
543			`if ((rnd_mode) == ROUNDING_UP) {`
544			`CS.w[0]++;`
545			`if (!CS.w[0])`
546			`CS.w[1]++;`
547			`}`
548
549			`}`
550			`#endif`
551			`#endif`
552
553			`#ifdef SET_STATUS_FLAGS`
554			`__set_status_flags (pfpsf, INEXACT_EXCEPTION);`
555			`#endif`
556			`get_BID128_fast (&res, 0, (exponent_q + DECIMAL_EXPONENT_BIAS_128) >> 1,`
557			`CS);`
558			`#ifdef UNCHANGED_BINARY_STATUS_FLAGS`
559			`(void) fesetexceptflag (&binaryflags, FE_ALL_FLAGS);`
560			`#endif`
561			`BID_RETURN (res);`
562
563
564			`}`