Dogcows Code - chaz/yoink/blob - src/Moof/cml/vector/vector_products.h

   1 /* -*- C++ -*- ------------------------------------------------------------
   2
   3 Copyright (c) 2007 Jesse Anders and Demian Nave http://cmldev.net/
   4
   5 The Configurable Math Library (CML) is distributed under the terms of the
   6 Boost Software License, v1.0 (see cml/LICENSE for details).
   7
   8  *-----------------------------------------------------------------------*/
   9 /** @file
  10  *  @brief Defines vector dot and outer products.
  11  *
  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.
  14  */
  15
  16 #ifndef vector_products_h
  17 #define vector_products_h
  18
  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>
  25
  26 /* This is used below to create a more meaningful compile-time error when
  27  * dot() is not provided with vector or VectorExpr arguments:
  28  */
  29 struct dot_expects_vector_args_error;
  30
  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:
  33  */
  34 struct perp_dot_expects_vector_args_error;
  35 struct perp_dot_expects_2D_vector_args_error;
  36
  37 /* This is used below to create a more meaningful compile-time error when
  38  * outer() is not provided with vector or VectorExpr arguments:
  39  */
  40 struct outer_expects_vector_args_error;
  41
  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:
  44  */
  45 struct cross_expects_vector_args_error;
  46 struct cross_expects_3D_vector_args_error;
  47
  48
  49 namespace cml {
  50 namespace detail {
  51
  52 template<typename LeftT, typename RightT>
  53 struct DotPromote
  54 {
  55     /* Shorthand: */
  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;
  60
  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;
  67 };
  68
  69 template<typename LeftT, typename RightT>
  70 struct CrossPromote
  71 {
  72     /* Shorthand: */
  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;
  77
  78     /* Deduce the matrix result type: */
  79     typedef typename et::VectorPromote<
  80         left_type,right_type>::temporary_type promoted_vector;
  81 };
  82
  83 template<typename LeftT, typename RightT>
  84 struct OuterPromote
  85 {
  86     /* Shorthand: */
  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;
  91
  92     /* Deduce the matrix result type: */
  93     typedef typename et::MatrixPromote<
  94         left_type,right_type>::temporary_type promoted_matrix;
  95 };
  96
  97 /** Construct a dot unroller for fixed-size arrays.
  98  *
  99  * @note This should only be called for vectors.
 100  *
 101  * @sa cml::dot
 102  */
 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)
 106 {
 107     /* Shorthand: */
 108     typedef DotPromote<LeftT,RightT> dot_helper;
 109
 110     /* Compile-type vector size check: */
 111     typedef typename et::GetCheckedSize<LeftT,RightT,fixed_size_tag>
 112         ::check_type check_sizes;
 113
 114     /* Get the fixed array size using the helper: */
 115     enum { Len = check_sizes::array_size };
 116
 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. */
 124
 125     /* Now, call the unroller: */
 126     return Unroller()(left,right);
 127 }
 128
 129 /** Use a loop to compute the dot product for dynamic arrays.
 130  *
 131  * @note This should only be called for vectors.
 132  *
 133  * @sa cml::dot
 134  */
 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)
 138 {
 139     /* Shorthand: */
 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;
 145
 146     /* Record the return type: */
 147     typedef typename dot_helper::promoted_scalar sum_type;
 148
 149     /* Verify expression sizes: */
 150     const size_t N = et::CheckedSize(left,right,dynamic_size_tag());
 151
 152     /* Initialize the sum. Left and right must be vector expressions, so
 153      * it's okay to use array notation here:
 154      */
 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
 160          * <op>= operators.
 161          */
 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.
 165          */
 166     }
 167     return sum;
 168 }
 169
 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);
 176 }
 177
 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));
 183 }
 184
 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);
 191 }
 192
 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));
 198 }
 199
 200 } // namespace detail
 201
 202
 203 /** Vector dot (inner) product implementation.
 204  */
 205 template<typename LeftT, typename RightT>
 206 inline typename detail::DotPromote<LeftT,RightT>::promoted_scalar
 207 dot(const LeftT& left, const RightT& right)
 208 {
 209     /* Shorthand: */
 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;
 217
 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
 223      * commas:
 224      */
 225
 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;
 230
 231     /* Call unroller: */
 232     return detail::UnrollDot(left,right,size_tag());
 233 }
 234
 235 /** perp_dot()
 236  */
 237 template<typename LeftT, typename RightT>
 238 inline typename detail::DotPromote<LeftT,RightT>::promoted_scalar
 239 perp_dot(const LeftT& left, const RightT& right)
 240 {
 241     /* Shorthand: */
 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;
 246
 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
 253      * commas.
 254      */
 255
 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());
 259
 260     /* Get result type: */
 261     typedef typename detail::DotPromote<
 262         LeftT,RightT>::promoted_scalar result_type;
 263
 264     /* Compute and return: */
 265     return result_type(left[0]*right[1]-left[1]*right[0]);
 266 }
 267
 268 template<typename LeftT, typename RightT>
 269 inline typename detail::CrossPromote<LeftT,RightT>::promoted_vector
 270 cross(const LeftT& left, const RightT& right)
 271 {
 272     /* Shorthand: */
 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;
 277
 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
 284      * commas.
 285      */
 286
 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());
 290
 291     /* Get result type: */
 292     typedef typename detail::CrossPromote<
 293         LeftT,RightT>::promoted_vector result_type;
 294
 295     /* Now, compute and return the cross product: */
 296     result_type result(
 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]
 300             );
 301     return result;
 302 }
 303
 304 /** Return the triple product of three 3D vectors.
 305  *
 306  * No checking is done here, as dot() and cross() will catch any size or
 307  * type errors.
 308  */
 309
 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
 313 >::promoted_scalar
 314 triple_product(const VecT_1& v1, const VecT_2& v2, const VecT_3& v3) {
 315     return dot(v1,cross(v2,v3));
 316 }
 317
 318 template<typename LeftT, typename RightT>
 319 inline typename detail::OuterPromote<LeftT,RightT>::promoted_matrix
 320 outer(const LeftT& left, const RightT& right)
 321 {
 322     /* Shorthand: */
 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;
 327
 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
 334      * commas.
 335      */
 336
 337     /* Create a matrix with the right size (resize() is a no-op for
 338      * fixed-size matrices):
 339      */
 340     typename detail::OuterPromote<LeftT,RightT>::promoted_matrix C;
 341     cml::et::detail::Resize(C, left.size(), right.size());
 342
 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
 348              * is okay here.
 349              */
 350         }
 351     }
 352
 353     return C;
 354 }
 355
 356 } // namespace cml
 357
 358 #endif
 359
 360 // -------------------------------------------------------------------------
 361 // vim:ft=cpp