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?)
The code is written in a way to use as little as fp ops as possible.
>
> 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?
These functions were not generally available when the code was
written. iirc these were only standarized in C99.
>
> Also, is "return x-x" just a fast/clever way to return NaN
> for x being Inf _or_ NaN, preferable to return nan("")?
Same as above.
-Otto
>
> 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 {
