#include <iostream>
#include <vector>
 
int main()
{
  // shape, stride & data of the matrix
 
  size_t shape  [] = { 2, 3, 4, 5};
  size_t strides[] = {60,20, 5, 1};
 
  std::vector<double> data(2*3*4*5);
 
  for ( size_t i = 0 ; i < data.size() ; ++i ) data[i] = 1.;
 
  // shape, stride & data (zero-initialized) of the reduced matrix
 
  size_t rshape  [] = { 2, 4, 5};
  size_t rstrides[] = {20, 5, 1};
 
  std::vector<double> rdata(2*4*5, 0.0);
 
  // compute reduction
 
//   for ( size_t a = 0 ; a < shape[0] ; ++a )
//     for ( size_t c = 0 ; c < shape[2] ; ++c )
//       for ( size_t d = 0 ; d < shape[3] ; ++d )
//         for ( size_t b = 0 ; b < shape[1] ; ++b )
//           rdata[ a*rstrides[0]                 + c*rstrides[1] + d*rstrides[2] ] += \
//           data [ a*strides [0] + b*strides [1] + c*strides [2] + d*strides [3] ];
 
  size_t cmp_axis = 1, axis_count = sizeof shape/ sizeof *shape;
  std::vector<size_t> adjusted_strides;
  //adjusted strides is basically same as strides
  //only difference being that the first element is the 
  //total number of elements in the n dim array.
 
  //The only reason to introduce this array was
  //so that I don't have to write any if-elses
  adjusted_strides.push_back(shape[0]*strides[0]);
  adjusted_strides.insert(adjusted_strides.end(), strides, strides + axis_count);
  for(size_t i = 0; i < data.size(); ++i) {
    size_t ni = i/adjusted_strides[cmp_axis]*adjusted_strides[cmp_axis+1] + i%adjusted_strides[cmp_axis+1];
    rdata[ni] += data[i];
  }
  // print resulting reduced matrix
 
  for ( size_t a = 0 ; a < rshape[0] ; ++a )
    for ( size_t b = 0 ; b < rshape[1] ; ++b )
      for ( size_t c = 0 ; c < rshape[2] ; ++c )
        std::cout << "(" << a << "," << b << "," << c << ") " << \
        rdata[ a*rstrides[0] + b*rstrides[1] + c*rstrides[2] ] << std::endl;
 
  return 0;
}

I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8dmVjdG9yPgoKaW50IG1haW4oKQp7CiAgLy8gc2hhcGUsIHN0cmlkZSAmIGRhdGEgb2YgdGhlIG1hdHJpeAoKICBzaXplX3Qgc2hhcGUgIFtdID0geyAyLCAzLCA0LCA1fTsKICBzaXplX3Qgc3RyaWRlc1tdID0gezYwLDIwLCA1LCAxfTsKCiAgc3RkOjp2ZWN0b3I8ZG91YmxlPiBkYXRhKDIqMyo0KjUpOwoKICBmb3IgKCBzaXplX3QgaSA9IDAgOyBpIDwgZGF0YS5zaXplKCkgOyArK2kgKSBkYXRhW2ldID0gMS47CgogIC8vIHNoYXBlLCBzdHJpZGUgJiBkYXRhICh6ZXJvLWluaXRpYWxpemVkKSBvZiB0aGUgcmVkdWNlZCBtYXRyaXgKCiAgc2l6ZV90IHJzaGFwZSAgW10gPSB7IDIsIDQsIDV9OwogIHNpemVfdCByc3RyaWRlc1tdID0gezIwLCA1LCAxfTsKCiAgc3RkOjp2ZWN0b3I8ZG91YmxlPiByZGF0YSgyKjQqNSwgMC4wKTsKCiAgLy8gY29tcHV0ZSByZWR1Y3Rpb24KCi8vICAgZm9yICggc2l6ZV90IGEgPSAwIDsgYSA8IHNoYXBlWzBdIDsgKythICkKLy8gICAgIGZvciAoIHNpemVfdCBjID0gMCA7IGMgPCBzaGFwZVsyXSA7ICsrYyApCi8vICAgICAgIGZvciAoIHNpemVfdCBkID0gMCA7IGQgPCBzaGFwZVszXSA7ICsrZCApCi8vICAgICAgICAgZm9yICggc2l6ZV90IGIgPSAwIDsgYiA8IHNoYXBlWzFdIDsgKytiICkKLy8gICAgICAgICAgIHJkYXRhWyBhKnJzdHJpZGVzWzBdICAgICAgICAgICAgICAgICArIGMqcnN0cmlkZXNbMV0gKyBkKnJzdHJpZGVzWzJdIF0gKz0gXAovLyAgICAgICAgICAgZGF0YSBbIGEqc3RyaWRlcyBbMF0gKyBiKnN0cmlkZXMgWzFdICsgYypzdHJpZGVzIFsyXSArIGQqc3RyaWRlcyBbM10gXTsKCiAgc2l6ZV90IGNtcF9heGlzID0gMSwgYXhpc19jb3VudCA9IHNpemVvZiBzaGFwZS8gc2l6ZW9mICpzaGFwZTsKICBzdGQ6OnZlY3RvcjxzaXplX3Q+IGFkanVzdGVkX3N0cmlkZXM7CiAgLy9hZGp1c3RlZCBzdHJpZGVzIGlzIGJhc2ljYWxseSBzYW1lIGFzIHN0cmlkZXMKICAvL29ubHkgZGlmZmVyZW5jZSBiZWluZyB0aGF0IHRoZSBmaXJzdCBlbGVtZW50IGlzIHRoZSAKICAvL3RvdGFsIG51bWJlciBvZiBlbGVtZW50cyBpbiB0aGUgbiBkaW0gYXJyYXkuCgogIC8vVGhlIG9ubHkgcmVhc29uIHRvIGludHJvZHVjZSB0aGlzIGFycmF5IHdhcwogIC8vc28gdGhhdCBJIGRvbid0IGhhdmUgdG8gd3JpdGUgYW55IGlmLWVsc2VzCiAgYWRqdXN0ZWRfc3RyaWRlcy5wdXNoX2JhY2soc2hhcGVbMF0qc3RyaWRlc1swXSk7CiAgYWRqdXN0ZWRfc3RyaWRlcy5pbnNlcnQoYWRqdXN0ZWRfc3RyaWRlcy5lbmQoKSwgc3RyaWRlcywgc3RyaWRlcyArIGF4aXNfY291bnQpOwogIGZvcihzaXplX3QgaSA9IDA7IGkgPCBkYXRhLnNpemUoKTsgKytpKSB7CiAgICBzaXplX3QgbmkgPSBpL2FkanVzdGVkX3N0cmlkZXNbY21wX2F4aXNdKmFkanVzdGVkX3N0cmlkZXNbY21wX2F4aXMrMV0gKyBpJWFkanVzdGVkX3N0cmlkZXNbY21wX2F4aXMrMV07CiAgICByZGF0YVtuaV0gKz0gZGF0YVtpXTsKICB9CiAgLy8gcHJpbnQgcmVzdWx0aW5nIHJlZHVjZWQgbWF0cml4CgogIGZvciAoIHNpemVfdCBhID0gMCA7IGEgPCByc2hhcGVbMF0gOyArK2EgKQogICAgZm9yICggc2l6ZV90IGIgPSAwIDsgYiA8IHJzaGFwZVsxXSA7ICsrYiApCiAgICAgIGZvciAoIHNpemVfdCBjID0gMCA7IGMgPCByc2hhcGVbMl0gOyArK2MgKQogICAgICAgIHN0ZDo6Y291dCA8PCAiKCIgPDwgYSA8PCAiLCIgPDwgYiA8PCAiLCIgPDwgYyA8PCAiKSAiIDw8IFwKICAgICAgICByZGF0YVsgYSpyc3RyaWRlc1swXSArIGIqcnN0cmlkZXNbMV0gKyBjKnJzdHJpZGVzWzJdIF0gPDwgc3RkOjplbmRsOwoKICByZXR1cm4gMDsKfQ==