Replace lrint(3), lrintf(3), llrint(3) and llrintf(3) implementations with

the same implementation that we're already using for lrintl(3) and llrintl(3).
The old implementations were derived from code from NetBSD that didn't pass
the lib/libm/msun/lrint_test regress test.  NetBSD replaced their
implementation with the FreeBSD implementation of this code which we were
already using for lrintl(3) and llrintl(3).  This fixes the regress test.

ok bluhm@, millert@

5 files changed, 72 insertions, 151 deletions

diff --git a/lib/libm/src/s_llrint.c b/lib/libm/src/s_llrint.c
index a8d74065b71..85176640490 100644
--- a/lib/libm/src/s_llrint.c
+++ b/lib/libm/src/s_llrint.c
@@ -1,14 +1,12 @@
-/*	$OpenBSD: s_llrint.c,v 1.6 2016/09/12 19:47:02 guenther Exp $	*/
-/* $NetBSD: llrint.c,v 1.2 2004/10/13 15:18:32 drochner Exp $ */
+/*	$OpenBSD: s_llrint.c,v 1.7 2021/10/14 21:30:00 kettenis Exp $	*/
 
 /*
- * Written by Matthias Drochner <drochner@NetBSD.org>.
- * Public domain.
+ * Written by Martynas Venckus.  Public domain
  */
 
-#define LRINTNAME llrint
-#define RESTYPE long long int
-#define RESTYPE_MIN LLONG_MIN
-#define RESTYPE_MAX LLONG_MAX
+#define type		double
+#define roundit		rint
+#define dtype		long long
+#define fn		llrint
 
 #include "s_lrint.c"
diff --git a/lib/libm/src/s_llrintf.c b/lib/libm/src/s_llrintf.c
index fdf763693d6..af69ba3fe48 100644
--- a/lib/libm/src/s_llrintf.c
+++ b/lib/libm/src/s_llrintf.c
@@ -1,14 +1,12 @@
-/*	$OpenBSD: s_llrintf.c,v 1.2 2006/09/25 22:16:48 kettenis Exp $	*/
-/* $NetBSD: llrintf.c,v 1.2 2004/10/13 15:18:32 drochner Exp $ */
+/*	$OpenBSD: s_llrintf.c,v 1.3 2021/10/14 21:30:00 kettenis Exp $	*/
 
 /*
- * Written by Matthias Drochner <drochner@NetBSD.org>.
- * Public domain.
+ * Written by Martynas Venckus.  Public domain
  */
 
-#define LRINTNAME llrintf
-#define RESTYPE long long int
-#define RESTYPE_MIN LLONG_MIN
-#define RESTYPE_MAX LLONG_MAX
+#define type		float
+#define roundit		rintf
+#define dtype		long long
+#define fn		llrintf
 
 #include "s_lrintf.c"
diff --git a/lib/libm/src/s_lrint.c b/lib/libm/src/s_lrint.c
index 865cb0e1b40..542ff1da7b7 100644
--- a/lib/libm/src/s_lrint.c
+++ b/lib/libm/src/s_lrint.c
@@ -1,9 +1,8 @@
-/*	$OpenBSD: s_lrint.c,v 1.11 2016/09/12 19:47:02 guenther Exp $	*/
-/* $NetBSD: lrint.c,v 1.3 2004/10/13 15:18:32 drochner Exp $ */
+/*	$OpenBSD: s_lrint.c,v 1.12 2021/10/14 21:30:00 kettenis Exp $	*/
 
 /*-
- * Copyright (c) 2004
- *	Matthias Drochner. All rights reserved.
+ * Copyright (c) 2005 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
@@ -27,75 +26,35 @@
  * SUCH DAMAGE.
  */
 
-#include <sys/types.h>
-#include <sys/limits.h>
-#include <float.h>
+#include <fenv.h>
 #include <math.h>
-#include <ieeefp.h>
-#include <machine/ieee.h>
 
-#include "math_private.h"
-
-#ifndef LRINTNAME
-#define LRINTNAME lrint
-#define RESTYPE long int
-#define RESTYPE_MIN LONG_MIN
-#define RESTYPE_MAX LONG_MAX
+#ifndef type
+#define type		double
+#define roundit		rint
+#define dtype		long
+#define fn		lrint
 #endif
 
-#define RESTYPE_BITS (sizeof(RESTYPE) * 8)
-
-static const double
-TWO52[2]={
-  4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
- -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
-};
-
-RESTYPE
-LRINTNAME(double x)
+/*
+ * C99 says we should not raise a spurious inexact exception when an
+ * invalid exception is raised.  Unfortunately, the set of inputs
+ * that overflows depends on the rounding mode when 'dtype' has more
+ * significant bits than 'type'.  Hence, we bend over backwards for the
+ * sake of correctness; an MD implementation could be more efficient.
+ */
+dtype
+fn(type x)
 {
-	u_int32_t i0, i1;
-	int e, s, shift;
-	RESTYPE res;
-
-	GET_HIGH_WORD(i0, x);
-	e = i0 >> DBL_FRACHBITS;
-	s = e >> DBL_EXPBITS;
-	e = (e & 0x7ff) - DBL_EXP_BIAS;
-
-	/* 1.0 x 2^31 (or 2^63) is already too large */
-	if (e >= (int)RESTYPE_BITS - 1)
-		return (s ? RESTYPE_MIN : RESTYPE_MAX); /* ??? unspecified */
-
-	/* >= 2^52 is already an exact integer */
-	if (e < DBL_FRACBITS) {
-		volatile double t = x;	/* clip extra precision */
-		/* round, using current direction */
-		t += TWO52[s];
-		t -= TWO52[s];
-		x = t;
-	}
-
-	EXTRACT_WORDS(i0, i1, x);
-	e = ((i0 >> DBL_FRACHBITS) & 0x7ff) - DBL_EXP_BIAS;
-	i0 &= 0xfffff;
-	i0 |= (1 << DBL_FRACHBITS);
-
-	if (e < 0)
-		return (0);
-
-	shift = e - DBL_FRACBITS;
-	if (shift >=0)
-		res = (shift < RESTYPE_BITS ? (RESTYPE)i1 << shift : 0);
-	else
-		res = (shift > -RESTYPE_BITS ? (RESTYPE)i1 >> -shift : 0);
-	shift += 32;
-	if (shift >=0)
-		res |= (shift < RESTYPE_BITS ? (RESTYPE)i0 << shift : 0);
-	else
-		res |= (shift > -RESTYPE_BITS ? (RESTYPE)i0 >> -shift : 0);
-
-	return (s ? -res : res);
+	fenv_t env;
+	dtype d;
+
+	feholdexcept(&env);
+	d = (dtype)roundit(x);
+	if (fetestexcept(FE_INVALID))
+		feclearexcept(FE_INEXACT);
+	feupdateenv(&env);
+	return (d);
 }
-DEF_STD(LRINTNAME);
-LDBL_MAYBE_CLONE(LRINTNAME);
+DEF_STD(fn);
+LDBL_MAYBE_CLONE(fn);
diff --git a/lib/libm/src/s_lrintf.c b/lib/libm/src/s_lrintf.c
index fae445483e6..a57ee8370bc 100644
--- a/lib/libm/src/s_lrintf.c
+++ b/lib/libm/src/s_lrintf.c
@@ -1,9 +1,8 @@
-/*	$OpenBSD: s_lrintf.c,v 1.6 2016/09/12 19:47:02 guenther Exp $	*/
-/* $NetBSD: lrintf.c,v 1.3 2004/10/13 15:18:32 drochner Exp $ */
+/*	$OpenBSD: s_lrintf.c,v 1.7 2021/10/14 21:30:00 kettenis Exp $	*/
 
 /*-
- * Copyright (c) 2004
- *	Matthias Drochner. All rights reserved.
+ * Copyright (c) 2005 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
@@ -27,67 +26,34 @@
  * SUCH DAMAGE.
  */
 
-#include <sys/types.h>
-#include <sys/limits.h>
+#include <fenv.h>
 #include <math.h>
-#include <ieeefp.h>
-#include <machine/ieee.h>
-#include "math_private.h"
 
-#ifndef LRINTNAME
-#define LRINTNAME lrintf
-#define RESTYPE long int
-#define RESTYPE_MIN LONG_MIN
-#define RESTYPE_MAX LONG_MAX
+#ifndef type
+#define type		float
+#define roundit		rintf
+#define dtype		long
+#define fn		lrintf
 #endif
 
-#define RESTYPE_BITS (sizeof(RESTYPE) * 8)
-
-static const float
-TWO23[2]={
-  8.3886080000e+06, /* 0x4b000000 */
- -8.3886080000e+06, /* 0xcb000000 */
-};
-
-RESTYPE
-LRINTNAME(float x)
+/*
+ * C99 says we should not raise a spurious inexact exception when an
+ * invalid exception is raised.  Unfortunately, the set of inputs
+ * that overflows depends on the rounding mode when 'dtype' has more
+ * significant bits than 'type'.  Hence, we bend over backwards for the
+ * sake of correctness; an MD implementation could be more efficient.
+ */
+dtype
+fn(type x)
 {
-	u_int32_t i0;
-	int e, s, shift;
-	RESTYPE res;
-
-	GET_FLOAT_WORD(i0, x);
-	e = i0 >> SNG_FRACBITS;
-	s = e >> SNG_EXPBITS;
-	e = (e & 0xff) - SNG_EXP_BIAS;
-
-	/* 1.0 x 2^31 (or 2^63) is already too large */
-	if (e >= (int)RESTYPE_BITS - 1)
-		return (s ? RESTYPE_MIN : RESTYPE_MAX); /* ??? unspecified */
-
-	/* >= 2^23 is already an exact integer */
-	if (e < SNG_FRACBITS) {
-		volatile float t = x;	/* clip extra precision */
-		/* round, using current direction */
-		t += TWO23[s];
-		t -= TWO23[s];
-		x = t;
-	}
-
-	GET_FLOAT_WORD(i0, x);
-	e = ((i0 >> SNG_FRACBITS) & 0xff) - SNG_EXP_BIAS;
-	i0 &= 0x7fffff;
-	i0 |= (1 << SNG_FRACBITS);
-
-	if (e < 0)
-		return (0);
-
-	shift = e - SNG_FRACBITS;
-	if (shift >=0)
-		res = (shift < RESTYPE_BITS ? (RESTYPE)i0 << shift : 0);
-	else
-		res = (shift > -RESTYPE_BITS ? (RESTYPE)i0 >> -shift : 0);
-
-	return (s ? -res : res);
+	fenv_t env;
+	dtype d;
+
+	feholdexcept(&env);
+	d = (dtype)roundit(x);
+	if (fetestexcept(FE_INVALID))
+		feclearexcept(FE_INEXACT);
+	feupdateenv(&env);
+	return (d);
 }
-DEF_STD(LRINTNAME);
+DEF_STD(fn);
diff --git a/lib/libm/src/s_lrintl.c b/lib/libm/src/s_lrintl.c
index 62b3155c2b5..d317b30fc14 100644
--- a/lib/libm/src/s_lrintl.c
+++ b/lib/libm/src/s_lrintl.c
@@ -1,4 +1,4 @@
-/*	$OpenBSD: s_lrintl.c,v 1.3 2019/03/15 05:42:38 kevlo Exp $	*/
+/*	$OpenBSD: s_lrintl.c,v 1.4 2021/10/14 21:30:00 kettenis Exp $	*/
 
 /*-
  * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
@@ -31,9 +31,9 @@
 
 #ifndef type
 #define type		long double
-#define	roundit		rintl
+#define roundit		rintl
 #define dtype		long
-#define	fn		lrintl
+#define fn		lrintl
 #endif
 
 /*