1 
/* * Mode: C; tabwidth: 8; indenttabsmode: nil; cbasicoffset: 4 * 
2 
* 
3 
* ***** BEGIN LICENSE BLOCK ***** 
4 
* Version: MPL 1.1/GPL 2.0/LGPL 2.1 
5 
* 
6 
* The contents of this file are subject to the Mozilla Public License Version 
7 
* 1.1 (the "License"); you may not use this file except in compliance with 
8 
* the License. You may obtain a copy of the License at 
9 
* http://www.mozilla.org/MPL/ 
10 
* 
11 
* Software distributed under the License is distributed on an "AS IS" basis, 
12 
* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License 
13 
* for the specific language governing rights and limitations under the 
14 
* License. 
15 
* 
16 
* The Original Code is Mozilla Communicator client code, released 
17 
* March 31, 1998. 
18 
* 
19 
* The Initial Developer of the Original Code is 
20 
* Sun Microsystems, Inc. 
21 
* Portions created by the Initial Developer are Copyright (C) 1998 
22 
* the Initial Developer. All Rights Reserved. 
23 
* 
24 
* Contributor(s): 
25 
* 
26 
* Alternatively, the contents of this file may be used under the terms of 
27 
* either of the GNU General Public License Version 2 or later (the "GPL"), 
28 
* or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"), 
29 
* in which case the provisions of the GPL or the LGPL are applicable instead 
30 
* of those above. If you wish to allow use of your version of this file only 
31 
* under the terms of either the GPL or the LGPL, and not to allow others to 
32 
* use your version of this file under the terms of the MPL, indicate your 
33 
* decision by deleting the provisions above and replace them with the notice 
34 
* and other provisions required by the GPL or the LGPL. If you do not delete 
35 
* the provisions above, a recipient may use your version of this file under 
36 
* the terms of any one of the MPL, the GPL or the LGPL. 
37 
* 
38 
* ***** END LICENSE BLOCK ***** */ 
39 

40 
/* @(#)e_log10.c 1.3 95/01/18 */ 
41 
/* 
42 
* ==================================================== 
43 
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 
44 
* 
45 
* Developed at SunSoft, a Sun Microsystems, Inc. business. 
46 
* Permission to use, copy, modify, and distribute this 
47 
* software is freely granted, provided that this notice 
48 
* is preserved. 
49 
* ==================================================== 
50 
*/ 
51 

52 
/* __ieee754_log10(x) 
53 
* Return the base 10 logarithm of x 
54 
* 
55 
* Method : 
56 
* Let log10_2hi = leading 40 bits of log10(2) and 
57 
* log10_2lo = log10(2)  log10_2hi, 
58 
* ivln10 = 1/log(10) rounded. 
59 
* Then 
60 
* n = ilogb(x), 
61 
* if(n<0) n = n+1; 
62 
* x = scalbn(x,n); 
63 
* log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) 
64 
* 
65 
* Note 1: 
66 
* To guarantee log10(10**n)=n, where 10**n is normal, the rounding 
67 
* mode must set to RoundtoNearest. 
68 
* Note 2: 
69 
* [1/log(10)] rounded to 53 bits has error .198 ulps; 
70 
* log10 is monotonic at all binary break points. 
71 
* 
72 
* Special cases: 
73 
* log10(x) is NaN with signal if x < 0; 
74 
* log10(+INF) is +INF with no signal; log10(0) is INF with signal; 
75 
* log10(NaN) is that NaN with no signal; 
76 
* log10(10**N) = N for N=0,1,...,22. 
77 
* 
78 
* Constants: 
79 
* The hexadecimal values are the intended ones for the following constants. 
80 
* The decimal values may be used, provided that the compiler will convert 
81 
* from decimal to binary accurately enough to produce the hexadecimal values 
82 
* shown. 
83 
*/ 
84 

85 
#include "fdlibm.h" 
86 

87 
#ifdef __STDC__ 
88 
static const double 
89 
#else 
90 
static double 
91 
#endif 
92 
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 
93 
ivln10 = 4.34294481903251816668e01, /* 0x3FDBCB7B, 0x1526E50E */ 
94 
log10_2hi = 3.01029995663611771306e01, /* 0x3FD34413, 0x509F6000 */ 
95 
log10_2lo = 3.69423907715893078616e13; /* 0x3D59FEF3, 0x11F12B36 */ 
96 

97 
static double zero = 0.0; 
98 

99 
#ifdef __STDC__ 
100 
double __ieee754_log10(double x) 
101 
#else 
102 
double __ieee754_log10(x) 
103 
double x; 
104 
#endif 
105 
{ 
106 
fd_twoints u; 
107 
double y,z; 
108 
int i,k,hx; 
109 
unsigned lx; 
110 

111 
u.d = x; 
112 
hx = __HI(u); /* high word of x */ 
113 
lx = __LO(u); /* low word of x */ 
114 

115 
k=0; 
116 
if (hx < 0x00100000) { /* x < 2**1022 */ 
117 
if (((hx&0x7fffffff)lx)==0) 
118 
return two54/zero; /* log(+0)=inf */ 
119 
if (hx<0) return (xx)/zero; /* log(#) = NaN */ 
120 
k = 54; x *= two54; /* subnormal number, scale up x */ 
121 
u.d = x; 
122 
hx = __HI(u); /* high word of x */ 
123 
} 
124 
if (hx >= 0x7ff00000) return x+x; 
125 
k += (hx>>20)1023; 
126 
i = ((unsigned)k&0x80000000)>>31; 
127 
hx = (hx&0x000fffff)((0x3ffi)<<20); 
128 
y = (double)(k+i); 
129 
u.d = x; 
130 
__HI(u) = hx; 
131 
x = u.d; 
132 
z = y*log10_2lo + ivln10*__ieee754_log(x); 
133 
return z+y*log10_2hi; 
134 
} 