Module Name: src
Committed By: joerg
Date: Tue Nov 19 19:24:34 UTC 2013
Modified Files:
src/distrib/sets/lists/comp: mi
src/lib/libm: Makefile
src/lib/libm/man: sqrt.3
src/lib/libm/src: k_standard.c math_private.h namespace.h s_cbrt.c
w_sqrt.c
src/tests/lib/libm: t_cbrt.c t_sqrt.c
Added Files:
src/lib/libm/src: e_sqrtl.c s_cbrtl.c w_sqrtl.c
Log Message:
Add cbrtl(3) and sqrtl(3), from FreeBSD.
To generate a diff of this commit:
cvs rdiff -u -r1.1861 -r1.1862 src/distrib/sets/lists/comp/mi
cvs rdiff -u -r1.149 -r1.150 src/lib/libm/Makefile
cvs rdiff -u -r1.13 -r1.14 src/lib/libm/man/sqrt.3
cvs rdiff -u -r0 -r1.1 src/lib/libm/src/e_sqrtl.c src/lib/libm/src/s_cbrtl.c \
src/lib/libm/src/w_sqrtl.c
cvs rdiff -u -r1.18 -r1.19 src/lib/libm/src/k_standard.c
cvs rdiff -u -r1.19 -r1.20 src/lib/libm/src/math_private.h
cvs rdiff -u -r1.9 -r1.10 src/lib/libm/src/namespace.h \
src/lib/libm/src/w_sqrt.c
cvs rdiff -u -r1.11 -r1.12 src/lib/libm/src/s_cbrt.c
cvs rdiff -u -r1.1 -r1.2 src/tests/lib/libm/t_cbrt.c
cvs rdiff -u -r1.3 -r1.4 src/tests/lib/libm/t_sqrt.c
Please note that diffs are not public domain; they are subject to the
copyright notices on the relevant files.
Modified files:
Index: src/distrib/sets/lists/comp/mi
diff -u src/distrib/sets/lists/comp/mi:1.1861 src/distrib/sets/lists/comp/mi:1.1862
--- src/distrib/sets/lists/comp/mi:1.1861 Sat Nov 16 13:01:38 2013
+++ src/distrib/sets/lists/comp/mi Tue Nov 19 19:24:34 2013
@@ -1,4 +1,4 @@
-# $NetBSD: mi,v 1.1861 2013/11/16 13:01:38 alnsn Exp $
+# $NetBSD: mi,v 1.1862 2013/11/19 19:24:34 joerg Exp $
#
# Note: don't delete entries from here - mark them as "obsolete" instead.
#
@@ -5666,6 +5666,7 @@
./usr/share/man/cat3/cbreak.0 comp-c-catman .cat
./usr/share/man/cat3/cbrt.0 comp-c-catman .cat
./usr/share/man/cat3/cbrtf.0 comp-c-catman .cat
+./usr/share/man/cat3/cbrtl.0 comp-c-catman .cat
./usr/share/man/cat3/ccos.0 comp-c-catman complex,.cat
./usr/share/man/cat3/ccosf.0 comp-c-catman complex,.cat
./usr/share/man/cat3/ccosh.0 comp-c-catman complex,.cat
@@ -8804,6 +8805,7 @@
./usr/share/man/cat3/sprintf.0 comp-c-catman .cat
./usr/share/man/cat3/sqrt.0 comp-c-catman .cat
./usr/share/man/cat3/sqrtf.0 comp-c-catman .cat
+./usr/share/man/cat3/sqrtl.0 comp-c-catman .cat
./usr/share/man/cat3/sradixsort.0 comp-c-catman .cat
./usr/share/man/cat3/srand.0 comp-c-catman .cat
./usr/share/man/cat3/srand48.0 comp-c-catman .cat
@@ -12269,6 +12271,7 @@
./usr/share/man/html3/cbreak.html comp-c-htmlman html
./usr/share/man/html3/cbrt.html comp-c-htmlman html
./usr/share/man/html3/cbrtf.html comp-c-htmlman html
+./usr/share/man/html3/cbrtl.html comp-c-htmlman html
./usr/share/man/html3/ccos.html comp-c-htmlman complex,html
./usr/share/man/html3/ccosf.html comp-c-htmlman complex,html
./usr/share/man/html3/ccosh.html comp-c-htmlman complex,html
@@ -15299,6 +15302,7 @@
./usr/share/man/html3/sprintf.html comp-c-htmlman html
./usr/share/man/html3/sqrt.html comp-c-htmlman html
./usr/share/man/html3/sqrtf.html comp-c-htmlman html
+./usr/share/man/html3/sqrtl.html comp-c-htmlman html
./usr/share/man/html3/sradixsort.html comp-c-htmlman html
./usr/share/man/html3/srand.html comp-c-htmlman html
./usr/share/man/html3/srand48.html comp-c-htmlman html
@@ -18696,6 +18700,7 @@
./usr/share/man/man3/cbreak.3 comp-c-man .man
./usr/share/man/man3/cbrt.3 comp-c-man .man
./usr/share/man/man3/cbrtf.3 comp-c-man .man
+./usr/share/man/man3/cbrtl.3 comp-c-man .man
./usr/share/man/man3/ccos.3 comp-c-man complex,.man
./usr/share/man/man3/ccosf.3 comp-c-man complex,.man
./usr/share/man/man3/ccosh.3 comp-c-man complex,.man
@@ -21830,6 +21835,7 @@
./usr/share/man/man3/sprintf.3 comp-c-man .man
./usr/share/man/man3/sqrt.3 comp-c-man .man
./usr/share/man/man3/sqrtf.3 comp-c-man .man
+./usr/share/man/man3/sqrtl.3 comp-c-man .man
./usr/share/man/man3/sradixsort.3 comp-c-man .man
./usr/share/man/man3/srand.3 comp-c-man .man
./usr/share/man/man3/srand48.3 comp-c-man .man
Index: src/lib/libm/Makefile
diff -u src/lib/libm/Makefile:1.149 src/lib/libm/Makefile:1.150
--- src/lib/libm/Makefile:1.149 Wed Nov 13 22:09:55 2013
+++ src/lib/libm/Makefile Tue Nov 19 19:24:33 2013
@@ -1,4 +1,4 @@
-# $NetBSD: Makefile,v 1.149 2013/11/13 22:09:55 joerg Exp $
+# $NetBSD: Makefile,v 1.150 2013/11/19 19:24:33 joerg Exp $
#
# @(#)Makefile 5.1beta 93/09/24
#
@@ -156,11 +156,11 @@ COMMON_SRCS+= b_exp.c b_log.c b_tgamma.c
e_j1.c e_j1f.c e_jn.c e_jnf.c e_lgamma_r.c e_lgammaf_r.c e_log.c \
e_log2.c e_log10.c e_log10f.c e_log2f.c e_logf.c e_pow.c e_powf.c \
e_rem_pio2.c e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c \
- e_scalbf.c e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c \
+ e_scalbf.c e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c e_sqrtl.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 \
ldbl_dummy.c \
- s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c \
+ s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_cbrt.c s_cbrtf.c s_cbrtl.c \
s_ceil.c s_ceilf.c s_ceill.c s_copysign.c s_copysignf.c s_copysignl.c \
s_cos.c s_cosf.c s_erf.c \
s_erff.c s_exp2.c s_exp2f.c s_expm1.c s_expm1f.c s_fabsf.c s_fabsl.c \
@@ -184,7 +184,7 @@ COMMON_SRCS+= b_exp.c b_log.c b_tgamma.c
w_lgammaf.c w_lgammaf_r.c w_log.c w_log10.c w_log10f.c w_log2.c \
w_log2f.c w_logf.c \
w_pow.c w_powf.c w_remainder.c w_remainderf.c w_scalb.c w_scalbf.c \
- w_sinh.c w_sinhf.c w_sqrt.c w_sqrtf.c \
+ w_sinh.c w_sinhf.c w_sqrt.c w_sqrtf.c w_sqrtl.c \
lrint.c lrintf.c llrint.c llrintf.c lround.c lroundf.c llround.c \
llroundf.c s_frexp.c s_frexpl.c s_modf.c \
s_fmax.c s_fmaxf.c s_fmaxl.c \
@@ -313,7 +313,8 @@ MLINKS+=scalbn.3 scalbnf.3 \
scalbn.3 scalbnl.3
MLINKS+=sin.3 sinf.3
MLINKS+=sin.3 sinhf.3
-MLINKS+=sqrt.3 sqrtf.3 sqrt.3 cbrt.3 sqrt.3 cbrtf.3
+MLINKS+=sqrt.3 sqrtf.3 sqrt.3 sqrtl.3 \
+ sqrt.3 cbrt.3 sqrt.3 cbrtf.3 sqrt.3 cbrtl.3
MLINKS+=tan.3 tanf.3
MLINKS+=tanh.3 tanhf.3
MLINKS+=round.3 roundf.3 \
Index: src/lib/libm/man/sqrt.3
diff -u src/lib/libm/man/sqrt.3:1.13 src/lib/libm/man/sqrt.3:1.14
--- src/lib/libm/man/sqrt.3:1.13 Thu Aug 7 16:44:49 2003
+++ src/lib/libm/man/sqrt.3 Tue Nov 19 19:24:33 2013
@@ -26,16 +26,18 @@
.\" SUCH DAMAGE.
.\"
.\" from: @(#)sqrt.3 6.4 (Berkeley) 5/6/91
-.\" $NetBSD: sqrt.3,v 1.13 2003/08/07 16:44:49 agc Exp $
+.\" $NetBSD: sqrt.3,v 1.14 2013/11/19 19:24:33 joerg Exp $
.\"
-.Dd May 6, 1991
+.Dd November 19, 2013
.Dt SQRT 3
.Os
.Sh NAME
.Nm cbrt ,
.Nm cbrtf ,
+.Nm cbrtl ,
.Nm sqrt ,
-.Nm sqrtf
+.Nm sqrtf ,
+.Nm sqrtl
.Nd cube root and square root functions
.Sh LIBRARY
.Lb libm
@@ -45,30 +47,37 @@
.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
.Fn sqrtf "float x"
+.Ft long double
+.Fn sqrtl "long double x"
.Sh DESCRIPTION
The
-.Fn cbrt
-and
+.Fn cbrt ,
.Fn cbrtf
+and
+.Fn cbrtl
functions compute
the cube root of
.Ar x .
.Pp
The
-.Fn sqrt
-and
+.Fn sqrt ,
.Fn sqrtf
+and
+.Fn sqrtl
functions compute
the non-negative square root of x.
.Sh RETURN VALUES
If x is negative,
-.Fn sqrt "x"
-and
+.Fn sqrt "x" ,
.Fn sqrtf "x"
+and
+.Fn sqrtl "x"
.\" POSIX_MODE
set the global variable
.Va errno
Index: src/lib/libm/src/k_standard.c
diff -u src/lib/libm/src/k_standard.c:1.18 src/lib/libm/src/k_standard.c:1.19
--- src/lib/libm/src/k_standard.c:1.18 Tue Nov 19 14:04:24 2013
+++ src/lib/libm/src/k_standard.c Tue Nov 19 19:24:34 2013
@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: k_standard.c,v 1.18 2013/11/19 14:04:24 joerg Exp $");
+__RCSID("$NetBSD: k_standard.c,v 1.19 2013/11/19 19:24:34 joerg Exp $");
#endif
#include "math.h"
@@ -514,9 +514,15 @@ __kernel_standard(double x, double y, in
break;
case 26:
case 126:
+ case 226:
/* sqrt(x<0) */
exc.type = DOMAIN;
- exc.name = type < 100 ? "sqrt" : "sqrtf";
+ if (type == 26)
+ exc.name = "sqrt";
+ else if (type == 126)
+ exc.name = "sqrtf";
+ else
+ exc.name = "sqrtl";
if (_LIB_VERSION == _SVID_)
exc.retval = zero;
else
Index: src/lib/libm/src/math_private.h
diff -u src/lib/libm/src/math_private.h:1.19 src/lib/libm/src/math_private.h:1.20
--- src/lib/libm/src/math_private.h:1.19 Tue Nov 12 16:48:39 2013
+++ src/lib/libm/src/math_private.h Tue Nov 19 19:24:34 2013
@@ -11,7 +11,7 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
- * $NetBSD: math_private.h,v 1.19 2013/11/12 16:48:39 joerg Exp $
+ * $NetBSD: math_private.h,v 1.20 2013/11/19 19:24:34 joerg Exp $
*/
#ifndef _MATH_PRIVATE_H_
@@ -301,6 +301,7 @@ extern int __kernel_rem_pio2f __P((flo
/* ieee style elementary long double functions */
extern long double __ieee754_fmodl(long double, long double);
+extern long double __ieee754_sqrtl(long double);
/*
* TRUNC() is a macro that sets the trailing 27 bits in the mantissa of an
Index: src/lib/libm/src/namespace.h
diff -u src/lib/libm/src/namespace.h:1.9 src/lib/libm/src/namespace.h:1.10
--- src/lib/libm/src/namespace.h:1.9 Wed Nov 13 12:58:11 2013
+++ src/lib/libm/src/namespace.h Tue Nov 19 19:24:34 2013
@@ -1,4 +1,4 @@
-/* $NetBSD: namespace.h,v 1.9 2013/11/13 12:58:11 joerg Exp $ */
+/* $NetBSD: namespace.h,v 1.10 2013/11/19 19:24:34 joerg Exp $ */
#define atan2 _atan2
#define atan2f _atan2f
@@ -46,6 +46,8 @@
#define scalblnf _scalblnf
#define scalblnl _scalblnl
+#define sqrtl _sqrtl
+#define cbrtl _cbrtl
#define ceill _ceill
#define floorl _floorl
#define roundl _roundl
Index: src/lib/libm/src/w_sqrt.c
diff -u src/lib/libm/src/w_sqrt.c:1.9 src/lib/libm/src/w_sqrt.c:1.10
--- src/lib/libm/src/w_sqrt.c:1.9 Sun May 26 22:02:03 2002
+++ src/lib/libm/src/w_sqrt.c Tue Nov 19 19:24:34 2013
@@ -12,16 +12,22 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: w_sqrt.c,v 1.9 2002/05/26 22:02:03 wiz Exp $");
+__RCSID("$NetBSD: w_sqrt.c,v 1.10 2013/11/19 19:24:34 joerg Exp $");
#endif
/*
* wrapper sqrt(x)
*/
+#include "namespace.h"
#include "math.h"
#include "math_private.h"
+#ifndef __HAVE_LONG_DOUBLE
+__strong_alias(_sqrtl, sqrt)
+__weak_alias(sqrtl, _sqrtl)
+#endif
+
double
sqrt(double x) /* wrapper sqrt */
{
Index: src/lib/libm/src/s_cbrt.c
diff -u src/lib/libm/src/s_cbrt.c:1.11 src/lib/libm/src/s_cbrt.c:1.12
--- src/lib/libm/src/s_cbrt.c:1.11 Sun May 26 22:01:54 2002
+++ src/lib/libm/src/s_cbrt.c Tue Nov 19 19:24:34 2013
@@ -12,12 +12,18 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
-__RCSID("$NetBSD: s_cbrt.c,v 1.11 2002/05/26 22:01:54 wiz Exp $");
+__RCSID("$NetBSD: s_cbrt.c,v 1.12 2013/11/19 19:24:34 joerg Exp $");
#endif
+#include "namespace.h"
#include "math.h"
#include "math_private.h"
+#ifndef __HAVE_LONG_DOUBLE
+__strong_alias(cbrt, _cbrtl)
+__weak_alias(cbrtl, _cbrtl)
+#endif
+
/* cbrt(x)
* Return cube root of x
*/
Index: src/tests/lib/libm/t_cbrt.c
diff -u src/tests/lib/libm/t_cbrt.c:1.1 src/tests/lib/libm/t_cbrt.c:1.2
--- src/tests/lib/libm/t_cbrt.c:1.1 Sun Oct 16 08:25:40 2011
+++ src/tests/lib/libm/t_cbrt.c Tue Nov 19 19:24:33 2013
@@ -1,4 +1,4 @@
-/* $NetBSD: t_cbrt.c,v 1.1 2011/10/16 08:25:40 jruoho Exp $ */
+/* $NetBSD: t_cbrt.c,v 1.2 2013/11/19 19:24:33 joerg Exp $ */
/*-
* Copyright (c) 2011 The NetBSD Foundation, Inc.
@@ -29,7 +29,7 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
-__RCSID("$NetBSD: t_cbrt.c,v 1.1 2011/10/16 08:25:40 jruoho Exp $");
+__RCSID("$NetBSD: t_cbrt.c,v 1.2 2013/11/19 19:24:33 joerg Exp $");
#include <atf-c.h>
#include <math.h>
@@ -261,6 +261,119 @@ ATF_TC_BODY(cbrtf_zero_pos, tc)
#endif
}
+/*
+ * cbrtl(3)
+ */
+ATF_TC(cbrtl_nan);
+ATF_TC_HEAD(cbrtl_nan, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cbrtl(NaN) == NaN");
+}
+
+ATF_TC_BODY(cbrtl_nan, tc)
+{
+#ifndef __vax__
+ const long double x = 0.0L / 0.0L;
+
+ ATF_CHECK(isnan(x) != 0);
+ ATF_CHECK(isnan(cbrtl(x)) != 0);
+#endif
+}
+
+ATF_TC(cbrtl_powl);
+ATF_TC_HEAD(cbrtl_powl, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cbrtl(3) vs. powl(3)");
+}
+
+ATF_TC_BODY(cbrtl_powl, tc)
+{
+#ifndef __vax__
+ const long double x[] = { 0.0, 0.005, 1.0, 99.0, 123.123, 9999.0 };
+ const long double eps = 1.0e-15;
+ long double y, z;
+ size_t i;
+
+ for (i = 0; i < __arraycount(x); i++) {
+
+ y = cbrtl(x[i]);
+ z = powl(x[i], 1.0 / 3.0);
+
+ if (fabsl(y - z) > eps * fabsl(1 + x[i]))
+ atf_tc_fail_nonfatal("cbrtl(%0.03Lf) != "
+ "powl(%0.03Lf, 1/3)\n", x[i], x[i]);
+ }
+#endif
+}
+
+ATF_TC(cbrtl_inf_neg);
+ATF_TC_HEAD(cbrtl_inf_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cbrtl(-Inf) == -Inf");
+}
+
+ATF_TC_BODY(cbrtl_inf_neg, tc)
+{
+#ifndef __vax__
+ const long double x = -1.0L / 0.0L;
+ long double y = cbrtl(x);
+
+ ATF_CHECK(isinf(y) != 0);
+ ATF_CHECK(signbit(y) != 0);
+#endif
+}
+
+ATF_TC(cbrtl_inf_pos);
+ATF_TC_HEAD(cbrtl_inf_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cbrtl(+Inf) == +Inf");
+}
+
+ATF_TC_BODY(cbrtl_inf_pos, tc)
+{
+#ifndef __vax__
+ const long double x = 1.0L / 0.0L;
+ long double y = cbrtl(x);
+
+ ATF_CHECK(isinf(y) != 0);
+ ATF_CHECK(signbit(y) == 0);
+#endif
+}
+
+ATF_TC(cbrtl_zero_neg);
+ATF_TC_HEAD(cbrtl_zero_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cbrtl(-0.0) == -0.0");
+}
+
+ATF_TC_BODY(cbrtl_zero_neg, tc)
+{
+#ifndef __vax__
+ const long double x = -0.0L;
+ long double y = cbrtl(x);
+
+ if (fabsl(y) > 0.0 || signbit(y) == 0)
+ atf_tc_fail_nonfatal("cbrtl(-0.0) != -0.0");
+#endif
+}
+
+ATF_TC(cbrtl_zero_pos);
+ATF_TC_HEAD(cbrtl_zero_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test cbrtl(+0.0) == +0.0");
+}
+
+ATF_TC_BODY(cbrtl_zero_pos, tc)
+{
+#ifndef __vax__
+ const long double x = 0.0L;
+ long double y = cbrtl(x);
+
+ if (fabsl(y) > 0.0 || signbit(y) != 0)
+ atf_tc_fail_nonfatal("cbrtl(+0.0) != +0.0");
+#endif
+}
+
ATF_TP_ADD_TCS(tp)
{
@@ -278,5 +391,12 @@ ATF_TP_ADD_TCS(tp)
ATF_TP_ADD_TC(tp, cbrtf_zero_neg);
ATF_TP_ADD_TC(tp, cbrtf_zero_pos);
+ ATF_TP_ADD_TC(tp, cbrtl_nan);
+ ATF_TP_ADD_TC(tp, cbrtl_powl);
+ ATF_TP_ADD_TC(tp, cbrtl_inf_neg);
+ ATF_TP_ADD_TC(tp, cbrtl_inf_pos);
+ ATF_TP_ADD_TC(tp, cbrtl_zero_neg);
+ ATF_TP_ADD_TC(tp, cbrtl_zero_pos);
+
return atf_no_error();
}
Index: src/tests/lib/libm/t_sqrt.c
diff -u src/tests/lib/libm/t_sqrt.c:1.3 src/tests/lib/libm/t_sqrt.c:1.4
--- src/tests/lib/libm/t_sqrt.c:1.3 Mon Feb 13 05:09:01 2012
+++ src/tests/lib/libm/t_sqrt.c Tue Nov 19 19:24:33 2013
@@ -1,4 +1,4 @@
-/* $NetBSD: t_sqrt.c,v 1.3 2012/02/13 05:09:01 jruoho Exp $ */
+/* $NetBSD: t_sqrt.c,v 1.4 2013/11/19 19:24:33 joerg Exp $ */
/*-
* Copyright (c) 2011 The NetBSD Foundation, Inc.
@@ -29,7 +29,7 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
-__RCSID("$NetBSD: t_sqrt.c,v 1.3 2012/02/13 05:09:01 jruoho Exp $");
+__RCSID("$NetBSD: t_sqrt.c,v 1.4 2013/11/19 19:24:33 joerg Exp $");
#include <atf-c.h>
#include <math.h>
@@ -259,6 +259,118 @@ ATF_TC_BODY(sqrtf_zero_pos, tc)
#endif
}
+/*
+ * sqrtl(3)
+ */
+ATF_TC(sqrtl_nan);
+ATF_TC_HEAD(sqrtl_nan, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sqrtl(NaN) == NaN");
+}
+
+ATF_TC_BODY(sqrtl_nan, tc)
+{
+#ifndef __vax__
+ const long double x = 0.0L / 0.0L;
+
+ ATF_CHECK(isnan(x) != 0);
+ ATF_CHECK(isnan(sqrtl(x)) != 0);
+#endif
+}
+
+ATF_TC(sqrtl_powl);
+ATF_TC_HEAD(sqrtl_powl, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sqrtl(3) vs. powl(3)");
+}
+
+ATF_TC_BODY(sqrtl_powl, tc)
+{
+#ifndef __vax__
+ const long double x[] = { 0.0, 0.005, 1.0, 99.0, 123.123, 9999.9999 };
+ const long double eps = 1.0e-30;
+ volatile long double y, z;
+ size_t i;
+
+ for (i = 0; i < __arraycount(x); i++) {
+
+ y = sqrtl(x[i]);
+ z = powl(x[i], 1.0 / 2.0);
+
+ if (fabsl(y - z) > eps)
+ atf_tc_fail_nonfatal("sqrtl(%0.03Lf) != "
+ "powl(%0.03Lf, 1/2)\n", x[i], x[i]);
+ }
+#endif
+}
+
+ATF_TC(sqrtl_inf_neg);
+ATF_TC_HEAD(sqrtl_inf_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sqrtl(-Inf) == NaN");
+}
+
+ATF_TC_BODY(sqrtl_inf_neg, tc)
+{
+#ifndef __vax__
+ const long double x = -1.0L / 0.0L;
+ long double y = sqrtl(x);
+
+ ATF_CHECK(isnan(y) != 0);
+#endif
+}
+
+ATF_TC(sqrtl_inf_pos);
+ATF_TC_HEAD(sqrtl_inf_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sqrtl(+Inf) == +Inf");
+}
+
+ATF_TC_BODY(sqrtl_inf_pos, tc)
+{
+#ifndef __vax__
+ const long double x = 1.0L / 0.0L;
+ long double y = sqrtl(x);
+
+ ATF_CHECK(isinf(y) != 0);
+ ATF_CHECK(signbit(y) == 0);
+#endif
+}
+
+ATF_TC(sqrtl_zero_neg);
+ATF_TC_HEAD(sqrtl_zero_neg, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sqrtl(-0.0) == -0.0");
+}
+
+ATF_TC_BODY(sqrtl_zero_neg, tc)
+{
+#ifndef __vax__
+ const long double x = -0.0L;
+ long double y = sqrtl(x);
+
+ if (fabsl(y) > 0.0 || signbit(y) == 0)
+ atf_tc_fail_nonfatal("sqrtl(-0.0) != -0.0");
+#endif
+}
+
+ATF_TC(sqrtl_zero_pos);
+ATF_TC_HEAD(sqrtl_zero_pos, tc)
+{
+ atf_tc_set_md_var(tc, "descr", "Test sqrtl(+0.0) == +0.0");
+}
+
+ATF_TC_BODY(sqrtl_zero_pos, tc)
+{
+#ifndef __vax__
+ const long double x = 0.0L;
+ long double y = sqrtl(x);
+
+ if (fabsl(y) > 0.0 || signbit(y) != 0)
+ atf_tc_fail_nonfatal("sqrtl(+0.0) != +0.0");
+#endif
+}
+
ATF_TP_ADD_TCS(tp)
{
@@ -276,5 +388,12 @@ ATF_TP_ADD_TCS(tp)
ATF_TP_ADD_TC(tp, sqrtf_zero_neg);
ATF_TP_ADD_TC(tp, sqrtf_zero_pos);
+ ATF_TP_ADD_TC(tp, sqrtl_nan);
+ ATF_TP_ADD_TC(tp, sqrtl_powl);
+ ATF_TP_ADD_TC(tp, sqrtl_inf_neg);
+ ATF_TP_ADD_TC(tp, sqrtl_inf_pos);
+ ATF_TP_ADD_TC(tp, sqrtl_zero_neg);
+ ATF_TP_ADD_TC(tp, sqrtl_zero_pos);
+
return atf_no_error();
}
Added files:
Index: src/lib/libm/src/e_sqrtl.c
diff -u /dev/null src/lib/libm/src/e_sqrtl.c:1.1
--- /dev/null Tue Nov 19 19:24:34 2013
+++ src/lib/libm/src/e_sqrtl.c Tue Nov 19 19:24:34 2013
@@ -0,0 +1,161 @@
+/*-
+ * Copyright (c) 2007 Steven G. Kargl
+ * 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 unmodified, 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 ``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 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>
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/e_sqrtl.c 176720 2008-03-02 01:47:58Z das $");
+#endif
+__RCSID("$NetBSD: e_sqrtl.c,v 1.1 2013/11/19 19:24:34 joerg Exp $");
+
+#include <machine/ieee.h>
+#include <fenv.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __HAVE_LONG_DOUBLE
+
+/* Return (x + ulp) for normal positive x. Assumes no overflow. */
+static inline long double
+inc(long double x)
+{
+ union ieee_ext_u ux = { .extu_ld = x, };
+
+ if (++ux.extu_fracl == 0) {
+ if (++ux.extu_frach == 0) {
+ ux.extu_exp++;
+ ux.extu_frach |= LDBL_NBIT;
+ }
+ }
+ return (ux.extu_ld);
+}
+
+/* Return (x - ulp) for normal positive x. Assumes no underflow. */
+static inline long double
+dec(long double x)
+{
+ union ieee_ext_u ux = { .extu_ld = x, };
+
+ if (ux.extu_fracl-- == 0) {
+ if (ux.extu_frach-- == LDBL_NBIT) {
+ ux.extu_exp--;
+ ux.extu_frach |= LDBL_NBIT;
+ }
+ }
+ return (ux.extu_ld);
+}
+
+/*
+ * This is slow, but simple and portable. You should use hardware sqrt
+ * if possible.
+ */
+
+long double
+__ieee754_sqrtl(long double x)
+{
+ union ieee_ext_u ux = { .extu_ld = x, };
+ int k, r;
+ long double lo, xn;
+ fenv_t env;
+
+ /* If x = NaN, then sqrt(x) = NaN. */
+ /* If x = Inf, then sqrt(x) = Inf. */
+ /* If x = -Inf, then sqrt(x) = NaN. */
+ if (ux.extu_exp == LDBL_MAX_EXP * 2 - 1)
+ return (x * x + x);
+
+ /* If x = +-0, then sqrt(x) = +-0. */
+ if ((ux.extu_frach | ux.extu_fracl | ux.extu_exp) == 0)
+ return (x);
+
+ /* If x < 0, then raise invalid and return NaN */
+ if (ux.extu_sign)
+ return ((x - x) / (x - x));
+
+ feholdexcept(&env);
+
+ if (ux.extu_exp == 0) {
+ /* Adjust subnormal numbers. */
+ ux.extu_ld *= 0x1.0p514;
+ k = -514;
+ } else {
+ k = 0;
+ }
+ /*
+ * ux.extu_ld is a normal number, so break it into ux.extu_ld = e*2^n where
+ * ux.extu_ld = (2*e)*2^2k for odd n and ux.extu_ld = (4*e)*2^2k for even n.
+ */
+ if ((ux.extu_exp - EXT_EXP_BIAS) & 1) { /* n is even. */
+ k += ux.extu_exp - EXT_EXP_BIAS - 1; /* 2k = n - 2. */
+ ux.extu_exp = EXT_EXP_BIAS + 1; /* ux.extu_ld in [2,4). */
+ } else {
+ k += ux.extu_exp - EXT_EXP_BIAS; /* 2k = n - 1. */
+ ux.extu_exp = EXT_EXP_BIAS; /* ux.extu_ld in [1,2). */
+ }
+
+ /*
+ * Newton's iteration.
+ * Split ux.extu_ld into a high and low part to achieve additional precision.
+ */
+ xn = sqrt(ux.extu_ld); /* 53-bit estimate of sqrtl(x). */
+#if LDBL_MANT_DIG > 100
+ xn = (xn + (ux.extu_ld / xn)) * 0.5; /* 106-bit estimate. */
+#endif
+ lo = ux.extu_ld;
+ ux.extu_fracl = 0; /* Zero out lower bits. */
+ lo = (lo - ux.extu_ld) / xn; /* Low bits divided by xn. */
+ xn = xn + (ux.extu_ld / xn); /* High portion of estimate. */
+ ux.extu_ld = xn + lo; /* Combine everything. */
+ ux.extu_exp += (k >> 1) - 1;
+
+ feclearexcept(FE_INEXACT);
+ r = fegetround();
+ fesetround(FE_TOWARDZERO); /* Set to round-toward-zero. */
+ xn = x / ux.extu_ld; /* Chopped quotient (inexact?). */
+
+ if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */
+ if (xn == ux.extu_ld) {
+ fesetenv(&env);
+ return (ux.extu_ld);
+ }
+ /* Round correctly for inputs like x = y**2 - ulp. */
+ xn = dec(xn); /* xn = xn - ulp. */
+ }
+
+ if (r == FE_TONEAREST) {
+ xn = inc(xn); /* xn = xn + ulp. */
+ } else if (r == FE_UPWARD) {
+ ux.extu_ld = inc(ux.extu_ld); /* ux.extu_ld = ux.extu_ld + ulp. */
+ xn = inc(xn); /* xn = xn + ulp. */
+ }
+ ux.extu_ld = ux.extu_ld + xn; /* Chopped sum. */
+ feupdateenv(&env); /* Restore env and raise inexact */
+ ux.extu_exp--;
+ return (ux.extu_ld);
+}
+
+#endif
Index: src/lib/libm/src/s_cbrtl.c
diff -u /dev/null src/lib/libm/src/s_cbrtl.c:1.1
--- /dev/null Tue Nov 19 19:24:34 2013
+++ src/lib/libm/src/s_cbrtl.c Tue Nov 19 19:24:34 2013
@@ -0,0 +1,145 @@
+/*-
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ * Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
+ *
+ * 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.
+ * ====================================================
+ *
+ * The argument reduction and testing for exceptional cases was
+ * written by Steven G. Kargl with input from Bruce D. Evans
+ * and David A. Schultz.
+ */
+
+#include <sys/cdefs.h>
+__RCSID("$NetBSD: s_cbrtl.c,v 1.1 2013/11/19 19:24:34 joerg Exp $");
+#if 0
+__FBSDID("$FreeBSD: head/lib/msun/src/s_cbrtl.c 238924 2012-07-30 21:58:28Z kargl $");
+#endif
+
+#include "namespace.h"
+#include <machine/ieee.h>
+#include <float.h>
+
+#include "math.h"
+#include "math_private.h"
+
+#ifdef __HAVE_LONG_DOUBLE
+__weak_alias(cbrtl, _cbrtl)
+
+static const unsigned
+ B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
+
+long double
+cbrtl(long double x)
+{
+ union ieee_ext_u ux, vx;
+ long double r, s, t, w;
+ double dr, dt, dx;
+ float ft, fx;
+ uint32_t hx;
+ int k;
+
+ ux.extu_ld = x;
+
+
+ /*
+ * If x = +-Inf, then cbrt(x) = +-Inf.
+ * If x = NaN, then cbrt(x) = NaN.
+ */
+ if (ux.extu_exp == EXT_EXP_INFNAN)
+ return (x + x);
+ if ((ux.extu_frach | ux.extu_fracl | ux.extu_exp) == 0)
+ return (x);
+
+ vx.extu_ld = 1;
+ vx.extu_ext.ext_sign = ux.extu_ext.ext_sign;
+ ux.extu_ext.ext_sign = 0;
+ if (ux.extu_exp == 0) {
+ /* Adjust subnormal numbers. */
+ ux.extu_ld *= 0x1.0p514;
+ k = ux.extu_exp - EXT_EXP_BIAS - 514;
+ } else {
+ k = ux.extu_exp - EXT_EXP_BIAS;
+ }
+
+ ux.extu_exp = EXT_EXP_BIAS;
+ x = ux.extu_ld;
+ switch (k % 3) {
+ case 1:
+ case -2:
+ x = 2*x;
+ k--;
+ break;
+ case 2:
+ case -1:
+ x = 4*x;
+ k -= 2;
+ break;
+ }
+ vx.extu_exp = EXT_EXP_BIAS + k / 3;
+
+ /*
+ * The following is the guts of s_cbrtf, with the handling of
+ * special values removed and extra care for accuracy not taken,
+ * but with most of the extra accuracy not discarded.
+ */
+
+ /* ~5-bit estimate: */
+ fx = x;
+ GET_FLOAT_WORD(hx, fx);
+ SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
+
+ /* ~16-bit estimate: */
+ dx = x;
+ dt = ft;
+ dr = dt * dt * dt;
+ dt = dt * (dx + dx + dr) / (dx + dr + dr);
+
+ /* ~47-bit estimate: */
+ dr = dt * dt * dt;
+ dt = dt * (dx + dx + dr) / (dx + dr + dr);
+
+#if LDBL_MANT_DIG == 64
+ /*
+ * dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
+ * Round it away from zero to 32 bits (32 so that t*t is exact, and
+ * away from zero for technical reasons).
+ */
+ volatile double vd2 = 0x1.0p32;
+ volatile double vd1 = 0x1.0p-31;
+ #define vd ((long double)vd2 + vd1)
+
+ t = dt + vd - 0x1.0p32;
+#elif LDBL_MANT_DIG == 113
+ /*
+ * Round dt away from zero to 47 bits. Since we don't trust the 47,
+ * add 2 47-bit ulps instead of 1 to round up. Rounding is slow and
+ * might be avoidable in this case, since on most machines dt will
+ * have been evaluated in 53-bit precision and the technical reasons
+ * for rounding up might not apply to either case in cbrtl() since
+ * dt is much more accurate than needed.
+ */
+ t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
+#else
+#error "Unsupported long double format"
+#endif
+
+ /*
+ * Final step Newton iteration to 64 or 113 bits with
+ * error < 0.667 ulps
+ */
+ s=t*t; /* t*t is exact */
+ r=x/s; /* error <= 0.5 ulps; |r| < |t| */
+ w=t+t; /* t+t is exact */
+ r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
+ t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
+
+ t *= vx.extu_ld;
+ return t;
+}
+
+#endif /* __HAVE_LONG_DOUBLE */
Index: src/lib/libm/src/w_sqrtl.c
diff -u /dev/null src/lib/libm/src/w_sqrtl.c:1.1
--- /dev/null Tue Nov 19 19:24:34 2013
+++ src/lib/libm/src/w_sqrtl.c Tue Nov 19 19:24:34 2013
@@ -0,0 +1,40 @@
+/* @(#)w_sqrt.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#include <sys/cdefs.h>
+__RCSID("$NetBSD: w_sqrtl.c,v 1.1 2013/11/19 19:24:34 joerg Exp $");
+
+/*
+ * wrapper sqrtl(x)
+ */
+
+#include "namespace.h"
+#include "math.h"
+#include "math_private.h"
+
+__weak_alias(sqrtl, _sqrtl)
+
+long double
+sqrtl(long double x) /* wrapper sqrtl */
+{
+#ifdef _IEEE_LIBM
+ return __ieee754_sqrtl(x);
+#else
+ long double z;
+ z = __ieee754_sqrtl(x);
+ if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
+ if(x<0.0) {
+ return __kernel_standard(x,x,226); /* sqrtl(negative) */
+ } else
+ return z;
+#endif
+}
