From owner-svn-src-head@FreeBSD.ORG Fri Oct 19 22:47:45 2012 Return-Path: Delivered-To: svn-src-head@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 6DD69430; Fri, 19 Oct 2012 22:47:45 +0000 (UTC) (envelope-from imp@FreeBSD.org) Received: from svn.freebsd.org (svn.freebsd.org [IPv6:2001:4f8:fff6::2c]) by mx1.freebsd.org (Postfix) with ESMTP id 3CF788FC08; Fri, 19 Oct 2012 22:47:45 +0000 (UTC) Received: from svn.freebsd.org (localhost [127.0.0.1]) by svn.freebsd.org (8.14.4/8.14.4) with ESMTP id q9JMlj1O093100; Fri, 19 Oct 2012 22:47:45 GMT (envelope-from imp@svn.freebsd.org) Received: (from imp@localhost) by svn.freebsd.org (8.14.4/8.14.4/Submit) id q9JMljL4093098; Fri, 19 Oct 2012 22:47:45 GMT (envelope-from imp@svn.freebsd.org) Message-Id: <201210192247.q9JMljL4093098@svn.freebsd.org> From: Warner Losh Date: Fri, 19 Oct 2012 22:47:45 +0000 (UTC) To: src-committers@freebsd.org, svn-src-all@freebsd.org, svn-src-head@freebsd.org Subject: svn commit: r241756 - head/lib/msun/src X-SVN-Group: head MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit X-BeenThere: svn-src-head@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: SVN commit messages for the src tree for head/-current List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 19 Oct 2012 22:47:45 -0000 Author: imp Date: Fri Oct 19 22:47:44 2012 New Revision: 241756 URL: http://svn.freebsd.org/changeset/base/241756 Log: Document the method used to compute expf. Taken from exp, with changes to reflect differences in computation between the two. Modified: head/lib/msun/src/e_expf.c Modified: head/lib/msun/src/e_expf.c ============================================================================== --- head/lib/msun/src/e_expf.c Fri Oct 19 22:46:48 2012 (r241755) +++ head/lib/msun/src/e_expf.c Fri Oct 19 22:47:44 2012 (r241756) @@ -21,6 +21,68 @@ __FBSDID("$FreeBSD$"); #include "math.h" #include "math_private.h" +/* __ieee754_expf + * Returns the exponential of x. + * + * Method + * 1. Argument reduction: + * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2. + * + * Here r will be represented as r = hi-lo for better + * accuracy. + * + * 2. Approximation of exp(r) by a special rational function on + * the interval [0,0.34658]: + * Write + * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... + * We use a special Remes algorithm on [0,0.34658] to generate + * a polynomial of degree 2 to approximate R. The maximum error + * of this polynomial approximation is bounded by 2**-27. In + * other words, + * R(z) ~ 2.0 + P1*z + P2*z*z + * (where z=r*r, and the values of P1 and P2 are listed below) + * and + * | 2 | -27 + * | 2.0+P1*z+P2*z - R(z) | <= 2 + * | | + * The computation of expf(r) thus becomes + * 2*r + * expf(r) = 1 + ------- + * R - r + * r*R1(r) + * = 1 + r + ----------- (for better accuracy) + * 2 - R1(r) + * where + * 2 4 + * R1(r) = r - (P1*r + P2*r) + * + * 3. Scale back to obtain expf(x): + * From step 1, we have + * expf(x) = 2^k * expf(r) + * + * Special cases: + * expf(INF) is INF, exp(NaN) is NaN; + * expf(-INF) is 0, and + * for finite argument, only exp(0)=1 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 0.5013 ulp (unit in the last place). + * + * Misc. info. + * For IEEE float + * if x > 8.8721679688e+01 then exp(x) overflow + * if x < -1.0397208405e+02 then exp(x) underflow + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ static const float one = 1.0, halF[2] = {0.5,-0.5,},