/*******************************************************************
 * Julius C demo case: Fractal Matrix Multiply
 * -- author Paolo D'Alberto 
 *
 *
 * This algorithm implements a version of matrix multiply:  
 * 1) it is for demostration only
 * 2) the algorithm implements a fractal algorithm with no standard layout 
 * 3) the algorithm has been simplified and therefore produces incorrect results
 * 4) the complexity is unaffected
 */


#define MAX 1500
#define fl(a) ((a>>1))
#define cl(a) ((a>>1)+((a&1)))


#define M0(r,c) 0
#define M1(r,c) (cl(r)*cl(c))
#define M2(r,c) ((cl(r)*c))
#define M3(r,c) (r*c-fl(r)*fl(c))

#define ADD000(c,a,b)  c         , a          ,  b
#define ADD012(c,a,b)  c         , a+M1(m,n)  ,  b+M2(n,p)
#define ADD101(c,a,b)  c+M1(m,p) , a          ,  b+M1(n,p)
#define ADD113(c,a,b)  c+M1(m,p) , a+M1(m,n)  ,  b+M3(n,p)
#define ADD220(c,a,b)  c+M2(m,p) , a+M2(m,n)  ,  b                    
#define ADD232(c,a,b)  c+M2(m,p) , a+M3(m,n)  ,  b+M2(n,p)
#define ADD321(c,a,b)  c+M3(m,p) , a+M2(m,n)  ,  b+M1(n,p)
#define ADD333(c,a,b)  c+M3(m,p) , a+M3(m,n)  ,  b+M3(n,p)

#define TYPE000(m,n,p)   cl(m), cl(n), cl(p)
#define TYPE012(m,n,p)   cl(m), fl(n), cl(p)
#define TYPE101(m,n,p)   cl(m), cl(n), fl(p)
#define TYPE113(m,n,p)   cl(m), fl(n), fl(p)
#define TYPE220(m,n,p)   fl(m), cl(n), cl(p)
#define TYPE232(m,n,p)   fl(m), fl(n), cl(p)
#define TYPE333(m,n,p)   fl(m), fl(n), fl(p)
#define TYPE321(m,n,p)   fl(m), cl(n), fl(p)

int C[MAX*MAX],B[MAX*MAX], A[MAX*MAX];

void mm(int *c, int *a, int *b,
	int m, int n, int p) {
  
  int i,j,k;

  for (i=0;i<m;i++) 
    for (j=0;j<p;j++)
      for (k=0;k<n; k++) 
	c[i*p+j] += a[i*n+k]*b[k*p+j];
  
}

#define LS 4

void fractal_matrix_multiply(int *c, int *a, int *b, int m, int n, int p) { 

  if (m <=LS &&  n<=LS && p<=LS)    
    mm(c,a,b,m,n,p);
  
  else {

    fractal_matrix_multiply(ADD000(c,a,b),TYPE000(m,n,p));  /* A0 += A0A0 */
    fractal_matrix_multiply(ADD101(c,a,b),TYPE101(m,n,p));  /* A1 += A0A1 */
    fractal_matrix_multiply(ADD220(c,a,b),TYPE220(m,n,p));  /* A2 += A2A0 */
    fractal_matrix_multiply(ADD321(c,a,b),TYPE321(m,n,p));  /* A3 += A2A1 */
    fractal_matrix_multiply(ADD333(c,a,b),TYPE333(m,n,p));  /* A3 += A3A3 */
    fractal_matrix_multiply(ADD113(c,a,b),TYPE113(m,n,p));  /* A1 += A1A3 */
    fractal_matrix_multiply(ADD232(c,a,b),TYPE232(m,n,p));  /* A2 += A3A2 */
    fractal_matrix_multiply(ADD012(c,a,b),TYPE012(m,n,p));  /* A0 += A1A2 */
  }
}

int main() { 

  int i;
  
  for (i=117; i<1000; i+=53)
    fractal_matrix_multiply(C,A,B,i,i,i);
  
  return 0;
}