Implement and document ceill().

4 files changed, 117 insertions, 8 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 61dfac0..07078af 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -80,8 +80,9 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
 	e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \
 	k_cos.c k_cosf.c k_rem_pio2.c k_rem_pio2f.c k_sin.c k_sinf.c \
 	k_standard.c k_tan.c k_tanf.c \
-	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c s_ceil.c \
-	s_ceilf.c s_copysign.c s_copysignf.c s_cos.c s_cosf.c s_erf.c s_erff.c \
+	s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c \
+	s_ceil.c s_ceilf.c s_ceill.c \
+	s_copysign.c s_copysignf.c s_cos.c s_cosf.c s_erf.c s_erff.c \
 	s_expm1.c s_expm1f.c s_fabsf.c s_fdim.c s_finite.c s_finitef.c \
 	s_floor.c s_floorf.c s_floorl.c s_fmax.c s_fmaxf.c s_fmaxl.c s_fmin.c \
 	s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
@@ -143,7 +144,7 @@ MLINKS+=asinh.3 asinhf.3
 MLINKS+=atan.3 atanf.3
 MLINKS+=atanh.3 atanhf.3
 MLINKS+=atan2.3 atan2f.3
-MLINKS+=ceil.3 ceilf.3
+MLINKS+=ceil.3 ceilf.3 ceil.3 ceill.3
 MLINKS+=cimag.3 cimagf.3 cimag.3 cimagl.3 \
 	cimag.3 conj.3 cimag.3 conjf.3 cimag.3 conjl.3 \
 	cimag.3 creal.3 cimag.3 crealf.3 cimag.3 creall.3
diff --git a/lib/msun/man/ceil.3 b/lib/msun/man/ceil.3
index 15b0ff4..153624e 100644
--- a/lib/msun/man/ceil.3
+++ b/lib/msun/man/ceil.3
@@ -32,12 +32,13 @@
 .\"     from: @(#)ceil.3	5.1 (Berkeley) 5/2/91
 .\" $FreeBSD$
 .\"
-.Dd March 10, 1994
+.Dd January 13, 2005
 .Dt CEIL 3
 .Os
 .Sh NAME
 .Nm ceil ,
-.Nm ceilf
+.Nm ceilf ,
+.Nm ceill
 .Nd smallest integral value greater than or equal to x
 .Sh LIBRARY
 .Lb libm
@@ -47,11 +48,15 @@
 .Fn ceil "double x"
 .Ft float
 .Fn ceilf "float x"
+.Ft long double
+.Fn ceill "long double x"
 .Sh DESCRIPTION
 The
-.Fn ceil
-and the
+.Fn ceil ,
+the
 .Fn ceilf
+and the
+.Fn ceill
 functions return the smallest integral value
 greater than or equal to
 .Fa x ,
@@ -70,3 +75,9 @@ The
 .Fn ceil
 function conforms to
 .St -isoC .
+The
+.Fn ceilf
+and the
+.Fn ceill
+functions conform to
+.St -isoC-99 .
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 5b6ea73..841402d 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -431,8 +431,8 @@ long double	atan2l(long double, long double);
 long double	atanhl(long double);
 long double	atanl(long double);
 long double	cbrtl(long double);
-long double	ceill(long double);
 #endif
+long double	ceill(long double);
 long double	copysignl(long double, long double);
 #if 0
 long double	coshl(long double);
diff --git a/lib/msun/src/s_ceill.c b/lib/msun/src/s_ceill.c
new file mode 100644
index 0000000..92f5033
--- /dev/null
+++ b/lib/msun/src/s_ceill.c
@@ -0,0 +1,97 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ *
+ * From: @(#)s_ceil.c 5.1 93/09/24
+ */
+
+#ifndef lint
+static char rcsid[] = "$FreeBSD$";
+#endif
+
+/*
+ * ceill(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ *	Bit twiddling.
+ * Exception:
+ *	Inexact flag raised if x not equal to ceill(x).
+ */
+
+#include <float.h>
+#include <math.h>
+#include <stdint.h>
+
+#include "fpmath.h"
+
+#ifdef LDBL_IMPLICIT_NBIT
+#define	MANH_SIZE	(LDBL_MANH_SIZE + 1)
+#define	INC_MANH(u, c)	do {					\
+	uint64_t o = u.bits.manh;				\
+	u.bits.manh += (c);					\
+	if (u.bits.manh < o)					\
+		u.bits.exp++;					\
+} while (0)
+#else
+#define	MANH_SIZE	LDBL_MANH_SIZE
+#define	INC_MANH(u, c)	do {					\
+	uint64_t o = u.bits.manh;				\
+	u.bits.manh += (c);					\
+	if (u.bits.manh < o) {					\
+		u.bits.exp++;					\
+		u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1);	\
+	}							\
+} while (0)
+#endif
+
+long double
+ceill(long double x)
+{
+	union IEEEl2bits u = { .e = x };
+	int e = u.bits.exp - LDBL_MAX_EXP + 1;
+
+	if (e < MANH_SIZE - 1) {
+		if (e < 0) {			/* raise inexact if x != 0 */
+			if (u.bits.exp > 0 || (u.bits.manh | u.bits.manl) != 0)
+				u.e = u.bits.sign ? 0.0 : 1.0;
+		} else {
+			uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
+			if (((u.bits.manh & m) | u.bits.manl) == 0)
+				return (x);	/* x is integral */
+			/* raise inexact flag */
+			if (!u.bits.sign) {
+#ifdef LDBL_IMPLICIT_NBIT
+				if (e == 0)
+					u.bits.exp++;
+				else
+#endif
+				INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
+			}
+			u.bits.manh &= ~m;
+			u.bits.manl = 0;
+		}
+	} else if (e < LDBL_MANT_DIG - 1) {
+		uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
+		if ((u.bits.manl & m) == 0)
+			return (x);	/* x is integral */
+		/* raise inexact flag */
+		if (!u.bits.sign) {
+			if (e == MANH_SIZE - 1)
+				INC_MANH(u, 1);
+			else {
+				uint64_t o = u.bits.manl;
+				u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
+				if (u.bits.manl < o)	/* got a carry */
+					INC_MANH(u, 1);
+			}
+		}
+		u.bits.manl &= ~m;
+	}
+	return (u.e);
+}