#ifndef BOOST_MULTI_ARRAY_NUMERIC_INNER_PRODUCT_LE060811_HPP
#define BOOST_MULTI_ARRAY_NUMERIC_INNER_PRODUCT_LE060811_HPP
#include <boost/multi_array.hpp>
#include <boost/assert.hpp>
#include <functional> //for std::multiplies, std::plus
#include <algorithm>
//#define MULTI_ARRAY_NUMERIC_INNER_PRODUCT_TRACE

namespace boost {
  template
  < typename BinaryFunctor
  >
  struct 
identity
  /**@brief
   *  static value returns identity value
   *  for functor, BinaryFunctor.
   */
;
  template
  < typename Operand
  >
  struct
identity
  < std::plus<Operand>
  >
{
        static
      Operand const
    value=Operand(0)
    ;  
};

  template
  < typename Operand
  >
  struct
identity
  < std::multiplies<Operand>
  >
{
        static
      Operand const
    value=Operand(1)
    ;  
};

  template
  < typename T
  , std::size_t NumDimsLeft
  , std::size_t NumDimsRight
  , typename OpOuter=std::plus<T>
  , typename OpInner=std::multiplies<T>
  >
  void
inner_product
  ( multi_array< T, NumDimsLeft>const& a_left
  , multi_array< T, NumDimsRight>const& a_right
  , multi_array< T, NumDimsLeft-1+NumDimsRight-1>& a_result
  , OpOuter op_outer=std::plus<T>()
    //used to calculate fold of results of op_inner.
  , OpInner op_inner=std::multiplies<T>()
    //Used to calculate op_inner(left_element,right_element)
    //where:
    //  left_element is some element from a_left
    //  right_element is some element from a_right
    //The result is then fed to op_outer.
  )
  /**@brief
   *  Calculates inner prodduct of a_left and a_right and puts
   *  result in a_result.
   **@reference:
   *  The method used is, essentially, that outlined
   *  in the pseudo-code in section 5.7 of:
   *    _An Apl Compiler_
   *    Timothy Budd
   *    Springer-Verlag 1988
   */
{
    auto left_shape=a_left.shape();
    auto length_common=left_shape[NumDimsLeft-1]
      //L in "Shape Phase" and "Value Phase"
      //code of @reference.
      ;
    { //Check arguments for proper shape.
      //
      //The following assertions correspond to the
      //following passage from @reference:
      //
      //  let L_1,L_2,...,L_n represent the
      //  size of the left argument, and
      //  R_1,R_2,...,R_m represent the size
      //  of the right argument.  It must be
      //  true that L_n == R_1.  The shape of
      //  the result is given by...
      //    L_1,L_2,...,L_n-1,R_2,...,R_m
      //
        auto right_shape=a_right.shape();
        BOOST_ASSERT(length_common == right_shape[0])
          //L_n == R_1
          ;
        auto result_shape=a_result.shape();
        bool eq_front_shapes =
          std::equal
          ( left_shape
          , left_shape+NumDimsLeft-2
          , result_shape
          );
        bool eq_back_shapes =
          std::equal
          ( right_shape+1
          , right_shape+NumDimsRight-1
          , result_shape+NumDimsLeft-1
          );
        BOOST_ASSERT( eq_front_shapes && eq_back_shapes)
          //the shape of the result is given by...
          //    L_1,L_2,...,L_n-1,R_2,...,R_m
          ;
    }
    
  //The "Shape Phase" code of @reference.
    auto stride_n_left=a_left.strides()[NumDimsLeft-1]
      //not in @reference.  Used to account for 
      //non-c_storage_order. (with c_storage_order, ==1).
      ;
    auto stride_1_right
      =a_right.strides()[0]
      //e_1 of @reference *when* used on rhs of offset_right.
      ;
    multi_array_types::size_type const space_1_right
      =a_right.num_elements()/length_common
      //e_1 of reference *when* used with div and mod.
      //Different from @reference to account for
      //non-c_storage_order.  This is
      //"the size of the right subexpression 
      //without the first dimension".
      ;
    auto data_left=a_left.data();
    auto data_right=a_right.data();
    auto data_result=a_result.data();
    
    auto size_result=a_result.num_elements();
    for
    ( unsigned offset_result=0//offset of @reference
    ; offset_result<size_result
    ; ++offset_result
    )
    { //the body of this loop is the "Value Phase" code of @reference.
        auto quot_rem=div(offset_result,space_1_right);
        int offset_left//offset'_L in @reference.
          = quot_rem.quot//offset div e_1 in @reference
          * length_common
          + (length_common-1)
          * stride_n_left //1 in @reference. Difference to account for
                          //non-c_storage_order.
          ;
        int offset_right//offset'_R in @reference.
          = quot_rem.rem//offset mod e_1 in @reference.
          + (length_common-1)
          * stride_1_right//e_1 in @reference
          ;
        typedef T elem_t;
        elem_t elem_result//called 'result' in @reference.
          =identity<OpOuter>::value//identity for op_outer functor.
          ;
      #ifdef MULTI_ARRAY_NUMERIC_INNER_PRODUCT_TRACE
        std::cout<<":offset_result="<<offset_result<<"\n";
        std::cout<<":elem_result="<<elem_result<<"\n";
        std::cout<<indent_buf_in;
      #endif
        for(int counter=length_common-1; 0<=counter; --counter)
        {
            elem_t elem_left//called 'value_L' in @reference.
              =data_left[offset_left];
            elem_t elem_right//called 'value_R' in @reference.
              =data_right[offset_right];
            elem_t elem_inner=op_inner(elem_left,elem_right);
            elem_result=op_outer(elem_inner,elem_result);
          #ifdef MULTI_ARRAY_NUMERIC_INNER_PRODUCT_TRACE
            std::cout<<":offset_left="<<offset_left<<"\n";
            std::cout<<":offset_right="<<offset_right<<"\n";
            std::cout<<":elem_left="<<elem_left<<"\n";
            std::cout<<":elem_right="<<elem_right<<"\n";
            std::cout<<":elem_result="<<elem_result<<"\n";
          #endif
            offset_left -=stride_n_left 
              //change from @reference 
              //(account non-c_storage_order).
              ;
            offset_right-=stride_1_right
              //same as @reference since 
              //stride_1_right corresponds to e_1 in @reference
              ;
        }
      #ifdef MULTI_ARRAY_NUMERIC_INNER_PRODUCT_TRACE
        std::cout<<indent_buf_out;
      #endif
        data_result[offset_result]=elem_result;
    }
}

} // namespace boost
#endif