summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorbde <bde@FreeBSD.org>2005-11-24 02:04:26 +0000
committerbde <bde@FreeBSD.org>2005-11-24 02:04:26 +0000
commitcaae9bf0814e6d393e79b98b8b66b456cf281a12 (patch)
tree65ffbfd0dcba12473103936c33020f757ea4c7f8
parentf9340d98845e83aca6b2febdfa5a1ca70a412a2e (diff)
downloadFreeBSD-src-caae9bf0814e6d393e79b98b8b66b456cf281a12.zip
FreeBSD-src-caae9bf0814e6d393e79b98b8b66b456cf281a12.tar.gz
Optimized by eliminating the special case for 0.67434 <= |x| < pi/4.
A single polynomial approximation for tan(x) works in infinite precision up to |x| < pi/2, but in finite precision, to restrict the accumulated roundoff error to < 1 ulp, |x| must be restricted to less than about sqrt(0.5/((1.5+1.5)/3)) ~= 0.707. We restricted it a bit more to give a safety margin including some slop for optimizations. Now that we use double precision for the calculations, the accumulated roundoff error is in double-precision ulps so it can easily be made almost 2**29 times smaller than a single-precision ulp. Near x = pi/4 its maximum is about 0.5+(1.5+1.5)*x**2/3 ~= 1.117 double-precision ulps. The minimax polynomial needs to be different to work for the larger interval. I didn't increase its degree the old degree is just large enough to keep the final error less than 1 ulp and increasing the degree would be a pessimization. The maximum error is now ~0.80 ulps instead of ~0.53 ulps. The speedup from this optimization for uniformly distributed args in [-2pi, 2pi] is 28-43% on athlons, depending on how badly gcc selected and scheduled the instructions in the old version. The old version has some int-to-float conversions that are apparently difficult to schedule well, but gcc-3.3 somehow did everything ~10 cycles or ~10% faster than gcc-3.4, with the difference especially large on AXPs. On A64s, the problem seems to be related to documented penalties for moving single precision data to undead xmm registers. With this version, the speed is cycles is almost independent of the athlon and gcc version despite the large differences in instruction selection to use the FPU on AXPs and SSE on A64s.
-rw-r--r--lib/msun/src/k_tanf.c23
1 files changed, 7 insertions, 16 deletions
diff --git a/lib/msun/src/k_tanf.c b/lib/msun/src/k_tanf.c
index 77bc7b71..263ec96 100644
--- a/lib/msun/src/k_tanf.c
+++ b/lib/msun/src/k_tanf.c
@@ -21,16 +21,15 @@ static char rcsid[] = "$FreeBSD$";
#include "math.h"
#include "math_private.h"
+/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
static const double
-pio4 = M_PI_4,
-/* |tan(x)/x - t(x)| < 2**-29.1 (~[-1.72e-09, 1.719e-09]). */
T[] = {
- 0x1555545f8b54d0.0p-54, /* 0.333333104424423432022 */
- 0x111160cdc2c9af.0p-55, /* 0.133342838734802765499 */
- 0x1b9097e5693cd0.0p-57, /* 0.0538375346701457369036 */
- 0x173b2333895b6f.0p-58, /* 0.0226865291791357691353 */
- 0x19fcb197e825ab.0p-60, /* 0.00634450313965243938713 */
- 0x1d5f3701b44a27.0p-60, /* 0.00717088210082520490646 */
+ 0x15554d3418c99f.0p-54, /* 0.333331395030791399758 */
+ 0x1112fd38999f72.0p-55, /* 0.133392002712976742718 */
+ 0x1b54c91d865afe.0p-57, /* 0.0533812378445670393523 */
+ 0x191df3908c33ce.0p-58, /* 0.0245283181166547278873 */
+ 0x185dadfcecf44e.0p-61, /* 0.00297435743359967304927 */
+ 0x1362b9bf971bcd.0p-59, /* 0.00946564784943673166728 */
};
#ifdef INLINE_KERNEL_TANF
@@ -44,10 +43,6 @@ __kernel_tandf(double x, int iy)
GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff;
- if(ix>=0x3f2ca140) { /* |x|>=0.67434 */
- if(hx<0) {x = -x;}
- x = pio4-x;
- }
z = x*x;
w = z*z;
/* Break x^5*(T[1]+x^2*T[2]+...) into
@@ -60,10 +55,6 @@ __kernel_tandf(double x, int iy)
r = z*s*(r+v);
r += T[0]*s;
w = x+r;
- if(ix>=0x3f2ca140) {
- v = (double)iy;
- return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
- }
if(iy==1) return w;
else return -1.0/w;
}
OpenPOWER on IntegriCloud