From owner-freebsd-standards@FreeBSD.ORG Sat Jun 25 23:40:12 2005 Return-Path: X-Original-To: freebsd-standards@hub.freebsd.org Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id C171216A41F for ; Sat, 25 Jun 2005 23:40:12 +0000 (GMT) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [216.136.204.21]) by mx1.FreeBSD.org (Postfix) with ESMTP id 657CF43D48 for ; Sat, 25 Jun 2005 23:40:12 +0000 (GMT) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (gnats@localhost [127.0.0.1]) by freefall.freebsd.org (8.13.3/8.13.3) with ESMTP id j5PNeCrb005985 for ; Sat, 25 Jun 2005 23:40:12 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j5PNeCqj005984; Sat, 25 Jun 2005 23:40:12 GMT (envelope-from gnats) Resent-Date: Sat, 25 Jun 2005 23:40:12 GMT Resent-Message-Id: <200506252340.j5PNeCqj005984@freefall.freebsd.org> Resent-From: FreeBSD-gnats-submit@FreeBSD.org (GNATS Filer) Resent-To: freebsd-standards@FreeBSD.org Resent-Reply-To: FreeBSD-gnats-submit@FreeBSD.org, "Steven G. Kargl" Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id 95B4216A41C for ; Sat, 25 Jun 2005 23:35:04 +0000 (GMT) (envelope-from kargl@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.FreeBSD.org (Postfix) with ESMTP id 6485243D49 for ; Sat, 25 Jun 2005 23:35:04 +0000 (GMT) (envelope-from kargl@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.13.4/8.13.4) with ESMTP id j5PNZ4Pr089560 for ; Sat, 25 Jun 2005 16:35:04 -0700 (PDT) (envelope-from kargl@troutmask.apl.washington.edu) Received: (from kargl@localhost) by troutmask.apl.washington.edu (8.13.4/8.13.1/Submit) id j5PNZ4Er089559; Sat, 25 Jun 2005 16:35:04 -0700 (PDT) (envelope-from kargl) Message-Id: <200506252335.j5PNZ4Er089559@troutmask.apl.washington.edu> Date: Sat, 25 Jun 2005 16:35:04 -0700 (PDT) From: "Steven G. Kargl" To: FreeBSD-gnats-submit@FreeBSD.org X-Send-Pr-Version: 3.113 Cc: Subject: standards/82654: C99 long double math functions are missing X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list Reply-To: "Steven G. Kargl" List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 25 Jun 2005 23:40:12 -0000 >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@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 +#include + +/* + * 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 + +#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 +#include + +#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 +#include + +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: