wrong_interval.cpp

10/24/2019 - 7:15 AM
#include <CGAL/Simple_cartesian.h>

#include <CGAL/Interval_nt.h>
#include <CGAL/leda_real.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

using EK = CGAL::Simple_cartesian<leda::real>;
using FK = CGAL::Simple_cartesian<CGAL::Interval_nt<true> >;

namespace CGAL {

template < typename EFT, typename FFT>
bool is_exact_in_interval(const EFT e, const FFT f)
{
  if(!is_finite(f.inf()) && !is_finite(f.sup()))
  {
//    std::cerr << "infinite interval..." << std::endl;
    return true; // not true if both are -inf or +inf, but doesn't matter we're never there with the data
  }

//  std::cout << "finiteness of the interval: " << is_finite(f.inf()) << " " << is_finite(f.sup()) << std::endl;

  if(!is_finite(f.inf()))
  {
//    std::cout << "infinite inf: " << (e <= EFT(f.sup())) << std::endl;
    return (e <= EFT(f.sup()));
  }
  else if(!is_finite(f.sup()))
  {
//    std::cout << "infinite sup: " << (EFT(f.inf()) <= e) << std::endl;
    return (EFT(f.inf()) <= e);
  }
  else
  {
//    std::cout << "general case: " << (EFT(f.inf()) <= e) << " " << (e <= EFT(f.sup())) << std::endl;
    return (EFT(f.inf()) <= e) && (e <= EFT(f.sup()));
  }
}

} // namespace CGAL

int main()
{
  const double px = -1.7976931348622648e+308;
  const double py =  5.3779408224704862e-299;
  const double pz =  2.5394974196240072e-321;

  const double qx =  5.3779407578037265e-299;
  const double qy =  8.7171284671510145e-121;
  const double qz =  5.3779407512700166e-299;

  const double rx =  5.3779414276905704e-299;
  const double ry =  5.3779407512682172e-299;
  const double rz =  7.1237266591579577e-307;

  const double sx =  9.4175561178482311e+298;
  const double sy =  5.3779428873390216e-299;
  const double sz =  5.3779407512681174e-299;

  EK::Point_3 ep = { px, py, pz };
  EK::Point_3 eq = { qx, qy, qz };
  EK::Point_3 er = { rx, ry, rz };
  EK::Point_3 es = { sx, sy, sz };

  FK::Point_3 fp = { px, py, pz };
  FK::Point_3 fq = { qx, qy, qz };
  FK::Point_3 fr = { rx, ry, rz };
  FK::Point_3 fs = { sx, sy, sz };

  CGAL_assertion(CGAL::is_exact_in_interval(ep.x(), fp.x()));
  CGAL_assertion(CGAL::is_exact_in_interval(ep.y(), fp.y()));
  CGAL_assertion(CGAL::is_exact_in_interval(ep.z(), fp.z()));
  CGAL_assertion(CGAL::is_exact_in_interval(eq.x(), fq.x()));
  CGAL_assertion(CGAL::is_exact_in_interval(eq.y(), fq.y()));
  CGAL_assertion(CGAL::is_exact_in_interval(eq.z(), fq.z()));
  CGAL_assertion(CGAL::is_exact_in_interval(er.x(), fr.x()));
  CGAL_assertion(CGAL::is_exact_in_interval(er.y(), fr.y()));
  CGAL_assertion(CGAL::is_exact_in_interval(er.z(), fr.z()));
  CGAL_assertion(CGAL::is_exact_in_interval(es.x(), fs.x()));
  CGAL_assertion(CGAL::is_exact_in_interval(es.y(), fs.y()));
  CGAL_assertion(CGAL::is_exact_in_interval(es.z(), fs.z()));

//    CGAL::Protect_FPU_rounding<true> protector;

  std::cout << "---- a00" << std::endl;
  EK::FT e_a00 = eq.x() - ep.x();
  FK::FT f_a00 = fq.x() - fp.x();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a00, f_a00));

  std::cout << "---- a10" << std::endl;
  EK::FT e_a10 = eq.y() - ep.y();
  FK::FT f_a10 = fq.y() - fp.y();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a10, f_a10));

  std::cout << "---- a20" << std::endl;
  EK::FT e_a20 = eq.z() - ep.z();
  FK::FT f_a20 = fq.z() - fp.z();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a20, f_a20));

  std::cout << "---- a01" << std::endl;
  EK::FT e_a01 = er.x() - ep.x();
  FK::FT f_a01 = fr.x() - fp.x();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a01, f_a01));

  std::cout << "---- a11" << std::endl;
  EK::FT e_a11 = er.y() - ep.y();
  FK::FT f_a11 = fr.y() - fp.y();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a11, f_a11));

  std::cout << "---- a21" << std::endl;
  EK::FT e_a21 = er.z() - ep.z();
  FK::FT f_a21 = fr.z() - fp.z();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a21, f_a21));

  std::cout << "---- a02" << std::endl;
  EK::FT e_a02 = es.x() - ep.x();
  FK::FT f_a02 = fs.x() - fp.x();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a02, f_a02));

  std::cout << "---- a12" << std::endl;
  EK::FT e_a12 = es.y() - ep.y();
  FK::FT f_a12 = fs.y() - fp.y();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a12, f_a12));

  std::cout << "---- a22" << std::endl;
  EK::FT e_a22 = es.z() - ep.z();
  FK::FT f_a22 = fs.z() - fp.z();
  CGAL_assertion(CGAL::is_exact_in_interval(e_a22, f_a22));

  std::cout << "---- m01" << std::endl;
  const EK::FT e_m01 = e_a00 * e_a11 - e_a10 * e_a01;
  const FK::FT f_m01 = f_a00 * f_a11 - f_a10 * f_a01;
  CGAL_assertion(CGAL::is_exact_in_interval(e_m01, f_m01));

  std::cout << "---- m02" << std::endl;
  const EK::FT e_m02 = e_a00 * e_a21 - e_a20 * e_a01;
  const FK::FT f_m02 = f_a00 * f_a21 - f_a20 * f_a01;
  CGAL_assertion(CGAL::is_exact_in_interval(e_m02, f_m02));

  std::cout << "---- m12" << std::endl;
  const EK::FT e_m12 = e_a10 * e_a21 - e_a20 * e_a11;
  const FK::FT f_m12 = f_a10 * f_a21 - f_a20 * f_a11;
  CGAL_assertion(CGAL::is_exact_in_interval(e_m12, f_m12));

  std::cout << "---- part_m012_1" << std::endl;
  const EK::FT e_m012_1 = e_m01 * e_a22;
  const FK::FT f_m012_1 = f_m01 * f_a22;
  CGAL_assertion(CGAL::is_exact_in_interval(e_m012_1, f_m012_1));

  std::cout << "---- part_m012_2" << std::endl;
  const EK::FT e_m012_2 = e_m02 * e_a12;
  const FK::FT f_m012_2 = f_m02 * f_a12;
  CGAL_assertion(CGAL::is_exact_in_interval(e_m012_2, f_m012_2));

  std::cout << "---- part_m012_3" << std::endl;
  const EK::FT e_m012_3 = e_m12 * e_a02;
  const FK::FT f_m012_3 = f_m12 * f_a02;
  std::cout << "e_m12: " << e_m12 << std::endl;
  std::cout << "f_m12: " << f_m12 << std::endl;
  std::cout << "e_a02: " << e_a02 << std::endl;
  std::cout << "f_a02: " << f_a02 << std::endl;
  std::cout << "e_m12 * e_a02: " << e_m12 * e_a02 << std::endl;
  std::cout << "f_m12 * f_a02: " << f_m12 * f_a02 << std::endl;
  CGAL_assertion(CGAL::is_exact_in_interval(e_m012_3, f_m012_3));

//  std::cout << "---- sum_l" << std::endl;
//  const EK::FT e_m012_l = e_m01 * e_a22 - e_m02 * e_a12;
//  const FK::FT f_m012_l = f_m01 * f_a22 - f_m02 * f_a12;
//  CGAL_assertion(CGAL::is_exact_in_interval(e_m012_l, f_m012_l));

//  std::cout << "---- sum_r" << std::endl;
//  const EK::FT e_m012_r = - e_m02 * e_a12 + e_m12 * e_a02;
//  const FK::FT f_m012_r = - f_m02 * f_a12 + f_m12 * f_a02;
//  CGAL_assertion(CGAL::is_exact_in_interval(e_m012_r, f_m012_r));

//  std::cout << "---- m012" << std::endl;
//  const EK::FT e_m012 = e_m01 * e_a22 - e_m02 * e_a12 + e_m12 * e_a02;
//  const FK::FT f_m012 = f_m01 * f_a22 - f_m02 * f_a12 + f_m12 * f_a02;
//  CGAL_assertion(CGAL::is_exact_in_interval(e_m012, f_m012));
}
Cacher is the code snippet organizer for pro developers

We empower you and your team to get more done, faster

wrong_interval.cpp
