--- /dev/null
+/* -*- C++ -*- ------------------------------------------------------------
+
+Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
+
+The Configurable Math Library (CML) is distributed under the terms of the
+Boost Software License, v1.0 (see cml/LICENSE for details).
+
+ *-----------------------------------------------------------------------*/
+/** @file
+ * @brief Defines vector dot and outer products.
+ *
+ * @todo Figure out if the source or destination size type should trigger
+ * unrolling. May need a per-compiler compile-time option for this.
+ */
+
+#ifndef vector_products_h
+#define vector_products_h
+
+#include <cml/core/cml_assert.h>
+#include <cml/et/scalar_promotions.h>
+#include <cml/et/size_checking.h>
+#include <cml/vector/vector_unroller.h>
+#include <cml/vector/vector_expr.h>
+#include <cml/matrix/matrix_expr.h>
+
+/* This is used below to create a more meaningful compile-time error when
+ * dot() is not provided with vector or VectorExpr arguments:
+ */
+struct dot_expects_vector_args_error;
+
+/* This is used below to create a more meaningful compile-time error when
+ * perp_dot() is not provided with 2D vector or VectorExpr arguments:
+ */
+struct perp_dot_expects_vector_args_error;
+struct perp_dot_expects_2D_vector_args_error;
+
+/* This is used below to create a more meaningful compile-time error when
+ * outer() is not provided with vector or VectorExpr arguments:
+ */
+struct outer_expects_vector_args_error;
+
+/* This is used below to create a more meaningful compile-time error when
+ * cross() is not provided with 3D vector or VectorExpr arguments:
+ */
+struct cross_expects_vector_args_error;
+struct cross_expects_3D_vector_args_error;
+
+
+namespace cml {
+namespace detail {
+
+template<typename LeftT, typename RightT>
+struct DotPromote
+{
+ /* Shorthand: */
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::value_type left_value;
+ typedef typename right_traits::value_type right_value;
+
+ /* Deduce the promoted scalar type: */
+ typedef et::OpMul<left_value, right_value> op_mul;
+ typedef typename et::OpAdd<
+ typename op_mul::value_type,
+ typename op_mul::value_type> op_add;
+ typedef typename op_add::value_type promoted_scalar;
+};
+
+template<typename LeftT, typename RightT>
+struct CrossPromote
+{
+ /* Shorthand: */
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::result_type left_type;
+ typedef typename right_traits::result_type right_type;
+
+ /* Deduce the matrix result type: */
+ typedef typename et::VectorPromote<
+ left_type,right_type>::temporary_type promoted_vector;
+};
+
+template<typename LeftT, typename RightT>
+struct OuterPromote
+{
+ /* Shorthand: */
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::result_type left_type;
+ typedef typename right_traits::result_type right_type;
+
+ /* Deduce the matrix result type: */
+ typedef typename et::MatrixPromote<
+ left_type,right_type>::temporary_type promoted_matrix;
+};
+
+/** Construct a dot unroller for fixed-size arrays.
+ *
+ * @note This should only be called for vectors.
+ *
+ * @sa cml::dot
+ */
+template<typename LeftT, typename RightT>
+inline typename DotPromote<LeftT,RightT>::promoted_scalar
+UnrollDot(const LeftT& left, const RightT& right, fixed_size_tag)
+{
+ /* Shorthand: */
+ typedef DotPromote<LeftT,RightT> dot_helper;
+
+ /* Compile-type vector size check: */
+ typedef typename et::GetCheckedSize<LeftT,RightT,fixed_size_tag>
+ ::check_type check_sizes;
+
+ /* Get the fixed array size using the helper: */
+ enum { Len = check_sizes::array_size };
+
+ /* Record the unroller type: */
+ typedef typename dot_helper::op_mul op_mul;
+ typedef typename dot_helper::op_add op_add;
+ typedef typename et::detail::VectorAccumulateUnroller<
+ op_add,op_mul,LeftT,RightT>::template
+ Eval<0, Len-1, (Len <= CML_VECTOR_DOT_UNROLL_LIMIT)> Unroller;
+ /* Note: Len is the array size, so Len-1 is the last element. */
+
+ /* Now, call the unroller: */
+ return Unroller()(left,right);
+}
+
+/** Use a loop to compute the dot product for dynamic arrays.
+ *
+ * @note This should only be called for vectors.
+ *
+ * @sa cml::dot
+ */
+template<typename LeftT, typename RightT>
+inline typename DotPromote<LeftT,RightT>::promoted_scalar
+UnrollDot(const LeftT& left, const RightT& right, dynamic_size_tag)
+{
+ /* Shorthand: */
+ typedef DotPromote<LeftT,RightT> dot_helper;
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename dot_helper::op_mul op_mul;
+ typedef typename dot_helper::op_add op_add;
+
+ /* Record the return type: */
+ typedef typename dot_helper::promoted_scalar sum_type;
+
+ /* Verify expression sizes: */
+ const size_t N = et::CheckedSize(left,right,dynamic_size_tag());
+
+ /* Initialize the sum. Left and right must be vector expressions, so
+ * it's okay to use array notation here:
+ */
+ sum_type sum(op_mul().apply(left[0],right[0]));
+ for(size_t i = 1; i < N; ++i) {
+ /* XXX This might not be optimized properly by some compilers.
+ * but to do anything else requires changing the requirements
+ * of a scalar operator, or requires defining a new class of scalar
+ * <op>= operators.
+ */
+ sum = op_add().apply(sum, op_mul().apply(left[i], right[i]));
+ /* Note: we don't need get(), since both arguments are required to
+ * be vector expressions.
+ */
+ }
+ return sum;
+}
+
+/** For cross(): compile-time check for a 3D vector. */
+template<typename VecT> inline void
+Require3D(const VecT&, fixed_size_tag) {
+ CML_STATIC_REQUIRE_M(
+ ((size_t)VecT::array_size == 3),
+ cross_expects_3D_vector_args_error);
+}
+
+/** For cross(): run-time check for a 3D vector. */
+template<typename VecT> inline void
+Require3D(const VecT& v, dynamic_size_tag) {
+ et::GetCheckedSize<VecT,VecT,dynamic_size_tag>()
+ .equal_or_fail(v.size(),size_t(3));
+}
+
+/** For perp_dot(): compile-time check for a 2D vector. */
+template<typename VecT> inline void
+Require2D(const VecT& v, fixed_size_tag) {
+ CML_STATIC_REQUIRE_M(
+ ((size_t)VecT::array_size == 2),
+ perp_dot_expects_2D_vector_args_error);
+}
+
+/** For perp_dot(): run-time check for a 2D vector. */
+template<typename VecT> inline void
+Require2D(const VecT& v, dynamic_size_tag) {
+ et::GetCheckedSize<VecT,VecT,dynamic_size_tag>()
+ .equal_or_fail(v.size(),size_t(2));
+}
+
+} // namespace detail
+
+
+/** Vector dot (inner) product implementation.
+ */
+template<typename LeftT, typename RightT>
+inline typename detail::DotPromote<LeftT,RightT>::promoted_scalar
+dot(const LeftT& left, const RightT& right)
+{
+ /* Shorthand: */
+ typedef detail::DotPromote<LeftT,RightT> dot_helper;
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::result_type left_type;
+ typedef typename right_traits::result_type right_type;
+ typedef typename left_traits::size_tag left_size;
+ typedef typename right_traits::size_tag right_size;
+
+ /* dot() requires vector expressions: */
+ CML_STATIC_REQUIRE_M(
+ (et::VectorExpressions<LeftT,RightT>::is_true),
+ dot_expects_vector_args_error);
+ /* Note: parens are required here so that the preprocessor ignores the
+ * commas:
+ */
+
+ /* Figure out the unroller to use (fixed or dynamic): */
+ typedef typename et::VectorPromote<
+ left_type, right_type>::temporary_type promoted_vector;
+ typedef typename promoted_vector::size_tag size_tag;
+
+ /* Call unroller: */
+ return detail::UnrollDot(left,right,size_tag());
+}
+
+/** perp_dot()
+ */
+template<typename LeftT, typename RightT>
+inline typename detail::DotPromote<LeftT,RightT>::promoted_scalar
+perp_dot(const LeftT& left, const RightT& right)
+{
+ /* Shorthand: */
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::result_tag left_result;
+ typedef typename right_traits::result_tag right_result;
+
+ /* perp_dot() requires vector expressions: */
+ CML_STATIC_REQUIRE_M(
+ (same_type<left_result, et::vector_result_tag>::is_true
+ && same_type<right_result, et::vector_result_tag>::is_true),
+ perp_dot_expects_vector_args_error);
+ /* Note: parens are required here so that the preprocessor ignores the
+ * commas.
+ */
+
+ /* Make sure arguments are 2D vectors: */
+ detail::Require2D(left, typename left_traits::size_tag());
+ detail::Require2D(right, typename right_traits::size_tag());
+
+ /* Get result type: */
+ typedef typename detail::DotPromote<
+ LeftT,RightT>::promoted_scalar result_type;
+
+ /* Compute and return: */
+ return result_type(left[0]*right[1]-left[1]*right[0]);
+}
+
+template<typename LeftT, typename RightT>
+inline typename detail::CrossPromote<LeftT,RightT>::promoted_vector
+cross(const LeftT& left, const RightT& right)
+{
+ /* Shorthand: */
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::result_tag left_result;
+ typedef typename right_traits::result_tag right_result;
+
+ /* outer() requires vector expressions: */
+ CML_STATIC_REQUIRE_M(
+ (same_type<left_result, et::vector_result_tag>::is_true
+ && same_type<right_result, et::vector_result_tag>::is_true),
+ cross_expects_vector_args_error);
+ /* Note: parens are required here so that the preprocessor ignores the
+ * commas.
+ */
+
+ /* Make sure arguments are 3D vectors: */
+ detail::Require3D(left, typename left_traits::size_tag());
+ detail::Require3D(right, typename right_traits::size_tag());
+
+ /* Get result type: */
+ typedef typename detail::CrossPromote<
+ LeftT,RightT>::promoted_vector result_type;
+
+ /* Now, compute and return the cross product: */
+ result_type result(
+ left[1]*right[2] - left[2]*right[1],
+ left[2]*right[0] - left[0]*right[2],
+ left[0]*right[1] - left[1]*right[0]
+ );
+ return result;
+}
+
+/** Return the triple product of three 3D vectors.
+ *
+ * No checking is done here, as dot() and cross() will catch any size or
+ * type errors.
+ */
+
+template < class VecT_1, class VecT_2, class VecT_3 >
+typename detail::DotPromote<
+ VecT_1, typename detail::CrossPromote< VecT_2, VecT_3 >::promoted_vector
+>::promoted_scalar
+triple_product(const VecT_1& v1, const VecT_2& v2, const VecT_3& v3) {
+ return dot(v1,cross(v2,v3));
+}
+
+template<typename LeftT, typename RightT>
+inline typename detail::OuterPromote<LeftT,RightT>::promoted_matrix
+outer(const LeftT& left, const RightT& right)
+{
+ /* Shorthand: */
+ typedef et::ExprTraits<LeftT> left_traits;
+ typedef et::ExprTraits<RightT> right_traits;
+ typedef typename left_traits::result_tag left_result;
+ typedef typename right_traits::result_tag right_result;
+
+ /* outer() requires vector expressions: */
+ CML_STATIC_REQUIRE_M(
+ (same_type<left_result, et::vector_result_tag>::is_true
+ && same_type<right_result, et::vector_result_tag>::is_true),
+ dot_expects_vector_args_error);
+ /* Note: parens are required here so that the preprocessor ignores the
+ * commas.
+ */
+
+ /* Create a matrix with the right size (resize() is a no-op for
+ * fixed-size matrices):
+ */
+ typename detail::OuterPromote<LeftT,RightT>::promoted_matrix C;
+ cml::et::detail::Resize(C, left.size(), right.size());
+
+ /* Now, compute the outer product: */
+ for(size_t i = 0; i < left.size(); ++i) {
+ for(size_t j = 0; j < right.size(); ++j) {
+ C(i,j) = left[i]*right[j];
+ /* Note: both arguments are vectors, so array notation
+ * is okay here.
+ */
+ }
+ }
+
+ return C;
+}
+
+} // namespace cml
+
+#endif
+
+// -------------------------------------------------------------------------
+// vim:ft=cpp
projects / chaz / yoink / blobdiff