#include "array.hfa"


// TODO: can I mini-sweeten the type overloads of -[[-,-,-]] beyond ptrdiff_t
// so that it works with for(i; z(-))

// PARK: refining the feature idea of equationally-defined operators, to include
// declaration:  foo a = b * c;
// meaning construct a to be b times c

forall( ztype(R), ztype(C), ztype(Mid) )
void mul( array( float, R, Mid) & srcL,
          array( float, Mid, C) & srcR,
          array( float, R, C ) & dst ) {
    for (ptrdiff_t r = 0; r < z(R); r++) for (ptrdiff_t c = 0; c < z(C); c++) {
        dst[[r,c]] = 0;
        for (ptrdiff_t i = 0; i < z(Mid); i++)
            dst[[r,c]] += srcL[[r,i]] * srcR[[i,c]];
    }
}