Re: libm-test.inc: Computing ulps near FP_ZERO.

[Brooks Moses <brooks at codesourcery dot com>]

Carlos O'Donell wrote, at 4/7/2013 7:33 AM:
> This is a very good point, one which I wasn't thinking about as
> I worked through the problem of computing ulps in glibc.
> 
> In this particular case I would argue the test suite is wrong and
> instead of expecting zero it should expect ulp(pi/2)/2, given
> that one of the expected input arguments is impossible to represent
> exactly. That makes sense.

Careful: It should expect a value in the range 0 +/- ulp(pi/2)/2, not
the particular value ulp(pi/2)/2.

> It still doesn't answer the question that I'm trying to raise here,
> which is how do we compute ulps near FP_ZERO?
> 
> Take another example, particularly the one from the LSB oliver
> testsuite which does w^z, where w = -1 -0*i and z = -1 -0*i.
> In this case both inputs are perfectly representable as floating
> point arguments but the answer is significantly off from the
> expected zero. In this case I can make no easy assumptions that
> the input arguments dictate that I adjust the answer away from
> -1 + 0*i in any appreciable way.

My claim that "if using ulp(0) = 0 gives you trouble you're doing
something wrong" applies here too.  :)

There are a couple of ways of looking at this computation.

The first is to say that any number equal to w+e, where e is a value
with magnitude less than ulp(w), will be represented as w [1] -- and so
there is this possible range of "true" inputs that could be represented
by our numerical input.  We could thus argue that the answer must be
correct for _some_ true input in that range (but not necessarily the one
exactly corresponding to the numerical input).  Thus, you get a
deviation in the output that's bounded by ulp(w).  This is effectively
what glibc's implementation is doing by converting w to polar
coordinates as the first step of the computation, and why the current
error and error bound are as they are.

The second way to look at this is to say that the implementation method
should recognize that the imaginary part of z, and one of the components
of w, are both zero -- and should as a result of that recognition make
sure that the relevant component of the output is also zero.  Here, the
error in that component is not going to be some tiny finite number; it's
going to be _exactly_ zero.

There is _no_ plausible implementation of this computation that will
give you an error in your "zero" component that is some tiny finite
number such as you suggest for ulp(0).  Either the result will be
exactly zero, or it will be something in the range of ulp(w).

[...]
> Were I to have used this solution to compute w^z one can immediately
> see that we face the same problem as cos(pi/2), in that pi is not
> representable and thus the result of sin(-0-pi) is ulp(pi)/2
> away from zero or ~2.2e-16 (0x1.0p-52).
> 
> I can't know this apriori though. The algorithm might change,
> and several other factors come into play here.

Right.  And the _real_ question is, is this algorithm good enough?  Is
it good enough to be computing (w + e1*ulp(w))^(z + e2*ulp(z)) for some
complex e1 and e2 where |e1| < 1/2 and |e2| < 1/2?

I posit no -- we want an algorithm of type 2.

There are a couple of ways to get one.  One of them is the brute-force
way; check for zeros and branch appropriately.  Glibc does this sort of
branching for other complex operations where the cost of branching is
notable, and here it's going to be minimal compared to the cost of
computing the sine and cosine.

The other way is this, which is probably faster than computing arg(t):

  r_w = |w|
  unit_w_real = w_real / r_w
  unit_w_imag = w_imag / r_w

  r_out = r_w ^ z_real
  unit_exp_zimag_real = cos(z_imag)
  unit_exp_zimag_imag = sin(z_imag)

  unit_out_real = unit_w_real * cos(z_imag) + unit_w_imag * sin(z_imag)
  unit_out_imag = -unit_w_real * sin(z_imag) + unit_w_imag * cos(z_imag)

  out_real = r_out * unit_out_real
  out_imag = r_out * unit_out_imag

Note that, if z_imag is zero, the computed cosine and sine will be
_exactly_ 0 and 1 respectively.  Further, if a component of w_real is
zero, the corresponding component of unit_w will be _exactly_ 0, and so
the resulting component of out will be _exactly_ 0.

[...]
> In the case of w^z, I see no way to make that argument. I must
> request that the test case expect -1 + 0*i, and in that case I must
> be able to compare ~0*i to +0*i and give some ulp value for this.
> 
> That is to say: How many ulps away from +0*i is ~0*i?

Well, no, you have to give some _allowable error_ value for this; ulps
are only useful as a way of getting to the allowable error value.

Here, I answer your question by saying that the answer should be exactly
zero, and ulps are thus irrelevant.  Or, if you want to go with the
customary "1 ulp" expression of the error bound, define ulp(0) as 0.

> We are back to this:
> 
> java says ulp(0x0.0p0) is 0x0.0p-1022, which is the next step to
> a normal fp value.
> 
> I expected ulp(0x0.0p0) to be 0x0.0p-1074, which is the next step
> to the smallest subnormal fp value.

...and neither of those represent a likely error magnitude for the w^z
computation.

- Brooks