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
17/** \internal \struct global_math_functions_filtering_base
18  *
19  * What it does:
20  * Defines a typedef 'type' as follows:
21  * - if type T has a member typedef Eigen_BaseClassForSpecializationOfGlobalMathFuncImpl, then
22  *   global_math_functions_filtering_base<T>::type is a typedef for it.
23  * - otherwise, global_math_functions_filtering_base<T>::type is a typedef for T.
24  *
25  * How it's used:
26  * To allow to defined the global math functions (like sin...) in certain cases, like the Array expressions.
27  * When you do sin(array1+array2), the object array1+array2 has a complicated expression type, all what you want to know
28  * is that it inherits ArrayBase. So we implement a partial specialization of sin_impl for ArrayBase<Derived>.
29  * So we must make sure to use sin_impl<ArrayBase<Derived> > and not sin_impl<Derived>, otherwise our partial specialization
30  * won't be used. How does sin know that? That's exactly what global_math_functions_filtering_base tells it.
31  *
32  * How it's implemented:
33  * SFINAE in the style of enable_if. Highly susceptible of breaking compilers. With GCC, it sure does work, but if you replace
34  * the typename dummy by an integer template parameter, it doesn't work anymore!
35  */
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>
222struct 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
231template<typename RealScalar>
232struct 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
240template<typename Scalar>
241struct abs2_retval
242{
243  typedef typename NumTraits<Scalar>::Real type;
244};
245
246/****************************************************************************
247* Implementation of norm1                                                *
248****************************************************************************/
249
250template<typename Scalar, bool IsComplex>
251struct 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
261template<typename Scalar>
262struct 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
271template<typename Scalar>
272struct norm1_impl : norm1_default_impl<Scalar, NumTraits<Scalar>::IsComplex> {};
273
274template<typename Scalar>
275struct norm1_retval
276{
277  typedef typename NumTraits<Scalar>::Real type;
278};
279
280/****************************************************************************
281* Implementation of hypot                                                *
282****************************************************************************/
283
284template<typename Scalar>
285struct 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 0;
298    RealScalar q = (min)(_x, _y);
299    RealScalar qp = q/p;
300    return p * sqrt(RealScalar(1) + qp*qp);
301  }
302};
303
304template<typename Scalar>
305struct hypot_retval
306{
307  typedef typename NumTraits<Scalar>::Real type;
308};
309
310/****************************************************************************
311* Implementation of cast                                                 *
312****************************************************************************/
313
314template<typename OldType, typename NewType>
315struct 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
325template<typename OldType, typename NewType>
326inline NewType cast(const OldType& x)
327{
328  return cast_impl<OldType, NewType>::run(x);
329}
330
331/****************************************************************************
332* Implementation of atanh2                                                *
333****************************************************************************/
334
335template<typename Scalar, bool IsInteger>
336struct 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
353template<typename Scalar>
354struct 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
363template<typename Scalar>
364struct atanh2_impl : atanh2_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
365
366template<typename Scalar>
367struct atanh2_retval
368{
369  typedef Scalar type;
370};
371
372/****************************************************************************
373* Implementation of pow                                                  *
374****************************************************************************/
375
376template<typename Scalar, bool IsInteger>
377struct 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
387template<typename Scalar>
388struct 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
406template<typename Scalar>
407struct pow_impl : pow_default_impl<Scalar, NumTraits<Scalar>::IsInteger> {};
408
409template<typename Scalar>
410struct pow_retval
411{
412  typedef Scalar type;
413};
414
415/****************************************************************************
416* Implementation of random                                               *
417****************************************************************************/
418
419template<typename Scalar,
420         bool IsComplex,
421         bool IsInteger>
422struct random_default_impl {};
423
424template<typename Scalar>
425struct random_impl : random_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
426
427template<typename Scalar>
428struct random_retval
429{
430  typedef Scalar type;
431};
432
433template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random(const Scalar& x, const Scalar& y);
434template<typename Scalar> inline EIGEN_MATHFUNC_RETVAL(random, Scalar) random();
435
436template<typename Scalar>
437struct 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
449enum {
450  floor_log2_terminate,
451  floor_log2_move_up,
452  floor_log2_move_down,
453  floor_log2_bogus
454};
455
456template<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
466template<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>
470struct floor_log2 {};
471
472template<unsigned int n, int lower, int upper>
473struct 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
478template<unsigned int n, int lower, int upper>
479struct 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
484template<unsigned int n, int lower, int upper>
485struct floor_log2<n, lower, upper, floor_log2_terminate>
486{
487  enum { value = (n >= ((unsigned int)(1) << (lower+1))) ? lower+1 : lower };
488};
489
490template<unsigned int n, int lower, int upper>
491struct floor_log2<n, lower, upper, floor_log2_bogus>
492{
493  // no value, error at compile time
494};
495
496template<typename Scalar>
497struct 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
521template<typename Scalar>
522struct 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
536template<typename Scalar>
537inline 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
542template<typename Scalar>
543inline 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
554namespace numext {
555
556template<typename Scalar>
557inline EIGEN_MATHFUNC_RETVAL(real, Scalar) real(const Scalar& x)
558{
559  return EIGEN_MATHFUNC_IMPL(real, Scalar)::run(x);
560}
561
562template<typename Scalar>
563inline 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
568template<typename Scalar>
569inline EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) real_ref(Scalar& x)
570{
571  return EIGEN_MATHFUNC_IMPL(real_ref, Scalar)::run(x);
572}
573
574template<typename Scalar>
575inline EIGEN_MATHFUNC_RETVAL(imag, Scalar) imag(const Scalar& x)
576{
577  return EIGEN_MATHFUNC_IMPL(imag, Scalar)::run(x);
578}
579
580template<typename Scalar>
581inline 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
586template<typename Scalar>
587inline EIGEN_MATHFUNC_RETVAL(imag_ref, Scalar) imag_ref(Scalar& x)
588{
589  return EIGEN_MATHFUNC_IMPL(imag_ref, Scalar)::run(x);
590}
591
592template<typename Scalar>
593inline EIGEN_MATHFUNC_RETVAL(conj, Scalar) conj(const Scalar& x)
594{
595  return EIGEN_MATHFUNC_IMPL(conj, Scalar)::run(x);
596}
597
598template<typename Scalar>
599inline EIGEN_MATHFUNC_RETVAL(abs2, Scalar) abs2(const Scalar& x)
600{
601  return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x);
602}
603
604template<typename Scalar>
605inline EIGEN_MATHFUNC_RETVAL(norm1, Scalar) norm1(const Scalar& x)
606{
607  return EIGEN_MATHFUNC_IMPL(norm1, Scalar)::run(x);
608}
609
610template<typename Scalar>
611inline 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
616template<typename Scalar>
617inline 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
622template<typename Scalar>
623inline 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.
630template<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
637namespace internal {
638
639/****************************************************************************
640* Implementation of fuzzy comparisons                                       *
641****************************************************************************/
642
643template<typename Scalar,
644         bool IsComplex,
645         bool IsInteger>
646struct scalar_fuzzy_default_impl {};
647
648template<typename Scalar>
649struct 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
670template<typename Scalar>
671struct 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
689template<typename Scalar>
690struct 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
705template<typename Scalar>
706struct scalar_fuzzy_impl : scalar_fuzzy_default_impl<Scalar, NumTraits<Scalar>::IsComplex, NumTraits<Scalar>::IsInteger> {};
707
708template<typename Scalar, typename OtherScalar>
709inline 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
715template<typename Scalar>
716inline 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
722template<typename Scalar>
723inline 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
733template<> struct random_impl<bool>
734{
735  static inline bool run()
736  {
737    return random<int>(0,1)==0 ? false : true;
738  }
739};
740
741template<> 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
769