]>
Commit | Line | Data |
---|---|---|
2aa62f2b | 1 | /** @file\r |
2 | Compute the base 10 logrithm of x.\r | |
3 | \r | |
4 | Copyright (c) 2010 - 2011, Intel Corporation. All rights reserved.<BR>\r | |
5 | This program and the accompanying materials are licensed and made available under\r | |
6 | the terms and conditions of the BSD License that accompanies this distribution.\r | |
7 | The full text of the license may be found at\r | |
8 | http://opensource.org/licenses/bsd-license.\r | |
9 | \r | |
10 | THE PROGRAM IS DISTRIBUTED UNDER THE BSD LICENSE ON AN "AS IS" BASIS,\r | |
11 | WITHOUT WARRANTIES OR REPRESENTATIONS OF ANY KIND, EITHER EXPRESS OR IMPLIED.\r | |
12 | \r | |
13 | * ====================================================\r | |
14 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.\r | |
15 | *\r | |
16 | * Developed at SunPro, a Sun Microsystems, Inc. business.\r | |
17 | * Permission to use, copy, modify, and distribute this\r | |
18 | * software is freely granted, provided that this notice\r | |
19 | * is preserved.\r | |
20 | * ====================================================\r | |
21 | \r | |
22 | e_log10.c 5.1 93/09/24\r | |
23 | NetBSD: e_log10.c,v 1.12 2002/05/26 22:01:51 wiz Exp\r | |
24 | **/\r | |
25 | #include <LibConfig.h>\r | |
26 | #include <sys/EfiCdefs.h>\r | |
27 | \r | |
28 | /* __ieee754_log10(x)\r | |
29 | * Return the base 10 logarithm of x\r | |
30 | *\r | |
31 | * Method :\r | |
32 | * Let log10_2hi = leading 40 bits of log10(2) and\r | |
33 | * log10_2lo = log10(2) - log10_2hi,\r | |
34 | * ivln10 = 1/log(10) rounded.\r | |
35 | * Then\r | |
36 | * n = ilogb(x),\r | |
37 | * if(n<0) n = n+1;\r | |
38 | * x = scalbn(x,-n);\r | |
39 | * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))\r | |
40 | *\r | |
41 | * Note 1:\r | |
42 | * To guarantee log10(10**n)=n, where 10**n is normal, the rounding\r | |
43 | * mode must set to Round-to-Nearest.\r | |
44 | * Note 2:\r | |
45 | * [1/log(10)] rounded to 53 bits has error .198 ulps;\r | |
46 | * log10 is monotonic at all binary break points.\r | |
47 | *\r | |
48 | * Special cases:\r | |
49 | * log10(x) is NaN with signal if x < 0;\r | |
50 | * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;\r | |
51 | * log10(NaN) is that NaN with no signal;\r | |
52 | * log10(10**N) = N for N=0,1,...,22.\r | |
53 | *\r | |
54 | * Constants:\r | |
55 | * The hexadecimal values are the intended ones for the following constants.\r | |
56 | * The decimal values may be used, provided that the compiler will convert\r | |
57 | * from decimal to binary accurately enough to produce the hexadecimal values\r | |
58 | * shown.\r | |
59 | */\r | |
60 | \r | |
61 | #include "math.h"\r | |
62 | #include "math_private.h"\r | |
63 | #include <errno.h>\r | |
64 | \r | |
65 | #if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */\r | |
66 | // potential divide by 0 -- near line 80, (x-x)/zero is on purpose\r | |
67 | #pragma warning ( disable : 4723 )\r | |
68 | #endif\r | |
69 | \r | |
70 | static const double\r | |
71 | two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */\r | |
72 | ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */\r | |
73 | log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */\r | |
74 | log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */\r | |
75 | \r | |
76 | static const double zero = 0.0;\r | |
77 | \r | |
78 | double\r | |
79 | __ieee754_log10(double x)\r | |
80 | {\r | |
81 | double y,z;\r | |
82 | int32_t i,k,hx;\r | |
83 | u_int32_t lx;\r | |
84 | \r | |
85 | EXTRACT_WORDS(hx,lx,x);\r | |
86 | \r | |
87 | k=0;\r | |
88 | if (hx < 0x00100000) { /* x < 2**-1022 */\r | |
89 | if (((hx&0x7fffffff)|lx)==0)\r | |
90 | return -two54/zero; /* log(+-0)=-inf */\r | |
91 | if (hx<0) {\r | |
92 | errno = EDOM;\r | |
93 | return (x-x)/zero; /* log(-#) = NaN */\r | |
94 | }\r | |
95 | k -= 54; x *= two54; /* subnormal number, scale up x */\r | |
96 | GET_HIGH_WORD(hx,x);\r | |
97 | }\r | |
98 | if (hx >= 0x7ff00000) return x+x;\r | |
99 | k += (hx>>20)-1023;\r | |
100 | i = ((u_int32_t)k&0x80000000)>>31;\r | |
101 | hx = (hx&0x000fffff)|((0x3ff-i)<<20);\r | |
102 | y = (double)(k+i);\r | |
103 | SET_HIGH_WORD(x,hx);\r | |
104 | z = y*log10_2lo + ivln10*__ieee754_log(x);\r | |
105 | return z+y*log10_2hi;\r | |
106 | }\r |