@MaximSmolskiy, thank you for the patch! I think, generally this is the right direction to modify the function, and in this particular example it helps, but can it be improved even further in the same patch?

Generally, when you have (K*a*b - K*a*c)/(N*K), it makes sense to reduce it to (b-c)*a/N, because the largest source of rounding errors is when you add or subtract two big numbers or when you divide one big number by another big one, and so the methodology to reduce rounding errors is to reduce (to a reasonable level) magnitude of values that you add/subtract.

In particular, the last part of d formula is:

double d = ... + (18 * a1 * a2 * a3 - 27 * a3 * a3) * (1./108);

which could be transformed to

double d = ... + (2 * a1 * a2 - 3 * a3) * (a3/12);

here you greatly reduce magnitude of the two subtracted items and thus probably achieve higher accuracy. Probably the first part of d formula can be simplified as well.

Numeric methods are always very tricky, and probably here we need to spend a little more time in order to improve accuracy in most cases, not only in some cases.

@vpisarev I grouped terms with same total degree

Looks good to me in general. It'll be great to add a comment with the original formulae and the reason why the parts are grouped.

Added brief comment