From owner-freebsd-standards@FreeBSD.ORG Sun Jun 19 17:32:33 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 B621416A41F; Sun, 19 Jun 2005 17:32:33 +0000 (GMT) (envelope-from matteo@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [216.136.204.21]) by mx1.FreeBSD.org (Postfix) with ESMTP id 8E6F643D1F; Sun, 19 Jun 2005 17:32:33 +0000 (GMT) (envelope-from matteo@FreeBSD.org) Received: from freefall.freebsd.org (matteo@localhost [127.0.0.1]) by freefall.freebsd.org (8.13.3/8.13.3) with ESMTP id j5JHWXC7055988; Sun, 19 Jun 2005 17:32:33 GMT (envelope-from matteo@freefall.freebsd.org) Received: (from matteo@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j5JHWXAG055984; Sun, 19 Jun 2005 17:32:33 GMT (envelope-from matteo) Date: Sun, 19 Jun 2005 17:32:33 GMT From: Matteo Riondato Message-Id: <200506191732.j5JHWXAG055984@freefall.freebsd.org> To: dotz@irc.pl, matteo@FreeBSD.org, freebsd-standards@FreeBSD.org Cc: Subject: Re: standards/60597: FreeBSD's /usr/include lacks of cpio.h X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 19 Jun 2005 17:32:33 -0000 Synopsis: FreeBSD's /usr/include lacks of cpio.h State-Changed-From-To: patched->closed State-Changed-By: matteo State-Changed-When: Sun Jun 19 17:31:18 GMT 2005 State-Changed-Why: stefanf@ said this will never be MFCed to RELENG_4 http://www.freebsd.org/cgi/query-pr.cgi?pr=60597 From owner-freebsd-standards@FreeBSD.ORG Mon Jun 20 11:02:08 2005 Return-Path: X-Original-To: freebsd-standards@freebsd.org Delivered-To: freebsd-standards@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id F35DA16A422 for ; Mon, 20 Jun 2005 11:02:07 +0000 (GMT) (envelope-from owner-bugmaster@freebsd.org) Received: from freefall.freebsd.org (freefall.freebsd.org [216.136.204.21]) by mx1.FreeBSD.org (Postfix) with ESMTP id D861143D49 for ; Mon, 20 Jun 2005 11:02:07 +0000 (GMT) (envelope-from owner-bugmaster@freebsd.org) Received: from freefall.freebsd.org (peter@localhost [127.0.0.1]) by freefall.freebsd.org (8.13.3/8.13.3) with ESMTP id j5KB271T011578 for ; Mon, 20 Jun 2005 11:02:07 GMT (envelope-from owner-bugmaster@freebsd.org) Received: (from peter@localhost) by freefall.freebsd.org (8.13.3/8.13.1/Submit) id j5KB250H011572 for freebsd-standards@freebsd.org; Mon, 20 Jun 2005 11:02:05 GMT (envelope-from owner-bugmaster@freebsd.org) Date: Mon, 20 Jun 2005 11:02:05 GMT Message-Id: <200506201102.j5KB250H011572@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: peter set sender to owner-bugmaster@freebsd.org using -f From: FreeBSD bugmaster To: freebsd-standards@FreeBSD.org Cc: Subject: Current problem reports assigned to you X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 20 Jun 2005 11:02:08 -0000 Current FreeBSD problem reports Critical problems Serious problems S Submitted Tracker Resp. Description ------------------------------------------------------------------------------- o [2001/03/05] bin/25542 standards /bin/sh: null char in quoted string p [2002/02/25] standards/35307standards standard include files are not standard c o [2002/12/13] kern/46239 standards posix semaphore implementation errors o [2003/04/21] standards/51209standards [PATCH] add a64l()/l64a/l64a_r functions p [2003/06/05] standards/52972standards /bin/sh arithmetic not POSIX compliant o [2003/07/12] standards/54410standards one-true-awk not POSIX compliant (no exte o [2004/01/01] standards/60772standards _Bool and bool should be unsigned o [2004/11/03] standards/73500standards 'set +o' in /bin/sh does not include unse o [2005/03/03] standards/78357standards getaddrinfo() doesn't appear to support A 9 problems total. Non-critical problems S Submitted Tracker Resp. Description ------------------------------------------------------------------------------- o [2000/09/24] bin/21519 standards sys/dir.h should be deprecated some more o [2001/01/16] bin/24390 standards Replacing old dir-symlinks when using /bi s [2001/01/24] standards/24590standards timezone function not compatible witn Sin s [2001/06/18] kern/28260 standards UIO_MAXIOV needs to be made public p [2001/11/20] standards/32126standards getopt(3) not Unix-98 conformant o [2002/02/27] misc/35381 standards incorrect floating-point display of large s [2002/03/19] standards/36076standards Implementation of POSIX fuser command o [2002/06/14] standards/39256standards [v]snprintf aren't POSIX-conformant for s o [2002/07/09] kern/40378 standards stdlib.h gives needless warnings with -an p [2002/08/12] standards/41576standards POSIX compliance of ln(1) o [2002/10/23] standards/44425standards getcwd() succeeds even if current dir has o [2002/12/09] standards/46119standards Priority problems for SCHED_OTHER using p o [2003/07/24] standards/54809standards pcvt deficits o [2003/07/25] standards/54833standards more pcvt deficits o [2003/07/25] standards/54839standards pcvt deficits o [2003/07/31] standards/55112standards glob.h, glob_t's gl_pathc should be "size o [2003/09/05] standards/56476standards cd9660 unicode support simple hack o [2003/10/29] standards/58676standards grantpt(3) alters storage used by ptsname s [2004/02/14] standards/62858standards malloc(0) not C99 compliant p [2004/02/21] standards/63173standards Patch to add getopt_long_only(3) to libc o [2004/03/29] kern/64875 standards [patch] add a system call: fdatasync() o [2004/05/07] standards/66357standards make POSIX conformance problem ('sh -e' & o [2004/05/11] standards/66531standards _gettemp uses a far smaller set of filena o [2004/08/22] standards/70813standards [PATCH] ls not Posix compliant o [2004/08/26] docs/70985 standards [patch] sh(1): incomplete documentation o o [2004/09/22] standards/72006standards floating point formating in non-C locales o [2005/03/20] standards/79055standards Add an IFS regression test for shells o [2005/03/20] standards/79056standards regex(3) regression tests o [2005/03/21] standards/79067standards /bin/sh should be more intelligent about o [2005/05/20] standards/81287standards [PATCH]: fingerd(8) might send a line not 30 problems total. From owner-freebsd-standards@FreeBSD.ORG Sat Jun 25 17:14:57 2005 Return-Path: X-Original-To: freebsd-standards@freebsd.org Delivered-To: freebsd-standards@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id D4E0716A41C for ; Sat, 25 Jun 2005 17:14:57 +0000 (GMT) (envelope-from sgk@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 B556743D1D for ; Sat, 25 Jun 2005 17:14:57 +0000 (GMT) (envelope-from sgk@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 j5PHEvYq087058 for ; Sat, 25 Jun 2005 10:14:57 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.13.4/8.13.1/Submit) id j5PHEv24087057 for freebsd-standards@freebsd.org; Sat, 25 Jun 2005 10:14:57 -0700 (PDT) (envelope-from sgk) Date: Sat, 25 Jun 2005 10:14:57 -0700 From: Steve Kargl To: freebsd-standards@freebsd.org Message-ID: <20050625171457.GA86993@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline User-Agent: Mutt/1.4.2.1i Subject: Implementation of logl(3) X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 25 Jun 2005 17:14:58 -0000 I've noticed that several of C99's long double math function are not implemented. In an effort to start filling in some of the holes, I've implemented logl(3). Before I proceed in looking at other functions, I would appreciate comments because I'm not a number theorist by training. The code uses a Taylor series expansion to evaluate the log() of the fractional portion of x = f 2**n where f = [0.5,1). Although log(f) is very smooth over the indicated range, my attempts to fit a simple polynomial to the curve have failed. The failures are probably due to my fitting routines which are in double not long double (actually they use Fortran's double precision). Note, if the following is teemed okay, then cbrtl, sqrtl, and the other logrithm functions can be implemented in a similar manner. -- steve /* ln(x) = ln(f * 2**n) = ln(f) + n * ln(2) where f = [0.5,1). Use a Taylor series about f0 = 0.75 to compute ln(f). ln(f) = ln(f0) + sum (1/n) * (-1)**(n-1) * ((f - f0)/f0)**n */ #include #define LOGL2 6.931471805599453094e-1L #define LOGLF0 -2.876820724517809274e-1L static long double zero = 0.0L; #define NUM 32 long double logl(long double x) { int i, n, s = 1; long double f, t, val; /* 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(+-0) = -Inf with signal */ if (isnan(x)) return (x+x); f = frexpl(x, &n); f = t = (f - 0.75L) / 0.75L; val = LOGLF0; for (i = 1; i < NUM; i++) { val += s * t / (long double) i; t *= f; s = -s; } val += n * LOGL2; return val; } 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: