#include <stdio.h>

// Large math problems running on massively parallel systems sometimes use a lot
// of upper triangular matrices.  Implemented naively, these waste 50% of the memory
// in the machine, which is not recoverable by virtual memory techniques because it
// is interspersed with data on each row.  By mapping the array elements into an
// array that is half the size and not actually storing the zeroes, we can do twice
// the computation in the same machine or use half as many machines in total.

// To implement a safety feature of making the zero-part of an upper triangular matrix
// read-only, we place all the zeroes in write-protected memory and cause a memory violation
// if the programmer attempts to write to them.  System dependent but relatively portable.
// Requires that you compile with the -Wno-discarded-qualifiers option.

// for the awkward case (an even-sized array bound):

//                         +--------/
//     row  0, 40 items -> |0      /
//     row  1, 39 items -> |      /
//     row 19, 21 items -> |     /
//     row 20, 20 items -> |----/ <------  cut and rotate here to form a rectangle.
//     row 21, 19 items -> |   /
//                         |  /
//     row 39,  1 item  -> | /
//     row 40,  0 items -> |/
//                         /


//    x y   x  y
//    0,0   39,0
//     +----/           |                     +--------/
//     |   /            | row  0, 40 items -> |0      /| <-- row 40, 0 items
//     |  / - 20,18     | row  1, 39 items -> |      /0| <-- row 39, 1 item
//     | /\             | row 19, 21 items -> |     /  | <-- row 21, 19 items
//     |/  19,19        | row 20, 20 items -> |    /???| <-- row 20, 20 items  half of row 20 is wasted...
//0,39 v                |                     ~~~~~~~~~~
//     |                |

// for odd-sized array bounds, there is no need for the wasted half-row marked '???' above...

// And once the mapping above is done, mirror the indexes in x to get a proper Upper Triangular Matrix which
// looks like this...
//     ____
//     \  |
//      \ |
//       \|
//

// Rather than store the underlying data in a 2D array, if we use a 1-D array,
// and map the indexes ourselves, it is possible to recover that final half-row...

// The implementation allows for the matrix elements to be any scalar type.

#define DECLARE_TRIANGULAR_MATRIX(type, name, bound, zero)			\
  type _##name[bound * (bound+1) / 2 + 1]; /* +1 is for a debugging tombstone */ \
  type *__##name(int x, int y) { \
    static const type Zero = zero; /* writing to the lower half of the matrix will segfault */ \
    x = (bound-1)-x; /* mirror */ \
    if (x+y >= bound) return &Zero; /* requires cc -Wno-discarded-qualifiers */ \
    if (y > bound/2) {x = (bound-1)-x; y = bound-y;} \
    return &_##name[y*bound+x]; /* standard mapping of x,y -> X */ \
  }
#define TRIANGULAR_MATRIX(name, x, y)  *__##name(x,y)


// ----------------------------------------------------------------------------------------


// Simulate 'int fred[11][11];' as an upper triangular matrix: 
#define ARRAYSIZE 11
DECLARE_TRIANGULAR_MATRIX(int, fred, ARRAYSIZE, 0)
#define fred(x, y) TRIANGULAR_MATRIX(fred, x, y)
/* unfortunately we can't #define fred[x][y] here ... In the Imp language which used () both
   for array indexes and procedure parameters, we could write a mapping function fred(x,y)
   which made the indirected function call indistinguishable from a normal array access.
   We attempt to do something similar here using macros, but C is not as cooperative. */



int main(int argc, char **argv) {
  int x,y, element;

  // treat as fully populated 2D array...
  for (y = 0; y < ARRAYSIZE; y++) {
    for (x = 0; x < ARRAYSIZE; x++) {
      if (ARRAYSIZE-x+y <= ARRAYSIZE) fred(x,y) = (x+1) * 100 + (y+1); // upper triangle test
    }
  }

  fprintf(stdout, "Upper Triangular matrix:\n\n");
  fprintf(stdout, "    ");
  for (x = 0; x < ARRAYSIZE; x++) fprintf(stdout, "%5d", x);
  fprintf(stdout, "\n    ");
  for (x = 0; x < ARRAYSIZE; x++) fprintf(stdout, "_____");
  fprintf(stdout, "\n");
  for (y = 0; y < ARRAYSIZE; y++) {
    fprintf(stdout, "%2d |", y);
    for (x = 0; x < ARRAYSIZE; x++) {
      element = fred(x,y);
      fprintf(stdout, "%5d", element);
      if (ARRAYSIZE-x+y <= ARRAYSIZE) { // upper triangle test
        if (element # (x+1) * 100 + (y+1)) {
          fflush(stdout); fprintf(stderr, "Mismatch! at %d,%d (%d # %d)\n", x, y, element, x * 100 + y);
        }
      } else if (element # 0) {
        fflush(stdout); fprintf(stderr, "Mismatch! at %d,%d (%d # 0)\n", x, y, element);
      }
    }
    fprintf(stdout, "\n");
  }

  return 0;
}


#ifdef IMP
begin
  integer name Force an error  { being uninitialised, Imp will raise a signal on access through the pointer... }

  const integer ARRAY SIZE = 11

  ! for compatibility with the C version we will base our array at 0
  ! even though in Imp we would probably normally base it at 1...

  integer array triangular fred(0 : ARRAY SIZE * (ARRAY SIZE+1)//2 + 1)
  ! Actually, in imp, we would more likely declare the array as a two dimensional array
  ! similar to: %integer %array triangular fred(0 : ARRAY SIZE, 0 : (ARRAY SIZE+1)//2)

  integer map fred(integer x, y)
    x = (ARRAY SIZE-1) - x
    result == Force an error if x+y >= bound
    if y > ARRAY SIZE//2 start
      x = (ARRAY SIZE-1) - x; y = ARRAY SIZE - y
    finish
    result == triangular fred(y*ARRAY SIZE + x)
    ! or if we had used a 2d array, ... %result == triangular fred(x, y)
  end

  integer x, y, element

  select output(0)

  ! treat as fully populated 2D array...
  for y = 0, 1, ARRAY SIZE-1 cycle
    for x = 0, 1, ARRAY SIZE-1 cycle
      if ARRAY SIZE-x+y <= ARRAY SIZE start
	fred(x,y) = (x+1) * 100 + (y+1)    { upper triangle test }
      finish
    repeat
  repeat

  print string("Upper Triangular matrix:".snl.snl)
  spaces(4)
  write(x, 5) for (x = 0; x < ARRAY SIZE; x++)
  newline; spaces(4);
  print string("_____") for x = 0, 1, ARRAY SIZE-1
  newline
  for y = 0, 1, ARRAY SIZE-1 cycle
    write(y, 2); print string(" |")
    for x = 0, 1, ARRAY SIZE-1 cycle
      element = fred(x,y)
      write(element, 5);
      if ARRAY SIZE-x+y <= ARRAY SIZE start   { upper triangle test }
        if element # (x+1) * 100 + (y+1) start
          print string("Mismatch! at ")
          write(x, 0); print symbol(','); write(y, 0)
	  print string(" ("); write(element, 0); print string(" # ")
	  write(x * 100 + y, 0); print symbol(')'); newline
        finish
      else if element # 0
        print string("Mismatch! at ")
        write(x, 0); print symbol(','); write(y, 0)
	print string(" ("); write(element, 0); print string(" # 0)".snl)
      finish
    repeat
    newline
  repeat
end of program
#endif