Logo ROOT   6.08/07
Reference Guide
log.h
Go to the documentation of this file.
1 /*
2  * log.h
3  * The basic idea is to exploit Pade polynomials.
4  * A lot of ideas were inspired by the cephes math library (by Stephen L. Moshier
5  * moshier@na-net.ornl.gov) as well as actual code.
6  * The Cephes library can be found here: http://www.netlib.org/cephes/
7  *
8  * Created on: Jun 23, 2012
9  * Author: Danilo Piparo, Thomas Hauth, Vincenzo Innocente
10  */
11 
12 /*
13  * VDT is free software: you can redistribute it and/or modify
14  * it under the terms of the GNU Lesser Public License as published by
15  * the Free Software Foundation, either version 3 of the License, or
16  * (at your option) any later version.
17  *
18  * This program is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU Lesser Public License for more details.
22  *
23  * You should have received a copy of the GNU Lesser Public License
24  * along with this program. If not, see <http://www.gnu.org/licenses/>.
25  */
26 
27 #ifndef LOG_H_
28 #define LOG_H_
29 
30 #include "vdtcore_common.h"
31 #include <limits>
32 
33 namespace vdt{
34 
35 // local namespace for the constants/functions which are necessary only here
36 namespace details{
37 
38 const double LOG_UPPER_LIMIT = 1e307;
39 const double LOG_LOWER_LIMIT = 0;
40 
41 const double SQRTH = 0.70710678118654752440;
42 
43 inline double get_log_px(const double x){
44  const double PX1log = 1.01875663804580931796E-4;
45  const double PX2log = 4.97494994976747001425E-1;
46  const double PX3log = 4.70579119878881725854E0;
47  const double PX4log = 1.44989225341610930846E1;
48  const double PX5log = 1.79368678507819816313E1;
49  const double PX6log = 7.70838733755885391666E0;
50 
51  double px = PX1log;
52  px *= x;
53  px += PX2log;
54  px *= x;
55  px += PX3log;
56  px *= x;
57  px += PX4log;
58  px *= x;
59  px += PX5log;
60  px *= x;
61  px += PX6log;
62  return px;
63 
64 }
65 
66 inline double get_log_qx(const double x){
67  const double QX1log = 1.12873587189167450590E1;
68  const double QX2log = 4.52279145837532221105E1;
69  const double QX3log = 8.29875266912776603211E1;
70  const double QX4log = 7.11544750618563894466E1;
71  const double QX5log = 2.31251620126765340583E1;
72 
73  double qx = x;
74  qx += QX1log;
75  qx *=x;
76  qx += QX2log;
77  qx *=x;
78  qx += QX3log;
79  qx *=x;
80  qx += QX4log;
81  qx *=x;
82  qx += QX5log;
83  return qx;
84 }
85 
86 }
87 
88 // Log double precision --------------------------------------------------------
89 inline double fast_log(double x){
90 
91  const double original_x = x;
92 
93  /* separate mantissa from exponent */
94  double fe;
95  x = details::getMantExponent(x,fe);
96 
97  // blending
98  x > details::SQRTH? fe+=1. : x+=x ;
99  x -= 1.0;
100 
101  /* rational form */
102  double px = details::get_log_px(x);
103 
104  //for the final formula
105  const double x2 = x*x;
106  px *= x;
107  px *= x2;
108 
109  const double qx = details::get_log_qx(x);
110 
111  double res = px / qx ;
112 
113  res -= fe * 2.121944400546905827679e-4;
114  res -= 0.5 * x2 ;
115 
116  res = x + res;
117  res += fe * 0.693359375;
118 
119  if (original_x > details::LOG_UPPER_LIMIT)
120  res = std::numeric_limits<double>::infinity();
121  if (original_x < details::LOG_LOWER_LIMIT) // THIS IS NAN!
122  res = - std::numeric_limits<double>::quiet_NaN();
123 
124  return res;
125 
126 }
127 
128 // Log single precision --------------------------------------------------------
129 
130 
131 
132 namespace details{
133 
134 const float LOGF_UPPER_LIMIT = MAXNUMF;
135 const float LOGF_LOWER_LIMIT = 0;
136 
137 const float PX1logf = 7.0376836292E-2f;
138 const float PX2logf = -1.1514610310E-1f;
139 const float PX3logf = 1.1676998740E-1f;
140 const float PX4logf = -1.2420140846E-1f;
141 const float PX5logf = 1.4249322787E-1f;
142 const float PX6logf = -1.6668057665E-1f;
143 const float PX7logf = 2.0000714765E-1f;
144 const float PX8logf = -2.4999993993E-1f;
145 const float PX9logf = 3.3333331174E-1f;
146 
147 inline float get_log_poly(const float x){
148  float y = x*PX1logf;
149  y += PX2logf;
150  y *= x;
151  y += PX3logf;
152  y *= x;
153  y += PX4logf;
154  y *= x;
155  y += PX5logf;
156  y *= x;
157  y += PX6logf;
158  y *= x;
159  y += PX7logf;
160  y *= x;
161  y += PX8logf;
162  y *= x;
163  y += PX9logf;
164  return y;
165 }
166 
167 const float SQRTHF = 0.707106781186547524f;
168 
169 }
170 
171 // Log single precision --------------------------------------------------------
172 inline float fast_logf( float x ) {
173 
174  const float original_x = x;
175 
176  float fe;
177  x = details::getMantExponentf( x, fe);
178 
179  x > details::SQRTHF? fe+=1.f : x+=x ;
180  x -= 1.0f;
181 
182  const float x2 = x*x;
183 
184  float res = details::get_log_poly(x);
185  res *= x2*x;
186 
187  res += -2.12194440e-4f * fe;
188  res += -0.5f * x2;
189 
190  res= x + res;
191 
192  res += 0.693359375f * fe;
193 
194  if (original_x > details::LOGF_UPPER_LIMIT)
195  res = std::numeric_limits<float>::infinity();
196  if (original_x < details::LOGF_LOWER_LIMIT)
197  res = -std::numeric_limits<float>::quiet_NaN();
198 
199  return res;
200 }
201 
202 
203 //------------------------------------------------------------------------------
204 
205 // void logv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
206 // void fast_logv(const uint32_t size, double const * __restrict__ iarray, double* __restrict__ oarray);
207 // void logfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
208 // void fast_logfv(const uint32_t size, float const * __restrict__ iarray, float* __restrict__ oarray);
209 
210 } //vdt namespace
211 
212 #endif /* LOG_H_ */
const double SQRTH
Definition: log.h:41
const float PX6logf
Definition: log.h:142
const float PX9logf
Definition: log.h:145
double get_log_px(const double x)
Definition: log.h:43
const float PX8logf
Definition: log.h:144
double get_log_qx(const double x)
Definition: log.h:66
const float MAXNUMF
const float LOGF_LOWER_LIMIT
Definition: log.h:135
const double LOG_UPPER_LIMIT
Definition: log.h:38
const float PX5logf
Definition: log.h:141
const float PX1logf
Definition: log.h:137
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
float getMantExponentf(const float x, float &fe)
Like frexp but vectorising and the exponent is a float.
float get_log_poly(const float x)
Definition: log.h:147
double fast_log(double x)
Definition: log.h:89
const float PX3logf
Definition: log.h:139
const double LOG_LOWER_LIMIT
Definition: log.h:39
const float PX7logf
Definition: log.h:143
const float PX4logf
Definition: log.h:140
double f(double x)
const float PX2logf
Definition: log.h:138
Double_t y[n]
Definition: legend1.C:17
Definition: asin.h:32
const float SQRTHF
Definition: log.h:167
const float LOGF_UPPER_LIMIT
Definition: log.h:134
float fast_logf(float x)
Definition: log.h:172
double getMantExponent(const double x, double &fe)
Like frexp but vectorising and the exponent is a double.