#include<stdio.h>
#include<math.h>
#include <stdlib.h>
#define PI 3.14159265
 
double ProducteEscalar(double u[], double v[], int dim);
 
double Norma(double u[], int dim);
 
int BaseR3(double u[], double v[], double w[], int dim);
 
double Angle(double u[], double v[], int dim);
 
double Determinant(double u[], double v[], double w[]);
 
int Perpendiculars(double u[], double v[], int dim);
 
double ModAngle(double u[], double v1[], double v2[], int dim);
 
void ProjeccioOrtogonal(double u[], double v[],double w[], int dim);
 
void ProducteVectorial(double u[], double v[],double w[]);
 
int tipus2d(int n, double v[n][3], double u[], unsigned int vmin[], unsigned int * pmin);
 
int projectacon(int n, double u[], double v[][3], double proj[][3]);
 
int tipus3d(int n, double v[n][3],unsigned int vmin[], unsigned int * pmin);
 
 
int main()
{
unsigned int * p = malloc(sizeof(*p));
double d,e,f;
int n=0,i=0;
FILE * fitxer;
fitxer = fopen("poliedric4.dat","r");
if ( fitxer == NULL )
{
fprintf (stderr , " ERROR : El fitxer ’recta.txt’ no es pot obrir ... \n\n");
free(p);
return -100;
}
while (EOF!= fscanf (fitxer , "%lf %lf %lf ", &d, &e, &f)) //No llegeix fitxers amb comes.
{
 n++;
}
fclose(fitxer);
double v[n][3];
unsigned int vmin[n];
fopen("poliedric4.dat","r");
while (EOF!= fscanf (fitxer , "%lf %lf %lf ", &d, &e, &f))
{
 v[i][0]=d;
 v[i][1]=e;
 v[i][2]=f;
 i++;
}
fclose(fitxer);
 
int a = tipus3d(n, v, vmin, &(*p));
 
FILE * meufitxer ;
meufitxer = fopen (" generadorsminimals.txt ", "w");
if ( meufitxer == NULL )
{
fprintf (stderr , " ERROR : El fitxer ’generadorsminimals.txt’ no es pot obrir ...\n\n");
free(p);
return 1;
}
 for(int k=0;k<(*p);k++)
 {
     fprintf(meufitxer, "%lf\t", v[vmin[k]][0]) ;
     fprintf(meufitxer, "%lf\t", v[vmin[k]][1]) ;
     fprintf(meufitxer, "%lf\t", v[vmin[k]][2]) ;
     fprintf(meufitxer, "\n") ;
 }
fclose(meufitxer);
 
  printf("S'ha llegit un con generat per %d vectors.", n);
  if(a==0)
  {
      printf("\nEl con és el {0}.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
    if(a==1)
  {
      printf("\nEl con és una semirecta.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
    if(a==2)
  {
      printf("\nEl con és una recta.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
    if(a==3)
  {
      printf("\nEl con és un angle del pla.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
    if(a==4)
  {
      printf("\nEl con és un semiplà.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
    if(a==5)
  {
      printf("\nEl con és un pla.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      return 100;
  }
  if(a==6)
  {
      printf("\nEl con és polièdric amb %d arestes.", *p);
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
      if(a==7)
  {
      printf("\nEl con és un angle dièdric.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
      if(a==8)
  {
      printf("\nEl con és un semiespai.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
      }
      free(p);
      return 100;
  }
  if(a==9)
  {
      printf("\nEl con és tot l'espai.");
      printf("\nUn sistema de generadors és:\n");
      for(int i=0;i<(*p);i++)
      {
          printf("\n");
          for(int j=0;j<3;j++)
          {
              printf("%lf\t", v[vmin[i]][j]);
          }
   }
   free(p);
   return 100;
  }
  free(p);
return -100;
}
int tipus3d(int n, double v[n][3],unsigned int vmin[], unsigned int * pmin)
{
    * pmin = 0;
    for (int i=0;i<100;i++) //Buidem el vector
    {
        vmin[100]=0;
    }
    unsigned int vcares[n];  //Llista de <= n-1 enters.
    for (int i=0;i<n-1;i++) //Buidem el vector
    {
        vcares[i]=0;
    }
    unsigned int wmin[4];  //Llista de <= 4 enters.
    for (int i=0;i<4;i++) //Buidem el vector
    {
        wmin[i]=0;
    }
    int i, j, m, tip;
    int repetit;
    double w[n][3];   //Llista de m <= n vectors.
    unsigned int * long_wmin  = malloc(sizeof(*long_wmin));
    unsigned int * long_vcares  = malloc(sizeof(*long_vcares));
    double perp[3];
 
    if(n == 0)  //Con zero.
      {
          free(long_wmin);
          free(long_vcares);
         return 0;
      }
    if(n == 1) //Semirecta.
      {
         vmin[(*pmin)]=0;
         (*pmin)++;
         free(long_wmin);
         free(long_vcares);
         return 1;
      }
    //A partir d'aquí  n >= 1
   m = projectacon(n, v[0], v, w);
   if(m == 0) // Projecció igual a punt: dimensió 1
   {
       for(i=0;i<n;i++) // Cal distingir si recta o semirecta
       {
           if(ProducteEscalar(v[0], v[i], 3)<0)
           {
               vmin[(*pmin)]=0; // Algun angle &#960;: recta
               (*pmin)++;
               vmin[(*pmin)]=i;
               (*pmin)++;
               free(long_wmin);
               free(long_vcares);
               return 2;
           }
       }
       vmin[(*pmin)]=0;
       (*pmin)++; // Cap angle &#960;: semirecta
        free(long_wmin);
        free(long_vcares);
       return 1;
   }
   tip = tipus2d(m, w, v[0], wmin, long_wmin);
   if(tip<3)
   {
       ProducteVectorial(v[0],w[0], perp); // Vector perpendicular
       free(long_wmin);
       free(long_vcares);
       return tipus2d(n, v, perp, vmin, &(*pmin));;
   }
   //A partir d'aquí el con és de dimensió 3.
   for(i=0;i<n;i++)
   {
       m = projectacon(n, v[i], v, w );
       tip = tipus2d(m, w, v[i], wmin, long_wmin);
       if(tip==3)       // projecció és angle; v[i] genera aresta
       {
           repetit = 0; // control d’aresta ja trobada
           for(j=0;j<(*pmin); j++)
           {
               if(Angle(v[i],v[vmin[j]], 3)==0)
               {
                   repetit = 1;
               }
           }
           if(repetit == 0)
           {
               vmin[(*pmin)]=i;
               (*pmin)++;
           }
       }
       if(tip == 4) // projecció és semiplà; v[i] està a cara
       {
           vcares[(*long_vcares)] = i;
           (*long_vcares)++;
       }
   }
   // En aquest punt vmin generen les arestes.
   if ((*pmin)>2)// con polièdric
   {
       free(long_wmin);
       free(long_vcares);
       return 6;
   }
   if ((*pmin)==2)// con dièdric; cal trobar dues cares
   {
       perp[0]=v[vmin[0]][0]; // vector aresta, per projectar
       perp[1]=v[vmin[0]][1];
       perp[2]=v[vmin[0]][2];
       vmin[(*pmin)]=vcares[0];  // primera aresta-cara
       (*pmin)++;
       ProjeccioOrtogonal(perp, v[vcares[0]], w[0], 3); // projecció aresta-cara
       for(i=1;i<(*long_vcares); i++)
       {
           ProjeccioOrtogonal(perp, v[vcares[i]], w[1], 3); // projecció
           if(Angle(w[0], w[1], 3)!=0) // és a l’altra cara
           {
               vmin[(*pmin)]=vcares[i];   // segona aresta-cara
               (*pmin)++;
               free(long_wmin);
               free(long_vcares);
               return 7;
           }
       }
   }
   if((*long_vcares)>0)// És un semiespai
   {
       for(i=0;i<m;i++)// Col.leccionem els vectors cara
       {
           w[i][0]=0;
           w[i][1]=0;
           w[i][2]=0;
       }
    j=n-1;             // i també un vector interior
 
    for(i=0;i<(*long_vcares);i++)
    {
        w[i][0]=v[vcares[i]][0];
        w[i][1]=v[vcares[i]][1];
        w[i][2]=v[vcares[i]][2];
 
        if((abs(Angle(w[0], w[i], 3)))>0 && (abs(Angle(w[0], w[i], 3))) < PI)
        {
            ProducteVectorial(w[0], w[i], perp);// Un vector ortogonal
        }
        if((j) == n-1 && vcares[i]>i)
        {
            i=j;
        }
    }
    tip=tipus2d(m, w, perp, wmin, long_wmin);//És la mateixa llargada per lin.152
    for(i=0;i<(*long_wmin);i++)
    {
        vmin[(*pmin)]=vcares[wmin[i]];
        (*pmin)++;
    }
    vmin[(*pmin)]=j;
    (*pmin)++;
    free(long_wmin);
    free(long_vcares);
    return 8;
   }
   for(i=0;i<m;i++)// Col.leccionem els vectors cara
       {
           w[i][0]=0;
           w[i][1]=0;
           w[i][2]=0;
       }
   for(i=1;i<n;i++)// Tots menys v[0]
   {
       w[i][0]=v[i][0];
       w[i][1]=v[i][1];
       w[i][2]=v[i][2];
   }
 i=0;
  while(tipus3d(n-1, w, vmin, pmin)<9)// No es pot treure v[i]
  {
      w[i][0]=v[i][0]; // Posar-lo, i ara no incloure v[i+1]
      w[i][1]=v[i][1];
      w[i][2]=v[i][2];
      i++;
      if(i==n)// No se’n pot treure cap!
         {
             vmin[(*pmin)]=n-1;
             (*pmin)++;
             free(long_wmin);
             free(long_vcares);
             return 9;
         }
         for (int i=0;i<4;i++) // Preparat per cridar tipus3d
    {
        vmin[i]=0;
    }
  }
for(j=1;j<(*pmin);j++)// Com que v[i] falta a w, corregir els índexs trobats recursivament
{
    vmin[j]=vmin[j]+1; // vmin s’ha calculat recursivament
}
free(long_wmin);
free(long_vcares);
return 9;
}
 
 
double ProducteEscalar(double u[], double v[], int dim)
{
	double PE = 0;
	for (int i=0; i<dim; i++)
       {
		PE += u[i]*v[i];
	   }
	return PE;
}
 
 
double Norma(double u[], int dim)
{
    double n = 0;
    for (int i=0; i<dim; i++)
      {
        n += u[i]*u[i];
      }
    return sqrt(n);
}
 
 
int BaseR3(double u[], double v[], double w[], int dim)
{
    //Si hi ha vectors linealment dependents retorna 0, si són linealment independents retorna 1.
  for(int i=0;i+1<dim;i++)
 {
 if(v[i]/u[i] == v[i+1]/u[i+1])
 {
  return 0;
 }
 if (v[i]/w[i] == v[i+1]/w[i+1])
 {
  return 0;
 }
  if(u[i]/w[i] == u[i+1]/w[i+1])
 {
  return 0;
 }
 }
    return 1;
}
 
 
double Angle(double u[], double v[], int dim)
{
    return acos(ProducteEscalar(u,v, dim)/(Norma(u, dim)*Norma(v, dim)));
}
 
 
double Determinant(double u[], double v[], double w[])
{
    // El calcularem per la regla de Sarrurs. Només funciona amb tres vectors de R3.
    return u[0]*v[1]*w[2] + u[1]*v[2]*w[0] + u[2]*v[0]*w[1] - u[2]*v[1]*w[0] - u[0]*v[2]*w[1] - u[1]*v[0]*w[2];;
}
 
 
int Perpendiculars(double u[], double v[], int dim)
{
   // Si són perpendiculars retorna 1, si no ho són retorna 0.
   // Com tots els vectors tenen origen en el 0, només cal comprovar que el producte escalar es 0, és a dir, que formen un angle recte.
   if(abs(ProducteEscalar(u, v, dim))<0.0003)//Revisar la tolerància!!
     {
      return 1;
     }
    else
    {
   return 0;
    }
}
 
 
double ModAngle(double u[], double v1[], double v2[], int dim)
{
    // Modificació de la funció angle. Calcula angle amb signe. Els vectors v_1 i v_2 han de ser perpendiculars a u.
    // No comprova si són perpendiculars.
    {
    int signe;
 
     if(Determinant(u, v1, v2)>=0) //Accepta el conveni: si el determinant és 0 té signe positiu.
     {
       signe = 1  ;
     }
     if(Determinant(u, v1, v2)<0)
     {
       signe = -1;
     }
    return (signe*acos(ProducteEscalar(v1,v2, dim)/(Norma(v1, dim)*Norma(v2, dim))));
    }
}
 
 
void ProjeccioOrtogonal(double u[], double v[],double w[], int dim)
{
    // Has de declarar amb anterioritat el vector w i la seva dimensió.
  for(int i=0;i<dim;i++)
{
    w[i]= v[i]-(ProducteEscalar(u,v,dim)/ProducteEscalar(u,u,dim))*u[i] ;
}
  return ;
}
 
 
void ProducteVectorial(double u[], double v[],double w[])
{
     // Has de declarar amb anterioritat el vector w i la seva dimensió.
     // Només funciona amb vectors a R3.
     w[0]= u[1]*v[2]-u[2]*v[1];
     w[1]= u[0]*v[2]-u[2]*v[0];
     w[2]= u[0]*v[1]-u[1]*v[0];
 
    return;
}
 
 
int tipus2d(int n, double v[n][3], double u[], unsigned int vmin[], unsigned int * pmin)
{
    * pmin = 0;
    for (int i=0;i<4;i++) //Buidem el vector
    {
        vmin[i]=0;
    }
    int i,j;
    int dr, esq;
    double angle[n][n];
    double anglemax;
 
      if(n == 0)  //Con zero.
        {
          return 0;
        }
      if(n == 1) //Semirecta.
        {
          vmin[(*pmin)]=0;
          (*pmin)++;
          return 1;
        }
 
    dr = 0;
    esq = 1;
    anglemax = ModAngle(u, v[dr], v[esq], 3);
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++) //Troba anglemax.
        {
            angle[i][j] = ModAngle(u, v[i], v[j], 3);
            if(angle[i][j]>anglemax)
            {
                dr = i;
                esq = j;
                anglemax = angle[dr][esq];
            }
        }
    }
    if(anglemax==0) //És semirecta.
    {
       vmin[(*pmin)] = 0;
       (*pmin)++;
       return 1;
    }
    if(anglemax < M_PI) //És un angle pla o un pla.
    {
        for(i=0;i<n;i++)
        {
            if(angle[dr][i]<0)
            {  vmin[(*pmin)]=dr;
               (*pmin)++;
               vmin[*pmin]=esq;
               (*pmin)++;
               vmin[(*pmin)]=i;
              (*pmin)++;
               return 5; //És un pla.
            }
 
        }
                //Si hem arribat aquí és un angle pla.
    vmin[(*pmin)]=dr;
    (*pmin)++;
    vmin[(*pmin)]=esq;
    (*pmin)++;
    return 3;
    }
    // Trobar si es una recta.
    for(i=0;i<n; i++)
    {
        if(angle[dr][i]<0)
        {
               vmin[(*pmin)]=dr;
               (*pmin)++;
               vmin[(*pmin)]=esq;
               (*pmin)++;
               vmin[(*pmin)]=i;
               (*pmin)++;
               return 5; //És un pla.
        }
    }
    for(i=0;i<n;i++)
    {
        if(abs(angle[dr][i])>0.000001 && abs(angle[dr][i]-M_PI)>0.000001)//És un semipla
        {
               vmin[(*pmin)]=dr;
               (*pmin)++;
               vmin[(*pmin)]=esq;
               (*pmin)++;
               vmin[(*pmin)]=i;
               (*pmin)++;
 
               return 4;
        }
    }
    //Si hem arribat aquí és una recta.
    vmin[(*pmin)]=dr;
    (*pmin)++;
    vmin[(*pmin)]=esq;
    (*pmin)++;
    return 2;
}
int projectacon(int n, double u[], double v[n][3], double proj[][3])
{
   int i,j=0;
   double w[n][3];
   int b = n; //Nombre d'enters amb projecció no nul·la
   for(i=0;i<n;i++)
     {
         ProjeccioOrtogonal(u,v[i],w[i],3);
     }
    for(i=0;i<n;i++)
     {
         if(v[i][0]/u[0]==v[i][1]/u[1] && v[i][0]/u[0]==v[i][2]/u[2] && v[i][1]/u[1]==v[i][2]/u[2])
         {
             w[i][0]=0;
             w[i][1]=0;
             w[i][2]=0;
         }
     }
     for(i=0;i<n;i++)
     {
         if(w[i][0]==0 && w[i][1]==0 && w[i][2]==0)
         {
             b -= 1;
         }
     }
     for(i=0;i<n;i++)
     {
         if(w[i][0]==0 && w[i][1]==0 && w[i][2]==0)
         {
             j++;
         }
         else
         {
             proj[i-j][0] = w[0][i];
             proj[i-j][1] = w[1][i];
             proj[i-j][2] = w[2][i];
         }
     }
return b;
}