`to_double` for non-trivial expressions that have an exact double evaluation

qwenger commented 3 months ago

Issue Details

Consider the program below. I get the following output running it:

direct: 1.0000000000000000 1.0000000000000000
double: 0.99999999999999989 1
direct: 2 0
double: 2 0
direct: 0 2
double: 0 2
direct: 2 2
double: 2 2
direct: 0 0
double: 0 0

The first vertex should be at the exact square center, (1,1).

I understand the limitations of converting to floating point values, but here the correct value does have an exact double representation. I also understand the difficulty of evaluating an expression, but the code for printing the value directly does a fine (or at least, better) job.

So:

Why doesn't to_double return the double value that is actually closest to the exact value?
Why does direct printing do better?
Is there a way to robustly get the "correct" double value in such cases? (other than stringifying/parsing...)

I've read https://github.com/CGAL/cgal/discussions/6000 but it doesn't seem to apply here, as CORE::Expr doesn't have an exact method.

The doc for to_double writes

Remark: In order to control the quality of approximation one has to resort to methods that are specific to NT. There are no general guarantees whatsoever.

What would then be the appropriate "specific methods" for CORE::Expr?

Source Code

#include <iostream>

#include <CGAL/Exact_predicates_exact_constructions_kernel_with_sqrt.h>
#include <CGAL/Segment_Delaunay_graph_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_policies_2.h>
#include <CGAL/Segment_Delaunay_graph_adaptation_traits_2.h>
#include <CGAL/Segment_Delaunay_graph_traits_2.h>
#include <CGAL/Voronoi_diagram_2.h>

using K = CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt;
using Point = CGAL::Point_2<K>;
using Gt = CGAL::Segment_Delaunay_graph_traits_2<K>;
using SDG2 = CGAL::Segment_Delaunay_graph_2<Gt>;
using AT = CGAL::Segment_Delaunay_graph_adaptation_traits_2<SDG2>;
using AP = CGAL::Segment_Delaunay_graph_degeneracy_removal_policy_2<SDG2>;
using VoronoiDiagram = CGAL::Voronoi_diagram_2<SDG2, AT, AP>;
using Site = AT::Site_2;

int main() {
    VoronoiDiagram vd;
    vd.insert(Site::construct_site_2({0.0, 0.0}, {2.0, 0.0}));
    vd.insert(Site::construct_site_2({2.0, 0.0}, {2.0, 2.0}));
    vd.insert(Site::construct_site_2({2.0, 2.0}, {0.0, 2.0}));
    vd.insert(Site::construct_site_2({0.0, 2.0}, {0.0, 0.0}));
    for (auto vit = vd.vertices_begin(); vit != vd.vertices_end(); ++vit) {
        std::cout << "direct: " << std::setprecision(17) << vit->point().x() << " " << vit->point().y() << std::endl;
        std::cout << "double: " << std::setprecision(17) << CGAL::to_double(vit->point().x()) << " " << CGAL::to_double(vit->point().y()) << std::endl;
    }
    return 0;
}

Environment

Operating system (Windows/Mac/Linux, 32/64 bits): Windows/Linux, 64 bits
Compiler: MSVC/GCC
Release or debug mode: debug
Specific flags used (if any):
CGAL version: 5.6.1
Boost version: 1.83.0
Other libraries versions if used (Eigen, TBB, etc.):

sloriot commented 3 months ago

In https://github.com/CGAL/cgal/pull/8179 with @afabri , we modified the behavior of the outstream operator for Core Expr and now you would get consistent results:

direct: 0.99999999999999989 1
double: 0.99999999999999989 1

Note that if you ask for more than 17 for the precision, you'll get the old behavior. Even if the outstream operator prints with more precision, if you want to go back into the world of double you will lose that precision. Maybe that would need internally to refine the interval around the real value but I'm not sure if that's already an option available.

qwenger commented 3 months ago

Thanks for the answer. So in the future it will be consistent, albeit exactly in the opposite way as what I would like.

if you want to go back into the world of double you will lose that precision

I'm not sure that I completely get this point. Yes in general doubles have limited precision, that is there are only finitely many values that can be represented exactly. But here the exact value is one of them, so doubles are not an excuse per se. It's the process of evaluating an approximation of the expression that returns the "wrong" double.

So, is there a way to get the "correct" one, even if it means sacrificing some performance?

And more generally, what does it mean to have "17 for the precision", when the 16th digit is incorrect? What are the precision guarantees, or is it really that there are no guarantees at all? If so, what would be the point of offering a conversion to double when the result could be arbitrarily far away from the truth? Or is it at least possible to find guarantees for the specific application at hand (Voronoi vertices positions)?

CGAL / cgal