QUADRIQUES ET REDUCTION DE GAUSS

//-------------------------------------------------
// GAUSS
//-------------------------------------------------
void gauss(double **M,int N)
{
int n,x,y,k;

n = N-1;
while(n > 0)
  {

  // si on n'a pas un terme carre ...
  if(CmpDouble0(M[n][n]))
    {
    int s; // indice du "swap"

    // ... on cherche un autre terme carre ...
    s = -1;
    for(k=0;k<n;k++)
      {
      if(!CmpDouble0(M[k][k]))
        {
        s = k;
        break;
        }
      }
    
    if(-1 == s)
      {
      // ... si on n'en trouver decidement pas, on regarde si on a des produits
      // croise entre la derniere variable et l'avant derniere, puis on elimine
      // les deux dernieres lignes et colonne :
      if(!CmpDouble0(M[n][n-1]))
        {
        // 
        // ( *   *   *   0   0    )
        // ( *   *   *   0   0    )
        // ( *   *   *   0   0    )
        // ( 0   0   0  1/2  0    )
        // ( 0   0   0   0  -1/2  )
        // 
        for(y=0;y<n-1;y++)
          {
          for(x=0;x<n-1;x++)
            {
            M[y][x] -= (M[n-1][x]*M[n][y]+M[n-1][y]*M[n][x])/M[n][n-1];
            }
          }
        M[n-1][n-1] = + 0.5;
        M[n-1][n  ] =   0.0;
        M[n  ][n-1] =   0.0;
        M[n  ][n  ] = - 0.5;

        // on diminue de 2 le rang
        n -= 2;
        }
      else
        {
        // sinon
        // 
        // ( 0   *   *   *   0  )
        // ( *   0   *   *   0  )
        // ( *   *   0   *   0  )
        // ( *   *   *   0   0  )
        // ( 0   0   0   0   0  )
        // 

        // on diminue de 1 le rang
        n -= 1;
        }
      }
    else
      {
      // ... si on en a trouve un, on l'echange et on recommence la boucle au meme rang
      // mais cette fois-ci avec un terme carre
      for(k=0;k<=n;k++) // ligne : L[n]  <->  L[s]
        {
        SwapDouble(&M[n][k],&M[s][k]);
        }
      for(k=0;k<=n;k++) // colonne : C[n]  <->  C[s]
        {
        SwapDouble(&M[k][n],&M[k][s]);
        }
      }
    }
  else
    {
    // ... on a donc un terme carre, on elimine la derniere ligne et colonne :
    // 
    // ( * * * * 0 )
    // ( * * * * 0 )
    // ( * * * * 0 )
    // ( * * * * 0 )
    // ( 0 0 0 0 ? )
    // 
    for(y=0;y<n;y++)
      {
      for(x=0;x<n;x++)
        {
        M[y][x] -= M[n][x] * M[y][n] / M[n][n];
        }
      }

    for(k=0;k<n;k++)
      {
      M[k][n] = M[n][k] = 0.;
      }

    // on passe au rang inferieur
    n --;
    }
  }
} // gauss()


//-------------------------------------------------
void sign(int *p,int *q,double **mat,int n)
{
double **diag;
int      k;

(*p) = 0;
(*q) = 0;

diag = CopyMat(mat,n);
gauss(diag,n);

for(k=0;k<n;k++)
  {
  if(diag[k][k] > MATH_EPSILON)
    {
    (*p) ++;
    }
  else if(diag[k][k] < -MATH_EPSILON)
    {
    (*q) ++;
    }
  }

DeleteMat(diag,n);
} // sign()


//-------------------------------------------------
void info(int pA,int qA,int pB,int qB)
{
if(4==pB && 0==qB && 3==pA && 0==qA)
  {
  printf("ensemble vide (quadrique propre)\n");
  }
else if(3==pB && 1==qB && 3==pA && 0==qA)
  {
  printf("ellipsoide (quadrique propre)\n");
  }
else if(3==pB && 1==qB && 2==pA && 1==qA)
  {
  printf("hyperboloide à deux nappes (quadrique propre)\n");
  }
else if(3==pB && 1==qB && 2==pA && 0==qA)
  {
  printf("paraboloide elliptique (quadrique propre)\n");
  }
else if(2==pB && 2==qB && 2==pA && 1==qA)
  {
  printf("hyperboloide à une nappe (quadrique propre)\n");
  }
else if(2==pB && 2==qB && 1==pA && 1==qA)
  {
  printf("paraboloide hyperbolique (quadrique propre)\n");
  }
else if(3==pB && 0==qB && 3==pA && 0==qA)
  {
  printf("point\n");
  }
else if(3==pB && 0==qB && 2==pA && 0==qA)
  {
  printf("ensemble vide\n");
  }
else if(2==pB && 1==qB && 2==pA && 1==qA)
  {
  printf("cone\n");
  }
else if(2==pB && 1==qB && 2==pA && 0==qA)
  {
  printf("cylindre elliptique\n");
  }
else if(2==pB && 1==qB && 1==pA && 1==qA)
  {
  printf("cylindre hyperbolique\n");
  }
else if(2==pB && 1==qB && 1==pA && 0==qA)
  {
  printf("cylindre parabolique\n");
  }
else if(2==pB && 0==qB && 2==pA && 0==qA)
  {
  printf("droite\n");
  }
else if(2==pB && 0==qB && 1==pA && 0==qA)
  {
  printf("ensemble vide\n");
  }
else if(1==pB && 1==qB && 1==pA && 1==qA)
  {
  printf("deux plans sécants\n");
  }
else if(1==pB && 1==qB && 1==pA && 0==qA)
  {
  printf("deux plans parallèles\n");
  }
else if(1==pB && 0==qB && 1==pA && 0==qA)
  {
  printf("plan double\n");
  }
else
  {
  printf("impossible\n");
  }
} // info()