set library 1;

complex matrix diag(complex vector v) {
  complex matrix m[#v];
  int i;
  for i=0 to #v-1 { m[i,i]=v[i]; }
  return m;
}

complex vector unitv(int d,int i) {
  complex vector v[d];
  v[i]=1;
  return v;
}

complex matrix unitm(int d) {
  complex matrix m[d];
  int i;
  for i=0 to d-1 { m[i,i]=1; }
  return m;
}

complex matrix trans(complex matrix m) {
  complex matrix mm[#m];
  int i; int j;
  for i=0 to #m-1 {
    for j=0 to #m-1 {
      mm[i,j]=m[j,i];
    }
  }
  return mm;
}

complex matrix adj(complex matrix m) {
  complex matrix mm[#m];
  int i; int j;
  for i=0 to #m-1 {
    for j=0 to #m-1 {
      mm[i,j]=conj(m[j,i]);
    }
  }
  return mm;
}

complex vector column(complex matrix m,int j) {
  int i;
  complex vector v[#m];
  for i=0 to #m-1 {
    v[i]=m[i,j];
  }
  return v;
}

complex vector row(complex matrix m,int i) {
  int j;
  complex vector v[#m];
  for j=0 to #m-1 {
    v[j]=m[i,j];
  }
  return v;
}

complex matrix submatrix(complex matrix m,int i,int j) {
  complex matrix mm[#m-1];
  int ii; int jj;
  for ii=0 to #m-2 {
    for jj=0 to #m-2 {
      mm[ii,jj]=m[ii+int(ii>=i),jj+int(jj>=j)];
    }
  }
  return mm;
}

complex det(complex matrix m) {
  int i;
  complex d;
  if #m==1 { return m[0,0]; }
  for i=0 to #m-1 {
    if m[i,0]!=0 { 
      d=d+(-1)^i*m[i,0]*det(submatrix(m,i,0)); 
    }
  }
  return d;
}

complex trace(complex matrix m) {
  int i;
  complex t;
  for i=0 to #m-1 { t=t+m[i,i]; }
  return t;
}

complex matrix inverse(complex matrix m) {
  int i; int j;
  complex matrix mm[#m];
  complex d=det(m);
  if d==0 { exit "matrix is singular"; }
  for i=0 to #m-1 {
    for j=0 to #m-1 {
      mm[i,j]=(-1)^(i+j)*det(submatrix(m,j,i))/d;
    }
  }
  return mm;
}

const epsilon=tensor3(0,0,0,0,0,1,0,-1,0,0,0,-1,0,0,0,1,0,0,0,1,0,-1,0,0,0,0,0);

complex vector cross(complex vector a,complex vector b) {
  return (epsilon*b)*a;
}

complex matrix outer(complex vector a,complex vector b) {
  int i; int j; 
  complex matrix m[#a];
  for i=0 to #a {
    for j=0 to #b {
      m[i,j]=a[i]*b[j];
    }
  }
  return m;
}
      
complex tensor3 outervm(complex vector v,complex matrix m) {
  int i; int j; int k;
  complex tensor3 t[#v];
  for i=0 to #v-1 {
    for j=0 to #m-1 {
      for k=0 to #m-1 {
        t[i,j,k]=v[i]*m[j,k];
      }
    }
  }
  return t;
}

complex tensor3 outermv(complex matrix m,complex vector v) {
  int i; int j; int k;
  complex tensor3 t[#m];
  for i=0 to #m-1 {
    for j=0 to #m-1 {
      for k=0 to #v-1 {
        t[i,j,k]=m[i,j]*v[k];
      }
    }
  }
  return t;
}

set library 0;