X-Git-Url: https://git.dogcows.com/gitweb?p=chaz%2Fyoink;a=blobdiff_plain;f=src%2Fcml%2Fvector%2Fvector_products.h;fp=src%2Fcml%2Fvector%2Fvector_products.h;h=89f5edcf66530fae11d2ac2f5fd1e504e743b05a;hp=0000000000000000000000000000000000000000;hb=0fffd0097d7b496454413e57b398c903ecc252e4;hpb=79becf045222f385da5a1b9eb79081f6f5266c86 diff --git a/src/cml/vector/vector_products.h b/src/cml/vector/vector_products.h new file mode 100644 index 0000000..89f5edc --- /dev/null +++ b/src/cml/vector/vector_products.h @@ -0,0 +1,361 @@ +/* -*- 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 +#include +#include +#include +#include +#include + +/* 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 +struct DotPromote +{ + /* Shorthand: */ + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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 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 +struct CrossPromote +{ + /* Shorthand: */ + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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 +struct OuterPromote +{ + /* Shorthand: */ + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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 +inline typename DotPromote::promoted_scalar +UnrollDot(const LeftT& left, const RightT& right, fixed_size_tag) +{ + /* Shorthand: */ + typedef DotPromote dot_helper; + + /* Compile-type vector size check: */ + typedef typename et::GetCheckedSize + ::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 +inline typename DotPromote::promoted_scalar +UnrollDot(const LeftT& left, const RightT& right, dynamic_size_tag) +{ + /* Shorthand: */ + typedef DotPromote dot_helper; + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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 + * = 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 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 inline void +Require3D(const VecT& v, dynamic_size_tag) { + et::GetCheckedSize() + .equal_or_fail(v.size(),size_t(3)); +} + +/** For perp_dot(): compile-time check for a 2D vector. */ +template 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 inline void +Require2D(const VecT& v, dynamic_size_tag) { + et::GetCheckedSize() + .equal_or_fail(v.size(),size_t(2)); +} + +} // namespace detail + + +/** Vector dot (inner) product implementation. + */ +template +inline typename detail::DotPromote::promoted_scalar +dot(const LeftT& left, const RightT& right) +{ + /* Shorthand: */ + typedef detail::DotPromote dot_helper; + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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::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 +inline typename detail::DotPromote::promoted_scalar +perp_dot(const LeftT& left, const RightT& right) +{ + /* Shorthand: */ + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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::is_true + && same_type::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 +inline typename detail::CrossPromote::promoted_vector +cross(const LeftT& left, const RightT& right) +{ + /* Shorthand: */ + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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::is_true + && same_type::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 +inline typename detail::OuterPromote::promoted_matrix +outer(const LeftT& left, const RightT& right) +{ + /* Shorthand: */ + typedef et::ExprTraits left_traits; + typedef et::ExprTraits 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::is_true + && same_type::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::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