RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
powfast.cpp
Go to the documentation of this file.
1 /*------------------------------------------------------------------------------
2 
3  HXA7241 General library.
4  Copyright (c) 2005-2007, Harrison Ainsworth / HXA7241.
5 
6  http://www.hxa7241.org/
7 
8 ------------------------------------------------------------------------------*/
9 
10 #include "powfast.hpp"
11 #include <math.h>
12 
14 
15 namespace {
16 
40 const float _2p23 = 8388608.0f;
41 
42 /*
43  * Initialize powFast lookup table.
44  *
45  * @pTable length must be 2 ^ precision
46  * @precision number of mantissa bits used, >= 0 and <= 18
47  */
48 void
49 powFastSetTable(unsigned int* const pTable, const unsigned int precision)
50 {
51  // step along table elements and x-axis positions
52  float zeroToOne = 1.0f / (static_cast<float>(1 << precision) * 2.0f);
53  for (int i = 0; i < (1 << precision); ++i) {
54  // make y-axis value for table element
55  const float f = (::powf(2.0f, zeroToOne) - 1.0f) * _2p23;
56  pTable[i] = static_cast<unsigned int>(f < _2p23 ? f : (_2p23 - 1.0f));
57  zeroToOne += 1.0f / static_cast<float>(1 << precision);
58  }
59 }
60 
61 /*
62  * Get pow (fast!).
63  *
64  * @val power to raise radix to
65  * @ilog2 one over log, to required radix, of two
66  * @pTable length must be 2 ^ precision
67  * @precision number of mantissa bits used, >= 0 and <= 18
68  */
69 inline float
70 powFastLookup(const float val, const float ilog2, unsigned int* const pTable,
71  const unsigned int precision)
72 
73 {
74  // build float bits
75  const int i = static_cast<int>((val * (_2p23 * ilog2)) + (127.0f * _2p23));
76 
77  // replace mantissa with lookup
78  const int it = (i & 0xFF800000) | pTable[(i & 0x7FFFFF) >> (23 - precision)];
79 
80  // convert bits to float
81  return *reinterpret_cast<const float*>(&it);
82 }
83 }
84 
86 
87 PowFast::PowFast(const unsigned int precision)
88  : precision_m(precision <= 18u ? precision : 18u)
89  , pTable_m(new unsigned int[1 << precision_m])
90 {
91  powFastSetTable(pTable_m, precision_m);
92 }
93 
95 {
96  delete[] pTable_m;
97 }
98 
99 float
100 PowFast::two(const float f) const
101 {
102  return powFastLookup(f, 1.0f, pTable_m, precision_m);
103 }
104 
105 float
106 PowFast::e(const float f) const
107 {
108  return powFastLookup(f, 1.44269504088896f, pTable_m, precision_m);
109 }
110 
111 float
112 PowFast::ten(const float f) const
113 {
114  return powFastLookup(f, 3.32192809488736f, pTable_m, precision_m);
115 }
116 
117 float
118 PowFast::r(const float logr, const float f) const
119 {
120  return powFastLookup(f, (logr * 1.44269504088896f), pTable_m, precision_m);
121 }
122 
123 unsigned int
125 {
126  return precision_m;
127 }
128 
130 const PowFast&
132 {
133  static const PowFast k(11);
134  return k;
135 }
136 
138 #ifdef TESTING
139 
140 #include <float.h>
141 #include <iomanip>
142 #include <ostream>
143 #include <stdlib.h>
144 
145 typedef unsigned int udword;
146 typedef int dword;
147 
148 namespace {
149 
150 /*
151  * rand [0,1) (may be coarsely quantized).
152  */
153 float
154 randFloat()
155 {
156  return static_cast<float>(static_cast<udword>(::rand())) /
157  static_cast<float>(static_cast<udword>(RAND_MAX) + 1);
158 }
159 }
160 
161 bool
162 test_PowFast(std::ostream* pOut, const bool isVerbose, const dword seed)
163 {
164  bool isOk = true;
165 
166  if (pOut)
167  *pOut << "[ test_PowFast ]\n\n";
168 
169  if (pOut) {
170  pOut->setf(std::ios_base::scientific, std::ios_base::floatfield);
171  pOut->precision(6);
172  }
173 
174  ::srand(static_cast<unsigned>(seed));
175 
177  {
178  const PowFast powFastAdj(11);
179 
180  const dword J_COUNT = 1000;
181 
182  float sumDifE = 0.0f;
183  float maxDifE = static_cast<float>(FLT_MIN);
184 
185  const dword E_START = -86;
186  const dword E_END = +88;
187  for (dword i = E_START; i < E_END; ++i) {
188  for (dword j = 0; j < J_COUNT; ++j) {
189  const float number = static_cast<float>(i) + randFloat();
190  const float pe = ::powf(2.71828182845905f, number);
191  const float pef = powFastAdj.e(number);
192  const float peDif = ::fabsf(pef - pe) / pe;
193  sumDifE += peDif;
194  maxDifE = (maxDifE >= peDif) ? maxDifE : peDif;
195 
196  if ((0 == j) && pOut && isVerbose) {
197  *pOut << std::showpos << std::fixed << std::setw(10) << number
198  << std::scientific << std::noshowpos;
199  *pOut << " E " << pef << " " << peDif << "\n";
200  }
201  }
202  }
203 
204  if (pOut && isVerbose)
205  *pOut << "\n";
206 
207  float sumDifT = 0.0f;
208  float maxDifT = static_cast<float>(FLT_MIN);
209  const dword T_START = -36;
210  const dword T_END = +38;
211 
212  for (dword i = T_START; i < T_END; ++i) {
213  for (dword j = 0; j < J_COUNT; ++j) {
214  const float number = static_cast<float>(i) + randFloat();
215  const float pt = ::powf(10.0f, number);
216  const float ptf = powFastAdj.ten(number);
217  const float ptDif = ::fabsf(ptf - pt) / pt;
218  sumDifT += ptDif;
219  maxDifT = (maxDifT >= ptDif) ? maxDifT : ptDif;
220 
221  if ((0 == j) && pOut && isVerbose) {
222  *pOut << std::showpos << std::fixed << std::setw(10) << number
223  << std::scientific << std::noshowpos;
224  *pOut << " T " << ptf << " " << ptDif << "\n";
225  }
226  }
227  }
228 
229  if (pOut && isVerbose)
230  *pOut << "\n";
231 
232  const float meanDifE = sumDifE / (static_cast<float>(E_END - E_START) *
233  static_cast<float>(J_COUNT));
234  const float meanDifT = sumDifT / (static_cast<float>(T_END - T_START) *
235  static_cast<float>(J_COUNT));
236 
237  if (pOut && isVerbose)
238  *pOut << "precision: " << powFastAdj.precision() << "\n";
239  if (pOut && isVerbose)
240  *pOut << "mean diff, E: " << meanDifE << " 10: " << meanDifT << "\n";
241  if (pOut && isVerbose)
242  *pOut << "max diff, E: " << maxDifE << " 10: " << maxDifT << "\n";
243  if (pOut && isVerbose)
244  *pOut << "\n";
245 
246  bool isOk_ = (meanDifE < 0.0001f) & (meanDifT < 0.0001f) &
247  (maxDifE < 0.0002f) & (maxDifT < 0.0002f);
248 
249  if (pOut)
250  *pOut << "adjustable : " << (isOk_ ? "--- succeeded" : "*** failed")
251  << "\n\n";
252 
253  isOk &= isOk_;
254  }
255 
256  if (pOut)
257  *pOut << (isOk ? "--- successfully" : "*** failurefully") << " completed "
258  << "\n\n";
259  if (pOut)
260  pOut->flush();
261 
262  return isOk;
263 }
264 
265 #endif // TESTING
float ten(float) const
Definition: powfast.cpp:112
unsigned int precision_m
Definition: powfast.hpp:57
unsigned int * pTable_m
Definition: powfast.hpp:58
float e(float) const
Definition: powfast.cpp:106
unsigned int precision() const
Definition: powfast.cpp:124
const PowFast & POWFAST()
default instance --------------------------------------------------------—
Definition: powfast.cpp:131
float r(float logr, float f) const
Definition: powfast.cpp:118
~PowFast()
Definition: powfast.cpp:94
int ilog2(int value_inp)
Definition: log2.cpp:9
float two(float) const
Definition: powfast.cpp:100
PowFast(unsigned int precision=11)
wrapper class -----------------------------------------------------------—
Definition: powfast.cpp:87