An attempt to fix sym_difference/get_turns failures on MinGW, QNX/QCC and PowerPC/GCC

[awulkiew]

This is an attempt to fix an issue observable on MinGW on Win, probably QCC on QNX and maybe GCC5.0 on PowerPC.

Consider a side calculation (side_by_triangle) for points ({13, 0.3}, {14, 0.4}, {20, 1}). Internally calculated value of a determinant is -2.2204460492503131e-016 on MSVC (there is no problem on the majority of tested compilers so I'm guessing that the value is similar for them). However for MinGW (I tested only this one from the "failing" ones mentioned above) the value is -2.7755575615628914e-016, but only if it's calculated as one expression (see the code). Maybe the rvalues are handled differently? Anyway if temporary values are explicitly defined and used the result is the same as the one calculated by the code compiled by MSVC (and probably the majority of compilers).

I'm not sure what's the cause, a bug, optimization or something else, I stopped digging at this stage. Maybe we should compare the assemblies and other compilers? Or maybe someone has an idea what's the cause? I also haven't checked the release build or various optimization parameters. So I don't find this workaround as a solid one. It's rather a test or a temporary workaround.

I think that such cases should be handled with more care. There is a class of issues related to the way how very small values of sides and cramer's rule are handled. If they were supported somehow more robustly by the intersection strategy then I'm guessing that this workaround wouldn't be needed and could be reverted.

[mkaravel]

Have you tried writing this as:

    return (rt(ux) * rt(vy)) - (rt(uy) * rt(vx));

It should give the same results as your alternative code.

[awulkiew]

Unfortunately, it gives the same result as the original. Btw, I'm compiling using MinGW-w64, it's also shipped with QtCreator for Windows, in case you'd like to play with it.

[barendgehrels]

Is it about overflow?
Anyway, I've no objection to this change, I'm OK with merging

[awulkiew]

I checked the dissassembly of this simple program:

double ux = 1;
double uy = 0.10000000000000003;
double vx = 7;
double vy = 0.69999999999999996;

double c1 = ux * vy;
double c2 = uy * vx;
double t0 = c1 - c2;
double t1 = ux * vy - uy * vx;

In MinGW-w64 4.9.1 for c1 and c2 the u_ and v_ variables are loaded, multiplied and stored. Then c1 and c2 are loaded subtracted with fsubl instruction and stored. However in the case of t1 after loading and multiplying the values are subtracted using fsubrp instruction.

On the other hand in MSVC 10 fsub and fsubp are used respectively.

According to the docs of FSUBR/FSUBRP/FISUBR--Reverse Subtract: These instructions perform the reverse operations of the FSUB, FSUBP, and FISUB instructions. They are provided to support more efficient coding.

Unfortunately, if O2 flag is passed into the MinGW compiler it optimizes all of the variables out and generates the result calculated by fsubr not taking into account the temporaries. MSVC with optimization enabled still generates the same result.

So this PR doesn't solve the issue. I'm going to close it and just add those test cases in get_turns test behind macro check for further investigation.