Eigen  3.2.10
 
Loading...
Searching...
No Matches
MathFunctions.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2006-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_MATHFUNCTIONS_H
11#define EIGEN_MATHFUNCTIONS_H
12
13namespace Eigen {
14
15namespace internal {
16
36
37template<typename T, typename dummy = void>
38struct global_math_functions_filtering_base
39{
40 typedef T type;
41};
42
43template<typename T> struct always_void { typedef void type; };
44
45template<typename T>
46struct global_math_functions_filtering_base
47 <T,
48 typename always_void<typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl>::type
49 >
50{
51 typedef typename T::Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl type;
52};
53
54#define EIGEN_MATHFUNC_IMPL(func, scalar) Eigen::internal::func##_impl<typename Eigen::internal::global_math_functions_filtering_base<scalar>::type>
55#define EIGEN_MATHFUNC_RETVAL(func, scalar) typename Eigen::internal::func##_retval<typename Eigen::internal::global_math_functions_filtering_base<scalar>::type>::type
56
57/****************************************************************************
58* Implementation of real *
59****************************************************************************/
60
61template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
62struct real_default_impl
63{
64 typedef typename NumTraits<Scalar>::Real RealScalar;
65 static inline RealScalar run(const Scalar& x)
66 {
67 return x;
68 }
69};
70
71template<typename Scalar>
72struct real_default_impl<Scalar,true>
73{
74 typedef typename NumTraits<Scalar>::Real RealScalar;
75 static inline RealScalar run(const Scalar& x)
76 {
77 using std::real;
78 return real(x);
79 }
80};
81
82template<typename Scalar> struct real_impl : real_default_impl<Scalar> {};
83
84template<typename Scalar>
85struct real_retval
86{
87 typedef typename NumTraits<Scalar>::Real type;
88};
89
90
91/****************************************************************************
92* Implementation of imag *
93****************************************************************************/
94
95template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
96struct imag_default_impl
97{
98 typedef typename NumTraits<Scalar>::Real RealScalar;
99 static inline RealScalar run(const Scalar&)
100 {
101 return RealScalar(0);
102 }
103};
104
105template<typename Scalar>
106struct imag_default_impl<Scalar,true>
107{
108 typedef typename NumTraits<Scalar>::Real RealScalar;
109 static inline RealScalar run(const Scalar& x)
110 {
111 using std::imag;
112 return imag(x);
113 }
114};
115
116template<typename Scalar> struct imag_impl : imag_default_impl<Scalar> {};
117
118template<typename Scalar>
119struct imag_retval
120{
121 typedef typename NumTraits<Scalar>::Real type;
122};
123
124/****************************************************************************
125* Implementation of real_ref *
126****************************************************************************/
127
128template<typename Scalar>
129struct real_ref_impl
130{
131 typedef typename NumTraits<Scalar>::Real RealScalar;
132 static inline RealScalar& run(Scalar& x)
133 {
134 return reinterpret_cast<RealScalar*>(&x)[0];
135 }
136 static inline const RealScalar& run(const Scalar& x)
137 {
138 return reinterpret_cast<const RealScalar*>(&x)[0];
139 }
140};
141
142template<typename Scalar>
143struct real_ref_retval
144{
145 typedef typename NumTraits<Scalar>::Real & type;
146};
147
148/****************************************************************************
149* Implementation of imag_ref *
150****************************************************************************/
151
152template<typename Scalar, bool IsComplex>
153struct imag_ref_default_impl
154{
155 typedef typename NumTraits<Scalar>::Real RealScalar;
156 static inline RealScalar& run(Scalar& x)
157 {
158 return reinterpret_cast<RealScalar*>(&x)[1];
159 }
160 static inline const RealScalar& run(const Scalar& x)
161 {
162 return reinterpret_cast<RealScalar*>(&x)[1];
163 }
164};
165
166template<typename Scalar>
167struct imag_ref_default_impl<Scalar, false>
168{
169 static inline Scalar run(Scalar&)
170 {
171 return Scalar(0);
172 }
173 static inline const Scalar run(const Scalar&)
174 {
175 return Scalar(0);
176 }
177};
178
179template<typename Scalar>
180struct imag_ref_impl : imag_ref_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
181
182template<typename Scalar>
183struct imag_ref_retval
184{
185 typedef typename NumTraits<Scalar>::Real & type;
186};
187
188/****************************************************************************
189* Implementation of conj *
190****************************************************************************/
191
192template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
193struct conj_impl
194{
195 static inline Scalar run(const Scalar& x)
196 {
197 return x;
198 }
199};
200
201template<typename Scalar>
202struct conj_impl<Scalar,true>
203{
204 static inline Scalar run(const Scalar& x)
205 {
206 using std::conj;
207 return conj(x);
208 }
209};
210
211template<typename Scalar>
212struct conj_retval
213{
214 typedef Scalar type;
215};
216
217/****************************************************************************
218* Implementation of abs2 *
219****************************************************************************/
220
221template<typename Scalar,bool IsComplex>
222struct abs2_impl_default
223{
224 typedef typename NumTraits<Scalar>::Real RealScalar;
225 static inline RealScalar run(const Scalar& x)
226 {
227 return x*x;
228 }
229};
230
231template<typename Scalar>
232struct abs2_impl_default<Scalar, true> // IsComplex
233{
234 typedef typename NumTraits<Scalar>::Real RealScalar;
235 static inline RealScalar run(const Scalar& x)
236 {
237 return real(x)*real(x) + imag(x)*imag(x);
238 }
239};
240
241template<typename Scalar>
242struct abs2_impl
243{
244 typedef typename NumTraits<Scalar>::Real RealScalar;
245 static inline RealScalar run(const Scalar& x)
246 {
247 return abs2_impl_default<Scalar,NumTraits<Scalar>::IsComplex>::run(x);
248 }
249};
250
251template<typename Scalar>
252struct abs2_retval
253{
254 typedef typename NumTraits<Scalar>::Real type;
255};
256
257/****************************************************************************
258* Implementation of norm1 *
259****************************************************************************/
260
261template<typename Scalar, bool IsComplex>
262struct norm1_default_impl
263{
264 typedef typename NumTraits<Scalar>::Real RealScalar;
265 static inline RealScalar run(const Scalar& x)
266 {
267 using std::abs;
268 return abs(real(x)) + abs(imag(x));
269 }
270};
271
272template<typename Scalar>
273struct norm1_default_impl<Scalar, false>
274{
275 static inline Scalar run(const Scalar& x)
276 {
277 using std::abs;
278 return abs(x);
279 }
280};
281
282template<typename Scalar>
283struct norm1_impl : norm1_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
284
285template<typename Scalar>
286struct norm1_retval
287{
288 typedef typename NumTraits<Scalar>::Real type;
289};
290
291/****************************************************************************
292* Implementation of hypot *
293****************************************************************************/
294
295template<typename Scalar>
296struct hypot_impl
297{
298 typedef typename NumTraits<Scalar>::Real RealScalar;
299 static inline RealScalar run(const Scalar& x, const Scalar& y)
300 {
301 using std::max;
302 using std::min;
303 using std::abs;
304 using std::sqrt;
305 RealScalar _x = abs(x);
306 RealScalar _y = abs(y);
307 RealScalar p = (max)(_x, _y);
308 if(p==RealScalar(0)) return RealScalar(0);
309 RealScalar q = (min)(_x, _y);
310 RealScalar qp = q/p;
311 return p * sqrt(RealScalar(1) + qp*qp);
312 }
313};
314
315template<typename Scalar>
316struct hypot_retval
317{
318 typedef typename NumTraits<Scalar>::Real type;
319};
320
321/****************************************************************************
322* Implementation of cast *
323****************************************************************************/
324
325template<typename OldType, typename NewType>
326struct cast_impl
327{
328 static inline NewType run(const OldType& x)
329 {
330 return static_cast<NewType>(x);
331 }
332};
333
334// here, for once, we're plainly returning NewType: we don't want cast to do weird things.
335
336template<typename OldType, typename NewType>
337inline NewType cast(const OldType& x)
338{
339 return cast_impl<OldType, NewType>::run(x);
340}
341
342/****************************************************************************
343* Implementation of atanh2 *
344****************************************************************************/
345
346template<typename Scalar, bool IsInteger>
347struct atanh2_default_impl
348{
349 typedef Scalar retval;
350 typedef typename NumTraits<Scalar>::Real RealScalar;
351 static inline Scalar run(const Scalar& x, const Scalar& y)
352 {
353 using std::abs;
354 using std::log;
355 using std::sqrt;
356 Scalar z = x / y;
357 if (y == Scalar(0) || abs(z) > sqrt(NumTraits<RealScalar>::epsilon()))
358 return RealScalar(0.5) * log((y + x) / (y - x));
359 else
360 return z + z*z*z / RealScalar(3);
361 }
362};
363
364template<typename Scalar>
365struct atanh2_default_impl<Scalar, true>
366{
367 static inline Scalar run(const Scalar&, const Scalar&)
368 {
369 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
370 return Scalar(0);
371 }
372};
373
374template<typename Scalar>
375struct atanh2_impl : atanh2_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
376
377template<typename Scalar>
378struct atanh2_retval
379{
380 typedef Scalar type;
381};
382
383/****************************************************************************
384* Implementation of pow *
385****************************************************************************/
386
387template<typename Scalar, bool IsInteger>
388struct pow_default_impl
389{
390 typedef Scalar retval;
391 static inline Scalar run(const Scalar& x, const Scalar& y)
392 {
393 using std::pow;
394 return pow(x, y);
395 }
396};
397
398template<typename Scalar>
399struct pow_default_impl<Scalar, true>
400{
401 static inline Scalar run(Scalar x, Scalar y)
402 {
403 Scalar res(1);
404 eigen_assert(!NumTraits<Scalar>::IsSigned || y >= 0);
405 if(y & 1) res *= x;
406 y >>= 1;
407 while(y)
408 {
409 x *= x;
410 if(y&1) res *= x;
411 y >>= 1;
412 }
413 return res;
414 }
415};
416
417template<typename Scalar>
418struct pow_impl : pow_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
419
420template<typename Scalar>
421struct pow_retval
422{
423 typedef Scalar type;
424};
425
426/****************************************************************************
427* Implementation of random *
428****************************************************************************/
429
430template<typename Scalar,
431 bool IsComplex,
432 bool IsInteger>
433struct random_default_impl {};
434
435template<typename Scalar>
436struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
437
438template<typename Scalar>
439struct random_retval
440{
441 typedef Scalar type;
442};
443
444template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y);
445template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random();
446
447template<typename Scalar>
448struct random_default_impl<Scalar, false, false>
449{
450 static inline Scalar run(const Scalar& x, const Scalar& y)
451 {
452 return x + (y-x) * Scalar(std::rand()) / Scalar(RAND_MAX);
453 }
454 static inline Scalar run()
455 {
456 return run(Scalar(NumTraits<Scalar>::IsSigned ? -1 : 0), Scalar(1));
457 }
458};
459
460enum {
461 floor_log2_terminate,
462 floor_log2_move_up,
463 floor_log2_move_down,
464 floor_log2_bogus
465};
466
467template<unsigned int n, int lower, int upper> struct floor_log2_selector
468{
469 enum { middle = (lower + upper) / 2,
470 value = (upper <= lower + 1) ? int(floor_log2_terminate)
471 : (n < (1 << middle)) ? int(floor_log2_move_down)
472 : (n==0) ? int(floor_log2_bogus)
473 : int(floor_log2_move_up)
474 };
475};
476
477template<unsigned int n,
478 int lower = 0,
479 int upper = sizeof(unsigned int) * CHAR_BIT - 1,
480 int selector = floor_log2_selector<n, lower, upper>::value>
481struct floor_log2 {};
482
483template<unsigned int n, int lower, int upper>
484struct floor_log2<n, lower, upper, floor_log2_move_down>
485{
486 enum { value = floor_log2<n, lower, floor_log2_selector<n, lower, upper>::middle>::value };
487};
488
489template<unsigned int n, int lower, int upper>
490struct floor_log2<n, lower, upper, floor_log2_move_up>
491{
492 enum { value = floor_log2<n, floor_log2_selector<n, lower, upper>::middle, upper>::value };
493};
494
495template<unsigned int n, int lower, int upper>
496struct floor_log2<n, lower, upper, floor_log2_terminate>
497{
498 enum { value = (n >= ((unsigned int)(1) << (lower+1))) ? lower+1 : lower };
499};
500
501template<unsigned int n, int lower, int upper>
502struct floor_log2<n, lower, upper, floor_log2_bogus>
503{
504 // no value, error at compile time
505};
506
507template<typename Scalar>
508struct random_default_impl<Scalar, false, true>
509{
510 static inline Scalar run(const Scalar& x, const Scalar& y)
511 {
512 typedef typename conditional<NumTraits<Scalar>::IsSigned,std::ptrdiff_t,std::size_t>::type ScalarX;
513 if(y<x)
514 return x;
515 // the following difference might overflow on a 32 bits system,
516 // but since y>=x the result converted to an unsigned long is still correct.
517 std::size_t range = ScalarX(y)-ScalarX(x);
518 std::size_t offset = 0;
519 // rejection sampling
520 std::size_t divisor = 1;
521 std::size_t multiplier = 1;
522 if(range<RAND_MAX) divisor = (std::size_t(RAND_MAX)+1)/(range+1);
523 else multiplier = 1 + range/(std::size_t(RAND_MAX)+1);
524 do {
525 offset = (std::size_t(std::rand()) * multiplier) / divisor;
526 } while (offset > range);
527 return Scalar(ScalarX(x) + offset);
528 }
529
530 static inline Scalar run()
531 {
532#ifdef EIGEN_MAKING_DOCS
533 return run(Scalar(NumTraits<Scalar>::IsSigned ? -10 : 0), Scalar(10));
534#else
535 enum { rand_bits = floor_log2<(unsigned int)(RAND_MAX)+1>::value,
536 scalar_bits = sizeof(Scalar) * CHAR_BIT,
537 shift = EIGEN_PLAIN_ENUM_MAX(0, int(rand_bits) - int(scalar_bits)),
538 offset = NumTraits<Scalar>::IsSigned ? (1 << (EIGEN_PLAIN_ENUM_MIN(rand_bits,scalar_bits)-1)) : 0
539 };
540 return Scalar((std::rand() >> shift) - offset);
541#endif
542 }
543};
544
545template<typename Scalar>
546struct random_default_impl<Scalar, true, false>
547{
548 static inline Scalar run(const Scalar& x, const Scalar& y)
549 {
550 return Scalar(random(real(x), real(y)),
551 random(imag(x), imag(y)));
552 }
553 static inline Scalar run()
554 {
555 typedef typename NumTraits<Scalar>::Real RealScalar;
556 return Scalar(random<RealScalar>(), random<RealScalar>());
557 }
558};
559
560template<typename Scalar>
561inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y)
562{
563 return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y);
564}
565
566template<typename Scalar>
567inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random()
568{
569 return EIGEN_MATHFUNC_IMPL(random, Scalar)::run();
570}
571
572} // end namespace internal
573
574/****************************************************************************
575* Generic math function *
576****************************************************************************/
577
578namespace numext {
579
580template<typename Scalar>
581inline EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar& x)
582{
583 return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x);
584}
585
586template<typename Scalar>
587inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar& x)
588{
589 return internal::real_ref_impl<Scalar>::run(x);
590}
591
592template<typename Scalar>
593inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) real_ref(Scalar& x)
594{
595 return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x);
596}
597
598template<typename Scalar>
599inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) imag(const Scalar& x)
600{
601 return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x);
602}
603
604template<typename Scalar>
605inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) >::type imag_ref(const Scalar& x)
606{
607 return internal::imag_ref_impl<Scalar>::run(x);
608}
609
610template<typename Scalar>
611inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) imag_ref(Scalar& x)
612{
613 return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x);
614}
615
616template<typename Scalar>
617inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) conj(const Scalar& x)
618{
619 return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x);
620}
621
622template<typename Scalar>
623inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) abs2(const Scalar& x)
624{
625 return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x);
626}
627
628template<typename Scalar>
629inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) norm1(const Scalar& x)
630{
631 return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x);
632}
633
634template<typename Scalar>
635inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) hypot(const Scalar& x, const Scalar& y)
636{
637 return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y);
638}
639
640template<typename Scalar>
641inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y)
642{
643 return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y);
644}
645
646template<typename Scalar>
647inline EIGEN_MATHFUNC_RETVAL(pow, Scalar) pow(const Scalar& x, const Scalar& y)
648{
649 return EIGEN_MATHFUNC_IMPL(pow, Scalar)::run(x, y);
650}
651
652// std::isfinite is non standard, so let's define our own version,
653// even though it is not very efficient.
654template<typename T> bool (isfinite)(const T& x)
655{
656 return x<NumTraits<T>::highest() && x>NumTraits<T>::lowest();
657}
658
659} // end namespace numext
660
661namespace internal {
662
663/****************************************************************************
664* Implementation of fuzzy comparisons *
665****************************************************************************/
666
667template<typename Scalar,
668 bool IsComplex,
669 bool IsInteger>
670struct scalar_fuzzy_default_impl {};
671
672template<typename Scalar>
673struct scalar_fuzzy_default_impl<Scalar, false, false>
674{
675 typedef typename NumTraits<Scalar>::Real RealScalar;
676 template<typename OtherScalar>
677 static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec)
678 {
679 using std::abs;
680 return abs(x) <= abs(y) * prec;
681 }
682 static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
683 {
684 using std::min;
685 using std::abs;
686 return abs(x - y) <= (min)(abs(x), abs(y)) * prec;
687 }
688 static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar& prec)
689 {
690 return x <= y || isApprox(x, y, prec);
691 }
692};
693
694template<typename Scalar>
695struct scalar_fuzzy_default_impl<Scalar, false, true>
696{
697 typedef typename NumTraits<Scalar>::Real RealScalar;
698 template<typename OtherScalar>
699 static inline bool isMuchSmallerThan(const Scalar& x, const Scalar&, const RealScalar&)
700 {
701 return x == Scalar(0);
702 }
703 static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar&)
704 {
705 return x == y;
706 }
707 static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar&)
708 {
709 return x <= y;
710 }
711};
712
713template<typename Scalar>
714struct scalar_fuzzy_default_impl<Scalar, true, false>
715{
716 typedef typename NumTraits<Scalar>::Real RealScalar;
717 template<typename OtherScalar>
718 static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec)
719 {
720 return numext::abs2(x) <= numext::abs2(y) * prec * prec;
721 }
722 static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
723 {
724 using std::min;
725 return numext::abs2(x - y) <= (min)(numext::abs2(x), numext::abs2(y)) * prec * prec;
726 }
727};
728
729template<typename Scalar>
730struct scalar_fuzzy_impl : scalar_fuzzy_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
731
732template<typename Scalar, typename OtherScalar>
733inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y,
734 const typename NumTraits<Scalar>::Real &precision = NumTraits<Scalar>::dummy_precision())
735{
736 return scalar_fuzzy_impl<Scalar>::template isMuchSmallerThan<OtherScalar>(x, y, precision);
737}
738
739template<typename Scalar>
740inline bool isApprox(const Scalar& x, const Scalar& y,
741 const typename NumTraits<Scalar>::Real &precision = NumTraits<Scalar>::dummy_precision())
742{
743 return scalar_fuzzy_impl<Scalar>::isApprox(x, y, precision);
744}
745
746template<typename Scalar>
747inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y,
748 const typename NumTraits<Scalar>::Real &precision = NumTraits<Scalar>::dummy_precision())
749{
750 return scalar_fuzzy_impl<Scalar>::isApproxOrLessThan(x, y, precision);
751}
752
753/******************************************
754*** The special case of the bool type ***
755******************************************/
756
757template<> struct random_impl<bool>
758{
759 static inline bool run()
760 {
761 return random<int>(0,1)==0 ? false : true;
762 }
763};
764
765template<> struct scalar_fuzzy_impl<bool>
766{
767 typedef bool RealScalar;
768
769 template<typename OtherScalar>
770 static inline bool isMuchSmallerThan(const bool& x, const bool&, const bool&)
771 {
772 return !x;
773 }
774
775 static inline bool isApprox(bool x, bool y, bool)
776 {
777 return x == y;
778 }
779
780 static inline bool isApproxOrLessThan(const bool& x, const bool& y, const bool&)
781 {
782 return (!x) || y;
783 }
784
785};
786
787
788} // end namespace internal
789
790} // end namespace Eigen
791
792#endif // EIGEN_MATHFUNCTIONS_H
Eigen::internal
Definition LDLT.h:18