Eigen  3.2.7
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 <[email protected]>
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 
13 namespace Eigen {
14 
15 namespace internal {
16 
37 template<typename T, typename dummy = void>
38 struct global_math_functions_filtering_base
39 {
40  typedef T type;
41 };
42 
43 template<typename T> struct always_void { typedef void type; };
44 
45 template<typename T>
46 struct 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 
61 template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
62 struct 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 
71 template<typename Scalar>
72 struct 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 
82 template<typename Scalar> struct real_impl : real_default_impl<Scalar> {};
83 
84 template<typename Scalar>
85 struct real_retval
86 {
87  typedef typename NumTraits<Scalar>::Real type;
88 };
89 
90 
91 /****************************************************************************
92 * Implementation of imag *
93 ****************************************************************************/
94 
95 template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
96 struct 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 
105 template<typename Scalar>
106 struct 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 
116 template<typename Scalar> struct imag_impl : imag_default_impl<Scalar> {};
117 
118 template<typename Scalar>
119 struct imag_retval
120 {
121  typedef typename NumTraits<Scalar>::Real type;
122 };
123 
124 /****************************************************************************
125 * Implementation of real_ref *
126 ****************************************************************************/
127 
128 template<typename Scalar>
129 struct 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 
142 template<typename Scalar>
143 struct real_ref_retval
144 {
145  typedef typename NumTraits<Scalar>::Real & type;
146 };
147 
148 /****************************************************************************
149 * Implementation of imag_ref *
150 ****************************************************************************/
151 
152 template<typename Scalar, bool IsComplex>
153 struct 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 
166 template<typename Scalar>
167 struct 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 
179 template<typename Scalar>
180 struct imag_ref_impl : imag_ref_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
181 
182 template<typename Scalar>
183 struct imag_ref_retval
184 {
185  typedef typename NumTraits<Scalar>::Real & type;
186 };
187 
188 /****************************************************************************
189 * Implementation of conj *
190 ****************************************************************************/
191 
192 template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex>
193 struct conj_impl
194 {
195  static inline Scalar run(const Scalar& x)
196  {
197  return x;
198  }
199 };
200 
201 template<typename Scalar>
202 struct 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 
211 template<typename Scalar>
212 struct conj_retval
213 {
214  typedef Scalar type;
215 };
216 
217 /****************************************************************************
218 * Implementation of abs2 *
219 ****************************************************************************/
220 
221 template<typename Scalar>
222 struct abs2_impl
223 {
224  typedef typename NumTraits<Scalar>::Real RealScalar;
225  static inline RealScalar run(const Scalar& x)
226  {
227  return x*x;
228  }
229 };
230 
231 template<typename RealScalar>
232 struct abs2_impl<std::complex<RealScalar> >
233 {
234  static inline RealScalar run(const std::complex<RealScalar>& x)
235  {
236  return real(x)*real(x) + imag(x)*imag(x);
237  }
238 };
239 
240 template<typename Scalar>
241 struct abs2_retval
242 {
243  typedef typename NumTraits<Scalar>::Real type;
244 };
245 
246 /****************************************************************************
247 * Implementation of norm1 *
248 ****************************************************************************/
249 
250 template<typename Scalar, bool IsComplex>
251 struct norm1_default_impl
252 {
253  typedef typename NumTraits<Scalar>::Real RealScalar;
254  static inline RealScalar run(const Scalar& x)
255  {
256  using std::abs;
257  return abs(real(x)) + abs(imag(x));
258  }
259 };
260 
261 template<typename Scalar>
262 struct norm1_default_impl<Scalar, false>
263 {
264  static inline Scalar run(const Scalar& x)
265  {
266  using std::abs;
267  return abs(x);
268  }
269 };
270 
271 template<typename Scalar>
272 struct norm1_impl : norm1_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
273 
274 template<typename Scalar>
275 struct norm1_retval
276 {
277  typedef typename NumTraits<Scalar>::Real type;
278 };
279 
280 /****************************************************************************
281 * Implementation of hypot *
282 ****************************************************************************/
283 
284 template<typename Scalar>
285 struct hypot_impl
286 {
287  typedef typename NumTraits<Scalar>::Real RealScalar;
288  static inline RealScalar run(const Scalar& x, const Scalar& y)
289  {
290  using std::max;
291  using std::min;
292  using std::abs;
293  using std::sqrt;
294  RealScalar _x = abs(x);
295  RealScalar _y = abs(y);
296  RealScalar p = (max)(_x, _y);
297  if(p==RealScalar(0)) return RealScalar(0);
298  RealScalar q = (min)(_x, _y);
299  RealScalar qp = q/p;
300  return p * sqrt(RealScalar(1) + qp*qp);
301  }
302 };
303 
304 template<typename Scalar>
305 struct hypot_retval
306 {
307  typedef typename NumTraits<Scalar>::Real type;
308 };
309 
310 /****************************************************************************
311 * Implementation of cast *
312 ****************************************************************************/
313 
314 template<typename OldType, typename NewType>
315 struct cast_impl
316 {
317  static inline NewType run(const OldType& x)
318  {
319  return static_cast<NewType>(x);
320  }
321 };
322 
323 // here, for once, we're plainly returning NewType: we don't want cast to do weird things.
324 
325 template<typename OldType, typename NewType>
326 inline NewType cast(const OldType& x)
327 {
328  return cast_impl<OldType, NewType>::run(x);
329 }
330 
331 /****************************************************************************
332 * Implementation of atanh2 *
333 ****************************************************************************/
334 
335 template<typename Scalar, bool IsInteger>
336 struct atanh2_default_impl
337 {
338  typedef Scalar retval;
339  typedef typename NumTraits<Scalar>::Real RealScalar;
340  static inline Scalar run(const Scalar& x, const Scalar& y)
341  {
342  using std::abs;
343  using std::log;
344  using std::sqrt;
345  Scalar z = x / y;
346  if (y == Scalar(0) || abs(z) > sqrt(NumTraits<RealScalar>::epsilon()))
347  return RealScalar(0.5) * log((y + x) / (y - x));
348  else
349  return z + z*z*z / RealScalar(3);
350  }
351 };
352 
353 template<typename Scalar>
354 struct atanh2_default_impl<Scalar, true>
355 {
356  static inline Scalar run(const Scalar&, const Scalar&)
357  {
358  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
359  return Scalar(0);
360  }
361 };
362 
363 template<typename Scalar>
364 struct atanh2_impl : atanh2_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
365 
366 template<typename Scalar>
367 struct atanh2_retval
368 {
369  typedef Scalar type;
370 };
371 
372 /****************************************************************************
373 * Implementation of pow *
374 ****************************************************************************/
375 
376 template<typename Scalar, bool IsInteger>
377 struct pow_default_impl
378 {
379  typedef Scalar retval;
380  static inline Scalar run(const Scalar& x, const Scalar& y)
381  {
382  using std::pow;
383  return pow(x, y);
384  }
385 };
386 
387 template<typename Scalar>
388 struct pow_default_impl<Scalar, true>
389 {
390  static inline Scalar run(Scalar x, Scalar y)
391  {
392  Scalar res(1);
393  eigen_assert(!NumTraits<Scalar>::IsSigned || y >= 0);
394  if(y & 1) res *= x;
395  y >>= 1;
396  while(y)
397  {
398  x *= x;
399  if(y&1) res *= x;
400  y >>= 1;
401  }
402  return res;
403  }
404 };
405 
406 template<typename Scalar>
407 struct pow_impl : pow_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
408 
409 template<typename Scalar>
410 struct pow_retval
411 {
412  typedef Scalar type;
413 };
414 
415 /****************************************************************************
416 * Implementation of random *
417 ****************************************************************************/
418 
419 template<typename Scalar,
420  bool IsComplex,
421  bool IsInteger>
422 struct random_default_impl {};
423 
424 template<typename Scalar>
425 struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
426 
427 template<typename Scalar>
428 struct random_retval
429 {
430  typedef Scalar type;
431 };
432 
433 template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y);
434 template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random();
435 
436 template<typename Scalar>
437 struct random_default_impl<Scalar, false, false>
438 {
439  static inline Scalar run(const Scalar& x, const Scalar& y)
440  {
441  return x + (y-x) * Scalar(std::rand()) / Scalar(RAND_MAX);
442  }
443  static inline Scalar run()
444  {
445  return run(Scalar(NumTraits<Scalar>::IsSigned ? -1 : 0), Scalar(1));
446  }
447 };
448 
449 enum {
450  floor_log2_terminate,
451  floor_log2_move_up,
452  floor_log2_move_down,
453  floor_log2_bogus
454 };
455 
456 template<unsigned int n, int lower, int upper> struct floor_log2_selector
457 {
458  enum { middle = (lower + upper) / 2,
459  value = (upper <= lower + 1) ? int(floor_log2_terminate)
460  : (n < (1 << middle)) ? int(floor_log2_move_down)
461  : (n==0) ? int(floor_log2_bogus)
462  : int(floor_log2_move_up)
463  };
464 };
465 
466 template<unsigned int n,
467  int lower = 0,
468  int upper = sizeof(unsigned int) * CHAR_BIT - 1,
469  int selector = floor_log2_selector<n, lower, upper>::value>
470 struct floor_log2 {};
471 
472 template<unsigned int n, int lower, int upper>
473 struct floor_log2<n, lower, upper, floor_log2_move_down>
474 {
475  enum { value = floor_log2<n, lower, floor_log2_selector<n, lower, upper>::middle>::value };
476 };
477 
478 template<unsigned int n, int lower, int upper>
479 struct floor_log2<n, lower, upper, floor_log2_move_up>
480 {
481  enum { value = floor_log2<n, floor_log2_selector<n, lower, upper>::middle, upper>::value };
482 };
483 
484 template<unsigned int n, int lower, int upper>
485 struct floor_log2<n, lower, upper, floor_log2_terminate>
486 {
487  enum { value = (n >= ((unsigned int)(1) << (lower+1))) ? lower+1 : lower };
488 };
489 
490 template<unsigned int n, int lower, int upper>
491 struct floor_log2<n, lower, upper, floor_log2_bogus>
492 {
493  // no value, error at compile time
494 };
495 
496 template<typename Scalar>
497 struct random_default_impl<Scalar, false, true>
498 {
499  typedef typename NumTraits<Scalar>::NonInteger NonInteger;
500 
501  static inline Scalar run(const Scalar& x, const Scalar& y)
502  {
503  return x + Scalar((NonInteger(y)-x+1) * std::rand() / (RAND_MAX + NonInteger(1)));
504  }
505 
506  static inline Scalar run()
507  {
508 #ifdef EIGEN_MAKING_DOCS
509  return run(Scalar(NumTraits<Scalar>::IsSigned ? -10 : 0), Scalar(10));
510 #else
511  enum { rand_bits = floor_log2<(unsigned int)(RAND_MAX)+1>::value,
512  scalar_bits = sizeof(Scalar) * CHAR_BIT,
513  shift = EIGEN_PLAIN_ENUM_MAX(0, int(rand_bits) - int(scalar_bits)),
514  offset = NumTraits<Scalar>::IsSigned ? (1 << (EIGEN_PLAIN_ENUM_MIN(rand_bits,scalar_bits)-1)) : 0
515  };
516  return Scalar((std::rand() >> shift) - offset);
517 #endif
518  }
519 };
520 
521 template<typename Scalar>
522 struct random_default_impl<Scalar, true, false>
523 {
524  static inline Scalar run(const Scalar& x, const Scalar& y)
525  {
526  return Scalar(random(real(x), real(y)),
527  random(imag(x), imag(y)));
528  }
529  static inline Scalar run()
530  {
531  typedef typename NumTraits<Scalar>::Real RealScalar;
532  return Scalar(random<RealScalar>(), random<RealScalar>());
533  }
534 };
535 
536 template<typename Scalar>
537 inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y)
538 {
539  return EIGEN_MATHFUNC_IMPL(random, Scalar)::run(x, y);
540 }
541 
542 template<typename Scalar>
543 inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random()
544 {
545  return EIGEN_MATHFUNC_IMPL(random, Scalar)::run();
546 }
547 
548 } // end namespace internal
549 
550 /****************************************************************************
551 * Generic math function *
552 ****************************************************************************/
553 
554 namespace numext {
555 
556 template<typename Scalar>
557 inline EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar& x)
558 {
559  return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x);
560 }
561 
562 template<typename Scalar>
563 inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar& x)
564 {
565  return internal::real_ref_impl<Scalar>::run(x);
566 }
567 
568 template<typename Scalar>
569 inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) real_ref(Scalar& x)
570 {
571  return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x);
572 }
573 
574 template<typename Scalar>
575 inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) imag(const Scalar& x)
576 {
577  return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x);
578 }
579 
580 template<typename Scalar>
581 inline typename internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) >::type imag_ref(const Scalar& x)
582 {
583  return internal::imag_ref_impl<Scalar>::run(x);
584 }
585 
586 template<typename Scalar>
587 inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) imag_ref(Scalar& x)
588 {
589  return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x);
590 }
591 
592 template<typename Scalar>
593 inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) conj(const Scalar& x)
594 {
595  return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x);
596 }
597 
598 template<typename Scalar>
599 inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) abs2(const Scalar& x)
600 {
601  return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x);
602 }
603 
604 template<typename Scalar>
605 inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) norm1(const Scalar& x)
606 {
607  return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x);
608 }
609 
610 template<typename Scalar>
611 inline EIGEN_MATHFUNC_RETVAL(hypot, Scalar) hypot(const Scalar& x, const Scalar& y)
612 {
613  return EIGEN_MATHFUNC_IMPL(hypot, Scalar)::run(x, y);
614 }
615 
616 template<typename Scalar>
617 inline EIGEN_MATHFUNC_RETVAL(atanh2, Scalar) atanh2(const Scalar& x, const Scalar& y)
618 {
619  return EIGEN_MATHFUNC_IMPL(atanh2, Scalar)::run(x, y);
620 }
621 
622 template<typename Scalar>
623 inline EIGEN_MATHFUNC_RETVAL(pow, Scalar) pow(const Scalar& x, const Scalar& y)
624 {
625  return EIGEN_MATHFUNC_IMPL(pow, Scalar)::run(x, y);
626 }
627 
628 // std::isfinite is non standard, so let's define our own version,
629 // even though it is not very efficient.
630 template<typename T> bool (isfinite)(const T& x)
631 {
632  return x<NumTraits<T>::highest() && x>NumTraits<T>::lowest();
633 }
634 
635 } // end namespace numext
636 
637 namespace internal {
638 
639 /****************************************************************************
640 * Implementation of fuzzy comparisons *
641 ****************************************************************************/
642 
643 template<typename Scalar,
644  bool IsComplex,
645  bool IsInteger>
646 struct scalar_fuzzy_default_impl {};
647 
648 template<typename Scalar>
649 struct scalar_fuzzy_default_impl<Scalar, false, false>
650 {
651  typedef typename NumTraits<Scalar>::Real RealScalar;
652  template<typename OtherScalar>
653  static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec)
654  {
655  using std::abs;
656  return abs(x) <= abs(y) * prec;
657  }
658  static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
659  {
660  using std::min;
661  using std::abs;
662  return abs(x - y) <= (min)(abs(x), abs(y)) * prec;
663  }
664  static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar& prec)
665  {
666  return x <= y || isApprox(x, y, prec);
667  }
668 };
669 
670 template<typename Scalar>
671 struct scalar_fuzzy_default_impl<Scalar, false, true>
672 {
673  typedef typename NumTraits<Scalar>::Real RealScalar;
674  template<typename OtherScalar>
675  static inline bool isMuchSmallerThan(const Scalar& x, const Scalar&, const RealScalar&)
676  {
677  return x == Scalar(0);
678  }
679  static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar&)
680  {
681  return x == y;
682  }
683  static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar&)
684  {
685  return x <= y;
686  }
687 };
688 
689 template<typename Scalar>
690 struct scalar_fuzzy_default_impl<Scalar, true, false>
691 {
692  typedef typename NumTraits<Scalar>::Real RealScalar;
693  template<typename OtherScalar>
694  static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec)
695  {
696  return numext::abs2(x) <= numext::abs2(y) * prec * prec;
697  }
698  static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
699  {
700  using std::min;
701  return numext::abs2(x - y) <= (min)(numext::abs2(x), numext::abs2(y)) * prec * prec;
702  }
703 };
704 
705 template<typename Scalar>
706 struct scalar_fuzzy_impl : scalar_fuzzy_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
707 
708 template<typename Scalar, typename OtherScalar>
709 inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y,
710  typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision())
711 {
712  return scalar_fuzzy_impl<Scalar>::template isMuchSmallerThan<OtherScalar>(x, y, precision);
713 }
714 
715 template<typename Scalar>
716 inline bool isApprox(const Scalar& x, const Scalar& y,
717  typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision())
718 {
719  return scalar_fuzzy_impl<Scalar>::isApprox(x, y, precision);
720 }
721 
722 template<typename Scalar>
723 inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y,
724  typename NumTraits<Scalar>::Real precision = NumTraits<Scalar>::dummy_precision())
725 {
726  return scalar_fuzzy_impl<Scalar>::isApproxOrLessThan(x, y, precision);
727 }
728 
729 /******************************************
730 *** The special case of the bool type ***
731 ******************************************/
732 
733 template<> struct random_impl<bool>
734 {
735  static inline bool run()
736  {
737  return random<int>(0,1)==0 ? false : true;
738  }
739 };
740 
741 template<> struct scalar_fuzzy_impl<bool>
742 {
743  typedef bool RealScalar;
744 
745  template<typename OtherScalar>
746  static inline bool isMuchSmallerThan(const bool& x, const bool&, const bool&)
747  {
748  return !x;
749  }
750 
751  static inline bool isApprox(bool x, bool y, bool)
752  {
753  return x == y;
754  }
755 
756  static inline bool isApproxOrLessThan(const bool& x, const bool& y, const bool&)
757  {
758  return (!x) || y;
759  }
760 
761 };
762 
763 
764 } // end namespace internal
765 
766 } // end namespace Eigen
767 
768 #endif // EIGEN_MATHFUNCTIONS_H