// R = return type
// A, B = matrix types
// a, b are names of corresponding matrices
// m, n, p = fixed sizes
// dyn_m, dyn_n, dyn_p = dynamic sizes
// Optionally define TRANS_A, TRANS_B to access those matrices transposed.
#ifdef TRANS_A
#undef TRANS_A
#define TRANS_A(x,y) y,x
#else
#define TRANS_A(x,y) x,y
#endif

#ifdef TRANS_B
#undef TRANS_B
#define TRANS_B(x,y) y,x
#else
#define TRANS_B(x,y) x,y
#endif

{
    if constexpr(R::type==mat_zero) return R{};
    else if constexpr(R::type==mat_id) return R{};
    else if constexpr(R::type==mat_scale) // SCALE_MATRIX
    {
        // Careful not to access out of bounds.
        if constexpr(R::m==0) return R{};
        else return R{diag(a,0)*diag(b,0)};
    }
    else if constexpr(R::type==mat_diag)
    {
        // If we are small enough to hard-code.
        if constexpr(R::m<=size_limit)
        {
            return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
                {
                    return R{(get<TRANS_A(i,i)>(a)*get<TRANS_B(i,i)>(b))...};
                }(a,b,std::make_integer_sequence<int,R::m>{});
        }
        else // Too big, generate dynamically.
        {
            R r;
            for(int i=0;i<r.m;i++)
                r.x(i)=diag(a,i)*diag(b,i);
            return r;
        }
    }
    else if constexpr(R::type==mat_sym)
    {
        // If we are small enough to hard-code.
        if constexpr(m<=size_limit)
        {
            if constexpr(A::is_diag) // no loop over k, since k=i
            {
                return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
                    {
                        return R{(get<sym_i(m,i),sym_i(m,i)>(a)*
                                get<TRANS_B(sym_i(m,i),sym_j(m,i))>(b))...};
                    }(a,b,std::make_integer_sequence<int,m*(m+1)/2>{});
            }
            else if constexpr(B::is_diag) // no loop over k, since k=j
            {
                return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
                    {
                        return R{(get<TRANS_A(sym_i(m,i),sym_j(m,i))>(a)*
                                get<sym_j(m,i),sym_j(m,i)>(b))...};
                    }(a,b,std::make_integer_sequence<int,m*(m+1)/2>{});
            }
            else if constexpr(p==0) return R{};
            else
            {
                return []<int... i,int... k>(auto&& a,auto&& b,std::integer_sequence<int,i...>,std::integer_sequence<int,k...>)
                    {
                        int j{};(void)j;
                        return R{(j=i,((a(TRANS_A(sym_idx(m,j).first,k))*b(TRANS_B(k,sym_idx(m,j).second)))+...))...};
                    }(a,b,std::make_integer_sequence<int,m*(m+1)/2>{},std::make_integer_sequence<int,p>{});
            }
        }
        else if constexpr(A::is_diag) // no loop over k, since k=i
        {
            R r;
            for(int j=0;j<m;j++)
                for(int i=j;i<m;i++)
                    r.Element_Lower(i,j)=diag(a,i)*b(TRANS_B(i,j));
            return r;
        }
        else if constexpr(B::is_diag) // no loop over k, since k=j
        {
            R r;
            for(int j=0;j<m;j++)
                for(int i=j;i<m;i++)
                    r.Element_Lower(i,j)=a(TRANS_A(i,j))*diag(b,j);
            return r;
        }
        else
        {
            R r;
            for(int j=0;j<m;j++)
                for(int i=j;i<m;i++)
                {
                    T z{};
                    for(int k=0;k<dyn_p;k++)
                        z+=a(TRANS_A(i,k))*b(TRANS_B(k,j));
                    r.Element_Lower(i,j)=z;
                }
            return r;
        }
    }
    else // mat_gen
    {
        // If we are small enough to hard-code.
        if constexpr(m>=0 && n>=0 && p>=0 && m<=size_limit && n<=size_limit && p<=size_limit)
        {
            if constexpr(A::is_diag) // no loop over k, since k=i
            {
                return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
                    {
                        return R{(get<i%m,i%m>(a)*get<TRANS_B(i%m,i/m)>(b))...};
                    }(a,b,std::make_integer_sequence<int,m*n>{});
            }
            else if constexpr(B::is_diag) // no loop over k, since k=j
            {
                return []<int... i>(auto&& a,auto&& b,std::integer_sequence<int,i...>)
                    {
                        return R{(get<TRANS_A(i%m,i/m)>(a)*get<i/m,i/m>(b))...};
                    }(a,b,std::make_integer_sequence<int,m*n>{});
            }
            else if constexpr(p==0) return R{};
            else
            {
                return []<int... i,int... k>(auto&& a,auto&& b,std::integer_sequence<int,i...>,std::integer_sequence<int,k...>)
                    {
                        int j{};
                        return R{(j=i,((a(TRANS_A(j%m,k))*b(TRANS_B(k,j/m)))+...))...};
                    }(a,b,std::make_integer_sequence<int,m*n>{},std::make_integer_sequence<int,p>{});
            }
        }
        else if constexpr(A::is_diag) // no loop over k, since k=i
        {
            R r((INITIAL_SIZE)dyn_m,(INITIAL_SIZE)dyn_n);
            for(int j=0;j<dyn_n;j++)
                for(int i=0;i<dyn_m;i++)
                    r(i,j)=diag(a,i)*b(TRANS_B(i,j));
            return r;
        }
        else if constexpr(B::is_diag) // no loop over k, since k=j
        {
            R r((INITIAL_SIZE)dyn_m,(INITIAL_SIZE)dyn_n);
            for(int j=0;j<dyn_n;j++)
                for(int i=0;i<dyn_m;i++)
                    r(i,j)=a(TRANS_A(i,j))*diag(b,j);
            return r;
        }
        else
        {
            R r((INITIAL_SIZE)dyn_m,(INITIAL_SIZE)dyn_n);
            for(int j=0;j<dyn_n;j++)
                for(int i=0;i<dyn_m;i++)
                {
                    T z{};
                    for(int k=0;k<dyn_p;k++)
                        z+=a(TRANS_A(i,k))*b(TRANS_B(k,j));
                    r(i,j)=z;
                }
            return r;
        }
    }
}
#undef TRANS_A
#undef TRANS_B
