FPU: Support platforms without round-towards-infinity

[mglisse]

The platform I have in mind is WebAssembly, which does not support non-default rounding modes. This platform may still be worth targeting, for small examples or demos where performance is not that important (or for applications where static filters are sufficient).

• First idea would be a software emulation of round-towards-infinity, or more roughly replace a +up b with nextafter(a +near b, +infinity). Looks doable, nextafter is not so cheap but probably still better than full rationals.
• Second idea is to disable intervals. Note that static filters are still ok (they use round-to-nearest), so we would likely want Epeck to be Simple_cartesian with static filters on top, and Epick predicates would try static filters and then go directly to the exact computation.
• I am not sure how well Boost supports it for now, which may cause other issues. GMP/MPFR also require unofficial patches IIRC, although the branch Number_types-boost_mp-glisse makes them not-so-necessary.

[lrineau]

@mglisse, @horasio Will you be at next developers meeting? That could be an opportunity to start and work on that topic.

[mglisse]

I am unlikely to be at the next meeting.

[sloriot]

@mglisse did you already start a branch for option 1, I'd be curious to see the slow down in practice.

[mglisse]

@sloriot no, I haven't experimented with anything like that. For a quick try, I guess you could just redefine CGAL_IA_ADD, etc.

[lrineau]

LLVM 8, just released, has now an official WebAssembly target:
http://releases.llvm.org/8.0.0/docs/ReleaseNotes.html#changes-to-the-webassembly-target

[pentacular]

I'm taking a look at option 1.

But I'm not sure I understand how Interval_nt.h is supposed to be working.

I see

  compare (const Interval_nt<Protected> & d, const Interval_nt<Protected> & e)
  {
    if (d.inf() > e.sup()) return LARGER;
    if (e.inf() > d.sup()) return SMALLER;
    if (e.inf() == d.sup() && d.inf() == e.sup()) return EQUAL;
    return Uncertain<Comparison_result>::indeterminate();
  }

These looks safe for round-to-nearest to me (please tell me if I am wrong)

    if (d.inf() > e.sup()) return LARGER;
    if (e.inf() > d.sup()) return SMALLER;

But I don't understand

    if (e.inf() == d.sup() && d.inf() == e.sup()) return EQUAL;

In is_same I see:

    return inf() == d.inf() && sup() == d.sup();

which does make sense -- why are we comparing inf with sub in compare, but inf with inf and sub with sub in is_same?

[sloriot]

Basically, the rounding mode only matters when you update the interval with operator+-/*.

[sloriot]

sqrt() and co. too btw

[mglisse]

We have Internal_protector where the rounding mode matters. Comparing inf and sup: this is an efficient way to say both is_same and is_point, we can only know for sure that 2 numbers are equal if their intervals are reduced to a point, and they are the same. Conceptually, Interval_nt is not meant to represent an interval, but to represent a number that we do not know exactly, we only know an interval containing it. And we try to go as far as possible with that information, throwing an exception when we don't know what to do (say if we try to branch on whether a number in [-2,3] is positive).

For option 1, I think you want to look at macros like CGAL_IA_ADD since that's the operation that needs rounding.

[pentacular]

Ok, that makes sense -- I was worried that I'd misunderstood inf and sup.

Looking at + I see

      return Interval_nt (-CGAL_IA_ADD(-a.inf(), -b.inf()), CGAL_IA_ADD(a.sup(), b.sup()));

Based on

// Note: When rounding is towards +infinity, to make an operation rounded
// towards -infinity, it's enough to take the opposite of some of the operand,
// and the opposite of the result (see operator+, operator*,...).

-CGAL_IA_ADD(-a.inf(), -b.inf()) rounding toward -inf

So this is really expressing

    return Interval_nt(round_down(a.inf() + b.inf()), round_up(a.sup() + b.sup()));

Which makes sense as it gives a conservative interval around the result.

But I think it also means that it needs to be fixed at a higher level than CGAL_IA_ADD, since CGAL_IA_ADD doesn't know which way you want to round the result.

Does that sound right to you?

[mglisse]

CGAL_IA_ADD doesn't know which way you want to round the result

Actually, I believe CGAL_IA_ADD always needs to round towards +infinity. And then we use tricks like -CGAL_IA_ADD(-x,-y) to compute an addition rounded towards -infinity without changing the rounding mode.
But if it helps, we could add some CGAL_IA_ADD_UP macro (actually, we should replace these macros with inline functions), or add a comment in FPU.h clarifying things.

[pentacular]

I think the key here is that -OP(-a, -b) only works with round towards +infinity, which we don't support here.
So we need to provide an alternative to that trick.

I think that CGAL_IA_ADD_UP and CGAL_IA_ADD_DOWN would be a suitable abstraction for this case.
Likewise CGAL_IA_MUL_UP and CGAL_IA_MUL_DOWN, etc.

I agree that inline functions would probably be a better choice than macros.

I think that I was a bit off earlier, and it should effectively be

    return Interval_nt(down(down(a.inf()) + down(b.inf())), up(up(a.sup()) + up(b.sup())));

Although maybe we can get away with

    return Interval_nt(down(a.inf()) + down(b.inf()), up(a.sup()) + up(b.sup()));

What do you think?

[mglisse]

CGAL_IA_ADD_DOWN(a,b) and -CGAL_IA_ADD_UP(-a,-b) still look similar to me.
If by up and down you mean a corresponding call to nextafter, then up(a.sup()+b.sup()) is right, it gives a number larger than the sum computed with infinite precision. The idea is to fake rounding towards +infinity using nextafter.

[pentacular]

Yes.

I see what you mean, but don't you want up(a.sup) + up(b.sup())?

We need the rounding to apply before the operation in order to get a conservative result, right?

[mglisse]

a.sup is an upper bound on a, b.sup is an upper bound on b, computing a.sup+b.sup with upwards rounding is an upper bound on a+b. Now, if n is x+y computed with the default to-nearest rounding and u is x+y computed with upwards rounding, either u is n or u is nextafter(n,+inf), and in all cases nextafter(n,+inf) is an upperbound on u.

[pentacular]

Let's imagine that both a and b were rounded down toward nearest, could we undershoot the result?

e.g., (in an imaginary system) 10.1 + 10.1, rounded down to 10 + 10, giving 20 which we round up to 20.1, but ought to actually be 20.3?

I think in general, at least it needs to be up(a) OP up(b) , at least with MUL and DIV -- but perhaps we could get away with it for ADD if we can be sure about it ...

[mglisse]

There is no such thing as "a and b were rounded down toward nearest". Let me use capital letters for the mathematical exact value (infinite precision). The invariant is that a.inf() <= A <= a.sup(). All we want is to preserve this invariant in operations. If A is 10.1, then there is no way a.sup() can be 10.

[pentacular]

Ok, providing a.inf() <= A <= a.sup() is initially true, then it should remain that way.

Do we need to fix the constructor to be conservative?
Otherwise a.inf() <= A <= a.sup() might not be really true.

Or is it sufficient that <= for double on that system is consistent about it?

[mglisse]

Do we need to fix the constructor to be conservative?

They already are. The original value is assumed to be exact, whether it is an int, double, Gmpq or whatnot, and the conversion produces an interval that contains the original value. For int or double, an interval with inf==sup is sufficient. For Gmpq, the interval may be larger.

[pentacular]

Hmm, I don't see a constructor for Interval_nt from Gmpq or any exact non-integer value.
I guess the conversion happens earlier in Filtered_predicate in C2A, but it's not obvious to me where we ensure it precisely contains the value, for a given rounding mode.

I've started testing the #define CGAL_IA_ADD(a,b) CGAL_IA_UP(CGAL_IA_FORCE_TO_DOUBLE((a)+CGAL_IA_STOP_CPROP(b))) approach, but it's not straight-forward to isolate the interval FPU requirements for testing.

What I've done so far is to split out a variant of Protect_FPU_rounding which enforces CGAL_FE_TONEAREST and have Interval_nt use this instead, and to force always setting and backing up mode to work around assumptions of inheriting the expected mode.

I've also been trying to use values like IA(-CGAL_IA_MAX_DOUBLE, CGAL_IA_MAX_DOUBLE) to selectively opt out of the approximation for particular operations, like div, sqrt, etc.

Can you let me know if that sounds reasonable?

I've been getting some exceptions still, so either I'm doing something fundamentally wrong, or I've missed something somewhere.

[mglisse]

Hmm, I don't see a constructor for Interval_nt from Gmpq or any exact non-integer value.

Look for To_interval in files for those types (Gmpq.h, mpq_class.h, etc).

[mglisse]

I tried this hack to see what happens

diff --git a/Number_types/include/CGAL/FPU.h b/Number_types/include/CGAL/FPU.h
index 6d5f02a3075..693a92279dd 100644
--- a/Number_types/include/CGAL/FPU.h
+++ b/Number_types/include/CGAL/FPU.h
@@ -117,7 +117,7 @@ extern "C" {
 #if defined(CGAL_HAS_SSE2) && (defined(__x86_64__) || defined(_M_X64))
 #  define CGAL_USE_SSE2 1
 #endif
-
+#undef CGAL_USE_SSE2
 #ifdef CGAL_CFG_DENORMALS_COMPILE_BUG
 double& get_static_minimin(); // Defined in Interval_arithmetic_impl.h
 #endif
@@ -347,16 +347,17 @@ inline double IA_bug_sqrt(double d)
 // With GCC, we can do slightly better : test with __builtin_constant_p()
 // that both arguments are constant before stopping one of them.
 // Use inline functions instead ?
-#define CGAL_IA_ADD(a,b) CGAL_IA_FORCE_TO_DOUBLE((a)+CGAL_IA_STOP_CPROP(b))
-#define CGAL_IA_SUB(a,b) CGAL_IA_FORCE_TO_DOUBLE(CGAL_IA_STOP_CPROP(a)-(b))
-#define CGAL_IA_MUL(a,b) CGAL_IA_FORCE_TO_DOUBLE(CGAL_IA_STOP_CPROP(a)*CGAL_IA_STOP_CPROP(b))
-#define CGAL_IA_DIV(a,b) CGAL_IA_FORCE_TO_DOUBLE(CGAL_IA_STOP_CPROP(a)/CGAL_IA_STOP_CPROP(b))
+inline double CGAL_IA_UP(double d){return nextafter(d,std::numeric_limits<double>::infinity());}
+#define CGAL_IA_ADD(a,b) CGAL_IA_UP((a)+CGAL_IA_STOP_CPROP(b))
+#define CGAL_IA_SUB(a,b) CGAL_IA_UP(CGAL_IA_STOP_CPROP(a)-(b))
+#define CGAL_IA_MUL(a,b) CGAL_IA_UP(CGAL_IA_STOP_CPROP(a)*CGAL_IA_STOP_CPROP(b))
+#define CGAL_IA_DIV(a,b) CGAL_IA_UP(CGAL_IA_STOP_CPROP(a)/CGAL_IA_STOP_CPROP(b))
 inline double CGAL_IA_SQUARE(double a){
   double b = CGAL_IA_STOP_CPROP(a); // only once
-  return CGAL_IA_FORCE_TO_DOUBLE(b*b);
+  return CGAL_IA_UP(b*b);
 }
 #define CGAL_IA_SQRT(a) \
-        CGAL_IA_FORCE_TO_DOUBLE(CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(a)))
+        CGAL_IA_UP(CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(a)))
 
 
 #if defined CGAL_SAFE_SSE2
@@ -476,6 +477,7 @@ inline
 FPU_CW_t
 FPU_get_cw (void)
 {
+    return CGAL_FE_TONEAREST;
     FPU_CW_t cw;
     CGAL_IA_GETFPCW(cw);
     return cw;
@@ -487,6 +489,7 @@ inline
 void
 FPU_set_cw (FPU_CW_t cw)
 {
+  return;
   CGAL_IA_SETFPCW(cw);
 }
 
@@ -494,6 +497,7 @@ inline
 FPU_CW_t
 FPU_get_and_set_cw (FPU_CW_t cw)
 {
+    return CGAL_FE_TONEAREST;
     FPU_CW_t old = FPU_get_cw();
     FPU_set_cw(cw);
     return old;
diff --git a/Number_types/include/CGAL/Interval_nt.h b/Number_types/include/CGAL/Interval_nt.h
index 8a51e7cd860..80e0d0b07e5 100644
--- a/Number_types/include/CGAL/Interval_nt.h
+++ b/Number_types/include/CGAL/Interval_nt.h
@@ -1148,7 +1148,7 @@ namespace INTERN_INTERVAL_NT {
     // - compute y = CGAL_IA_SQUARE(x)
     // - if y==d.inf() use x, else use -CGAL_IA_SUB(CGAL_IA_MIN_DOUBLE,x)
     FPU_set_cw(CGAL_FE_DOWNWARD);
-    double i = (d.inf() > 0.0) ? CGAL_IA_SQRT(d.inf()) : 0.0;
+    double i = (d.inf() > 0.0) ? nextafter(std::sqrt(d.inf()), 0.) : 0.0;
     FPU_set_cw(CGAL_FE_UPWARD);
 #endif
     return Interval_nt<Protected>(i, CGAL_IA_SQRT(d.sup()));

The tests Interval_nt(_new).cpp fail in several places that assume that the computation of intervals is as tight as possible, which is obviously not going to be the case here, but that's something to handle in the tests. The precise tests need to be skipped when using sloppy intervals, and new relaxed tests need to be implemented.

[pentacular]

Yes, that's pretty much what I have -- thanks for the reference. :)

However, I find that when I try to run corefinement_mesh_union.cpp with an Epeck kernel with these modifications, it starts to fail with

terminate called after throwing an instance of 'CGAL::Assertion_exception'
  what():  CGAL ERROR: assertion violation!
Expr: s != LEFT_TURN
File: ./CGAL/Triangulation_2.h
Line: 933
Aborted

I've added CGAL_ALWAYS_ROUND_TO_NEAREST to control the modifications, and with this set false it works as expected -- so the basic kernel configuration seems sound.

Which leads me to two theories.

(a) Rounding mode requirements outside of Interval_nt are being violated, or
(b) There is an error in the operations within Interval_nt.

To test (a) I've built out a second Protect_FPU_Rounding_2 class, and forced them to always set and reset on each operation.

This should mean that the rest of the system continues rounding as per usual, and can be incrementally opted in to CGAL_ALWAYS_ROUND_TO_NEAREST.

To test (b) I've changed the Interval_nt operations to return return Interval_nt<Protected>(-CGAL_IA_MAX_DOUBLE, CGAL_IA_MAX_DOUBLE);

I hope this means that it falls back to the exact kernel for these operations, allowing us to effectively turn them on and off.

Does my understanding in (b) seem correct?

[pentacular]

Having written that out, I realized that I can more directly test the Interval_nt ops by simply forcing the rounding mode inside each.

This let me bisect the problem down to here:

#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
      auto r1 = -CGAL_IA_DIV(-a.inf(), aa);
      FPU_set_cw(CGAL_FE_UPWARD);
      auto r2 = CGAL_IA_DIV(a.sup(), bb); // <- here
      FPU_set_cw(CGAL_FE_TONEAREST);
      return IA(r1, r2);
#else
      return IA(-CGAL_IA_DIV(-a.inf(), aa), CGAL_IA_DIV(a.sup(), bb));
#endif

So I guess CGAL_IA_DIV requires a bit more thought, although I'm not sure why the issue doesn't manifest on the other CGAL_IA_DIV paths -- perhaps they're just not sufficiently exercised.

I'll take a closer look at it tomorrow -- but thought you might have some insight. :)

[mglisse]

"corefinement_mesh_union.cpp with an Epeck kernel": that doesn't even compile on master. It would be easier to understand your comments if they referred to some public branch, showed a command line, etc.

[pentacular]

Hmm, I can't see any reason that corefinement ought not work with Epeck -- isn't this essential for being able to produce non-self-intersecting meshes for incremental operation?

Here is what I am using.

#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
// #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/corefinement.h>
#include <fstream>
// typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Exact_predicates_exact_constructions_kernel K;
typedef CGAL::Surface_mesh<K::Point_3>             Mesh;
namespace PMP = CGAL::Polygon_mesh_processing;
int main(int argc, char* argv[])
{
  const char* filename1 = (argc > 1) ? argv[1] : "data/blobby.off";
  const char* filename2 = (argc > 2) ? argv[2] : "data/eight.off";
  std::ifstream input(filename1);
  Mesh mesh1, mesh2;
  if (!input || !(input >> mesh1))
  {
    std::cerr << "First mesh is not a valid off file." << std::endl;
    return 1;
  }
  input.close();
  input.open(filename2);
  if (!input || !(input >> mesh2))
  {
    std::cerr << "Second mesh is not a valid off file." << std::endl;
    return 1;
  }
  Mesh out;
  bool valid_union = PMP::corefine_and_compute_union(mesh1,mesh2, out);
  if (valid_union)
  {
    std::cout << "Union was successfully computed\n";
    std::ofstream output("union.off");
    output.precision(17);
    output << out;
    return 0;
  }
  std::cout << "Union could not be computed\n";
  return 1;
}

The build environment I have for wasm is not an easy fit to how cgal is normally set up, so it's not easy to share at the moment.

--

Having investigated a bit further, it turns out that in some case cases a / b needs to be bumped upward by two steps in round-to-nearest to match what happens in round-upward.

Having bumped the result up twice in this case, the example successfully completes.

I'm not sure why this is so, but it should be straight-forward to investigate now that we have particular values.

If you're interested, one set of example values is:

a.sup: 3.2991353107230831e-06
bb: 3.9193289838502824e-05

      auto r1 = -CGAL_IA_DIV(-a.inf(), aa);
      auto r2a = CGAL_IA_DIV(a.sup(), bb);
      auto r2au = CGAL_IA_UP(r2a);
      FPU_set_cw(CGAL_FE_UPWARD); 
      auto r2b = CGAL_IA_DIV(a.sup(), bb);
      FPU_set_cw(CGAL_FE_TONEAREST);

r2a: 0.084176024118344586 // round-to-nearest bumped up once
r2a2: 0.0841760241183446 // round-to-nearest bumped up twice
r2b: 0.0841760241183446 // round-up

Perhaps it is due to something like excess precision.

--

What I still don't understand is why returning largest() from INTERN_INTERVAL_NT's / operator doesn't cause it to fall back gracefully to the exact operation.

Would you expect returning largest() to cause it to fall back to the exact operation, or am I confused about how intervals are supposed to work?

[mglisse]

Hmm, I can't see any reason that corefinement ought not work with Epeck

When I tried yesterday, it failed because of a call to sqrt, which doesn't work with Epeck.

[mglisse]

Here is what I am using.

This test passes here with the patch I posted.

[lrineau]

The build environment I have for wasm is not an easy fit to how cgal is normally set up, so it's not easy to share at the moment.

Before you start running that modified CGAL code in the WASM environment, you should first debug it in the usual host environment (probably using an x86_64 processor).

[pentacular]

Good advice. :)

I rebuilt CGAL and the dependencies and made sure that the compilation options were aligned between the x86 and wasm builds.
I am not entirely sure why, but the differences I was seeing under x86 disappeared.

Having done that, I still had issues under wasm but was able to produce readable backtraces in the browser and track them down.

I am still testing the results, but they seem plausible so far.
I'll clean up the code and tests and prepare a PR.

There will need to be some additional work to eliminate the dependency on catching Uncertain_conversion_exception in Filtered_predicate.h, as this produces major overhead, but it is not necessary for correctness and can be handled independently.

Thanks for all of your help -- CGAL internals made a good deal more sense to me now.

[mglisse]

eliminate the dependency on catching Uncertain_conversion_exception in Filtered_predicate.h, as this produces major overhead, but it is not necessary for correctness and can be handled independently.

Er, catching this exception there is a very important part of CGAL's filtering strategy, getting rid of it doesn't seem trivial. I understand that c++ exceptions are still a second class feature in wasm (https://github.com/WebAssembly/exception-handling/ indicates this may change eventually), but it looks like emscripten can handle them, at some cost.

[pentacular]

Yes -- it would require a different mechanism to coordinate the filters.

The cost is pretty huge -- the wasm becomes about 160% of the size (3.2M to 5.2M) and the latency is only about 65% that with simple gmpq.

Anyhow, we can worry about that independently of this. :)

[pentacular]

I've put together a PR, but suspect it needs restructuring.

In particular I'm not sure what branch I should be targeting, and I'm not sure about the test organization.

[pentacular]

I'm re-opening this issue as I've found a possible problem.

Using the following test program (under Surface_mesh/test/Surface_mesh) I've found some unexpected behavior with CGAL_ALWAYS_ROUND_TO_NEAREST.

#define CGAL_ALWAYS_ROUND_TO_NEAREST
  
#include <CGAL/Surface_mesh/Surface_mesh.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Polygon_mesh_processing/transform.h>

#include <cstring>
#include <iostream>
#include <fstream>

typedef CGAL::Exact_predicates_exact_constructions_kernel     Kernel;
typedef Kernel::Point_3                                       Point;
typedef Kernel::Aff_transformation_3                          Transformation;

typedef CGAL::Surface_mesh<Point>                             SMesh;
typedef boost::graph_traits<SMesh>::vertex_descriptor         vertex_descriptor;
typedef boost::graph_traits<SMesh>::face_descriptor           face_descriptor;

int main()
{
  std::ifstream in("sm_off_io_input.off");
  SMesh mesh;
  CGAL::read_OFF(in, mesh);

  std::cerr << "before: " << std::setprecision(20) << mesh << std::endl;

  Transformation identity(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0);
  CGAL::Polygon_mesh_processing::transform(identity, mesh, CGAL::parameters::all_default());

  std::cerr << "after: " << std::setprecision(20) << mesh << std::endl;

  return 0;
}

This produces the following output for me with CGAL_ALWAYS_ROUND_TO_NEAREST true.

# Output of a CGAL tool
OFF
8 12 0

# 8 vertices
# ------------------------------------------


-0.9000000000000000222 -0.9000000000000000222 -0.9000000000000000222
-1.000000000000000222 1.000000000000000222 -1.000000000000000222
1.000000000000000222 1.000000000000000222 -1.000000000000000222
1.000000000000000222 -1.000000000000000222 -1.000000000000000222
-1.000000000000000222 -1.000000000000000222 1.000000000000000222
-1.000000000000000222 1.000000000000000222 1.000000000000000222
1.000000000000000222 1.000000000000000222 1.000000000000000222
1.000000000000000222 -1.000000000000000222 1.000000000000000222

# 12 facets
# ------------------------------------------

3  0 1 3
3  3 1 2
3  0 4 1
3  1 4 5
3  3 2 7
3  7 2 6
3  4 0 3
3  7 4 3
3  6 4 7
3  6 5 4
3  1 5 6
3  2 1 6


# End of OFF #

But with CGAL_ALWAYS_ROUND_TO_NEAREST false the values are unchanged.

This shouldn't be happening, right? :)

[lrineau]

What is in the file sm_off_io_input.off?

Note that you use an ASCII precision of 20, but double have between 15 or 16 digits of precision (53 bits of precision). For example, those two values, as double, are actually equal:

• 1.000000000000000222, and
• 1.0.

Probably what you see is a strange effect of the algorithm that converts from binary to ASCII for floating point values, with ROUND_TO_NEAREST.

[lrineau]

For double, std::numeric_limits<double>::max_digits10 is equal to 17. That is the precision that ensures ASCII inputs/outputs for double values are lossless.

[pentacular]

The effect remains with setprecision(17).

Surface_mesh/test/Surface_mesh/sm_off_io_input.off used here contains

OFF
8 12 0
-0.9 -0.9 -0.9
-1 1 -1
1 1 -1
1 -1 -1
-1 -1 1
-1 1 1
1 1 1
1 -1 1
3  0 1 3
3  3 1 2
3  0 4 1
3  1 4 5
3  3 2 7
3  7 2 6
3  4 0 3
3  7 4 3
3  6 4 7
3  6 5 4
3  1 5 6
3  2 1 6

So I think this may be a real deviation.

[mglisse]

Sounds like what happens when we convert an interval (from Lazy_exact_nt) to double?

[lrineau]

Sounds like what happens when we convert an interval (from Lazy_exact_nt) to double?

... You are right! The kernel is CGAL::Exact_predicates_exact_constructions_kernel, so in memory the coordinates are not double.

[pentacular]

Ok, so is the problem that it's not calling exact() before converting to double?

[pentacular]

Ok, so I've confirmed that calling exact() before converting to double addresses the problem.

I think this means that it is due to converting an interval to a double, and this is not an issue with CGAL_ALWAYS_ROUND_TO_NEAREST.

Thanks.