Index: libcfa/src/containers/array.hfa
===================================================================
--- libcfa/src/containers/array.hfa	(revision 3eb55f9806679f91fcafa686a33e6c8fdb6a1639)
+++ libcfa/src/containers/array.hfa	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
@@ -25,4 +25,12 @@
 
     Timmed & ?[?]( arpk(Zn, S, Timmed, Tbase) & a, ptrdiff_t i ) {
+        return (Timmed &) a.strides[i];
+    }
+
+    Timmed & ?[?]( arpk(Zn, S, Timmed, Tbase) & a, int i ) {
+        return (Timmed &) a.strides[i];
+    }
+
+    Timmed & ?[?]( arpk(Zn, S, Timmed, Tbase) & a, size_t i ) {
         return (Timmed &) a.strides[i];
     }
@@ -80,44 +88,33 @@
 // Core -[[-,-,-]] operator
 
+#ifdef TRY_BROKEN_DESIRED_MD_SUBSCRIPT
+
 // Desired form.  One definition with recursion on IxBC (worked until Jan 2021, see trac #__TODO__)
-// forall( TA &, TB &, TC &, IxAB, IxBC ... | { TB & ?[?]( TA &, IxAB ); TC & ?[?]( TB &, IxBC ); } )
-// TC & ?[?]( TA & this, IxAB ab, IxBC bc ) {
-//     return this[ab][bc];
-// }
 
-// Workaround form.  Listing all possibilities up to 4 dims.
-forall( TA &, TB &, IxAB | { TB & ?[?]( TA &, IxAB ); }
-            , TC &, IxBC | { TC & ?[?]( TB &, IxBC ); } )
+forall( TA &, TB &, TC &, IxAB, IxBC ... | { TB & ?[?]( TA &, IxAB ); TC & ?[?]( TB &, IxBC ); } )
 TC & ?[?]( TA & this, IxAB ab, IxBC bc ) {
     return this[ab][bc];
 }
-forall( TA &, TB &, IxAB | { TB & ?[?]( TA &, IxAB ); }
-            , TC &, IxBC | { TC & ?[?]( TB &, IxBC ); }
-            , TD &, IxCD | { TD & ?[?]( TC &, IxCD ); } )
-TD & ?[?]( TA & this, IxAB ab, IxBC bc, IxCD cd ) {
-    return this[ab][bc][cd];
-}
-forall( TA &, TB &, IxAB | { TB & ?[?]( TA &, IxAB ); }
-            , TC &, IxBC | { TC & ?[?]( TB &, IxBC ); }
-            , TD &, IxCD | { TD & ?[?]( TC &, IxCD ); }
-            , TE &, IxDE | { TE & ?[?]( TD &, IxDE ); } )
-TE & ?[?]( TA & this, IxAB ab, IxBC bc, IxCD cd, IxDE de ) {
-    return this[ab][bc][cd][de];
+
+#else 
+
+// Workaround form.  Listing all possibilities up to 4 dims.
+
+forall( TA &, TB &, TC &, IxAB_0, IxBC | { TB & ?[?]( TA &, IxAB_0 ); TC & ?[?]( TB &, IxBC ); } )
+TC & ?[?]( TA & this, IxAB_0 ab, IxBC bc ) {
+    return this[ab][bc];
 }
 
-// Adapters for "indexed by ptrdiff_t" implies "indexed by [this other integral type]"
-// Work around restriction that assertions underlying -[[-,-,-]] must match excatly
-forall( C &, E & | { E & ?[?]( C &, ptrdiff_t ); } ) {
+forall( TA &, TB &, TC &, IxAB_0, IxAB_1, IxBC | { TB & ?[?]( TA &, IxAB_0, IxAB_1 ); TC & ?[?]( TB &, IxBC ); } )
+TC & ?[?]( TA & this, IxAB_0 ab0, IxAB_1 ab1, IxBC bc ) {
+    return this[[ab0,ab1]][bc];
+}
 
-    // Targeted to support:  for( i; z(N) ) ... a[[ ..., i, ... ]]
-    E & ?[?]( C & this, size_t i ) {
-        return this[ (ptrdiff_t) i ];
-    }
+forall( TA &, TB &, TC &, IxAB_0, IxAB_1, IxAB_2, IxBC | { TB & ?[?]( TA &, IxAB_0, IxAB_1, IxAB_2 ); TC & ?[?]( TB &, IxBC ); } )
+TC & ?[?]( TA & this, IxAB_0 ab0, IxAB_1 ab1, IxAB_2 ab2, IxBC bc ) {
+    return this[[ab0,ab1,ab2]][bc];
+}
 
-    // Targeted to support:  for( i; 5 ) ... a[[ ..., i, ... ]]
-    E & ?[?]( C & this, int i ) {
-        return this[ (ptrdiff_t) i ];
-    }
-}
+#endif
 
 //
Index: tests/array-container/.expect/array-basic.x64.txt
===================================================================
--- tests/array-container/.expect/array-basic.x64.txt	(revision 3eb55f9806679f91fcafa686a33e6c8fdb6a1639)
+++ tests/array-container/.expect/array-basic.x64.txt	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
@@ -1,6 +1,10 @@
 expect Ws             = 7.060606
 result Ws [][][][] lo = 7.060606
+result Ws [,,,]    lo = 7.060606
 result Ws [][][][] hi = 7.060606
+result Ws [,,,]    hi = 7.060606
 expect Xs             = 8.150808
 result Xs [][][][] lo = 8.150808
+result Xs [,,,]    lo = 8.150808
 result Xs [][][][] hi = 8.150808
+result Xs [,,,]    hi = 8.150808
Index: tests/array-container/.expect/array-md-sbscr-cases.x64.txt
===================================================================
--- tests/array-container/.expect/array-md-sbscr-cases.x64.txt	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
+++ tests/array-container/.expect/array-md-sbscr-cases.x64.txt	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
@@ -0,0 +1,1 @@
+done
Index: tests/array-container/array-basic.cfa
===================================================================
--- tests/array-container/array-basic.cfa	(revision 3eb55f9806679f91fcafa686a33e6c8fdb6a1639)
+++ tests/array-container/array-basic.cfa	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
@@ -105,14 +105,12 @@
     printf("result Ws [][][][] lo = %f\n", result);
 
-    // fixme: -[[-,-,-,-]] not working
-    // result = total1d_low( wxyz[[all, slice_ix, slice_ix, slice_ix]] );
-    // printf("result Ws [,,,]    lo = %f\n", result);
+    result = total1d_low( wxyz[[all, slice_ix, slice_ix, slice_ix]] );
+    printf("result Ws [,,,]    lo = %f\n", result);
 
     result = total1d_hi( wxyz[all][slice_ix][slice_ix][slice_ix] );
     printf("result Ws [][][][] hi = %f\n", result);
 
-    // fixme: -[[-,-,-,-]] not working
-    // result = total1d_hi( wxyz[[all, slice_ix, slice_ix, slice_ix]] );
-    // printf("result Ws [,,,]    hi = %f\n", result);
+    result = total1d_hi( wxyz[[all, slice_ix, slice_ix, slice_ix]] );
+    printf("result Ws [,,,]    hi = %f\n", result);
 
     // summing across X, with w=y=z=1
@@ -126,14 +124,12 @@
     printf("result Xs [][][][] lo = %f\n", result);
 
-    // fixme: -[[-,-,-,-]] not working
-    // result = total1d_low( wxyz[[slice_ix, all, slice_ix, slice_ix]] );
-    // printf("result Xs [,,,]    lo = %f\n", result);
+    result = total1d_low( wxyz[[slice_ix, all, slice_ix, slice_ix]] );
+    printf("result Xs [,,,]    lo = %f\n", result);
 
     result = total1d_hi( wxyz[slice_ix][all][slice_ix][slice_ix] );    
     printf("result Xs [][][][] hi = %f\n", result);
 
-    // fixme: -[[-,-,-,-]] not working
-    // result = total1d_hi( wxyz[[slice_ix, all, slice_ix, slice_ix]] );
-    // printf("result Xs [,,,]    hi = %f\n", result);
+    result = total1d_hi( wxyz[[slice_ix, all, slice_ix, slice_ix]] );
+    printf("result Xs [,,,]    hi = %f\n", result);
 
 }
Index: tests/array-container/array-md-sbscr-cases.cfa
===================================================================
--- tests/array-container/array-md-sbscr-cases.cfa	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
+++ tests/array-container/array-md-sbscr-cases.cfa	(revision 63a4b92a1b215265addfc287907ae92a2316f258)
@@ -0,0 +1,285 @@
+#include <containers/array.hfa>
+
+#include <assert.h>
+
+float getMagicNumber( ptrdiff_t w, ptrdiff_t x, ptrdiff_t y, ptrdiff_t z ) {
+
+    assert( 0 <= w && w < 3 );
+    assert( 0 <= x && x < 4 );
+    assert( 0 <= y && y < 5 );
+    assert( 0 <= z && z < 6 );
+
+    float ww = (2.0f \ w) / 1.0f;
+    float xx = (2.0f \ x) / 100.0f;
+    float yy = (2.0f \ y) / 10000.0f;
+    float Nz = (2.0f \ z) / 1000000.0f;
+
+    return ww+xx+yy+Nz;
+}
+
+forall( ztype(Nw), ztype(Nx), ztype(Ny), ztype(Nz) )
+void fillHelloData( array( float, Nw, Nx, Ny, Nz ) & wxyz ) {
+    for (w; z(Nw))
+    for (x; z(Nx))
+    for (y; z(Ny))
+    for (z; z(Nz))
+        wxyz[w][x][y][z] = getMagicNumber(w, x, y, z);
+}
+
+forall( ztype(Zn)
+      , S & | sized(S)
+      )
+float total1d_low( arpk(Zn, S, float, float ) & a ) {
+    float total = 0.0f;
+    for (i; z(Zn))
+        total += a[i];
+    return total;
+}
+
+forall( A & | ar(A, float) )
+float total1d_hi( A & a ) {
+    float total = 0.0f;
+    for (i; a`len)
+        total += a[i];
+    return total;
+}
+
+// Tests all the ways to split dimensions into CFA-supported chunks, by the only order that C supports: coarsest to finest stride.
+forall( ztype(Nw), ztype(Nx), ztype(Ny), ztype(Nz) )
+void test_inOrderSplits( tag(Nw), tag(Nx), tag(Ny), tag(Nz) ) {
+
+    array( float, Nw, Nx, Ny, Nz ) wxyz;
+    fillHelloData(wxyz);
+
+    ptrdiff_t iw = 2, ix = 3, iy=4, iz=5;
+
+    const float valExpected = getMagicNumber(iw, ix, iy, iz);
+    assert( wxyz[iw][ix][iy][iz] == valExpected );
+
+    // order wxyz, natural split (4-0 or 0-4, no intermediate to declare)
+
+    assert(( wxyz[[iw, ix, iy, iz]] == valExpected ));
+
+    // order wxyz, unnatural split 1-3  (three ways declared)
+
+    typeof( wxyz[iw] ) xyz1 = wxyz[iw];
+    assert(( xyz1[[ix, iy, iz]]  == valExpected ));
+
+    typeof( wxyz[iw] ) xyz2;
+    &xyz2 = &wxyz[iw];
+    assert(( xyz2[[ix, iy, iz]] == valExpected ));
+
+    assert(( wxyz[iw][[ix, iy, iz]] == valExpected ));
+
+    // order wxyz, unnatural split 2-2  (three ways declared)
+
+    typeof( wxyz[[iw, ix]] ) yz1 = wxyz[[iw,ix]];
+    assert(( yz1[[iy, iz]]  == valExpected ));
+
+    typeof( wxyz[[iw, ix]] ) yz2;
+    &yz2 = &wxyz[[iw, ix]];
+    assert(( yz2[[iy, iz]]  == valExpected ));
+
+    assert(( wxyz[[iw, ix]][[iy, iz]] == valExpected ));
+
+    // order wxyz, unnatural split 3-1  (three ways declared)
+
+    typeof( wxyz[[iw, ix, iy]] ) z1 = wxyz[[iw, ix, iy]];
+    assert(( z1[iz]  == valExpected ));
+
+    typeof( wxyz[[iw, ix, iy]] ) z2;
+    &z2 = &wxyz[[iw, ix, iy]];
+    assert(( z2[iz] == valExpected ));
+
+    assert(( wxyz[[iw, ix, iy]][iz] == valExpected ));
+}
+
+// All orders that skip a single dimension, each in its most natural split.
+forall( ztype(Nw), ztype(Nx), ztype(Ny), ztype(Nz) )
+void test_skipSingle( tag(Nw), tag(Nx), tag(Ny), tag(Nz) ) {
+
+    array( float, Nw, Nx, Ny, Nz ) wxyz;
+    fillHelloData(wxyz);
+
+    ptrdiff_t iw = 2, ix = 3, iy=4, iz=5;
+
+    const float valExpected = getMagicNumber(iw, ix, iy, iz);
+    assert( wxyz[iw][ix][iy][iz] == valExpected );
+
+
+    // order wxyz (no intermediates to declare)
+
+    assert(( wxyz[[iw  , ix  , iy  , iz  ]]       == valExpected ));
+    assert(( wxyz[[iw-1, ix  , iy  , iz  ]]       != valExpected ));
+
+    // order xyzw: *xyz, w
+
+    assert(( wxyz[[all , ix  , iy  , iz  ]][iw  ] == valExpected ));
+    assert(( wxyz[[all , ix-1, iy  , iz  ]][iw  ] != valExpected ));
+    assert(( wxyz[[all , ix  , iy  , iz  ]][iw-1] != valExpected ));
+
+    // order wyzx: w*yz, x
+
+    assert(( wxyz[[iw  , all , iy  , iz  ]][ix  ] == valExpected ));
+    assert(( wxyz[[iw  , all , iy-1, iz  ]][ix  ] != valExpected ));
+    assert(( wxyz[[iw  , all , iy  , iz  ]][ix-1] != valExpected ));
+
+    // order wxzy: wx*z, y
+
+    assert(( wxyz[[iw  , ix  , all , iz  ]][iy  ] == valExpected ));
+    assert(( wxyz[[iw  , ix  , all , iz-1]][iy  ] != valExpected ));
+    assert(( wxyz[[iw  , ix  , all , iz  ]][iy-1] != valExpected ));
+}
+
+
+// The comments specify a covering set of orders, each in its most natural split.
+// Covering means that each edge on the lattice of dimesnions-provided is used.
+// Natural split means the arity of every -[[-,...]] tuple equals the dimensionality of its "this" operand, then that the fewest "all" subscripts are given.
+// The commented-out test code shows cases that don't work.  We wish all the comment-coverd cases worked.
+forall( ztype(Nw), ztype(Nx), ztype(Ny), ztype(Nz) )
+void test_latticeCoverage( tag(Nw), tag(Nx), tag(Ny), tag(Nz) ) {
+
+    array( float, Nw, Nx, Ny, Nz ) wxyz;
+    fillHelloData(wxyz);
+
+    ptrdiff_t iw = 2, ix = 3, iy=4, iz=5;
+
+    const float valExpected = getMagicNumber(iw, ix, iy, iz);
+    assert( wxyz[iw][ix][iy][iz] == valExpected );
+
+
+    // order wxyz (no intermediates to declare)
+
+    assert(( wxyz[[iw, ix, iy, iz]] == valExpected ));
+
+    {
+        // order wyxz: w*y*, xz
+        assert( wxyz[iw][all][iy][all] [ix][iz] == valExpected );
+
+        typeof( wxyz[[iw, all, iy, all]] ) xz1 = wxyz[[iw, all, iy, all]];
+        assert(( xz1[[ix, iz]]  == valExpected ));
+
+        typeof( wxyz[[iw, all, iy, all]] ) xz2;
+        &xz2 = &wxyz[[iw, all, iy, all]];
+        assert(( xz2[[ix, iz]]  == valExpected ));
+
+        assert(( wxyz[[iw  , all, iy  , all]][[ix  , iz  ]] == valExpected ));
+        assert(( wxyz[[iw-1, all, iy  , all]][[ix  , iz  ]] != valExpected ));
+        assert(( wxyz[[iw  , all, iy-1, all]][[ix  , iz  ]] != valExpected ));
+        assert(( wxyz[[iw  , all, iy  , all]][[ix-1, iz  ]] != valExpected ));
+        assert(( wxyz[[iw  , all, iy  , all]][[ix  , iz-1]] != valExpected ));
+    }
+    {
+        // order wzxy: w**z, xy
+        assert( wxyz[iw][all][all][iz] [ix][iy] == valExpected );
+
+        // typeof( wxyz[[iw, all, all, iz]] ) xy1 = wxyz[[iw, all, all, iz]];
+        // assert(( xy1[[ix, iy]]  == valExpected ));
+
+        // typeof(  wxyz[[iw, all, all, iz]] ) xy2;
+        // &xy2 = &wxyz[[iw, all, all, iz]];
+        // assert(( xy2[[ix, iy]]  == valExpected ));
+
+        // assert(( wxyz[[iw  , all, all, iz  ]][[ix  , iy  ]] == valExpected ));
+        // assert(( wxyz[[iw-1, all, all, iz  ]][[ix  , iy  ]] != valExpected ));
+        // assert(( wxyz[[iw  , all, all, iz-1]][[ix  , iy  ]] != valExpected ));
+        // assert(( wxyz[[iw  , all, all, iz  ]][[ix-1, iy  ]] != valExpected ));
+        // assert(( wxyz[[iw  , all, all, iz  ]][[ix  , iy-1]] != valExpected ));
+    }
+    {
+        // order xywz: *xy*, wz
+        assert( wxyz[all][ix][iy][all] [iw][iz] == valExpected );
+
+        typeof( wxyz[[all, ix, iy, all]] ) wz1 = wxyz[[all, ix, iy, all]];
+        assert(( wz1[[iw, iz]]  == valExpected ));
+
+        assert(( wxyz[[all  , ix, iy  , all]][[iw  , iz  ]] == valExpected ));
+    }
+    {
+        // order xzwy: *x*z, wy
+        assert( wxyz[all][ix][all][iz] [iw][iy] == valExpected );
+
+        // assert(( wxyz[[all , ix  , all , iz  ]][[iw  , iy  ]] == valExpected ));
+    }
+    {
+        // order yzwx: **yz, wx
+        assert( wxyz[all][all][iy][iz] [iw][ix] == valExpected );
+
+        // assert(( wxyz[[all , all , iy  , iz  ]][[iw  , ix  ]] == valExpected ));
+    }
+    {
+        // order xwzy: *x**, w*z, y
+        assert( wxyz[all][ix][all][all] [iw][all][iz] [iy] == valExpected );
+
+        typeof( wxyz[all][ix][all][all] ) wyz_workaround = wxyz[[all , ix , all  , all  ]];
+        typeof( wyz_workaround[iw][all][iz] ) y_workaround = wyz_workaround[[iw , all , iz  ]];
+        assert( y_workaround[iy] == valExpected );
+
+        // assert(( wxyz[[all , ix , all  , all  ]][[iw  , all , iz  ]][iy  ] == valExpected ));
+    }
+    {
+        // order ywzx: **y*, w*z, x
+    }
+    {
+        // order zwyx: ***z, w*y, x
+    }
+    {
+        // order yxzw: **y*, *xz, w
+    }
+    {
+        // order zxyw: ***z, *xy, w
+    }
+    {
+        // order zyxw: ***z, **y, *x, w
+    }
+}
+
+forall( ztype(Nw), ztype(Nx), ztype(Ny), ztype(Nz) )
+void test_numSubscrTypeCompatibility( tag(Nw), tag(Nx), tag(Ny), tag(Nz) ) {
+
+    array( float, Nw, Nx, Ny, Nz ) wxyz;
+    fillHelloData(wxyz);
+
+    const float valExpected = getMagicNumber(2, 3, 4, 5);
+    assert(( wxyz [2] [3] [4] [5]  == valExpected ));
+    assert(( wxyz[[2,  3]][4] [5]  == valExpected ));
+    assert(( wxyz [2][[3,  4]][5]  == valExpected ));
+    assert(( wxyz [2] [3][[4,  5]] == valExpected ));
+    assert(( wxyz[[2,  3,  4]][5]  == valExpected ));
+    assert(( wxyz [2][[3,  4,  5]] == valExpected ));
+    assert(( wxyz[[2,  3,  4,  5]] == valExpected ));
+
+    for ( i; z(Nw) ) {
+        assert(( wxyz[[ i, 3, 4, 5 ]] == getMagicNumber(i, 3, 4, 5) ));
+    }
+
+    for ( i; z(Nx) ) {
+        assert(( wxyz[[ 2, i, 4, 5 ]] == getMagicNumber(2, i, 4, 5) ));
+    }
+
+    for ( i; z(Ny) ) {
+        assert(( wxyz[[ 2, 3, i, 5 ]] == getMagicNumber(2, 3, i, 5) ));
+    }
+
+    for ( i; z(Nz) ) {
+        assert(( wxyz[[ 2, 3, 4, i ]] == getMagicNumber(2, 3, 4, i) ));
+    }
+
+    for ( i; z(Nw) ) {
+        assert(( wxyz[[ i, all, 4, 5 ]][3] == getMagicNumber(i, 3, 4, 5) ));
+    }
+
+    for ( i; z(Nw) ) {
+        assert(( wxyz[[ all, 3, 4, 5 ]][i] == getMagicNumber(i, 3, 4, 5) ));
+    }
+}
+
+const size_t  KW = 3,  KX = 4,  KY = 5,  KZ = 6;
+
+int main() {
+
+    test_inOrderSplits  ( ztag(KW), ztag(KX), ztag(KY), ztag(KZ) );
+    test_skipSingle     ( ztag(KW), ztag(KX), ztag(KY), ztag(KZ) );
+    test_latticeCoverage( ztag(KW), ztag(KX), ztag(KY), ztag(KZ) );
+    printf("done\n");
+}