1 /* -*- C++ -*- ------------------------------------------------------------
3 Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
5 The Configurable Math Library (CML) is distributed under the terms of the
6 Boost Software License, v1.0 (see cml/LICENSE for details).
8 *-----------------------------------------------------------------------*/
10 * @brief Defines vector dot and outer products.
12 * @todo Figure out if the source or destination size type should trigger
13 * unrolling. May need a per-compiler compile-time option for this.
16 #ifndef vector_products_h
17 #define vector_products_h
19 #include <cml/core/cml_assert.h>
20 #include <cml/et/scalar_promotions.h>
21 #include <cml/et/size_checking.h>
22 #include <cml/vector/vector_unroller.h>
23 #include <cml/vector/vector_expr.h>
24 #include <cml/matrix/matrix_expr.h>
26 /* This is used below to create a more meaningful compile-time error when
27 * dot() is not provided with vector or VectorExpr arguments:
29 struct dot_expects_vector_args_error
;
31 /* This is used below to create a more meaningful compile-time error when
32 * perp_dot() is not provided with 2D vector or VectorExpr arguments:
34 struct perp_dot_expects_vector_args_error
;
35 struct perp_dot_expects_2D_vector_args_error
;
37 /* This is used below to create a more meaningful compile-time error when
38 * outer() is not provided with vector or VectorExpr arguments:
40 struct outer_expects_vector_args_error
;
42 /* This is used below to create a more meaningful compile-time error when
43 * cross() is not provided with 3D vector or VectorExpr arguments:
45 struct cross_expects_vector_args_error
;
46 struct cross_expects_3D_vector_args_error
;
52 template<typename LeftT
, typename RightT
>
56 typedef et::ExprTraits
<LeftT
> left_traits
;
57 typedef et::ExprTraits
<RightT
> right_traits
;
58 typedef typename
left_traits::value_type left_value
;
59 typedef typename
right_traits::value_type right_value
;
61 /* Deduce the promoted scalar type: */
62 typedef et::OpMul
<left_value
, right_value
> op_mul
;
63 typedef typename
et::OpAdd
<
64 typename
op_mul::value_type
,
65 typename
op_mul::value_type
> op_add
;
66 typedef typename
op_add::value_type promoted_scalar
;
69 template<typename LeftT
, typename RightT
>
73 typedef et::ExprTraits
<LeftT
> left_traits
;
74 typedef et::ExprTraits
<RightT
> right_traits
;
75 typedef typename
left_traits::result_type left_type
;
76 typedef typename
right_traits::result_type right_type
;
78 /* Deduce the matrix result type: */
79 typedef typename
et::VectorPromote
<
80 left_type
,right_type
>::temporary_type promoted_vector
;
83 template<typename LeftT
, typename RightT
>
87 typedef et::ExprTraits
<LeftT
> left_traits
;
88 typedef et::ExprTraits
<RightT
> right_traits
;
89 typedef typename
left_traits::result_type left_type
;
90 typedef typename
right_traits::result_type right_type
;
92 /* Deduce the matrix result type: */
93 typedef typename
et::MatrixPromote
<
94 left_type
,right_type
>::temporary_type promoted_matrix
;
97 /** Construct a dot unroller for fixed-size arrays.
99 * @note This should only be called for vectors.
103 template<typename LeftT
, typename RightT
>
104 inline typename DotPromote
<LeftT
,RightT
>::promoted_scalar
105 UnrollDot(const LeftT
& left
, const RightT
& right
, fixed_size_tag
)
108 typedef DotPromote
<LeftT
,RightT
> dot_helper
;
110 /* Compile-type vector size check: */
111 typedef typename
et::GetCheckedSize
<LeftT
,RightT
,fixed_size_tag
>
112 ::check_type check_sizes
;
114 /* Get the fixed array size using the helper: */
115 enum { Len
= check_sizes::array_size
};
117 /* Record the unroller type: */
118 typedef typename
dot_helper::op_mul op_mul
;
119 typedef typename
dot_helper::op_add op_add
;
120 typedef typename
et::detail::VectorAccumulateUnroller
<
121 op_add
,op_mul
,LeftT
,RightT
>::template
122 Eval
<0, Len
-1, (Len
<= CML_VECTOR_DOT_UNROLL_LIMIT
)> Unroller
;
123 /* Note: Len is the array size, so Len-1 is the last element. */
125 /* Now, call the unroller: */
126 return Unroller()(left
,right
);
129 /** Use a loop to compute the dot product for dynamic arrays.
131 * @note This should only be called for vectors.
135 template<typename LeftT
, typename RightT
>
136 inline typename DotPromote
<LeftT
,RightT
>::promoted_scalar
137 UnrollDot(const LeftT
& left
, const RightT
& right
, dynamic_size_tag
)
140 typedef DotPromote
<LeftT
,RightT
> dot_helper
;
141 typedef et::ExprTraits
<LeftT
> left_traits
;
142 typedef et::ExprTraits
<RightT
> right_traits
;
143 typedef typename
dot_helper::op_mul op_mul
;
144 typedef typename
dot_helper::op_add op_add
;
146 /* Record the return type: */
147 typedef typename
dot_helper::promoted_scalar sum_type
;
149 /* Verify expression sizes: */
150 const size_t N
= et::CheckedSize(left
,right
,dynamic_size_tag());
152 /* Initialize the sum. Left and right must be vector expressions, so
153 * it's okay to use array notation here:
155 sum_type
sum(op_mul().apply(left
[0],right
[0]));
156 for(size_t i
= 1; i
< N
; ++i
) {
157 /* XXX This might not be optimized properly by some compilers.
158 * but to do anything else requires changing the requirements
159 * of a scalar operator, or requires defining a new class of scalar
162 sum
= op_add().apply(sum
, op_mul().apply(left
[i
], right
[i
]));
163 /* Note: we don't need get(), since both arguments are required to
164 * be vector expressions.
170 /** For cross(): compile-time check for a 3D vector. */
171 template<typename VecT
> inline void
172 Require3D(const VecT
&, fixed_size_tag
) {
173 CML_STATIC_REQUIRE_M(
174 ((size_t)VecT::array_size
== 3),
175 cross_expects_3D_vector_args_error
);
178 /** For cross(): run-time check for a 3D vector. */
179 template<typename VecT
> inline void
180 Require3D(const VecT
& v
, dynamic_size_tag
) {
181 et::GetCheckedSize
<VecT
,VecT
,dynamic_size_tag
>()
182 .equal_or_fail(v
.size(),size_t(3));
185 /** For perp_dot(): compile-time check for a 2D vector. */
186 template<typename VecT
> inline void
187 Require2D(const VecT
& v
, fixed_size_tag
) {
188 CML_STATIC_REQUIRE_M(
189 ((size_t)VecT::array_size
== 2),
190 perp_dot_expects_2D_vector_args_error
);
193 /** For perp_dot(): run-time check for a 2D vector. */
194 template<typename VecT
> inline void
195 Require2D(const VecT
& v
, dynamic_size_tag
) {
196 et::GetCheckedSize
<VecT
,VecT
,dynamic_size_tag
>()
197 .equal_or_fail(v
.size(),size_t(2));
200 } // namespace detail
203 /** Vector dot (inner) product implementation.
205 template<typename LeftT
, typename RightT
>
206 inline typename
detail::DotPromote
<LeftT
,RightT
>::promoted_scalar
207 dot(const LeftT
& left
, const RightT
& right
)
210 typedef detail::DotPromote
<LeftT
,RightT
> dot_helper
;
211 typedef et::ExprTraits
<LeftT
> left_traits
;
212 typedef et::ExprTraits
<RightT
> right_traits
;
213 typedef typename
left_traits::result_type left_type
;
214 typedef typename
right_traits::result_type right_type
;
215 typedef typename
left_traits::size_tag left_size
;
216 typedef typename
right_traits::size_tag right_size
;
218 /* dot() requires vector expressions: */
219 CML_STATIC_REQUIRE_M(
220 (et::VectorExpressions
<LeftT
,RightT
>::is_true
),
221 dot_expects_vector_args_error
);
222 /* Note: parens are required here so that the preprocessor ignores the
226 /* Figure out the unroller to use (fixed or dynamic): */
227 typedef typename
et::VectorPromote
<
228 left_type
, right_type
>::temporary_type promoted_vector
;
229 typedef typename
promoted_vector::size_tag size_tag
;
232 return detail::UnrollDot(left
,right
,size_tag());
237 template<typename LeftT
, typename RightT
>
238 inline typename
detail::DotPromote
<LeftT
,RightT
>::promoted_scalar
239 perp_dot(const LeftT
& left
, const RightT
& right
)
242 typedef et::ExprTraits
<LeftT
> left_traits
;
243 typedef et::ExprTraits
<RightT
> right_traits
;
244 typedef typename
left_traits::result_tag left_result
;
245 typedef typename
right_traits::result_tag right_result
;
247 /* perp_dot() requires vector expressions: */
248 CML_STATIC_REQUIRE_M(
249 (same_type
<left_result
, et::vector_result_tag
>::is_true
250 && same_type
<right_result
, et::vector_result_tag
>::is_true
),
251 perp_dot_expects_vector_args_error
);
252 /* Note: parens are required here so that the preprocessor ignores the
256 /* Make sure arguments are 2D vectors: */
257 detail::Require2D(left
, typename
left_traits::size_tag());
258 detail::Require2D(right
, typename
right_traits::size_tag());
260 /* Get result type: */
261 typedef typename
detail::DotPromote
<
262 LeftT
,RightT
>::promoted_scalar result_type
;
264 /* Compute and return: */
265 return result_type(left
[0]*right
[1]-left
[1]*right
[0]);
268 template<typename LeftT
, typename RightT
>
269 inline typename
detail::CrossPromote
<LeftT
,RightT
>::promoted_vector
270 cross(const LeftT
& left
, const RightT
& right
)
273 typedef et::ExprTraits
<LeftT
> left_traits
;
274 typedef et::ExprTraits
<RightT
> right_traits
;
275 typedef typename
left_traits::result_tag left_result
;
276 typedef typename
right_traits::result_tag right_result
;
278 /* outer() requires vector expressions: */
279 CML_STATIC_REQUIRE_M(
280 (same_type
<left_result
, et::vector_result_tag
>::is_true
281 && same_type
<right_result
, et::vector_result_tag
>::is_true
),
282 cross_expects_vector_args_error
);
283 /* Note: parens are required here so that the preprocessor ignores the
287 /* Make sure arguments are 3D vectors: */
288 detail::Require3D(left
, typename
left_traits::size_tag());
289 detail::Require3D(right
, typename
right_traits::size_tag());
291 /* Get result type: */
292 typedef typename
detail::CrossPromote
<
293 LeftT
,RightT
>::promoted_vector result_type
;
295 /* Now, compute and return the cross product: */
297 left
[1]*right
[2] - left
[2]*right
[1],
298 left
[2]*right
[0] - left
[0]*right
[2],
299 left
[0]*right
[1] - left
[1]*right
[0]
304 /** Return the triple product of three 3D vectors.
306 * No checking is done here, as dot() and cross() will catch any size or
310 template < class VecT_1
, class VecT_2
, class VecT_3
>
311 typename
detail::DotPromote
<
312 VecT_1
, typename
detail::CrossPromote
< VecT_2
, VecT_3
>::promoted_vector
314 triple_product(const VecT_1
& v1
, const VecT_2
& v2
, const VecT_3
& v3
) {
315 return dot(v1
,cross(v2
,v3
));
318 template<typename LeftT
, typename RightT
>
319 inline typename
detail::OuterPromote
<LeftT
,RightT
>::promoted_matrix
320 outer(const LeftT
& left
, const RightT
& right
)
323 typedef et::ExprTraits
<LeftT
> left_traits
;
324 typedef et::ExprTraits
<RightT
> right_traits
;
325 typedef typename
left_traits::result_tag left_result
;
326 typedef typename
right_traits::result_tag right_result
;
328 /* outer() requires vector expressions: */
329 CML_STATIC_REQUIRE_M(
330 (same_type
<left_result
, et::vector_result_tag
>::is_true
331 && same_type
<right_result
, et::vector_result_tag
>::is_true
),
332 dot_expects_vector_args_error
);
333 /* Note: parens are required here so that the preprocessor ignores the
337 /* Create a matrix with the right size (resize() is a no-op for
338 * fixed-size matrices):
340 typename
detail::OuterPromote
<LeftT
,RightT
>::promoted_matrix C
;
341 cml::et::detail::Resize(C
, left
.size(), right
.size());
343 /* Now, compute the outer product: */
344 for(size_t i
= 0; i
< left
.size(); ++i
) {
345 for(size_t j
= 0; j
< right
.size(); ++j
) {
346 C(i
,j
) = left
[i
]*right
[j
];
347 /* Note: both arguments are vectors, so array notation
360 // -------------------------------------------------------------------------