Subversion Repositories openrisc_me

/* 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);

}