On Sep 26 13:27:52, [email protected] wrote: > I picked sin() as an example while trying to walk though the implementation > of (pieces of) libm. If someone has the time for it, I have some questions. > > I understand the implementation originaly stems from Sun's libm of 1993. > (As does that of FreeBSD and NetBSD.)
And Android https://android.googlesource.com/platform/external/fdlibm/ and Julia http://openlibm.org/, apparently. > It has been tweaked over the years, > but the code that actually computes the values is mostly untouched. > > s_sin.c normalizes the argument to [-pi/4, +pi/4]. > This is how |x| <= pi/4 is tested: > > GET_HIGH_WORD(ix,x); > ix &= 0x7fffffff; > if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0); > > Why is it done like that? Is it faster or more portable > or in any way better that comparing the value itself to M_PI_4? > (Or was it, in 1993?) Hm, is 0x3fe921fb the nearest float <= pi/4? s_sinf.c uses 0x3f490fd8 for that test. Jan > For Inf and NaN arguments, NaN is returned as follows: > > /* sin(Inf or NaN) is NaN */ > else if (ix>=0x7ff00000) return x-x; > > Again, why the integer conversion? Would > > else if (isinf(x) || isnan(x)) > > be slower? Less portable? > > Also, is "return x-x" just a fast/clever way to return NaN > for x being Inf _or_ NaN, preferable to return nan("")? > > I have other questions in k_sin.c and elsewhere, but I'll stop here. > Please see my naive diff below. I tried this on amd64, macppc and armv7; > regress/lib/libm/fpaccuracy results in zero difference. The only > other thing I tried is "sox -n out.raw synth 10 sin" which also gives > the same result. (That's hardly "testing" of course.) > > Surely I am missing something elementary. > I'll be glad for any insight. > > Thank you for your time > > Jan > > > Index: s_sin.c > =================================================================== > RCS file: /cvs/src/lib/libm/src/s_sin.c,v > retrieving revision 1.11 > diff -u -p -r1.11 s_sin.c > --- s_sin.c 12 Sep 2016 19:47:02 -0000 1.11 > +++ s_sin.c 26 Sep 2017 08:54:03 -0000 > @@ -50,17 +50,12 @@ double > sin(double x) > { > double y[2],z=0.0; > - int32_t n, ix; > + int32_t n; > > - /* High word of x. */ > - GET_HIGH_WORD(ix,x); > - > - /* |x| ~< pi/4 */ > - ix &= 0x7fffffff; > - if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0); > - > - /* sin(Inf or NaN) is NaN */ > - else if (ix>=0x7ff00000) return x-x; > + if (-M_PI_4 < x && x < M_PI_4) > + return __kernel_sin(x,z,0); > + else if (isinf(x) || isnan(x)) > + return nan(""); > > /* argument reduction needed */ > else { >
