#include <boost/multi_array/io/print.hpp>
#include <boost/iterator/counting_iterator.hpp>
#include <vector>
#include <utility>
#include <functional>

  template
  < typename Left
  , typename Right
  >
  struct
make_value_pair
  {
          typedef
        std::pair< Left, Right>
      result_type
      ;
        result_type
      operator()(Left a_left, Right a_right)
      {
          std::pair< Left, Right> a_pair(a_left, a_right);
          return a_pair;
      }
  };
  
using namespace boost;

  template< typename T>
  std::ostream&
operator<<
  ( std::ostream& sout
  , std::pair<T*,T*> ip
  )
  {
      sout<<"{ ";
      for( unsigned i=0; ip.first < ip.second; ++i,++ip.first)
      {
          if(0<i)sout<<", ";
          sout<<*ip.first;
      }
      sout<<"}\n";
      return sout;
  }
  
  template< typename Left, typename Right>
  std::ostream&
operator<<
  ( std::ostream& sout
  , std::pair<Left,Right> const& a_pair
  )
  {
      sout<<"pair( "<<a_pair.first<<", "<<a_pair.second<<")";
      return sout;
  }
  
  template
  < typename TypeLeft
  , std::size_t NumDimsLeft
  , typename TypeRight
  , std::size_t NumDimsRight
  , typename OpBinary
  >
  multi_array
  < typename std::result_of<OpBinary( TypeLeft&&, TypeRight&&)>::type
  , NumDimsLeft+NumDimsRight
  >
outer_product
  ( multi_array< TypeLeft, NumDimsLeft>const& a_left
  , multi_array< TypeRight, NumDimsRight>const& a_right
  , OpBinary op_binary
  )
  /**@brief
   *  This code is based on that in section:
   *    4.5 Outer Product
   *  from book:
   *    _An APL Compiler_
   *    Timothy Budd
   *    Springer-Verlag 1988
   */
  {
      typedef typename std::result_of<OpBinary( TypeLeft, TypeRight)>::type type_result;
      typedef multi_array_types::size_type size_type;
      std::size_t const num_dims_result=NumDimsLeft+NumDimsRight;
      std::vector<size_type> shape_result( num_dims_result, 0);
      auto shape_left=a_left.shape();
      auto shape_right=a_right.shape();
      std::copy
        ( shape_right
        , shape_right+NumDimsRight
        , std::copy
          ( shape_left
          , shape_left+NumDimsLeft
          , shape_result.begin()
          )
        );
      multi_array< type_result, num_dims_result> a_result(shape_result);
      auto data_left=a_left.data();
      auto data_right=a_right.data();
      auto data_result=a_result.data();
      size_type space_left=a_left.num_elements();
      size_type space_right=a_right.num_elements();
      size_type offset_result=0;
      for( size_type offset_left=0; offset_left<space_left; ++offset_left)
      {
          for( size_type offset_right=0; offset_right<space_right; ++offset_right)
          {
              data_result[offset_result++]
                = op_binary
                  ( data_left [offset_left ]
                  , data_right[offset_right]
                  );
          }
      }
      return a_result;
  }

    typedef 
  long long 
elem_t
  ;

int main(void)
{
    iostreams::indent_scoped_ostreambuf<char> indent_outbuf(std::cout,2);
    std::size_t const dim_left=1;
    std::size_t const dim_right=1;
    typedef std::vector<multi_array_types::size_type> shape_t;
    shape_t shape_left({2});
    shape_t shape_right({2});
    {
        multi_array<elem_t,dim_left> a_left( shape_left);
        elem_t const left0=100;
        a_left.assign
          ( counting_iterator<elem_t>(left0)
          , counting_iterator<elem_t>(left0+a_left.num_elements())
          );
        std::cout<<":a_left=\n"<<a_left<<"\n";
        std::cout<<":a_left.data=\n"
          <<std::make_pair
            ( a_left.data()
            , a_left.data()+a_left.num_elements()
            )
          <<"\n";
        multi_array<elem_t,dim_right> a_right( shape_right);
        elem_t const right0=1000;
        a_right.assign
          ( counting_iterator<elem_t>(right0)
          , counting_iterator<elem_t>(right0+a_left.num_elements())
          );
        std::cout<<":a_right=\n"<<a_right<<"\n";
        std::cout<<":a_right.data=\n"
          <<std::make_pair
            ( a_right.data()
            , a_right.data()+a_right.num_elements()
            )
          <<"\n";
        auto a_result
          = outer_product
            ( a_left
            , a_right
            , make_value_pair<elem_t,elem_t>()
            );
        std::cout<<":a_result=\n"<<a_result<<"\n";
        std::cout<<":a_result.data=\n"
          <<std::make_pair
            ( a_result.data()
            , a_result.data()+a_result.num_elements()
            )
          <<"\n";
    }
    return 0;
}