This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: [PATCH] Remove slow paths from log


On Fri, 2 Feb 2018, Wilco Dijkstra wrote:

> Interestingly removing the slow paths makes hardly any difference in practice:
> the worst case error is still ~0.5ULP, and exp(log(x)) shows identical results
> before/after on many millions of random cases.  All GLIBC math tests pass on
> AArch64 and x64 with no change in ULP error.

Could you please give the error analysis that justifies that the worst 
case error is still ~0.5ulp?

I don't expect a direct error analysis of the fast case code, though you 
can do one if you want.  Rather, in such cases, an indirect justification 
of the following form (which supposes that the existing code does have an 
error bound of 0.5ulp) is probably a lot easier to produce:

The fast case produces a result of the form HIGH + LOW, with error bound 
E, and checks that HIGH + (LOW + E) and HIGH + (LOW - E) round to the same 
value to nearest - if they don't, it uses extra precision.  That is, on 
the assumption that the existing code is correct, E is an error bound for 
HIGH + LOW.  And assuming that E is less than 0.5ulp of the correct 
result, that means we don't actually need the error calculation, and can 
get sufficiently good results just from the fast path.

I suspect that analyses of this form will allow eliminating most if not 
all of the slow paths.  But it's necessary in each case where removal of a 
slow path is produced to say what E is, what the minimum magnitude of the 
correct result is, and so how small E must be in terms of the correct 
result, to justify that the slow path can be removed.

> -  /* End stage I, case abs(x-1) < 0.03 */
> -  if ((y = b + (c + b * E2)) == b + (c - b * E2))
> -    return y;

> +  return b + (c + b * E2);

And in such a case, the analysis will imply that the b * E2 part of 
computing the result is not needed - it's actually computing one bound on 
the possible correctly rounded result - and you should just return b + c.

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]