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

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 TypeResult
  >
  multi_array< TypeResult, NumDimsLeft+NumDimsRight>
outer_product
  ( multi_array< TypeLeft, NumDimsLeft>const& a_left
  , multi_array< TypeRight, NumDimsRight>const& a_right
  , TypeResult (*op_binary)( TypeLeft, TypeRight)
  )
  /**@brief
   *  This code is based on that in section:
   *    4.5 Outer Product
   *  from book:
   *    _An APL Compiler_
   *    Timothy Budd
   *    Springer-Verlag 1988
   */
  {
      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< TypeResult, 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 size_right=a_right.num_elements();
      size_type space_result=a_result.num_elements();
      for( size_type offset_result=0; offset_result<space_result; ++offset_result)
      {
          auto quot_rem=div( offset_result, size_right);
          data_result[offset_result]
            =op_binary( data_left[quot_rem.quot], data_right[quot_rem.rem]);
      }
      return a_result;
  }

typedef long long elem_t;

  std::pair<elem_t,elem_t>
mk_elem_t_pair(elem_t a_left, elem_t a_right)
{
    std::pair< elem_t, elem_t> a_pair(a_left, a_right);
    return a_pair;
}    


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, mk_elem_t_pair);
        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;
}