SPAR/AIModule/BWTA/vendors/CGAL/CGAL/CORE/BigFloat.h

Go to the documentation of this file.

/****************************************************************************
 * Core Library Version 1.7, August 2004
 * Copyright (c) 1995-2004 Exact Computation Project
 * All rights reserved.
 *
 * This file is part of CORE (http://cs.nyu.edu/exact/core/); you may
 * redistribute it under the terms of the Q Public License version 1.0.
 * See the file LICENSE.QPL distributed with CORE.
 *
 * Licensees holding a valid commercial license may use this file in
 * accordance with the commercial license agreement provided with the
 * software.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * File: BigFloat.h
 * Synopsis: 
 *              An implementation of BigFloat numbers with error bounds.
 * 
 * Written by 
 *       Chee Yap <yap@cs.nyu.edu>
 *       Chen Li <chenli@cs.nyu.edu>
 *       Zilin Du <zilin@cs.nyu.edu>
 *
 * WWW URL: http://cs.nyu.edu/exact/
 * Email: exact@cs.nyu.edu
 *
 * $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.5-branch/Core/include/CGAL/CORE/BigFloat.h $
 * $Id: BigFloat.h 37563 2007-03-27 13:28:17Z afabri $
 ***************************************************************************/

#ifndef _CORE_BIGFLOAT_H_
#define _CORE_BIGFLOAT_H_

#include <CGAL/CORE/BigFloatRep.h>

CORE_BEGIN_NAMESPACE

class Expr;

typedef RCImpl<BigFloatRep> RCBigFloat;
class BigFloat : public RCBigFloat {
public:


  BigFloat() : RCBigFloat(new BigFloatRep()) {}
  BigFloat(short i) : RCBigFloat(new BigFloatRep(i)) {}
  BigFloat(float i) : RCBigFloat(new BigFloatRep(i)) {}
  BigFloat(int i) : RCBigFloat(new BigFloatRep(i)) {}
  BigFloat(long l) : RCBigFloat(new BigFloatRep(l)) {}
  BigFloat(double d) : RCBigFloat(new BigFloatRep(d)) {}
  BigFloat(const char* s) : RCBigFloat(new BigFloatRep(s)) {}
  BigFloat(const std::string& s) : RCBigFloat(new BigFloatRep(s)) {}

  //     This is a hack because in Sturm, we need to approximate any
  //     coefficient type NT to a BigFloat, and it would complain if we
  //     do not have this method explicitly:
  BigFloat(int& i, const extLong& /*r*/, const extLong& /*a*/)
      : RCBigFloat(new BigFloatRep(i)) {}
  BigFloat(long& x, const extLong& /*r*/, const extLong& /*a*/)
      : RCBigFloat(new BigFloatRep(x)) {}
  BigFloat(const BigInt& I, unsigned long er, long ex)
      : RCBigFloat(new BigFloatRep(I, er, ex)) {}
  BigFloat(const BigInt& I, long ex)
      : RCBigFloat(new BigFloatRep(I, ex)) {}
  BigFloat(const BigInt& I)
      : RCBigFloat(new BigFloatRep(I)) {}
  BigFloat(const BigRat& R, const extLong& r = defRelPrec,
           const extLong& a = defAbsPrec)
      : RCBigFloat(new BigFloatRep()) {
    rep->approx(R, r, a);
  }

  // REMARK: it is somewhat against our principles to have BigFloat
  // know about Expr, but BigFloat has a special role in our system!
  // ===============================
  explicit BigFloat(const Expr& E, const extLong& r = defRelPrec,
           const extLong& a = defAbsPrec);

  //Dummy
  explicit BigFloat(const BigFloat& E, const extLong& ,
           const extLong&): RCBigFloat(E) {
    rep->incRef();
  }

  explicit BigFloat(BigFloatRep* r) : RCBigFloat(new BigFloatRep()) {
    rep = r;
  }



  BigFloat(const BigFloat& rhs) : RCBigFloat(rhs) {
    rep->incRef();
  }
  BigFloat& operator=(const BigFloat& rhs) {
    if (this != &rhs) {
      rep->decRef();
      rep = rhs.rep;
      rep->incRef();
    }
    return *this;
  }
  ~BigFloat() {
    rep->decRef();
  }



  BigFloat& operator+= (const BigFloat& x) {
    BigFloat z;
    z.rep->add(*rep, *x.rep);
    *this = z;
    return *this;
  }
  BigFloat& operator-= (const BigFloat& x) {
    BigFloat z;
    z.rep->sub(*rep, *x.rep);
    *this = z;
    return *this;
  }
  BigFloat& operator*= (const BigFloat& x) {
    BigFloat z;
    z.rep->mul(*rep, *x.rep);
    *this = z;
    return *this;
  }
  BigFloat& operator/= (const BigFloat& x) {
    BigFloat z;
    z.rep->div(*rep, *x.rep, defBFdivRelPrec);
    *this = z;
    return *this;
  }



  BigFloat operator+() const {
    return BigFloat(*this);
  }
  BigFloat operator-() const {
    return BigFloat(-rep->m, rep->err, rep->exp);
  }



  void fromString(const char* s, const extLong& p=defBigFloatInputDigits) {
    rep->fromString(s, p);
  }
  std::string toString(long prec=defBigFloatOutputDigits, bool sci=false) const {
    return rep->toString(prec, sci);
  }
  std::string str() const {
    return toString();
  }



  int intValue() const {
    return static_cast<int>(rep->toLong());
  }
  long longValue() const {
    long l = rep->toLong();
    if ((l == LONG_MAX) || (l == LONG_MIN))
      return l; // return the overflown value.
    if ((sign() < 0) && (cmp(BigFloat(l)) != 0)) {
      // a negative value not exactly rounded.
      l--; // rounded to floor.
    }
    return l;
  }
  float floatValue() const {
    return static_cast<float>(rep->toDouble());
  }
  double doubleValue() const {
    return rep->toDouble();
  }
  BigInt BigIntValue() const {
    return rep->toBigInt();
  }
  BigRat BigRatValue() const {
    return rep->BigRatize();
  }



  static bool hasExactDivision() {
    return false;
  }
  
  //CONSTANTS
  static const BigFloat& getZero();
  static const BigFloat& getOne();


  int sign() const {
    assert((err() == 0 && m() == 0) || !(isZeroIn()));
    return rep->signM();
  }

  bool isZeroIn() const {
    return rep->isZeroIn();
  }
  BigFloat abs() const {
    return (sign()>0) ? +(*this) : -(*this);
  }
  int cmp(const BigFloat& x) const {
    return rep->compareMExp(*x.rep);
  }

  const BigInt& m() const {
    return rep->m;
  }
  unsigned long err() const {
    return rep->err;
  }
  long exp() const {
    return rep->exp;
  }


  bool isExact() const {
    return rep->err == 0;
  }

  BigFloat& makeExact() {
    makeCopy();
    rep->err =0;
    return *this;
  }

  BigFloat& makeCeilExact() {
    makeCopy();
    rep->m += rep->err;
    rep->err =0;
    return *this;
  }

  BigFloat& makeFloorExact() {
    makeCopy();
    rep->m -= rep->err;
    rep->err =0;
    return *this;
  }

  BigFloat& makeInexact() {
    makeCopy();
    rep->err =1;
    return *this;
  }

  extLong lMSB() const {
    return rep->lMSB();
  }
  extLong uMSB() const {
    return rep->uMSB();
  }
  extLong MSB() const {
    return rep->MSB();
  }

  extLong flrLgErr() const {
    return rep->flrLgErr();
  }
  extLong clLgErr() const {
    return rep->clLgErr();
  }

  BigFloat div(const BigFloat& x, const extLong& r) const {
    BigFloat y;
    y.rep->div(*rep, *x.rep, r);
    return y;
  }
  BigFloat div2() const {
    BigFloat y;
    y.rep->div2(*rep);
    return y;
  }

  BigFloat sqrt(const extLong& a) const {
    BigFloat x;
    x.rep->sqrt(*rep, a);
    return x;
  }
  BigFloat sqrt(const extLong& a, const BigFloat& init) const {
    BigFloat x;
    x.rep->sqrt(*rep, a, init);
    return x;
  }



  void approx(const BigInt& I, const extLong& r, const extLong& a) {
    makeCopy();
    rep->trunc(I, r, a);
  }
  void approx(const BigFloat& B, const extLong& r, const extLong& a) {
    makeCopy();
    rep->approx(*B.rep, r, a);
  }
  void approx(const BigRat& R, const extLong& r, const extLong& a) {
    makeCopy();
    rep->approx(R, r, a);
  }
  void dump() const {
    rep->dump();
  }

  static BigFloat exp2(int e) {
    return BigFloat(BigFloatRep::exp2(e));
  }

}; // class BigFloat




void readFromFile(BigFloat& bf, std::istream& in, long maxLength = 0);
void writeToFile(const BigFloat& bf, std::ostream& in,
                int base=10, int charsPerLine=80);

inline std::ostream& operator<< (std::ostream& o, const BigFloat& x) {
  x.getRep().operator<<(o);
  return o;
}
inline std::istream& operator>> (std::istream& i, BigFloat& x) {
  x.makeCopy();
  x.getRep().operator>>(i);
  return i;
}

inline BigFloat operator+ (const BigFloat& x, const BigFloat& y) {
  BigFloat z;
  z.getRep().add(x.getRep(), y.getRep());
  return z;
}
inline BigFloat operator- (const BigFloat& x, const BigFloat& y) {
  BigFloat z;
  z.getRep().sub(x.getRep(), y.getRep());
  return z;
}
inline BigFloat operator* (const BigFloat& x, const BigFloat& y) {
  BigFloat z;
  z.getRep().mul(x.getRep(), y.getRep());
  return z;
}
inline BigFloat operator/ (const BigFloat& x, const BigFloat& y) {
  BigFloat z;
  z.getRep().div(x.getRep(),y.getRep(),defBFdivRelPrec);
  return z;
}

inline bool operator== (const BigFloat& x, const BigFloat& y) {
  return x.cmp(y) == 0;
}
inline bool operator!= (const BigFloat& x, const BigFloat& y) {
  return x.cmp(y) != 0;
}
inline bool operator>= (const BigFloat& x, const BigFloat& y) {
  return x.cmp(y) >= 0;
}
inline bool operator> (const BigFloat& x, const BigFloat& y) {
  return x.cmp(y) > 0;
}
inline bool operator<= (const BigFloat& x, const BigFloat& y) {
  return x.cmp(y) <= 0;
}
inline bool operator< (const BigFloat& x, const BigFloat& y) {
  return x.cmp(y) < 0;
}

inline int sign(const BigFloat& x) {
  return x.sign();
}
inline BigFloat div2(const BigFloat& x){
  return x.div2();
}
inline BigFloat abs(const BigFloat& x) {
  return x.abs();
}
inline int cmp(const BigFloat& x, const BigFloat& y) {
  return x.cmp(y);
}
BigFloat pow(const BigFloat&, unsigned long);
inline BigFloat power(const BigFloat& x, unsigned long p) {
  return pow(x, p);
}
BigFloat root(const BigFloat&, unsigned long k, const extLong&, const BigFloat&);
inline BigFloat root(const BigFloat& x, unsigned long k) {
  return root(x, k, defBFsqrtAbsPrec, x);
}

inline BigFloat sqrt(const BigFloat& x) {
  return x.sqrt(defBFsqrtAbsPrec);
}

inline BigFloat centerize(const BigFloat& a, const BigFloat& b) {
  BigFloat z;
  z.getRep().centerize(a.getRep(), b.getRep());
  return z;
}

inline long minStar(long m, long n) {
  if (m*n <= 0) return 0;
  if (m>0) 
    return core_min(m, n);
  else 
    return core_max(m, n);
}



inline bool isDivisible(const BigFloat& a, const BigFloat& b) {
  // assert: a and b are exact BigFloats.
  if (sign(a.m()) == 0) return true;
  if (sign(b.m()) == 0) return false;
  unsigned long bin_a = getBinExpo(a.m());
  unsigned long bin_b = getBinExpo(b.m());
  
  BigInt m_a = a.m() >> bin_a;
  BigInt m_b = b.m() >> bin_b;
  long e_a = bin_a + BigFloatRep::bits(a.exp());
  long e_b = bin_b + BigFloatRep::bits(b.exp());
  long dx = minStar(e_a, e_b);

  return isDivisible(m_a, m_b) && (dx == e_b); 
}

inline bool isDivisible(double x, double y) {
  //Are these exact?
  return isDivisible(BigFloat(x), BigFloat(y)); 
}


// Chee (8/1/2004)   The definition of div_exact(x,y) 
//   ensure that Polynomials<NT> works with NT=BigFloat and NT=double.
//Bug: We should first normalize the mantissas of the Bigfloats and
//then do the BigInt division. For e.g. 1 can be written as 2^{14}*2{-14}.
//Now if we divide 2 by one using this representation of one and without
// normalizing it then we get zero.
inline BigFloat div_exact(const BigFloat& x, const BigFloat& y) {
  BigInt z;
  assert (isDivisible(x,y));
  unsigned long bin_x = getBinExpo(x.m());
  unsigned long bin_y = getBinExpo(y.m());

  BigInt m_x = x.m() >> bin_x;
  BigInt m_y = y.m() >> bin_y;
  long e_x = bin_x + BigFloatRep::bits(x.exp());
  long e_y = bin_y + BigFloatRep::bits(y.exp());
  //Since y divides x, e_y = minstar(e_x, e_y)
  z = div_exact(m_x, m_y);

  //  mpz_divexact(z.get_mp(), x.m().get_mp(), y.m().get_mp()); THIS WAS THE BUG
  // assert: x.exp() - y.exp() does not under- or over-flow.
  return BigFloat(z, e_x - e_y);  
}

inline BigFloat div_exact(double x, double y) {
  return div_exact(BigFloat(x), BigFloat(y));
}
// Remark: there is another notion of "exact division" for BigFloats,
//      and that is to make the division return an "exact" BigFloat
//      i.e., err()=0.  

inline BigFloat gcd(const BigFloat& a, const BigFloat& b) {
  if (sign(a.m()) == 0) return core_abs(b);
  if (sign(b.m()) == 0) return core_abs(a);

  BigInt r;
  long dx;
  unsigned long bin_a = getBinExpo(a.m());
  unsigned long bin_b = getBinExpo(b.m());

/* THE FOLLOWING IS ALTERNATIVE CODE, for GCD using base B=2^{14}:
 *std::cout << "bin_a=" << bin_a << ",bin_b=" << bin_b << std::endl;
  std::cout << "a.exp()=" << a.exp() << ",b.exp()=" << b.exp() << std::endl;
  long chunk_a = BigFloatRep::chunkFloor(bin_a);
  long chunk_b = BigFloatRep::chunkFloor(bin_b);
  BigInt m_a = BigFloatRep::chunkShift(a.m(), chunk_a);
  BigInt m_b = BigFloatRep::chunkShift(b.m(), chunk_b);

  r = gcd(m_a, m_b);
  dx = minStar(chunk_a + a.exp(), chunk_b + b.exp());
*/
  BigInt m_a = a.m() >> bin_a;
  BigInt m_b = b.m() >> bin_b;
  r = gcd(m_a, m_b);
  dx = minStar(bin_a + BigFloatRep::bits(a.exp()),
                  bin_b + BigFloatRep::bits(b.exp()));

  long chunks = BigFloatRep::chunkFloor(dx);
  r <<= (dx - BigFloatRep::bits(chunks));
  dx = chunks;

  return BigFloat(r, 0, dx);
}

// Not needed for now:
// inline void div_rem(BigFloat& q, BigFloat& r,
//      const BigFloat& a, const BigFloat& b) {
  //q.makeCopy();
  //r.makeCopy();
  //mpz_tdiv_qr(q.get_mp(), r.get_mp(), a.get_mp(), b.get_mp());
//}//


// constructor BigRat from BigFloat
inline BigRat::BigRat(const BigFloat& f) : RCBigRat(new BigRatRep()){
  *this = f.BigRatValue();
}
CORE_END_NAMESPACE
#endif // _CORE_BIGFLOAT_H_