// student4D.java   compile    javac -cp . student4D.java
//                  execute    java  -cp . student4D

public class student4D
{
  double A4[] = {1.0, 1.0, 1.0, 1.0};
  double B4[] = {1.0, 1.0, 3.0, -5.0};
  double C4[] = {1.0, -1.0, 0.0, 0.0};
  double D4[] = new double[4];
  double detM4;
  double dotp4;
  double M4[][] = {{1.0, 1.0, 1.0, 1.0},
      	           {1.0, 1.0, 1.0, 1.0},  // A  all 1.0
		   {1.0, 2.0, 4.0, 5.0},
		   {10.0, 21.0, 32.0, 43.0}};
  double M4u[][]= {{1.0, 0.0, 0.0, 0.0},
		   {0.0, 1.0, 0.0, 0.0},  // A  unit vectors
		   {0.0, 0.0, 1.0, 0.0},
		   {0.0, 0.0, 0.0, 1.0}};
  double err1, err2, err3, err4;
  double B4l;
    
  public student4D()
  {
    System.out.println("student4D.java running");

    prtvec(A4, "A4");
    prtvec(B4, "B4");
    prtvec(C4, "C4");
    dotp4 = dot(A4, B4);
    System.out.println("dot(A4, B4)="+dotp4);
    System.out.println(" ");
    dotp4 = dot(A4, C4);
    System.out.println("dot(A4, C4)="+dotp4);
    System.out.println(" ");
    dotp4 = dot(B4, C4);
    System.out.println("dot(B4, C4)="+dotp4);
    System.out.println(" ");

    cross(A4, B4, C4, D4);
    System.out.println("cross4(A4, B4, C4, D4)=");    
    prtvec(D4, "D4");
    System.out.println("check orthoganal");
    err1 = Math.abs(dot(A4,D4));
    System.out.println("dot(A4, D4)="+err1);
    err2 = Math.abs(dot(B4,D4));
    System.out.println("dot(B4, D4)="+err2);
    err3 = Math.abs(dot(C4,D4));
    System.out.println("dot(C4, D4)="+err3);
    System.out.println("cross(A4, B4, C4, D4) err="+(err1+err2+err3));    

    System.out.println("student4D.java finished");
  }

  // utility functions

  public double length(double A[])
  {             //  return length of vector A
    int n = A.length;
    double len = 0.0;
    for(int i=0; i<n; i++) { len = len + A[i]*A[i]; }
    return Math.sqrt(len);  
  } // end length
    
  public void normalize(double A[], double B[])
  {                    // vector A is direction of B with length 1.0
    int n = A.length;
    double len = length(A);
    for(int i=0; i<n; i++) { B[i] = A[i]/len; }
  } // end normalize

  public double dot(double A[], double B[])
  {
    return A[0]*B[0] + A[1]*B[1] + A[2]*B[2] + A[3]*B[3];
  } // end dot

  public void prtvec(double V[], String name)
  {
    int n = V.length;	
    for(int i=0; i<n; i++)
    {
      System.out.println(name+"["+i+"]="+V[i]);
    }
    System.out.println(" ");
  } // end prtvec

  public void prtmat(double A[][], String name)
  {
    int n = A.length;	
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        System.out.println(name+"["+i+"]["+j+"]="+A[i][j]);
      }
    }
    System.out.println(" ");
  } // end prtmat

  public void cross(double A[], double B[], double C[], double D[])
    {                       //  vector D orthogonal to A and B and C
    double A3[][] = new double[3][3]; // (n-1)*(n-1) matrix  n times
    int n=4; // 4 = A.length = B.length = C.length = D.length
    int i, j, k, ii, jj;
    int swap[] = { 1, 2, 3, 0 }; // D subscripts
    for(k=0; k<n; k++) // D subscripts
    {
      // System.out.println("k="+k);
      for(i=0; i<n-1; i++) // i=0 A, i=1 B, i=2 C
      {
	jj = 0;
	for(j=0; j<n; j++)
        {
          if(j==k) continue; // skip row k
	  jj = jj +1;
          if(i==0) A3[i][jj-1] = A[j];
          if(i==1) A3[i][jj-1] = B[j];
          if(i==2) A3[i][jj-1] = C[j];
	  // if(i==0) { System.out.println("j="+(j));}
        } // j
      } // i
      // matprt(A3);
      D[k] = det(A3);
      if(k%2==1) {D[k] = - D[k];}
    } // k
  } // end cross

  double det(final double A[][])
  {
    int n=A.length;
    double D = 1.0;                 // determinant
    double B[][]=new double[n][n];  // working matrix
    int row[]=new int[n];           // row interchange indicies
    int hold , I_pivot;             // pivot indicies
    double pivot;                   // pivot element value
    double abs_pivot;

    if(A[0].length!=n)
    {
      System.out.println("Error in det, inconsistent array sizes.");
    }
    // build working matrix
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
        B[i][j]=A[i][j];
    // set up row interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k]= k;
    }
    // begin main reduction loop
    for(int k=0; k<n-1; k++)
    {
      // find largest element for pivot
      pivot = B[row[k]][k];
      abs_pivot = Math.abs(pivot);
      I_pivot = k;
      for(int i=k; i<n; i++)
      {
        if( Math.abs(B[row[i]][k]) > abs_pivot )
        {
          I_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = Math.abs(pivot);
        }
      }
      // have pivot, interchange row indicies
      if(I_pivot != k)
      {
        hold = row[k];
        row[k] = row[I_pivot];
        row[I_pivot] = hold;
        D = - D;
      }
      // check for near singular
      if(abs_pivot < 1.0E-20)
      {
        return abs_pivot;
      }
      else
      {
        D = D * pivot;
        // reduce about pivot
        for(int j=k+1; j<n; j++)
        {
          B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
        }
        //  inner reduction loop
        for(int i=0; i<n; i++)
        {
          if(i != k)
          {
            for(int j=k+1; j<n; j++)
            {
              B[row[i]][j] = B[row[i]][j] - B[row[i]][k]* B[row[k]][j];
            }
          }
        }
      }
      //  finished inner reduction
    }
    // end of main reduction loop
    return D * B[row[n-1]][n-1];
  } // end det

    
  public static void main (String[] args)
  {
    new student4D();
  }
} // end student4D.java