On Tue, Sep 26, 2017 at 01:27:52PM +0200, Jan Stary 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.) 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?)
Even on modern amd64s integer arithmetics and bitwise operations are
faster (and more precise in many cases) than floating point
equivalents. Recently I did few measurements on "real life" signal
processing code (intel i5-2500K).
I don't know if this is true for above condition and even if it was,
I'm not sure the optimization would be visible in the execution time
of programs using sinf() a lot.
Hex constants would be "useful" to avoid rounding errors in
decimal->binary conversions done during compilation. This could be
avoided by using hex float literals, but these didn't exist 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;
interestingly, 0x7f800001 is also NaN but doesn't pass the test.
AFAICS, to implement what the comment says, the test would be:
else if (ix >= 0x7f800000) 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("")?
>
Above function could also just set x to 0x7fffffff.
For sure nan(3) and isnan(3) were not common in 1993, so this might be
for portability. Furthermore, we used to support non-IEEE-754 machines
which may need this, I don't known.
