avoid calling exactinit

[mlivesu]

Hi,

in case you are interested, I managed to statically define the error bounds used by the Shewchuk predicates using standard C defines, so that it is no longer necessary to compute the machine epsilon in the way exactinit() does. It's not a big deal, but allows people to directly use the exact version of the predicates without having to call exactinit at the beginning of the program. The initialization is as follows:

#include <float.h>
#define SHEWCHUK_EPSILON DBL_EPSILON/2;

const REAL resulterrbound = ( 3.0 + 8.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL ccwerrboundA = ( 3.0 + 16.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL ccwerrboundB = ( 2.0 + 12.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL ccwerrboundC = ( 9.0 + 64.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON * SHEWCHUK_EPSILON;

const REAL o3derrboundA = ( 7.0 + 56.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL o3derrboundB = ( 3.0 + 28.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL o3derrboundC = (26.0 + 288.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON * SHEWCHUK_EPSILON;

const REAL iccerrboundA = (10.0 + 96.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL iccerrboundB = ( 4.0 + 48.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL iccerrboundC = (44.0 + 576.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON * SHEWCHUK_EPSILON;

const REAL isperrboundA = (16.0 + 224.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL isperrboundB = ( 5.0 + 72.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON;
const REAL isperrboundC = (71.0 + 1408.0 * SHEWCHUK_EPSILON) * SHEWCHUK_EPSILON * SHEWCHUK_EPSILON;

const REAL splitter = (1 << (DBL_MANT_DIG/2+1)) + 1.0;

Here is the full C file containing the (modified) Shewchuk predicates

cheers,
cino

[alecjacobson]

This seems great. Do you mind preparing a PR?

[jdumas]

Nice. Does it work with both float and double?

[mlivesu]

I extensively tested the double version, never tried with floats. However, everything is based on DBL_EPSILON and DBL_MANT_DIG, which exist also with the FLT_ prefix, so it should be fine.

I wouldn't be worried about the epsilon: for the way it is defined it just can't be wrong. In case things don't work, I bet it's because of that /2+1 in

const REAL splitter = (1 << (DBL_MANT_DIG/2+1)) + 1.0;

which is there because this is the code that matches the splitter value obtained with the original exactinit(), although

const REAL splitter = (1 << (DBL_MANT_DIG/2)) + 1.0;

seems more correct to me.

@alecjacobson this float/double version of the predicates you have requires more than a trivial copy/paste from my library. I will prepare a PR asaic

[jdumas]

I think we should also add a unit test on the libigl side to make sure the values match (I trust your code, but having a unit test for it is better). In any case, triangle also has a copy of exactinit() and might benefit from the same trick.

On a side note, people might still want to use their own version of triangle.c/predicates.c, so we should still call an empty exactinit() (using Mayer's singleton trick) to be conservative. But still, that's neat!

[qnzhou]

Just to be extra cautious, I am curious if different rounding mode could give rise to different values of SHEWCHUK_EPSILON using Shewchuk's original implementation?