1 /*
2 ================================================================================
3 PROJECT:
4
5 Eddy C++ Utilities Project
6
7 CONTENTS:
8
9 Implementation of class Math.
10
11 NOTES:
12
13 See notes of Math.hpp.
14
15 PROGRAMMERS:
16
17 John Eddy (jpeddy@sandia.gov) (JE)
18
19 ORGANIZATION:
20
21 Sandia National Laboratories
22
23 COPYRIGHT:
24
25 See the LICENSE file in the top level JEGA directory.
26
27 VERSION:
28
29 1.0.0
30
31 CHANGES:
32
33 Tue Jun 03 10:30:43 2003 - Original Version (JE)
34
35 ================================================================================
36 */
37
38
39
40 /*
41 ================================================================================
42 Document This File
43 ================================================================================
44 */
45 /** \file
46 * \brief Contains the implementation of the Math class.
47 */
48
49
50
51
52 /*
53 ================================================================================
54 Includes
55 ================================================================================
56 */
57 #include <cmath>
58 #include <cfloat>
59 #include <cstddef>
60 #include <cstdlib>
61 #include <cassert>
62 #include "../include/Math.hpp"
63
64
65
66
67
68 /*
69 ================================================================================
70 Namespace Using Directives
71 ================================================================================
72 */
73 using namespace std;
74
75
76
77
78
79
80
81
82 /*
83 ================================================================================
84 Begin Namespace
85 ================================================================================
86 */
87 namespace eddy {
88 namespace utilities {
89
90
91
92
93
94
95
96 /*
97 ================================================================================
98 Static Member Data Definitions
99 ================================================================================
100 */
101
102
103
104
105
106
107
108
109 /*
110 ================================================================================
111 Mutators
112 ================================================================================
113 */
114
115
116
117
118
119
120
121
122 /*
123 ================================================================================
124 Accessors
125 ================================================================================
126 */
127
128
129
130
131
132
133
134
135 /*
136 ================================================================================
137 Public Methods
138 ================================================================================
139 */
140
141 double
Floor(double arg)142 Math::Floor(
143 double arg
144 )
145 {
146 return floor(arg);
147
148 } // Math::Floor
149
150 float
Floor(float arg)151 Math::Floor(
152 float arg
153 )
154 {
155 return static_cast<float>(floor(static_cast<double>(arg)));
156
157 } // Math::Floor
158
159 double
Ceil(double arg)160 Math::Ceil(
161 double arg
162 )
163 {
164 return ceil(arg);
165
166 } // Math::Ceil
167
168 float
Ceil(float arg)169 Math::Ceil(
170 float arg
171 )
172 {
173 return static_cast<float>(ceil(static_cast<double>(arg)));
174
175 } // Math::Ceil
176
177 double
Exp(double arg)178 Math::Exp(
179 double arg
180 )
181 {
182 return exp(arg);
183
184 } // Math::Exp
185
186 float
Exp(float arg)187 Math::Exp(
188 float arg
189 )
190 {
191 return static_cast<float>(exp(static_cast<double>(arg)));
192
193 } // Math::Exp
194
195 double
Round(double arg,int prec)196 Math::Round(
197 double arg,
198 int prec
199 )
200 {
201 double conv = prec == 0.0 ? 1.0 : pow(10.0, static_cast<double>(prec));
202 return (arg > 0.0) ?
203 floor(arg * conv + 0.5) / conv :
204 ceil (arg * conv - 0.5) / conv;
205
206 } // Math::Round
207
208 float
Round(float arg,int prec)209 Math::Round(
210 float arg,
211 int prec
212 )
213 {
214 float conv = prec == 0.0f ? 1.0f : pow(10.0f, static_cast<float>(prec));
215 return (arg > 0.0) ?
216 floor(arg * conv + 0.5f) / conv :
217 ceil (arg * conv - 0.5f) / conv;
218
219 } // Math::Round
220
221 double
Truncate(double arg,int prec)222 Math::Truncate(
223 double arg,
224 int prec
225 )
226 {
227 if(prec == 0) return static_cast<double>(static_cast<int>(arg));
228 double conv = pow(10.0, static_cast<double>(prec));
229 return (arg > 0.0) ?
230 floor(arg * conv) / conv :
231 ceil (arg * conv) / conv;
232
233 } // Math::Truncate
234
235 float
Truncate(float arg,int prec)236 Math::Truncate(
237 float arg,
238 int prec
239 )
240 {
241 if(prec == 0) return static_cast<float>(static_cast<int>(arg));
242 float conv = pow(10.0f, static_cast<float>(prec));
243 return (arg > 0.0) ?
244 floor(arg * conv) / conv :
245 ceil (arg * conv) / conv;
246
247 } // Math::Truncate
248
249 double
Sqrt(double arg)250 Math::Sqrt(
251 double arg
252 )
253 {
254 return sqrt(arg);
255
256 } // Math::Sqrt
257
258 float
Sqrt(float arg)259 Math::Sqrt(
260 float arg
261 )
262 {
263 return static_cast<float>(sqrt(static_cast<double>(arg)));
264
265 } // Math::Sqrt
266
267 double
Log(double arg,double base)268 Math::Log(
269 double arg,
270 double base
271 )
272 {
273 return (base==10.0) ? log10(arg) : (log10(arg) / log10(base));
274
275 } // Math::Log
276
277 float
Log(float arg,float base)278 Math::Log(
279 float arg,
280 float base
281 )
282 {
283 return (base==10.0f) ? Log10(arg) : (Log10(arg) / Log10(base));
284
285 } // Math::Log
286
287 double
Log10(double arg)288 Math::Log10(
289 double arg
290 )
291 {
292 return log10(arg);
293
294 } // Math::Log10
295
296 float
Log10(float arg)297 Math::Log10(
298 float arg
299 )
300 {
301 return static_cast<float>(log10(static_cast<double>(arg)));
302
303 } // Math::Log10
304
305 double
Ln(double arg)306 Math::Ln(
307 double arg
308 )
309 {
310 return log(arg);
311
312 } // Math::Ln
313
314 float
Ln(float arg)315 Math::Ln(
316 float arg
317 )
318 {
319 return static_cast<float>(log(static_cast<double>(arg)));
320
321 } // Math::Ln
322
323 double
Pow(double arg,double exp)324 Math::Pow(
325 double arg,
326 double exp
327 )
328 {
329 return pow(arg, exp);
330
331 } // Math::Pow
332
333 float
Pow(float arg,float exp)334 Math::Pow(
335 float arg,
336 float exp
337 )
338 {
339 return static_cast<float>(
340 pow(static_cast<double>(arg), static_cast<double>(exp))
341 );
342
343 } // Math::Pow
344
345 double
Abs(double arg)346 Math::Abs(
347 double arg
348 )
349 {
350 return fabs(arg);
351
352 } // Math::Abs
353
354 int
Abs(int arg)355 Math::Abs(
356 int arg
357 )
358 {
359 return abs(arg);
360
361 } // Math::Abs
362
363 long
Abs(long arg)364 Math::Abs(
365 long arg
366 )
367 {
368 return labs(arg);
369
370 } // Math::Abs
371
372 float
Abs(float arg)373 Math::Abs(
374 float arg
375 )
376 {
377 return static_cast<float>(fabs(static_cast<double>(arg)));
378
379 } // Math::Abs
380
381 double
Sin(double arg)382 Math::Sin(
383 double arg
384 )
385 {
386 return sin(arg);
387
388 } // Math::Sin
389
390 float
Sin(float arg)391 Math::Sin(
392 float arg
393 )
394 {
395 return static_cast<float>(sin(static_cast<double>(arg)));
396
397 } // Math::Sin
398
399 double
Cos(double arg)400 Math::Cos(
401 double arg
402 )
403 {
404 return cos(arg);
405
406 } // Math::Cos
407
408 float
Cos(float arg)409 Math::Cos(
410 float arg
411 )
412 {
413 return static_cast<float>(cos(static_cast<double>(arg)));
414
415 } // Math::Cos
416
417 double
Tan(double arg)418 Math::Tan(
419 double arg
420 )
421 {
422 return tan(arg);
423
424 } // Math::Tan
425
426 float
Tan(float arg)427 Math::Tan(
428 float arg
429 )
430 {
431 return static_cast<float>(tan(static_cast<double>(arg)));
432
433 } // Math::Tan
434
435 double
ArcSin(double arg)436 Math::ArcSin(
437 double arg
438 )
439 {
440 return asin(arg);
441
442 } // Math::ArcSin
443
444 float
ArcSin(float arg)445 Math::ArcSin(
446 float arg
447 )
448 {
449 return static_cast<float>(asin(static_cast<double>(arg)));
450
451 } // Math::ArcSin
452
453 double
ArcCos(double arg)454 Math::ArcCos(
455 double arg
456 )
457 {
458 return acos(arg);
459
460 } // Math::ArcCos
461
462 float
ArcCos(float arg)463 Math::ArcCos(
464 float arg
465 )
466 {
467 return static_cast<float>(acos(static_cast<double>(arg)));
468
469 } // Math::ArcCos
470
471 double
ArcTan(double arg)472 Math::ArcTan(
473 double arg
474 )
475 {
476 return atan(arg);
477
478 } // Math::ArcTan
479
480 float
ArcTan(float arg)481 Math::ArcTan(
482 float arg
483 )
484 {
485 return static_cast<float>(atan(static_cast<double>(arg)));
486
487 } // Math::ArcTan
488
489 double
LambertW(const double arg)490 Math::LambertW(
491 const double arg
492 )
493 {
494 if(0.0 == arg) return 0.0;
495
496 const double eps = 4.0e-16;
497 const double em1 = 0.3678794411714423215955237701614608;
498
499 assert(arg>=-em1);
500
501 double p,e,t,w;
502
503 if(arg < (-em1+1e-4))
504 {
505 // series near -em1 in sqrt(q)
506 double q = arg+em1;
507 double r=Sqrt(q);
508 double q2=q*q;
509 double q3=q2*q;
510
511 return
512 -1.0
513 +2.331643981597124203363536062168*r
514 -1.812187885639363490240191647568*q
515 +1.936631114492359755363277457668*r*q
516 -2.353551201881614516821543561516*q2
517 +3.066858901050631912893148922704*r*q2
518 -4.175335600258177138854984177460*q3
519 +5.858023729874774148815053846119*r*q3
520 -8.401032217523977370984161688514*q3*q; // error approx 1e-16
521 }
522
523 /* initial approx for iteration... */
524 if(arg < 1.0)
525 {
526 /* series near 0 */
527 p = Sqrt(2.0 * (2.7182818284590452353602874713526625*arg + 1.0));
528 w = -1.0 + p * (1.0 + p *
529 (-0.333333333333333333333 + p * 0.152777777777777777777777)
530 );
531 }
532 else w=Ln(arg); /* asymptotic */
533
534 if(arg > 3.0) w -= Ln(w); /* useful? */
535
536 for (unsigned int i=0; i<10; ++i)
537 {
538 /* Halley iteration */
539 e = Exp(w);
540 t = w * e - arg;
541 p = w + 1.0;
542 t /= e * p - 0.5 * (p + 1.0) * t / p;
543 w -= t;
544 if (Abs(t) < eps * (1.0 + Abs(w))) return w; /* rel-abs error */
545 }
546
547 assert(false && "Lambert W - No convergence.");
548 return 0.0f;
549
550 } // Math::LambertW
551
552 /*
553 ================================================================================
554 Subclass Visible Methods
555 ================================================================================
556 */
557
558
559
560
561
562
563
564
565 /*
566 ================================================================================
567 Subclass Overridable Methods
568 ================================================================================
569 */
570
571
572
573
574
575
576
577
578 /*
579 ================================================================================
580 Private Methods
581 ================================================================================
582 */
583
584
585
586
587
588
589
590
591 /*
592 ================================================================================
593 Structors
594 ================================================================================
595 */
596
597
598
599
600
601
602
603
604 /*
605 ================================================================================
606 End Namespace
607 ================================================================================
608 */
609 } // namespace utilities
610 } // namespace eddy
611