Add C99's nearbyint{,f}() functions as wrappers around rint().

These trivial implementations are about 25 times slower than
rint{,f}() on x86 due to the FP environment save/restore.
They should eventually be redone in terms of fegetround() and
bit fiddling.

4 files changed, 94 insertions, 6 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 0d10696..9dc02e3 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -86,7 +86,8 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 	s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
 	s_isnanf.c s_ldexpf.c s_lib_version.c s_log1p.c \
 	s_log1pf.c s_logb.c s_logbf.c s_matherr.c s_modff.c \
-	s_nextafter.c s_nextafterf.c s_rint.c s_rintf.c s_round.c s_roundf.c \
+	s_nearbyint.c s_nextafter.c s_nextafterf.c \
+	s_rint.c s_rintf.c s_round.c s_roundf.c \
 	s_scalbln.c s_scalbn.c s_scalbnf.c \
 	s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c s_tan.c \
 	s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c \
@@ -170,7 +171,7 @@ MLINKS+=ieee_test.3 significand.3 ieee_test.3 significandf.3
 MLINKS+=j0.3 j1.3 j0.3 jn.3 j0.3 y0.3 j0.3 y1.3 j0.3 y1f.3 j0.3 yn.3
 MLINKS+=j0.3 j0f.3 j0.3 j1f.3 j0.3 jnf.3 j0.3 y0f.3 j0.3 ynf.3
 MLINKS+=lgamma.3 gamma.3 lgamma.3 gammaf.3 lgamma.3 lgammaf.3 lgamma.3 tgamma.3
-MLINKS+=rint.3 rintf.3
+MLINKS+=rint.3 rintf.3 rint.3 nearbyint.3 rint.3 nearbyintf.3
 MLINKS+=round.3 roundf.3
 MLINKS+=sin.3 sinf.3
 MLINKS+=sinh.3 sinhf.3
diff --git a/lib/msun/man/rint.3 b/lib/msun/man/rint.3
index 3af9110..09f7e19 100644
--- a/lib/msun/man/rint.3
+++ b/lib/msun/man/rint.3
@@ -32,10 +32,12 @@
 .\"     from: @(#)rint.3	5.1 (Berkeley) 5/2/91
 .\" $FreeBSD$
 .\"
-.Dd March 10, 1994
+.Dd July 5, 2004
 .Dt RINT 3
 .Os
 .Sh NAME
+.Nm nearbyint ,
+.Nm nearbyintf ,
 .Nm rint ,
 .Nm rintf
 .Nd round to integral value in floating-point format
@@ -44,6 +46,10 @@
 .Sh SYNOPSIS
 .In math.h
 .Ft double
+.Fn nearbyint "double x"
+.Ft float
+.Fn nearbyintf "float x"
+.Ft double
 .Fn rint "double x"
 .Ft float
 .Fn rintf "float x"
@@ -52,19 +58,44 @@ The
 .Fn rint
 and the
 .Fn rintf
-functions return the integral value (represented as a double or float
-precision number) nearest to
+functions return the integral value nearest to
 .Fa x
 according to the prevailing rounding mode.
+These functions raise an inexact exception when the original argument
+is not an exact integer.
+.Pp
+The
+.Fn nearbyint
+and
+.Fn nearbyintf
+functions perform the same operation, except that they do not raise
+an inexact exception.
 .Sh SEE ALSO
 .Xr abs 3 ,
 .Xr ceil 3 ,
 .Xr fabs 3 ,
+.Xr fenv 3 ,
 .Xr floor 3 ,
 .Xr ieee 3 ,
-.Xr math 3
+.Xr math 3 ,
+.Xr round 3
+.Sh STANDARDS
+The
+.Fn nearbyint ,
+.Fn nearbyintf ,
+.Fn rint ,
+and
+.Fn rintf
+functions conform to
+.St -isoC-99 .
 .Sh HISTORY
 A
 .Fn rint
 function appeared in
 .At v6 .
+The
+.Fn nearbyint
+and
+.Fn nearbyintf
+functions appeared in
+.Fx 5.3 .
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index f666039..2da3de7 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -215,6 +215,7 @@ int	ilogb(double);
 double	lgamma(double);
 double	log1p(double) __pure2;
 double	logb(double) __pure2;
+double	nearbyint(double) __pure2;
 double	nextafter(double, double);
 double	remainder(double, double);
 double	rint(double) __pure2;
@@ -311,6 +312,7 @@ float	atanhf(float);
 float	cbrtf(float) __pure2;
 float	logbf(float) __pure2;
 float	copysignf(float, float) __pure2;
+float	nearbyintf(float) __pure2;
 float	nextafterf(float, float);
 float	remainderf(float, float);
 float	rintf(float);
diff --git a/lib/msun/src/s_nearbyint.c b/lib/msun/src/s_nearbyint.c
new file mode 100644
index 0000000..a075d2f
--- /dev/null
+++ b/lib/msun/src/s_nearbyint.c
@@ -0,0 +1,54 @@
+/*-
+ * Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <fenv.h>
+#include <math.h>
+
+/*
+ * We save and restore the floating-point environment to avoid raising
+ * an inexact exception.  We can get away with using fesetenv()
+ * instead of feclearexcept()/feupdateenv() to restore the environment
+ * because the only exception defined for rint() is overflow, and
+ * rounding can't overflow as long as emax >= p.
+ */
+#define	DECL(type, fn, rint)	\
+type				\
+fn(type x)			\
+{				\
+	type ret;		\
+	fenv_t env;		\
+				\
+	fegetenv(&env);		\
+	ret = rint(x);		\
+	fesetenv(&env);		\
+	return (ret);		\
+}
+
+DECL(double, nearbyint, rint)
+DECL(float, nearbyintf, rintf)