Re: [PATCH] Optimized generic expf and exp2f

[Szabolcs Nagy <szabolcs dot nagy at arm dot com>]

On 05/09/17 21:45, Joseph Myers wrote:
> On Tue, 5 Sep 2017, Szabolcs Nagy wrote:
> 
>> Error checks are inline and actual error handling is
>> in separate functions that are tail called and expected
>> to be reusable when other math functions are implemented
>> without wrappers. (expf, __expf, __ieee754_expf symbols
>> are aliases, _LIB_VERSION is not checked, errno is set
>> unconditionally according to POSIX rules.)
> 
> The _LIB_VERSION and matherr handling is part of the ABI for the existing 
> symbol versions of expf and exp2f.  Thus, if you wish to use integrated 
> error handling that does not handle _LIB_VERSION and matherr the same way 
> as the existing code, you need a new symbol version (with the old version 
> then probably keeping the existing wrapper but having it wrap the new 
> implementation - duplicate errno setting for existing binaries should be 
> fine).
> 

i tried to do this but it became complicated dealing with
targets that have their own implementation.

so my plan now is to first post the new code with the
wrapper around it (so there is redundant errno setting),
then in a separate patch try to bump the symbol version
and move to errno-only-wrapper (on all targets for expf
and exp2f) and then finally remove the wrapper for the
new generic code.

it is not yet clear to me how to do the errno-only-wrapper,
since the existing wrapper-template machinery does not
work on a per function basis, but that code would be nice
to reuse, some guidance on that would be helpful.

>> One issue with the argument reduction is that it only
>> works right with nearest rounded rint, otherwise the
>> interval is [0,2c] or [-2c,0] instead of [-c,c]. The
>> polynomial is optimized for [-c,c] but it has sufficent
>> extra precision that it gives acceptable results on
>> [-2c,2c] too assuming users are less interested in
>> non-nearest rounded precision, however this means some
>> glibc ulp error limits will need adjustment.
> 
> By ulp error limits do you mean the entries in libm-test-ulps, or the 
> global max_valid_error limit in libm-test-support.c which no 
> libm-test-ulps entry is allowed to exceed?
> 

limb-test-ulps.

expf failures on x86_64 (no-fma, lrint using +-shift):

testing float (without inline functions)
Failure: Test: exp_downward (-0x1p-20)
Result:
 is:          9.99998986e-01   0x1.ffffdep-1
 should be:   9.99999046e-01   0x1.ffffe0p-1
 difference:  5.96046447e-08   0x1.000000p-24
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: exp_downward (0x5.8b90b8p+4)
Result:
 is:          3.40279831e+38   0x1.ffff06p+127
 should be:   3.40279851e+38   0x1.ffff08p+127
 difference:  2.02824096e+31   0x1.000000p+104
 ulp       :  1.0000
 max.ulp   :  0.0000
Maximal error of `exp_downward'
 is      : 1 ulp
 accepted: 0 ulp
Failure: Test: exp_towardzero (-0x1p-20)
Result:
 is:          9.99998986e-01   0x1.ffffdep-1
 should be:   9.99999046e-01   0x1.ffffe0p-1
 difference:  5.96046447e-08   0x1.000000p-24
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: exp_towardzero (0x5.8b90b8p+4)
Result:
 is:          3.40279831e+38   0x1.ffff06p+127
 should be:   3.40279851e+38   0x1.ffff08p+127
 difference:  2.02824096e+31   0x1.000000p+104
 ulp       :  1.0000
 max.ulp   :  0.0000
Maximal error of `exp_towardzero'
 is      : 1 ulp
 accepted: 0 ulp