1 //===========================================================================
5 // Part of the standard mathematical function library
7 //===========================================================================
8 //####ECOSGPLCOPYRIGHTBEGIN####
9 // -------------------------------------------
10 // This file is part of eCos, the Embedded Configurable Operating System.
11 // Copyright (C) 1998, 1999, 2000, 2001, 2002 Red Hat, Inc.
13 // eCos is free software; you can redistribute it and/or modify it under
14 // the terms of the GNU General Public License as published by the Free
15 // Software Foundation; either version 2 or (at your option) any later version.
17 // eCos is distributed in the hope that it will be useful, but WITHOUT ANY
18 // WARRANTY; without even the implied warranty of MERCHANTABILITY or
19 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 // You should have received a copy of the GNU General Public License along
23 // with eCos; if not, write to the Free Software Foundation, Inc.,
24 // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
26 // As a special exception, if other files instantiate templates or use macros
27 // or inline functions from this file, or you compile this file and link it
28 // with other works to produce a work based on this file, this file does not
29 // by itself cause the resulting work to be covered by the GNU General Public
30 // License. However the source code for this file must still be made available
31 // in accordance with section (3) of the GNU General Public License.
33 // This exception does not invalidate any other reasons why a work based on
34 // this file might be covered by the GNU General Public License.
36 // Alternative licenses for eCos may be arranged by contacting Red Hat, Inc.
37 // at http://sources.redhat.com/ecos/ecos-license/
38 // -------------------------------------------
39 //####ECOSGPLCOPYRIGHTEND####
40 //===========================================================================
41 //#####DESCRIPTIONBEGIN####
43 // Author(s): jlarmour
44 // Contributors: jlarmour
50 //####DESCRIPTIONEND####
52 //===========================================================================
56 #include <pkgconf/libm.h> // Configuration header
58 // Include the Math library?
61 // Derived from code with the following copyright
64 /* @(#)e_log.c 1.3 95/01/18 */
66 * ====================================================
67 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
69 * Developed at SunSoft, a Sun Microsystems, Inc. business.
70 * Permission to use, copy, modify, and distribute this
71 * software is freely granted, provided that this notice
73 * ====================================================
77 * Return the logrithm of x
80 * 1. Argument Reduction: find k and f such that
82 * where sqrt(2)/2 < 1+f < sqrt(2) .
84 * 2. Approximation of log(1+f).
85 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
86 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
88 * We use a special Reme algorithm on [0,0.1716] to generate
89 * a polynomial of degree 14 to approximate R The maximum error
90 * of this polynomial approximation is bounded by 2**-58.45. In
93 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
94 * (the values of Lg1 to Lg7 are listed in the program)
97 * | Lg1*s +...+Lg7*s - R(z) | <= 2
99 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
100 * In order to guarantee error in log below 1ulp, we compute log
102 * log(1+f) = f - s*(f - R) (if f is not too large)
103 * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
105 * 3. Finally, log(x) = k*ln2 + log(1+f).
106 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
107 * Here ln2 is split into two floating point number:
109 * where n*ln2_hi is always exact for |n| < 2000.
112 * log(x) is NaN with signal if x < 0 (including -INF) ;
113 * log(+INF) is +INF; log(0) is -INF with signal;
114 * log(NaN) is that NaN with no signal.
117 * according to an error analysis, the error is always less than
118 * 1 ulp (unit in the last place).
121 * The hexadecimal values are the intended ones for the following
122 * constants. The decimal values may be used, provided that the
123 * compiler will convert from decimal to binary accurately enough
124 * to produce the hexadecimal values shown.
127 #include "mathincl/fdlibm.h"
130 ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
131 ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
132 two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
133 Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
134 Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
135 Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
136 Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
137 Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
138 Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
139 Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
141 static double zero = 0.0;
143 double __ieee754_log(double x)
145 double hfsq,f,s,z,R,w,t1,t2,dk;
149 hx = CYG_LIBM_HI(x); /* high word of x */
150 lx = CYG_LIBM_LO(x); /* low word of x */
153 if (hx < 0x00100000) { /* x < 2**-1022 */
154 if (((hx&0x7fffffff)|lx)==0)
155 return -two54/zero; /* log(+-0)=-inf */
156 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
157 k -= 54; x *= two54; /* subnormal number, scale up x */
158 hx = CYG_LIBM_HI(x); /* high word of x */
160 if (hx >= 0x7ff00000) return x+x;
163 i = (hx+0x95f64)&0x100000;
164 CYG_LIBM_HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
167 if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
169 if(k==0) return zero;
172 return dk*ln2_hi+dk*ln2_lo;
175 R = f*f*(0.5-0.33333333333333333*f);
179 return dk*ln2_hi-((R-dk*ln2_lo)-f);
188 t1= w*(Lg2+w*(Lg4+w*Lg6));
189 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
194 if(k==0) return f-(hfsq-s*(hfsq+R)); else
195 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
197 if(k==0) return f-s*(f-R); else
198 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
202 #endif // ifdef CYGPKG_LIBM