#define PROG_NAME "test_indexing"
#define PROG_DESC "test of {indexing.h}"
#define PROG_VERS "1.0"

/* Last edited on 2009-09-01 18:56:47 by stolfi */ 
/* Created on 2005-02-14 (or earlier) by J. Stolfi, UNICAMP */

#define test_indexing_COPYRIGHT \
  "Copyright © 2005  by the State University of Campinas (UNICAMP)"

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <assert.h>

#include <indexing.h>
#include <affirm.h>
#include <jsrandom.h>

#define bug(FMT_AND_ARGS...) \
  do { \
    fprintf(stderr, "%s:%d: ", __FILE__, __LINE__); \
    fprintf(stderr, FMT_AND_ARGS); \
    fprintf(stderr, "\n"); \
    exit(1); \
  } while(0)

#define MAXDIM 6

typedef struct desc_t
  { ix_dim_t d;
    ix_size_t z[MAXDIM];
    ix_pos_t b;
    ix_step_t s[MAXDIM];
    int nel;     /* Number of elements allocated for  {el}. */
    uint8_t *el; /* Element storage area. */
  } desc_t;

desc_t *make_desc(ix_dim_t d, ix_size_t sz[]);
  /* Creates a descriptor with axis count {d} and
    size vector {sz[0..d-1]}.  Some indices {i}
    may have virtual replication ({st[i] == 0}
    with {sz[i] > 1}). */

desc_t *make_packed_desc(ix_dim_t d, ix_size_t sz[], ix_order_t ixor);
  /* Creates a descriptor with axis count {d} and
    size vector {sz[0..d-1]}, with steps consistent
    with a packed element arrangement with index
    order {L}, and no virtual replication. */

void do_test_suite(int nt);
  /* Repeats the full test suite {nt} times. */

void extend_size_vector(ix_dim_t dA, ix_size_t szA[], ix_dim_t dB, ix_size_t szB[]);
  /* Copies {szA[0..dA-1]} to {szB[0..dA-1]} and fills { {szB[dA..dB-1]} with '1'. */

void test_size_ops(ix_dim_t d);
/* void test_new_array(desc_t *A, ix_size_t *sz); */
void test_addressing(desc_t *A);
void test_enum(desc_t *A);
void test_packed(desc_t *A);
/* void test_indices(desc_t *A); */
void test_crop(desc_t *A);
void test_subsample(desc_t *A);
void test_swap_indices(desc_t *A);
void test_flip_indices(desc_t *A);
void test_flip(desc_t *A);
void test_diagonal(desc_t *A);
void test_chop(desc_t *A);
void test_replicate(desc_t *A);
void test_reduce(void);

void print_desc(FILE *wr, char *pf, desc_t *A, char *sf);
void print_indices(FILE *wr, char *pf, ix_dim_t d, ix_index_t ix[], char *sf);

/* void dump_storage(FILE *wr, void *el, int nw, ix_nbits_t bpw); */
void check_sizes(ix_dim_t d, ix_size_t sza[], ix_size_t szb[]);
void check_max_indices(ix_dim_t d, ix_index_t ixa[], ix_index_t ixb[]);
void check_max_indices_sz(ix_dim_t d, ix_index_t ixa[], ix_size_t szb[]);

int main (int argn, char **argv)
{
  
  /* General limit and data type tests: */
  demand((ix_pos_NONE < 0) || (ix_pos_NONE > ix_MAX_POS), "invalid ix_pos_NONE");
  do_test_suite(30);
  return 0;
}
 
void do_test_suite(int nt)
{
  int it;
  ix_size_t sz[MAXDIM];
  for (it = 0; it < nt; it++)
    { 
      ix_dim_t d = (it % MAXDIM);
      test_size_ops(d);

      /* An array for testing: */
      sz[0] = abrandom(0, 2) + abrandom(0, 2); 
      sz[1] = abrandom(0, 5) + abrandom(0, 5); 
      sz[2] = 1; 
      sz[3] = abrandom(0, 9) + abrandom(0, 9); 
      sz[4] = abrandom(0, 7) + abrandom(0, 7);
      sz[5] = abrandom(0, 3) + abrandom(0, 3);
      assert(MAXDIM == 6);

      desc_t *A = make_desc(d, sz);

      /* test_new_array(A, sz, bps, bpw); */
      test_packed(A);
      test_enum(A);
      test_addressing(A);
      test_crop(A);
      test_subsample(A);
      test_swap_indices(A);
      test_flip_indices(A);
      test_flip(A);
      test_diagonal(A);
      test_chop(A);
      test_replicate(A);
      test_reduce();
    }
}

desc_t *make_desc(ix_dim_t d, ix_size_t sz[])
{
  fprintf(stderr, "Creating a descriptor (d = %d)...\n", d);
  desc_t *A = notnull(malloc(sizeof(desc_t)), "no mem");
  ix_axis_t iz = abrandom(0,2*d+1);     /* The index that will get zero step. */
  A->d = d;
  int i;
  int64_t npos = 1;
  for (i = 0; i < d; i++)
    { A->z[i] = sz[i];
      if ((i == iz) || (sz[i] <= 1))
        { /* Use virtual replication on this index: */
          A->s[i] = 0;
        }
      else
        { /* Use distinct elements on this index: */
          A->s[i] = (sz[i] == 1 ? 0 : npos);
        }
      /* Update {npos} if axis is neither trivial nor replicated: */
      if (A->z[i] == 0)
        { npos = 0; }
      else if ((A->z[i] > 1) && (A->s[i] != 0))
        { npos *= sz[i]; }
    }
  /* If the array is empty, all steps and the base must be zero: */
  if (npos == 0)
    { for (i = 0; i < d; i++) { A->s[i] = 0; } 
      A->b = 0;
    }
  else
    { /* A random base: */
      A->b = 1000*abrandom(1,9);
    }
  if (! ix_parms_are_valid(A->d, A->z, A->b, A->s, TRUE))
    { print_desc(stderr, "A = { ", A, " }\n");
      bug("make_desc creates an invalid desc");
      assert(FALSE);
    }
  return A;
}

desc_t *make_packed_desc(ix_dim_t d, ix_size_t sz[], ix_order_t ixor)
{
  desc_t *A = notnull(malloc(sizeof(desc_t)), "no mem");
  A->d = d;
  ix_sizes_assign(d, A->z, sz);
  ix_packed_steps(d, sz, ixor, A->s);
  A->b = 0;
  if (! ix_parms_are_valid(A->d, A->z, A->b, A->s, FALSE))
    { print_desc(stderr, "A = { ", A, " }\n");
      bug("make_desc creates an invalid desc");
      assert(FALSE);
    }
  return A;
}

void test_size_ops(ix_dim_t d)
{
  fprintf(stderr, "Checking ix_sizes_assign, ix_sizes_shrink, ix_sizes_expand (d = %d)...\n", d);
  bool_t ch_r, ch_g; 

  /* Size vectors for tests: */
  ix_size_t sza[MAXDIM];
  ix_size_t szb[MAXDIM];
  ix_size_t szc[MAXDIM];
  ix_size_t szd[MAXDIM];
  
  /* Initialize {sza,szb,szc} with random values: */
  auto ix_size_t rnd_sz(int a, int b, int c, int i); 
    /* A funct of {a,b,c}, and {i}. */

  ix_size_t rnd_sz(int a, int b, int c, int i) 
    { return a + ((b*i*i) % c); }
  
  int i;

  /* Test {ix_size_assign}: */
  /* fprintf(stderr, "  checking {ix_size_assign}\n"); */
  for (i = 0; i < MAXDIM; i++) { sza[i] = szc[i] = rnd_sz(17,  8, 13, i); }
  for (i = 0; i < MAXDIM; i++) { szb[i] = rnd_sz(19, 10, 17, i); }
  ix_sizes_assign(d, sza, szb);
  for (i = 0; i < MAXDIM; i++)
    { ix_size_t g = (i < d ? szb[i] : szc[i]); 
      if (sza[i] != g) 
        { bug("ix_size_assign mismatch i = %d sza = %llu should be %llu", i, sza[i], g); }
    }
    
  /* Test {ix_size_min}: */
  /* fprintf(stderr, "  checking {ix_size_min}\n"); */
  for (i = 0; i < MAXDIM; i++) { sza[i] = rnd_sz(17,  8, 13, i); }
  ix_sizes_assign(MAXDIM, szb, sza); /* Save current value of {sza} in {szb}. */
  ch_r = ix_sizes_shrink(d, sza, szc);
  ch_g = FALSE;
  for (i = 0; i < MAXDIM; i++)
    { ix_size_t g = (i < d ? (szb[i] < szc[i] ? szb[i] : szc[i]) : szb[i]); 
      if (sza[i] != g) 
        { bug("ix_size_min mismatch i = %d sza = %llu should be %llu", i, sza[i], g); }
      if ((i < d) && (sza[i] != szb[i])) { ch_g = TRUE; }
    }
  if (ch_r != ch_g)
    { bug("ix_size_min wrong return value"); }
    
  /* Test {ix_size_max}: */
  /* fprintf(stderr, "  checking {ix_size_max}\n"); */
  for (i = 0; i < MAXDIM; i++) { sza[i] = rnd_sz(17,  8, 13, i); }
  ix_sizes_assign(MAXDIM, szb, sza); /* Save current value of {sza} in {szb}. */
  ch_r = ix_sizes_expand(d, sza, szd);
  ch_g = FALSE;
  for (i = 0; i < MAXDIM; i++)
    { ix_size_t g = (i < d ? (szb[i] > szd[i] ? szb[i] : szd[i]) : szb[i]); 
      if (sza[i] != g) 
        { bug("ix_size_max mismatch i = %d sza = %llu should be %llu", i, sza[i], g); }
      if ((i < d) && (sza[i] != szb[i])) { ch_g = TRUE; }
    }
  if (ch_r != ch_g)
    { bug("ix_size_max wrong return value"); }
}

void test_addressing(desc_t *A)
{
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_position, ix_position_safe, setting and checking all elems (d = %d nt = %llu)...\n", A->d, nt);
  /* Clear elements: */
  int i;
  for (i = 0; i < A->nel; i++) { A->el[i] = 0; }
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  ix_index_t ixA[MAXDIM];

  /* Allocate a large enough boolean array: */
  int64_t npos = 1;
  { int i; for (i = 0; i < A->d; i++) { npos *= A->z[i]; } }
  int64_t nel = A->b + npos + 1000;
  uint8_t *el = (uint8_t *)notnull(malloc(nel*sizeof(uint8_t)), "no mem");
  
  /* Get the size vector, extend with 1's: */
  ix_size_t sz[MAXDIM];
  extend_size_vector(A->d, A->z, MAXDIM, sz);

  /* Reduce the size to 1 whenever the step is 0: */
  for (i = 0; i < A->d; i++) 
    { if ((A->z[i] > 1) && (A->s[i] == 0)) { sz[i] = 1; } }

  assert(MAXDIM == 6);
  for (ixA[5] = 0; ixA[5] < sz[5]; ixA[5]++)
    for (ixA[4] = 0; ixA[4] < sz[4]; ixA[4]++)
      for (ixA[3] = 0; ixA[3] < sz[3]; ixA[3]++)
        for (ixA[2] = 0; ixA[2] < sz[2]; ixA[2]++)
          for (ixA[1] = 0; ixA[1] < sz[1]; ixA[1]++)
            for (ixA[0] = 0; ixA[0] < sz[0]; ixA[0]++)
              { 
                /* Compute the expected address: */
                ix_pos_t pos_exp = A->b;
                for (i = 0; i < A->d; i++) { pos_exp += A->s[i]*ixA[i]; }
                
                /* Compare with library procs: */
                if (! ix_is_valid(A->d, ixA, A->z)) { bug("not ix_is_valid"); } 
                ix_pos_t pos = ix_position(A->d, ixA, A->b, A->s);
                if (pos != pos_exp) 
                  { print_desc(stderr, "A = { ", A, " }\n");
                    print_indices(stderr, "ixA = [ ", A->d, ixA, " ]\n");
                    bug("ix_position mismatch pos = %llu pos_exp = %llu", pos, pos_exp);
                  } 
                ix_pos_t pos_safe = ix_position_safe(A->d, ixA, A->z, A->b, A->s);
                if (pos != pos_safe) 
                  { print_desc(stderr, "A = { ", A, " }\n");
                    print_indices(stderr, "ixA = [ ", A->d, ixA, " ]\n");
                    bug("ix_position_safe mismatch pos = %llu pos_safe = %llu", pos, pos_safe);
                  } 
                uint8_t smp_exp = 1 + (pos % 253);
                el[pos] = smp_exp;
              }

  fprintf(stderr, "Re-checking the value of all samples (d = %d)...\n", A->d);
  assert(MAXDIM == 6);
  for (ixA[5] = 0; ixA[5] < sz[5]; ixA[5]++)
    for (ixA[4] = 0; ixA[4] < sz[4]; ixA[4]++)
      for (ixA[3] = 0; ixA[3] < sz[3]; ixA[3]++)
        for (ixA[2] = 0; ixA[2] < sz[2]; ixA[2]++)
          for (ixA[1] = 0; ixA[1] < sz[1]; ixA[1]++)
            for (ixA[0] = 0; ixA[0] < sz[0]; ixA[0]++)
              { 
                ix_pos_t pos = ix_position(A->d, ixA, A->b, A->s);
                uint8_t smp_exp = 1 + (pos % 253);
                uint8_t smp = el[pos];
                if (smp != smp_exp) 
                  { bug("set/get sample mismatch smp = %u smp_exp = %u", smp, smp_exp); } 
            }
}

void test_enum(desc_t *A)
{
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_enum (d = %d nt = %llu)...\n", A->d, nt);
  
  ix_dim_t d = A->d;

  /* Get the size vector, extend with 1's: */
  ix_size_t sz[MAXDIM];
  extend_size_vector(A->d, A->z, MAXDIM, sz);

  /* Fabricate the other two arrays: */
  desc_t *B = make_desc(d, sz);
  desc_t *C = make_desc(d, sz);
  bool_t verbose = FALSE;
  if (verbose) 
    { print_desc(stderr, "A = { ", A, " }\n");
      print_desc(stderr, "B = { ", B, " }\n");
      print_desc(stderr, "C = { ", C, " }\n");
    }
    
  ix_index_t ix[MAXDIM];

  /* Allocate a vector to store all positions of {A}: */
  ix_pos_count_t npos = ix_num_tuples(d, sz);
  ix_pos_t *posA = notnull(malloc(npos*sizeof(ix_pos_t)), "no mem");
  ix_pos_t *posB = notnull(malloc(npos*sizeof(ix_pos_t)), "no mem");
  ix_pos_t *posC = notnull(malloc(npos*sizeof(ix_pos_t)), "no mem");
  int64_t kpos;     /* Current index into {posA}. */

  auto void check_fw( ix_index_t ix[], ix_pos_t pA, ix_pos_t pB, ix_pos_t pC );
  auto void check_bw( ix_index_t ix[], ix_pos_t pA, ix_pos_t pB, ix_pos_t pC );
    /* These are procedures that can be passed as the {op} argument of 
      {ix_enum}. They check whether the given positions {pA,pB,pC}
      match {posA[kpos],posB[kpos],posC[kpos]}.  The procedure 
      {check_fw} increments {kpos} after the test,
      while {check_bw} decrements it before the test. */

  /* Enumerate all positions by hand in "F" order: */
  kpos = 0;
  assert(MAXDIM == 6);
  for (ix[5] = 0; ix[5] < sz[5]; ix[5]++)
    for (ix[4] = 0; ix[4] < sz[4]; ix[4]++)
      for (ix[3] = 0; ix[3] < sz[3]; ix[3]++)
        for (ix[2] = 0; ix[2] < sz[2]; ix[2]++)
          for (ix[1] = 0; ix[1] < sz[1]; ix[1]++)
            for (ix[0] = 0; ix[0] < sz[0]; ix[0]++)
              { posA[kpos] = ix_position(A->d, ix, A->b, A->s);
                posB[kpos] = ix_position(B->d, ix, B->b, B->s);
                posC[kpos] = ix_position(C->d, ix, C->b, C->s);
                kpos++;
              }
  assert(kpos == npos);
  if (verbose) { fprintf(stderr, "  counted %llu tuples\n", npos); }
  
  /* Check forward enumeration, Fortran order: */
  kpos = 0;
  ix_enum(check_fw, d, sz, ix_order_F, FALSE, A->b, A->s, B->b, B->s, C->b, C->s);
  assert(kpos == npos);
  
  /* Check backward enumeration, Fortran order: */
  kpos = npos;
  ix_enum(check_bw, d, sz, ix_order_F, TRUE,  A->b, A->s, B->b, B->s, C->b, C->s);
  assert(kpos == 0);
  
  /* Enumerate all positions by hand in "L" order: */
  kpos = 0;
  assert(MAXDIM == 6);
  for (ix[0] = 0; ix[0] < sz[0]; ix[0]++)
    for (ix[1] = 0; ix[1] < sz[1]; ix[1]++)
      for (ix[2] = 0; ix[2] < sz[2]; ix[2]++)
        for (ix[3] = 0; ix[3] < sz[3]; ix[3]++)
          for (ix[4] = 0; ix[4] < sz[4]; ix[4]++)
            for (ix[5] = 0; ix[5] < sz[5]; ix[5]++)
              { posA[kpos] = ix_position(A->d, ix, A->b, A->s);
                posB[kpos] = ix_position(B->d, ix, B->b, B->s);
                posC[kpos] = ix_position(C->d, ix, C->b, C->s);
                kpos++;
              }
  assert(kpos == npos);
  
  /* Check forward enumeration, C/Pascal order: */
  kpos = 0;
  ix_enum(check_fw, d, sz, ix_order_L, FALSE, A->b, A->s, B->b, B->s, C->b, C->s);
  assert(kpos == npos);
  
  /* Check backward enumeration, C/Pascal order: */
  kpos = npos;
  ix_enum(check_bw, d, sz, ix_order_L, TRUE,  A->b, A->s, B->b, B->s, C->b, C->s);
  assert(kpos == 0);

  /* IMPLEMENTATION OF LOCAL PROCS */

  void check_fw( ix_index_t ix[], ix_pos_t pA, ix_pos_t pB, ix_pos_t pC )
    {
      assert(kpos >= 0);
      assert(kpos < npos);
      assert(pA == posA[kpos]);
      assert(pB == posB[kpos]);
      assert(pC == posC[kpos]);
      kpos++;
    }
  
  void check_bw( ix_index_t ix[], ix_pos_t pA, ix_pos_t pB, ix_pos_t pC )
    {
      kpos--;
      assert(kpos >= 0);
      assert(kpos < npos);
      assert(pA == posA[kpos]);
      assert(pB == posB[kpos]);
      assert(pC == posC[kpos]);
    }
}

void test_packed(desc_t *A)
{
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_packed_XXX ops (d = %d nt = %llu)...\n", A->d, nt);
  
  /* Get number of indices of {A}: */
  ix_dim_t d = A->d;
  int i;

  /* Get the size vector, extend with 1's: */
  ix_size_t sz[MAXDIM];
  extend_size_vector(A->d, A->z, MAXDIM, sz);
  
  /* Fabricate a *packed* array (can't use {A}, but use its size): */
  desc_t *F = make_packed_desc(d, A->z, ix_order_F);
  desc_t *L = make_packed_desc(d, A->z, ix_order_L);
  bool_t verbose = FALSE;
  if (verbose) 
    { print_desc(stderr, "A = { ", A, " }\n");
      print_desc(stderr, "F = { ", F, " }\n");
      print_desc(stderr, "L = { ", L, " }\n");
    }
    
  ix_index_t ix[MAXDIM], jx[MAXDIM];

  /* Allocate a vector to store all positions of {B}: */
  ix_pos_count_t npos = ix_num_tuples(d, sz);
  int64_t kpos;     /* Current index into {posA}. */

  /* Enumerate all tuples by hand in "F" order, check {F} positions: */
  kpos = 0;
  assert(MAXDIM == 6);
  for (ix[5] = 0; ix[5] < sz[5]; ix[5]++)
    for (ix[4] = 0; ix[4] < sz[4]; ix[4]++)
      for (ix[3] = 0; ix[3] < sz[3]; ix[3]++)
        for (ix[2] = 0; ix[2] < sz[2]; ix[2]++)
          for (ix[1] = 0; ix[1] < sz[1]; ix[1]++)
            for (ix[0] = 0; ix[0] < sz[0]; ix[0]++)
              { ix_pos_t posF1 = ix_packed_position(d, ix, sz, ix_order_F);
                ix_pos_t posF2 = ix_position(d, ix, F->b, F->s);
                if (verbose)
                  { fprintf(stderr, "  kpos = %llu\n", kpos);  
                    fprintf(stderr, "  posF1 = %llu\n", posF1); 
                    fprintf(stderr, "  posF2 = %llu\n", posF2);
                  }
                assert(posF1 == kpos);
                assert(posF2 == posF1);
                ix_packed_indices(d, kpos, sz, ix_order_F, jx);
                for (i = 0; i < d; i++) { assert(ix[i] == jx[i]); }
                kpos++;
              }
  assert(kpos == npos);

  /* Enumerate all tuples by hand in "L" order, check {L} positions: */
  kpos = 0;
  assert(MAXDIM == 6);
  for (ix[0] = 0; ix[0] < sz[0]; ix[0]++)
    for (ix[1] = 0; ix[1] < sz[1]; ix[1]++)
      for (ix[2] = 0; ix[2] < sz[2]; ix[2]++)
        for (ix[3] = 0; ix[3] < sz[3]; ix[3]++)
          for (ix[4] = 0; ix[4] < sz[4]; ix[4]++)
            for (ix[5] = 0; ix[5] < sz[5]; ix[5]++)
              { ix_pos_t posL1 = ix_packed_position(d, ix, sz, ix_order_L);
                ix_pos_t posL2 = ix_position(d, ix, L->b, L->s);
                assert(posL1 == kpos);
                assert(posL2 == posL1);
                ix_packed_indices(d, kpos, sz, ix_order_L, jx);
                for (i = 0; i < d; i++) { assert(ix[i] == jx[i]); }
                kpos++;
              }
  assert(kpos == npos);

}

void test_crop(desc_t *A)
{
  if (A->d == 0) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_crop (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  ix_index_t sh[MAXDIM]; /* Shift of cropped array along each axis. */
  ix_size_t sz[MAXDIM]; /* Size of cropped array along each axis. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {sh} and {sz} for full array: */
      { int i; for (i = 0; i < MAXDIM; i++) { sh[i] = 0; } }
      extend_size_vector(A->d, A->z, MAXDIM, sz);

      /* Start with the standard array: */
      desc_t B = (*A);
      if (verbose) { print_desc(stderr, "B = { ", &B, " }\n"); } 

      /* Perform first cropping: */
      ix_axis_t ax1 = (trial % A->d);
      ix_size_t skip1 = (trial & 3)*sz[ax1]/8;
      ix_size_t keep1 = (trial/4 & 3)*(sz[ax1] - skip1)/8;
      if (keep1 < 0) { keep1 = 0; }
      sh[ax1] += skip1; sz[ax1] = keep1;
      ix_crop(B.d, B.z, &B.b, B.s, ax1, skip1, keep1);
      if (verbose) 
        { fprintf(stderr, "crop(%d, %llu, %llu)", ax1, skip1, keep1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { print_desc(stderr, "A = { ", A, " }\n");
          fprintf(stderr, "crop(%d, %llu, %llu)", ax1, skip1, keep1);
          print_desc(stderr, "B = { ", &B, " }\n");
          fprintf(stderr, "crop(%d, %llu, %llu)", ax1, skip1, keep1);
          bug("ix_crop returns an invalid desc");
        }
      /* Check whether the size of {B} is correct: */
      check_sizes(B.d, B.z, sz);
      
      /* Perform second cropping (possibly on same axis): */ 
      ix_axis_t ax2 = ((trial / A->d) % A->d);
      ix_size_t skip2 = (trial/2 & 3)*sz[ax2]/8;
      ix_size_t keep2 = (trial/5 & 3)*(sz[ax2] - skip2)/8;
      if (keep2 < 0) { keep2 = 0; }
      sh[ax2] += skip2; sz[ax2] = keep2;
      ix_crop(B.d, B.z, &B.b, B.s, ax2, skip2, keep2);
      if (verbose) 
        { fprintf(stderr, "crop(%d, %llu, %llu)", ax2, skip2, keep2);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_crop returns an invalid desc"); }
      
      /* Check whether the size of {B} is correct: */
      check_sizes(B.d, B.z, sz);

      /* Now check coincidence of the two arrays: */

      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);
      for (ixB[5] = 0; ixB[5] < sz[5]; ixB[5]++)
        for (ixB[4] = 0; ixB[4] < sz[4]; ixB[4]++)
          for (ixB[3] = 0; ixB[3] < sz[3]; ixB[3]++)
            for (ixB[2] = 0; ixB[2] < sz[2]; ixB[2]++)
              for (ixB[1] = 0; ixB[1] < sz[1]; ixB[1]++)
                for (ixB[0] = 0; ixB[0] < sz[0]; ixB[0]++)
                  { 
                    ix_assign(A->d, ixA, ixB);
                    ix_shift(A->d, ixA, sh);
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("crop position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}

void test_subsample(desc_t *A) 
{
  if (A->d == 0) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_subsample (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  ix_size_t st[MAXDIM]; /* Subsampling step along each axis. */
  ix_size_t sz[MAXDIM]; /* Size of subsampled array along each axis. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {st} and {sz} for full array: */
      { int i; for (i = 0; i < MAXDIM; i++) { st[i] = 1; } }
      extend_size_vector(A->d, A->z, MAXDIM, sz);
      /* Start with the standard array: */
      desc_t B = (*A);
      
      /* Perform first subsampling: */
      ix_axis_t ax1 = (trial % A->d);
      ix_size_t step1 = 1 + (trial & 3)*(sz[ax1] - 1)/8;
      ix_size_t keep1 = (sz[ax1] + step1 - 1)/step1;
      st[ax1] *= step1; sz[ax1] = keep1;
      ix_subsample(B.d, B.z, &(B.b), B.s, ax1, step1);
      if (verbose) 
        { fprintf(stderr, "subsample(%d, %llu)", ax1, step1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE)) 
        { bug("ix_flip returns an invalid desc"); }
      
      /* Perform second subsampling (possibly on same axis): */ 
      ix_axis_t ax2 = ((trial / A->d) % A->d);
      ix_size_t step2 = 1 + (trial/2 & 3)*(sz[ax2] - 1)/8;
      ix_size_t keep2 = (sz[ax2] + step2 - 1)/step2;
      st[ax2] *= step2; sz[ax2] = keep2;
      ix_subsample(B.d, B.z, &(B.b), B.s, ax2, step2);
      if (verbose) 
        { fprintf(stderr, "subsample(%d, %llu)", ax2, step2);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE)) 
        { bug("ix_subsample returns an invalid desc"); }
      
      check_sizes(B.d, B.z, sz); 
      /* Now check coincidence of the two arrays: */
      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);
      for (ixB[5] = 0; ixB[5] < sz[5]; ixB[5]++)
        for (ixB[4] = 0; ixB[4] < sz[4]; ixB[4]++)
          for (ixB[3] = 0; ixB[3] < sz[3]; ixB[3]++)
            for (ixB[2] = 0; ixB[2] < sz[2]; ixB[2]++)
              for (ixB[1] = 0; ixB[1] < sz[1]; ixB[1]++)
                for (ixB[0] = 0; ixB[0] < sz[0]; ixB[0]++)
                  { 
                    ix_axis_t ax;
                    for (ax = 0; ax < A->d; ax++) { ixA[ax] = ixB[ax]*st[ax]; }
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("subsample position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}

void test_swap_indices(desc_t *A)
{
  if (A->d == 0) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_swap_indices (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  ix_axis_t tr[MAXDIM]; /* Axis {i} of {B} is axis {tr[i]} of {A}. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {tr} for full array: */
      { int i; for (i = 0; i < MAXDIM; i++) { tr[i] = i; } }
      /* Start with the standard array: */
      desc_t B = (*A);

      /* Perform first tranpose: */
      ix_axis_t ax1 = (trial % A->d), bx1 = (trial/2 % A->d);
      { ix_axis_t tmp = tr[ax1]; tr[ax1] = tr[bx1]; tr[bx1] = tmp; }
      ix_swap_indices(B.d, B.z, &(B.b), B.s, ax1, bx1);
      if (verbose) 
        { fprintf(stderr, "swap_indices(%d, %d)", ax1, bx1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_flip returns an invalid desc"); }

      /* Perform second tranpose (possibly on same axis): */ 
      ix_axis_t ax2 = ((trial/A->d) % A->d), bx2 = ((trial/10) % A->d);
      { ix_axis_t tmp = tr[ax2]; tr[ax2] = tr[bx2]; tr[bx2] = tmp; }
      ix_swap_indices(B.d, B.z, &(B.b), B.s, ax2, bx2);
      if (verbose) 
        { fprintf(stderr, "swap_indices(%d, %d)", ax2, bx2);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_swap_indices returns an invalid desc"); }
      
      /* Check whether the size and increments of {B} is correct: */
      { int i;
        for (i = 0; i < A->d; i++)
          { if (B.z[i] != A->z[tr[i]])
              { bug("B.z[%d] = %llu  A.z[%d] = %llu", i, B.z[i], tr[i], A->z[tr[i]]); }
            if (B.s[0] != A->s[tr[0]])
              { bug("B.s[%d] = %llu  A.s[%d] = %llu", i, B.s[i], tr[i], A->s[tr[i]]); }
          }
      }
      /* Get {A}'s size vector, extend with 1's: */
      ix_size_t sz[MAXDIM];
      extend_size_vector(A->d, A->z, MAXDIM, sz);

      /* Now check coincidence of the two arrays: */
      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);

      for (ixA[5] = 0; ixA[5] < sz[5]; ixA[5]++)
        for (ixA[4] = 0; ixA[4] < sz[4]; ixA[4]++)
          for (ixA[3] = 0; ixA[3] < sz[3]; ixA[3]++)
            for (ixA[2] = 0; ixA[2] < sz[2]; ixA[2]++)
              for (ixA[1] = 0; ixA[1] < sz[1]; ixA[1]++)
                for (ixA[0] = 0; ixA[0] < sz[0]; ixA[0]++)
                  { 
                    ix_axis_t ax;
                    for (ax = 0; ax < A->d; ax++) { ixB[ax] = ixA[tr[ax]]; }
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("swap_indices position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}

void test_flip_indices(desc_t *A)
{
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "NOT Checking ix_flip_indices (d = %d nt = %llu)...\n", A->d, nt);
}

void test_flip(desc_t *A)
{
  if (A->d == 0) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_flip (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  bool_t fp[MAXDIM]; /* Tells whether each axis was flipped or not. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {fp} for unflipped array: */
      { int i; for (i = 0; i < MAXDIM; i++) { fp[i] = FALSE; } }
      /* Start with the standard array: */
      desc_t B = (*A);

      /* Perform first flip: */
      ix_axis_t ax1 = (trial % A->d);
      fp[ax1] = ! fp[ax1];
      ix_flip(B.d, B.z, &(B.b), B.s, ax1);
      if (verbose) 
        { fprintf(stderr, "flip(%d)", ax1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_flip returns an invalid desc"); }

      /* Perform second flip (possibly on same axis): */ 
      ix_axis_t ax2 = ((trial / A->d) % A->d);
      fp[ax2] = ! fp[ax2];
      ix_flip(B.d, B.z, &(B.b), B.s, ax2);
      if (verbose) 
        { fprintf(stderr, "flip(%d)", ax2);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_flip returns an invalid desc"); }

      /* Check whether the size of {B} is correct: */
      check_sizes(B.d, B.z, A->z);
      
      /* Check whether the increments of {B} are correct: */
      { int i;
        for (i = 0; i < A->d; i++)
          { if (B.s[i] != A->s[i]*(fp[i]?-1:+1))
              { bug("step[%d] = %llu", i, B.s[i]); }
          }
      }
      
      /* Get {B}'s size vector, extend with 1's: */
      ix_size_t sz[MAXDIM];
      extend_size_vector(B.d, B.z, MAXDIM, sz);

      /* Now check coincidence of the two arrays: */
      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);
      for (ixB[5] = 0; ixB[5] < sz[5]; ixB[5]++)
        for (ixB[4] = 0; ixB[4] < sz[4]; ixB[4]++)
          for (ixB[3] = 0; ixB[3] < sz[3]; ixB[3]++)
            for (ixB[2] = 0; ixB[2] < sz[2]; ixB[2]++)
              for (ixB[1] = 0; ixB[1] < sz[1]; ixB[1]++)
                for (ixB[0] = 0; ixB[0] < sz[0]; ixB[0]++)
                  { 
                    ix_axis_t ax;
                    for (ax = 0; ax < A->d; ax++) 
                      { ixA[ax] = (fp[ax] ? A->z[ax] - 1 - ixB[ax] : ixB[ax]); }
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("flip position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}

void test_diagonal(desc_t *A)
{
  if (A->d < 2) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_diagonal (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  ix_size_t sz[MAXDIM]; /* Size of diagonalized array along each axis. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {sz} for full array: */
      extend_size_vector(A->d, A->z, MAXDIM, sz);
      /* Start with the standard array: */
      desc_t B = (*A);

      /* Choose any two axes: */
      ix_axis_t ax1 = (trial % A->d), bx1 = (trial/A->d % A->d);
      if (bx1 == ax1) { bx1 = (ax1 + 1) % A->d; }
      /* Make it square: */
      ix_size_t szmin = (sz[ax1] < sz[bx1] ? sz[ax1] : sz[bx1]);
      ix_crop(B.d, B.z, &(B.b), B.s, ax1, 0, szmin); sz[ax1] = szmin;
      ix_crop(B.d, B.z, &(B.b), B.s, bx1, 0, szmin); sz[bx1] = szmin;
      /* Diagonalize array: */
      ix_diagonal(B.d, B.z, &(B.b), B.s, ax1, bx1);
      /* Compute expected sizes of {1}: */
      if (sz[bx1] < sz[ax1]) { sz[ax1] = sz[bx1]; }
      if (bx1 != ax1) { sz[bx1] = 1; }
      if (verbose) 
        { fprintf(stderr, "diagonal(%d, %d)", ax1, bx1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_diagonal returns an invalid desc"); }

      /* Check whether the size of {1} is correct: */
      check_sizes(B.d, B.z, sz);
      /* Now check coincidence of the two arrays: */
      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);
      for (ixB[5] = 0; ixB[5] < sz[5]; ixB[5]++)
        for (ixB[4] = 0; ixB[4] < sz[4]; ixB[4]++)
          for (ixB[3] = 0; ixB[3] < sz[3]; ixB[3]++)
            for (ixB[2] = 0; ixB[2] < sz[2]; ixB[2]++)
              for (ixB[1] = 0; ixB[1] < sz[1]; ixB[1]++)
                for (ixB[0] = 0; ixB[0] < sz[0]; ixB[0]++)
                  { 
                    ix_assign(A->d, ixA, ixB);
                    ixA[bx1] = ixA[ax1];
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("diagonal position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}

void test_chop(desc_t *A)
{
  if (A->d < 2) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_chop (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  ix_size_t sz[MAXDIM]; /* Size of chopped array along each axis. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {sz} for full array: */
      extend_size_vector(A->d, A->z, MAXDIM, sz);
      /* Start with the standard array: */
      desc_t B = (*A);

      /* Choose two DISTINCT axes: */
      ix_axis_t ax1 = (trial % A->d);
      ix_axis_t bx1 = (trial/A->d % (A->d - 1));
      if (bx1 >= ax1) { bx1++; }
      /* Choose a nonzero chunk size: */
      ix_size_t chunksz1 = 1 + (trial/3 & 3)*((A->z[ax1] <= 1 ? 1 : A->z[ax1]) + 7)/8;
      /* Ensure integral number of chunks: */
      ix_size_t szround = (sz[ax1]/chunksz1)*chunksz1;
      ix_crop(B.d, B.z, &(B.b), B.s, ax1, 0, szround); sz[ax1] = szround;
      /* Ensure that the argument of {chop} is trivial along axis {bx1}: */
      if (sz[bx1] == 0) { /* Skip this trial: */ continue; }
      ix_crop(B.d, B.z, &(B.b), B.s, bx1, 0, 1); sz[bx1] = 1;
      /* Chop array: */
      ix_chop(B.d, B.z, &(B.b), B.s, ax1, chunksz1, bx1);
      /* Compute expected sizes of {1}: */
      sz[bx1] = sz[ax1]/chunksz1; sz[ax1] = chunksz1;
      if (verbose) 
        { fprintf(stderr, "chop(%d, %lld, %d)", ax1, chunksz1, bx1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_chop returns an invalid desc"); }

      /* Check whether the size of {1} is correct: */
      check_sizes(B.d, B.z, sz);
      /* Now check coincidence of the two arrays: */
      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);
      for (ixB[5] = 0; ixB[5] < sz[5]; ixB[5]++)
        for (ixB[4] = 0; ixB[4] < sz[4]; ixB[4]++)
          for (ixB[3] = 0; ixB[3] < sz[3]; ixB[3]++)
            for (ixB[2] = 0; ixB[2] < sz[2]; ixB[2]++)
              for (ixB[1] = 0; ixB[1] < sz[1]; ixB[1]++)
                for (ixB[0] = 0; ixB[0] < sz[0]; ixB[0]++)
                  { 
                    ix_assign(A->d, ixA, ixB);
                    ixA[ax1] = ixA[ax1] + chunksz1 * ixB[bx1];
                    ixA[bx1] = 0;
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("chop position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}
  
void test_replicate(desc_t *A)
{
  if (A->d == 0) { return; }
  ix_pos_count_t nt = ix_num_tuples(A->d, A->z);
  fprintf(stderr, "Checking ix_replicate (d = %d nt = %llu)...\n", A->d, nt);
  bool_t verbose = FALSE;
  if (verbose) { print_desc(stderr, "A = { ", A, " }\n"); } 
  int trial, ntrials = 8*A->d + 1;
  ix_size_t sz[MAXDIM]; /* Size of replicated array along each axis. */
  for (trial = 0; trial < ntrials; trial++)
    { /* Initialze {sz} for full array: */
      extend_size_vector(A->d, A->z, MAXDIM, sz);
      /* Start with the standard array: */
      desc_t B = (*A);
      
      /* Choose axis for first replication: */
      ix_axis_t ax1 = (trial % A->d);
      if (sz[ax1] == 0) { continue; }
      /* Choose a positive replication factor: */
      ix_size_t rep1 = 1 + (trial/3 & 3)*A->d;
      /* Make trivial along axis {ax1}: */
      ix_size_t skip1 = (B.z[ax1]-1)/2;
      ix_crop(B.d, B.z, &(B.b), B.s, ax1, skip1, 1); sz[ax1] = 1;
      /* Replicate along axis {ax1}: */
      ix_replicate(B.d, B.z, &(B.b), B.s, ax1, rep1);
      sz[ax1] = rep1;
      if (verbose) 
        { fprintf(stderr, "chop(%d,%llu,1)+replicate(%d,%lld)", ax1, skip1, ax1, rep1);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_replicate returns an invalid desc"); }
      
      /* Choose axis for second replication (may be {ax1}): */
      ix_axis_t ax2 = (trial/A->d % A->d);
      if (sz[ax2] == 0) { continue; }
      /* Choose a replication factor: */
      ix_size_t rep2 = 1 + (trial/8 & 3)*A->d;
      /* Make trivial along axis {ax2}: */
      ix_size_t skip2 = (B.z[ax2]-1)/2;
      ix_crop(B.d, B.z, &(B.b), B.s, ax2, skip2, 1); sz[ax2] = 1;
      /* Replicate along axis {ax2}: */
      ix_replicate(B.d, B.z, &(B.b), B.s, ax2, rep2);
      sz[ax2] = rep2;
      if (verbose) 
        { fprintf(stderr, "chop(%d,%llu,1)+replicate(%d,%lld)", ax2, skip2, ax2, rep2);
          print_desc(stderr, " = { ", &B, " }\n");
        } 
      /* Check validity of descriptor: */
      if (! ix_parms_are_valid(B.d, B.z, B.b, B.s, FALSE))
        { bug("ix_replicate returns an invalid desc"); }

      /* Check whether the size of {1} is correct: */
      check_sizes(B.d, B.z, sz);
      /* Now check coincidence of the two arrays: */
      ix_index_t ixA[MAXDIM], ixB[MAXDIM];
      assert(MAXDIM == 6);
      for (ixB[5] = 0; ixB[5] < sz[5]; ixB[5]++)
        for (ixB[4] = 0; ixB[4] < sz[4]; ixB[4]++)
          for (ixB[3] = 0; ixB[3] < sz[3]; ixB[3]++)
            for (ixB[2] = 0; ixB[2] < sz[2]; ixB[2]++)
              for (ixB[1] = 0; ixB[1] < sz[1]; ixB[1]++)
                for (ixB[0] = 0; ixB[0] < sz[0]; ixB[0]++)
                  { 
                    ix_assign(A->d, ixA, ixB);
                    ixA[ax1] = skip1;
                    /* The second crop changes {b} only if {ax1 != ax2}: */
                    if (ax2 != ax1) { ixA[ax2] = skip2; }
                    ix_pos_t posB = ix_position(B.d, ixB, B.b, B.s);
                    ix_pos_t posA = ix_position(A->d, ixA, A->b, A->s);
                    if (posB != posA)
                      { bug("replicate position error posB = %llu posA = %llu", posB, posA); } 
                  }
    }
}
  
void test_reduce(void)
{
  fprintf(stderr, "Checking ix_reduce...\n");
  bool_t verbose = FALSE;
  /* Precomputed examples: */
  ix_size_t N = 6;
  ix_index_t iMin = -3;
  ix_index_t iMax = 15;
  char *tseq[ix_reduction_LAST + 1];
  tseq[ix_reduction_SINGLE] = "*,*,*,0,1,2,3,4,5,*,*,*,*,*,*,*,*,*,*,";
  tseq[ix_reduction_EXTEND] = "0,0,0,0,1,2,3,4,5,5,5,5,5,5,5,5,5,5,5,";
  tseq[ix_reduction_REPEAT] = "3,4,5,0,1,2,3,4,5,0,1,2,3,4,5,0,1,2,3,";
  tseq[ix_reduction_MIRROR] = "2,1,0,0,1,2,3,4,5,5,4,3,2,1,0,0,1,2,3,";
  tseq[ix_reduction_PXMIRR] = "3,2,1,0,1,2,3,4,5,4,3,2,1,0,1,2,3,4,5,";

  if (verbose) { fprintf(stderr, "N = %lld\n", N); } 
  ix_reduction_t red;
  for (red = ix_reduction_FIRST; red <= ix_reduction_LAST; red++)
    { if (verbose) { fprintf(stderr, "red = %d (%c)\n", red, "SERMP"[red]); } 
      ix_index_t i;
      for (i = iMin; i <= iMax; i++)
        { char cr = tseq[red][2*(i-iMin)];
          ix_index_t jExp = (cr == '*' ? -1 : cr - '0');
          ix_index_t j = ix_reduce(i, N, red);
          if (verbose) { fprintf(stderr, "%3lld --> %3lld  (should be %3lld)\n", i, j, jExp); } 
          demand(((red == ix_reduction_SINGLE) && (j == -1)) || ((j >= 0) && (j < N)), "invalid {ix_reduce} result");
          demand(j == jExp, "error in {ix_reduce}");
        }
    }
    
}

void extend_size_vector(ix_dim_t dA, ix_size_t szA[], ix_dim_t dB, ix_size_t szB[])
{
  int i;
  for (i = 0; i < dB; i++) { szB[i] = (i < dA ? szA[i] : 1); }
}

void check_max_indices(ix_dim_t d, ix_index_t ixa[], ix_index_t ixb[])
{
  ix_axis_t i;
  for (i = 0; i < d; i++)
    { if (ixa[i] != ixb[i]) 
        { bug("(ixa[%d] = %lld ixb[%d] = %lld", i, ixa[i], i, ixb[i]); }
    }
}
  
void check_max_indices_sz(ix_dim_t d, ix_index_t ixa[], ix_size_t szb[])
{
  ix_axis_t i;
  for (i = 0; i < d; i++)
    { ix_index_t ixbi = szb[i]-1;
      if (ixa[i] != ixbi) 
        { bug("(ixa[%d]-1) = %lld ixb[%d] = %lld", i, ixa[i], i, ixbi); }
    }
}
  
void check_sizes(ix_dim_t d, ix_size_t sza[], ix_size_t szb[])
{
  ix_axis_t i;
  for (i = 0; i < d; i++)
    { if (sza[i] != szb[i]) 
        { bug("sza[%d] = %llu szb[%d] = %llu", i, sza[i], i, szb[i]); }
    }
}

void print_desc(FILE *wr, char *pf, desc_t *A, char *sf)
  { fprintf(wr, "%s", pf); 
    int i;
    for (i = 0; i < A->d; i++)
      { fprintf(wr, " %llu(×%lld)", A->z[i], A->s[i]); }
    fprintf(wr, 
      " b = %llu tuples = %llu positions = %llu", 
      A->b, 
      ix_num_tuples(A->d,A->z), 
      ix_num_positions(A->d,A->z,A->s)
    );
    fprintf(wr, "%s", sf);
  }

void print_indices(FILE *wr, char *pf, ix_dim_t d, ix_index_t ix[], char *sf)
  { fprintf(wr, "%s", pf); 
    int i;
    for (i = 0; i < d; i++)
      { fprintf(wr, " %2lld", ix[i]); }
    fprintf(wr, "%s", sf);
  }