Implement the long double version for the cube root function, cbrtl.

The algorithm uses Newton's iterations with a crude estimate of the
cube root to converge to a result.

Reviewed by:	bde
Approved by:	das

6 files changed, 32 insertions, 18 deletions

diff --git a/lib/msun/Makefile b/lib/msun/Makefile
index 2d7d5cc..bb99e45 100644
--- a/lib/msun/Makefile
+++ b/lib/msun/Makefile
@@ -93,8 +93,8 @@ COMMON_SRCS+=	s_copysignl.c s_fabsl.c s_llrintl.c s_lrintl.c s_modfl.c
 COMMON_SRCS+=	e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
 	e_hypotl.c e_remainderl.c e_sqrtl.c \
 	invtrig.c k_cosl.c k_sinl.c k_tanl.c \
-	s_atanl.c s_ceill.c s_cosl.c s_cprojl.c s_csqrtl.c s_exp2l.c \
-	s_floorl.c s_fmal.c \
+	s_atanl.c s_cbrtl.c s_ceill.c s_cosl.c s_cprojl.c \
+	s_csqrtl.c s_exp2l.c s_floorl.c s_fmal.c \
 	s_frexpl.c s_logbl.c s_nanl.c s_nextafterl.c s_nexttoward.c \
 	s_remquol.c s_rintl.c s_scalbnl.c \
 	s_sinl.c s_tanl.c s_truncl.c w_cabsl.c
@@ -197,7 +197,8 @@ MLINKS+=scalbn.3 scalbln.3 scalbn.3 scalblnf.3 scalbn.3 scalblnl.3
 MLINKS+=scalbn.3 scalbnf.3 scalbn.3 scalbnl.3
 MLINKS+=sin.3 sinf.3 sin.3 sinl.3
 MLINKS+=sinh.3 sinhf.3
-MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 sqrtf.3 sqrt.3 sqrtl.3
+MLINKS+=sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3 sqrt.3 sqrtf.3 \
+	sqrt.3 sqrtl.3
 MLINKS+=tan.3 tanf.3 tan.3 tanl.3
 MLINKS+=tanh.3 tanhf.3
 MLINKS+=trunc.3 truncf.3 trunc.3 truncl.3
diff --git a/lib/msun/Symbol.map b/lib/msun/Symbol.map
index cb60ecf..bb53ac3 100644
--- a/lib/msun/Symbol.map
+++ b/lib/msun/Symbol.map
@@ -222,6 +222,7 @@ FBSD_1.1 {
 /* First added in 9.0-CURRENT */
 FBSD_1.2 {
 	__isnanf;
+	cbrtl;
 	cexp;
 	cexpf;
 	log2;
diff --git a/lib/msun/man/math.3 b/lib/msun/man/math.3
index 39aa8ab..9ffe080 100644
--- a/lib/msun/man/math.3
+++ b/lib/msun/man/math.3
@@ -220,12 +220,11 @@ and
 .Vt long double
 values, were written for or imported into subsequent versions of FreeBSD.
 .Sh BUGS
-The
-.Fn cbrt
-function and many of the transcendental functions
-are not available in their
+Some of the
 .Vt "long double"
-variants.
+math functions in
+.St -isoC-99
+are not available.
 .Pp
 Many of the routines to compute transcendental functions produce
 inaccurate results in other than the default rounding mode.
diff --git a/lib/msun/man/sqrt.3 b/lib/msun/man/sqrt.3
index 71460fa..e48f98b 100644
--- a/lib/msun/man/sqrt.3
+++ b/lib/msun/man/sqrt.3
@@ -28,12 +28,13 @@
 .\"     from: @(#)sqrt.3	6.4 (Berkeley) 5/6/91
 .\" $FreeBSD$
 .\"
-.Dd March 1, 2008
+.Dd March 5, 2011
 .Dt SQRT 3
 .Os
 .Sh NAME
 .Nm cbrt ,
 .Nm cbrtf ,
+.Nm cbrtl ,
 .Nm sqrt ,
 .Nm sqrtf ,
 .Nm sqrtl
@@ -46,6 +47,8 @@
 .Fn cbrt "double x"
 .Ft float
 .Fn cbrtf "float x"
+.Ft long double
+.Fn cbrtl "long double x"
 .Ft double
 .Fn sqrt "double x"
 .Ft float
@@ -54,9 +57,10 @@
 .Fn sqrtl "long double x"
 .Sh DESCRIPTION
 The
-.Fn cbrt
-and the
-.Fn cbrtf
+.Fn cbrt ,
+.Fn cbrtf ,
+and
+.Fn cbrtl
 functions compute
 the cube root of
 .Ar x .
@@ -67,12 +71,14 @@ The
 and
 .Fn sqrtl
 functions compute the
-non-negative square root of x.
+non-negative square root of
+.Ar x .
 .Sh RETURN VALUES
 The
-.Fn cbrt
-and the
-.Fn cbrtf
+.Fn cbrt ,
+.Fn cbrtf ,
+and
+.Fn cbrtl
 functions return the requested cube root.
 The
 .Fn sqrt ,
@@ -94,6 +100,7 @@ raises an invalid exception and causes an \*(Na to be returned
 The
 .Fn cbrt ,
 .Fn cbrtf ,
+.Fn cbrtl ,
 .Fn sqrt ,
 .Fn sqrtf ,
 and
@@ -109,3 +116,7 @@ The
 .Fn sqrtl
 function appeared in
 .Fx 8.0 .
+The
+.Fn cbrtl
+function appeared in
+.Fx 9.0 .
diff --git a/lib/msun/src/math.h b/lib/msun/src/math.h
index 9aa2484..8ad13ed 100644
--- a/lib/msun/src/math.h
+++ b/lib/msun/src/math.h
@@ -411,9 +411,7 @@ long double	atan2l(long double, long double);
 long double	atanhl(long double);
 #endif
 long double	atanl(long double);
-#if 0
 long double	cbrtl(long double);
-#endif
 long double	ceill(long double);
 long double	copysignl(long double, long double) __pure2;
 #if 0
diff --git a/lib/msun/src/s_cbrt.c b/lib/msun/src/s_cbrt.c
index 27ecd97..910f75b 100644
--- a/lib/msun/src/s_cbrt.c
+++ b/lib/msun/src/s_cbrt.c
@@ -111,3 +111,7 @@ cbrt(double x)
 
 	return(t);
 }
+
+#if (LDBL_MANT_DIG == 53)
+__weak_reference(cbrt, cbrtl);
+#endif