Are predicates used in CDT safe?

[msokalski]

It seams to me that predicates library can be easily fooled by flipping sign and/or swapping x,y coordinates of a point to form 4 exactly concyclic set of points. Results are not very optimistic ... any thoughts? Do you reproduce same behaviour?

// example where adaptive differs from exact predicate
// adaptive seams to be wrong (I think rising xxx_errbound could fix this)

ax	-1.4048103911580639e-11	const double
ay	5.3534412196824579e-12	const double
bx	1.4048103911580639e-11	const double
by	-5.3534412196824579e-12	const double
cx	-1.4048103911580639e-11	const double
cy	-5.3534412196824579e-12	const double
dx	-5.3534412196824579e-12	const double
dy	-1.4048103911580639e-11	const double

predicates::adaptive::incircle<double> returned	1.5785423408550432e-75	double
predicates::exact::incircle<double> returned	0.0000000000000000	double
crude_xa::incircle<double> returned 		0.0000000000000000	double

// ---------------------------------------------------
// example where adaptive & exact predicates seams to be both wrong!

ax	-9.1428964336885931e-10	const double
ay	-5.2844166345548637e-19	const double
bx	-5.2844166345548637e-19	const double
by	-9.1428964336885931e-10	const double
cx	5.2844166345548637e-19	const double
cy	-9.1428964336885931e-10	const double
dx	-5.2844166345548637e-19	const double
dy	9.1428964336885931e-10	const double

predicates::adaptive::incircle<double> returned	1.2392294127517818e-119	double
predicates::exact::incircle<double> returned	1.2392294127517818e-119	double
crude_xa::incircle<double> returned 		0.0000000000000000	double

EDIT:
adding a more representative set of 4 points (both adaptive and exact predicates seams to be wrong here):

	double ax = -1.00;
	double ay = -0.01;
	double bx = -0.01;
	double by = -1.00;
	double cx =  0.01;
	double cy = -1.00;
	double dx = -0.01;
	double dy =  1.00;

I've narrowed the problem to the Dekker's product function. Forcing use of std::fma() solves the problem, probably at the cost of performance drop on platforms not having fma in hardware.

[artem-ogre]

O, horror! You are right. I managed to track the bug down, see the recent issue. Could you please check if the fix works for you?

[msokalski]

Yea fixed.
I made same thing a bit less invasive, just before #include "predicates.h" I've added:

#undef FP_FAST_FMAF
#define FP_FAST_FMAF
#undef FP_FAST_FMA
#define FP_FAST_FMA
#undef FP_FAST_FMAL
#define FP_FAST_FMAL

[artem-ogre]

Yes, I just see no reason to do Dekker product instead of FMA. Even if FMA is not done with a special instruction, it should not be much slower.

[ramcdona]

I've found a reason to turn off hardware FMA.

When hardware FMA is used on X86 binaries, those files will not run properly on Apple Silicon under Rosetta2 as Rosetta does not implement :

Rosetta translates all x86_64 instructions, but it doesn’t support the execution of some newer instruction sets and processor features, such as AVX, AVX2, and AVX512 vector instructions. If you include these newer instructions in your code, execute them only after verifying that they are available. For example, to determine if AVX512 vector instructions are available, use the sysctlbyname function to check the hw.optional.avx512f attribute.

[artem-ogre]

@ramcdona
As I understand if the code is compiled for x86_64 no special instructions (including hardware FMA) will be generated. One needs to opt-in using, e.g., -march and co.

So that’s not really a problem to do anything about, what do you think?

[msokalski]

Funny it didn't come up any sooner 🤯
Thanks for quick response, I was in a panic for a while.

[ramcdona]

I found this post

That discusses potential drawbacks to Dekker's Algorithm as a software FMA implementation.

This makes me think that I might be able to disable the FMA algorithm simply via compiler options.

If not, then I think it may be possible to dig into one of the C librarys and find a software FMA implementation there that we can pull out and use here optionally.

[artem-ogre]

@ramcdona
Also, CDT should compile with c++98 and boost. There is no way of using std::fma without making c++11 a requirement.

[ramcdona]

I didn't realize you were targeting C++98. Good on ya. I'm being forced to move from 11 to 17 because a new library I need requires it (for no good reason). Ugh...

[ramcdona]

@msokalski Looks like CDT has fixed the issue in predicates.h -- you should consider updating to their version.

[msokalski]

I'm glad to hear you've found the solution. Please confirm if replacing the file works fine on all platforms you're targeting.

[artem-ogre]

@ramcdona
Hehe, thanks, I appreciate it :) I also find the situation with different standards quite broken. On a bright side of things, at least you can upgrade. I would love to use at least c++11, but there are some forsaken Microsoft platforms we still target. It’s good for karma to take the pain as a lib author, so that the users don’t have to, but c++11 happened more than a decade ago…