// cross_product6d.c  vector cross product for 6 dimensions   
//          from 2 to 5 vectors input
// 
// Cross product F = A x B x C x D x E 
// F orthogonal to A, B, C, D, and E
//
// 6 dimensions x,y,z,p,q,r unit vectors,  F = |f1, f2, f3, f4, f5, f6|
//
//        | x  y  z  p  q  r|  f1 =  det|not row 1, not col 1|
//        |a1 a2 a3 a4 a5 a6|  f2 = -det|not row 1, not col 2|
// F = det|b1 b2 b3 b4 b5 b6|  f3 =  det|not row 1, not col 3|
//        |c1 c2 c3 c4 c5 c6|  f4 = -det|not row 1, not col 4|
//        |d1 d2 d3 d4 d5 d6|  f5 =  det|not row 1, not col 5|
//        |e1 e2 e3 e4 e5 e6|  f6 =  det|not row 1, not col 6|
//
//          |a2 a3 a4 a5 a6|           |a1 a3 a4 a5 a6|
//          |b2 b3 b4 b5 b6|           |b1 b3 b4 b5 b6|
//  f1 = det|c2 c3 c4 c5 c6|  f2 = -det|c1 c3 c4 c5 c6|  ...
//          |d2 d3 d4 d5 d6|           |d1 d3 d4 d5 d6|
//          |e2 e3 e4 e5 e5|           |e1 e3 e4 e5 e6|
//
//          |a1 a2 a3 a4 a6|           |a1 a2 a3 a4 a5|
//          |b1 b2 b3 b4 b6|           |b1 b2 b3 b4 b5|
//  f5 = det|c1 c2 c3 c4 c6|  f6 = -det|c1 c2 c3 c4 c5|
//          |d1 d2 d3 d4 d6|           |d1 d2 d3 d4 d5|
//          |e1 e2 e3 e4 e6|           |e1 e2 e3 e4 e5|
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "determinant.h"

static void cross6(double M6[], double F[])
{
  double M5[5][5]; // (n-1)*(n-1) matrix for determinant
  int n=6;
  int i, j, k, ii, jj;

  for(k=0; k<n; k++)
  {
    for(i=1; i<n; i++) // skip top row  all 1.0
    {
      jj = 0;
      for(j=0; j<n; j++)
      {
        if(j==k) continue; // skip row k
        jj = jj +1;
        // printf("k=%1d, i=%1d, j=%1d, ii=%1d, jj=%1d, ix=%d \n",
	//       k, i, j, i-1, jj-1, i*n+j);
        M5[i-1][jj-1] = M6[i*n+j];
      } // j
    } // i
    // for(i=0; i<n-1; i++) printf("M5 |%f  %f  %f  %f  %f| \n",
    //		            M5[i][0], M5[i][1], M5[i][2], M5[i][3], M5[i][4]); 
    F[k] = determinant(n-1, M5);
    if(k%2==1) F[k] = - F[k];
    // printf("k=%d, det=%e \n", k, F[k]);
  } // k
} // end cross6

static void cross6_2(double A[], double B[], double F[]) // for two vectors
{
  int i, n = 6;
  double M6[36];
  for(i=0; i<n; i++)
  {
    M6[i] = 1.0; // must build complete matrix
    M6[  n+i] = 1.0;
    if(i==0) M6[  n+i] = 0.0; // fill with linearly indpendent vectors
    M6[2*n+i] = 1.0;
    if(i==1) M6[2*n+i] = 0.0;
    M6[3*n+i] = 1.0;
    if(i==2) M6[3*n+i] = 0.0;
    M6[4*n+i] = A[i]; // only two given
    M6[5*n+i] = B[i];
  }
  cross6(M6, F); // compute orthogonal vector
} // end cross6_2

static void cross6_3(double A[], double B[], double C[], 
                     double F[]) // for three vectors
{
  int i, n = 6;
  double M6[36];
  for(i=0; i<n; i++)
  {
    M6[i] = 1.0; // must build complete matrix
    M6[  n+i] = 1.0;
    if(i==0) M6[  n+i] = 0.0; // fill with linearly independent vectors
    M6[2*n+i] = 1.0;
    if(i==1) M6[2*n+i] = 0.0;
    M6[3*n+i] = A[i]; // only three given
    M6[4*n+i] = B[i];
    M6[5*n+i] = C[i];
  }
  cross6(M6, F); // compute orthogonal vector
} // end cross6_3

static void cross6_4(double A[], double B[], double C[], 
                     double D[], double F[]) // for four vectors
{
  int i, n = 6;
  double M6[36];
  for(i=0; i<n; i++)
  {
    M6[i] = 1.0; // must build complete matrix
    M6[  n+i] = 1.0;
    if(i==0) M6[  n+i] = 0.0; // fill with linearly independent vectors
    M6[2*n+i] = A[i]; // only four given
    M6[3*n+i] = B[i];
    M6[4*n+i] = C[i];
    M6[5*n+i] = D[i];
  }
  cross6(M6, F); // compute orthogonal vector
} // end cross6_4

static void cross6_5(double A[], double B[], double C[],
                     double D[], double E[], double F[]) // for five vectors
{
  int i, n = 6;
  double M6[36];
  for(i=0; i<n; i++)
  {
    M6[i] = 1.0; // must build complete matrix
    M6[  n+i] = A[i];
    M6[2*n+i] = B[i];
    M6[3*n+i] = C[i];
    M6[4*n+i] = D[i];
    M6[5*n+i] = E[i];
  }
  cross6(M6, F); // compute orthogonal vector
} // end cross6_5

static double dotn(int n, double A[], double B[])
{
  double sum = 0.0;
  int i;
  for(i=0; i<n; i++) sum = sum + A[i]*B[i];
  return sum;
} // end dotn

int main()
{
  double A[] = {2.0, 3.0, 5.0, 8.0, 12.0, 1.0}; // random test input vectors
  double B[] = {1.0, 4.0, 6.0, 9.0, 13.0, 2.0};
  double C[] = {3.0, 5.0, 14.0, 8.0, 16.0, 1.0};
  double D[] = {6.0, 9.0, 8.0, 15.0, 19.0, 2.0};
  double E[] = {5.0, 6.0, 7.0, 12.0, 11.0, 1.0};
  double F[6]; // computed orthogonal vector

  printf("cross_product6d.c  six dimensional vectors in R6 \n");
  printf("vectors are orthogonal if dot product is zero \n");
  printf("cos(angle between A and B) = (A dot B)/(|A|*|B|) \n");

  printf("A= |%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         A[0], A[1], A[2], A[3], A[4], A[5]);
  printf("B= |%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         B[0], B[1], B[2], B[3], B[4], B[5]);
  printf("C= |%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         C[0], C[1], C[2], C[3], C[4], C[5]);
  printf("D= |%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         D[0], D[1], D[2], D[3], D[4], D[5]);
  printf("E= |%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         E[0], E[1], E[2], E[3], E[4], E[5]);
  printf(" \n");

  printf("compute F vector orthogonal to A, B, C, D, E \n");
  cross6_5(A, B, C, D, E, F);
  printf("F=A x B x C x D x E=|%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         F[0], F[1], F[2], F[3], F[4], F[5]);
  printf("check dotn(6, A, F)=%8.2f \n", dotn(6, A, F));
  printf("check dotn(6, B, F)=%8.2f \n", dotn(6, B, F));
  printf("check dotn(6, C, F)=%8.2f \n", dotn(6, C, F));
  printf("check dotn(6, D, F)=%8.2f \n", dotn(6, D, F));
  printf("check dotn(6, E, F)=%8.2f \n", dotn(6, E, F));
  printf(" \n");

  printf("compute F vector orthogonal to A, B, C, D, not unique \n");
  cross6_4(A, B, C, D, F);
  printf("F=A x B x C x D =|%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         F[0], F[1], F[2], F[3], F[4], F[5]);
  printf("check dotn(6, A, F)=%8.2f \n", dotn(6, A, F));
  printf("check dotn(6, B, F)=%8.2f \n", dotn(6, B, F));
  printf("check dotn(6, C, F)=%8.2f \n", dotn(6, C, F));
  printf("check dotn(6, D, F)=%8.2f \n", dotn(6, D, F));
  printf(" \n");

  printf("compute F vector orthogonal to A, B, C, not unique \n");
  cross6_3(A, B, C, F);
  printf("F=A x B x C =|%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         F[0], F[1], F[2], F[3], F[4], F[5]);
  printf("check dotn(6, A, F)=%8.2f \n", dotn(6, A, F));
  printf("check dotn(6, B, F)=%8.2f \n", dotn(6, B, F));
  printf("check dotn(6, C, F)=%8.2f \n", dotn(6, C, F));
  printf(" \n");

  printf("compute F vector orthogonal to A, B, not unique \n");
  cross6_2(A, B, F);
  printf("F=A x B =|%8.2f, %8.2f, %8.2f, %8.2f, %8.2f, %8.2f| \n",
         F[0], F[1], F[2], F[3], F[4], F[5]);
  printf("check dotn(6, A, F)=%8.2f \n", dotn(6, A, F));
  printf("check dotn(6, B, F)=%8.2f \n", dotn(6, B, F));
  printf(" \n");


  printf("cross_producr6d.c finished \n");
  return 0;
} // end cross_product6d.c