Unexpected results with a specific rounding mode

First of all, sorry for the bad title: /

I am trying to reproduce the results of the calculation of the eigenvalues ​​of a tridiagonal symmetric matrix. I define the “upper and lower bounds” of some values, rounding to plus and minus infinity, respectively.

Instead of changing the rounding mode every time I just use a “trick”: fl⁻ (y) = -fl⁺ (-y), where fl⁻ (y) is the value of y when using minus infinity, the rounding mode and fl⁺ (y) are y value when using round to plus infinity mode. So, I have the following code snippet in C:

fesetround(FE_UPWARD);
first =  - (-d[i] + x);
second =  ( - (( e[i-1]*e[i-1] ) / a_inf ));
a_inf = first + second;

first = d[i] - x;
second = - ( ( e[i-1]*e[i-1] ) / a_sup );
a_sup = first + second;

and it works fine, except for one example where a_inf gives me the correct result, but a_sup gives the wrong result, although both the first and second variables seem to have the same value.

However, if I like it:

  fesetround(FE_UPWARD);
  first =  - (-d[i] + x);
  second =  ( - (( e[i-1]*e[i-1] ) / a_inf ));

  fesetround(FE_DOWNWARD);
  first =  - (-d[i] + x);
  second =  ( - (( e[i-1]*e[i-1] ) / a_sup ));

I get the right results. So, if I use the fl⁻ (y) = -fl⁺ (-y) trick, I get the correct results, if I change the rounding mode and use the original expression, I get the wrong results. Any idea why?

In both cases, the variables of the first and second values ​​are as follows:

first  1.031250000000000e+07,  second -1.031250000000000e+07
first  1.031250000000000e+07,  second -1.031250000000000e+07

And the correct values ​​for a_inf and a_sup are -1.862645149230957e-09 and + 1.862645149230957e-09, respectively, but in the first case a_sup = 0, which is incorrect

What I suppose this is happening is some kind of catastrophic cancellation, but I don’t know how to solve it in this case ...

Thanks in advance!