#ifndef DOTPROD_HPP
#define DOTPROD_HPP

#include <cassert>
#include <iterator>
#include "dim.hpp"

template<typename A, typename B>
auto dotprod(const A& a, const B& b) ->
      decltype(
	 dim(a), dim(b),
	 std::begin(a) + 3,
	 std::begin(b) + 3,
	 *std::begin(a) + *std::begin(a) * *std::begin(b)
      ) {
   using T = decltype(*std::begin(a) * *std::begin(b));
   auto len = dim(a);
   using Index = decltype(len);
   assert(len = dim(b));
   T sum{};
   for (Index i = 0; i < len; ++i) {
      sum += *(std::begin(a) + i) * *(std::begin(b) + i);
   }
   return sum;
}

#endif