Replace lrint(3), lrintf(3), llrint(3) and llrintf(3) implementations with
authorkettenis <kettenis@openbsd.org>
Thu, 14 Oct 2021 21:30:00 +0000 (21:30 +0000)
committerkettenis <kettenis@openbsd.org>
Thu, 14 Oct 2021 21:30:00 +0000 (21:30 +0000)
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@

lib/libm/src/s_llrint.c
lib/libm/src/s_llrintf.c
lib/libm/src/s_lrint.c
lib/libm/src/s_lrintf.c
lib/libm/src/s_lrintl.c

index a8d7406..8517664 100644 (file)
@@ -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"
index fdf7636..af69ba3 100644 (file)
@@ -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"
index 865cb0e..542ff1d 100644 (file)
@@ -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
  * 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);
index fae4454..a57ee83 100644 (file)
@@ -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
  * 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);
index 62b3155..d317b30 100644 (file)
@@ -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
 
 /*