2009-03-25 Craig Howland <howland@LGSInnovations.com>

	* libc/include/math.h:  (llround, llroundf): Declare.
	* libm/common/s_llround.c: New file, implementing llround().
	* libm/common/sf_llround.c: New file, implementing llroundf().
	* libm/common/sf_lround.c: Remove spurious cast in _DOUBLE_IS_32BITS
	version of function.
	* libm/common/sf_lrint.c: Ditto.
	* libm/common/sf_logb.c:  Corrected return for subnormal argument
	by replacing existing function with a version created from sf_ilogb.c.
	* libm/common/s_logb.c: Ditto, except starting point s_ilogb.c.  Also
	added documentation for logb() and logbf().
	* libm/common/s_signbit.c:  Add signbit() documentation.
	* libm/common/s_log2.c: Update return values to match what w_log2.c has,
	since log2 uses log(); add note about being derived instead of direct.
	* libm/common/sf_fma.c: Add casts to attempt to get correct results,
	as well as comments pointing out problems with the implementation.
	* libm/common/s_fma.c: Add fma() and fmaf() documentation.
	* libm/common/sf_remquo.c: Incorrect quotient returns for large values
	corrected by discarding existing function and replacing with Sun
	verion, with some enhancements.
	* libm/common/s_remquo.c: Ditto.  Add remquo() and remquof()
	documentation.
	* libm/common/s_fmax.c: Add fmax() and fmaxf() documentation.
	* libm/common/s_fmin.c: Add fmin() and fminf() documentation.
	* libm/common/s_fdim.c: Return NAN for NAN arg, add fdim() and fdimf()
	documentation.
	* libm/common/sf_fdim.c: Return NAN for NAN arg, HUGE_VALF for inf arg.
	* libm/common/s_trunc.c: Add trunc() and truncf() documentation.
	* libm/common/s_rint.c: Add rint() and rintf() documentation.
	* libm/common/s_round.c: Add round() and roundf() documentation.
	* libm/common/s_scalbn.c: Add scalbln() and scalblnf() documentation.
	* libm/common/s_infinity.c: Add infinity() and infinityf()
	documentation.
	* libm/common/s_lround.c: Add lround(), lroundf(), llround(), and
	llroundf() documentation.
	* libm/common/s_lrint.c: Add lrint(), lrintf(), llrint(), and llrintf()
	documentation.
	* libm/common/isgreater.c: New file for documenting math.h function-like
	macros isgreater(), isgreaterequal(), isless(), islessequal(),
	islessgreater(), and isunordered().
	* libm/common/s_isnan.c: Add documentation for function-like macros
	fpclassify(), isfinite(), isinf(), isnan(), and isnormal().
	* libm/common/s_nearbyint.c: Add nearbyint() and nearbyintf()
	documentation.
	* libm/common/Makefile.am: Add s_llround.c (src); sf_llround.c (fsrc);
	s_fdim.def, s_fma.def, s_fmax.def, s_fmin.def,
	s_logb.def, s_lrint.def, s_lround.def, s_nearbyint.def, s_remquo.def,
        s_rint.def, s_round.def, s_signbit.def, s_trunc.def, and
        isgreater.def (chobj);
	re-name all existing chew files (chobj) to match source file base
	names (put in underscores), delete all special targets for chew files
	(leaving all to be generated by rule).
	* libm/common/Makefile.in: regenerate.
	* libm/math/w_exp2.c: Add "base 2" to documentation description (and
	delete TRAD_SYNOPSIS).
	* libm/math/w_gamma.c: Add tgamma() and tgammaf() documentation, along
	with some history behind the function names.
	* libm/math/math.tex: Add includes for newly-added documentation (see
	.def additions to common/Makefile.am and math/Makefile.am in this
	ChangeLog list), adjusted existing .def file names to match source file
	base names (added underscores); add mention of HUGE_VALF; rename
	"Version of library" section to "Error Handling" and add some text
	about floating-point exception; added section "Standards Compliance And
	Portability".
	* libm/math/Makefile.am: Add w_exp2.def (chobj);
	re-name all existing chew files (chobj) to match source file base
	names, delete all special targets for chew files (leaving all to be
	generated by rule).
	* libm/math/Makefile.in: regenerated
	* doc/makedoc.c: Change silent ignoring of commands < 5 characters
	to a failure when reading macro file for commands < 4 characters;
	add -v (verbose) option for printing some debugging information;
	get rid of spurious translation of "@*" to "*" (no source files used @*,
	so no existing doc pages were affected); clean up some compiler
	warnings.
	* doc/doc.str: add BUGS and SEEALSO sections (to match texi2pod.pl
	which has them); Remove ITEM command (redundant with makedoc built-in
	"o", not used in any present source file so nothing is lost, anyway).
	* HOWTO: New file to hold information for maintainers regarding how
	to do things.  Initial sections on documentation and ELIX levels.

1 files changed, 86 insertions, 18 deletions

diff --git a/newlib/libm/common/s_logb.c b/newlib/libm/common/s_logb.c
index dbddd1994..1be4cddb3 100644
--- a/newlib/libm/common/s_logb.c
+++ b/newlib/libm/common/s_logb.c
@@ -1,5 +1,5 @@
-
-/* @(#)s_logb.c 5.1 93/09/24 */
+/* 2009 for Newlib:  Sun's s_ilogb.c converted to be s_logb.c.  */
+/* @(#)s_ilogb.c 5.1 93/09/24 */
 /*
  * ====================================================
  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
@@ -10,33 +10,101 @@
  * is preserved.
  * ====================================================
  */
-
 /*
- * double logb(x)
- * IEEE 754 logb. Included to pass IEEE test suite. Not recommend.
- * Use ilogb instead.
+FUNCTION
+       <<logb>>, <<logbf>>--get exponent of floating-point number
+INDEX
+	logb
+INDEX
+	logbf
+
+ANSI_SYNOPSIS
+	#include <math.h>
+        double logb(double <[x]>);
+        float logbf(float <[x]>);
+
+DESCRIPTION
+The <<logb>> functions extract the exponent of <[x]>, as a signed integer value
+in floating-point format.  If <[x]> is subnormal it is treated as though it were
+normalized; thus, for positive finite <[x]>,
+@ifnottex
+1 <= (<[x]> * FLT_RADIX to the power (-logb(<[x]>))) < FLT_RADIX.
+@end ifnottex
+@tex
+$1 \leq ( x \cdot FLT\_RADIX ^ {-logb(x)} ) < FLT\_RADIX$.
+@end tex
+A domain error may occur if the argument is zero.
+In this floating-point implementation, FLT_RADIX is 2.  Which also means
+that for finite <[x]>, <<logb>>(<[x]>) = <<floor>>(<<log2>>(<<fabs>>(<[x]>))).
+
+All nonzero, normal numbers can be described as
+@ifnottex
+<[m]> * 2**<[p]>, where 1.0 <= <[m]> < 2.0.
+@end ifnottex
+@tex
+$m \cdot 2^p$, where $1.0 \leq m < 2.0$.
+@end tex
+The <<logb>> functions examine the argument <[x]>, and return <[p]>.
+The <<frexp>> functions are similar to the <<logb>> functions, but
+returning <[m]> adjusted to the interval [.5, 1) or 0, and <[p]>+1.
+
+RETURNS
+@comment Formatting note:  "$@" forces a new line
+When <[x]> is:@*
++inf or -inf, +inf is returned;@*
+NaN, NaN is returned;@*
+0, -inf is returned, and the divide-by-zero exception is raised;@*
+otherwise, the <<logb>> functions return the signed exponent of <[x]>.
+
+PORTABILITY
+ANSI C, POSIX
+
+SEEALSO
+frexp, ilogb
+*/
+
+/* double logb(double x)
+ * return the binary exponent of non-zero x
+ * logb(0) = -inf, raise divide-by-zero floating point exception
+ * logb(+inf|-inf) = +inf (no signal is raised)
+ * logb(NaN) = NaN (no signal is raised)
+ * Per C99 recommendation, a NaN argument is returned unchanged.
  */
 
 #include "fdlibm.h"
 
 #ifndef _DOUBLE_IS_32BITS
 
+double
 #ifdef __STDC__
-	double logb(double x)
+logb(double x)
 #else
-	double logb(x)
-	double x;
+logb(x)
+double x;
 #endif
 {
-	__int32_t lx,ix;
-	EXTRACT_WORDS(ix,lx,x);
-	ix &= 0x7fffffff;			/* high |x| */
-	if((ix|lx)==0) return -1.0/fabs(x);
-	if(ix>=0x7ff00000) return x*x;
-	if((ix>>=20)==0) 			/* IEEE 754 logb */
-		return -1022.0; 
-	else
-		return (double) (ix-1023); 
+	__int32_t hx,lx,ix;
+
+	EXTRACT_WORDS(hx,lx,x);
+	hx &= 0x7fffffff;		/* high |x| */
+	if(hx<0x00100000) {		/* 0 or subnormal */
+	    if((hx|lx)==0)  {
+		double  xx;
+		/* arg==0:  return -inf and raise divide-by-zero exception */
+		INSERT_WORDS(xx,hx,lx);	/* +0.0 */
+		return -1./xx;	/* logb(0) = -inf */
+		}
+	    else			/* subnormal x */
+		if(hx==0) {
+		    for (ix = -1043; lx>0; lx<<=1) ix -=1;
+		} else {
+		    for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
+		}
+	    return (double) ix;
+	}
+	else if (hx<0x7ff00000) return (hx>>20)-1023;	/* normal # */
+	else if (hx>0x7ff00000 || lx)  return x;	/* x==NaN */
+	else  return HUGE_VAL;	/* x==inf (+ or -) */
 }
 
 #endif /* _DOUBLE_IS_32BITS */