From owner-freebsd-standards@FreeBSD.ORG Mon Nov 15 05:06:18 2010 Return-Path: Delivered-To: standards@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id D8E87106566B for ; Mon, 15 Nov 2010 05:06:18 +0000 (UTC) (envelope-from yanegomi@gmail.com) Received: from mail-wy0-f182.google.com (mail-wy0-f182.google.com [74.125.82.182]) by mx1.freebsd.org (Postfix) with ESMTP id 63B248FC1C for ; Mon, 15 Nov 2010 05:06:18 +0000 (UTC) Received: by wyb36 with SMTP id 36so2157289wyb.13 for ; Sun, 14 Nov 2010 21:06:17 -0800 (PST) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:sender:received :in-reply-to:references:date:x-google-sender-auth:message-id:subject :from:to:cc:content-type:content-transfer-encoding; bh=DT4TGd1eCF49c24Hl8oxeidApMRwe0OIFhimZZlBEgY=; b=BKzpvqiYxEiPO9oR8cetpjDckSqGwhE79F5I/TWA+zBdBSuRx6SgSf+IjQfS8W7jnT Rr0vSHljtYIiCVB+7OHq7vKuR13vKGp2ijjyGKq5BLqADCgLQatDAOewjzAr9ciJgbgo vv9dY2xtiVNpUt41S+tCcfcF92VYbtZUdmJ30= DomainKey-Signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:sender:in-reply-to:references:date :x-google-sender-auth:message-id:subject:from:to:cc:content-type :content-transfer-encoding; b=pxu4tQE3Kq/6JjZN48YG4MSnwxP/l2YWloYdk39d+TbP9srgJg/MfVoXIhqXCRlc/9 EiiWeV3yh9Kkj9bbhHokClhm3tFA78uFUznRoOHpywbvhU+eCWnxx6KZY2EOgELcMV+l +1n7wksgrWU1rLoc/v8aqHM5nQ6rEppa9Bhm4= MIME-Version: 1.0 Received: by 10.216.46.200 with SMTP id r50mr6760479web.45.1289797577203; Sun, 14 Nov 2010 21:06:17 -0800 (PST) Sender: yanegomi@gmail.com Received: by 10.216.198.27 with HTTP; Sun, 14 Nov 2010 21:06:17 -0800 (PST) In-Reply-To: <20101109015704.L1243@besplex.bde.org> References: <20101109015704.L1243@besplex.bde.org> Date: Sun, 14 Nov 2010 21:06:17 -0800 X-Google-Sender-Auth: gcxNRZAKs13UDkQjiwuNg7jSHdU Message-ID: From: Garrett Cooper To: Bruce Evans Content-Type: text/plain; charset=ISO-8859-1 Content-Transfer-Encoding: quoted-printable Cc: standards@freebsd.org Subject: Re: Question about non-__BSD_VISIBLE guarded CLOCK_* constants 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, 15 Nov 2010 05:06:18 -0000 On Mon, Nov 8, 2010 at 7:59 AM, Bruce Evans wrote: > On Sun, 7 Nov 2010, Garrett Cooper wrote: > >> =A0 None of the following constants in time.h are guarded by >> __BSD_VISIBLE, __FreeBSD__, etc, even though other sections of the >> file are blocked off that way (and the comments suggest that they're >> FreeBSD-specific). I was wondering why that's the case... > > It is because the other sections were written by standards-aware persons. > >> #define CLOCK_UPTIME =A0 =A05 =A0 =A0 =A0 =A0 =A0 =A0 =A0 /* FreeBSD-spe= cific. */ >> #define CLOCK_UPTIME_PRECISE =A0 =A07 =A0 =A0 =A0 /* FreeBSD-specific. *= / >> #define CLOCK_UPTIME_FAST =A0 =A0 =A0 8 =A0 =A0 =A0 /* FreeBSD-specific.= */ >> #define CLOCK_REALTIME_PRECISE =A09 =A0 =A0 =A0 /* FreeBSD-specific. */ >> #define CLOCK_REALTIME_FAST =A0 =A0 10 =A0 =A0 =A0/* FreeBSD-specific. *= / >> #define CLOCK_MONOTONIC_PRECISE 11 =A0 =A0 =A0/* FreeBSD-specific. */ >> #define CLOCK_MONOTONIC_FAST =A0 =A012 =A0 =A0 =A0/* FreeBSD-specific. *= / >> #define CLOCK_SECOND =A0 =A013 =A0 =A0 =A0 =A0 =A0 =A0 =A0/* FreeBSD-spe= cific. */ > > Defining these without using __BSD_VISBLE is just a style bug, since the > CLOCK_ prefix is reserved in , but FreeBSD normally ifdefs > extensions that use reserved prefixes anyway (e.g., FreeBSD-specific > signal numbers at least used to be ifdefed to a fault in . > > The definitions themselves have lots of style bugs: > - space instead of tab after #define (all tabs were corrupted in the mail= , > =A0except these are corrupt in the file) > - macro values don't line up, although they are indented with tabs > - comment duplicated ad nauseum > - there is now an id 14. =A0It is missing the comment, and seems to be > =A0undocumented. > >> =A0 Also, they're blocked off by #if !defined(CLOCK_REALTIME) && >> __POSIX_VISIBLE >=3D 200112 , which doesn't seem to make sense, given >> that it's an "advanced realtime" feature, according to POSIX 2008. > > Re-quoting everything to get more context. > > % /* These macros are also in sys/time.h. */ > > is considerably more broken than here. =A0It defines all of > these macros unconditionally, and also includes in the !_KERNEL > case. =A0Thus including the POSIX header gives massive names= pace > pollution. > > % #if !defined(CLOCK_REALTIME) && __POSIX_VISIBLE >=3D 200112 > % #define CLOCK_REALTIME =A0 =A0 =A0 =A00 > % #ifdef __BSD_VISIBLE > % #define CLOCK_VIRTUAL 1 > % #define CLOCK_PROF =A0 =A02 > % #endif > % #define CLOCK_MONOTONIC =A0 =A0 =A0 4 > /* These macros are also in sys/time.h. */ > #define CLOCK_MONOTONIC 4 > #define CLOCK_UPTIME =A0 =A05 =A0 =A0 =A0 =A0 =A0 =A0 =A0 /* FreeBSD-spec= ific. */ > #define ... > #endif /* !defined(CLOCK_REALTIME) && __POSIX_VISIBLE >=3D 200112 */ > > The CLOCK_REALTIME part of this does nothing except partially break > automatic checking that the defines here are the same as the ones in > : > - if is not included before here, then CLOCK_REALTIME should= n't > =A0be defined yet, so it has no effect in this ifdef > - if is =A0 =A0 included before here, then CLOCK_REALTIME > =A0is defined now, so its effect in this ifdef is to prevent repeating if= defs > =A0that should be the same. =A0Repeated ifdefs are a very small part of h= eader > =A0bloat, so it is hardly worth avoiding them, and not avoiding them give= s the > =A0feature of checking that they really are repeated. =A0I think the repe= titions > =A0are all indentical, including their style bugs and whitespace that wou= ld > =A0not be checked. > - if is not included before here, but is included later, the= n > =A0CLOCK_REALTIME etc. gets declared here, and since has no > =A0similar ifdefs, it gets defined again later. =A0Thus the automatic che= cking > =A0is only partially broken. > > As you noticed, the __POSIX_VISIBLE part of this has the wrong version > number > at best. =A0CLOCK_REALTIME dates from POSIX.4 in ~1994 (it is in my POSIX= .4 > book published in 1995, and in the 1996 POSIX standard). =A0FreeBSD doesn= 't > try to track old versions of POSIX very carefully, except original POSIX, > but it normally makes features that appeared in in-between POSIXes visibl= e > by default by putting them under __BSD_VISIBLE. > > So the correct ifdefs here aresomething like __BSD_VISIBLE || > POSIX_VISIBLE >=3D 199309 (check the latter) for the whole block, and > __BSD_VISIBLE for the sub-block consisting of FreeBSD extensions. > Since the organization is poor (historical, giving random alphabetical an= d > numericl order), there are actually several sub-blocks: > - CLOCK_MONOTONIC seems to be new in POSIX.1-2001. =A0Thus the current if= def > =A0is partially correct for it. =A0The rest of the file has careful ifdef= s for > =A01993 POSIX, so perhaps it is worth being equally careful for > CLOCK_REALTIME. > - the BSD extensions CLOCK_VIRTUAL and CLOCK_PROF are already in a > =A0__BSD_VISIBLE block. =A0Then there seems to be space for expansion (a > =A0whole 1 id =3D 3). =A0Then there is the whole 1 new id =3D 4 for 2001 = POSIX. > =A0Then there is a FreeBSD-specific block (ids 5-14). =A0Too much bloat f= or me. Looks like I got the isolation working properly: # gcc -D_BSD_SOURCE -o ~gcooper/test_clock_uptime ~gcooper/test_clock_uptim= e.c # ~gcooper/test_clock_uptime 100031 661918304 # gcc -D_POSIX_C_SOURCE=3D200809 -D_POSIX_SOURCE -o ~gcooper/test_clock_uptime ~gcooper/test_clock_uptime.c /home/gcooper/test_clock_uptime.c: In function 'main': /home/gcooper/test_clock_uptime.c:10: error: 'CLOCK_UPTIME' undeclared (first use in this function) /home/gcooper/test_clock_uptime.c:10: error: (Each undeclared identifier is reported only once /home/gcooper/test_clock_uptime.c:10: error: for each function it appears i= n.) Removing time.h looks like a mini-project in and of itself, as I've had to go fix a handful of files in libc that were looking for clock or timer definitions that were defined in time.h according to POSIX, found a few bugs elsewhere when running buildworld, etc. I'll be submitting patches for all of these items eventually, but it looks like the breakage is potentially quite large. Thanks, -Garrett From owner-freebsd-standards@FreeBSD.ORG Mon Nov 15 11:07:05 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 4CAF9106566B for ; Mon, 15 Nov 2010 11:07:05 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 394738FC22 for ; Mon, 15 Nov 2010 11:07:05 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.4/8.14.4) with ESMTP id oAFB75Cf086430 for ; Mon, 15 Nov 2010 11:07:05 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.4/8.14.4/Submit) id oAFB74dJ086428 for freebsd-standards@FreeBSD.org; Mon, 15 Nov 2010 11:07:04 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 15 Nov 2010 11:07:04 GMT Message-Id: <201011151107.oAFB74dJ086428@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-standards@FreeBSD.org Cc: Subject: Current problem reports assigned to freebsd-standards@FreeBSD.org 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, 15 Nov 2010 11:07:05 -0000 Note: to view an individual PR, use: http://www.freebsd.org/cgi/query-pr.cgi?pr=(number). The following is a listing of current problems submitted by FreeBSD users. These represent problem reports covering all versions including experimental development code and obsolete releases. S Tracker Resp. Description -------------------------------------------------------------------------------- o stand/151316 standards lib/libc/string/strerror.c r1.9 breaks POSIX o stand/150093 standards C++ std::locale support is broken a stand/149980 standards [libc] [patch] negative value integer to nanosleep(2) o stand/147210 standards xmmintrin.h and cstdlib conflicts with each other with p stand/145517 standards POSIX getline() missing o stand/144231 standards bind/connect/sendto too strict about sockaddr length o stand/143358 standards [libm] nearbyint(3) raises spurious inexact exception o stand/142803 standards j0 Bessel function inaccurate near zeros of the functi s stand/141705 standards [libc] [request] libc lacks cexp (and friends) o stand/130067 standards Wrong numeric limits in system headers? o stand/124860 standards flockfile(3) doesn't work when the memory has been exh o stand/123688 standards POSIX standard changes in unistd.h and grp.h o stand/121921 standards [patch] Add leap second support to at(1), atrun(8) o stand/116826 standards [patch] sh support for POSIX character classes o stand/116477 standards rm(1): rm behaves unexpectedly when using -r and relat o bin/116413 standards incorrect getconf(1) handling of unsigned constants gi o stand/116081 standards make does not work with the directive sinclude o stand/114633 standards /etc/rc.subr: line 511: omits a quotation mark: "force p stand/107561 standards [libc] [patch] [request] Missing SUS function tcgetsid o stand/104743 standards [headers] [patch] Wrong values for _POSIX_ minimal lim o stand/100017 standards [Patch] Add fuser(1) functionality to fstat(1) o stand/96236 standards [patch] [posix] sed(1) incorrectly describes a functio o stand/96016 standards [headers] clock_getres et al should be in o stand/94729 standards [libc] fcntl() throws undocumented ENOTTY o kern/93705 standards [headers] [patch] ENODATA and EGREGIOUS (for glibc com o stand/92362 standards [headers] [patch] Missing SIGPOLL in kernel headers a stand/86484 standards [patch] mkfifo(1) uses wrong permissions o stand/83845 standards [libm] [patch] add log2() and log2f() support for libm o stand/82654 standards C99 long double math functions are missing o stand/81287 standards [patch] fingerd(8) might send a line not ending in CRL a stand/80293 standards sysconf() does not support well-defined unistd values o stand/79056 standards [feature request] [atch] regex(3) regression tests o stand/70813 standards [patch] ls(1) not Posix compliant o stand/66357 standards make POSIX conformance problem ('sh -e' & '+' command- s kern/64875 standards [libc] [patch] [request] add a system call: fdatasync( s stand/62858 standards malloc(0) not C99 compliant o stand/56476 standards [patch] cd9660 unicode support simple hack o stand/54410 standards one-true-awk not POSIX compliant (no extended REs) o stand/46119 standards Priority problems for SCHED_OTHER using pthreads o stand/44365 standards [headers] [patch] [request] introduce ulong and unchar a stand/41576 standards ln(1): replacing old dir-symlinks o stand/39256 standards snprintf/vsnprintf aren't POSIX-conformant for strings o kern/27835 standards [libc] execve() doesn't conform to execve(2) spec in s a docs/26003 standards getgroups(2) lists NGROUPS_MAX but not syslimits.h s stand/24590 standards timezone function not compatible witn Single Unix Spec o stand/21519 standards sys/dir.h should be deprecated some more s bin/14925 standards getsubopt isn't poisonous enough 47 problems total. From owner-freebsd-standards@FreeBSD.ORG Fri Nov 19 23:10:09 2010 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 6D79C106566C for ; Fri, 19 Nov 2010 23:10:09 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 2A52D8FC16 for ; Fri, 19 Nov 2010 23:10:09 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.4/8.14.4) with ESMTP id oAJNA9Us035519 for ; Fri, 19 Nov 2010 23:10:09 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.4/8.14.4/Submit) id oAJNA9t5035518; Fri, 19 Nov 2010 23:10:09 GMT (envelope-from gnats) Resent-Date: Fri, 19 Nov 2010 23:10:09 GMT Resent-Message-Id: <201011192310.oAJNA9t5035518@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 [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 8594E106564A for ; Fri, 19 Nov 2010 23:07:19 +0000 (UTC) (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 626058FC14 for ; Fri, 19 Nov 2010 23:07:19 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.4/8.14.4) with ESMTP id oAJN7JUY030126 for ; Fri, 19 Nov 2010 15:07:19 -0800 (PST) (envelope-from kargl@troutmask.apl.washington.edu) Received: (from kargl@localhost) by troutmask.apl.washington.edu (8.14.4/8.14.4/Submit) id oAJN7JPK030125; Fri, 19 Nov 2010 15:07:19 -0800 (PST) (envelope-from kargl) Message-Id: <201011192307.oAJN7JPK030125@troutmask.apl.washington.edu> Date: Fri, 19 Nov 2010 15:07:19 -0800 (PST) From: "Steven G. Kargl" To: FreeBSD-gnats-submit@FreeBSD.org X-Send-Pr-Version: 3.113 Cc: Subject: standards/152415: [libm] implementation of expl() 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: Fri, 19 Nov 2010 23:10:09 -0000 >Number: 152415 >Category: standards >Synopsis: [libm] implementation of expl() >Confidential: no >Severity: non-critical >Priority: low >Responsible: freebsd-standards >State: open >Quarter: >Keywords: >Date-Required: >Class: change-request >Submitter-Id: current-users >Arrival-Date: Fri Nov 19 23:10:08 UTC 2010 >Closed-Date: >Last-Modified: >Originator: Steven G. Kargl >Release: FreeBSD 9.0-CURRENT amd64 >Organization: apl/uw >Environment: System: FreeBSD troutmask.apl.washington.edu 9.0-CURRENT FreeBSD 9.0-CURRENT #6 r213691M: Tue Oct 26 11:01:30 PDT 2010 kargl@troutmask.apl.washington.edu:/usr/obj/usr/src/sys/SPEW amd64 >Description: FreeBSD currently lacks a long double version of the exponential function, also known as expl(). An algorithm for both float and double versions of this function, expf() and exp(), has been given by Tang in PTP Tang, "Table-driven implementation of the exponential function in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, 144-157 (1989). I have implemented Tang's algorithm for expf(), exp(), and expl(). My testing on amd64 has found Tang Libm Range Max ULP sec Max ULP sec # of calls expf: [-18.0:88.72] 0.525668 1.11 0.901775 1.33 23124111 exp: [-37.0:709.7] 0.526008 1.47 0.895724 1.60 23000001 expl: [-44.0:11356] 0.506063 2.81 23000001 where "Range" is the range over which the function will return a finite normal result. The 'sec' column is the self-seconds reported by gprof for the "# of calls" to the function. Thus, for over 23M function calls expf() and exp() based on Tang are more accurate and actually faster than the corresponding functions in libm. In the included code, the functions are named t_expf(), t_exp(), and t_expl(), which permits inclusion into libm for ease of testing. If these implementations are to supercede the versions available in libm, then rename the functions accordingly. >How-To-Repeat: >Fix: /*- * Copyright (c) 2009-2010 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. */ /* * Single precision implemenation of exp(x) as outlined in * * PTP Tang, "Table-driven implementation of the exponential function * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, * 144-157 (1989). */ #include #include #include "math_private.h" /* * Tang states that THRESHOLD_1 is 341*log(2), which is about 236.36. But, * the constant given in the Appendix is 0x435c6bba, which is about 220.42. */ #define THRESHOLD_1 0x1.d8b9f4p+7f /* 341 * log(2) */ #define THRESHOLD_2 0x1.000000p-25f /* 2**(-25) */ #define INV_L 0x1.715476p+5f /* 32 / log(2) */ #define L1 0x1.62e400p-6f /* L1+L2=log(2)/32 */ #define L2 0x1.7f7d1cp-25f /* L1+L2=log(2)/32 */ /* Remes polynomial coefficients. */ #define A1 0x1.000088p-1f #define A2 0x1.5555d8p-3f /* Table of values for 2**(j/32) with j = 0, ..., 31. */ static const struct { float head; float tail; } s[32] = { 0x1.00000p+0f, 0x0.000000p+00f, 0x1.059b0p+0f, 0x1.a62b0ap-21f, 0x1.0b558p+0f, 0x1.b3e624p-22f, 0x1.11300p+0f, 0x1.d0125cp-20f, 0x1.172b8p+0f, 0x1.e3ea8cp-23f, 0x1.1d480p+0f, 0x1.cc5a2ep-18f, 0x1.23878p+0f, 0x1.373ab2p-19f, 0x1.29e98p+0f, 0x1.7d47f8p-18f, 0x1.306f8p+0f, 0x1.828c6ep-18f, 0x1.371a0p+0f, 0x1.cdceaap-18f, 0x1.3dea0p+0f, 0x1.93048ep-18f, 0x1.44e08p+0f, 0x1.818624p-22f, 0x1.4bfd8p+0f, 0x1.6a9b16p-19f, 0x1.53428p+0f, 0x1.ab4ea8p-19f, 0x1.5ab00p+0f, 0x1.f75216p-18f, 0x1.62478p+0f, 0x1.ac0e96p-18f, 0x1.6a098p+0f, 0x1.999fcep-18f, 0x1.71f70p+0f, 0x1.7a3b18p-18f, 0x1.7a110p+0f, 0x1.1cfac0p-18f, 0x1.82588p+0f, 0x1.994ccep-20f, 0x1.8ace0p+0f, 0x1.508aa8p-18f, 0x1.93730p+0f, 0x1.ec3372p-18f, 0x1.9c490p+0f, 0x1.82a3f0p-20f, 0x1.a5500p+0f, 0x1.d91f12p-19f, 0x1.ae898p+0f, 0x1.e656b4p-18f, 0x1.b7f70p+0f, 0x1.bcbed8p-18f, 0x1.c1998p+0f, 0x1.eec2aap-19f, 0x1.cb720p+0f, 0x1.b9df20p-21f, 0x1.d5818p+0f, 0x1.b9f74ap-21f, 0x1.dfc90p+0f, 0x1.ccdee6p-18f, 0x1.ea4a8p+0f, 0x1.e8a924p-18f, 0x1.f5070p+0f, 0x1.96db92p-18f }; static const float tiny = 0x1.p-100f; /* * From Tang's paper: INTRND rounds a floating-point number to the nearest * integer in the manner prescribed by the IEEE standard. It is crucial that * the default round-to-nearest mode, not any other rounding mode, is in * effect here. */ #define INTRND(x) roundf((x)) float t_expf(float x) { int32_t xsb; u_int32_t hx; int n, n1, n2; float ax, r1, r2, r, p, q, t; GET_FLOAT_WORD(hx,x); /* Determine the sign bit of x. */ xsb = (hx >> 31) & 1; /* Set high word to |x|. */ hx &= 0x7fffffff; /* Test for NaN and +-Inf. */ if (hx >= 0x7f800000) { /* exp(NaN) = NaN. */ if (hx > 0x7f800000) return (x + x); /* exp(+inf) = inf (no exception!) or exp(-inf) = 0. */ return (xsb ? 0 : x); } SET_FLOAT_WORD(ax, hx); if (ax > THRESHOLD_1) return (x > 0 ? (ax / (ax - ax)) : (0 + tiny * tiny)); if (ax < THRESHOLD_2) return (1 + x); n = INTRND(x * INV_L); n2 = n % 32; if (n2 < 0) n2 += 32; /* Tang's j. */ n1 = n - n2; if (n >= 512 || n <= -512) r1 = (x - n1 * L1) - n2 * L1; else r1 = x - n * L1; r2 = - n * L2; n1 /= 32; /* Tang's m. */ r = r1 + r2; q = r * r * (A1 + r * A2); p = r1 + (r2 + q); t = s[n2].head + s[n2].tail; t = s[n2].head + (s[n2].tail + t * p); return (ldexpf(t, n1)); } /*- * Copyright (c) 2009-2010 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. */ /* * Double precision implemenation of exp(x) as outlined in * * PTP Tang, "Table-driven implementation of the exponential function * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, * 144-157 (1989). * */ #include #include "math.h" #include "math_private.h" #define THRESHOLD_1 0x1.c4474e1726454p+10 /* 2610 * log(2) */ #define THRESHOLD_2 0x1.0000000000000p-54 /* 2**(-54) */ #define INV_L 0x1.71547652b82fep+5 /* 32 / log(2) */ #define L1 0x1.62e42fef00000p-6 /* L1+L2=log(2)/32 */ #define L2 0x1.473de6af278ecp-39 /* L1+L2=log(2)/32 */ /* Remes polynomial coefficients. */ #define A1 0x1.0000000000000p-1 #define A2 0x1.5555555548f7cp-3 #define A3 0x1.5555555545d4ep-5 #define A4 0x1.11115b7aa905ep-7 #define A5 0x1.6c1728d739764p-10 static const struct { double head; double tail; } s[32] = { 0x1p+0, 0x0p+0, 0x1.059b0d315854p+0, 0x1.a1d73e2a475b4p-47, 0x1.0b5586cf9890p+0, 0x1.ec5317256e308p-49, 0x1.11301d0125b4p+0, 0x1.0a4ebbf1aed93p-48, 0x1.172b83c7d514p+0, 0x1.d6e6fbe462876p-47, 0x1.1d4873168b98p+0, 0x1.53c02dc0144c8p-47, 0x1.2387a6e75620p+0, 0x1.c3360fd6d8e0bp-47, 0x1.29e9df51fdecp+0, 0x1.09612e8afad12p-47, 0x1.306fe0a31b70p+0, 0x1.52de8d5a46306p-48, 0x1.371a7373aa9cp+0, 0x1.54e28aa05e8a9p-49, 0x1.3dea64c12340p+0, 0x1.11ada0911f09fp-47, 0x1.44e086061890p+0, 0x1.68189b7a04ef8p-47, 0x1.4bfdad5362a0p+0, 0x1.38ea1cbd7f621p-47, 0x1.5342b569d4f8p+0, 0x1.df0a83c49d86ap-52, 0x1.5ab07dd48540p+0, 0x1.4ac64980a8c8fp-47, 0x1.6247eb03a558p+0, 0x1.2c7c3e81bf4b7p-50, 0x1.6a09e667f3bcp+0, 0x1.921165f626cddp-49, 0x1.71f75e8ec5f4p+0, 0x1.9ee91b8797785p-47, 0x1.7a11473eb018p+0, 0x1.b5f54408fdb37p-50, 0x1.82589994cce0p+0, 0x1.28acf88afab35p-48, 0x1.8ace5422aa0cp+0, 0x1.b5ba7c55a192dp-48, 0x1.93737b0cdc5cp+0, 0x1.27a280e1f92a0p-47, 0x1.9c49182a3f08p+0, 0x1.01c7c46b071f3p-48, 0x1.a5503b23e254p+0, 0x1.c8b424491caf8p-48, 0x1.ae89f995ad38p+0, 0x1.6af439a68bb99p-47, 0x1.b7f76f2fb5e4p+0, 0x1.baa9ec206ad4fp-50, 0x1.c199bdd85528p+0, 0x1.c2220cb12a092p-48, 0x1.cb720dcef904p+0, 0x1.48a81e5e8f4a5p-47, 0x1.d5818dcfba48p+0, 0x1.c976816bad9b8p-50, 0x1.dfc97337b9b4p+0, 0x1.eb968cac39ed3p-48, 0x1.ea4afa2a490cp+0, 0x1.9858f73a18f5ep-48, 0x1.f50765b6e454p+0, 0x1.9d3e12dd8a18bp-54 }; static const double tiny = 0x1.p-1000; /* * From Tang's paper: INTRND rounds a floating-point number to the nearest * integer in the manner prescribed by the IEEE standard. It is crucial that * the default round-to-nearest mode, not any other rounding mode, is in * effect here. */ #define INTRND(x) round((x)) double t_exp(double x) { u_int32_t hx, lx, xsb; int n, n1, n2; double ax, r1, r2, r, p, q, t; EXTRACT_WORDS(hx, lx, x); /* Detemine the sign bit of x. */ xsb = (hx >> 31) & 1; /* Set high word of |x| */ hx &= 0x7fffffff; /* Test for NaN and +-Inf. */ if (hx >= 0x7ff00000) { /* exp(NaN) = NaN. */ if (((hx & 0xfffff) | lx) != 0) return (x + x); /* exp(+inf) = inf (no exception!) or exp(-inf) = 0. */ return (xsb ? 0 : x); } ax = x; SET_HIGH_WORD(ax, hx); if (ax > THRESHOLD_1) return (x > 0 ?(ax / (ax - ax)) : (0 + tiny * tiny)); if (ax < THRESHOLD_2) return (1 + x); n = INTRND(x * INV_L); n2 = n % 32; if (n2 < 0) n2 += 32; /* Tang's j. */ n1 = n - n2; if (n >= 512 || n <= -512) r1 = (x - n1 * L1) - n2 * L1; else r1 = x - n * L1; r2 = - n * L2; n1 /= 32; /* Tang's m. */ r = r1 + r2; q = r * r * (A1 + r * (A2 + r * (A3 + r * (A4 + r * A5)))); p = r1 + (r2 + q); t = s[n2].head + s[n2].tail; t = s[n2].head + (s[n2].tail + t * p); return (ldexp(t, n1)); } /*- * Copyright (c) 2009-2010 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. */ /* * Implemenation of expl(x) in Intel 80-bit long double format. * The significand is 64 bits. This is based on * * PTP Tang, "Table-driven implementation of the exponential function * in IEEE floating-point arithmetic," ACM Trans. Math. Soft., 15, * 144-157 (1989). * * where the 32 table entries has been expanded to NUM (see below). * */ #include #include "math.h" #include "fpmath.h" #define BIAS (LDBL_MAX_EXP - 1) /* * The range of overflow to +Inf and of underflow to +0 is given by * 41021 * log(2). 41021 is obtained from -(-16382-64+1)+24576 where * -16382 is the minimum exponent, 64 is the precision and 24576 is * a bias adjustment given in IEEE 754 as 3*2^(e-2) with e being the * width of the exponent in bits. */ static const long double THRESHOLD_1 = 2.84335904937495165381e+04L, /* 41021 * log(2) */ THRESHOLD_2 = 2.71050543121376108502e-20L, /* 2**(-65) */ L1 = 5.4152123481245872938e-03L, /* L1+L2=log(2)/128 */ L2 = -1.4563974013471118145e-17L, /* L1+L2=log(2)/128 */ INV_L = 1.84664965233787316146e+02L; /* 128 / log(2) */ /* * The polynomial coefficients are the initial rational approximation * of the minimax polynomial coefficients determined by a Remes algorithm * coded in a much higher working precision. For details, see * * N. Brisebarre, et al, "Computing Machine-Efficient Polynomial * Approximations," ACM Tran. Math. Soft., 32, 236-256 (2006). */ static const long double P1 = 5.000000000000000000e-01L, P2 = 1.666666666666666666e-01L, P3 = 4.166666666666666677e-02L, P4 = 8.333333333333284767e-03L, P5 = 1.388888888835006555e-03L, P6 = 1.984127249252972877e-04L, P7 = 2.480718929383586384e-05L; /* * 2^(i/NUM) for i in [0,NUM] is represented by two values where the * first 47-bits of the significand is stored in head and the next 53 * bits are in tail. */ #define NUM 128 static const struct { double head; double tail; } s[NUM] = { 0x1p+0, 0x0p+0, 0x1.0163da9fb330p+0, 0x1.ab6c25335719bp-47, 0x1.02c9a3e77804p+0, 0x1.07737be56527cp-47, 0x1.04315e86e7f8p+0, 0x1.2f5ce3e688369p-50, 0x1.059b0d315854p+0, 0x1.a1d73e2a475b4p-47, 0x1.0706b29ddf6cp+0, 0x1.dc6dc403a9d88p-48, 0x1.0874518759bcp+0, 0x1.01186be4bb285p-49, 0x1.09e3ecac6f38p+0, 0x1.a290f03062c27p-51, 0x1.0b5586cf9890p+0, 0x1.ec5317256e308p-49, 0x1.0cc922b7247cp+0, 0x1.ba03db82dc49fp-47, 0x1.0e3ec32d3d18p+0, 0x1.10103a1727c58p-47, 0x1.0fb66affed30p+0, 0x1.af232091dd8a1p-48, 0x1.11301d0125b4p+0, 0x1.0a4ebbf1aed93p-48, 0x1.12abdc06c31cp+0, 0x1.7f72575a649adp-49, 0x1.1429aaea92dcp+0, 0x1.fb34101943b26p-48, 0x1.15a98c8a58e4p+0, 0x1.12480d573dd56p-48, 0x1.172b83c7d514p+0, 0x1.d6e6fbe462876p-47, 0x1.18af9388c8dcp+0, 0x1.4dddfb85cd1e1p-47, 0x1.1a35beb6fcb4p+0, 0x1.a9e5b4c7b4969p-47, 0x1.1bbe084045ccp+0, 0x1.39ab1e72b4428p-48, 0x1.1d4873168b98p+0, 0x1.53c02dc0144c8p-47, 0x1.1ed5022fcd90p+0, 0x1.cb8819ff61122p-48, 0x1.2063b88628ccp+0, 0x1.63b8eeb029509p-48, 0x1.21f49917ddc8p+0, 0x1.62552fd29294cp-48, 0x1.2387a6e75620p+0, 0x1.c3360fd6d8e0bp-47, 0x1.251ce4fb2a60p+0, 0x1.f9ac155bef4f5p-47, 0x1.26b4565e27ccp+0, 0x1.d257a673281d4p-48, 0x1.284dfe1f5638p+0, 0x1.2d9e2b9e07941p-53, 0x1.29e9df51fdecp+0, 0x1.09612e8afad12p-47, 0x1.2b87fd0dad98p+0, 0x1.ffbbd48ca71f9p-49, 0x1.2d285a6e4030p+0, 0x1.680123aa6da0fp-49, 0x1.2ecafa93e2f4p+0, 0x1.611ca0f45d524p-48, 0x1.306fe0a31b70p+0, 0x1.52de8d5a46306p-48, 0x1.32170fc4cd80p+0, 0x1.89a9ce78e1804p-47, 0x1.33c08b26416cp+0, 0x1.fa64e43086cb3p-47, 0x1.356c55f929fcp+0, 0x1.864a311a3b1bap-47, 0x1.371a7373aa9cp+0, 0x1.54e28aa05e8a9p-49, 0x1.38cae6d05d84p+0, 0x1.2c2d4e586cdf7p-47, 0x1.3a7db34e59fcp+0, 0x1.b750de494cf05p-47, 0x1.3c32dc313a8cp+0, 0x1.242000f9145acp-47, 0x1.3dea64c12340p+0, 0x1.11ada0911f09fp-47, 0x1.3fa4504ac800p+0, 0x1.ba0bf701aa418p-48, 0x1.4160a21f72e0p+0, 0x1.4fc2192dc79eep-47, 0x1.431f5d950a88p+0, 0x1.6dc704439410dp-48, 0x1.44e086061890p+0, 0x1.68189b7a04ef8p-47, 0x1.46a41ed1d004p+0, 0x1.772512f45922ap-48, 0x1.486a2b5c13ccp+0, 0x1.013c1a3b69063p-48, 0x1.4a32af0d7d3cp+0, 0x1.e672d8bcf46f9p-48, 0x1.4bfdad5362a0p+0, 0x1.38ea1cbd7f621p-47, 0x1.4dcb299fddd0p+0, 0x1.ac766dde353c2p-49, 0x1.4f9b2769d2c8p+0, 0x1.35699ec5b4d50p-47, 0x1.516daa2cf664p+0, 0x1.c112f52c84d82p-52, 0x1.5342b569d4f8p+0, 0x1.df0a83c49d86ap-52, 0x1.551a4ca5d920p+0, 0x1.d8a5d8c40486ap-49, 0x1.56f4736b527cp+0, 0x1.a66ecb004764fp-48, 0x1.58d12d497c7cp+0, 0x1.e9295e15b9a1ep-47, 0x1.5ab07dd48540p+0, 0x1.4ac64980a8c8fp-47, 0x1.5c9268a59468p+0, 0x1.b80e258dc0b4cp-47, 0x1.5e76f15ad214p+0, 0x1.0dd37c9840733p-49, 0x1.605e1b976dc0p+0, 0x1.160edeb25490ep-49, 0x1.6247eb03a558p+0, 0x1.2c7c3e81bf4b7p-50, 0x1.6434634ccc30p+0, 0x1.fc76f8714c4eep-48, 0x1.662388255220p+0, 0x1.24893ecf14dc8p-47, 0x1.68155d44ca94p+0, 0x1.9840e2b913dd0p-47, 0x1.6a09e667f3bcp+0, 0x1.921165f626cddp-49, 0x1.6c012750bda8p+0, 0x1.f76bb54cc007ap-47, 0x1.6dfb23c651a0p+0, 0x1.779107165f0dep-47, 0x1.6ff7df951948p+0, 0x1.e7c3f0da79f11p-51, 0x1.71f75e8ec5f4p+0, 0x1.9ee91b8797785p-47, 0x1.73f9a48a5814p+0, 0x1.9deae4d273456p-47, 0x1.75feb564267cp+0, 0x1.17edd35467491p-49, 0x1.780694fde5d0p+0, 0x1.fb0cd7014042cp-47, 0x1.7a11473eb018p+0, 0x1.b5f54408fdb37p-50, 0x1.7c1ed0130c10p+0, 0x1.93e2499a22c9cp-47, 0x1.7e2f336cf4e4p+0, 0x1.1082e815d0abdp-47, 0x1.80427543e1a0p+0, 0x1.1b60de67649a3p-48, 0x1.82589994cce0p+0, 0x1.28acf88afab35p-48, 0x1.8471a4623c78p+0, 0x1.667297b5cbe32p-47, 0x1.868d99b4492cp+0, 0x1.640720ec85613p-47, 0x1.88ac7d98a668p+0, 0x1.966530bcdf2d5p-48, 0x1.8ace5422aa0cp+0, 0x1.b5ba7c55a192dp-48, 0x1.8cf3216b5448p+0, 0x1.7de55439a2c39p-49, 0x1.8f1ae9915770p+0, 0x1.b15cc13a2e397p-47, 0x1.9145b0b91ffcp+0, 0x1.622986d1a7daep-50, 0x1.93737b0cdc5cp+0, 0x1.27a280e1f92a0p-47, 0x1.95a44cbc8520p+0, 0x1.dd36906d2b420p-49, 0x1.97d829fde4e4p+0, 0x1.f173d241f23d1p-49, 0x1.9a0f170ca078p+0, 0x1.cdd1884dc6234p-47, 0x1.9c49182a3f08p+0, 0x1.01c7c46b071f3p-48, 0x1.9e86319e3230p+0, 0x1.18c12653c7326p-47, 0x1.a0c667b5de54p+0, 0x1.2594d6d45c656p-47, 0x1.a309bec4a2d0p+0, 0x1.9ac60b8fbb86dp-47, 0x1.a5503b23e254p+0, 0x1.c8b424491caf8p-48, 0x1.a799e1330b34p+0, 0x1.86f2dfb2b158fp-48, 0x1.a9e6b5579fd8p+0, 0x1.fa1f5921deffap-47, 0x1.ac36bbfd3f34p+0, 0x1.ce06dcb351893p-47, 0x1.ae89f995ad38p+0, 0x1.6af439a68bb99p-47, 0x1.b0e07298db64p+0, 0x1.2c8421566fe38p-47, 0x1.b33a2b84f15cp+0, 0x1.d7b5fe873decap-47, 0x1.b59728de5590p+0, 0x1.cc71c40888b24p-47, 0x1.b7f76f2fb5e4p+0, 0x1.baa9ec206ad4fp-50, 0x1.ba5b030a1064p+0, 0x1.30819678d5eb7p-49, 0x1.bcc1e904bc1cp+0, 0x1.2247ba0f45b3dp-48, 0x1.bf2c25bd71e0p+0, 0x1.10811ae04a31cp-49, 0x1.c199bdd85528p+0, 0x1.c2220cb12a092p-48, 0x1.c40ab5fffd04p+0, 0x1.d368a6fc1078cp-47, 0x1.c67f12e57d14p+0, 0x1.694426ffa41e5p-49, 0x1.c8f6d9406e78p+0, 0x1.a88d65e24402ep-47, 0x1.cb720dcef904p+0, 0x1.48a81e5e8f4a5p-47, 0x1.cdf0b555dc3cp+0, 0x1.ce227c4ac7d63p-47, 0x1.d072d4a07894p+0, 0x1.dc68791790d0bp-47, 0x1.d2f87080d89cp+0, 0x1.8c56f091cc4f5p-47, 0x1.d5818dcfba48p+0, 0x1.c976816bad9b8p-50, 0x1.d80e316c9838p+0, 0x1.7bb84f9d04880p-48, 0x1.da9e603db328p+0, 0x1.5c2300696db53p-50, 0x1.dd321f301b44p+0, 0x1.025b4aef1e032p-47, 0x1.dfc97337b9b4p+0, 0x1.eb968cac39ed3p-48, 0x1.e264614f5a10p+0, 0x1.45093b0fd0bd7p-47, 0x1.e502ee78b3fcp+0, 0x1.b139e8980a9cdp-47, 0x1.e7a51fbc74c8p+0, 0x1.a5aa4594191bcp-51, 0x1.ea4afa2a490cp+0, 0x1.9858f73a18f5ep-48, 0x1.ecf482d8e67cp+0, 0x1.846d81897dca5p-47, 0x1.efa1bee615a0p+0, 0x1.3bb8fe90d496dp-47, 0x1.f252b376bba8p+0, 0x1.74e8696fc3639p-48, 0x1.f50765b6e454p+0, 0x1.9d3e12dd8a18bp-54, 0x1.f7bfdad9cbe0p+0, 0x1.38913b4bfe72cp-48, 0x1.fa7c1819e90cp+0, 0x1.82e90a7e74b26p-48, 0x1.fd3c22b8f71cp+0, 0x1.884badd25995ep-47 }; static const long double tiny = 0x1.p-10000L; /* * From Tang's paper: INTRND rounds a floating-point number to the nearest * integer in the manner prescribed by the IEEE standard. It is crucial that * the default round-to-nearest mode, not any other rounding mode, is in * effect here. */ #define INTRND(x) roundl((x)) long double t_expl(long double x) { union IEEEl2bits z; int k, n, n1, n2, sgn; long double r1, r2, r, p, q, t; z.e = x; sgn = z.bits.sign; z.bits.sign = 0; /* * If x = +Inf, then exp(x) = Inf. * If x = -Inf, then exp(x) = 0. * If x = NaN, then exp(x) = NaN. */ if (z.bits.exp == LDBL_MAX_EXP + BIAS) { mask_nbit_l(z); if (!(z.bits.manh | z.bits.manl)) return (sgn ? 0 : x); return (x + x); } if (z.e > THRESHOLD_1) return (x > 0 ? (z.e / (z.e - z.e)) : (0 + tiny * tiny)); if (z.e < THRESHOLD_2) return (1 + x); n = INTRND(x * INV_L); n2 = n % NUM; if (n2 < 0) n2 += NUM; /* Tang's j. */ n1 = n - n2; if (n >= 512 || n <= -512) r1 = (x - n1 * L1) - n2 * L1; else r1 = x - n * L1; r2 = - n * L2; n1 /= NUM; /* Tang's m. */ r = r1 + r2; q = r * r * (P1 + r * (P2 + r * (P3 + r * (P4 + r * (P5 + r * (P6 + r * P7)))))); p = r1 + (r2 + q); t = (long double)s[n2].tail + s[n2].head; t = (long double)s[n2].head + (s[n2].tail + t * p); return (ldexpl(t, n1)); } >Release-Note: >Audit-Trail: >Unformatted: From owner-freebsd-standards@FreeBSD.ORG Sat Nov 20 01:42:41 2010 Return-Path: Delivered-To: freebsd-standards@FreeBSD.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 52C94106564A; Sat, 20 Nov 2010 01:42:41 +0000 (UTC) (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 1767B8FC13; Sat, 20 Nov 2010 01:42:41 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.4/8.14.4) with ESMTP id oAK1ge4B030830; Fri, 19 Nov 2010 17:42:40 -0800 (PST) (envelope-from kargl@troutmask.apl.washington.edu) Received: (from kargl@localhost) by troutmask.apl.washington.edu (8.14.4/8.14.4/Submit) id oAK1geT0030829; Fri, 19 Nov 2010 17:42:40 -0800 (PST) (envelope-from kargl) From: "Steven G. Kargl" Message-Id: <201011200142.oAK1geT0030829@troutmask.apl.washington.edu> In-Reply-To: <201011192310.oAJNA8Ae035514@freefall.freebsd.org> To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Date: Fri, 19 Nov 2010 17:42:40 -0800 (PST) X-Mailer: ELM [version 2.4ME+ PL124d (25)] MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset="US-ASCII" Cc: Subject: Re: standards/152415: [libm] implementation of expl() 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, 20 Nov 2010 01:42:41 -0000 Forgot to mention that this is an ld80 implementation. I have no access to 128 bit long double hardware, so I have no intentions of writing code for ld128. -- Steve http://troutmask.apl.washington.edu/~kargl/ From owner-freebsd-standards@FreeBSD.ORG Sat Nov 20 01:50:10 2010 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 942341065771 for ; Sat, 20 Nov 2010 01:50:10 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 685A78FC08 for ; Sat, 20 Nov 2010 01:50:10 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.4/8.14.4) with ESMTP id oAK1oAoF099130 for ; Sat, 20 Nov 2010 01:50:10 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.4/8.14.4/Submit) id oAK1oAcg099129; Sat, 20 Nov 2010 01:50:10 GMT (envelope-from gnats) Date: Sat, 20 Nov 2010 01:50:10 GMT Message-Id: <201011200150.oAK1oAcg099129@freefall.freebsd.org> To: freebsd-standards@FreeBSD.org From: "Steven G. Kargl" Cc: Subject: Re: standards/152415: [libm] implementation of expl() 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, 20 Nov 2010 01:50:10 -0000 The following reply was made to PR standards/152415; it has been noted by GNATS. From: "Steven G. Kargl" To: FreeBSD-gnats-submit@FreeBSD.org, freebsd-standards@FreeBSD.org Cc: Subject: Re: standards/152415: [libm] implementation of expl() Date: Fri, 19 Nov 2010 17:42:40 -0800 (PST) Forgot to mention that this is an ld80 implementation. I have no access to 128 bit long double hardware, so I have no intentions of writing code for ld128. -- Steve http://troutmask.apl.washington.edu/~kargl/