| Age | Commit message (Collapse) | Author | Lines |
|---|
|
Handle when after reduction |y| > pi/4+tiny. This happens in directed
rounding modes because the fast round to int code does not give the
nearest integer. In such cases the reduction may not be symmetric
between x and -x so e.g. cos(x)==cos(-x) may not hold (but polynomial
evaluation is not symmetric either with directed rounding so fixing
that would require more changes with bigger performance impact).
The fix only adds two predictable branches in nearest rounding mode,
simple ubenchmark does not show relevant performance regression in
nearest rounding mode.
The code could be improved: e.g reducing the medium size threshold
such that two step reduction is enough instead of three, and the
single precision case can avoid the issue by doing the round to int
differently, but this fix was kept minimal.
|
|
the idiomatic rounding of x is
n = x + toint - toint;
where toint is either 1/EPSILON (x is non-negative) or 1.5/EPSILON
(x may be negative and nearest rounding mode is assumed) and EPSILON is
according to the evaluation precision (the type of toint is not very
important, because single precision float can represent the 1/EPSILON of
ieee binary128).
in case of FLT_EVAL_METHOD!=0 this avoids a useless store to double or
float precision, and the long double code became cleaner with
1/LDBL_EPSILON instead of ifdefs for toint.
__rem_pio2f and __rem_pio2 functions slightly changed semantics:
on i386 a double-rounding is avoided so close to half-way cases may
get evaluated differently eg. as sin(pi/4-eps) instead of cos(pi/4+eps)
|
|
ld128 support was added to internal kernel functions (__cosl, __sinl,
__tanl, __rem_pio2l) from freebsd (not tested, but should be a good
start for when ld128 arch arrives)
__rem_pio2l had some code cleanup, the freebsd ld128 code seems to
gather the results of a large reduction with precision loss (fixed
the bug but a todo comment was added for later investigation)
the old copyright was removed from the non-kernel wrapper functions
(cosl, sinl, sincosl, tanl) since these are trivial and the interesting
parts and comments had been already rewritten.
|
|
zero, one, two, half are replaced by const literals
The policy was to use the f suffix for float consts (1.0f),
but don't use suffix for long double consts (these consts
can be exactly represented as double).
|
|
|
