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.) 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?)
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 {
