SPAR/AIModule/BWTA/vendors/CGAL/CGAL/CORE/Filter.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: Filter.h
 * Synopsis:
 *      This is a simple filtered floating point number,
 *      represented by the main class, FilterFp.
 *      based on the Burnikel-Funke-Schirra (BFS) filter scheme.
 *      We do not use IEEE exception mechanism here.
 *      It is used by the Expr class.
 * 
 * Written by 
 *       Zilin Du <zilin@cs.nyu.edu>
 *       Chee Yap <yap@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/Filter.h $
 * $Id: Filter.h 45042 2008-08-20 11:14:00Z spion $
 ***************************************************************************/

#ifndef _CORE_FILTER_H_
#define _CORE_FILTER_H_

#include <CGAL/config.h>
#include <CGAL/CORE/Real.h>
#include <cmath>

#if !defined CGAL_CFG_NO_CPP0X_ISFINITE
  #define finite(x)     std::isfinite(x)
#elif defined (_MSC_VER) || defined (__MINGW32__) // add support for MinGW
  #define finite(x)     _finite(x)
  #define ilogb(x)      (int)_logb(x)
#endif

#if defined(sun) || defined(__sun)
  #include <ieeefp.h>
#endif

CORE_BEGIN_NAMESPACE

const int POWTWO_26 = (1 << 26);  

class filteredFp {
  double fpVal;         // approximate double value for some "real value"
  double maxAbs;        // if (|fpVal| > maxAbs * ind * 2^{-53}) then
  int ind;              // sign of value is sign(fpVal).  Else, don't know.
  // REFERENCE: Burnikel, Funke, Schirra (BFS) filter
  // Chee: in isOK(), you used the test "|fpVal| >= maxAbs * ind * 2^{-53}" 
  // which seems to be correct (i.e., not |fpVal| > maxAbs * ind * 2^{-53})
public:


  filteredFp (double val = 0.0)
      : fpVal(val), maxAbs(core_abs(val)), ind(0) {}
  filteredFp (double val, double m, int i)
      : fpVal(val), maxAbs(m), ind(i) {}


  filteredFp (const Real & value) : fpVal(0.0), maxAbs(0.0), ind(0) {
    if (value != CORE_REAL_ZERO) {
      ind = 1;
      fpVal = value.doubleValue();
      if (value.MSB() <= -1075)
        maxAbs = 1;
      else      
        maxAbs = core_abs(fpVal); // NaN are propagated correctly by core_abs.
    }
  }



  double getValue() const {
    return fpVal;
  }
  bool isOK() const {
    return (fpFilterFlag  && // To disable filter
            finite(fpVal) && // Test for infinite and NaNs
            (core_abs(fpVal) >= maxAbs*ind*CORE_EPS));
  }

  int sign() const {
#ifdef CORE_DEBUG
    assert(isOK());
#endif
    if (fpVal == 0.0)
      return 0;
    else
      return fpVal > 0.0 ? 1: -1;
  }

  extLong lMSB() const {
    return extLong(ilogb(core_abs(fpVal)-maxAbs*ind*CORE_EPS));
  }
  extLong uMSB() const {
    return extLong(ilogb(core_abs(fpVal)+maxAbs*ind*CORE_EPS));
  }



  filteredFp operator -() const {
    return filteredFp(-fpVal, maxAbs, ind);
  }
  filteredFp operator+ (const filteredFp& x) const {
    return filteredFp(fpVal+x.fpVal, maxAbs+x.maxAbs, 1+core_max(ind, x.ind));
  }
  filteredFp operator- (const filteredFp& x) const {
    return filteredFp(fpVal-x.fpVal, maxAbs+x.maxAbs, 1+core_max(ind, x.ind));
  }
  filteredFp operator* (const filteredFp& x) const {
    return filteredFp(fpVal*x.fpVal, maxAbs*x.maxAbs+DBL_MIN, 1+ind+x.ind);
  }
  filteredFp operator/ (const filteredFp& x) const {
    if (x.fpVal == 0.0)
      core_error("possible zero divisor!", __FILE__, __LINE__, false);
    double xxx = core_abs(x.fpVal) / x.maxAbs - (x.ind+1)*CORE_EPS + DBL_MIN;
    if (xxx > 0) {
      double val =  fpVal / x.fpVal;
      double maxVal = ( core_abs(val) + maxAbs / x.maxAbs) / xxx + DBL_MIN;
      return filteredFp(val, maxVal, 1 + core_max(ind, x.ind + 1));
    } else
      return filteredFp(getDoubleInfty(), 0.0, 0);
  }
  filteredFp sqrt () const {
    if (fpVal < 0.0)
      core_error("possible negative sqrt!", __FILE__, __LINE__, false);
    if (fpVal > 0.0) {
      double val = std::sqrt(fpVal);
      return filteredFp(val,  ( maxAbs / fpVal ) * val, 1 + ind);
    } else
      return filteredFp(0.0, std::sqrt(maxAbs) * POWTWO_26, 1 + ind);
  }

  void dump (std::ostream&os) const {
    os << "Filter=[fpVal=" << fpVal << ",maxAbs=" << maxAbs << ",ind=" << ind << "]";
  }

  static double getDoubleInfty() {
    static double d = DBL_MAX;
    return 2*d;
  }
}; //filteredFp class

inline std::ostream & operator<< (std::ostream & os, const filteredFp& fp) {
  fp.dump(os);
  return os;
}

CORE_END_NAMESPACE
#endif // _CORE_FILTER_H_