Minimal example of indeterminacy of CGAL::to_double ?

[alecjacobson]

I have a scenario where I would like the following to hold:

Given two CGAL::Epeck::FT numbers, a and b and two corresponding doubles c = cast_to_double(a) and d = cast_to_double(b), if a==b then I require that c==d .

I realize that the space of doubles is finite and there are (infinitely) many rationals who share the same closest float, so it's OK if cast_to_double()'s loss of precision implies the reverse (i.e., c==d ⇏ a==b) does not hold. Nevertheless, cast_to_double should be a good rounding of rational to float.

The first candidate for cast_to_double() is simply CGAL::to_double(), but in the documentation for CGAL::to_double there is a note which is giving me pause:

Note that in general, the value returned is not guaranteed to be the same when called several times on the same number. For example, if NT is a lazy number type (such as an instance of CGAL::Lazy_exact_nt), the double approximation returned might be affected by an exact computation internally triggered (that might have improved the double approximation).

Is there a minimal example which demonstrates the issue discussed in this note?

Does this imply that CGAL::to_double(x) will not meet my criteria above?

Would (the much slower) CGAL::to_double(x.exact())?

[alecjacobson]

Here's an attempt at a minimal example, but it'd be great to know if what's going on here is really matching that note.

I construct two equivalent numbers with disparate expressions a=1/3 and b=1-2/3. Calling CGAL::to_double() creates unequal doubles. Calling CGAL::to_double(.exact()) (I guess!) triggers some internal evaluations which improve a and b and now c==d. As hidden side-effects (const doesn't appear to make a difference here), those internal evaluations have been cached so now even just CGAL::to_double() will cause c==d.

Is my understanding correct?

Upon further scrutiny, I'm more confused.

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <cstdio>

int main()
{
  CGAL::Epeck::FT a = CGAL::Epeck::FT(1)/ CGAL::Epeck::FT(3);
  CGAL::Epeck::FT b = CGAL::Epeck::FT(1) - CGAL::Epeck::FT(2)/ CGAL::Epeck::FT(3);
  double c1,d1,c2,d2,c3,d3;
  c1 = CGAL::to_double(a);
  d1 = CGAL::to_double(b);
  printf("a%sb → c%sd\n",a==b?"=":"≠",c1==d1?"=":"≠");
  c2 = CGAL::to_double(a.exact());
  d2 = CGAL::to_double(b.exact());
  printf("a%sb → c%sd\n",a==b?"=":"≠",c2==d2?"=":"≠");
  c3 = CGAL::to_double(a);
  d3 = CGAL::to_double(b);
  printf("a%sb → c%sd\n",a==b?"=":"≠",c3==d3?"=":"≠");
  printf("c₁%sc₂\n",c1==c2?"=":"≠");
  printf("c₂%sc₃\n",c2==c3?"=":"≠");
  printf("c₁%sc₃\n",c1==c3?"=":"≠");
  printf("d₁%sd₂\n",d1==d2?"=":"≠");
  printf("d₂%sd₃\n",d2==d3?"=":"≠");
  printf("d₁%sd₃\n",d1==d3?"=":"≠");
}

this will output:

a=b → c≠d
a=b → c=d
a=b → c=d
c₁≠c₂
c₂≠c₃
c₁=c₃
d₁=d₂
d₂≠d₃
d₁≠d₃

A lingering confusion even if the understand above is correct: why is c₂≠c₃? If a's internal representation has improved. Why would it not reproduce the same double?

[afabri]

The Exact_predicates_exact_constructions_kernel has as number type Lazy_exact_nt which stores an approximation as Interval_nt and the arithmetic expression as a directed acyclic graph. For the first expression the root is a node storing the interval, the Division functor and two pointers to children for FT(1) and for FT(3). Similar for the other expression.
The intervals of a and b are different because the error propagation of interval arithmetic is different for IA(1) / IA(3) and IA(1) / (IA(1) - IA(2)).
The function to_double() produces the middle of the interval which explains that they differ for different intervals.

When you call a == b the intervals are compared. If the result can be computed just based on the intervals there is an early exit. If not, and that is the case here, first a.exact() and b.exact() are called. This evaluates the stored expression with arbitrary precision arithmetic, performs the comparision in arbitrary precision arithmetic, and updates the interval (which will be tighter now). So in the call

printf("a%sb → c%sd\n",a==b?"=":"≠",c1==d1?"=":"≠");

the values of c1 and d1 get computed before the a==b test, that is, for the numbers in question, the to_exact() in the next line is already done.

[MaelRL]

The second part of the comment above is specific to the actual example code that contains an explicit comparison a == b, which forces an exact computation. However, in general you hope to avoid doing this narrowing of the approximate representation (i.e., the intervals) for speed reasons. Hence you cannot guarantee a == b ==> to_double(a) == to_double(b) for random EPECK::FTs in the middle of your code because all you usually have are intervals that could be very different yet represent the same value, but you cannot know that before computing the exact value (which is what the code in cast_to_double in the comment from @lrineau below does).

[alecjacobson]

Thanks, Andreas and Mael! This is super helpful. When you write "produces the middle of the interval" does that mean that the interval is computed in double and then to_double returns something like (i0+i1)/2 ?

[MaelRL]

The bounds of the interval are always doubles, and the to_double() for the interval is (inf + sup) * 0.5.

[lrineau]

c2 = CGAL::to_double(a.exact());
c3 = CGAL::to_double(a);

a.exact() returns the exact type inside a (of type CGAL::Epeck::FT::ET), so c2 is constructed by calling CGAL::to_double(CGAL::Epeck::FT::ET), whereas c3 is constructed by calling CGAL::to_double(CGAL::Epeck::FT). That explains why c2 and c3 differ.

Your function cast_to_double could be implemented that way:

template <typename Lazy_FT>
double cast_to_double(const Lazy_FT& x) {
  x.exact();
  return CGAL::to_double(x);
}

[alecjacobson]

Excellent. This is precisely the reassuring answer I needed.