/*******************************************************************************
+
+  LEDA 3.5
+
+  _bigfloat.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/

/*
 * This file has been automatically generated from "bigfloat.w" by CTANGLE
 * (Version 3.1 [km2c])
 */


#include <LEDA/bigfloat.h>
#include <LEDA/floatingpoint.h>


long bigfloat::global_prec = 53;
rounding_modes bigfloat::round_mode = TO_NEAREST;
bool bigfloat::dbool = true;
bigfloat::output_modes bigfloat::output_mode = bigfloat::DEC_OUT;

const bigfloat bigfloat::pZero(bigfloat::PZERO_VAL);
const bigfloat bigfloat::nZero(bigfloat::NZERO_VAL);
const bigfloat bigfloat::pInf(bigfloat::PINF_VAL);
const bigfloat bigfloat::nInf(bigfloat::NINF_VAL);
const bigfloat bigfloat::NaN(bigfloat::NAN_VAL);

const bigfloat bigfloat::zero = bigfloat::pZero;
const bigfloat bigfloat::one = bigfloat(1, 0);



const integer integer_zero = 0;
const integer integer_one = 1;
const double pow52 = power_two(52);



const long bin_maxlen = 10000;



long max(long l1, long l2)
{
  return (l1 > l2 ? l1 : l2);
}
double logtwo(double d)
{
  return log(d) / log(2.0);
}
long sign(long i)
{
  if (i == 0)
    return 0;
  else
    return (i > 0 ? 1 : -1);
}
bool is_small(const integer & x)
{
  if (!x.islong())
    return false;
  long x_long = x.tolong();
  if (x_long < 0)
    x_long = -x_long;
  return !(x_long & 0xf000000);
}


bool isnum(char ch)
{
  if ((ch >= '0') && (ch <= '9'))
    return true;
  else
    return false;
}




void bigfloat::normalize(void)
{
  if (special != bigfloat::NOT_VAL)
    return;

  if ((::sign(significant) == 0) && (special == bigfloat::NOT_VAL))
    special = bigfloat::PZERO_VAL;

  long z = significant.zeros();
  if (z > 0) {
    significant = significant >> z;
    exponent += z;
  }
  precision = significant.length();
}



void cut(integer & b, long prec)
{
  b = (b >> (b.length() - prec));
}


void binout(ostream & os, integer b)
{
  char temp[bin_maxlen];
  long count = 0;
  do {
    temp[count++] = (char) (b.tolong() & 1) + '0';
    b = b >> 1;
  }
  while (b.tolong());
  for (long i = count - 1; i >= 0; i--)
    os << temp[i];
}


bigfloat Powl(const bigfloat & x, long n, long prec = 1,
   rounding_modes mode = EXACT)
{
  bigfloat z = 1, y = x;
  long n_prefix = n;

  while (n_prefix > 0) {
    if (n_prefix % 2)
      z = mul(z, y, prec, mode);
    n_prefix = n_prefix / 2;
    y = mul(y, y, prec, mode);
  }
  return z;
}


void decimal_output(ostream & os, bigfloat b, long prec,
   rounding_modes mode = TO_NEAREST)
{
  if (!b.get_exponent().islong())
    error_handler(1, "decimal_output: not implemented for large exponents");
  if (prec == 0)
    error_handler(1, "decimal_output: prec has to be bigger than 0!");
  long dd = 10;

  if (sign(b) < 0) {
    b = -b;
    os << "-";
  }




  long log2_b = b.get_significant_length() + b.get_exponent().tolong();
  bigfloat b_rem(b.get_significant(), -b.get_significant_length());
  long log10_b = (long) ceil(log10(2.0) * log2_b + log10(to_double(b_rem)));




  long diff = prec - log10_b;
  long digits = (long) ceil(logtwo((double) dd) * diff) + log2_b;

  bigfloat b_shift;
  integer significant;
  if (diff >= 0)
    b_shift = mul(b, Powl(dd, diff, digits, mode), digits, mode);
  else
    b_shift = div(b, Powl(dd, -diff, digits, mode), digits, mode);

  significant = to_integer(b_shift, mode);
  string str = significant.tostring();




  int len = str.length();
  int end = len - 1;
  int tail = end;
  while (tail && str[tail] == '0')
    tail--;
  if (tail < end)
    str = str.del(tail + 1, end);
  if (tail > 0)
    os << str.insert(1, ".");
  else
    os << str;

  if (log10_b != 1)
    os << "E" << len - 1 - diff;
/* before exponent log10_b - 1 was not stable */
/* (numerical errors near negative exponents of 10) */


}


bigfloat outofchar(char *rep, long prec = 0)
{
  integer sig = 0, exp = 0;
  bigfloat result, pow;
  long dd = 10;
  long frac_length = 0;

  int s;

  s = 1;
  if (rep[0] == '-') {
    s = -1;
    rep++;
  }
  else if (rep[0] == '+')
    rep++;


  int sign = s;

  while (isnum(*rep)) {
    sig = sig * dd + (*(rep++) - '0');
  }



  if (*rep == '.') {
    rep++;
    while (isnum(*rep)) {
      sig = sig * dd + (*(rep++) - '0');
      frac_length++;
    }
  }




  if (*rep == 'E') {
    rep++;

    s = 1;
    if (rep[0] == '-') {
      s = -1;
      rep++;
    }
    else if (rep[0] == '+')
      rep++;


    while (isnum(*rep))
      exp = exp * dd + (*(rep++) - '0');
    if (s == -1)
      exp = -exp;
  }





  long l = sig.length() + 2;
  long exp_diff = exp.tolong() - frac_length;


  if (prec <= 0)
    prec = l;

  pow = Powl(dd, (exp_diff > 0 ? exp_diff : -exp_diff), 1, EXACT);

  if (exp_diff > 0)
    result = mul(sig, pow, prec, TO_NEAREST);
  else
    result = div(sig, pow, prec, TO_NEAREST);

  if (sign == 1)
    return result;
  else
    return -result;
}



/* member functions */


long bigfloat::set_precision(long p)
{
  long tmp = global_prec;
  global_prec = p;
  return tmp;
}
rounding_modes bigfloat::set_rounding_mode(rounding_modes m)
{
  rounding_modes tmp = round_mode;
  round_mode = m;
  return tmp;
}
bigfloat::output_modes bigfloat::set_output_mode(output_modes o_mode)
{
  output_modes tmp = output_mode;
  output_mode = o_mode;
  return tmp;
}





long bigfloat::get_significant_length(void) const
{
  return precision;
}
integer bigfloat::get_exponent(void) const
{
  return exponent;
}
integer bigfloat::get_significant(void) const
{
  return significant;
}



bool isNaN(const bigfloat & x)
{
  return (x.special == bigfloat::NAN_VAL);
}
bool isnInf(const bigfloat & x)
{
  return (x.special == bigfloat::NINF_VAL);
}
bool ispInf(const bigfloat & x)
{
  return (x.special == bigfloat::PINF_VAL);
}
bool isnZero(const bigfloat & x)
{
  return (x.special == bigfloat::NZERO_VAL);
}
bool ispZero(const bigfloat & x)
{
  return (x.special == bigfloat::PZERO_VAL);
}
bool isZero(const bigfloat & x)
{
  return ((x.special == bigfloat::PZERO_VAL) || (x.special == bigfloat::NZERO_VAL));
}
bool isInf(const bigfloat & x)
{
  return ((x.special == bigfloat::PINF_VAL) || (x.special == bigfloat::NINF_VAL));
}



bigfloat pow2(const integer & p)
{
  return bigfloat(1, p);
}
integer ceil(const bigfloat & x)
{
  return to_integer(x, TO_P_INF);
}
integer floor(const bigfloat & x)
{
  return to_integer(x, TO_N_INF);
}





bigfloat::bigfloat()
{
  special = bigfloat::PZERO_VAL;
  exponent = significant = integer_zero;
  precision = 0;
}
bigfloat::bigfloat(const integer & s, const integer & e)
{
  if (sign(s) != 0) {
    significant = s;
    exponent = e;
    special = bigfloat::NOT_VAL;
    precision = s.length();
  }
  else {
    special = bigfloat::PZERO_VAL;
    exponent = significant = integer_zero;
    precision = 0;
  }
}
bigfloat::bigfloat(special_values sp)
{
  special = sp;
}
bigfloat::bigfloat(const integer & a)
{
  special = bigfloat::NOT_VAL;
  significant = a;
  exponent = integer_zero;
  precision = significant.length();

  if (sign(significant) == 0)
    special = bigfloat::PZERO_VAL;



}
bigfloat::bigfloat(long a)
{
  special = bigfloat::NOT_VAL;
  significant = a;
  exponent = integer_zero;
  precision = significant.length();

  if (sign(significant) == 0)
    special = bigfloat::PZERO_VAL;



}
bigfloat::bigfloat(int a)
{
  special = bigfloat::NOT_VAL;
  significant = a;
  exponent = integer_zero;
  precision = significant.length();

  if (sign(significant) == 0)
    special = bigfloat::PZERO_VAL;



}


bigfloat::bigfloat(double d)
{
  long a = (long) d;
  if (a == d) {

    special = bigfloat::NOT_VAL;
    significant = a;
    exponent = integer_zero;
    precision = significant.length();

    if (sign(significant) == 0)
      special = bigfloat::PZERO_VAL;



    return;
  }

  bool is_normalized = is_not_underflowed(d);
  if (!is_normalized)
    d *= pow52;

  int sign_bit;
  long exp_11;
  unsigned long most_sig_20, least_sig_32;
  read_parts(d, sign_bit, exp_11, most_sig_20, least_sig_32);


  if (!exp_11 || exp_11 == 2047) {
    if (d == 0)
      special = (sign_bit) ? bigfloat::NZERO_VAL : bigfloat::PZERO_VAL;
    if (is_infinite(d))
      special = (sign_bit) ? bigfloat::NINF_VAL : bigfloat::PINF_VAL;
    if (is_nan(d))
      special = bigfloat::NAN_VAL;
    return;
  }





  special = bigfloat::NOT_VAL;

  most_sig_20 |= 0x00100000;
  significant = integer(most_sig_20);
  significant = significant << 32;
  significant = significant + integer(least_sig_32);
  if (sign_bit)
    significant = -significant;



  exp_11 -= 1075;
  if (!is_normalized)
    exp_11 -= 52;
  exponent = exp_11;




  precision = significant.length();
}



long sign(const bigfloat & x)
{
  switch (x.special) {
 case bigfloat::NOT_VAL:
	return::sign(x.significant);
      case bigfloat::PZERO_VAL:
      case bigfloat::NZERO_VAL:
	return 0;
      case bigfloat::PINF_VAL:
	return 1;
      case bigfloat::NINF_VAL:
	return -1;
      case bigfloat::NAN_VAL:
	error_handler(1, "sign: NaN has no sign");
	return 1;
  }
  return 0;
}
bigfloat abs(const bigfloat & x)
{
  if (isZero(x))
    return bigfloat::pZero;
  else
    if (sign(x) > 0)
      return  x;
    else
      return -x;
}
integer ilog2(const bigfloat & x)
{
  if (x.precision != 1 + x.significant.zeros())
    return x.precision + x.exponent;
  else
    return (x.precision - 1) + x.exponent;
}



long sign_of_special_value(const bigfloat & x)
{
  if (x.special == bigfloat::NAN_VAL)
    error_handler(1, "sign_of_special_value: want a special value but not NaN");
  if (x.special == bigfloat::NOT_VAL)
    return sign(x);
  if (x.special == bigfloat::PZERO_VAL || x.special == bigfloat::PINF_VAL)
    return 1;
  else
    return -1;
}



bigfloat round(const bigfloat & y, long digits, rounding_modes mode, bool &is_exact, long bias)
{
  if (isSpecial(y) || (mode == EXACT)) {
    is_exact = true;
    return y;
  }
  if ((y.precision <= digits) && (bias == 0)) {
    is_exact = true;
    return y;
  }

  long precision = y.precision - y.significant.zeros();
  bool need_rounding = (digits < precision);
  integer significant = y.significant;
  integer exponent = y.exponent;

  if (need_rounding) {
    switch (mode) {
	case TO_NEAREST:

	  bool is_odd;
	  cut(significant, digits + 1);
	  is_odd = significant.tolong() & 1;
	  cut(significant, digits);
	  if (is_odd) {
	/* $x_2$ starts with a one */
	    if ((precision > digits + 1)
	      || (::sign(bias) == ::sign(significant))
	      || ((bias == 0) && (significant.tolong() & 1))
	      ) {
	  /* $x_1+sign(x_1)$ is the rounded value */
	      if (::sign(significant) > 0)
		significant++;
	      else
		significant--;
	    }
	  }


	  break;
	case TO_ZERO:

	  cut(significant, digits);

	  break;
	case TO_P_INF:


	  cut(significant, digits);
	  if (::sign(significant) > 0)
	    significant++;

	  break;
	case TO_N_INF:


	  cut(significant, digits);
	  if (::sign(significant) < 0)
	    significant--;

	  break;
	case TO_INF:

	  cut(significant, digits);
	  if (::sign(significant) > 0)
	    significant++;
	  else
	    significant--;

	  break;
	case EXACT:
      /* never reached */
	  break;
    }
    exponent += (y.precision - digits);
  }
  else if (bias != 0) {


    long shift = digits - y.precision;
    long bf_sign = ::sign(significant);

    switch (mode) {
	case TO_ZERO:
	  if (::sign(bias) != bf_sign) {

	    if (shift >= 0)
	      significant = significant << shift;
	    else
	      significant = significant >> -shift;
	    exponent -= shift;


	    if (bf_sign == 1)
	      significant--;
	    else
	      significant++;
	  }
	  break;
	case TO_INF:
	  if (::sign(bias) == bf_sign) {

	    if (shift >= 0)
	      significant = significant << shift;
	    else
	      significant = significant >> -shift;
	    exponent -= shift;


	    if (bf_sign == 1)
	      significant++;
	    else
	      significant--;
	  }
	  break;
	case TO_P_INF:
	  if (::sign(bias) == 1) {

	    if (shift >= 0)
	      significant = significant << shift;
	    else
	      significant = significant >> -shift;
	    exponent -= shift;


	    significant++;
	  }
	  break;
	case TO_N_INF:
	  if (::sign(bias) == -1) {

	    if (shift >= 0)
	      significant = significant << shift;
	    else
	      significant = significant >> -shift;
	    exponent -= shift;


	    significant--;
	  }
	  break;
      /* in the cases |TO_NEAREST| and |EXACT| we have to do nothing */
	case TO_NEAREST:
	case EXACT:;
    }


  }
  is_exact = !need_rounding && (bias == 0);

  if (significant.length() > digits) {
    long z = significant.zeros();
    significant = significant >> z;
    exponent += z;
  }

  return bigfloat(significant, exponent);
}


integer to_integer(bigfloat x, rounding_modes rmode)
{
  if ((isNaN(x)) || (isInf(x)))
    error_handler(1,
      "to_integer : special values cannot be converted to integer");

  if (isZero(x))
    return integer_zero;

  if (!is_small(x.exponent))
    error_handler(1, "error to_integer: number too large");

  long x_exponent = x.exponent.tolong();
  if (x.precision + x_exponent >= 1) {
    long length = max(1, x.precision + x_exponent);
    x = round(x, length, rmode);
    return (x.significant << x.exponent.tolong());
  }
  else {
    switch (rmode) {
	case TO_NEAREST:
	case EXACT:
	  if (x.precision + x_exponent < 0)
	    return integer_zero;
	  else
	    return sign(x);

	case TO_ZERO:
	  return integer_zero;

	case TO_INF:
	  return sign(x);

	case TO_P_INF:
	  if (sign(x) > 0)
	    return integer_one;
	  else
	    return integer_zero;

	case TO_N_INF:
	  if (sign(x) > 0)
	    return integer_zero;
	  else
	    return -integer_one;

    }
  }
  return 0;

}



double to_double(const bigfloat & x)
{


  if (isSpecial(x)) {
    if (ispZero(x))
      return pZero_double;
    if (isnZero(x))
      return nZero_double;
    if (isNaN(x))
      return NaN_double;
    if (ispInf(x))
      return pInf_double;
    if (isnInf(x))
      return nInf_double;
  }






  bigfloat y = round(x, 53, TO_NEAREST);
  double significant = y.significant.to_double();
  long exponent;
  if (is_small(y.exponent))
    exponent = y.exponent.tolong();
  else
    exponent = sign(x) * 2000;
  if (exponent > 1023)
    return sign(x) * pInf_double;
  if (exponent >= -1022)
    return power_two(exponent) * significant;
  long logx = exponent + y.precision;
  if (logx >= -1074)
    return double_min * (power_two(1022 + exponent) * significant);
  else
    return sign(x) * pZero_double;
}



bigfloat operator + (const bigfloat & a, const bigfloat & b) {
  return add(a, b);
}

bigfloat operator - (const bigfloat & a, const bigfloat & b) {
  return sub(a, b);
}

bigfloat operator *(const bigfloat & a, const bigfloat & b) {
  return mul(a, b);
}

bigfloat operator / (const bigfloat & a, const bigfloat & b) {
  return div(a, b);
}

bigfloat operator += (bigfloat & a, const bigfloat & b) {
  a = add(a, b);
  return a;
}

bigfloat operator -= (bigfloat & a, const bigfloat & b) {
  a = sub(a, b);
  return a;
}

bigfloat operator *= (bigfloat & a, const bigfloat & b) {
  a = mul(a, b);
  return a;
}

bigfloat operator /= (bigfloat & a, const bigfloat & b) {
  a = div(a, b);
  return a;
}



bigfloat add(const bigfloat & x, const bigfloat & y, long prec, rounding_modes mode,
   bool &is_exact)
{


  if (isSpecial(x) || isSpecial(y)) {
    if (isNaN(x) || isNaN(y))
      return bigfloat::NaN;

    if (isZero(x) && isZero(y)) {
      if (isnZero(x) && isnZero(y))
	return bigfloat::nZero;
      else
	return bigfloat::pZero;
    }

    if (isZero(x))
      return y;
    if (isZero(y))
      return x;

    if (isInf(x) && isInf(y)) {
      if (sign_of_special_value(x) == sign_of_special_value(y))
	return x;
      else
	return bigfloat::NaN;
    }
    if (isInf(x))
      return x;

/* it is obvious that |y| has to be $\infty$ */
    return y;
  }




  bigfloat sum;
  int bias;

  integer exp_diff = x.exponent - y.exponent;
  bool no_exact_sum, x_is_bigger;
  long exp_diff_long = exp_diff.tolong();
  if (is_small(exp_diff_long)) {
    long log_diff = exp_diff_long + (x.precision - y.precision);
    x_is_bigger = (log_diff >= 0);
    if (x_is_bigger)
      no_exact_sum = (log_diff >= max(prec + 1, x.precision));
    else
      no_exact_sum = (-log_diff >= max(prec + 1, y.precision));

  }
  else {
    integer log_diff = exp_diff + (x.precision - y.precision);
    x_is_bigger = (log_diff >= 0);
    if (x_is_bigger)
      no_exact_sum = (log_diff >= max(prec + 1, x.precision));
    else
      no_exact_sum = (-log_diff >= max(prec + 1, y.precision));
  }
  if ((mode != EXACT) && (no_exact_sum)) {
    if (x_is_bigger) {
      sum = x;
      bias = sign(y);
    }
    else {
      sum = y;
      bias = sign(x);
    }
  }
  else {
    if (!exp_diff.islong())
      error_handler(1, "bigfloat::add(): exponential difference too large");

    if (sign(exp_diff) >= 0) {
      sum = y;
      sum.significant += x.significant << exp_diff.tolong();
    }
    else {
      sum = x;
      sum.significant += y.significant << (-exp_diff.tolong());
    }
    if (sign(sum.significant) == 0)
      sum.special = bigfloat::PZERO_VAL;
    sum.precision = sum.significant.length();
    bias = 0;
  }




  return round(sum, prec, mode, is_exact, bias);
}



bigfloat sub(const bigfloat & a, const bigfloat & b, long prec,
   rounding_modes mode, bool &is_exact)
{
  return add(a, -b, prec, mode, is_exact);
}



bigfloat operator - (const bigfloat & a) {
  if (isSpecial(a)) {
    if (ispZero(a))
      return bigfloat::nZero;
    if (isnZero(a))
      return bigfloat::pZero;
    if (ispInf(a))
      return bigfloat::nInf;
    if (isnInf(a))
      return bigfloat::pInf;
    return bigfloat::NaN;
  }
  return bigfloat(-a.significant, a.exponent);
}




bigfloat mul(const bigfloat & a, const bigfloat & b, long prec,
   rounding_modes mode, bool &is_exact)
{


  long sign_result;

  if ((isSpecial(a)) || (isSpecial(b))) {
    if ((isNaN(a)) || (isNaN(b)))
      return bigfloat(bigfloat::NAN_VAL);
    if ((isZero(a) && isInf(b)) || (isInf(a) && isZero(b)))
      return bigfloat(bigfloat::NAN_VAL);

    sign_result = sign_of_special_value(a) * sign_of_special_value(b);

    if (isZero(a) || isZero(b)) {
      if (sign_result == 1)
	return bigfloat(bigfloat::PZERO_VAL);
      else
	return bigfloat(bigfloat::NZERO_VAL);
    }
    if (isInf(a) || isInf(b)) {
      if (sign_result == 1)
	return bigfloat(bigfloat::PINF_VAL);
      else
	return bigfloat(bigfloat::NINF_VAL);
    }
  }



  bigfloat result(a.significant * b.significant, a.exponent + b.exponent);
  return round(result, prec, mode, is_exact);
}



bigfloat div(const bigfloat & aa, const bigfloat & b, long prec,
   rounding_modes mode, bool &is_exact)
{
  bigfloat result;
  long bias;

  bigfloat a = aa;


  if ((isSpecial(a)) || (isSpecial(b))) {
    long sign_result;

    if ((isNaN(a)) || (isNaN(b)))
      return bigfloat(bigfloat::NAN_VAL);
    if (((isZero(a)) && (isZero(b))) || ((isInf(a)) && (isInf(b))))
      return bigfloat(bigfloat::NAN_VAL);
    sign_result = sign_of_special_value(a) * sign_of_special_value(b);
    if ((isInf(a)) || (isZero(b))) {
      if (sign_result == 1)
	return bigfloat(bigfloat::PINF_VAL);
      else
	return bigfloat(bigfloat::NINF_VAL);
    }

/* it is clear that |isZero(a)||(isInf(b)| */

    if (sign_result == 1)
      return bigfloat(bigfloat::PZERO_VAL);
    else
      return bigfloat(bigfloat::NZERO_VAL);
  }





  long d = prec + b.significant.length() - a.significant.length() + 1;
  if (d > 0) {
    a.significant = a.significant << d;
    a.exponent -= d;
  }



  result.special = bigfloat::NOT_VAL;
  result.significant = a.significant / b.significant;

  result.exponent = a.exponent - b.exponent;
  result.precision = result.significant.length();

  integer R = a.significant - b.significant * result.significant;
  bias = sign(R) * sign(b.significant);
  if (sign(R) != 0)
    is_exact = false;
  if (mode == EXACT)
    mode = TO_NEAREST;



  return round(result, prec, mode, is_exact, bias);
}



bigfloat sqrt(const bigfloat & a, long prec, rounding_modes mode,
   bool &is_exact, const bigfloat & start)
{
  bigfloat result;
  bigfloat start_value;
  bool exact = true;


  if (isZero(a))
    return a;
  if (sign(a) < 0) {
    is_exact = false;
    return bigfloat::NAN_VAL;
  }
  if (isSpecial(a)) {
    is_exact = false;
    return a;
  }



  if (isSpecial(start) || start <= 0) {
    bigfloat tmp = bigfloat(a.significant, -a.precision);
    double d = to_double(tmp);
    integer log_a = a.precision + a.exponent;
    if (log_a.tolong() & 1)
      d *= 2;
    start_value = sqrt(d);
    start_value.exponent += log_a >> 1;
  }
  else
    start_value = start;



  bigfloat it;			/* iteration value */
  bigfloat ne;			/* next iteration value */

  bigfloat square, sum;
  ne = start_value;
  long t_prec = ne.precision;
  long p_bound = prec + 1;

  while (1 == 1) {
    exact = true;
    t_prec *= 2;
    if (t_prec > p_bound)
      t_prec = p_bound;
    it = ne;
    square = mul(it, it, t_prec, EXACT);
    sum = add(square, a, t_prec, EXACT);
    ne = div(sum, it, t_prec, TO_P_INF, exact);
    ne.exponent--;
    if ((t_prec == p_bound) || exact)
      if (ne == it)
	break;
  }
  result = ne;
  if (!exact)
    is_exact = false;



  if (exact)
    result = round(result, prec, mode, is_exact);
  else {
    switch (mode) {
	case TO_NEAREST:
	  result = round(result, prec, TO_ZERO, is_exact);
	  break;
	case EXACT:
	  error_handler(1, "bigfloat::sqrt : the EXACT-rounding mode is not provided!");
	case TO_P_INF:
	case TO_INF:
	  result = round(result, prec, mode, is_exact);
	  break;
	case TO_N_INF:
	case TO_ZERO:
	  result = round(result, prec, mode, is_exact, -1);
    }
    is_exact = false;
  }


  return result;
}


bigfloat sqrt_d(const bigfloat & a, long prec, int d)
{
  if (d < 2)
    error_handler(1, "sqrt_d called with d < 2");
  bigfloat result;
  bigfloat start_value;

  if (isZero(a))
    return a;
  if (sign(a) < 0) {
    return bigfloat::NAN_VAL;
  }
  if (isSpecial(a)) {
    return a;
  }



  bigfloat tmp = bigfloat(a.significant, -a.precision);
  double dbl = to_double(tmp);
  integer log_a = a.precision + a.exponent;
  long log_a_long = log_a.tolong();
  if (log_a_long % d)
    dbl *= power_two(log_a_long % d);
  start_value = pow(dbl, 1 / double (d));
  start_value.exponent += log_a / d;



  bigfloat it;			/* iteration value */
  bigfloat ne;			/* next iteration value */

  bigfloat D = d, D_minus_one = d - 1;
  bigfloat pow_it;
  long logd = (int) ceil(log10(double(d)) / log10(2.0));
  ne = start_value;
  long t_prec = 51;
  long p_bound = prec + 1;
  bool it_has_full_precision = false;

  long old_prec = bigfloat::set_precision(53);
  rounding_modes old_mode = bigfloat::set_rounding_mode(TO_INF);
  while (1 == 1) {
    if (t_prec == p_bound)
      it_has_full_precision = true;
    t_prec *= 2;
    if (t_prec > p_bound)
      t_prec = p_bound;
    bigfloat::set_precision(t_prec + 5);
    it = ne;
    pow_it = Powl(it, d - 1, t_prec + logd + 5, TO_INF);
    ne = (it * D_minus_one + a / pow_it) / D;
    ne = round(ne, prec + 2, TO_INF);
    if (t_prec == p_bound)
      if (ne >= it)
	if (it_has_full_precision)
	  break;
  }
  bigfloat::set_precision(old_prec);
  bigfloat::set_rounding_mode(old_mode);
  result = round(it, prec + 1, TO_NEAREST);



#ifdef DEBUG_BIGFLOAT
  if (abs(Powl(result, d, prec + 20) - a) / (d * a) > pow2(-prec))
    error_handler(1, "error sqrt_d: precision not achieved");
#endif
  return result;
}



bool operator == (const bigfloat & a, const bigfloat & b) {

  if (isSpecial(a) || isSpecial(b)) {
    if (isZero(a) && isZero(b))
      return true;
    if (isNaN(a) || isNaN(b))
      error_handler(1, "bigfloat::operator == : NaN case occurred");
    return (a.special == b.special);
  }




  integer exp_diff = a.exponent - b.exponent;
  long prec_diff = a.precision - b.precision;
  integer log_diff = exp_diff + prec_diff;
  if (sign(log_diff) != 0)
    return false;
  if (!exp_diff.islong())
    return false;
  long a_shift = exp_diff.tolong();
  if (a_shift > 0)
    return ((a.significant << a_shift) == b.significant);
  if (a_shift < 0)
    return (a.significant == (b.significant << -a_shift));
  return (a.significant == b.significant);
}


bool operator > (const bigfloat & a, const bigfloat & b) {

  if (isSpecial(a) || isSpecial(b)) {
    if (isNaN(a) || isNaN(b))
      error_handler(1, "bigfloat::operator > : NaN case occurred!");
    if (isZero(a) && isZero(b))
      return false;
    if (ispInf(a))
      return !ispInf(b);
    if (isnInf(b))
      return !isnInf(a);
    return (sign(a) > sign(b));
  }




  int sign_a = sign(a.significant);
  int sign_b = sign(b.significant);

  if (sign_a != sign_b)
    return (sign_a > sign_b);

  integer exp_diff = a.exponent - b.exponent;
  integer log_diff = exp_diff + (a.precision - b.precision);
  int sign_log_diff = sign(log_diff);
  if (sign_log_diff != 0)
    return (sign_a * sign_log_diff > 0);

  int sign_exp_diff = sign(exp_diff);
  if (sign_exp_diff >= 0) {
    integer sig = a.significant << (exp_diff.tolong());
    return (sig > b.significant);
  }
  integer sig = b.significant << ((-exp_diff).tolong());
  return (a.significant > sig);
}



bool operator != (const bigfloat & a, const bigfloat & b) {
  return !(a == b);
}
bool operator >= (const bigfloat & a, const bigfloat & b) {
  return ((a > b) || (a == b));
}
bool operator < (const bigfloat & a, const bigfloat & b) {
  return (!(a >= b));
}
bool operator <= (const bigfloat & a, const bigfloat & b) {
  return (!(a > b));
}



ostream & operator << (ostream & os, const bigfloat & b) {
  if (isSpecial(b)) {
    if (isNaN(b))
      return os << "NaN";
    if (ispInf(b))
      return os << "+Inf";
    if (isnInf(b))
      return os << "-Inf";
    if (ispZero(b))
      return os << "0";
    if (isnZero(b))
      return os << "-0";
  }
  int sign_b = sign(b.significant);
  if (bigfloat::output_mode == bigfloat::BIN_OUT) {
    if (sign_b < 0)
      os << "-";
    os << "0.";
    if (sign_b >= 0)
      binout(os, b.significant);
    else
      binout(os, -b.significant);
    os << "E";
    if (b.exponent < 0)
      os << "-";
    else
      os << "+";
    binout(os, b.exponent);
  }
  if (bigfloat::output_mode == bigfloat::DEC_OUT) {
    if (b == to_integer(b))
      os << to_integer(b);
    else {
      bigfloat bb = b;
      bb.normalize();
      long prec = (long) floor(log10(2.0) * bb.significant.length());
      if (prec < 7)
	prec = 7;
      decimal_output(os, b, prec);
    }
  }

  return os;
}


istream & operator >> (istream & is, bigfloat & b) {
  char tmp[bin_maxlen];
  is >> tmp;
  b = outofchar(tmp);
  return is;
}
