/*******************************************************************************
+
+  LEDA 3.5
+
+  _real.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 "real.w" by CTANGLE
 * (Version 3.1 [km2c])
 */

#include <LEDA/real.h>
#include <math.h>
#include <LEDA/floatingpoint.h>

const double MaxError = power_two(915);
const double Infinity = pInf_double;
const double eps = power_two(-53);
const double two_eps = power_two(-52);
const double one_minus_two_eps = 1 - power_two(-52);
const double correction = 1 + 8 * eps;
const double MinDbl = double_min;
const double twoMinDbl = 2 * double_min;



int isdouble(const integer & x)
{
  return (x == integer(to_double(bigfloat(x))));
}
int isdouble(const bigfloat & x)
{
  return (x == bigfloat(to_double(x)));
}


const integer zero_integer(0);


bigfloat ldexp_bf(const bigfloat & x, const integer & k)
{
  return bigfloat(x.get_significant(), x.get_exponent() + k);
}
integer Maximum(integer x, integer y)
{
  return (x > y ? x : y);
}




struct sep_bound;
class real_rep {
  friend class real;

  struct sep_bound {
    integer U, L;
     sep_bound() {
    }
  };



  struct app_bf {
    bigfloat NUM_BF;
    bigfloat ERROR_BF;
     app_bf() {
    }
     app_bf(const bigfloat & x, const bigfloat & error_x = bigfloat::pZero) {
      NUM_BF = x;
      ERROR_BF = error_x;
    }
  };



  enum constructortype {
    DOUBLE = 0, BIGFLOAT = 1, NEGATION = 2,
    SQUAREROOT = 3, ADDITION = 4, SUBTRACTION = 5, MULTIPLICATION = 6,
  DIVISION = 7, ROOT = 8};



  enum statustype {
  EMPTY = 0, CLEARED = 1, VISITED = 2};



  double NUM;
  double ERROR;
  app_bf *APP_BF_PTR;
  constructortype CON;
  real_rep *OP1, *OP2;
  int count;
  int d;
  long status;
  sep_bound *sep_bound_ptr;

  real_rep();
  real_rep(double);
  real_rep(const bigfloat &);

  ~real_rep();

  LEDA_MEMORY(real_rep)
    void init_app_bf();
  friend double to_double(const real &);
  friend bigfloat to_bigfloat(const real &);
  void adjust_dbl();
  void compute_op(long);
  void compute_approximation(long);
  int sign_with_separation_bound(const integer &, bool = false);
  void guarantee_bound_two_to(const integer &);
  inline bigfloat & num_bf() {
    return APP_BF_PTR->NUM_BF;
  }
  inline bigfloat & error_bf() {
    return APP_BF_PTR->ERROR_BF;
  }
  inline bool exact() {
    return isZero(error_bf());
  }
  inline integer num_exp() {
    return ilog2(num_bf());
  }
  inline integer err_exp() {
    return ilog2(error_bf());
  }


  void compute_concrete_parameters();



  void compute_parameters();


  void estimate_degree(integer & D, long &k, bool &contains_divisions);


  void clear_visited_marks();


  void compute_concrete_system_bound(integer & concrete_system_bound);


  void improve_add(const integer &);
  void improve_sub(const integer &);
  void improve_mul(const integer &);
  void improve_div(const integer &);
  void improve_sqrt(const integer &);
  void improve_root(const integer &);



  void compute_mul_error(bigfloat &);
  void compute_div_error(bigfloat &);
  void compute_sqrt_error(bigfloat &);
  void compute_root_error(bigfloat &);




  friend ostream & operator << (ostream &, const real &);

  friend bool operator < (const real &, const real &);
  friend bool operator <= (const real &, const real &);
  friend bool operator > (const real &, const real &);
  friend bool operator >= (const real &, const real &);
  friend bool operator == (const real &, const real &);
  friend bool operator != (const real &, const real &);

  friend real operator + (const real &, const real &);
  friend real operator - (const real &, const real &);
  friend real operator *(const real &, const real &);
  friend real operator / (const real &, const real &);
  friend real operator - (const real &);
  friend real sqrt(const real &);
  friend real root(const real &, int);
};



real::real()
{
  is_double = true;
  value = 0;
  PTR = nil;
}
real::real(int x)
{
  is_double = true;
  value = x;
  PTR = nil;
}
real::real(double x)
{
  is_double = true;
  value = x;
  PTR = nil;
}
real::real(const bigfloat & x)
{
  if (is_double = isdouble(x)) {
    value = ::to_double(x);
    PTR = nil;
  }
  else
    PTR = new real_rep(x);
}
real::real(const integer & x)
{
  if (is_double = isdouble(x)) {
    value = x.todouble();
    PTR = nil;
  }
  else
    PTR = new real_rep(bigfloat(x));
}
real::real(const rational & x)
{
  real OP1 = *new real_rep(x.numerator());
  real OP2 = *new real_rep(x.denominator());
  PTR = nil;
  *this = OP1 / OP2;
}
real::real(real_rep & x_rep)
{
  is_double = false;
  PTR = &x_rep;
}


real::real(const real & x)
{
  if (is_double = x.is_double) {
    value = x.value;
    PTR = nil;
  }
  else {
    x.PTR->count++;
    PTR = x.PTR;
  }
}
real & real::operator = (const real & x) {

/* we avoid problems by ruling out the assignment |x = x;| */
  if (this == &x)
    return (*this);

/*
 * first the actions that must anyway be performed, in the respective cases
 * of |is_double|
 */
  if (is_double = x.is_double)
    value = x.value;
  else
    x.PTR->count++;

/*
 * second, an existing reference in PTR must always be deleted; if |x.PTR ==
 * PTR != nil| and |is_double == false|, this leaves the number of references
 * to |*PTR| unchanged, because we increased |x.PTR| before, and for
 * |is_double == true| we in fact lose one reference because |PTR| is set to
 * |nil| at the end
 */
  if (PTR)
    if (--PTR->count == 0)
      delete PTR;

/* last, we set |this->PTR| */
  if (is_double)
    PTR = nil;
  else
    PTR = x.PTR;

  return (*this);
}


real_rep::real_rep()
{
  APP_BF_PTR = nil;
  sep_bound_ptr = nil;
  count = 1;
  status = 0;
}
real_rep::real_rep(double x)
{
  APP_BF_PTR = nil;
  sep_bound_ptr = nil;
  OP1 = OP2 = nil;
  CON = DOUBLE;
  NUM = x;
  ERROR = 0;
  count = 1;
  status = 0;
}
real_rep::real_rep(const bigfloat & x)
{
  APP_BF_PTR = new app_bf(x, 0);
  sep_bound_ptr = nil;
  OP1 = OP2 = nil;
  CON = BIGFLOAT;
  adjust_dbl();
  count = 1;
  status = 0;
}


real::~real()
{
  if (PTR)
    if (--PTR->count == 0)
      delete PTR;
}
real_rep::~real_rep()
{
  if (APP_BF_PTR)
    delete APP_BF_PTR;
  if (sep_bound_ptr)
    delete sep_bound_ptr;
  if (OP1 && (--OP1->count == 0))
    delete OP1;
  if (OP2 && (--OP2->count == 0))
    delete OP2;
}




void real_rep::init_app_bf()
{

  if (APP_BF_PTR)
    return;
  if (OP1)
    OP1->init_app_bf();
  if (OP2)
    OP2->init_app_bf();
  if (is_finite(ERROR)) {
    bigfloat error = ERROR * float_abs(NUM) * correction;
    bigfloat approximation = NUM;
    APP_BF_PTR = new app_bf(approximation, error);
    return;
  }
  APP_BF_PTR = new app_bf();
  compute_op(53);
}


void real_rep::adjust_dbl()
{
  double n_head = to_double(num_bf());
  double e_head = to_double(error_bf());
  if (n_head != 0) {
    NUM = n_head;
    ERROR = eps + Max(MinDbl, e_head) / float_abs(n_head);
  }
  else {
    NUM = e_head + 2 * MinDbl;
    ERROR = 2;
  }
  ERROR = ERROR * correction;
  if ((ERROR >= MaxError) || (is_infinite(n_head))) {
    NUM = 1;
    ERROR = Infinity;
  }
}


void real::adjust() const
{
  if (is_double)
    return;
  if (PTR->exact() && isdouble(PTR->num_bf())) {
    ((real &) * this).is_double = true;
    ((real &) * this).value = ::to_double(PTR->num_bf());
    return;
  }
  PTR->adjust_dbl();
}



void decimal_output(ostream & os, bigfloat b, long dec_prec,
   rounding_modes mode = TO_NEAREST);

ostream & operator << (ostream & out, const real & x) {
  if (x.get_double_error() == 0) {
    double d = to_double(x);
    if (d == floor(d))
      out << "[" << d << "," << d << "]";
    return out;
  }
  x.PTR->init_app_bf();
  if (x.get_bigfloat_error() == 0) {
    bigfloat b = to_bigfloat(x);
    bigfloat I = to_integer(b);
    if (b == I)
      out << "[" << I << "," << I << "]";
    return out;
  }
  bigfloat low_bound, upp_bound;
  long prec = (ilog2(x.PTR->num_bf()) - ilog2(x.PTR->error_bf())).tolong();
  if (prec < 53)
    prec = 53;
  low_bound =
    sub(x.to_bigfloat(), x.PTR->error_bf(), prec, TO_N_INF);
  upp_bound =
    add(x.to_bigfloat(), x.PTR->error_bf(), prec, TO_P_INF);
  long dec_prec = (long) floor(log10(2.0) * prec);
  out << "[";
  decimal_output(out, low_bound, dec_prec, TO_N_INF);
  out << ",";
  decimal_output(out, upp_bound, dec_prec, TO_P_INF);
  out << "]";
  return out;
}
istream & operator >> (istream & in, real & x) {
  double x_num;
  in >> x_num;
  x = real(x_num);
  return in;
}


double to_double(const real & x)
{
  if (x.is_double)
    return x.value;
  if (x.PTR->APP_BF_PTR)
    x.PTR->adjust_dbl();
  return x.PTR->NUM;
}
double real::to_double() const
{
  return::to_double(*this);
}
double real::get_double_error() const
{
  if (is_double)
    return 0;
  if (PTR->APP_BF_PTR)
    PTR->adjust_dbl();
  return PTR->ERROR;
}
bigfloat to_bigfloat(const real & x)
{
  if (x.is_double)
    return x.value;
  else {
    x.PTR->init_app_bf();
    return x.PTR->num_bf();
  }
}
bigfloat real::to_bigfloat() const
{
  return::to_bigfloat(*this);
}
bigfloat real::get_bigfloat_error() const
{
  if (is_double)
    return 0;
  else {
    PTR->init_app_bf();
    return PTR->error_bf();
  }
}





real operator + (const real & x, const real & y) {

  double NUM_result;
  if (x.is_double && y.is_double) {
    NUM_result = x.value + y.value;
    if ((NUM_result - x.value == y.value)
      && ((NUM_result - y.value) == x.value)
      )
      return NUM_result;
  }



  if (!x.PTR)
    ((real &) x).PTR = new real_rep(x.value);
  if (!y.PTR)
    ((real &) y).PTR = new real_rep(y.value);
  double &x_num = x.PTR->NUM;
  double &y_num = y.PTR->NUM;


  real_rep & z_rep = *new real_rep();

  z_rep.OP1 = x.PTR;
  z_rep.OP2 = y.PTR;
  z_rep.OP1->count++;
  z_rep.OP2->count++;
  double &z_num = z_rep.NUM;


  z_rep.CON = real_rep::ADDITION;
  z_num = x_num + y_num;

  if (is_not_underflowed(z_num)) {
    double term_x = x_num / z_num;
    clear_sign_bit(term_x);
    double term_y = y_num / z_num;
    clear_sign_bit(term_y);
    z_rep.ERROR = term_x * x.PTR->ERROR + term_y * y.PTR->ERROR + eps;
  }
  else {

    z_num = float_abs(x_num) * x.PTR->ERROR + float_abs(y_num) * y.PTR->ERROR + 4 * MinDbl;
    z_rep.ERROR = 2;
  }
  z_rep.ERROR *= correction;



  if (!((is_finite(z_num)) && (z_rep.ERROR < MaxError))) {
    z_rep.ERROR = Infinity;
    z_num = 1;
  }


  return z_rep;
}
real operator - (const real & x, const real & y) {

  double NUM_result;
  if (x.is_double && y.is_double) {
    NUM_result = x.value - y.value;
    if ((NUM_result + y.value == x.value)
      && (x.value - NUM_result == y.value)
      )
      return NUM_result;
  }




  if (!x.PTR)
    ((real &) x).PTR = new real_rep(x.value);
  if (!y.PTR)
    ((real &) y).PTR = new real_rep(y.value);
  double &x_num = x.PTR->NUM;
  double &y_num = y.PTR->NUM;


  real_rep & z_rep = *new real_rep();

  z_rep.OP1 = x.PTR;
  z_rep.OP2 = y.PTR;
  z_rep.OP1->count++;
  z_rep.OP2->count++;
  double &z_num = z_rep.NUM;


  z_rep.CON = real_rep::SUBTRACTION;
  z_num = x_num - y_num;

  if (is_not_underflowed(z_num)) {
    double term_x = x_num / z_num;
    clear_sign_bit(term_x);
    double term_y = y_num / z_num;
    clear_sign_bit(term_y);
    z_rep.ERROR = term_x * x.PTR->ERROR + term_y * y.PTR->ERROR + eps;
  }
  else {

    z_num = float_abs(x_num) * x.PTR->ERROR + float_abs(y_num) * y.PTR->ERROR + 4 * MinDbl;
    z_rep.ERROR = 2;
  }
  z_rep.ERROR *= correction;



  if (!((is_finite(z_num)) && (z_rep.ERROR < MaxError))) {
    z_rep.ERROR = Infinity;
    z_num = 1;
  }


  return z_rep;
}

real operator *(const real & x, const real & y) {

  double NUM_result;
  if (x.is_double && y.is_double) {
    if (is_single_precision(x.value) && is_single_precision(y.value)) {
      NUM_result = x.value * y.value;
      if (is_finite(NUM_result) && is_not_underflowed(NUM_result))
	return NUM_result;
    }
  }




  if (!x.PTR)
    ((real &) x).PTR = new real_rep(x.value);
  if (!y.PTR)
    ((real &) y).PTR = new real_rep(y.value);
  double &x_num = x.PTR->NUM;
  double &y_num = y.PTR->NUM;


  real_rep & z_rep = *new real_rep();

  z_rep.OP1 = x.PTR;
  z_rep.OP2 = y.PTR;
  z_rep.OP1->count++;
  z_rep.OP2->count++;
  double &z_num = z_rep.NUM;


  z_rep.CON = real_rep::MULTIPLICATION;
  z_num = x_num * y_num;

  double &qx = x.PTR->ERROR;
  double &qy = y.PTR->ERROR;
  if (is_not_underflowed(z_num)) {
    z_rep.ERROR = qx + qy + qx * qy + eps;
  }
  else {
    if ((x_num != 0) && (y_num != 0)) {
      z_num = MinDbl * (1 + qx) * (1 + qy);
      z_rep.ERROR = 2;
    }
    else {
      delete & z_rep;
      return 0;
    }
  }
  z_rep.ERROR *= correction;



  if (!((is_finite(z_num)) && (z_rep.ERROR < MaxError))) {
    z_rep.ERROR = Infinity;
    z_num = 1;
  }


  return z_rep;
}


real operator / (const real & x, const real & y) {

  if (y.is_double) {
    if (y.value == 0)
      error_handler(1, "real::operator/:Division by zero");
  }
  else if ((y.PTR->NUM == 0) && is_finite(y.PTR->ERROR))
    error_handler(1, "real::operator/:Division by zero");




  double NUM_result;
  if (x.is_double && y.is_double)
    if (is_single_precision(y.value)) {
      NUM_result = x.value / y.value;
      if (is_single_precision(NUM_result)
	&& NUM_result * y.value == x.value
	)
	if (is_finite(NUM_result) && is_not_underflowed(NUM_result))
	  return NUM_result;
    }



  if (!x.PTR)
    ((real &) x).PTR = new real_rep(x.value);
  if (!y.PTR)
    ((real &) y).PTR = new real_rep(y.value);
  double &x_num = x.PTR->NUM;
  double &y_num = y.PTR->NUM;


  real_rep & z_rep = *new real_rep();

  z_rep.OP1 = x.PTR;
  z_rep.OP2 = y.PTR;
  z_rep.OP1->count++;
  z_rep.OP2->count++;
  double &z_num = z_rep.NUM;


  z_rep.CON = real_rep::DIVISION;
  z_num = x_num / y_num;

  if (y.PTR->ERROR >= 1) {
    z_rep.ERROR = Infinity;
    z_num = 1;
    return z_rep;
  }
  if (is_not_underflowed(z_num))
    z_rep.ERROR = (x.PTR->ERROR + y.PTR->ERROR) / (1 - y.PTR->ERROR)
      + eps;
  else {
    if (x_num != 0) {
      z_num = MinDbl * (1 + x.PTR->ERROR) / (1 - y.PTR->ERROR);
      z_rep.ERROR = 2;
    }
    else {
      delete & z_rep;
      return 0;
    }
  }
  z_rep.ERROR *= correction;




  if (!((is_finite(z_num)) && (z_rep.ERROR < MaxError))) {
    z_rep.ERROR = Infinity;
    z_num = 1;
  }


  return z_rep;
}


real sqrt(const real & x)
{

  if ((x.is_double) && (x.value < 0))
    error_handler(1, "real::operator sqrt:negative argument");




  double NUM_result;
  if (x.is_double) {
    NUM_result = sqrt(x.value);
    if (is_single_precision(NUM_result))
      return NUM_result;
  }



  if (!x.PTR)
    ((real &) x).PTR = new real_rep(x.value);
  real_rep & z_rep = *new real_rep();
  z_rep.OP1 = x.PTR;
  z_rep.OP1->count++;
  z_rep.OP2 = nil;
  z_rep.CON = real_rep::SQUAREROOT;
  if (x.PTR->ERROR >= 1) {
    z_rep.ERROR = Infinity;
    z_rep.NUM = MinDbl;
    return real(z_rep);
  }
  if (x.PTR->NUM < 0)
    error_handler(1, "real::operator sqrt:negative argument");
  z_rep.NUM = sqrt(x.PTR->NUM);

  z_rep.ERROR = x.PTR->ERROR + eps;
  z_rep.ERROR *= correction;



  return z_rep;
}


real root(const real & x, int d)
{

  if ((x.is_double) && (x.value < 0))
    error_handler(1, "real::operator root:negative argument");



  if (!x.PTR)
    ((real &) x).PTR = new real_rep(x.value);
  real_rep & z_rep = *new real_rep();
  z_rep.OP1 = x.PTR;
  z_rep.OP1->count++;
  z_rep.OP2 = nil;
  z_rep.CON = real_rep::ROOT;
  z_rep.d = d;
  if (x.PTR->ERROR >= 1) {
    z_rep.ERROR = Infinity;
    z_rep.NUM = MinDbl;
    return real(z_rep);
  }
  if (x.PTR->NUM < 0)
    error_handler(1, "real::operator root:negative argument");
  z_rep.NUM = pow(x.PTR->NUM, 1 / double (d));

  double delta;
  if (x.PTR->NUM > 1)
    delta = two_eps * (1 + get_unbiased_exponent(x.PTR->NUM));
  else
    delta = two_eps;
  if (delta < 1)
    z_rep.ERROR = delta + x.PTR->ERROR / ((1 - x.PTR->ERROR) * d * (1 - delta));
  else
    z_rep.ERROR = Infinity;
  z_rep.ERROR *= correction;



  return z_rep;
}


real operator - (const real & x) {
  if (x.is_double)
    return -x.value;
  real_rep & z_rep = *new real_rep();
  z_rep.CON = real_rep::NEGATION;
  z_rep.OP1 = x.PTR;
  z_rep.OP1->count++;
  z_rep.OP2 = nil;
  z_rep.NUM = -x.PTR->NUM;
  z_rep.ERROR = x.PTR->ERROR;
  return z_rep;
}


real operator += (real & x, const real & y) {
  x = x + y;
  return x;
}
real operator -= (real & x, const real & y) {
  x = x - y;
  return x;
}
real operator *= (real & x, const real & y) {
  x = x * y;
  return x;
}




int real::sign() const
{
  if (is_double)
    return float_sign(value);
  else
    return PTR->sign_with_separation_bound(zero_integer, false);
}
int real::sign(const integer & q) const
{
  if (is_double)
    return float_sign(value);
  else
    return PTR->sign_with_separation_bound(q, true);
}
int real::sign(long p) const
{
  if (is_double)
    return float_sign(value);
  else
    return PTR->sign_with_separation_bound(p, true);
}


void real_rep::compute_concrete_parameters()
{
  sep_bound *ptr = sep_bound_ptr;
  sep_bound *ptr1 = nil;
  if (OP1)
    ptr1 = OP1->sep_bound_ptr;
  sep_bound *ptr2 = nil;
  if (OP2)
    ptr2 = OP2->sep_bound_ptr;

  switch (CON) {
      case DOUBLE:
      case BIGFLOAT:
	{
	  bigfloat & num = APP_BF_PTR->NUM_BF;
	  ptr->L = zero_integer;
	  if (isZero(num))
	    ptr->U = zero_integer;
	  else {
	    if (sign(num.get_exponent()) >= 0)
	      ptr->U = ilog2(num) + 1;
	    else {
	      ptr->U = num.get_significant_length() + 1;
	      ptr->L = -num.get_exponent();
	    }
	  }
	  break;
	}
      case NEGATION:
	ptr->U = ptr1->U;
	ptr->L = ptr1->L;
	break;
      case MULTIPLICATION:
	ptr->U = ptr1->U + ptr2->U;
	ptr->L = ptr1->L + ptr2->L;
	break;
      case DIVISION:
	ptr->U = ptr1->U + ptr2->L;
	ptr->L = ptr1->L + ptr2->U;
	break;
      case ADDITION:
      case SUBTRACTION:
	ptr->U = Max(ptr1->U + ptr2->L, ptr2->U + ptr1->L) + 1;
	ptr->L = ptr1->L + ptr2->L;
	break;
      case SQUAREROOT:
	ptr->U = ptr1->U;
	ptr->L = ptr1->L;
	if (!ptr->U.iszero())
	  ptr->U = (ptr->U >> 1) + 1;
	if (!ptr->L.iszero())
	  ptr->L = (ptr->L >> 1) + 1;
	break;
      case ROOT:
	ptr->U = ptr1->U;
	ptr->L = ptr1->L;
	if (!ptr->U.iszero())
	  ptr->U = ptr->U / d + 1;
	if (!ptr->L.iszero())
	  ptr->L = ptr->L / d + 1;
  }
}




void real_rep::compute_parameters()
{
  if (status == VISITED)
    return;
  if (OP1)
    OP1->compute_parameters();
  if (OP2)
    OP2->compute_parameters();
  if (!sep_bound_ptr)
    sep_bound_ptr = new sep_bound();
  compute_concrete_parameters();
  status = VISITED;
}



void real_rep::estimate_degree(integer & D, long &k, bool &contains_divisions)
{
  if (status == VISITED)
    return;
  bool contains_divisions_1 = false, contains_divisions_2 = false;
  if (OP1)
    OP1->estimate_degree(D, k, contains_divisions_1);
  if (OP2)
    OP2->estimate_degree(D, k, contains_divisions_2);
  contains_divisions = contains_divisions_1 || contains_divisions_2;
  if (CON == SQUAREROOT) {
    if (contains_divisions)
      k += 2;
    else
      k++;
  }
  if (CON == ROOT) {
    if (contains_divisions)
      D *= (d * d);
    else
      D *= d;
  }
  if (CON == DIVISION)
    contains_divisions = true;
  status = VISITED;
}


void real_rep::clear_visited_marks()
{
  if (status == CLEARED)
    return;
  if (OP1)
    OP1->clear_visited_marks();
  if (OP2)
    OP2->clear_visited_marks();
  status = CLEARED;
}


void real_rep::compute_concrete_system_bound(integer & system_bound)
{
  clear_visited_marks();
  compute_parameters();
  clear_visited_marks();
  long number_of_sqrt_operations = 0;
  integer degree = 1;
  bool contains_divisions = false;
  estimate_degree(degree, number_of_sqrt_operations, contains_divisions);
  degree *= power_two(number_of_sqrt_operations);
  system_bound =
    sep_bound_ptr->U * (1 - degree) - sep_bound_ptr->L;
}



int real_rep::sign_with_separation_bound(const integer & user_bound, bool supported)
{

  if (is_finite(ERROR) && (ERROR < 1 || NUM == 0))
    return float_sign(NUM);
  init_app_bf();
  if
    (exact() || (!isZero(num_bf()) && (num_exp() > err_exp())))
    return sign(num_bf());



  bool user_bound_is_reached = false;
  integer relative_precision = 26;
  integer absolute_precision;
  do {
    relative_precision = relative_precision * 2;
    absolute_precision = err_exp() - relative_precision;
    guarantee_bound_two_to(absolute_precision);
    user_bound_is_reached = supported && (1 + absolute_precision < -user_bound);
  }
  while (!(user_bound_is_reached ||
      (exact() || (!isZero(num_bf()) && (num_exp() > err_exp())))
    ));


  if
    (exact() || (!isZero(num_bf()) && (num_exp() > err_exp())))
    return sign(num_bf());

  NUM = ERROR = 0;
  num_bf() = error_bf() = bigfloat::pZero;

  {
    if (OP1)
      if (--OP1->count == 0)
	delete OP1;
    if (OP2)
      if (--OP2->count == 0)
	delete OP2;
    OP1 = OP2 = nil;
    CON = BIGFLOAT;
  }






  return 0;
}


bool operator == (const real & x, const real & y) {
  double x_error, y_error;
  bool is_easy;
  if (x.is_double && y.is_double)
    is_easy = true;
  else {
    if (x.is_double)
      x_error = 0;
    else {
      ((real &) x).value = x.PTR->NUM;
      x_error = float_abs(x.value) * x.PTR->ERROR;
    }
    if (y.is_double)
      y_error = 0;
    else {
      ((real &) y).value = y.PTR->NUM;
      y_error = float_abs(y.value) * y.PTR->ERROR;
    }
    is_easy =
      (float_abs(x.value - y.value) > (x_error + y_error + twoMinDbl) * correction);
  }


  if (is_easy)
    return (x.value == y.value);
  return ((y - x).sign() == 0);
}

bool operator != (const real & x, const real & y) {
  double x_error, y_error;
  bool is_easy;
  if (x.is_double && y.is_double)
    is_easy = true;
  else {
    if (x.is_double)
      x_error = 0;
    else {
      ((real &) x).value = x.PTR->NUM;
      x_error = float_abs(x.value) * x.PTR->ERROR;
    }
    if (y.is_double)
      y_error = 0;
    else {
      ((real &) y).value = y.PTR->NUM;
      y_error = float_abs(y.value) * y.PTR->ERROR;
    }
    is_easy =
      (float_abs(x.value - y.value) > (x_error + y_error + twoMinDbl) * correction);
  }

  ;
  if (is_easy)
    return (x.value != y.value);
  return ((y - x).sign() != 0);
}

bool operator < (const real & x, const real & y) {
  double x_error, y_error;
  bool is_easy;
  if (x.is_double && y.is_double)
    is_easy = true;
  else {
    if (x.is_double)
      x_error = 0;
    else {
      ((real &) x).value = x.PTR->NUM;
      x_error = float_abs(x.value) * x.PTR->ERROR;
    }
    if (y.is_double)
      y_error = 0;
    else {
      ((real &) y).value = y.PTR->NUM;
      y_error = float_abs(y.value) * y.PTR->ERROR;
    }
    is_easy =
      (float_abs(x.value - y.value) > (x_error + y_error + twoMinDbl) * correction);
  }

  ;
  if (is_easy)
    return (x.value < y.value);
  return ((y - x).sign() > 0);
}

bool operator > (const real & x, const real & y) {
  double x_error, y_error;
  bool is_easy;
  if (x.is_double && y.is_double)
    is_easy = true;
  else {
    if (x.is_double)
      x_error = 0;
    else {
      ((real &) x).value = x.PTR->NUM;
      x_error = float_abs(x.value) * x.PTR->ERROR;
    }
    if (y.is_double)
      y_error = 0;
    else {
      ((real &) y).value = y.PTR->NUM;
      y_error = float_abs(y.value) * y.PTR->ERROR;
    }
    is_easy =
      (float_abs(x.value - y.value) > (x_error + y_error + twoMinDbl) * correction);
  }

  ;
  if (is_easy)
    return (x.value > y.value);
  return ((y - x).sign() < 0);
}

bool operator <= (const real & x, const real & y) {
  double x_error, y_error;
  bool is_easy;
  if (x.is_double && y.is_double)
    is_easy = true;
  else {
    if (x.is_double)
      x_error = 0;
    else {
      ((real &) x).value = x.PTR->NUM;
      x_error = float_abs(x.value) * x.PTR->ERROR;
    }
    if (y.is_double)
      y_error = 0;
    else {
      ((real &) y).value = y.PTR->NUM;
      y_error = float_abs(y.value) * y.PTR->ERROR;
    }
    is_easy =
      (float_abs(x.value - y.value) > (x_error + y_error + twoMinDbl) * correction);
  }

  ;
  if (is_easy)
    return (x.value <= y.value);
  return ((y - x).sign() >= 0);
}

bool operator >= (const real & x, const real & y) {
  double x_error, y_error;
  bool is_easy;
  if (x.is_double && y.is_double)
    is_easy = true;
  else {
    if (x.is_double)
      x_error = 0;
    else {
      ((real &) x).value = x.PTR->NUM;
      x_error = float_abs(x.value) * x.PTR->ERROR;
    }
    if (y.is_double)
      y_error = 0;
    else {
      ((real &) y).value = y.PTR->NUM;
      y_error = float_abs(y.value) * y.PTR->ERROR;
    }
    is_easy =
      (float_abs(x.value - y.value) > (x_error + y_error + twoMinDbl) * correction);
  }

  ;
  if (is_easy)
    return (x.value >= y.value);
  return ((y - x).sign() <= 0);
}



void real::improve_approximation_to(const integer & p) const
{
  if (is_double)
    return;
  PTR->init_app_bf();
  PTR->guarantee_bound_two_to(-p);
  adjust();
}


void real_rep::guarantee_bound_two_to(const integer & e)
{
  bool current_precision_is_sufficient =
  isZero(error_bf()) || (err_exp() <= e);
  if (current_precision_is_sufficient)
    return;

  switch (CON) {
      case DOUBLE:
      case BIGFLOAT:
	return;
      case NEGATION:
	OP1->guarantee_bound_two_to(e);
	num_bf() = -OP1->num_bf();
	error_bf() = OP1->error_bf();
	break;
      case ADDITION:
	improve_add(e);
	break;
      case SUBTRACTION:
	improve_sub(e);
	break;
      case MULTIPLICATION:
	improve_mul(e);
	break;
      case DIVISION:
	improve_div(e);
	break;
      case SQUAREROOT:
	improve_sqrt(e);
	break;
      case ROOT:
	improve_root(e);
	break;
  }

  if (exact()) {
    if (OP1)
      if (--OP1->count == 0)
	delete OP1;
    if (OP2)
      if (--OP2->count == 0)
	delete OP2;
    OP1 = OP2 = nil;
    CON = BIGFLOAT;
  }





  integer system_bound;
  if (!
    (exact() || (!isZero(num_bf()) && (num_exp() > err_exp())))
    ) {
    compute_concrete_system_bound(system_bound);
    bool system_bound_is_reached = (e + 1 < system_bound);
    if (system_bound_is_reached) {

      NUM = ERROR = 0;
      num_bf() = error_bf() = bigfloat::pZero;

      {
	if (OP1)
	  if (--OP1->count == 0)
	    delete OP1;
	if (OP2)
	  if (--OP2->count == 0)
	    delete OP2;
	OP1 = OP2 = nil;
	CON = BIGFLOAT;
      }






      sep_bound_ptr->U = sep_bound_ptr->L = zero_integer;
    }
  }

  adjust_dbl();
}


void real_rep::compute_op(long p)
{
  bool isexact;
  bigfloat operation_error, operand_error;
  switch (CON) {
      case DOUBLE:
      case BIGFLOAT:
	return;
      case NEGATION:
	num_bf() = -OP1->num_bf();
	error_bf() = OP1->error_bf();
	break;
      case ADDITION:
	num_bf() = add(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);

	operand_error = add(OP1->error_bf(), OP2->error_bf(), 16, TO_INF);



	break;
      case SUBTRACTION:
	num_bf() = sub(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);

	operand_error = add(OP1->error_bf(), OP2->error_bf(), 16, TO_INF);



	break;
      case MULTIPLICATION:
	num_bf() = mul(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);
	compute_mul_error(operand_error);
	break;
      case DIVISION:
	if (OP2->sign_with_separation_bound(zero_integer, false) == 0)
	  error_handler(1, "compute_op:division by zero");
	num_bf() = div(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);
	compute_div_error(operand_error);
	break;
      case SQUAREROOT:
	if (OP1->sign_with_separation_bound(zero_integer, false) < 0)
	  error_handler(1, "compute_op:sqrt of negative number");
	num_bf() = sqrt(OP1->num_bf(), p, TO_NEAREST, isexact);
	compute_sqrt_error(operand_error);
	break;
      case ROOT:
	if (OP1->sign_with_separation_bound(zero_integer, false) < 0)
	  error_handler(1, "compute_op:root of negative number");
	num_bf() = sqrt_d(OP1->num_bf(), p, d);
	isexact = false;
	compute_root_error(operand_error);
	break;
  }
  if (!isexact)
    operation_error = ldexp_bf(num_bf(), -p);
  error_bf() = operation_error + operand_error;
  if (exact()) {
    if (OP1)
      if (--OP1->count == 0)
	delete OP1;
    if (OP2)
      if (--OP2->count == 0)
	delete OP2;
    OP1 = OP2 = nil;
    CON = BIGFLOAT;
  }




  adjust_dbl();
}


void real_rep::improve_add(const integer & e)
{
  long p;

  OP1->guarantee_bound_two_to(e - 2);
  OP2->guarantee_bound_two_to(e - 2);
  integer log_epsilon = -1;
  if (!isZero(OP1->num_bf()) && !isZero(OP2->num_bf()))
    log_epsilon = e - Maximum(OP1->num_exp(), OP2->num_exp()) - 2;

  if (!log_epsilon.islong()) {
    error_handler(1, "sorry:improve:computation to expensive");
    p = 10000;			/* just to silence compiler warnings */
  }
  else
    p = -log_epsilon.tolong();
  if (p < 2)
    p = 2;




  bool isexact;
  num_bf() = add(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);

  if (isexact && OP1->exact() && OP2->exact())
    error_bf() = bigfloat::pZero;
  else
    error_bf() = pow2(e);



}

void real_rep::improve_sub(const integer & e)
{
  long p;

  OP1->guarantee_bound_two_to(e - 2);
  OP2->guarantee_bound_two_to(e - 2);
  integer log_epsilon = -1;
  if (!isZero(OP1->num_bf()) && !isZero(OP2->num_bf()))
    log_epsilon = e - Maximum(OP1->num_exp(), OP2->num_exp()) - 2;

  if (!log_epsilon.islong()) {
    error_handler(1, "sorry:improve:computation to expensive");
    p = 10000;			/* just to silence compiler warnings */
  }
  else
    p = -log_epsilon.tolong();
  if (p < 2)
    p = 2;




  bool isexact;
  num_bf() = sub(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);

  if (isexact && OP1->exact() && OP2->exact())
    error_bf() = bigfloat::pZero;
  else
    error_bf() = pow2(e);



}


void real_rep::compute_mul_error(bigfloat & operand_error)
{
  bigfloat x = OP1->num_bf();
  bigfloat y = OP2->num_bf();
  bigfloat error_x = mul(abs(round(x, 16, TO_INF)), OP2->error_bf(), 32, TO_INF);
  bigfloat error_y = mul(abs(round(y, 16, TO_INF)), OP1->error_bf(), 32, TO_INF);
  operand_error = add(error_x, error_y, 16, TO_INF);
}


void real_rep::improve_mul(const integer & e)
{
  bigfloat y_high = add(abs(OP2->num_bf()), OP2->error_bf(), 32, TO_P_INF);
  integer ex = e - ilog2(y_high) - 2;
  OP1->guarantee_bound_two_to(ex);
  integer log_epsilon = -1;
  if (!isZero(OP1->num_bf())) {
    integer ey = e - OP1->num_exp() - 2;
    OP2->guarantee_bound_two_to(ey);
    log_epsilon = e - (OP1->num_exp() + OP2->num_exp()) - 1;
  }
  bool isexact;
  long p;

  if (!log_epsilon.islong()) {
    error_handler(1, "sorry:improve:computation to expensive");
    p = 10000;			/* just to silence compiler warnings */
  }
  else
    p = -log_epsilon.tolong();
  if (p < 2)
    p = 2;


  num_bf() = mul(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);

  if (isexact && OP1->exact() && OP2->exact())
    error_bf() = bigfloat::pZero;
  else
    error_bf() = pow2(e);



}



void real_rep::compute_div_error(bigfloat & operand_error)
{
  bigfloat y_low = sub(abs(OP2->num_bf()), OP2->error_bf(), 32, TO_N_INF);
  bigfloat x = OP1->num_bf();
  bigfloat z = mul(abs(round(x, 16)), OP2->error_bf(), 32, TO_INF);
  z = ldexp_bf(z, 1 - OP2->num_exp());
  operand_error = add(z, abs(OP1->error_bf()), 16, TO_INF);
  operand_error = ldexp_bf(operand_error, 1 - ilog2(y_low));
}


void real_rep::improve_div(const integer & e)
{
  bigfloat y_low = sub(abs(OP2->num_bf()), OP2->error_bf(), 32, TO_N_INF);
  integer ex = e + ilog2(y_low) - 3;
  OP1->guarantee_bound_two_to(ex);
  integer log_epsilon = -1;
  if (!isZero(OP1->num_bf())) {
    integer ey = ex + OP2->num_exp() - OP1->num_exp() - 3;
    OP2->guarantee_bound_two_to(ey);
    log_epsilon = e + OP2->num_exp() - OP1->num_exp() - 2;
  }
  bool isexact;
  long p;

  if (!log_epsilon.islong()) {
    error_handler(1, "sorry:improve:computation to expensive");
    p = 10000;			/* just to silence compiler warnings */
  }
  else
    p = -log_epsilon.tolong();
  if (p < 2)
    p = 2;


  num_bf() = div(OP1->num_bf(), OP2->num_bf(), p, TO_NEAREST, isexact);

  if (isexact && OP1->exact() && OP2->exact())
    error_bf() = bigfloat::pZero;
  else
    error_bf() = pow2(e);



}


void real_rep::compute_sqrt_error(bigfloat & operand_error)
{
  integer log_sqrt = (num_exp() + 1) >> 1;
  if (!OP1->exact())
    operand_error = ldexp_bf(OP1->error_bf(), -log_sqrt);
}



void real_rep::improve_sqrt(const integer & e)
{
  bool isexact;
  bigfloat x_low = sub(OP1->num_bf(), OP1->error_bf(), 32, TO_N_INF);
  integer ex = e + (ilog2(x_low) - 1) / 2 - 1;
  OP1->guarantee_bound_two_to(ex);
  bigfloat x_high = add(OP1->num_bf(), OP1->error_bf(), 32, TO_P_INF);
  integer log_epsilon = e - (ilog2(x_high) + 1) / 2 - 1;
  long p;

  if (!log_epsilon.islong()) {
    error_handler(1, "sorry:improve:computation to expensive");
    p = 10000;			/* just to silence compiler warnings */
  }
  else
    p = -log_epsilon.tolong();
  if (p < 2)
    p = 2;


  num_bf() = sqrt(OP1->num_bf(), p, TO_NEAREST, isexact, num_bf());
  if (isexact && OP1->exact())
    error_bf() = bigfloat::pZero;
  else
    error_bf() = pow2(e);
}


void real_rep::compute_root_error(bigfloat & operand_error)
{
  integer log_sqrt = (num_exp() * (1 - d)) / d + 1;
  operand_error = ldexp_bf(OP1->error_bf(), log_sqrt - ilog2(d) + 1);
}



void real_rep::improve_root(const integer & e)
{
  bigfloat x_low = sub(OP1->num_bf(), OP1->error_bf(), 32, TO_N_INF);
  integer ex = e + ilog2(d) + ((ilog2(x_low) - 1) * (d - 1)) / d - 2;
  OP1->guarantee_bound_two_to(ex);
  bigfloat x_high = add(OP1->num_bf(), OP1->error_bf(), 32, TO_P_INF);
  integer log_epsilon = e - (ilog2(x_high) - 1) / d - 2;
  long p;

  if (!log_epsilon.islong()) {
    error_handler(1, "sorry:improve:computation to expensive");
    p = 10000;			/* just to silence compiler warnings */
  }
  else
    p = -log_epsilon.tolong();
  if (p < 2)
    p = 2;


  num_bf() = sqrt_d(OP1->num_bf(), p, d);
  error_bf() = pow2(e);
}




void real::guarantee_relative_error(long p) const
{
  if (is_double)
    return;
  if (PTR->sign_with_separation_bound(zero_integer) == 0)
    return;
  PTR->init_app_bf();
  bigfloat num_low = abs(PTR->num_bf()) - PTR->error_bf();
  PTR->guarantee_bound_two_to(ilog2(num_low) - p);
  adjust();
}


void real::compute_with_precision(long p) const
{
  if (is_double)
    return;
  PTR->init_app_bf();
  PTR->clear_visited_marks();
  PTR->compute_approximation(p);
  adjust();
}


void real_rep::compute_approximation(long p)
{
  if (status == VISITED)
    return;
  if (OP1)
    OP1->compute_approximation(p);
  if (OP2)
    OP2->compute_approximation(p);
  compute_op(p);
  status = VISITED;
}



int sign(const real & x)
{
  return x.sign();
}
real abs(const real & x)
{
  if (x.sign() < 0)
    return -x;
  else
    return x;
}
real sqr(const real & x)
{
  return x * x;
}
real dist(const real & x, const real & y)
{
  return sqrt(x * x + y * y);
}
real powi(const real & x, int n)
{
  real y = x, z = 1;
  int n_prefix = n;
  while (n_prefix > 0) {
    if (n_prefix % 2)
      z = z * y;
    n_prefix = n_prefix / 2;
    y = y * y;
  }
  return z;
}
