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, 194 insertions, 25 deletions

diff --git a/newlib/libm/common/s_remquo.c b/newlib/libm/common/s_remquo.c
index bd52f39fc..9e18ab132 100644
--- a/newlib/libm/common/s_remquo.c
+++ b/newlib/libm/common/s_remquo.c
@@ -1,39 +1,208 @@
-/* Copyright (C) 2002 by  Red Hat, Incorporated. All rights reserved.
+/* Adapted for Newlib, 2009.  (Allow for int < 32 bits; return *quo=0 during
+ * errors to make test scripts easier.)  */
+/* @(#)e_fmod.c 1.3 95/01/18 */
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  *
- * Permission to use, copy, modify, and distribute this software
- * is freely granted, provided that this notice is preserved.
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice 
+ * is preserved.
+ * ====================================================
  */
+/*
+FUNCTION
+<<remquo>>, <<remquof>>--remainder and part of quotient
+INDEX
+	remquo
+INDEX
+	remquof
 
-#include "fdlibm.h"
+ANSI_SYNOPSIS
+	#include <math.h>
+	double remquo(double <[x]>, double <[y]>, int *<[quo]>);
+	float remquof(float <[x]>, float <[y]>, int *<[quo]>);
+
+DESCRIPTION
+The <<remquo>> functions compute the same remainder as the <<remainder>>
+functions; this value is in the range -<[y]>/2 ... +<[y]>/2.  In the object
+pointed to by <<quo>> they store a value whose sign is the sign of <<x>>/<<y>>
+and whose magnitude is congruent modulo 2**n to the magnitude of the integral
+quotient of <<x>>/<<y>>.  (That is, <<quo>> is given the n lsbs of the
+quotient, not counting the sign.)  This implementation uses n=31 if int is 32
+bits or more, otherwise, n is 1 less than the width of int.
+
+For example:
+.	remquo(-29.0, 3.0, &<[quo]>)
+returns -1.0 and sets <[quo]>=10, and
+.	remquo(-98307.0, 3.0, &<[quo]>)
+returns -0.0 and sets <[quo]>=-32769, although for 16-bit int, <[quo]>=-1.  In
+the latter case, the actual quotient of -(32769=0x8001) is reduced to -1
+because of the 15-bit limitation for the quotient.
+
+RETURNS
+When either argument is NaN, NaN is returned.  If <[y]> is 0 or <[x]> is
+infinite (and neither is NaN), a domain error occurs (i.e. the "invalid"
+floating point exception is raised or errno is set to EDOM), and NaN is
+returned.
+Otherwise, the <<remquo>> functions return <[x]> REM <[y]>.
 
-#ifndef _DOUBLE_IS_32BITS
+BUGS
+IEEE754-2008 calls for <<remquo>>(subnormal, inf) to cause the "underflow"
+floating-point exception.  This implementation does not.
 
-#ifdef __STDC__
-	double remquo(double x, double y, int *quo)	/* wrapper remquo */
-#else
-	double remquo(x,y,quo)			/* wrapper remquo */
-	double x,y;
-        int *quo;
+PORTABILITY
+C99, POSIX.
+
+*/
+
+#include <limits.h>
+#include <math.h>
+#include "fdlibm.h"
+
+/* For quotient, return either all 31 bits that can from calculation (using
+ * int32_t), or as many as can fit into an int that is smaller than 32 bits.  */
+#if INT_MAX > 0x7FFFFFFFL
+  #define QUO_MASK 0x7FFFFFFF
+# else
+  #define QUO_MASK INT_MAX
 #endif
+
+static const double Zero[] = {0.0, -0.0,};
+
+/*
+ * Return the IEEE remainder and set *quo to the last n bits of the
+ * quotient, rounded to the nearest integer.  We choose n=31--if that many fit--
+ * because we wind up computing all the integer bits of the quotient anyway as
+ * a side-effect of computing the remainder by the shift and subtract
+ * method.  In practice, this is far more bits than are needed to use
+ * remquo in reduction algorithms.
+ */
+double
+remquo(double x, double y, int *quo)
 {
-        int signx, signy, signres;
-        int mswx;
-        int mswy;
-        double x_over_y;
+	__int32_t n,hx,hy,hz,ix,iy,sx,i;
+	__uint32_t lx,ly,lz,q,sxy;
 
-        GET_HIGH_WORD(mswx, x);
-        GET_HIGH_WORD(mswy, y);
+	EXTRACT_WORDS(hx,lx,x);
+	EXTRACT_WORDS(hy,ly,y);
+	sxy = (hx ^ hy) & 0x80000000;
+	sx = hx&0x80000000;		/* sign of x */
+	hx ^=sx;		/* |x| */
+	hy &= 0x7fffffff;	/* |y| */
 
-        signx = (mswx & 0x80000000) >> 31;
-        signy = (mswy & 0x80000000) >> 31;
+    /* purge off exception values */
+	if((hy|ly)==0||(hx>=0x7ff00000)||	/* y=0,or x not finite */
+	  ((hy|((ly|-ly)>>31))>0x7ff00000))  {	/* or y is NaN */
+	    *quo = 0;	/* Not necessary, but return consistent value */
+	    return (x*y)/(x*y);
+	}
+	if(hx<=hy) {
+	    if((hx<hy)||(lx<ly)) {
+		q = 0;
+		goto fixup;	/* |x|<|y| return x or x-y */
+	    }
+	    if(lx==ly) {
+		*quo = 1;
+		return Zero[(__uint32_t)sx>>31];	/* |x|=|y| return x*0*/
+	    }
+	}
 
-        signres = (signx ^ signy) ? -1 : 1;
+    /* determine ix = ilogb(x) */
+	if(hx<0x00100000) {	/* subnormal x */
+	    if(hx==0) {
+		for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
+	    } else {
+		for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
+	    }
+	} else ix = (hx>>20)-1023;
 
-        x_over_y = fabs(x / y);
+    /* determine iy = ilogb(y) */
+	if(hy<0x00100000) {	/* subnormal y */
+	    if(hy==0) {
+		for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
+	    } else {
+		for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
+	    }
+	} else iy = (hy>>20)-1023;
 
-        *quo = signres * (lrint(x_over_y) & 0x7f);
+    /* set up {hx,lx}, {hy,ly} and align y to x */
+	if(ix >= -1022) 
+	    hx = 0x00100000|(0x000fffff&hx);
+	else {		/* subnormal x, shift x to normal */
+	    n = -1022-ix;
+	    if(n<=31) {
+	        hx = (hx<<n)|(lx>>(32-n));
+	        lx <<= n;
+	    } else {
+		hx = lx<<(n-32);
+		lx = 0;
+	    }
+	}
+	if(iy >= -1022) 
+	    hy = 0x00100000|(0x000fffff&hy);
+	else {		/* subnormal y, shift y to normal */
+	    n = -1022-iy;
+	    if(n<=31) {
+	        hy = (hy<<n)|(ly>>(32-n));
+	        ly <<= n;
+	    } else {
+		hy = ly<<(n-32);
+		ly = 0;
+	    }
+	}
 
-        return remainder(x,y);
-}
+    /* fix point fmod */
+	n = ix - iy;
+	q = 0;
+	while(n--) {
+	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	    if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
+	    else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;}
+	    q <<= 1;
+	}
+	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
+	if(hz>=0) {hx=hz;lx=lz;q++;}
 
-#endif /* defined(_DOUBLE_IS_32BITS) */
+    /* convert back to floating value and restore the sign */
+	if((hx|lx)==0) {			/* return sign(x)*0 */
+	    q &= QUO_MASK;
+	    *quo = (sxy ? -q : q);
+	    return Zero[(__uint32_t)sx>>31];
+	}
+	while(hx<0x00100000) {		/* normalize x */
+	    hx = hx+hx+(lx>>31); lx = lx+lx;
+	    iy -= 1;
+	}
+	if(iy>= -1022) {	/* normalize output */
+	    hx = ((hx-0x00100000)|((iy+1023)<<20));
+	} else {		/* subnormal output */
+	    n = -1022 - iy;
+	    if(n<=20) {
+		lx = (lx>>n)|((__uint32_t)hx<<(32-n));
+		hx >>= n;
+	    } else if (n<=31) {
+		lx = (hx<<(32-n))|(lx>>n); hx = sx;
+	    } else {
+		lx = hx>>(n-32); hx = sx;
+	    }
+	}
+fixup:
+	INSERT_WORDS(x,hx,lx);
+	y = fabs(y);
+	if (y < 0x1p-1021) {
+	    if (x+x>y || (x+x==y && (q & 1))) {
+		q++;
+		x-=y;
+	    }
+	} else if (x>0.5*y || (x==0.5*y && (q & 1))) {
+	    q++;
+	    x-=y;
+	}
+	GET_HIGH_WORD(hx,x);
+	SET_HIGH_WORD(x,hx^sx);
+	q &= QUO_MASK;
+	*quo = (sxy ? -q : q);
+	return x;
+}