standards/82654: C99 long double math functions are missing
Steven G. Kargl
kargl at troutmask.apl.washington.edu
Sat Jun 25 23:40:12 GMT 2005
>Number: 82654
>Category: standards
>Synopsis: C99 long double math functions are missing
>Confidential: no
>Severity: serious
>Priority: medium
>Responsible: freebsd-standards
>State: open
>Quarter:
>Keywords:
>Date-Required:
>Class: sw-bug
>Submitter-Id: current-users
>Arrival-Date: Sat Jun 25 23:40:12 GMT 2005
>Closed-Date:
>Last-Modified:
>Originator: Steven G. Kargl
>Release: FreeBSD 6.0-CURRENT amd64
>Organization:
APL/UW
>Environment:
System: FreeBSD troutmask.apl.washington.edu 6.0-CURRENT FreeBSD 6.0-CURRENT #1: Thu Jun 16 15:47:33 PDT 2005 kargl at troutmask.apl.washington.edu:/usr/obj/usr/src/sys/SPEW amd64
>Description:
Most of the long double math functions as prescribed by C99 are
missing.
>How-To-Repeat:
>Fix:
The enclosed patch implements logl(), log10l(), sqrtl(), and cbrtl().
I'm sure someone will want bit twiddling or assembly code, but the
c code works on both i386 and amd64.
diff -rNu msun.orig/Makefile msun/Makefile
--- msun.orig/Makefile Sat Jun 25 16:02:58 2005
+++ msun/Makefile Sat Jun 25 16:07:48 2005
@@ -38,7 +38,7 @@
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_tan.c k_tanf.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_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_fdim.c \
@@ -48,13 +48,15 @@
s_fminf.c s_fminl.c s_frexp.c s_frexpf.c s_ilogb.c s_ilogbf.c \
s_ilogbl.c s_isfinite.c s_isnan.c s_isnormal.c \
s_llrint.c s_llrintf.c s_llround.c s_llroundf.c s_llroundl.c \
- s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_lrint.c s_lrintf.c \
+ s_log1p.c s_log1pf.c s_logb.c s_logbf.c s_logl.c s_log10l.c \
+ s_lrint.c s_lrintf.c \
s_lround.c s_lroundf.c s_lroundl.c s_modff.c \
s_nearbyint.c s_nextafter.c s_nextafterf.c \
s_nexttowardf.c s_remquo.c s_remquof.c \
s_rint.c s_rintf.c s_round.c s_roundf.c s_roundl.c \
s_scalbln.c s_scalbn.c s_scalbnf.c s_signbit.c \
- s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c s_tan.c \
+ s_signgam.c s_significand.c s_significandf.c s_sin.c s_sinf.c \
+ s_sqrtl.c s_tan.c \
s_tanf.c s_tanh.c s_tanhf.c s_trunc.c s_truncf.c s_truncl.c \
w_cabs.c w_cabsf.c w_drem.c w_dremf.c
diff -rNu msun.orig/man/exp.3 msun/man/exp.3
--- msun.orig/man/exp.3 Sat Jun 25 16:02:58 2005
+++ msun/man/exp.3 Sat Jun 25 16:03:53 2005
@@ -45,8 +45,10 @@
.Nm expm1f ,
.Nm log ,
.Nm logf ,
+.Nm logl ,
.Nm log10 ,
.Nm log10f ,
+.Nm log10l ,
.Nm log1p ,
.Nm log1pf ,
.Nm pow ,
@@ -72,10 +74,14 @@
.Fn log "double x"
.Ft float
.Fn logf "float x"
+.Ft "long double"
+.Fn logl "long double x"
.Ft double
.Fn log10 "double x"
.Ft float
.Fn log10f "float x"
+.Ft "long double"
+.Fn log10l "long double x"
.Ft double
.Fn log1p "double x"
.Ft float
@@ -109,16 +115,18 @@
.Fa x .
.Pp
The
-.Fn log
-and the
-.Fn logf
+.Fn log ,
+.Fn logf ,
+and
+.Fn logl
functions compute the value of the natural logarithm of argument
.Fa x .
.Pp
The
-.Fn log10
-and the
-.Fn log10f
+.Fn log10 ,
+.Fn log10f ,
+and
+.Fn log10l
functions compute the value of the logarithm of argument
.Fa x
to base 10.
diff -rNu msun.orig/man/sqrt.3 msun/man/sqrt.3
--- msun.orig/man/sqrt.3 Sat Jun 25 16:02:58 2005
+++ msun/man/sqrt.3 Sat Jun 25 16:03:59 2005
@@ -49,35 +49,43 @@
.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 the
-.Fn cbrtf
+.Fn cbrt ,
+.Fn cbrtf ,
+and
+.Fn cbrtl
functions compute
the cube root of
.Ar x .
.Pp
The
-.Fn sqrt
-and the
-.Fn sqrtf
-functions compute the
-non-negative square root of x.
+.Fn sqrt ,
+.Fn sqrtf ,
+and
+.Fn sqrtl
+functions compute the 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
-and the
-.Fn sqrtf
+.Fn sqrt ,
+.Fn sqrtf ,
+and
+.Fn sqrtl
functions return the requested square root
unless an error occurs.
An attempt to take the
diff -rNu msun.orig/src/math.h msun/src/math.h
--- msun.orig/src/math.h Sat Jun 25 16:02:58 2005
+++ msun/src/math.h Sat Jun 25 16:19:42 2005
@@ -397,8 +397,8 @@
long double atan2l(long double, long double);
long double atanhl(long double);
long double atanl(long double);
-long double cbrtl(long double);
#endif
+long double cbrtl(long double);
long double ceill(long double);
long double copysignl(long double, long double) __pure2;
#if 0
@@ -430,12 +430,14 @@
long long llrintl(long double);
#endif
long long llroundl(long double);
-#if 0
long double log10l(long double);
+#if 0
long double log1pl(long double);
long double log2l(long double);
long double logbl(long double);
+#endif
long double logl(long double);
+#if 0
long lrintl(long double);
#endif
long lroundl(long double);
@@ -460,7 +462,9 @@
#if 0
long double sinhl(long double);
long double sinl(long double);
+#endif
long double sqrtl(long double);
+#if 0
long double tanhl(long double);
long double tanl(long double);
long double tgammal(long double);
diff -rNu msun.orig/src/s_cbrtl.c msun/src/s_cbrtl.c
--- msun.orig/src/s_cbrtl.c Wed Dec 31 16:00:00 1969
+++ msun/src/s_cbrtl.c Sat Jun 25 16:05:17 2005
@@ -0,0 +1,84 @@
+/*-
+ * Copyright (c) 2005, 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 <math.h>
+#include <float.h>
+
+/*
+ * Compute the cube root of a long double argument by decomposing the
+ * the argument into its base 2 fraction and exponent.
+ *
+ * x**(1/3) = (f 2**n)**(1/3)
+ * = f**(1/3) * 2**(n/3) (mod(n,3) = 0)
+ * = f**(1/3) * 2**(1/3) * 2**[(n-1)/3] (mod(n,3) = 1)
+ * = f**(1/3) * 2**(2/3) * 2**[(n-2)/3] (mod(n,3) = 2)
+ *
+ * where f = [1/2,1).
+ *
+ * Use Newton's rule to compute f**(1/3)
+ * y_(k+1) = (x / y_k**2 + 2*y_k) / 3
+ * where k is the iteration count.
+ */
+
+long double cbrtl(long double x) {
+
+ int i, n, s = 1;
+ long double f, yk, oyk;
+
+ /* cbrtl(x) = NaN or +-Inf. */
+ if (isinf(x) || isnan(x))
+ return (x+x);
+
+ /* Save the sign. */
+ if (x < 0.L) {
+ s = -1;
+ x = -x;
+ }
+
+ f = frexpl(x, &n);
+
+ yk = oyk = f;
+ while(1) {
+ yk = (f / yk / yk + 2.0L * yk) / 3.0L;
+ if (fabsl(oyk - yk) < LDBL_EPSILON) break;
+ oyk = yk;
+ }
+
+ switch (n%3) {
+ case 0:
+ yk = ldexpl(yk, n / 3);
+ break;
+ case 1:
+ yk = ldexpl(yk, (n-1) / 3) * 1.259921049894873165L;
+ break;
+ case 2:
+ yk = ldexpl(yk, (n-2) / 3) * 1.587401051968199475L;
+ break;
+ }
+
+ return (s * yk);
+
+}
diff -rNu msun.orig/src/s_log10l.c msun/src/s_log10l.c
--- msun.orig/src/s_log10l.c Wed Dec 31 16:00:00 1969
+++ msun/src/s_log10l.c Sat Jun 25 16:05:34 2005
@@ -0,0 +1,67 @@
+/*-
+ * Copyright (c) 2005, 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.
+ */
+
+/*
+ * Compute the base 10 logrithm of the argument x by .
+ *
+ * log10(x) = log(x) / log(10).
+ *
+ * where log(x) is computed via logl(x).
+ */
+
+#include <math.h>
+
+#define LOG10L2 3.010299956639811952e-1L
+#define LOGL10 2.302585092994045684L
+
+static long double zero = 0.0L;
+
+long double log10l(long double x) {
+
+ int n;
+ long double f, val;
+
+ /* log10(x) = sNaN if x < 0. */
+ if (x < 0.0L)
+ return (x - x) / (x - x);
+
+ /* log10(+Inf) = +Inf */
+ if (isinf(x))
+ return x*x+x;
+
+ /* log10(+-0) = -Inf with signal */
+ if (x == 0.0L)
+ return (- 1.0L / zero);
+
+ /* log10(NaN) = NaN */
+ if (isnan(x))
+ return x;
+
+ val = logl(x) / LOGL10;
+
+ return val;
+
+}
diff -rNu msun.orig/src/s_logl.c msun/src/s_logl.c
--- msun.orig/src/s_logl.c Wed Dec 31 16:00:00 1969
+++ msun/src/s_logl.c Sat Jun 25 16:05:52 2005
@@ -0,0 +1,102 @@
+/*-
+ * Copyright (c) 2005, 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.
+ */
+
+/*
+ * Compute the natural logrithm of x by decomposing x into its base 2
+ * representation.
+ *
+ * log(x) = log(f * 2**n)
+ * = log(f) + n * log(2)
+ *
+ * where f = [1/2,1).
+ *
+ * Use a Taylor's series about f0 = 0.75 to compute log(f).
+ *
+ * log(f) = log(f0) + sum{(1/n) * (-1)**(n-1) * ((f - f0)/f0)**n}.
+ */
+
+#include <math.h>
+#include <float.h>
+
+#define LOGL2 6.931471805599453094e-1L
+#define LOGLF0 -2.876820724517809274e-1L
+
+static long double zero = 0.0L;
+
+long double logl(long double x) {
+
+ int i, n;
+ long double f, t, val, oval;
+
+ /* log(x) = sNaN if x < 0. */
+ if (x < 0.0L)
+ return (x - x) / (x - x);
+
+ /* log(+Inf) = +Inf */
+ if (isinf(x))
+ return x*x+x;
+
+ /* log(+-0) = -Inf with signal */
+ if (x == 0.0L)
+ return (- 1.0L / zero);
+
+ /* log(NaN) = NaN with signal */
+ if (isnan(x))
+ return (x+x);
+
+ /*
+ * Special case for values near log(1). Use the series expansion
+ * for log(1+x) = x - (1/2) * x**2 + (1/3) * x**3 + ...
+ */
+ if (0.95L < x && x < 1.05L) {
+
+ f = t = x - 1.0L;
+ oval = val = 0.L;
+
+ for (i = 1; ; i++) {
+ val += t / (long double) i;
+ t *= -f;
+ if (fabsl(oval - val) < LDBL_EPSILON) break;
+ oval = val;
+ }
+
+ return (val);
+ }
+
+ f = frexpl(x, &n);
+ f = t = (f - 0.75L) / 0.75L;
+
+ oval = val = LOGLF0;
+ for (i = 1; ; i++) {
+ val += t / (long double) i;
+ t *= -f;
+ if (fabsl(oval-val) < LDBL_EPSILON) break;
+ oval = val;
+ }
+
+ return (val + n * LOGL2);
+
+}
diff -rNu msun.orig/src/s_sqrtl.c msun/src/s_sqrtl.c
--- msun.orig/src/s_sqrtl.c Wed Dec 31 16:00:00 1969
+++ msun/src/s_sqrtl.c Sat Jun 25 16:06:06 2005
@@ -0,0 +1,77 @@
+/*-
+ * Copyright (c) 2005, 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.
+ */
+
+/*
+ * Compute the sqrt root of a long double argument by decomposing the
+ * the argument into its base 2 fraction and exponent.
+ *
+ * x**(1/2) = (f 2**n)**(1/2)
+ * = f**(1/2) * 2**(n/2) (n even)
+ * = f**(1/2) * 2**(1/2) * 2**[(n-1)/2] (n odd)
+ * where f = [1/2,1).
+ *
+ * Use Newton's rule to compute f
+ * y_(k+1) = (1/2) * (x / y_k + y_k)
+ * where k is the iteration count.
+ */
+
+#include <math.h>
+#include <float.h>
+
+long double sqrtl(long double x) {
+
+ int n;
+ long double f, yk, oyk;
+
+ /* sqrt(NaN) = NaN, sqrt(+Inf) = +Inf, or sqrt(-Inf) = sNaN */
+ if (isnan(x) || isinf(x))
+ return x*x+x;
+
+ /* sqrt(+-0) = 0 */
+ if (x == 0.0L)
+ return fabsl(x);
+
+ /* sqrt(x), x < 0 */
+ if (x < 0.0L)
+ return (x - x) / (x - x);
+
+ f = frexpl(x, &n);
+
+ yk = oyk = f;
+ while(1) {
+ yk = 0.5L * (f / yk + yk);
+ if (fabsl(oyk - yk) < LDBL_EPSILON) break;
+ oyk = yk;
+ }
+
+ if (n%2 == 0)
+ yk = ldexpl(yk, n / 2);
+ else
+ yk = 1.41421356237309505L * ldexpl(yk,(n - 1) / 2);
+
+ return yk;
+
+}
>Release-Note:
>Audit-Trail:
>Unformatted:
More information about the freebsd-standards
mailing list