Skip site navigation (1)Skip section navigation (2)
Date:      Wed, 4 Feb 2015 04:37:14 +1100 (EST)
From:      Bruce Evans <brde@optusnet.com.au>
To:        Steve Kargl <sgk@troutmask.apl.washington.edu>
Cc:        svn-src-head@freebsd.org, svn-src-all@freebsd.org, "Pedro F. Giffuni" <pfg@freebsd.org>, src-committers@freebsd.org, Bruce Evans <brde@optusnet.com.au>
Subject:   Re: svn commit: r278154 - head/lib/msun/src
Message-ID:  <20150204035337.Y2718@besplex.bde.org>
In-Reply-To: <20150203163329.GA12706@troutmask.apl.washington.edu>
References:  <201502031417.t13EHU9g074687@svn.freebsd.org> <20150204012614.T1972@besplex.bde.org> <20150203163329.GA12706@troutmask.apl.washington.edu>

next in thread | previous in thread | raw e-mail | index | archive | help
On Tue, 3 Feb 2015, Steve Kargl wrote:

> On Wed, Feb 04, 2015 at 02:38:40AM +1100, Bruce Evans wrote:
>> On Tue, 3 Feb 2015, Pedro F. Giffuni wrote:
>>
>>> Log:
>>>  Reduce confusion in scalbnl() family of functions.
>>>
>>>  The changes unrelated to the bug in r277948 made
>>>  the code very difficult to understand to both
>>>  coverity and regular humans so take a step back
>>>  to something that is much easier to understand
>>>  for both and follows better the original code.
>>>
>>>  CID:	1267992, 1267993, 1267994
>>>  Discussed with:	kargl
>>
>> You mean, take a step backwards to something that is harder to understand.
>
> Well, the correct fix should have been to ONLY fix the typo,
> and leave the code rewrite for a second commit.
>
> Index: s_scalbln.c
> ===================================================================
> --- s_scalbln.c (revision 276768)
> +++ s_scalbln.c (working copy)
> @@ -72,5 +72,5 @@
>                else
>                        in = INT_MIN;
>        }
> -       return (scalbnl(x, (int)n));
> +       return (scalbnl(x, in));
> }
>
> But, that's water under the bridge.

Indeed.

> You forgot to include a diff.  Here's one untested attempt at
> addressing your concerns.
>
> Index: s_scalbln.c
> ===================================================================
> --- s_scalbln.c	(revision 276768)
> +++ s_scalbln.c	(working copy)
> @@ -27,21 +27,23 @@
> #include <sys/cdefs.h>
> __FBSDID("$FreeBSD$");
>
> -#include <limits.h>
> +#include <float.h>
> #include <math.h>
>
> +#define	FLT_LARGE	FLT_MAX_EXP - FLT_MIN_EXP + FLT_MANT_DIG
> +#define	FLT_SMALL	FLT_MIN_EXP - FLT_MAX_EXP
> +#define	DBL_LARGE	DBL_MAX_EXP - DBL_MIN_EXP + DBL_MANT_DIG
> +#define	DBL_SMALL	DBL_MIN_EXP - DBL_MAX_EXP
> +#define	LDBL_LARGE	LDBL_MAX_EXP - LDBL_MIN_EXP + LDBL_MANT_DIG
> +#define	LDBL_SMALL	LDBL_MIN_EXP - LDBL_MAX_EXP
> +

This is a lot of code and bugs (missing parentheses) for negative gain.
To work, it depends on all of the above expressions having results that
can be represented as an int.  There is no guarantee of this in general.
In fact, 'ln' variants of the functions only exist because this is not
possible in general -- they exist because int is not large enough to
represent the limits.  This also means that expressions like DBL_MAX_EXP
cannot have type int in general.  In fact, it is not even clear that they
have integer type.  In n869.txt, I could only find:

X       [#8] The  values  given  in  the  following  list  shall  be
X       replaced by implementation-defined constant expressions with
X       values that are greater  or  equal  in  magnitude  (absolute
X       value) to those shown, with the same sign:

This is the implementation, so it can guarantee anything it wants and
assume things that other parts of the implementation do.  Other parts
use only 32-bit ints (larger would also work) and make all the
relevant constants have type plain int.

However, we can show that the above values are within certain fixed
limits just as easily as we can show that they are within the INT32
limits, and then hard-code the fixed limits.  E,g, +-16384 might work,
and +-32768 surely works when the exponent field is limited to 8 bits.
Double a few more times for safety.  Hard-coding a limit related to
an assumed maximum supported width for the exponent field is easier
to get right that the above expressions.  I think the above are quite
buggy, and so are my +-16384 and +-32768.  The full exponent range
plus more for denormals is needed in both directions.  You only have
it for the LARGE direction.  This can be hard-coded for a 15-bit
exponent field as +-32768 +- extras for denormals, -+ a couple for
special exponents.  This needs strictly larger than 16-bit ints for
the denormals only.  Make it +-65536 to avoid having to count carefully
to avoid off-by-1 errors.  Double a few times for larger exponent fields,
or just add assertions to detect increases of the maxima, or just hard-
code the limits as the int32_t extrema.

> double
> scalbln (double x, long n)
> {
> 	int in;
>
> -	in = (int)n;
> -	if (in != n) {
> -		if (n > 0)
> -			in = INT_MAX;
> -		else
> -			in = INT_MIN;
> -	}
> +	in = n > DBL_LARGE ? DBL_LARGE : n < DBL_SMALL ? DBL_SMALL : n;
> +

New style bug: extra blank line.  I think fdlibm doesn't have many of these,
and this file is supposed to be in KNF which has none.

> 	return (scalbn(x, in));
> }

The code can't be further simlified by avoiding the local variable and the
assignment to it, since that would make the single line in the function too
long.

Bruce



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20150204035337.Y2718>