More precise abstractions of trigonometric functions using c-stubs

[FelixKrayer]

These changes allow for a more precise abstraction of the trigonometric functions on monotonic intervals.
Similar to #1253, c-stubs are used for the calculation of the resulting interval, once an interval is identified to be monotonic.
The option to enable/disable this option can be removed, if we trust, that the implementation of the trigonometric functions used in the analyzer is the same as the one used in the program to be analyzed.

You can find a rough explanation in this PDF:
Explanation_draft.pdf

• add more detailed explanation

[github-advanced-security]

Semgrep found more than 10 potential problems in the proposed changes. Check the Files changed tab for more details.

[FelixKrayer]

I have added a first draft of an explanation that hopefully helps to understand the computations in the floatDomain.

[stilscher]

Thank you for this addition! I took a look at the conceptual description and think it is a nice idea. Since I am quite busy right now, I might have to postpone the review for another two weeks.

[stilscher]

I had a few questions regarding the rounding modes. Apart from that I think it is a nice and generic idea to get a good abstraction for trigonometric computations in case the interval of the input variable is not too imprecise. I think a few more comments with some selected information from the pdf draft could improve the clarity of the code.

[FelixKrayer]

I addressed the comments, fixed some rounding modes and added some comments to the code.
Do I have to worry about the Semgrep findings semgrep.big_int_z ? I am not sure how to understand and approach these.
Also do you have any idea how to approach the arbitrary constants for approximating pi? Do we just have to leave them like this for the lack of a better approach?

[stilscher]

Semgrep warns because there are several usages of functions from Big_int_Z. Instead of using Big_int_Z which is just a wrapper for Z functions from Z should be called directly. See #992 and #1317. You can avoid the warning by substitute e.g. Big_int_Z.float_of_big_int with Z.to_float. https://github.com/goblint/analyzer/blob/master/.semgrep/zarith.yml lists the corresponding semgrep rules and how to fix the warnings.

[stilscher]

For approximating pi, @michael-schwarz and me discussed, that one could write a C-stub for M_PI and round it for a lower and upper bound in case it is not exactly representable. I added commit that implements this idea.

[FelixKrayer]

Thanks for fixing the semgrep issues. I have removed the phase shift and changed the case checks instead

[stilscher]

And there now seems to be an issue with GobView because there is no JavaScript implementation for the stub I added for pi.

[FelixKrayer]

And there now seems to be an issue with GobView because there is no JavaScript implementation for the stub I added for pi.

Would it work to do something like
#ifndef M_PI
#define M_PI -1.
#endif
in the c-stubs file and then check if Float_t.pi <= 0. and use the arbitrary under/overapprox constants in this case?

[sim642]

There's something our stubs don't account for: in C, fesetround may fail, but we're not checking/asserting its return value. Maybe adding such assertions would reveal the mode not changing/not being available on certain CPUs.

[stilscher]

So the problem seems to be that the rounding mode only specifies how floating point numbers are rounded per instruction. cos is however usually approximated with different polynomials. The implementation can vary depending on the system. The rounding mode only influences how the single instructions within this approximation are rounded. It is not safe to assume, that rounding cos(x) with rounding mode down is necessarily smaller than cos x with rounding mode up.

#include <math.h>
#include <float.h>
#include <fenv.h>
#include <stdio.h>

int main()
{
    int old_roundingmode = fegetround();

    fesetround(FE_DOWNWARD);
    printf("%.20f\n", cos(-1.0));

    fesetround(FE_UPWARD);
    printf("%.20f\n", cos(-1.0));

    fesetround(old_roundingmode);
}

On the Intel Mac this outputs:
0.54030230586813965398
0.54030230586813798866

(fesetround returns 0 which indicates that the change of the rounding mode is successful)

[stilscher]

This not only affects the trigonometric functions but also other more complex mathematical functions such as sqrt. A proposal of @michael-schwarz and me would be to make this feature configurable and compute the cos for the upper bound with both, Up and Down as rounding mode and take the greater of both (and similarly for the lower bound). This is of course not completely safe as there is no guarantee that the cos implementation is monotonic in between, but would allow to retain precision for the float values and it probably is a proper approximation in most cases.
Another option would be to add some error to account for these rounding effects, but it seems difficult to find an error bound for different implementations of the trigonometric functions that often choose different approximations for certain input intervals.

[sim642]

I might have mentioned it somewhere at some point, but there's also crlibm. Maybe we'll have to turn to something like that.

[michael-schwarz]

@sim642 How would crlibm help us here? It correctly rounds w.r.t. to the mathematical definition and not w.r.t. some arbitrarily (mis-)behaving implementation.

Am I missing something here?

[sim642]

CRlibm, an efficient and proven correctly-rounded mathematical library

That at least sounds like it should provide that, but I haven't looked further into it.

[FelixKrayer]

If we can assume that for any operation op(<implementation>, <rounding mode>, <argument>), it holds that op(CRlibm, up, x) >= op(<any impl>, <any RM>, x) (and vice versa for down), we could use CRlibm for a safe upper (& lower) bound.
However, looking at the weirdly behaving cos(), I am not sure if that assumption holds here.

[stilscher]

A proposal of @michael-schwarz and me would be to make this feature configurable and compute the cos for the upper bound with both, Up and Down as rounding mode and take the greater of both (and similarly for the lower bound). This is of course not completely safe as there is no guarantee that the cos implementation is monotonic in between, but would allow to retain precision for the float values and it probably is a proper approximation in most cases.

At our meeting today, we discussed that we would like to implement the above idea, make this feature optional (turned off by default) and comment/warn that it can be unsound.

[FelixKrayer]

At our meeting today, we discussed that we would like to implement the above idea, make this feature optional (turned off by default) and comment/warn that it can be unsound.

I implemented the option math_fun_eval_cstub to toggle this functionality. Also all previous function applications of fun down x are now guarded by min (fun up x) (fun down x) and analogously applications of fun up x are guarded by a max. This is done for all calls to trigonometric functions as well as for sqrt

However, thinking about it, I am wondering if we should include fun nearest x into the min/max. Nearest can round up or down depending on what is closer to the actual outcome. Assuming a function calculates (a op b) / (y op x), the result of this calculated with nearest can be greater than the max of up and down. This can happen if the numerator is rounded up and the denominator is rounded down, while for up and down both are rounded the same.
An example for this would be a = y = 1.; b = x = 0.000000001; in (a - b) / (y + x);, where for up and down a result less than 1 is calculated while for nearest the result is exactly 1 (reproduce with roundtest.c).

I have not yet thought about if there is a similar issue for rounding mode ToZero

[stilscher]

The failing indentation check is unrelated.