#include <stdlib.h>
#include <math.h>
#include <iostream>
#include "SuperTensor.h"

// The following macro defines a family of functions that work with 3D
// arrays with the forms
//
//     TYPE SNAMENorm(   TYPE supertensor[2][2][2][2]);
//     TYPE SNAMEMax(    TYPE * supertensor, int cubes, int slices, int rows, int cols);
//     TYPE SNAMEMin(    int cubes, int slices, int rows, int cols, TYPE * supertensor);
//     void SNAMEScale(  TYPE supertensor[3][3][3][3]);
//     void SNAMEFloor(  TYPE * array,  int cubes, int slices, int rows, int cols, TYPE floor);
//     void SNAMECeil(   int slices, int cubes, int slices, int rows, int cols, TYPE * array, TYPE ceil);
//     void SNAMELUSplit(TYPE in[2][2][2][2], TYPE lower[2][2][2][2], TYPE upper[2][2][2][2]);
//
// for any specified type TYPE (for example: short, unsigned int, long
// long, etc.) with given short name SNAME (for example: short, uint,
// longLong, etc.).  The macro is then expanded for the given
// TYPE/SNAME pairs.  The resulting functions are for testing numpy
// interfaces, respectively, for:
//
//  * 4D input arrays, hard-coded length
//  * 4D input arrays
//  * 4D input arrays, data last
//  * 4D in-place arrays, hard-coded lengths
//  * 4D in-place arrays
//  * 4D in-place arrays, data last
//  * 4D argout arrays, hard-coded length
//
#define TEST_FUNCS(TYPE, SNAME) \
\
TYPE SNAME ## Norm(TYPE supertensor[2][2][2][2]) {	     \
  double result = 0;				     \
  for (int l=0; l<2; ++l)			     \
    for (int k=0; k<2; ++k)			     \
      for (int j=0; j<2; ++j)			     \
        for (int i=0; i<2; ++i)			     \
          result += supertensor[l][k][j][i] * supertensor[l][k][j][i]; \
  return (TYPE)sqrt(result/16);			     \
}						     \
\
TYPE SNAME ## Max(TYPE * supertensor, int cubes, int slices, int rows, int cols) { \
  int i, j, k, l, index;						\
  TYPE result = supertensor[0];					\
  for (l=0; l<cubes; ++l) {					\
    for (k=0; k<slices; ++k) {					\
      for (j=0; j<rows; ++j) {					\
        for (i=0; i<cols; ++i) {					\
          index = l*slices*rows*cols + k*rows*cols + j*cols + i;			\
          if (supertensor[index] > result) result = supertensor[index];	\
        }                           \
      }								\
    }								\
  }								\
  return result;						\
}								\
\
TYPE SNAME ## Min(int cubes, int slices, int rows, int cols, TYPE * supertensor) {	\
  int i, j, k, l, index;						\
  TYPE result = supertensor[0];					\
  for (l=0; l<cubes; ++l) {					\
    for (k=0; k<slices; ++k) {					\
      for (j=0; j<rows; ++j) {					\
        for (i=0; i<cols; ++i) {					\
          index = l*slices*rows*cols + k*rows*cols + j*cols + i;			\
          if (supertensor[index] < result) result = supertensor[index];	\
        }                       \
      }								\
    }								\
  }								\
  return result;						\
}								\
\
void SNAME ## Scale(TYPE array[3][3][3][3], TYPE val) { \
  for (int l=0; l<3; ++l)			     \
    for (int k=0; k<3; ++k)			     \
      for (int j=0; j<3; ++j)			     \
        for (int i=0; i<3; ++i)			     \
          array[l][k][j][i] *= val;			     \
}						     \
\
void SNAME ## Floor(TYPE * array, int cubes, int slices, int rows, int cols, TYPE floor) { \
  int i, j, k, l, index;							     \
  for (l=0; l<cubes; ++l) {			     \
    for (k=0; k<slices; ++k) {						     \
      for (j=0; j<rows; ++j) {						     \
        for (i=0; i<cols; ++i) {						     \
          index = l*slices*rows*cols + k*rows*cols + j*cols + i;			\
          if (array[index] < floor) array[index] = floor;			     \
        }									     \
      }                                         \
    }									     \
  }									     \
}									     \
\
void SNAME ## Ceil(int cubes, int slices, int rows, int cols, TYPE * array, TYPE ceil) { \
  int i, j, k, l, index;							   \
  for (l=0; l<cubes; ++l) {						   \
    for (k=0; k<slices; ++k) {						   \
      for (j=0; j<rows; ++j) {						   \
        for (i=0; i<cols; ++i) {						   \
          index = l*slices*rows*cols + k*rows*cols + j*cols + i;			\
          if (array[index] > ceil) array[index] = ceil;			   \
        }                                   \
      }									   \
    }									   \
  }									   \
}									   \
\
void SNAME ## LUSplit(TYPE supertensor[2][2][2][2], TYPE lower[2][2][2][2], \
		      TYPE upper[2][2][2][2]) {			 \
  int sum;							 \
  for (int l=0; l<2; ++l) {					 \
    for (int k=0; k<2; ++k) {					 \
      for (int j=0; j<2; ++j) {					 \
        for (int i=0; i<2; ++i) {					 \
          sum = i + j + k + l;					 \
          if (sum < 2) {						 \
            lower[l][k][j][i] = supertensor[l][k][j][i];			 \
            upper[l][k][j][i] = 0;					 \
          } else {						 \
            upper[l][k][j][i] = supertensor[l][k][j][i];			 \
            lower[l][k][j][i] = 0;					 \
          }							 \
        }                           \
      }								 \
    }								 \
  }								 \
}

TEST_FUNCS(signed char       , schar    )
TEST_FUNCS(unsigned char     , uchar    )
TEST_FUNCS(short             , short    )
TEST_FUNCS(unsigned short    , ushort   )
TEST_FUNCS(int               , int      )
TEST_FUNCS(unsigned int      , uint     )
TEST_FUNCS(long              , long     )
TEST_FUNCS(unsigned long     , ulong    )
TEST_FUNCS(long long         , longLong )
TEST_FUNCS(unsigned long long, ulongLong)
TEST_FUNCS(float             , float    )
TEST_FUNCS(double            , double   )