Round-to-nearest support for intervals #5510

pentacular · 2021-03-03T11:32:08Z

Summary of Changes

This adds support for using round-to-nearest for intervals, as an alternative to round-to-infinity, suitable for use with WASM.

Note: I'm a bit unclear on how this ought to be best organized for CGAL -- please let me know what would fit better.

Release Management

Affected package(s): Number_types
Issue(s) solved (if any): fix FPU: Support platforms without round-towards-infinity #2791

sloriot · 2021-03-03T11:56:37Z

Thanks for the PR, this is something that was on my TODO list but I did not find time to do it.

We should add some cmake machinery so that the there is an option to allow to not change the rounding mode + do not add the -frounding-math for g++ when ON.

sloriot · 2021-03-03T12:02:20Z

Number_types/include/CGAL/FPU.h

-#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))
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST


IMO, if we are adding an extra case, we should simply assume that we do not know what the rounding mode is and not rely on to-nearest. I'm in favor of something like CGAL_NO_FPU_ROUNDING_MODE_CHANGE

I guess, this implies more changes but that should be the way to go.

It's an interesting point.

I think it's difficult to predict correct behavior with all existing and future rounding modes.

So I'd avoid a claim like CGAL_NO_FPU_ROUNDING_MODE_CHANGE.

But something like CGAL_CONSERVATIVE_INTERVALS might describe what's actually happening in this change, and then you could decide if that covers whatever floating point modes you care about or not.

I also notice that there are systems other than interval which care about floating point mode, and it's not clear to me that there will be a generic fix for all of them.

In which case CGAL_ALWAYS_ROUND_TO_NEAREST would be the right kind of flag to tell them about the particular environment they're operating in -- maybe some would need a CGAL_ALWAYS_ROUND_UPWARD, etc

I don't think so. I think as long as you use nextafter() for the upper and -lower bound of the interval you are safe side no matter how the value was rounded

For intervals, certainly -- but CGAL_NO_FPU_ROUNDING_MODE_CHANGE seems to make a more universal claim that everything in CGAL will always work without changing rounding mode.

AFAIK, this is the only place where is is needed.

Static filters do require round-to-nearest AFAIK (the bounds would be worse with arbitrary rounding), not sure if that's what you are talking about.

Hmm, I see a number of cases that aren't clear to me.

Bounding_volumes/include/CGAL/Approximate_min_ellipsoid_d.h -- toward-zero and upward
Apollonius_graph_2/include/CGAL/Apollonius_graph_2/compare_quadratic.h -- upward
Number_types/include/CGAL/leda_bigfloat.h -- upward
Mesh_3/include/CGAL/Mesh_3/Robust_intersection_traits_3.h -- to nearest
Skin_surface_3/include/CGAL/Skin_surface_base_3.h -- to nearest
Straight_skeleton_2/include/CGAL/Straight_skeleton_2/Straight_skeleton_builder_traits_2_aux.h -- to nearest
Polynomial/test/Polynomial/Polynomial_traits_d_Residue.cpp -- upward
Polynomial/test/Polynomial/Polynomial_using_core.cpp -- upward
Polynomial/test/Polynomial/resultant.cpp -- upward
Polynomial/test/Polynomial/Polynomial_using_leda.cpp -- upward
Convex_hull_3/include/CGAL/convex_hull_3.h

Some of these don't seem to be part of filters, e.g., Approximate_min_ellipsoid.

My suggestion is to require explicit opt-in to any forced rounding mode.

i.e. If we use CGAL_ALWAYS_ROUND_TO_NEAREST then we can have CGAL_IA_SETFPCW throw an exception (or a compile-time error) when something tries to use a different mode, which will require systems to explicitly opt into CGAL_ALWAYS_AROUND_TO_NEAREST.

This should minimize surprising behavior, and allow gradual adoption, and allow alternate implementation for particular modes.

On the other hand, CGAL_NO_FPU_ROUNDING_MODE_CHANGE seems like it might require code to dynamically dispatch on the current rounding mode in the worst case, since it would need to operate in any existing or future rounding mode.

In any case, I'm operating from a position of ignorance, so maybe I'm being overly conservative -- but please make sure you're confident that you've considered all of the possibilities. :)

sloriot · 2021-03-03T12:03:24Z

Number_types/include/CGAL/FPU.h

 template <bool Protected = true> struct Protect_FPU_rounding;

 template <>
 struct Protect_FPU_rounding<true>
 {
-  Protect_FPU_rounding(FPU_CW_t r = CGAL_FE_UPWARD)
+  Protect_FPU_rounding(FPU_CW_t r = CGAL_FE_PROTECTED)


The protector should not do anything in the new mode. The partial specialization for the protected mode could be hidden behind the macro.

I think we should consider testing that we are using the expected rounding mode, and error otherwise.

mglisse · 2021-03-03T12:11:04Z

Number_types/include/CGAL/Interval_nt.h

    // TODO: Alternative for computing CGAL_IA_SQRT_DOWN(d.inf()) exactly
    // without changing the rounding mode:
    // - compute x = CGAL_IA_SQRT(d.inf())
    // - compute y = CGAL_IA_SQUARE(x)
    // - if y==d.inf() use x, else use -CGAL_IA_SUB(CGAL_IA_MIN_DOUBLE,x)
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST


#elif is usually nicer than nested #if.

I'm happy to change it, but was concerned that the context of the comment about the alternative approach would become confused.

I think this is obsoleted by the abstraction into IA_sqrt_toward_zero.

mglisse · 2021-03-03T12:14:09Z

Number_types/include/CGAL/FPU.h

@@ -347,16 +346,32 @@ 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))


You dropped CGAL_IA_FORCE_TO_DOUBLE for the normal case, that looks wrong.

Yes, I think you're right -- I didn't look deep enough into all of the cases it handles.

How about putting that into the other branch of CGAL_IA_UP?

Or is that getting too overloaded?

Putting it in the other branch of CGAL_IA_UP sounds fine.

mglisse · 2021-03-03T21:46:14Z

Number_types/include/CGAL/FPU.h

-#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))
+#ifdef CGAL_ALWAYS_ROUND_TO_NEAREST
+inline double CGAL_IA_UP(double d)


Since it isn't a macro and it is already in namespace CGAL, you don't need to prefix it with CGAL and write it in all caps. I see functions called IA_opacify or IA_force_to_double nearby.

Thanks -- I've changed to IA_up.

Please let me know if you have a better name.

mglisse · 2021-03-03T21:50:06Z

Number_types/include/CGAL/FPU.h

 }
 #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)))


It may be slightly confusing that in the new mode CGAL_IA_SQRT means up(sqrt) while in the old mode, it means sqrt with dynamic rounding, and it is actually used once with downward rounding (the other operations are only ever used with upward rounding).

I've abstracted IA_sqrt_up and IA_sqrt_toward_zero.

It's a little awkward because of scoping.

Let me know what you think.

mglisse · 2021-03-04T23:20:56Z

You mentioned benchmarks in the issue, could you post some benchmark results, explaining exactly what you tested? It would be interesting to know how this compares to current intervals, and how it compares to not using intervals at all.

pentacular · 2021-03-05T00:28:48Z

I just did a simple test which produced a sphere, then subtracted a half size sphere incrementally around the equator at each 360 / 7 degrees. https://imgur.com/a/ucZu3cl

The time taken with this approach was about 65% of the time taken without intervals, which is a significant improvement, but not as much as hoped.

It looks like the overhead of supporting exceptions was the main reason that it wasn't significantly faster.

I looked into the possibility of replacing exceptions within CGAL, but it seems like it may be infeasible.

I was concerned about this, but it looks like llvm has added support for native wasm exceptions, and these should become available via emscripten in a few months -- hopefully with this it should become much faster.

I haven't benchmarked this change in x86, but can do so if you like -- under wasm I don't think we can produce a meaningful comparison.

mglisse · 2021-03-05T06:49:55Z

Yes, I was thinking of a bench on x86_64. It would allow a comparison with the current code, and it would give some idea of what to expect compared to no-intervals once exceptions are properly supported on wasm.

pentacular · 2021-03-05T08:00:25Z

I'm probably missing something completely obvious, but I don't see any obvious way to run a benchmark suite.
I do see per package benchmarks, but I don't see suites for those either.

Should I run some ad hoc benchmarks (if so, which packages do you think would be most useful)?
Or is there some benchmark documentation that I should have found? :)

mglisse · 2021-03-05T08:08:58Z

Indeed, there is no benchmark suite. Some of the per package benchmarks are outdated and don't even compile...
I think anything that spends a lot of time in Interval_nt would do.

lrineau · 2021-03-05T08:12:01Z

Indeed, there is no benchmark suite. Some of the per package benchmarks are outdated and don't even compile...
I think anything that spends a lot of time in Interval_nt would do.

I think in CGAL we should treat benchmarks directories as demos directories:

compiled by the daily test suite,
but not run.

That would at least force us to update the code of the benchmarks in case of breaking changes.

(Cc @CGAL/geometryfactory)

pentacular · 2021-03-05T09:40:28Z

Well, I'm having trouble finding a significant difference.

I couldn't find a reasonable benchmark, so I adapted corefinement_mesh_union to run

  for (int i = 0; i < 2000; i++) {
    Mesh working_mesh1(mesh1);
    Mesh working_mesh2(mesh2);
    Mesh out;
    valid_union = PMP::corefine_and_compute_union(working_mesh1, working_mesh2, out);
  }

time ./corefinement_mesh_union
Union was successfully computed

real 0m23.875s
user 0m23.869s
sys 0m0.002s

time ./corefinement_mesh_union_nearest
Union was successfully computed

real 0m23.757s
user 0m23.735s
sys 0m0.026s

time ./corefinement_mesh_union
Union was successfully computed

real 0m23.048s
user 0m23.009s
sys 0m0.036s

time ./corefinement_mesh_union_nearest
Union was successfully computed

real 0m23.030s
user 0m23.012s
sys 0m0.022s

I then went and broke FPU.h inside a CGAL_ALWAYS_ROUND_TO_NEAREST just to make sure corefinement_mesh_union_nearest was using it.

Let me know if you have a better idea for measuring it, but it doesn't seem like there's much of a hit.

mglisse · 2021-03-05T14:01:53Z

Well, I'm having trouble finding a significant difference.

I am not sure how you tested exactly, but by default, we have CGAL_USE_SSE2 which enables a different code path that avoids CGAL_IA_ADD, etc. By the way, the combination CGAL_USE_SSE2 and CGAL_ALWAYS_ROUND_TO_NEAREST is likely broken in your branch.

Anyway, running Triangulation3/simple_2 100000 compiled with -DCGAL_NO_STATIC_FILTERS -DCGAL_NO_STRUCTURAL_FILTERING so the time is actually spent in intervals, I get .52s for normal CGAL, 1.24s if I disable CGAL_USE_SSE2, 5.1s if I define CGAL_ALWAYS_ROUND_TO_NEAREST, 7.7s if I edit Filtered_predicate to make it go straight to Mpzf, and 26s if I disable Mpzf so it uses rationals.

So it may still help to filter rationals, but the difference is significant.

pentacular · 2021-03-05T14:34:16Z

Ah, you're absolutely right.

time ./corefinement_mesh_union_nearest
Union was successfully computed

real 0m30.595s
user 0m30.512s
sys 0m0.038s

$ time ./corefinement_mesh_union

Union was successfully computed

real 0m24.609s
user 0m23.882s
sys 0m0.070s

$ time ./corefinement_mesh_union_nearest
Union was successfully computed

real 0m32.082s
user 0m31.778s
sys 0m0.013s

$ time ./corefinement_mesh_union

Union was successfully computed

real 0m23.572s
user 0m23.348s
sys 0m0.023s

So, CGAL_ALWAYS_ROUND_TO_NEAREST seems to take about 130% as long for this.

That makes much more sense.

Thanks. :)

mglisse · 2021-03-18T13:27:17Z

Is it possible that this PR is responsible for the failure of the testsuite today ?

Could be, since it touches those files, although it wasn't supposed to affect that target.

By the way, is there any non-x86_64 platform in the testsuite? At least one arm (32 or 64) or x86 (32) would be nice to check that the path without CGAL_USE_SSE2 is not broken.

lrineau · 2021-03-18T14:37:31Z

Is it possible that this PR is responsible for the failure of the testsuite today ?

Could be, since it touches those files, although it wasn't supposed to affect that target.

By the way, is there any non-x86_64 platform in the testsuite? At least one arm (32 or 64) or x86 (32) would be nice to check that the path without CGAL_USE_SSE2 is not broken.

We have a Linux Fedora running 32 bits, in the testsuite, but it seems to show up quite late in the day.

We also have a new Apple M1 machine (ARM64) running MacOS BigSur. But that one seems to be off for about a week.

I will check with @maxGimeno...

sbrian23 · 2021-03-18T15:01:46Z

I updated master from upstream and merged into a sub-branch to test.

On my machine, at least, both Interval_nt and Interval_nt_nearest pass.
So I guess it is either due to an interaction from another branch merged in for testing, or due to build options.
Or I've done something wrong in verification. :)

Would it be difficult to run the test suite with just this branch patched in?

danston · 2021-03-19T10:03:21Z

I have tested this PR in isolation on my x86_64 Apple machine and I reproduce the error from the test suite. So it is not caused by merging with other branches. Both tests bench_interval and Interval_nt fail as in the test suite. Change of the build options does not help either. The error comes from the files FPU.h and Interval_nt.h. I am not familiar with the modifications, but the error seems to come from the part that removes these lines:

#define CGAL_IA_SQRT(a) \
        CGAL_IA_FORCE_TO_DOUBLE(CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(a)))

in the file FPU.h (lines 370-380). By naively adding this macro back and undoing the corresponding modifications in the Interval_nt.h, I fixed the errors in both tests.

It is not a solution obviously but rather an idea of where to look to fix the actual error. Please let me know if more testing on my platform should be done.

mglisse · 2021-03-19T11:15:38Z

Number_types/include/CGAL/FPU.h

+  return (d > 0.0) ? nextafter(std::sqrt(d), 0.) : 0.0;
+#else
+  FPU_set_cw(CGAL_FE_DOWNWARD);
+  double i = (d > 0.0) ? CGAL_BUG_SQRT(CGAL_IA_STOP_CPROP(d)) : 0.0;


Did you lose CGAL_IA_FORCE_TO_DOUBLE somewhere?

Yes, I think so.

I went over the substitutions manually and have included them in the code as temporary notes.

I found one CGAL_IA_FORCE_TO_DOUBLE that looks like it was omitted.

sbrian23 · 2021-03-20T10:22:35Z

Thanks for isolating the issue.

I went over the code and I've put in a tentative fix, and added some temporary notes.

Please let me know how it goes on the x86_64 Apple target.

danston · 2021-03-22T10:04:04Z

Ok, now, it seems to work locally on my x86_64 machine. Both tests succeed.

maxGimeno · 2021-04-06T12:47:08Z

successfully tested in https://cgal.geometryfactory.com/CGAL/testsuite/results-5.3-Ic-98.shtml

mglisse · 2021-04-27T01:15:34Z

With -DCGAL_ALWAYS_ROUND_TO_NEAREST on master, the test Interval_nt_new.cpp fails because of _test_algebraic_structure.h assert( CGAL_NTS square( AS (23)) == AS (23*23) );, and indeed I don't see how it could work with nextafter (adding possibly would help, but there are many other strict tests around there). Is that test supposed to be skipped? Or is CGAL_ALWAYS_ROUND_TO_NEAREST considered experimental enough that it does not need to pass the testsuite?

sbrian23 · 2021-04-27T04:50:11Z

I guess the test expectations should consider CGAL_ALWAYS_ROUND_TO_NEAREST?

lrineau · 2021-04-27T12:17:30Z

With -DCGAL_ALWAYS_ROUND_TO_NEAREST on master, the test Interval_nt_new.cpp fails because of _test_algebraic_structure.h assert( CGAL_NTS square( AS (23)) == AS (23*23) );, and indeed I don't see how it could work with nextafter (adding possibly would help, but there are many other strict tests around there). Is that test supposed to be skipped? Or is CGAL_ALWAYS_ROUND_TO_NEAREST considered experimental enough that it does not need to pass the testsuite?

I think we should move forward and try to make CGAL_ALWAYS_ROUND_TO_NEAREST pass the testsuite. @maxGimeno What could be the right platform to test that macro?

maxGimeno · 2021-04-27T12:42:18Z

What about Fedora-rawhide ? It tests almost everything and it is usually all green.

pentacular added 8 commits March 2, 2021 23:28

Patch for intervals on wasm

6d048b9

wip

7f8f7e2

Add tests

d3d24fa

Clean up

ecda444

A little bit more

605a446

Revert PMP example

85c6b3b

Cleanup

4a67067

Revert out-of-scope fix

62d1f6b

lrineau added CLA/CAA The contributor must sign a CLA or a CAA so that the contribution can be merged into CGAL Enhancement Pkg::Number_types labels Mar 3, 2021

lrineau requested a review from mglisse March 3, 2021 11:44

sloriot reviewed Mar 3, 2021

View reviewed changes

mglisse reviewed Mar 3, 2021

View reviewed changes

Keep CGAL_IA_FORCE_TO_DOUBLE for the normal case

71bfcce

mglisse reviewed Mar 3, 2021

View reviewed changes

Clean up sqrt.

e545fcf

Prevent CGAL_USE_SSE2 with CGAL_ALWAYS_ROUND_TO_NEAREST

0f573cd

maxGimeno added Tests failing and removed Under Testing labels Mar 19, 2021

mglisse reviewed Mar 19, 2021

View reviewed changes

pentacular added 3 commits March 20, 2021 17:57

Merge branch 'master' into wasm

b05f606

Work through expansions to find test failure

97ae3cb

Clean up whitespace

6b75fed

maxGimeno removed the Tests failing label Mar 22, 2021

maxGimeno added the Under Testing label Mar 24, 2021

Clean up temporary notes

c26cb00

maxGimeno added the Tested label Apr 6, 2021

lrineau added rm only: ready for master For the release team only: that indicates that a PR is about to be merged in 'master' and removed Ready to be tested Under Testing labels Apr 6, 2021

lrineau merged commit d46dcec into CGAL:master Apr 6, 2021

lrineau removed the rm only: ready for master For the release team only: that indicates that a PR is about to be merged in 'master' label Apr 7, 2021

mglisse mentioned this pull request May 20, 2021

slow nextafter #5719

Open

lrineau added Merged_in_5.3-beta1 Merged_in_5.3 labels Jul 2, 2021

sloriot mentioned this pull request Oct 20, 2021

About -frounding-math in GCC #3180

Closed

Round-to-nearest support for intervals #5510

Round-to-nearest support for intervals #5510

Conversation

pentacular commented Mar 3, 2021

Summary of Changes

Release Management

sloriot commented Mar 3, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mglisse commented Mar 4, 2021

pentacular commented Mar 5, 2021

mglisse commented Mar 5, 2021

pentacular commented Mar 5, 2021

mglisse commented Mar 5, 2021

lrineau commented Mar 5, 2021

pentacular commented Mar 5, 2021

mglisse commented Mar 5, 2021 • edited Loading

pentacular commented Mar 5, 2021

mglisse commented Mar 18, 2021

lrineau commented Mar 18, 2021

sbrian23 commented Mar 18, 2021

danston commented Mar 19, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

sbrian23 commented Mar 20, 2021

danston commented Mar 22, 2021

maxGimeno commented Apr 6, 2021

mglisse commented Apr 27, 2021

sbrian23 commented Apr 27, 2021

lrineau commented Apr 27, 2021

maxGimeno commented Apr 27, 2021

mglisse commented Mar 5, 2021 •

edited

Loading

danston commented Mar 19, 2021 •

edited

Loading