CGAL users discussion list

["Sebastien Loriot (GeometryFactory)" <>]

Sebastian Kapfer wrote:
> Hallo!
>
> > I was able to reproduce that on one of my machines, with your last code+data, 
> > with -O1, but not with -O0. It is probably a compiler bug. I do not understand 
> > it well.
>
> Weird.  I tried with GCC 4.4.3 which is the latest version I can get
> my hands at without compiling GCC from source.  It's still buggy.
>
> I also tried Intel C++ a minute ago:
>
> > $ icpc --version
> > icpc (ICC) 11.0 20090609
> > Copyright (C) 1985-2009 Intel Corporation.  All rights reserved.
>
> Same message, aborts in line 267, even with -O0 and -fp-model strict.
>
> Hm. :-(
The problem is due to floating point rounding.
Depending of the optimization done by the compiler an expression
was not evaluated the same way.
For convenience, I made a fix which solve the "problem" on my machine 
(patch file attached).

Also have a look at this FAQ entry:
http://www.cgal.org/FAQ.html#inexact_NT

S.
// Copyright (c) 2002 Utrecht University (The Netherlands).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// 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.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.7-branch/Spatial_searching/include/CGAL/Euclidean_distance.h $
// $Id: Euclidean_distance.h 37183 2007-03-17 09:21:06Z afabri $
// 
//
// Author(s)     : Hans Tangelder (<>)


#ifndef CGAL_EUCLIDEAN_DISTANCE_H
#define CGAL_EUCLIDEAN_DISTANCE_H
#include <CGAL/Kd_tree_rectangle.h>

namespace CGAL {



  template <class SearchTraits>
  class Euclidean_distance {

    public:

    typedef typename SearchTraits::FT    FT;
    typedef typename SearchTraits::Point_d Point_d;
    typedef Point_d Query_item;
    
    public:

    	// default constructor
    	Euclidean_distance() {}


	inline FT transformed_distance(const Query_item& q, const Point_d& p) const {
	        FT distance = FT(0);
		typename SearchTraits::Construct_cartesian_const_iterator_d construct_it;
                typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
		  qe = construct_it(q,1), pit = construct_it(p);
		for(; qit != qe; qit++, pit++){
		  distance += ((*qit)-(*pit))*((*qit)-(*pit));
		}
        	return distance;
	}


	inline FT min_distance_to_rectangle(const Query_item& q,
					    const Kd_tree_rectangle<SearchTraits>& r) const {
		FT distance = FT(0);
		typename SearchTraits::Construct_cartesian_const_iterator_d construct_it;
                typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
		  qe = construct_it(q,1);
		for(unsigned int i = 0;qit != qe; i++, qit++){
		  if((*qit) < r.min_coord(i))
				distance += 
				(r.min_coord(i)-(*qit))*(r.min_coord(i)-(*qit));
		  else if ((*qit) > r.max_coord(i))
				distance +=  
				((*qit)-r.max_coord(i))*((*qit)-r.max_coord(i));
			
		}
		return distance;
	}

	inline FT max_distance_to_rectangle(const Query_item& q,
					     const Kd_tree_rectangle<SearchTraits>& r) const {
		FT distance=FT(0);
		typename SearchTraits::Construct_cartesian_const_iterator_d construct_it;
                typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
		  qe = construct_it(q,1);
		for(unsigned int i = 0;qit != qe; i++, qit++){
				if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0))
					distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit));
				else
					distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i));
		};
		return distance;
	}

	inline FT new_distance(FT dist, FT old_off, FT new_off,
			       int /* cutting_dimension */)  const {
		
		FT new_dist = dist + (new_off*new_off - old_off*old_off);
                return new_dist;
	}

  inline FT transformed_distance(FT d) const {
		return d*d;
	}

  inline FT inverse_of_transformed_distance(FT d) const {
		return CGAL::sqrt(d);
	}

  }; // class Euclidean_distance

} // namespace CGAL
#endif // EUCLIDEAN_DISTANCE_H