svn commit: r278154 - head/lib/msun/src

[Steve Kargl]

On Wed, Feb 04, 2015 at 04:37:14AM +1100, Bruce Evans wrote:
> On Tue, 3 Feb 2015, Steve Kargl wrote:
> 
> > -#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.

My crystal ball isn't function too well lately.

> 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.

On all architecture that FreeBSD currently supports, there is a guarantee
to the extent that src/sys/${MACHINE}/include/_limits.h shows

#define	__INT_MAX	0x7fffffff
#define	__INT_MIN	(-0x7fffffff - 1)

> 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.

Given that the widest field width in IEEE754-2008 for
binary128 is 15, which gives a maximum unbiased exponent
of 2**15, then using +-65536 should cover our needs.  The
functions scalbn[fl] can then deal with values of n 
that may result in exponents that are too small or too 
large.  Untested diff below.

-- 
Steve

Index: s_scalbln.c
===================================================================
--- s_scalbln.c	(revision 278165)
+++ s_scalbln.c	(working copy)
@@ -27,17 +27,17 @@
 #include <sys/cdefs.h>
 __FBSDID("$FreeBSD$");
 
-#include <limits.h>
 #include <math.h>
 
+#define	LARGE	65536
+#define	SMALL	-65536
+
 double
 scalbln (double x, long n)
 {
 	int in;
 
-	in = (int)n;
-	if (in != n)
-		in = (n > 0) ? INT_MAX: INT_MIN;
+	in = n > LARGE ? LARGE : n < SMALL ? SMALL : n;
 	return (scalbn(x, in));
 }
 
@@ -46,9 +46,7 @@
 {
 	int in;
 
-	in = (int)n;
-	if (in != n)
-		in = (n > 0) ? INT_MAX: INT_MIN;
+	in = n > LARGE ? LARGE : n < SMALL ? SMALL : n;
 	return (scalbnf(x, in));
 }
 
@@ -57,8 +55,6 @@
 {
 	int in;
 
-	in = (int)n;
-	if (in != n)
-		in = (n > 0) ? INT_MAX: INT_MIN;
+	in = n > LARGE ? LARGE : n < SMALL ? SMALL : n;
 	return (scalbnl(x, in));
 }