Date: Wed, 2 Jun 2004 17:51:23 +1000 (EST) From: Bruce Evans <bde@zeta.org.au> To: David Schultz <das@freebsd.org> Cc: freebsd-i386@freebsd.org Subject: Re: i386/67469: src/lib/msun/i387/s_tan.S gives incorrect results for large inputs Message-ID: <20040602172105.T23521@gamplex.bde.org> In-Reply-To: <200406012251.i51MpkkU024224@VARK.homeunix.com> References: <200406012251.i51MpkkU024224@VARK.homeunix.com>
next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, 1 Jun 2004, David Schultz wrote: > >Description: > src/lib/msun/i387/s_tan.S returns wildly inaccuate results when > its input has a large magnitude (>> 2*pi). For example: > > input s_tan.S k_tan.c > 1.776524190754802e+269 1.773388446261095e+16 -1.367233274980565e+01 > 1.182891728897420e+57 -1.9314539773999572e-01 1.0020569035866138e+03 > 2.303439778835110e+202 2.8465460220132694e+00 3.5686329695133922e+00 > > I suspect that this is caused in the modular reduction phase, which is > used for inputs with magnitude greater than 2^63. In these cases, the > inputs are taken mod 2*pi, but the double precision representation of > pi isn't precise enough to obtain a meaningful result for this > computation. All of the i387 trig functions have this bug. fldpi actually gives a value of pi with extended precision (64 bits), and I believe fprem uses extended precision and there is some magic to give a few more bits of precision (66 altogether?). However, fdlibm is much more precise than this. It uses a 1584-bit value for 2/pi to get 113 bits of precision for range reduction in the double precision case. (The single precision case is not so precise and is broken anyway (off by more than pi/2 for many values between 32768 and 65536 according to my notes). It is also very slow.) > >How-To-Repeat: > > >Fix: > One solution might involve doing the reduction using > __ieee754_rem_pio4() instead of the fprem1 instruction. > Unfortunately, since this uses pi/2 as the modulus, it > is necessary to apply the identity tan(x + pi/2) = -1/tan(x) > for odd-phase inputs. I tried this for the first example > input above and got an answer that was off by 2 ulps. Close, > but the MI implementation gets within 1 ulp. I'm not sure what > kind of correction (if any) should be used here. I think the complete fdlibm version of tan(), sin() and cos() should be used for large args. "large" could be classified by failure of the first fptan. Or just return the TLOSS error for very large args. Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20040602172105.T23521>