1331722Seadler/* 21573Srgrimes * Copyright (c) 1985, 1993 31573Srgrimes * The Regents of the University of California. All rights reserved. 41573Srgrimes * 51573Srgrimes * Redistribution and use in source and binary forms, with or without 61573Srgrimes * modification, are permitted provided that the following conditions 71573Srgrimes * are met: 81573Srgrimes * 1. Redistributions of source code must retain the above copyright 91573Srgrimes * notice, this list of conditions and the following disclaimer. 101573Srgrimes * 2. Redistributions in binary form must reproduce the above copyright 111573Srgrimes * notice, this list of conditions and the following disclaimer in the 121573Srgrimes * documentation and/or other materials provided with the distribution. 131573Srgrimes * 3. All advertising materials mentioning features or use of this software 141573Srgrimes * must display the following acknowledgement: 151573Srgrimes * This product includes software developed by the University of 161573Srgrimes * California, Berkeley and its contributors. 171573Srgrimes * 4. Neither the name of the University nor the names of its contributors 181573Srgrimes * may be used to endorse or promote products derived from this software 191573Srgrimes * without specific prior written permission. 201573Srgrimes * 211573Srgrimes * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 221573Srgrimes * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 231573Srgrimes * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 241573Srgrimes * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 251573Srgrimes * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 261573Srgrimes * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 271573Srgrimes * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 281573Srgrimes * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 291573Srgrimes * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 301573Srgrimes * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 311573Srgrimes * SUCH DAMAGE. 321573Srgrimes */ 331573Srgrimes 34176449Sdas/* @(#)exp.c 8.1 (Berkeley) 6/4/93 */ 3592887Sobrien#include <sys/cdefs.h> 3692887Sobrien__FBSDID("$FreeBSD$"); 371573Srgrimes 3892887Sobrien 391573Srgrimes/* EXP(X) 401573Srgrimes * RETURN THE EXPONENTIAL OF X 411573Srgrimes * DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS) 428870Srgrimes * CODED IN C BY K.C. NG, 1/19/85; 431573Srgrimes * REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86. 441573Srgrimes * 451573Srgrimes * Required system supported functions: 468870Srgrimes * scalb(x,n) 478870Srgrimes * copysign(x,y) 481573Srgrimes * finite(x) 491573Srgrimes * 501573Srgrimes * Method: 518870Srgrimes * 1. Argument Reduction: given the input x, find r and integer k such 521573Srgrimes * that 538870Srgrimes * x = k*ln2 + r, |r| <= 0.5*ln2 . 541573Srgrimes * r will be represented as r := z+c for better accuracy. 551573Srgrimes * 568870Srgrimes * 2. Compute exp(r) by 571573Srgrimes * 581573Srgrimes * exp(r) = 1 + r + r*R1/(2-R1), 591573Srgrimes * where 601573Srgrimes * R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))). 611573Srgrimes * 621573Srgrimes * 3. exp(x) = 2^k * exp(r) . 631573Srgrimes * 641573Srgrimes * Special cases: 651573Srgrimes * exp(INF) is INF, exp(NaN) is NaN; 661573Srgrimes * exp(-INF)= 0; 671573Srgrimes * for finite argument, only exp(0)=1 is exact. 681573Srgrimes * 691573Srgrimes * Accuracy: 701573Srgrimes * exp(x) returns the exponential of x nearly rounded. In a test run 711573Srgrimes * with 1,156,000 random arguments on a VAX, the maximum observed 721573Srgrimes * error was 0.869 ulps (units in the last place). 731573Srgrimes */ 741573Srgrimes 751573Srgrimes#include "mathimpl.h" 761573Srgrimes 77226414Sdasstatic const double p1 = 0x1.555555555553ep-3; 78226414Sdasstatic const double p2 = -0x1.6c16c16bebd93p-9; 79226414Sdasstatic const double p3 = 0x1.1566aaf25de2cp-14; 80226414Sdasstatic const double p4 = -0x1.bbd41c5d26bf1p-20; 81226414Sdasstatic const double p5 = 0x1.6376972bea4d0p-25; 82226414Sdasstatic const double ln2hi = 0x1.62e42fee00000p-1; 83226414Sdasstatic const double ln2lo = 0x1.a39ef35793c76p-33; 84226414Sdasstatic const double lnhuge = 0x1.6602b15b7ecf2p9; 85226414Sdasstatic const double lntiny = -0x1.77af8ebeae354p9; 86226414Sdasstatic const double invln2 = 0x1.71547652b82fep0; 871573Srgrimes 8893211Sbde#if 0 891573Srgrimesdouble exp(x) 901573Srgrimesdouble x; 911573Srgrimes{ 921573Srgrimes double z,hi,lo,c; 931573Srgrimes int k; 941573Srgrimes 951573Srgrimes#if !defined(vax)&&!defined(tahoe) 961573Srgrimes if(x!=x) return(x); /* x is NaN */ 971573Srgrimes#endif /* !defined(vax)&&!defined(tahoe) */ 981573Srgrimes if( x <= lnhuge ) { 991573Srgrimes if( x >= lntiny ) { 1001573Srgrimes 1011573Srgrimes /* argument reduction : x --> x - k*ln2 */ 1021573Srgrimes 1031573Srgrimes k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */ 1041573Srgrimes 1051573Srgrimes /* express x-k*ln2 as hi-lo and let x=hi-lo rounded */ 1061573Srgrimes 1071573Srgrimes hi=x-k*ln2hi; 1081573Srgrimes x=hi-(lo=k*ln2lo); 1091573Srgrimes 1101573Srgrimes /* return 2^k*[1+x+x*c/(2+c)] */ 1111573Srgrimes z=x*x; 1121573Srgrimes c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); 1131573Srgrimes return scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k); 1141573Srgrimes 1151573Srgrimes } 1161573Srgrimes /* end of x > lntiny */ 1171573Srgrimes 1188870Srgrimes else 1191573Srgrimes /* exp(-big#) underflows to zero */ 1201573Srgrimes if(finite(x)) return(scalb(1.0,-5000)); 1211573Srgrimes 1221573Srgrimes /* exp(-INF) is zero */ 1231573Srgrimes else return(0.0); 1241573Srgrimes } 1251573Srgrimes /* end of x < lnhuge */ 1261573Srgrimes 1278870Srgrimes else 1281573Srgrimes /* exp(INF) is INF, exp(+big#) overflows to INF */ 1291573Srgrimes return( finite(x) ? scalb(1.0,5000) : x); 1301573Srgrimes} 13193211Sbde#endif 1321573Srgrimes 1331573Srgrimes/* returns exp(r = x + c) for |c| < |x| with no overlap. */ 1341573Srgrimes 1351573Srgrimesdouble __exp__D(x, c) 1361573Srgrimesdouble x, c; 1371573Srgrimes{ 138138924Sdas double z,hi,lo; 1391573Srgrimes int k; 1401573Srgrimes 141138924Sdas if (x != x) /* x is NaN */ 142138924Sdas return(x); 1431573Srgrimes if ( x <= lnhuge ) { 1441573Srgrimes if ( x >= lntiny ) { 1451573Srgrimes 1461573Srgrimes /* argument reduction : x --> x - k*ln2 */ 1471573Srgrimes z = invln2*x; 1481573Srgrimes k = z + copysign(.5, x); 1491573Srgrimes 1501573Srgrimes /* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */ 1511573Srgrimes 1521573Srgrimes hi=(x-k*ln2hi); /* Exact. */ 1531573Srgrimes x= hi - (lo = k*ln2lo-c); 1541573Srgrimes /* return 2^k*[1+x+x*c/(2+c)] */ 1551573Srgrimes z=x*x; 1561573Srgrimes c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5)))); 1571573Srgrimes c = (x*c)/(2.0-c); 1581573Srgrimes 1591573Srgrimes return scalb(1.+(hi-(lo - c)), k); 1601573Srgrimes } 1611573Srgrimes /* end of x > lntiny */ 1621573Srgrimes 1638870Srgrimes else 1641573Srgrimes /* exp(-big#) underflows to zero */ 1651573Srgrimes if(finite(x)) return(scalb(1.0,-5000)); 1661573Srgrimes 1671573Srgrimes /* exp(-INF) is zero */ 1681573Srgrimes else return(0.0); 1691573Srgrimes } 1701573Srgrimes /* end of x < lnhuge */ 1711573Srgrimes 1728870Srgrimes else 1731573Srgrimes /* exp(INF) is INF, exp(+big#) overflows to INF */ 1741573Srgrimes return( finite(x) ? scalb(1.0,5000) : x); 1751573Srgrimes} 176