/****************************************************************************
 * Core Library Version 1.7, August 2004
 * Copyright (c) 1995-2004 Exact Computation Project
 * All rights reserved.
 *
 * This file is part of CGAL (www.cgal.org).
 *
 * File: BigFloat.cpp
 * Synopsis:
 *       BigFloat numbers with error bounds
 *
 *       EXACTNESS PROPERTY:
 *       ==================
 *       For BigFloats that are exact (i.e., error=0),
 *       addition/subtraction and multiplication return the
 *       exact result (i.e., error=0).  We also introduce the operation
 *       div2(), which simply divides a BigFloat by 2,
 *       but this again preserves exactness.  Such exactness
 *       properties are used in our Newton iteration/Sturm Sequences.
 *
 * Written by
 *       Chee Yap <yap@cs.nyu.edu>
 *       Chen Li <chenli@cs.nyu.edu>
 *       Zilin Du <zilin@cs.nyu.edu>
 *
 * WWW URL: https://cs.nyu.edu/exact/
 * Email: exact@cs.nyu.edu
 *
 * $URL: https://github.com/CGAL/cgal/blob/v6.1.1/CGAL_Core/include/CGAL/CORE/BigFloat_impl.h $
 * $Id: include/CGAL/CORE/BigFloat_impl.h 08b27d3db14 $
 * SPDX-License-Identifier: LGPL-3.0-or-later
 ***************************************************************************/

#ifdef CGAL_HEADER_ONLY
#define CGAL_INLINE_FUNCTION inline
#else
#define CGAL_INLINE_FUNCTION
#endif

#include <CGAL/disable_warnings.h>

#include <string>
#include <ctype.h>
#include <CGAL/CORE/BigFloat.h>
#include <CGAL/CORE/Expr.h>
#include <CGAL/tss.h>
#include <sstream>

namespace CORE {


////////////////////////////////////////////////////////////
// Misc Helper Functions
////////////////////////////////////////////////////////////

CGAL_INLINE_FUNCTION
BigInt FiveTo(unsigned long exp) {
  if (exp == 0)
    return BigInt(1);
  else if (exp == 1)
    return BigInt(5);
  else {
    BigInt x = FiveTo(exp / 2);

    x = x * x;

    if (exp & 1)
      x *= 5;

    return x;
  }
}

////////////////////////////////////////////////////////////
//  class BigFloat
////////////////////////////////////////////////////////////

// STATIC BIGFLOAT CONSTANTS
// ZERO
CGAL_INLINE_FUNCTION
const BigFloat& BigFloat::getZero() {
  init_CORE();
  CGAL_STATIC_THREAD_LOCAL_VARIABLE(BigFloat, Zero,0);
  return Zero;
}
// ONE
CGAL_INLINE_FUNCTION
const BigFloat& BigFloat::getOne() {
  init_CORE();
  CGAL_STATIC_THREAD_LOCAL_VARIABLE(BigFloat, One,1);
  return One;
}

// A special constructor for BigFloat from Expr
// -- this method is somewhat of an anomaly (we normally do not expect
//    BigFloats to know about Expr).
CGAL_INLINE_FUNCTION
BigFloat::BigFloat(const Expr& E, const extLong& r, const extLong& a)
  : RCBigFloat(new BigFloatRep()) {
  *this = E.approx(r, a).BigFloatValue(); // lazy implementation, any other way?
}

////////////////////////////////////////////////////////////
//  class BigFloatRep
////////////////////////////////////////////////////////////

CGAL_INLINE_FUNCTION
BigFloatRep::BigFloatRep(double d) : m(0), err(0), exp(0) {
  if (d != 0.0) {
    int isNegative = 0;

    if (d < 0.0) {
      isNegative = 1;
      d          = - d;
    }

    int    binExp;
    double f = frexp(d, &binExp);

    exp = chunkFloor(binExp);

    long s = binExp - bits(exp);

    long   stop = 0;
    double intPart;

    // convert f into a BigInt
    while (f != 0.0 && stop < DBL_MAX_CHUNK) {
      f =   ldexp(f, (int)CHUNK_BIT);
      f =   modf(f, &intPart);
      m <<= CHUNK_BIT;
      m +=  (long)intPart;
      exp--;
      stop++;
    }
#ifdef CGAL_CORE_DEBUG
    CGAL_assertion (s >= 0);
#endif

    if (s)
      m <<= s;
    if (isNegative)
      negate(m);
  }
}//BigFloatRep constructor

//  approximation
CGAL_INLINE_FUNCTION
void BigFloatRep::trunc(const BigInt& I, const extLong& r, const extLong& a) {
  if (sign(I)) {
    long tr = chunkFloor((- r + bitLength(I)).asLong());
    long ta = chunkFloor(- a.asLong());
    long t;

    if (r.isInfty() || a.isTiny())
      t = ta;
    else if (a.isInfty())
      t = tr;
    else
      t = ta < tr ? tr : ta;

    if (t > 0) {  // BigInt remainder;
      m   = chunkShift(I, - t);
      err = 1;
      exp = t;
    } else { //  t <= 0
      m   = I;
      err = 0;
      exp = 0;
    }
  } else {//  I == 0
    m   = 0;
    err = 0;
    exp = 0;
  }
}

CGAL_INLINE_FUNCTION
void BigFloatRep :: truncM(const BigFloatRep& B, const extLong& r, const extLong& a) {
  if (sign(B.m)) {
    long tr = chunkFloor((- 1 - r + bitLength(B.m)).asLong());
    long ta = chunkFloor(- 1 - a.asLong()) - B.exp;
    long t;

    if (r.isInfty() || a.isTiny())
      t = ta;
    else if (a.isInfty())
      t = tr;
    else
      t = ta < tr ? tr : ta;

    if (t >= chunkCeil(clLg(B.err))) {
      m   = chunkShift(B.m, - t);
      err = 2;
      exp = B.exp + t;
    } else //  t < chunkCeil(clLg(B.err))
      CGAL_error_msg("BigFloat error: truncM called with stricter precision than current error.");
  } else {//  B.m == 0
    long t = chunkFloor(- a.asLong()) - B.exp;

    if (t >= chunkCeil(clLg(B.err))) {
      m   = 0;
      err = 1;
      exp = B.exp + t;
    } else //  t < chunkCeil(clLg(B.err))
      CGAL_error_msg("BigFloat error: truncM called with stricter precision than current error.");
  }
}

// This is the main approximation function
// REMARK: would be useful to have a self-modifying version
//                 of this function (e.g., for Newton).
CGAL_INLINE_FUNCTION
void BigFloatRep::approx(const BigFloatRep& B,
              const extLong& r, const extLong& a) {
  if (B.err) {
    if (static_cast<std::size_t>(1 + clLg(B.err)) <= bitLength(B.m))
      truncM(B, r + 1, a);
    else //  1 + clLg(B.err) > lg(B.m)
      truncM(B, CORE_posInfty, a);
  } else {//  B.err == 0
    trunc(B.m, r, a - bits(B.exp));
    exp += B.exp;
  }
  // Call normalization globally     -- IP 10/9/98
  normal();
}

CGAL_INLINE_FUNCTION
void BigFloatRep::div(const BigInt& N, const BigInt& D,
              const extLong& r, const extLong& a) {
  if (sign(D)) {
    if (sign(N)) {
      long tr = chunkFloor((- r + bitLength(N) - bitLength(D) - 1).asLong());
      long ta = chunkFloor(- a.asLong());

      if (r.isInfty() || a.isTiny())
        exp = ta;
      else if (a.isInfty())
        exp = tr;
      else
        exp = ta < tr ? tr : ta;

      BigInt remainder;

      // divide(chunkShift(N, - exp), D, m, remainder);
      div_rem(m, remainder, chunkShift(N, - exp), D);

      if (exp <= 0 && sign(remainder) == 0)
        err = 0;
      else
        err = 1;
    } else {//  N == 0
      m   = 0;
      err = 0;
      exp = 0;
    }
  } else //  D == 0
    CGAL_error_msg( "BigFloat error: zero divisor.");

  // Call normalization globally     -- IP 10/9/98
  normal();
}//div

//  error-normalization
CGAL_INLINE_FUNCTION
void BigFloatRep::normal() {
  long le = flrLg(err);

  if (le >= CHUNK_BIT + 2) { // so we do not carry more than 16 = CHUNK_BIT + 2
                             // bits of error
    long f = chunkFloor(--le); // f is roughly equal to floor(le/CHUNK_BIT)
    long bits_f = bits(f);   // f chunks will have bits_f many bits
#ifdef CGAL_CORE_DEBUG
    CGAL_assertion (bits_f >= 0);
#endif

    m   >>= bits_f;  // reduce mantissa by bits_f many bits
    err >>= bits_f;  // same for err
    err +=  2;       // why 2?
    exp +=  f;
  }
  if (err == 0)      // unlikely, if err += 2 above
    eliminateTrailingZeroes();
}

// bigNormal(err)
//     convert a bigInt error value (=err) into an error that fits into
//     a long number.  This is done by
//     by increasing the exponent, and corresponding decrease
//     in the bit lengths of the mantissa and error.
//
CGAL_INLINE_FUNCTION
void BigFloatRep::bigNormal(BigInt& bigErr) {
  long le = bitLength(bigErr);

  if (le < CHUNK_BIT + 2) {
    err = ulongValue(bigErr);
  } else {
    long f = chunkFloor(--le);
    long bits_f = bits(f);
#ifdef CGAL_CORE_DEBUG
    CGAL_assertion(bits_f >= 0);
#endif

    m      >>= bits_f;
    bigErr >>= bits_f;
    err    = ulongValue(bigErr) + 2; // you need to add "2" because "1" comes
                    // from truncation error in the mantissa, and another
                // "1" comes from the truncation error in the bigErr.
                // (But there is danger of overflow...)
    exp    += f;
  }

  if (err == 0)
    eliminateTrailingZeroes();
}

// ARITHMETIC:
//  Addition
CGAL_INLINE_FUNCTION
void BigFloatRep::add(const BigFloatRep& x, const BigFloatRep& y) {
  long expDiff = x.exp - y.exp;

  if (expDiff > 0) {//  x.exp > y.exp
    if (!x.err) {
      m   = chunkShift(x.m, expDiff) + y.m;
      err = y.err;
      exp = y.exp;
    } else {//  x.err > 0
      m   = x.m + chunkShift(y.m, - expDiff); // negative shift!
      err = x.err + 5; // To account for y.err (but why 5?)
      exp = x.exp;     //
      // normal();
    }
  } else if (!expDiff) {//  x.exp == y.exp
    m   = x.m + y.m;
    err = x.err + y.err;
    exp = x.exp;
    // normal();
  } else {//  x.exp < y.exp
    if (!y.err) {
      m   = x.m + chunkShift(y.m, - expDiff);
      err = x.err;
      exp = x.exp;
    } else {//  y.err > 0
      m   = chunkShift(x.m, expDiff) + y.m;
      err = y.err + 5;
      exp = y.exp;
      // normal();
    }
  }
  // Call normalization globally     -- IP 10/9/98
  normal();
}

//  Subtraction
CGAL_INLINE_FUNCTION
void BigFloatRep::sub(const BigFloatRep& x, const BigFloatRep& y) {
  long expDiff = x.exp - y.exp;

  if (expDiff > 0) {//  x.exp > y.exp
    if (!x.err) {
      m   = chunkShift(x.m, expDiff) - y.m;
      err = y.err;
      exp = y.exp;
    } else {//  x.err > 0
      m   = x.m - chunkShift(y.m, - expDiff);
      err = x.err + 5;
      exp = x.exp;
      // normal();
    }
  } else if (!expDiff) {
    m   = x.m - y.m;
    err = x.err + y.err;
    exp = x.exp;
    // normal();
  } else { //  x.exp < y.exp
    if (!y.err) {
      m   = x.m - chunkShift(y.m, - expDiff);
      err = x.err;
      exp = x.exp;
    } else {//  y.err > 0
      m   = chunkShift(x.m, expDiff) - y.m;
      err = y.err + 5;
      exp = y.exp;
      // normal();
    }
  }
  // Call normalization globally     -- IP 10/9/98
  normal();
}

CGAL_INLINE_FUNCTION
void BigFloatRep::mul(const BigFloatRep& x, const BigFloatRep& y) {
  m = x.m * y.m;
  exp = x.exp + y.exp;
  // compute error (new code, much faster. Zilin Du, Nov 2003)
  if (x.err == 0 && y.err == 0) {
    err = 0;
    eliminateTrailingZeroes();
  } else {
    BigInt bigErr(0);
    if (y.err != 0)
      bigErr += abs(x.m)*y.err;
    if (x.err != 0)
      bigErr += abs(y.m)*x.err;
    if (x.err !=0 && y.err != 0)
      bigErr += x.err*y.err;
    bigNormal(bigErr);
  }
}
// BigFloat div2 will half the value of x, exactly with NO error
//         REMARK: should generalize this to dividing by any power of 2
//         We need this in our use of BigFloats to maintain isolation
//         intervals (e.g., in Sturm sequences)        --Chee/Vikram 4/2003
//
CGAL_INLINE_FUNCTION
void BigFloatRep :: div2(const BigFloatRep& x) {
  if (isEven(x.m)) {
    m = (x.m >> 1);
    exp = x.exp ;
  } else {
    m = (x.m << static_cast<unsigned long>(CHUNK_BIT-1));
    exp = x.exp -1;
  }
}

// Converts a BigFloat interval into one BigFloat with almost same error bound
// This routine ignores the errors in inputs a and b.
// But you cannot really ignore them since, they are taken into account
// when you compute "r.sub(a,b)"...
CGAL_INLINE_FUNCTION
void BigFloatRep::centerize(const BigFloatRep& a, const BigFloatRep& b) {
  if ((a.m == b.m) && (a.err == b.err) && (a.exp == b.exp)) {
    m = a.m;
    err = a.err;
    exp = a.exp;
    return;
  }

  BigFloatRep r;
  r.sub(a, b);
  r.div2(r);

  //setup mantissa and exponent, but not error bits
  // But this already sets the error bits? Chee
  add(a,b);
  div2(*this);
  // error bits = ceil ( B^{-exp}*|a-b|/2 )

  // bug fixed: possible overflow on converting
  // Zilin & Vikram, 08/24/04
  // err = 1 + longValue(chunkShift(r.m, r.exp - exp));
  BigInt E = chunkShift(r.m, r.exp - exp);
  E = abs(E);
  bigNormal(E);
}

// BigFloat Division, computing x/y:
//      Unlike +,-,*, this one takes a relative precision bound R
//        Note that R is only used when x and y are error-free!
//        (This remark means that we may be less efficient than we could be)
//
//          Assert( R>0  && R< CORE_Infty )
//
CGAL_INLINE_FUNCTION
void BigFloatRep :: div(const BigFloatRep& x, const BigFloatRep& y,
                        const extLong& R) {
  if (!y.isZeroIn()) { //  y.m > y.err, so we are not dividing by 0
    if (!x.err && !y.err) {
      if (R < 0 || R.isInfty()) //Oct 9, 2002: fixed major bug! [Zilin/Chee]
        div(x.m, y.m, get_static_defBFdivRelPrec(), CORE_posInfty);
      else
        div(x.m, y.m, R, CORE_posInfty);
      exp += x.exp - y.exp; // chen: adjust exp.
    } else {//  x.err > 0 or y.err > 0
      BigInt bigErr, errRemainder;

      if (x.isZeroIn()) { //  x.m <= x.err
        m   = 0;
        exp = x.exp - y.exp;

        div_rem(bigErr, errRemainder, abs(x.m) + static_cast<long>(x.err),
                abs(y.m) - static_cast<long>(y.err));
      } else { //  x.m > x.err
        long lx = bitLength(x.m);
        long ly = bitLength(y.m);
        long r;

        if (!x.err) //  x.err == 0 and y.err > 0
          r = ly + 2;
        else if(!y.err)  //  x.err > 0 and y.err == 0
          r = lx + 2;
        else  //  x.err > 0 and y.err > 0
          r = lx < ly ? lx + 2: ly + 2;

        long   t = chunkFloor(- r + lx - ly - 1);
        BigInt remainder;

        div_rem(m, remainder, chunkShift(x.m, - t), y.m);
        exp = t + x.exp - y.exp;

        long delta = ((t > 0) ? 2 : 0);

        // Chen Li: 9/9/99
        // here again, it use ">>" operator with a negative
        // right operand. So the result is not well defined.
        // Erroneous code:
        //   divide(abs(remainder) + (static_cast<long>(x.err) >> bits(t))
        //                     + delta + static_cast<long>(y.err) * abs(m),
        //                     abs(y.m) - static_cast<long>(y.err),
        //                     bigErr,
        //                     errRemainder);
        // New code:
        BigInt errx_over_Bexp = x.err;
        long bits_Bexp = bits(t);
        if (bits_Bexp >= 0) {
          errx_over_Bexp >>= bits_Bexp;
        } else {
          errx_over_Bexp <<= (-bits_Bexp);
        }

        // divide(abs(remainder) + errx_over_Bexp
        //        + delta + static_cast<long>(y.err) * abs(m),
        //        abs(y.m) - static_cast<long>(y.err),
        //        bigErr,
        //        errRemainder);
        div_rem(bigErr, errRemainder,
                abs(remainder) + errx_over_Bexp + delta + static_cast<long>(y.err) * abs(m),
                abs(y.m) - static_cast<long>(y.err));
      }

      if (sign(errRemainder))
        ++bigErr;

      bigNormal(bigErr);
    }
  } else {//  y.m <= y.err
    CGAL_error_msg("BigFloat error: possible zero divisor.");
  }

  // Call normalization globally     -- IP 10/9/98
  // normal(); -- Chen: after calling bigNormal, this call is redundant.
}// BigFloatRep::div

//  squareroot for BigInt argument, without initial approximation
//  sqrt(x,a) computes sqrt of x to absolute precision a.
//      -- this is where Newton is applied
//      -- this is called by BigFloatRep::sqrt(BigFloat, extLong)
CGAL_INLINE_FUNCTION
void BigFloatRep::sqrt(const BigInt& x, const extLong& a) {
  sqrt(x, a, BigFloat(x, 0, 0));
} // sqrt(BigInt x, extLong a) , without initial approx

//  sqrt(x,a,A) where
//      x = bigInt whose sqrt is to be computed
//      a = absolute precision bound
//      A = initial approximation in BigFloat
//  -- this is where Newton is applied
//  -- it is called by BigFloatRep::sqrt(BigFloatRep, extLong, BigFloat)
CGAL_INLINE_FUNCTION
void BigFloatRep::sqrt(const BigInt& x, const extLong& a, const BigFloat& A) {
  if (sign(x) == 0) {
    m = 0;
    err = 0;
    exp = 0;
  } else if (x == 1) {
    m = 1;
    err = 0;
    exp = 0;
  } else  {// main case
    // here is where we use the initial approximation
    m = A.m();
    err = 0;
    exp = A.exp();

    BigFloatRep q, z;
    extLong     aa;
    // need this to make sure that in case the
    // initial approximation A is less than sqrt(x)
    // then Newton iteration will still proceed at
    // least one step.
    bool firstTime = true;
    for (;;) {
      aa    = a - bits(exp);
      q.div(x, m, CORE_posInfty, aa);
      q.err = 0;
      q.exp -= exp;

      z.sub(*this, q);  // this=current approximation, so z = this - q
      /*if (sign(z.m) <= 0 || z.MSB() < - a)  // justification: see Koji's
          break;                              // thesis (p. 28) which states
                                              // that we can exit when
                                              // " (*this) <= q + 2**(-a)"
      */
      // The preceding code is replaced by what follows:
      if (z.MSB() < -a)
        break;
      if (sign(z.m) <= 0) {
        if (firstTime)
          firstTime = false;
        else
          break;
      }

      z.add(*this, q);
      // Chen Li: a bug fixed here.
      //      m   = z.m >> 1;
      //      err = 0;
      //      exp = z.exp;
      if ((z.m > 1) && isEven(z.m)) {
        m = z.m >> 1;               // exact division by 2
        err = 0;
        exp = z.exp;
      } else {                      // need to shift left before division by 2
        m = chunkShift(z.m, 1) >> 1;
        err = 0;
        exp = z.exp - 1;
      }//else
    }//for
  }//else
} // sqrt of BigInt, with initial approx

// MAIN ENTRY INTO SQRT FUNCTION (BIGFLOAT ARGUMENT, WITHOUT INITIAL APPROX)
CGAL_INLINE_FUNCTION
void BigFloatRep::sqrt(const BigFloatRep& x, const extLong& a) {
  sqrt(x, a, BigFloat(x.m, 0, x.exp));
} //sqrt(BigFloat, extLong a)

// MAIN ENTRY INTO SQRT FUNCTION (BIGFLOAT ARGUMENT WITH INITIAL APPROXIMATION)
CGAL_INLINE_FUNCTION
void BigFloatRep::sqrt(const BigFloatRep& x, const extLong& a, const BigFloat& A) {
  // This computes the sqrt of x to absolute precision a, starting with
  // the initial approximation A
  if (sign(x.m) >= 0) {          //  x.m >= 0
    int delta = x.exp & 1;    // delta=0 if x.exp is even, otherwise delta=1

    if (x.isZeroIn()) {         //  x.m <= x.err
      m = 0;
      if (!x.err)
        err = 0;
      else {                  //  x.err > 0
        err = (long)(std::sqrt((double)x.err));
        err++;
        err <<= 1;
        if (delta)
          err <<= HALF_CHUNK_BIT;
      }
      exp = x.exp >> 1;
      normal();
    } else {
      long aExp = A.exp() - (x.exp >> 1);
      BigFloat AA( chunkShift(A.m(), delta), 0, aExp);

      if (!x.err) {             //  x.m > x.err = 0 (ERROR FREE CASE)
        BigFloatRep z;
        extLong ppp;
        if (a.isInfty())        //Oct 9, 2002: fixed major bug! [Zilin/Chee]
          ppp = get_static_defBFsqrtAbsPrec();
        else
          ppp = a + EXTLONG_ONE;
        extLong absp  = ppp + bits(x.exp >> 1);

        z.sqrt(chunkShift(x.m, delta), absp, AA); // call sqrt(BigInt, a, AA)

        long p = (absp + bits(z.exp)).asLong();

        // Next, normalize the error:
        if (p <= 0) {
          m = z.m;
          // Chen Li: a bug fixed
          //    BigInt bigErr = 1 << (-p);
          BigInt bigErr(1);
          bigErr = bigErr << static_cast<unsigned long>(-p);
          exp = z.exp + (x.exp >> 1);
          bigNormal(bigErr);
        } else {                  //  p > 0
          m = chunkShift(z.m, chunkCeil(p));
          long r = CHUNK_BIT - 1 - (p + CHUNK_BIT - 1) % CHUNK_BIT;
#ifdef CGAL_CORE_DEBUG
          CGAL_assertion(r >= 0);
#endif

          err = 1 >> r;
          exp = - chunkCeil(ppp.asLong());
          normal();
        }
      } else {                      //  x.m > x.err > 0 (mantissa has error)
        BigFloatRep z;
        extLong absp=-flrLg(x.err)+bitLength(x.m)-(bits(delta) >> 1)+EXTLONG_FOUR;

        z.sqrt(chunkShift(x.m, delta), absp, AA);

        long qqq = - 1 + (bitLength(x.m) >> 1) - delta * HALF_CHUNK_BIT;
        long qq  = qqq - clLg(x.err);
        long q   = qq + bits(z.exp);

        if (q <= 0) {
          m = z.m;
          long   qqqq   = - qqq - bits(z.exp);
          // Chen Li (09/08/99), a bug fixed here:
          //        BigInt bigErr = x.err << - qqqq;
          // when (-qqqq) is negative, the result is not correct.
          // how "<<" and ">>" process negative second operand is
          // not well defined. Seems it just take it as a unsigned
          // integer and extract the last few bits.
          // x.err is a long number which easily overflows.
          // From page 22 of Koji's paper, I think the exponent is
          // wrong here. So I rewrote it as:
          BigInt bigErr = x.err;
          if (qqqq >= 0) {
            bigErr <<= qqqq;
          } else {
            bigErr >>= (-qqqq);
            ++bigErr; // we need to keep its ceiling.
          }

          exp = z.exp + (x.exp >> 1);
          bigNormal(bigErr);
        } else {         //  q > 0
          m = chunkShift(z.m, chunkCeil(q));
          long r = CHUNK_BIT - 1 - (q + CHUNK_BIT - 1) % CHUNK_BIT;
#ifdef CGAL_CORE_DEBUG
          CGAL_assertion(r >= 0);
#endif

          err = 1 >> r;
          exp = (x.exp >> 1) - chunkCeil(qq);
          normal();
        }
      }  // end of case with error in mantissa
    }//else
  } else
    CGAL_error_msg("BigFloat error: squareroot called with negative operand.");
} //sqrt with initial approximation

//  compareMExp(x)
//    returns  1 if *this > x
//             0 if *this = x,
//            -1 if *this < x,
//
//      Main comparison method for BigFloat
//      This is called by BigFloat::compare()
//      BE CAREFUL:  The error bits are ignored!
//      Need another version if we want to take care of error bits

CGAL_INLINE_FUNCTION
int BigFloatRep :: compareMExp(const BigFloatRep& x) const {
  int st = sign(m);
  int sx = sign(x.m);

  if (st > sx)
    return 1;
  else if (st == 0 && sx == 0)
    return 0;
  else if (st < sx)
    return - 1;
  else { //  need to compare m && exp
    long expDiff = exp - x.exp;

    if (expDiff > 0) //  exp > x.exp
      return cmp(chunkShift(m, expDiff), x.m);
    else if (!expDiff)
      return cmp(m, x.m);
    else  //  exp < x.exp
      return cmp(m, chunkShift(x.m, - expDiff));
  }
}

// 3/6/2000:
// This is a private function used by BigFloatRep::operator<<
// to get the exact value
// of floor(log10(M * 2^ e)) where E is an initial guess.
// We will return the correct E which satisfies
//              10^E <= M * 2^e < 10^{E+1}
// But we convert this into
//              mm <= M < 10.mm

CGAL_INLINE_FUNCTION
long BigFloatRep :: adjustE( long E, BigInt M, long ee) const {
  if (M<0)
    M=-M;
  BigInt mm(1);
  if (ee > 0)
    M = (M<<static_cast<unsigned long>(ee));
  else
    mm = (mm << static_cast<unsigned long>(-ee));
  if (E > 0)
    mm *= (FiveTo(E)<< static_cast<unsigned long>(E));
  else
    M *= (FiveTo(-E) << static_cast<unsigned long>(-E));

  if (M < mm) {
    do {
      E--;
      M *= 10;
    } while (M < mm);
  } else if (M >= 10*mm) {
    mm *= 10;
    do {
      E++;
      mm *= 10;
    } while (M >= mm);
  }
  return E;
}

CGAL_INLINE_FUNCTION
BigFloatRep::DecimalOutput
BigFloatRep::toDecimal(unsigned int width, bool Scientific) const {
  BigFloatRep::DecimalOutput decOut;                // to be returned
  if (err > 0) {
    decOut.isExact = false;
  } else { // err == 0
    decOut.isExact = true;
  }

  if (err > 0 && err >= abs(m)) {
    // if err is larger than mantissa, sign and significant values
    // can not be determined.
    CGAL_CORE_warning_msg(true, "BigFloat error: Error is too big!");
    decOut.rep = "0.0e0";          // error is too big
    decOut.isScientific = false;
    decOut.noSignificant = 0;
    decOut.errorCode = 1;          // sign of this number is unknown
    return decOut;
  }

  decOut.sign = sign(m);
  decOut.errorCode = 0;

  BigInt M(m);                  // temporary mantissa
  long lm = bitLength(M);       // binary length of mantissa
  long e2 = bits(exp);          // binary shift length represented by exponent
  long le = clLg(err);          // binary length of err
  if (le == -1)
    le = 0;

  long L10 = 0;
  if (M != 0) {
    L10 = (long)std::floor((lm + e2) / lgTenM);
    L10 = adjustE(L10, m, e2);     // L10: floor[log10(M 2^(e2))], M != 0
  } else {
    L10 = 0;
  }
  // Convention: in the positional format, when the output is
  // the following string of 8 characters:
  //             (d0, d1, d2, d3, ".", d4, d5, d6, d7)
  // then the decimal point is said to be in the 4th position.
  // E.g., (d0, ".", d1, d2) has the decimal point in the 1st position.
  // The value of L10 says that the decimal point of output should be at
  // the (L10 + 1)st position. This is
  // true regardingless of whether M = 0 or not. For zero, we output
  // {0.0*} so L10=0.  In general, the |value| is less than 10
  // if and only if L10 is 0 and the
  // decimal point is in the 1st place.  Note that L10 is defined even if
  // the output is an integer (in which case it does not physically appear
  // but conceptually terminates the sequence of digits).

  // First, get the decimal representation of (m * B^(exp)).
  if (e2 < 0) {
    M *= FiveTo(-e2); // M = x * 10^(-e2)
  } else if (e2 > 0) {
    M <<= e2;         // M = x * 2^(e2)
  }

  std::string decRep = M.convert_to<std::string>();
  // Determine the "significant part" of this string, i.e. the part which
  // is guaranteed to be correct in the presence of error,
  // except that the last digit which might be subject to +/- 1.

  if (err != 0) {     // valid = number of significant digits
    unsigned long valid = floorlg10(m) - (long)std::floor(std::log10(float(err)));
    if (decRep.length() > valid) {
      decRep.erase(valid);
    }
  }

  // All the digits in decM are correct, except the last one might
  // subject to an error +/- 1.

  if ((decRep[0] == '+') || (decRep[0] == '-')) {
    decRep.erase(0, 1);
  }

  // Second, make choice between positional representation
  // and scientific notation.  Use scientific notation when:
  // 0) if scientific notation flag is on
  // 1) err * B^exp >= 1, the error contribute to the integral part.
  // 2) (1 + L10) >= width, there is not have enough places to hold the
  //    positional representation, not including decimal point.
  // 3) The distance between the first significant digit and decimal
  //    point is too large for the width limit. This is equivalent to
  //            Either ((L10 >= 0 and (L10 + 1) > width))
  //            Or  ((L10 < 0) and (-L10 + 1) > width).

  if (Scientific ||
      ((err > 0) && (le + e2) >= 0) ||          // if err*B^exp >= 1
      ((L10 >= 0) && (L10 + 1 >= (long)width )) ||
      ((L10 < 0) && (-L10 + 1 > (long)width ))) {
    // use scientific notation
    decRep = round(decRep, L10, width);
    decOut.noSignificant = width;
    decRep.insert(1, ".");
    if (L10 != 0) {
      decRep += 'e';
      if (L10 > 0) {
        decRep += '+';
      } else { // L10 < 0
        decRep += '-';
      }
      std::ostringstream oss;
      oss << labs(L10);
      decRep += oss.str();
      decOut.isScientific = true;
    }
  } else {
    // use conventional positional notation.
    if (L10 >= 0) { // x >= 1 or x == 0 and L10 + 1 <= width
      // round when necessary
      if (decRep.length() > width ) {
        decRep = round(decRep, L10, width );
        if (decRep.length() > width ) {
          // overflow happens! use scientific notation
          return toDecimal(width, true);
        }
      }
      decOut.noSignificant = static_cast<int>(decRep.length());
      if (L10 + 1 < (long)width ) {
        decRep.insert(L10 + 1, ".");
      } else { // L10 + 1 == width
        // do nothing
      }
    } else { // L10 < 0, 0 < x < 1
      // (-L10) leading zeroes, including one to the left of decimal dot
      // need to be added in beginning.
      decRep = std::string(-L10, '0') + decRep;
      // then round when necessary
      if (decRep.length() > width ) {
        decRep = round(decRep, L10, width );
        // cannot overflow since there are L10 leading zeroes.
      }
      decOut.noSignificant = static_cast<int>(decRep.length() - (-L10));
      decRep.insert(1, ".");
    }
    decOut.isScientific = false;
  }
#ifdef CGAL_CORE_DEBUG
  CGAL_assertion(decOut.noSignificant >= 0);
#endif

  decOut.rep = decRep;
  return decOut;
}//toDecimal

CGAL_INLINE_FUNCTION
std::string BigFloatRep::round(std::string inRep, long& L10, unsigned int width) const {
  // round inRep so that the length would not exceed width.
  if (inRep.length() <= width)
    return inRep;

  int i = width; // < length
  bool carry = false;

  if ((inRep[i] >= '5') && (inRep[i] <= '9')) {
    carry = true;
    i--;
    while ((i >= 0) && carry) {
      if (carry) {
        inRep[i] ++;
        if (inRep[i] > '9') {
          inRep[i] = '0';
          carry = true;
        } else {
          carry = false;
        }
      }
      i-- ;
    }

    if ((i < 0) && carry) { // overflow
      inRep.insert(inRep.begin(), '1');
      L10 ++;
      width ++;
    }
  }

  return inRep.substr(0, width);
}//round(string,width)


// This function fromString(str, prec) is similar to the
//      constructor Real(char * str, extLong prec)
// See the file Real.cc for the differences

CGAL_INLINE_FUNCTION
void BigFloatRep :: fromString(const char *str, extLong prec ) {
  // NOTE: prec defaults to get_static_defBigFloatInputDigits() (see BigFloat.h)
  // check that prec is not INFTY
  if (prec.isInfty())
    CGAL_error_msg("BigFloat error: infinite precision not allowed");

  const char *e = strchr(str, 'e');
  int dot = 0;
  long e10 = 0;
  if (e != nullptr)
    e10 = atol(e+1);    // e10 is decimal precision of the input string
  // i.e., input is A/10^{e10}.
  else {
    e = str + strlen(str);
#ifdef CGAL_CORE_DEBUG
    CGAL_assertion(*e == '\0');
#endif

  }

  const char *p = str;
  if (*p == '-' || *p == '+')
    p++;
  m = 0;
  exp = 0;

  for (; p < e; p++) {
    if (*p == '.') {
      dot = 1;
      continue;
    }
    m = m * 10 + (*p - '0');
    if (dot)
      e10--;
  }

  BigInt one = 1;
  long t = (e10 < 0) ? -e10 : e10;
  BigInt ten = FiveTo(t) * (one << static_cast<unsigned long>(t));

  // HERE IS WHERE WE USE THE SYSTEM CONSTANT
  //           defBigFloatInputDigits
  // Note: this constant is rather similar to defInputDigits which
  //     is used by Real and Expr for controlling
  //     input accuracy.  The difference is that defInputDigits can
  //     be CORE_INFTY, but defBigFloatInputDigits must be finite.

  if (e10 < 0)
    div(m, ten, CORE_posInfty, 4 * prec);
  else
    m *= ten;
  if (*str == '-')
    m = -m;
}//BigFloatRep::fromString

CGAL_INLINE_FUNCTION
std::istream& BigFloatRep :: operator >>(std::istream& i) {
  int size = 20;
  std::string str;
  str.reserve(size);
  char c;
  int d = 0, e = 0, s = 0;
  // d=1 means dot is found
  // e=1 means 'e' or 'E' is found
  //  int done = 0;

  // Chen Li: fixed a bug, the original statement is
  //  for (i.get(c); c == ' '; i.get(c));
  // use isspace instead of testing c == ' ', since it must also
  // skip tab, catridge/return, etc.
  // Change to:
  //  int status;
  do {
    i.get(c);
  } while (isspace(c)); /* loop if met end-of-file, or
                               char read in is white-space. */
  // Chen Li, "if (c == EOF)" is unsafe since c is of char type and
  // EOF is of int type with a negative value -1
  if (i.eof()) {
    i.clear(std::ios::eofbit | std::ios::failbit);
    return i;
  }

  // the current content in "c" should be the first non-whitespace char
  if (c == '-' || c == '+') {
    str += c;
    i.get(c);
  }

  for (; isdigit(c) || (!d && c=='.') ||
       (!e && ((c=='e') || (c=='E'))) || (!s && (c=='-' || c=='+')); i.get(c)) {
    if (!e && (c == '-' || c == '+'))
      break;
    // Chen Li: put one more rule to prohibite input like
    //  xxxx.xxxe+xxx.xxx:
    if (e && (c == '.'))
      break;

    str += c;
    if (c == '.')
      d = 1;
    // Chen Li: fix a bug -- the sign of exponent can not happen before
    // the character "e" appears! It must follow the "e" actually.
    //    if (e || c == '-' || c == '+') s = 1;
    if (e)
      s = 1;
    if ((c == 'e') || (c=='E'))
      e = 1;
  }

  i.putback(c);
  fromString(str.c_str());
  return i;
}//operator >>


// BigFloatRep::toDouble()
//      converts the BigFloat to a machine double
//      This is a dangerous function as the method
//      is silent when it does not fit into a machine double!
// ToDo: fix this by return a machine NaN, +/- Infinity, +/- 0,
//      when appropriate.
//      Return NaN when error is larger than mantissa
//      Return +/- Infinity if BigFloat is too big
//      Return +/- 0 if BigFloat is too small
#ifdef _MSC_VER
#pragma warning(disable: 4723)
#endif
CGAL_INLINE_FUNCTION
double BigFloatRep :: toDouble() const {
  if (m == 0)
    return (sign(m) * 0.0);

  long e2 = bits(exp);
  long le = clLg(err);  // if err=0, le will be -1
  if (le == -1)
    le = 0;

  BigInt M = m >> static_cast<unsigned long>(le);// remove error bits in mantissa

  // Below, we want to return NaN by computing 0.0/0.0.
  // To avoid compiler warnings about divide by zero, we do this:

  double foolCompilerZero;
  foolCompilerZero = 0.0;

  // COMMENT: we should directly store the
  //    special IEEE values NaN, +/-Infinity, +/-0 in the code!!

  if (M == 0)
    return ( 0.0/foolCompilerZero ) ; // return NaN

  e2 += le;             // adjust exponent corresponding to error bits

  int len = bitLength(M) - 53;  // this is positive if M is too large

  if (len > 0) {
    M >>= len;
    e2 += len;
  }

  double tt = doubleValue(M);

  int ee = e2 + bitLength(M) - 1; // real exponent.

  if (ee >= 1024)       // overflow!
    return (  sign(m)/foolCompilerZero  );      // return a signed infinity

  if (ee <= -1075)      // underflow!
    // NOTE: if (-52 < ee <= 0) get denormalized number
    return ( sign(m) * 0.0 );  // return signed zero.

  // Execute this loop if e2 < 0;
  for (int i = 0; i > e2; i--)
    tt /= 2;

  // Execute this loop if e2 > 0;
  for (int j = 0; j < e2; j++)
    tt *= 2;

  return tt;
}//toDouble
#ifdef _MSC_VER
#pragma warning(default: 4723)
#endif
CGAL_INLINE_FUNCTION
BigInt BigFloatRep::toBigInt() const {
  long e2 = bits(exp);
  long le = clLg(err);
  if (le == -1)
    le = 0;
#ifdef CGAL_CORE_DEBUG
  CGAL_assertion (le >= 0);
#endif

  BigInt M = m >> static_cast<unsigned long>(le); // discard the contaminated bits.
  e2 += le;           // adjust the exponent

  if (e2 < 0)
    return M >> static_cast<unsigned long>(-e2);
  else if (e2 > 0)
    return M << static_cast<unsigned long>(e2);
  else
    return M;
}

CGAL_INLINE_FUNCTION
long BigFloatRep :: toLong() const {
  // convert a BigFloat to a long integer, rounded toward -\infty.
  long e2 = bits(exp);
  long le = clLg(err);
#ifdef CGAL_CORE_DEBUG
  CGAL_assertion (le >= 0);
#endif

  BigInt M = m >> static_cast<unsigned long>(le); // discard the contaminated bits.
  e2 += le;           // adjust the exponent
  long t;
  if (e2 < 0)
    t = ulongValue(M >> static_cast<unsigned long>(-e2));
  else if (e2 > 0)
    t = ulongValue(M << static_cast<unsigned long>(e2));
  else
    t = ulongValue(M);
  // t = M.as_long();
  // Note: as_long() will return LONG_MAX in case of overflow.

  return t;
}

// pow(r,n) function for BigFloat
// Note: power(r,n) calls pow(r,n)
CGAL_INLINE_FUNCTION
BigFloat pow(const BigFloat& r, unsigned long n) {
  if (n == 0)
    return BigFloat(1);
  else if (n == 1)
    return r;
  else {
    BigFloat x = r;
    while ((n % 2) == 0) { // n is even
      x *= x;
      n >>= 1;
    }
    BigFloat u = x;
    while (true) {
      n >>= 1;
      if (n == 0)
        return u;
      x *= x;
      if ((n % 2) == 1) // n is odd
        u *= x;
    }
    //return u; // unreachable
  }
}//pow

// experimental
CGAL_INLINE_FUNCTION
BigFloat root(const BigFloat& x, unsigned long k,
         const extLong& a, const BigFloat& A) {
  if (x.sign() == 0) {
    return BigFloat(0);
  } else if (x == 1) {
    return BigFloat(1);
  } else  {
    BigFloat q, del, zz;
    BigFloat z = A;
    BigFloat bk = long(k);
    for (; ;) {
      zz = pow(z, k-1);
      q = x.div(zz, a);
      q.makeExact();
      del = z - q;
      del.makeExact();
      if (del.MSB() < -a)
        break;
      z = ((bk-1)*z + q).div(bk, a);
          // newton's iteration: z_{n+1}=((k-1)z_n+x/z_n^{k-1})/k
      z.makeExact();
    }
    return z;
  }
}//root

  CORE_MEMORY_IMPL(BigFloatRep)

} //namespace CORE

#include <CGAL/enable_warnings.h>