// R = return type
// A, B, C = matrix types (may be NO_MATRIX)
// a, b, c are names of corresponding matrices
// s is not indexed (usually a scalar)
// dyn_m, dyn_n = dynamic size (only used for R=MATRIX/MATRIX_MXN/MATRIX_VIEW)
// #define FUNCTION(aij,bij,cij,s)
// Optionally define TRANS_A, TRANS_B, TRANS_C 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

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

{
    constexpr int ac=access_flags(A::type)&access_flags(B::type)&access_flags(C::type);
    if constexpr(R::type==mat_zero) return R{};
    else if constexpr(R::type==mat_id) return R{};
    else if constexpr(R::type==mat_scale)
    {
        // Careful not to access out of bounds.
        if constexpr(R::m==0) return R{};
        else return R{FUNCTION(diag(a,0),diag(b,0),diag(c,0),s)};
    }
    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,auto&& c,auto&& s,std::integer_sequence<int,i...>)
                {
                    return R{(FUNCTION(
                                (get<TRANS_A(i,i)>(a)),
                                (get<TRANS_B(i,i)>(b)),
                                (get<TRANS_C(i,i)>(c)),s))...};
                }(a,b,c,s,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)=FUNCTION(diag(a,i),diag(b,i),diag(c,i),s);
            return r;
        }
    }
    else if constexpr(R::type==mat_sym)
    {
        // If we are small enough to hard-code.
        if constexpr(R::m<=size_limit)
        {
            return []<int... i>(auto&& a,auto&& b,auto&& c,auto&& s,std::integer_sequence<int,i...>)
                {
                    return R{(FUNCTION(
                                (get<TRANS_A(sym_i(R::m,i),sym_j(R::m,i))>(a)),
                                (get<TRANS_B(sym_i(R::m,i),sym_j(R::m,i))>(b)),
                                (get<TRANS_C(sym_i(R::m,i),sym_j(R::m,i))>(c)),s))...};
                }(a,b,c,s,std::make_integer_sequence<int,R::m*(R::m+1)/2>{});
        }
        else if constexpr(ac&1) // we can use lower
        {
            R r;
            for(int j=0;j<R::m;j++)
                for(int i=j;i<R::m;i++)
                    r.Element_Lower(i,j)=FUNCTION(lower(a,TRANS_A(i,j)),lower(b,TRANS_B(i,j)),lower(c,TRANS_C(i,j)),s);
            return r;
        }
        else // must use diag/lower_
        {
            R r;
            for(int j=0;j<R::m;j++)
                for(int i=j+1;i<R::m;i++)
                    r.Element_Lower(i,j)=FUNCTION(lower_(a,TRANS_A(i,j)),lower_(b,TRANS_B(i,j)),lower_(c,TRANS_C(i,j)),s);
            for(int i=0;i<R::m;i++)
                r.Element_Lower(i,i)=FUNCTION(diag(a,i),diag(b,i),diag(c,i),s);
            return r;
        }
    }
    else // mat_gen
    {
        // If we are small enough to hard-code.
        if constexpr(R::fixed_m>=0 && R::fixed_m<=size_limit && R::fixed_n<=size_limit)
        {
            return []<int... i>(auto&& a,auto&& b,auto&& c,auto&& s,std::integer_sequence<int,i...>)
                {
                    return R{(FUNCTION(
                                (get<TRANS_A(i%R::m,i/R::m)>(a)),
                                (get<TRANS_B(i%R::m,i/R::m)>(b)),
                                (get<TRANS_C(i%R::m,i/R::m)>(c)),s))...};
                }(a,b,c,s,std::make_integer_sequence<int,R::m*R::n>{});
        }
        else if constexpr(ac&2) // we can use get
        {
            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)=FUNCTION(get(a,TRANS_A(i,j)),get(b,TRANS_B(i,j)),get(c,TRANS_C(i,j)),s);
            return r;
        }
        else if constexpr(ac&1) // we can use lower/upper
        {
            assert(dyn_m==dyn_n);
            R r((INITIAL_SIZE)dyn_m,(INITIAL_SIZE)dyn_n);
            for(int j=0;j<dyn_n;j++)
            {
                for(int i=0;i<j;i++)
                    r(i,j)=FUNCTION(upper(a,TRANS_A(i,j)),upper(b,TRANS_B(i,j)),upper(c,TRANS_C(i,j)),s);
                for(int i=j;i<dyn_m;i++)
                    r(i,j)=FUNCTION(lower(a,TRANS_A(i,j)),lower(b,TRANS_B(i,j)),lower(c,TRANS_C(i,j)),s);
            }
            return r;
        }
        else // must use diag/lower_
        {
            assert(dyn_m==dyn_n);
            R r((INITIAL_SIZE)dyn_m,(INITIAL_SIZE)dyn_n);
            for(int j=0;j<dyn_n;j++)
            {
                for(int i=0;i<j;i++)
                    r(i,j)=FUNCTION(upper_(a,TRANS_A(i,j)),upper_(b,TRANS_B(i,j)),upper_(c,TRANS_C(i,j)),s);
                r(j,j)=FUNCTION(diag(a,j),diag(b,j),diag(c,j),s);
                for(int i=j+1;i<dyn_m;i++)
                    r(i,j)=FUNCTION(lower_(a,TRANS_A(i,j)),lower_(b,TRANS_B(i,j)),lower_(c,TRANS_C(i,j)),s);
            }
            return r;
        }
    }
}
#undef FUNCTION
#undef TRANS_A
#undef TRANS_B
#undef TRANS_C
