1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
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_FUNCTORS_H
11#define EIGEN_FUNCTORS_H
12
13namespace Eigen {
14
15namespace internal {
16
17// associative functors:
18
19/** \internal
20  * \brief Template functor to compute the sum of two scalars
21  *
22  * \sa class CwiseBinaryOp, MatrixBase::operator+, class VectorwiseOp, MatrixBase::sum()
23  */
24template<typename Scalar> struct scalar_sum_op {
25  EIGEN_EMPTY_STRUCT_CTOR(scalar_sum_op)
26  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; }
27  template<typename Packet>
28  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
29  { return internal::padd(a,b); }
30  template<typename Packet>
31  EIGEN_STRONG_INLINE const Scalar predux(const Packet& a) const
32  { return internal::predux(a); }
33};
34template<typename Scalar>
35struct functor_traits<scalar_sum_op<Scalar> > {
36  enum {
37    Cost = NumTraits<Scalar>::AddCost,
38    PacketAccess = packet_traits<Scalar>::HasAdd
39  };
40};
41
42/** \internal
43  * \brief Template functor to compute the product of two scalars
44  *
45  * \sa class CwiseBinaryOp, Cwise::operator*(), class VectorwiseOp, MatrixBase::redux()
46  */
47template<typename LhsScalar,typename RhsScalar> struct scalar_product_op {
48  enum {
49    // TODO vectorize mixed product
50    Vectorizable = is_same<LhsScalar,RhsScalar>::value && packet_traits<LhsScalar>::HasMul && packet_traits<RhsScalar>::HasMul
51  };
52  typedef typename scalar_product_traits<LhsScalar,RhsScalar>::ReturnType result_type;
53  EIGEN_EMPTY_STRUCT_CTOR(scalar_product_op)
54  EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return a * b; }
55  template<typename Packet>
56  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
57  { return internal::pmul(a,b); }
58  template<typename Packet>
59  EIGEN_STRONG_INLINE const result_type predux(const Packet& a) const
60  { return internal::predux_mul(a); }
61};
62template<typename LhsScalar,typename RhsScalar>
63struct functor_traits<scalar_product_op<LhsScalar,RhsScalar> > {
64  enum {
65    Cost = (NumTraits<LhsScalar>::MulCost + NumTraits<RhsScalar>::MulCost)/2, // rough estimate!
66    PacketAccess = scalar_product_op<LhsScalar,RhsScalar>::Vectorizable
67  };
68};
69
70/** \internal
71  * \brief Template functor to compute the conjugate product of two scalars
72  *
73  * This is a short cut for conj(x) * y which is needed for optimization purpose; in Eigen2 support mode, this becomes x * conj(y)
74  */
75template<typename LhsScalar,typename RhsScalar> struct scalar_conj_product_op {
76
77  enum {
78    Conj = NumTraits<LhsScalar>::IsComplex
79  };
80
81  typedef typename scalar_product_traits<LhsScalar,RhsScalar>::ReturnType result_type;
82
83  EIGEN_EMPTY_STRUCT_CTOR(scalar_conj_product_op)
84  EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const
85  { return conj_helper<LhsScalar,RhsScalar,Conj,false>().pmul(a,b); }
86
87  template<typename Packet>
88  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
89  { return conj_helper<Packet,Packet,Conj,false>().pmul(a,b); }
90};
91template<typename LhsScalar,typename RhsScalar>
92struct functor_traits<scalar_conj_product_op<LhsScalar,RhsScalar> > {
93  enum {
94    Cost = NumTraits<LhsScalar>::MulCost,
95    PacketAccess = internal::is_same<LhsScalar, RhsScalar>::value && packet_traits<LhsScalar>::HasMul
96  };
97};
98
99/** \internal
100  * \brief Template functor to compute the min of two scalars
101  *
102  * \sa class CwiseBinaryOp, MatrixBase::cwiseMin, class VectorwiseOp, MatrixBase::minCoeff()
103  */
104template<typename Scalar> struct scalar_min_op {
105  EIGEN_EMPTY_STRUCT_CTOR(scalar_min_op)
106  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { using std::min; return (min)(a, b); }
107  template<typename Packet>
108  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
109  { return internal::pmin(a,b); }
110  template<typename Packet>
111  EIGEN_STRONG_INLINE const Scalar predux(const Packet& a) const
112  { return internal::predux_min(a); }
113};
114template<typename Scalar>
115struct functor_traits<scalar_min_op<Scalar> > {
116  enum {
117    Cost = NumTraits<Scalar>::AddCost,
118    PacketAccess = packet_traits<Scalar>::HasMin
119  };
120};
121
122/** \internal
123  * \brief Template functor to compute the max of two scalars
124  *
125  * \sa class CwiseBinaryOp, MatrixBase::cwiseMax, class VectorwiseOp, MatrixBase::maxCoeff()
126  */
127template<typename Scalar> struct scalar_max_op {
128  EIGEN_EMPTY_STRUCT_CTOR(scalar_max_op)
129  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { using std::max; return (max)(a, b); }
130  template<typename Packet>
131  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
132  { return internal::pmax(a,b); }
133  template<typename Packet>
134  EIGEN_STRONG_INLINE const Scalar predux(const Packet& a) const
135  { return internal::predux_max(a); }
136};
137template<typename Scalar>
138struct functor_traits<scalar_max_op<Scalar> > {
139  enum {
140    Cost = NumTraits<Scalar>::AddCost,
141    PacketAccess = packet_traits<Scalar>::HasMax
142  };
143};
144
145/** \internal
146  * \brief Template functor to compute the hypot of two scalars
147  *
148  * \sa MatrixBase::stableNorm(), class Redux
149  */
150template<typename Scalar> struct scalar_hypot_op {
151  EIGEN_EMPTY_STRUCT_CTOR(scalar_hypot_op)
152//   typedef typename NumTraits<Scalar>::Real result_type;
153  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const
154  {
155    using std::max;
156    using std::min;
157    using std::sqrt;
158    Scalar p = (max)(_x, _y);
159    Scalar q = (min)(_x, _y);
160    Scalar qp = q/p;
161    return p * sqrt(Scalar(1) + qp*qp);
162  }
163};
164template<typename Scalar>
165struct functor_traits<scalar_hypot_op<Scalar> > {
166  enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess=0 };
167};
168
169/** \internal
170  * \brief Template functor to compute the pow of two scalars
171  */
172template<typename Scalar, typename OtherScalar> struct scalar_binary_pow_op {
173  EIGEN_EMPTY_STRUCT_CTOR(scalar_binary_pow_op)
174  inline Scalar operator() (const Scalar& a, const OtherScalar& b) const { return numext::pow(a, b); }
175};
176template<typename Scalar, typename OtherScalar>
177struct functor_traits<scalar_binary_pow_op<Scalar,OtherScalar> > {
178  enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = false };
179};
180
181// other binary functors:
182
183/** \internal
184  * \brief Template functor to compute the difference of two scalars
185  *
186  * \sa class CwiseBinaryOp, MatrixBase::operator-
187  */
188template<typename Scalar> struct scalar_difference_op {
189  EIGEN_EMPTY_STRUCT_CTOR(scalar_difference_op)
190  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a - b; }
191  template<typename Packet>
192  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
193  { return internal::psub(a,b); }
194};
195template<typename Scalar>
196struct functor_traits<scalar_difference_op<Scalar> > {
197  enum {
198    Cost = NumTraits<Scalar>::AddCost,
199    PacketAccess = packet_traits<Scalar>::HasSub
200  };
201};
202
203/** \internal
204  * \brief Template functor to compute the quotient of two scalars
205  *
206  * \sa class CwiseBinaryOp, Cwise::operator/()
207  */
208template<typename LhsScalar,typename RhsScalar> struct scalar_quotient_op {
209  enum {
210    // TODO vectorize mixed product
211    Vectorizable = is_same<LhsScalar,RhsScalar>::value && packet_traits<LhsScalar>::HasDiv && packet_traits<RhsScalar>::HasDiv
212  };
213  typedef typename scalar_product_traits<LhsScalar,RhsScalar>::ReturnType result_type;
214  EIGEN_EMPTY_STRUCT_CTOR(scalar_quotient_op)
215  EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return a / b; }
216  template<typename Packet>
217  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
218  { return internal::pdiv(a,b); }
219};
220template<typename LhsScalar,typename RhsScalar>
221struct functor_traits<scalar_quotient_op<LhsScalar,RhsScalar> > {
222  enum {
223    Cost = (NumTraits<LhsScalar>::MulCost + NumTraits<RhsScalar>::MulCost), // rough estimate!
224    PacketAccess = scalar_quotient_op<LhsScalar,RhsScalar>::Vectorizable
225  };
226};
227
228
229
230/** \internal
231  * \brief Template functor to compute the and of two booleans
232  *
233  * \sa class CwiseBinaryOp, ArrayBase::operator&&
234  */
235struct scalar_boolean_and_op {
236  EIGEN_EMPTY_STRUCT_CTOR(scalar_boolean_and_op)
237  EIGEN_STRONG_INLINE bool operator() (const bool& a, const bool& b) const { return a && b; }
238};
239template<> struct functor_traits<scalar_boolean_and_op> {
240  enum {
241    Cost = NumTraits<bool>::AddCost,
242    PacketAccess = false
243  };
244};
245
246/** \internal
247  * \brief Template functor to compute the or of two booleans
248  *
249  * \sa class CwiseBinaryOp, ArrayBase::operator||
250  */
251struct scalar_boolean_or_op {
252  EIGEN_EMPTY_STRUCT_CTOR(scalar_boolean_or_op)
253  EIGEN_STRONG_INLINE bool operator() (const bool& a, const bool& b) const { return a || b; }
254};
255template<> struct functor_traits<scalar_boolean_or_op> {
256  enum {
257    Cost = NumTraits<bool>::AddCost,
258    PacketAccess = false
259  };
260};
261
262// unary functors:
263
264/** \internal
265  * \brief Template functor to compute the opposite of a scalar
266  *
267  * \sa class CwiseUnaryOp, MatrixBase::operator-
268  */
269template<typename Scalar> struct scalar_opposite_op {
270  EIGEN_EMPTY_STRUCT_CTOR(scalar_opposite_op)
271  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { return -a; }
272  template<typename Packet>
273  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
274  { return internal::pnegate(a); }
275};
276template<typename Scalar>
277struct functor_traits<scalar_opposite_op<Scalar> >
278{ enum {
279    Cost = NumTraits<Scalar>::AddCost,
280    PacketAccess = packet_traits<Scalar>::HasNegate };
281};
282
283/** \internal
284  * \brief Template functor to compute the absolute value of a scalar
285  *
286  * \sa class CwiseUnaryOp, Cwise::abs
287  */
288template<typename Scalar> struct scalar_abs_op {
289  EIGEN_EMPTY_STRUCT_CTOR(scalar_abs_op)
290  typedef typename NumTraits<Scalar>::Real result_type;
291  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { using std::abs; return abs(a); }
292  template<typename Packet>
293  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
294  { return internal::pabs(a); }
295};
296template<typename Scalar>
297struct functor_traits<scalar_abs_op<Scalar> >
298{
299  enum {
300    Cost = NumTraits<Scalar>::AddCost,
301    PacketAccess = packet_traits<Scalar>::HasAbs
302  };
303};
304
305/** \internal
306  * \brief Template functor to compute the squared absolute value of a scalar
307  *
308  * \sa class CwiseUnaryOp, Cwise::abs2
309  */
310template<typename Scalar> struct scalar_abs2_op {
311  EIGEN_EMPTY_STRUCT_CTOR(scalar_abs2_op)
312  typedef typename NumTraits<Scalar>::Real result_type;
313  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return numext::abs2(a); }
314  template<typename Packet>
315  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
316  { return internal::pmul(a,a); }
317};
318template<typename Scalar>
319struct functor_traits<scalar_abs2_op<Scalar> >
320{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasAbs2 }; };
321
322/** \internal
323  * \brief Template functor to compute the conjugate of a complex value
324  *
325  * \sa class CwiseUnaryOp, MatrixBase::conjugate()
326  */
327template<typename Scalar> struct scalar_conjugate_op {
328  EIGEN_EMPTY_STRUCT_CTOR(scalar_conjugate_op)
329  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { using numext::conj; return conj(a); }
330  template<typename Packet>
331  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const { return internal::pconj(a); }
332};
333template<typename Scalar>
334struct functor_traits<scalar_conjugate_op<Scalar> >
335{
336  enum {
337    Cost = NumTraits<Scalar>::IsComplex ? NumTraits<Scalar>::AddCost : 0,
338    PacketAccess = packet_traits<Scalar>::HasConj
339  };
340};
341
342/** \internal
343  * \brief Template functor to cast a scalar to another type
344  *
345  * \sa class CwiseUnaryOp, MatrixBase::cast()
346  */
347template<typename Scalar, typename NewType>
348struct scalar_cast_op {
349  EIGEN_EMPTY_STRUCT_CTOR(scalar_cast_op)
350  typedef NewType result_type;
351  EIGEN_STRONG_INLINE const NewType operator() (const Scalar& a) const { return cast<Scalar, NewType>(a); }
352};
353template<typename Scalar, typename NewType>
354struct functor_traits<scalar_cast_op<Scalar,NewType> >
355{ enum { Cost = is_same<Scalar, NewType>::value ? 0 : NumTraits<NewType>::AddCost, PacketAccess = false }; };
356
357/** \internal
358  * \brief Template functor to extract the real part of a complex
359  *
360  * \sa class CwiseUnaryOp, MatrixBase::real()
361  */
362template<typename Scalar>
363struct scalar_real_op {
364  EIGEN_EMPTY_STRUCT_CTOR(scalar_real_op)
365  typedef typename NumTraits<Scalar>::Real result_type;
366  EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return numext::real(a); }
367};
368template<typename Scalar>
369struct functor_traits<scalar_real_op<Scalar> >
370{ enum { Cost = 0, PacketAccess = false }; };
371
372/** \internal
373  * \brief Template functor to extract the imaginary part of a complex
374  *
375  * \sa class CwiseUnaryOp, MatrixBase::imag()
376  */
377template<typename Scalar>
378struct scalar_imag_op {
379  EIGEN_EMPTY_STRUCT_CTOR(scalar_imag_op)
380  typedef typename NumTraits<Scalar>::Real result_type;
381  EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return numext::imag(a); }
382};
383template<typename Scalar>
384struct functor_traits<scalar_imag_op<Scalar> >
385{ enum { Cost = 0, PacketAccess = false }; };
386
387/** \internal
388  * \brief Template functor to extract the real part of a complex as a reference
389  *
390  * \sa class CwiseUnaryOp, MatrixBase::real()
391  */
392template<typename Scalar>
393struct scalar_real_ref_op {
394  EIGEN_EMPTY_STRUCT_CTOR(scalar_real_ref_op)
395  typedef typename NumTraits<Scalar>::Real result_type;
396  EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return numext::real_ref(*const_cast<Scalar*>(&a)); }
397};
398template<typename Scalar>
399struct functor_traits<scalar_real_ref_op<Scalar> >
400{ enum { Cost = 0, PacketAccess = false }; };
401
402/** \internal
403  * \brief Template functor to extract the imaginary part of a complex as a reference
404  *
405  * \sa class CwiseUnaryOp, MatrixBase::imag()
406  */
407template<typename Scalar>
408struct scalar_imag_ref_op {
409  EIGEN_EMPTY_STRUCT_CTOR(scalar_imag_ref_op)
410  typedef typename NumTraits<Scalar>::Real result_type;
411  EIGEN_STRONG_INLINE result_type& operator() (const Scalar& a) const { return numext::imag_ref(*const_cast<Scalar*>(&a)); }
412};
413template<typename Scalar>
414struct functor_traits<scalar_imag_ref_op<Scalar> >
415{ enum { Cost = 0, PacketAccess = false }; };
416
417/** \internal
418  *
419  * \brief Template functor to compute the exponential of a scalar
420  *
421  * \sa class CwiseUnaryOp, Cwise::exp()
422  */
423template<typename Scalar> struct scalar_exp_op {
424  EIGEN_EMPTY_STRUCT_CTOR(scalar_exp_op)
425  inline const Scalar operator() (const Scalar& a) const { using std::exp; return exp(a); }
426  typedef typename packet_traits<Scalar>::type Packet;
427  inline Packet packetOp(const Packet& a) const { return internal::pexp(a); }
428};
429template<typename Scalar>
430struct functor_traits<scalar_exp_op<Scalar> >
431{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasExp }; };
432
433/** \internal
434  *
435  * \brief Template functor to compute the logarithm of a scalar
436  *
437  * \sa class CwiseUnaryOp, Cwise::log()
438  */
439template<typename Scalar> struct scalar_log_op {
440  EIGEN_EMPTY_STRUCT_CTOR(scalar_log_op)
441  inline const Scalar operator() (const Scalar& a) const { using std::log; return log(a); }
442  typedef typename packet_traits<Scalar>::type Packet;
443  inline Packet packetOp(const Packet& a) const { return internal::plog(a); }
444};
445template<typename Scalar>
446struct functor_traits<scalar_log_op<Scalar> >
447{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasLog }; };
448
449/** \internal
450  * \brief Template functor to multiply a scalar by a fixed other one
451  *
452  * \sa class CwiseUnaryOp, MatrixBase::operator*, MatrixBase::operator/
453  */
454/* NOTE why doing the pset1() in packetOp *is* an optimization ?
455 * indeed it seems better to declare m_other as a Packet and do the pset1() once
456 * in the constructor. However, in practice:
457 *  - GCC does not like m_other as a Packet and generate a load every time it needs it
458 *  - on the other hand GCC is able to moves the pset1() outside the loop :)
459 *  - simpler code ;)
460 * (ICC and gcc 4.4 seems to perform well in both cases, the issue is visible with y = a*x + b*y)
461 */
462template<typename Scalar>
463struct scalar_multiple_op {
464  typedef typename packet_traits<Scalar>::type Packet;
465  // FIXME default copy constructors seems bugged with std::complex<>
466  EIGEN_STRONG_INLINE scalar_multiple_op(const scalar_multiple_op& other) : m_other(other.m_other) { }
467  EIGEN_STRONG_INLINE scalar_multiple_op(const Scalar& other) : m_other(other) { }
468  EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a * m_other; }
469  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
470  { return internal::pmul(a, pset1<Packet>(m_other)); }
471  typename add_const_on_value_type<typename NumTraits<Scalar>::Nested>::type m_other;
472};
473template<typename Scalar>
474struct functor_traits<scalar_multiple_op<Scalar> >
475{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasMul }; };
476
477template<typename Scalar1, typename Scalar2>
478struct scalar_multiple2_op {
479  typedef typename scalar_product_traits<Scalar1,Scalar2>::ReturnType result_type;
480  EIGEN_STRONG_INLINE scalar_multiple2_op(const scalar_multiple2_op& other) : m_other(other.m_other) { }
481  EIGEN_STRONG_INLINE scalar_multiple2_op(const Scalar2& other) : m_other(other) { }
482  EIGEN_STRONG_INLINE result_type operator() (const Scalar1& a) const { return a * m_other; }
483  typename add_const_on_value_type<typename NumTraits<Scalar2>::Nested>::type m_other;
484};
485template<typename Scalar1,typename Scalar2>
486struct functor_traits<scalar_multiple2_op<Scalar1,Scalar2> >
487{ enum { Cost = NumTraits<Scalar1>::MulCost, PacketAccess = false }; };
488
489/** \internal
490  * \brief Template functor to divide a scalar by a fixed other one
491  *
492  * This functor is used to implement the quotient of a matrix by
493  * a scalar where the scalar type is not necessarily a floating point type.
494  *
495  * \sa class CwiseUnaryOp, MatrixBase::operator/
496  */
497template<typename Scalar>
498struct scalar_quotient1_op {
499  typedef typename packet_traits<Scalar>::type Packet;
500  // FIXME default copy constructors seems bugged with std::complex<>
501  EIGEN_STRONG_INLINE scalar_quotient1_op(const scalar_quotient1_op& other) : m_other(other.m_other) { }
502  EIGEN_STRONG_INLINE scalar_quotient1_op(const Scalar& other) : m_other(other) {}
503  EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; }
504  EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const
505  { return internal::pdiv(a, pset1<Packet>(m_other)); }
506  typename add_const_on_value_type<typename NumTraits<Scalar>::Nested>::type m_other;
507};
508template<typename Scalar>
509struct functor_traits<scalar_quotient1_op<Scalar> >
510{ enum { Cost = 2 * NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasDiv }; };
511
512// nullary functors
513
514template<typename Scalar>
515struct scalar_constant_op {
516  typedef typename packet_traits<Scalar>::type Packet;
517  EIGEN_STRONG_INLINE scalar_constant_op(const scalar_constant_op& other) : m_other(other.m_other) { }
518  EIGEN_STRONG_INLINE scalar_constant_op(const Scalar& other) : m_other(other) { }
519  template<typename Index>
520  EIGEN_STRONG_INLINE const Scalar operator() (Index, Index = 0) const { return m_other; }
521  template<typename Index>
522  EIGEN_STRONG_INLINE const Packet packetOp(Index, Index = 0) const { return internal::pset1<Packet>(m_other); }
523  const Scalar m_other;
524};
525template<typename Scalar>
526struct functor_traits<scalar_constant_op<Scalar> >
527// FIXME replace this packet test by a safe one
528{ enum { Cost = 1, PacketAccess = packet_traits<Scalar>::Vectorizable, IsRepeatable = true }; };
529
530template<typename Scalar> struct scalar_identity_op {
531  EIGEN_EMPTY_STRUCT_CTOR(scalar_identity_op)
532  template<typename Index>
533  EIGEN_STRONG_INLINE const Scalar operator() (Index row, Index col) const { return row==col ? Scalar(1) : Scalar(0); }
534};
535template<typename Scalar>
536struct functor_traits<scalar_identity_op<Scalar> >
537{ enum { Cost = NumTraits<Scalar>::AddCost, PacketAccess = false, IsRepeatable = true }; };
538
539template <typename Scalar, bool RandomAccess> struct linspaced_op_impl;
540
541// linear access for packet ops:
542// 1) initialization
543//   base = [low, ..., low] + ([step, ..., step] * [-size, ..., 0])
544// 2) each step (where size is 1 for coeff access or PacketSize for packet access)
545//   base += [size*step, ..., size*step]
546//
547// TODO: Perhaps it's better to initialize lazily (so not in the constructor but in packetOp)
548//       in order to avoid the padd() in operator() ?
549template <typename Scalar>
550struct linspaced_op_impl<Scalar,false>
551{
552  typedef typename packet_traits<Scalar>::type Packet;
553
554  linspaced_op_impl(const Scalar& low, const Scalar& step) :
555  m_low(low), m_step(step),
556  m_packetStep(pset1<Packet>(packet_traits<Scalar>::size*step)),
557  m_base(padd(pset1<Packet>(low), pmul(pset1<Packet>(step),plset<Scalar>(-packet_traits<Scalar>::size)))) {}
558
559  template<typename Index>
560  EIGEN_STRONG_INLINE const Scalar operator() (Index i) const
561  {
562    m_base = padd(m_base, pset1<Packet>(m_step));
563    return m_low+Scalar(i)*m_step;
564  }
565
566  template<typename Index>
567  EIGEN_STRONG_INLINE const Packet packetOp(Index) const { return m_base = padd(m_base,m_packetStep); }
568
569  const Scalar m_low;
570  const Scalar m_step;
571  const Packet m_packetStep;
572  mutable Packet m_base;
573};
574
575// random access for packet ops:
576// 1) each step
577//   [low, ..., low] + ( [step, ..., step] * ( [i, ..., i] + [0, ..., size] ) )
578template <typename Scalar>
579struct linspaced_op_impl<Scalar,true>
580{
581  typedef typename packet_traits<Scalar>::type Packet;
582
583  linspaced_op_impl(const Scalar& low, const Scalar& step) :
584  m_low(low), m_step(step),
585  m_lowPacket(pset1<Packet>(m_low)), m_stepPacket(pset1<Packet>(m_step)), m_interPacket(plset<Scalar>(0)) {}
586
587  template<typename Index>
588  EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return m_low+i*m_step; }
589
590  template<typename Index>
591  EIGEN_STRONG_INLINE const Packet packetOp(Index i) const
592  { return internal::padd(m_lowPacket, pmul(m_stepPacket, padd(pset1<Packet>(Scalar(i)),m_interPacket))); }
593
594  const Scalar m_low;
595  const Scalar m_step;
596  const Packet m_lowPacket;
597  const Packet m_stepPacket;
598  const Packet m_interPacket;
599};
600
601// ----- Linspace functor ----------------------------------------------------------------
602
603// Forward declaration (we default to random access which does not really give
604// us a speed gain when using packet access but it allows to use the functor in
605// nested expressions).
606template <typename Scalar, bool RandomAccess = true> struct linspaced_op;
607template <typename Scalar, bool RandomAccess> struct functor_traits< linspaced_op<Scalar,RandomAccess> >
608{ enum { Cost = 1, PacketAccess = packet_traits<Scalar>::HasSetLinear, IsRepeatable = true }; };
609template <typename Scalar, bool RandomAccess> struct linspaced_op
610{
611  typedef typename packet_traits<Scalar>::type Packet;
612  linspaced_op(const Scalar& low, const Scalar& high, DenseIndex num_steps) : impl((num_steps==1 ? high : low), (num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1))) {}
613
614  template<typename Index>
615  EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); }
616
617  // We need this function when assigning e.g. a RowVectorXd to a MatrixXd since
618  // there row==0 and col is used for the actual iteration.
619  template<typename Index>
620  EIGEN_STRONG_INLINE const Scalar operator() (Index row, Index col) const
621  {
622    eigen_assert(col==0 || row==0);
623    return impl(col + row);
624  }
625
626  template<typename Index>
627  EIGEN_STRONG_INLINE const Packet packetOp(Index i) const { return impl.packetOp(i); }
628
629  // We need this function when assigning e.g. a RowVectorXd to a MatrixXd since
630  // there row==0 and col is used for the actual iteration.
631  template<typename Index>
632  EIGEN_STRONG_INLINE const Packet packetOp(Index row, Index col) const
633  {
634    eigen_assert(col==0 || row==0);
635    return impl.packetOp(col + row);
636  }
637
638  // This proxy object handles the actual required temporaries, the different
639  // implementations (random vs. sequential access) as well as the
640  // correct piping to size 2/4 packet operations.
641  const linspaced_op_impl<Scalar,RandomAccess> impl;
642};
643
644// all functors allow linear access, except scalar_identity_op. So we fix here a quick meta
645// to indicate whether a functor allows linear access, just always answering 'yes' except for
646// scalar_identity_op.
647// FIXME move this to functor_traits adding a functor_default
648template<typename Functor> struct functor_has_linear_access { enum { ret = 1 }; };
649template<typename Scalar> struct functor_has_linear_access<scalar_identity_op<Scalar> > { enum { ret = 0 }; };
650
651// In Eigen, any binary op (Product, CwiseBinaryOp) require the Lhs and Rhs to have the same scalar type, except for multiplication
652// where the mixing of different types is handled by scalar_product_traits
653// In particular, real * complex<real> is allowed.
654// FIXME move this to functor_traits adding a functor_default
655template<typename Functor> struct functor_is_product_like { enum { ret = 0 }; };
656template<typename LhsScalar,typename RhsScalar> struct functor_is_product_like<scalar_product_op<LhsScalar,RhsScalar> > { enum { ret = 1 }; };
657template<typename LhsScalar,typename RhsScalar> struct functor_is_product_like<scalar_conj_product_op<LhsScalar,RhsScalar> > { enum { ret = 1 }; };
658template<typename LhsScalar,typename RhsScalar> struct functor_is_product_like<scalar_quotient_op<LhsScalar,RhsScalar> > { enum { ret = 1 }; };
659
660
661/** \internal
662  * \brief Template functor to add a scalar to a fixed other one
663  * \sa class CwiseUnaryOp, Array::operator+
664  */
665/* If you wonder why doing the pset1() in packetOp() is an optimization check scalar_multiple_op */
666template<typename Scalar>
667struct scalar_add_op {
668  typedef typename packet_traits<Scalar>::type Packet;
669  // FIXME default copy constructors seems bugged with std::complex<>
670  inline scalar_add_op(const scalar_add_op& other) : m_other(other.m_other) { }
671  inline scalar_add_op(const Scalar& other) : m_other(other) { }
672  inline Scalar operator() (const Scalar& a) const { return a + m_other; }
673  inline const Packet packetOp(const Packet& a) const
674  { return internal::padd(a, pset1<Packet>(m_other)); }
675  const Scalar m_other;
676};
677template<typename Scalar>
678struct functor_traits<scalar_add_op<Scalar> >
679{ enum { Cost = NumTraits<Scalar>::AddCost, PacketAccess = packet_traits<Scalar>::HasAdd }; };
680
681/** \internal
682  * \brief Template functor to compute the square root of a scalar
683  * \sa class CwiseUnaryOp, Cwise::sqrt()
684  */
685template<typename Scalar> struct scalar_sqrt_op {
686  EIGEN_EMPTY_STRUCT_CTOR(scalar_sqrt_op)
687  inline const Scalar operator() (const Scalar& a) const { using std::sqrt; return sqrt(a); }
688  typedef typename packet_traits<Scalar>::type Packet;
689  inline Packet packetOp(const Packet& a) const { return internal::psqrt(a); }
690};
691template<typename Scalar>
692struct functor_traits<scalar_sqrt_op<Scalar> >
693{ enum {
694    Cost = 5 * NumTraits<Scalar>::MulCost,
695    PacketAccess = packet_traits<Scalar>::HasSqrt
696  };
697};
698
699/** \internal
700  * \brief Template functor to compute the cosine of a scalar
701  * \sa class CwiseUnaryOp, ArrayBase::cos()
702  */
703template<typename Scalar> struct scalar_cos_op {
704  EIGEN_EMPTY_STRUCT_CTOR(scalar_cos_op)
705  inline Scalar operator() (const Scalar& a) const { using std::cos; return cos(a); }
706  typedef typename packet_traits<Scalar>::type Packet;
707  inline Packet packetOp(const Packet& a) const { return internal::pcos(a); }
708};
709template<typename Scalar>
710struct functor_traits<scalar_cos_op<Scalar> >
711{
712  enum {
713    Cost = 5 * NumTraits<Scalar>::MulCost,
714    PacketAccess = packet_traits<Scalar>::HasCos
715  };
716};
717
718/** \internal
719  * \brief Template functor to compute the sine of a scalar
720  * \sa class CwiseUnaryOp, ArrayBase::sin()
721  */
722template<typename Scalar> struct scalar_sin_op {
723  EIGEN_EMPTY_STRUCT_CTOR(scalar_sin_op)
724  inline const Scalar operator() (const Scalar& a) const { using std::sin; return sin(a); }
725  typedef typename packet_traits<Scalar>::type Packet;
726  inline Packet packetOp(const Packet& a) const { return internal::psin(a); }
727};
728template<typename Scalar>
729struct functor_traits<scalar_sin_op<Scalar> >
730{
731  enum {
732    Cost = 5 * NumTraits<Scalar>::MulCost,
733    PacketAccess = packet_traits<Scalar>::HasSin
734  };
735};
736
737
738/** \internal
739  * \brief Template functor to compute the tan of a scalar
740  * \sa class CwiseUnaryOp, ArrayBase::tan()
741  */
742template<typename Scalar> struct scalar_tan_op {
743  EIGEN_EMPTY_STRUCT_CTOR(scalar_tan_op)
744  inline const Scalar operator() (const Scalar& a) const { using std::tan; return tan(a); }
745  typedef typename packet_traits<Scalar>::type Packet;
746  inline Packet packetOp(const Packet& a) const { return internal::ptan(a); }
747};
748template<typename Scalar>
749struct functor_traits<scalar_tan_op<Scalar> >
750{
751  enum {
752    Cost = 5 * NumTraits<Scalar>::MulCost,
753    PacketAccess = packet_traits<Scalar>::HasTan
754  };
755};
756
757/** \internal
758  * \brief Template functor to compute the arc cosine of a scalar
759  * \sa class CwiseUnaryOp, ArrayBase::acos()
760  */
761template<typename Scalar> struct scalar_acos_op {
762  EIGEN_EMPTY_STRUCT_CTOR(scalar_acos_op)
763  inline const Scalar operator() (const Scalar& a) const { using std::acos; return acos(a); }
764  typedef typename packet_traits<Scalar>::type Packet;
765  inline Packet packetOp(const Packet& a) const { return internal::pacos(a); }
766};
767template<typename Scalar>
768struct functor_traits<scalar_acos_op<Scalar> >
769{
770  enum {
771    Cost = 5 * NumTraits<Scalar>::MulCost,
772    PacketAccess = packet_traits<Scalar>::HasACos
773  };
774};
775
776/** \internal
777  * \brief Template functor to compute the arc sine of a scalar
778  * \sa class CwiseUnaryOp, ArrayBase::asin()
779  */
780template<typename Scalar> struct scalar_asin_op {
781  EIGEN_EMPTY_STRUCT_CTOR(scalar_asin_op)
782  inline const Scalar operator() (const Scalar& a) const { using std::asin; return asin(a); }
783  typedef typename packet_traits<Scalar>::type Packet;
784  inline Packet packetOp(const Packet& a) const { return internal::pasin(a); }
785};
786template<typename Scalar>
787struct functor_traits<scalar_asin_op<Scalar> >
788{
789  enum {
790    Cost = 5 * NumTraits<Scalar>::MulCost,
791    PacketAccess = packet_traits<Scalar>::HasASin
792  };
793};
794
795/** \internal
796  * \brief Template functor to raise a scalar to a power
797  * \sa class CwiseUnaryOp, Cwise::pow
798  */
799template<typename Scalar>
800struct scalar_pow_op {
801  // FIXME default copy constructors seems bugged with std::complex<>
802  inline scalar_pow_op(const scalar_pow_op& other) : m_exponent(other.m_exponent) { }
803  inline scalar_pow_op(const Scalar& exponent) : m_exponent(exponent) {}
804  inline Scalar operator() (const Scalar& a) const { return numext::pow(a, m_exponent); }
805  const Scalar m_exponent;
806};
807template<typename Scalar>
808struct functor_traits<scalar_pow_op<Scalar> >
809{ enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = false }; };
810
811/** \internal
812  * \brief Template functor to compute the quotient between a scalar and array entries.
813  * \sa class CwiseUnaryOp, Cwise::inverse()
814  */
815template<typename Scalar>
816struct scalar_inverse_mult_op {
817  scalar_inverse_mult_op(const Scalar& other) : m_other(other) {}
818  inline Scalar operator() (const Scalar& a) const { return m_other / a; }
819  template<typename Packet>
820  inline const Packet packetOp(const Packet& a) const
821  { return internal::pdiv(pset1<Packet>(m_other),a); }
822  Scalar m_other;
823};
824
825/** \internal
826  * \brief Template functor to compute the inverse of a scalar
827  * \sa class CwiseUnaryOp, Cwise::inverse()
828  */
829template<typename Scalar>
830struct scalar_inverse_op {
831  EIGEN_EMPTY_STRUCT_CTOR(scalar_inverse_op)
832  inline Scalar operator() (const Scalar& a) const { return Scalar(1)/a; }
833  template<typename Packet>
834  inline const Packet packetOp(const Packet& a) const
835  { return internal::pdiv(pset1<Packet>(Scalar(1)),a); }
836};
837template<typename Scalar>
838struct functor_traits<scalar_inverse_op<Scalar> >
839{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasDiv }; };
840
841/** \internal
842  * \brief Template functor to compute the square of a scalar
843  * \sa class CwiseUnaryOp, Cwise::square()
844  */
845template<typename Scalar>
846struct scalar_square_op {
847  EIGEN_EMPTY_STRUCT_CTOR(scalar_square_op)
848  inline Scalar operator() (const Scalar& a) const { return a*a; }
849  template<typename Packet>
850  inline const Packet packetOp(const Packet& a) const
851  { return internal::pmul(a,a); }
852};
853template<typename Scalar>
854struct functor_traits<scalar_square_op<Scalar> >
855{ enum { Cost = NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasMul }; };
856
857/** \internal
858  * \brief Template functor to compute the cube of a scalar
859  * \sa class CwiseUnaryOp, Cwise::cube()
860  */
861template<typename Scalar>
862struct scalar_cube_op {
863  EIGEN_EMPTY_STRUCT_CTOR(scalar_cube_op)
864  inline Scalar operator() (const Scalar& a) const { return a*a*a; }
865  template<typename Packet>
866  inline const Packet packetOp(const Packet& a) const
867  { return internal::pmul(a,pmul(a,a)); }
868};
869template<typename Scalar>
870struct functor_traits<scalar_cube_op<Scalar> >
871{ enum { Cost = 2*NumTraits<Scalar>::MulCost, PacketAccess = packet_traits<Scalar>::HasMul }; };
872
873// default functor traits for STL functors:
874
875template<typename T>
876struct functor_traits<std::multiplies<T> >
877{ enum { Cost = NumTraits<T>::MulCost, PacketAccess = false }; };
878
879template<typename T>
880struct functor_traits<std::divides<T> >
881{ enum { Cost = NumTraits<T>::MulCost, PacketAccess = false }; };
882
883template<typename T>
884struct functor_traits<std::plus<T> >
885{ enum { Cost = NumTraits<T>::AddCost, PacketAccess = false }; };
886
887template<typename T>
888struct functor_traits<std::minus<T> >
889{ enum { Cost = NumTraits<T>::AddCost, PacketAccess = false }; };
890
891template<typename T>
892struct functor_traits<std::negate<T> >
893{ enum { Cost = NumTraits<T>::AddCost, PacketAccess = false }; };
894
895template<typename T>
896struct functor_traits<std::logical_or<T> >
897{ enum { Cost = 1, PacketAccess = false }; };
898
899template<typename T>
900struct functor_traits<std::logical_and<T> >
901{ enum { Cost = 1, PacketAccess = false }; };
902
903template<typename T>
904struct functor_traits<std::logical_not<T> >
905{ enum { Cost = 1, PacketAccess = false }; };
906
907template<typename T>
908struct functor_traits<std::greater<T> >
909{ enum { Cost = 1, PacketAccess = false }; };
910
911template<typename T>
912struct functor_traits<std::less<T> >
913{ enum { Cost = 1, PacketAccess = false }; };
914
915template<typename T>
916struct functor_traits<std::greater_equal<T> >
917{ enum { Cost = 1, PacketAccess = false }; };
918
919template<typename T>
920struct functor_traits<std::less_equal<T> >
921{ enum { Cost = 1, PacketAccess = false }; };
922
923template<typename T>
924struct functor_traits<std::equal_to<T> >
925{ enum { Cost = 1, PacketAccess = false }; };
926
927template<typename T>
928struct functor_traits<std::not_equal_to<T> >
929{ enum { Cost = 1, PacketAccess = false }; };
930
931template<typename T>
932struct functor_traits<std::binder2nd<T> >
933{ enum { Cost = functor_traits<T>::Cost, PacketAccess = false }; };
934
935template<typename T>
936struct functor_traits<std::binder1st<T> >
937{ enum { Cost = functor_traits<T>::Cost, PacketAccess = false }; };
938
939template<typename T>
940struct functor_traits<std::unary_negate<T> >
941{ enum { Cost = 1 + functor_traits<T>::Cost, PacketAccess = false }; };
942
943template<typename T>
944struct functor_traits<std::binary_negate<T> >
945{ enum { Cost = 1 + functor_traits<T>::Cost, PacketAccess = false }; };
946
947#ifdef EIGEN_STDEXT_SUPPORT
948
949template<typename T0,typename T1>
950struct functor_traits<std::project1st<T0,T1> >
951{ enum { Cost = 0, PacketAccess = false }; };
952
953template<typename T0,typename T1>
954struct functor_traits<std::project2nd<T0,T1> >
955{ enum { Cost = 0, PacketAccess = false }; };
956
957template<typename T0,typename T1>
958struct functor_traits<std::select2nd<std::pair<T0,T1> > >
959{ enum { Cost = 0, PacketAccess = false }; };
960
961template<typename T0,typename T1>
962struct functor_traits<std::select1st<std::pair<T0,T1> > >
963{ enum { Cost = 0, PacketAccess = false }; };
964
965template<typename T0,typename T1>
966struct functor_traits<std::unary_compose<T0,T1> >
967{ enum { Cost = functor_traits<T0>::Cost + functor_traits<T1>::Cost, PacketAccess = false }; };
968
969template<typename T0,typename T1,typename T2>
970struct functor_traits<std::binary_compose<T0,T1,T2> >
971{ enum { Cost = functor_traits<T0>::Cost + functor_traits<T1>::Cost + functor_traits<T2>::Cost, PacketAccess = false }; };
972
973#endif // EIGEN_STDEXT_SUPPORT
974
975// allow to add new functors and specializations of functor_traits from outside Eigen.
976// this macro is really needed because functor_traits must be specialized after it is declared but before it is used...
977#ifdef EIGEN_FUNCTORS_PLUGIN
978#include EIGEN_FUNCTORS_PLUGIN
979#endif
980
981} // end namespace internal
982
983} // end namespace Eigen
984
985#endif // EIGEN_FUNCTORS_H
986