Index: libcfa/src/Makefile.am
===================================================================
--- libcfa/src/Makefile.am	(revision 7ee3c873b9e89a6ba296bf3372897c315fee1b41)
+++ libcfa/src/Makefile.am	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
@@ -56,4 +56,5 @@
 	bits/queue.hfa \
 	bits/sequence.hfa \
+	containers/array.hfa \
 	concurrency/iofwd.hfa \
 	containers/list.hfa \
Index: libcfa/src/containers/array.hfa
===================================================================
--- libcfa/src/containers/array.hfa	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
+++ libcfa/src/containers/array.hfa	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
@@ -0,0 +1,149 @@
+
+
+// a type whose size is n
+#define Z(n) char[n]
+
+// the inverse of Z(-)
+#define z(Zn) sizeof(Zn)
+
+// if you're expecting a Z(n), say so, by asking for a ztype, instead of dtype or otype
+#define ztype(Zn) Zn & | sized(Zn)
+
+forall( T & ) struct tag {};
+#define ttag(T) ((tag(T)){})
+#define ztag(n) ttag(Z(n))
+
+
+//
+// Single-dim array sruct (with explicit packing and atom)
+//
+
+forall( ztype(Zn), ztype(S), Timmed &, Tbase & ) {
+    struct arpk {
+        S strides[z(Zn)];
+    };
+
+    Timmed & ?[?]( arpk(Zn, S, Timmed, Tbase) & a, ptrdiff_t i ) {
+        return (Timmed &) a.strides[i];
+    }
+
+    size_t ?`len( arpk(Zn, S, Timmed, Tbase) & a ) {
+        return z(Zn);
+    }
+
+    // workaround #226 (and array relevance thereof demonstrated in mike102/otype-slow-ndims.cfa)
+    void ?{}( arpk(Zn, S, Timmed, Tbase) & this ) {
+        void ?{}( S (&inner)[z(Zn)] ) {}
+        ?{}(this.strides);
+    }
+    void ^?{}( arpk(Zn, S, Timmed, Tbase) & this ) {
+        void ^?{}( S (&inner)[z(Zn)] ) {}
+        ^?{}(this.strides);
+    }
+}
+
+//
+// Sugar for declaring array structure instances
+//
+
+forall( Te )
+Te mkar_( tag(Te) ) {}
+
+forall( ztype(Zn), ZTags ... , Trslt &, Tatom & | { Trslt mkar_( tag(Tatom), ZTags ); } )
+arpk(Zn, Trslt, Trslt, Tatom) mkar_( tag(Tatom), tag(Zn), ZTags ) {}
+
+// based on https://stackoverflow.com/questions/1872220/is-it-possible-to-iterate-over-arguments-in-variadic-macros
+
+    // Make a FOREACH macro
+    #define FE_0(WHAT)
+    #define FE_1(WHAT, X) WHAT(X) 
+    #define FE_2(WHAT, X, ...) WHAT(X)FE_1(WHAT, __VA_ARGS__)
+    #define FE_3(WHAT, X, ...) WHAT(X)FE_2(WHAT, __VA_ARGS__)
+    #define FE_4(WHAT, X, ...) WHAT(X)FE_3(WHAT, __VA_ARGS__)
+    #define FE_5(WHAT, X, ...) WHAT(X)FE_4(WHAT, __VA_ARGS__)
+    //... repeat as needed
+
+    #define GET_MACRO(_0,_1,_2,_3,_4,_5,NAME,...) NAME 
+    #define FOR_EACH(action,...) \
+    GET_MACRO(_0,__VA_ARGS__,FE_5,FE_4,FE_3,FE_2,FE_1,FE_0)(action,__VA_ARGS__)
+
+#define COMMA_ttag(X) , ttag(X)
+#define array( TE, ...) typeof( mkar_( ttag(TE)  FOR_EACH( COMMA_ttag, __VA_ARGS__ ) ) )
+
+#define COMMA_ztag(X) , ztag(X)
+#define zarray( TE, ...) typeof( mkar_( ttag(TE)  FOR_EACH( COMMA_ztag, __VA_ARGS__ ) ) )
+
+//
+// Sugar for multidimensional indexing
+//
+
+// Core -[[-,-,-]] operator
+
+// 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 ); } )
+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];
+}
+
+// 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 ); } ) {
+
+    // Targeted to support:  for( i; z(N) ) ... a[[ ..., i, ... ]]
+    E & ?[?]( C & this, size_t i ) {
+        return this[ (ptrdiff_t) i ];
+    }
+
+    // Targeted to support:  for( i; 5 ) ... a[[ ..., i, ... ]]
+    E & ?[?]( C & this, int i ) {
+        return this[ (ptrdiff_t) i ];
+    }
+}
+
+//
+// Rotation
+//
+
+// Base
+forall( ztype(Zq), ztype(Sq), Tbase & )
+tag(arpk(Zq, Sq, Tbase, Tbase)) enq_( tag(Tbase), tag(Zq), tag(Sq), tag(Tbase) ) {}
+
+// Rec
+forall( ztype(Zq), ztype(Sq), ztype(Z), ztype(S), recq &, recr &, Tbase & | { tag(recr) enq_( tag(Tbase), tag(Zq), tag(Sq), tag(recq) ); } )
+tag(arpk(Z, S, recr, Tbase)) enq_( tag(Tbase), tag(Zq), tag(Sq), tag(arpk(Z, S, recq, Tbase)) ) {}
+
+// Wrapper
+struct all_t {} all;
+forall( ztype(Z), ztype(S), Te &, result &, Tbase & | { tag(result) enq_( tag(Tbase), tag(Z), tag(S), tag(Te) ); } )
+result & ?[?]( arpk(Z, S, Te, Tbase) & this, all_t ) {
+    return (result&) this;
+}
+
+//
+// Trait of array or slice
+//
+
+trait ar(A &, Tv &) {
+    Tv& ?[?]( A&, ptrdiff_t );
+    size_t ?`len( A& );
+};
Index: tests/array-container/.expect/array-basic.txt
===================================================================
--- tests/array-container/.expect/array-basic.txt	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
+++ tests/array-container/.expect/array-basic.txt	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
@@ -0,0 +1,6 @@
+expect Ws             = 7.060606
+result Ws [][][][] lo = 7.060606
+result Ws [][][][] hi = 7.060606
+expect Xs             = 8.150808
+result Xs [][][][] lo = 8.150808
+result Xs [][][][] hi = 8.150808
Index: tests/array-container/array-basic.cfa
===================================================================
--- tests/array-container/array-basic.cfa	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
+++ tests/array-container/array-basic.cfa	(revision c7625e0bd9985f4ba9e0f9199ca99d157c60506d)
@@ -0,0 +1,147 @@
+#include <containers/array.hfa>
+
+//
+// Type-theory demo (success is "it compiles")
+//
+
+forall( ztype(Nx), ztype(Ny), ztype(Nz) )
+void typesTest( tag(Nx), tag(Ny), tag(Nz) ) {
+
+    array( float, Nx, Ny, Nz )  xyz;
+
+    // numeric subscripts
+    ptrdiff_t ix = 1, iy = 2, iz = 3;
+    array( float, Ny, Nz ) & yz = xyz[ix];
+    array( float, Nz ) & z = xyz[ix][iy];
+    float & val = xyz[ix][iy][iz];
+
+    // deferral subscript
+    typeof( xyz[all] ) yzx = xyz[all];
+
+    // longform check that -[all] gives intended type (user doesn't write this)
+    arpk( Ny
+        , array( float, Nz)
+        , arpk( Nz
+              , float
+              , arpk( Nx
+                    , array( float, Ny, Nz )
+                    , float
+                    , float
+                    )
+              , float
+              )
+        , float
+        )
+        & yzx_long = xyz[all];
+    & yzx_long = & yzx;
+    & yzx  = & yzx_long;
+}
+
+//
+// Runtime demo
+//
+
+#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;
+}
+
+forall( ztype(Nw), ztype(Nx), ztype(Ny), ztype(Nz) )
+void runtimeTest( tag(Nw), tag(Nx), tag(Ny), tag(Nz) ) {
+
+    array( float, Nw, Nx, Ny, Nz ) wxyz;
+    fillHelloData(wxyz);
+
+    float expect, result;
+    ptrdiff_t slice_ix = 1;
+
+    // summing across W, with x=y=z=1
+
+    expect = 0;
+    for (i; z(Nw))
+        expect += getMagicNumber( i, slice_ix, slice_ix, slice_ix );
+    printf("expect Ws             = %f\n", expect);
+
+    result = total1d_low( wxyz[all][slice_ix][slice_ix][slice_ix] );
+    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_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);
+
+    // summing across X, with w=y=z=1
+
+    expect = 0;
+    for (i; z(Nx))
+        expect += getMagicNumber( slice_ix, i, slice_ix, slice_ix );
+    printf("expect Xs             = %f\n", expect);
+
+    result = total1d_low( wxyz[slice_ix][all][slice_ix][slice_ix] );    
+    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_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);
+
+}
+
+const size_t  KW = 3,  KX = 4,  KY = 5,  KZ = 6;
+
+int main() {
+
+    typesTest  (           ztag(KX), ztag(KY), ztag(KZ) );
+    runtimeTest( ztag(KW), ztag(KX), ztag(KY), ztag(KZ) );
+}