/***************************************
*                                      *
*   Copyright (c) 1998 Jean-Eric Pin   *
*   All rights reserved.               *
*                                      *
*   TAB = 2 spaces                     *
*                                      *
***************************************/

/*-------------------------------------------------------------------
 * Dclasses.c   Jean-Eric Pin 07/12/96
 *-------------------------------------------------------------------
 * Probleme majeur : les procedures recursives provoquent des erreurs
 * 28 : "Stack has moved into application heap"
 * La solution sur Mac 68k est SetApplLimit(GetApplLimit() - 16384)
 * pour augmenter la pile de 16k. Cependant, pour des graphes tres
 * grands, la pile augmente trop vite (80 octets par noeud semble-t-il)
 * d'ou la decision des procedures iteratives.
 *-------------------------------------------------------------------
 * Algorithme de calcul des D-classes. Principe general :
 * Les R-classes sont les composantes fortement connexes du graphe de Cayley droit
 * Les L-classes sont les composantes fortement connexes du graphe de Cayley gauche
 * On calcule les numeros de R-classes en utilisant l'algorithme de Tarjan 
 * (voir ci-dessous). On calcule ensuite les numeros de L-classes par un parcours
 * du graphe de Cayley gauche. On gere simultanement un tableau TableRD qui donne le
 * numero de D-classe de chaque R-classe. S'il y a par exemple 7 R-classes telles
 * que les R-classes 1 et 4, d'une part, 2, 3 et 6 d'autre part, soient dans la meme
 * D-classe, le tableau final aura la forme suivante :
 * ----------------------------
 * | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
 * -----------------------------
 * | 1 | 2 | 2 | 1 | 3 | 2 | 4 |
 * -----------------------------
 *
 * Au depart, le tableau est initialise a 0.
 * ----------------------------
 * | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
 * -----------------------------
 * | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
 * -----------------------------
 *
 * Lors du calcul des L-classes, on regarde si le sommet courant s est point d'entree
 * d'une L-classe. Si c'est le cas, on regarde le unsigned long de D-classe d = TableRD[R(s)]
 * ou R(s) est le numero de R-classe de s. L'idee est de propager ce numero d dans toute
 * la L-classe de s. Si d = 0, on commence une nouvelle D-classe. On incremente alors
 * le nombre courant de D-classes, et on attribue a d cette nouvelle valeur. Il reste
 * alors a effectuer le parcours de la L-classe.
 *
 */     

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "Main.h"
#include "Memoire.h"
#include "Initialisation.h"
#include "DclassesIteratif.h"

extern info *Table;
extern info2 *Table2;
extern unsigned long *TableIdempotents;
extern unsigned short NbLettres, PossedeUnNeutre, TypeCalcul;
extern unsigned long CalculsEffectues;
extern unsigned long NbElements, NbRclasses, NbLclasses, NbDclasses, NbHclasses, NbIdempotents,
              NbLclassesMinimales, TailleRIdealMinimal, MemoireAllouee;
elementpile *Pile;
unsigned long *PileComposanteConnexe;
unsigned long *TableRD;
 
static unsigned short EstElementRegulier;

/***************************************************************************
*
* ParcoursRAttache. Determine les points d'attache du graphe de Cayley 
* puis les composantes fortement connexes du graphe, ce qui donne les R-classes.
* On utilise l'algorithme de Tarjan.
* 
* s = &Table[s], sa_ = &Table[sa], as = &Table[as]
*
* On effectue un parcours en profondeur d'abord.
* Le point d'attache de s est le minimum de son unsigned long de parcours, des points
* d'attache de ses fils dans l'arbre de parcours et des extremites des arcs transverses
* ou de retour d'origine s.
* Le drapeau R_VISITE est mis a 1 lors du premier passage. Cette information est
* utilisee pour connaitre la nature d'un arc :
* - Arc de descente de s vers sa : (Table2[s].NumeroDeParcours < (*sa_).NumeroDeParcours) 
* - Arc retour de s vers sa : (R_VISITE = 0) en sa
*
* A l'issue de ParcoursRAttache, les points d'entree dans une R-classe sont
* ceux egaux a leur point d'attache.
*
* Gestion de la pile : chaque niveau de pile contient une lettre et la
* valeur cumulee du mot de pile situe dessous. Par exemple, si la pile est : 
* --------------------------------
* | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
* --------------------------------
* | a | b | a | b | b | a | b | a |
* --------------------------------
*
* le niveau 7 contiendra aussi l'adresse de l'element correspondant a ababbab. 
* A la fin du parcours, les champs Pile[i].Lettre sont revenus a 0.
*
***************************************************************************/

void ParcoursRclassesIteratif(void)
{
  register lettre a;
  register long i = 0;
  register unsigned long sa, s, sk;
  unsigned long nPassage = 0;
  unsigned long HauteurPile = 0;
  
  NbRclasses = 0;
  TailleRIdealMinimal = 0;
  Pile[0].Adresse = IDENTITE;
  PileComposanteConnexe[0] = 0;    /* Sentinelle */
  do
  {
    s = Pile[i].Adresse;
    if ((Table2[s].Statut & R_VISITE) == 0)
    {
      Table2[s].PointDAttache = Table2[s].NumeroDeParcours = ++nPassage;  
       PileComposanteConnexe[++HauteurPile] = s;    /* On empile s sur la pile de la composante connexe */
      Table2[s].Statut |= R_VISITE;                /* Premiere visite dans s */
    }
    a = Pile[i].Lettre;
    sa = Table[s].Produits[a].D & ~EST_REDUIT;
    if ((Table2[sa].Statut & R_VISITE) == 0)        /* Si possible, on descend en profondeur. */
      Pile[++i].Adresse = sa;
    else if (a < (NbLettres - 1))        /* Sinon, on explore les autres branches */
    {
      if ((Table2[s].NumeroDeParcours < Table2[sa].NumeroDeParcours) && (Table2[sa].PointDAttache < Table2[s].PointDAttache))  /* Arc de descente de s vers sa */
        Table2[s].PointDAttache = Table2[sa].PointDAttache;
      else if (((Table2[sa].Statut & R_DEJA_CALCULE) == 0) && (Table2[sa].NumeroDeParcours <  Table2[s].PointDAttache))  /* Arc de retour de s vers sa */
        Table2[s].PointDAttache = Table2[sa].NumeroDeParcours;
      Pile[i].Lettre++;
    }
    else                            /* Plus rien a explorer a ce niveau : on met a jour et on remonte */ 
    {
      if ((Table2[s].NumeroDeParcours < Table2[sa].NumeroDeParcours) && (Table2[sa].PointDAttache < Table2[s].PointDAttache))  /* Arc de descente de s vers sa */
        Table2[s].PointDAttache = Table2[sa].PointDAttache;
      else if (((Table2[sa].Statut & R_DEJA_CALCULE) == 0) && (Table2[sa].NumeroDeParcours <  Table2[s].PointDAttache))  /* Arc de retour de s vers sa */
        Table2[s].PointDAttache = Table2[sa].NumeroDeParcours;
      if (Table2[s].NumeroDeParcours == Table2[s].PointDAttache)  /* Si s est egal a son point d'attache, on a un point d'entree de la R-classe */
      {
        NbRclasses++;
        if (NbRclasses == 1)    /* Ideal minimal */
          do
          {
            TailleRIdealMinimal++;               
            sk = PileComposanteConnexe[HauteurPile--];
            Table2[sk].Rclasse = NbRclasses;        /* On depile : on parcours la composante connexe issue de s, ce qui donne une R-classe */
            Table2[sk].Statut |= R_DEJA_CALCULE;    /* Le calcul de la R-classe de s est acheve. On peut eliminer ce sommet. */
          }
          while (sk != s);
        else
          do
          {
            sk = PileComposanteConnexe[HauteurPile--];
            Table2[sk].Rclasse = NbRclasses;        /* On depile : on parcours la composante connexe issue de s, ce qui donne une R-classe */
            Table2[sk].Statut |= R_DEJA_CALCULE;    /* Le calcul de la R-classe de s est acheve. On peut eliminer ce sommet. */
          }
          while (sk != s);
      }
      Pile[i--].Lettre = 0;
    }
  }
  while (i != -1);  
}

/*   if (Table[s].Statut & IDEMPOTENT)        Mise a jour du statut 
      EstElementRegulier = 1;
    if (EstElementRegulier)    
      Table[s].Statut |= REGULIER; 
    if (Table2[s].NumeroDeParcours == Table2[s].PointDAttache)   Lorsque l'un sommet est egal a son point d'attache, on commence une nouvelle R-classe 
      EstElementRegulier = 0; 
    if (((Table[s].Statut & REGULIER) == 0) && ((Table[Table2[s].PointDAttache].Statut & REGULIER) != 0)
      Table[s].Statut |= REGULIER; */



/****************************************************
*
* CalculRclasses. OK
*
****************************************************/

void CalculRclassesIteratif(void)
{
  Table2 = AlloueMemoireTable2();
  Table2[IDENTITE].xOmega = IDENTITE;
  Pile = AlloueMemoirePile();
  PileComposanteConnexe = AlloueMemoirePileComposanteConnexe();
  InitPile(Pile);
  InitDrapeaux();
  EstElementRegulier = 0;
  ParcoursRclassesIteratif();
  CalculsEffectues |= CALCUL_R_CLASSES;    /* On prend note : le calcul des R-classes est termine */
}

/***************************************************************************
*
* ParcoursLAttache. Determine les points d'attache du graphe de Cayley a gauche
*
***************************************************************************/

void ParcoursLclassesIteratif(void)
{
  register lettre a;
  register long i = 0;
  register unsigned long as, s, sk;
  unsigned long nPassage = 0;
  unsigned long HauteurPile = 0;
  unsigned long Dclasse;
  
  if (NbLettres == 0)    /* Monoide trivial */
  {
    NbLclasses = 1;
    Table2[IDENTITE].Lclasse = 1;
    return;
  }
  NbLclasses = 0;  
  NbDclasses = 0;  
  NbLclassesMinimales = 0;  
  Pile[0].Adresse = IDENTITE;
  PileComposanteConnexe[0] = 0;    /* Sentinelle */
  do
  {
    s = Pile[i].Adresse;
    if ((Table2[s].Statut & L_VISITE) == 0)
    {
      Table2[s].PointDAttache = Table2[s].NumeroDeParcours = ++nPassage;  
      PileComposanteConnexe[++HauteurPile] = s;    /* On empile s sur la pile de la composante connexe */
      Table2[s].Statut |= L_VISITE;        /* Premiere visite dans s */
    }
    a = Pile[i].Lettre;
    as = Table[s].Produits[a].G;
    if ((Table2[as].Statut & L_VISITE) == 0)        /* Si possible, on descend en profondeur. */
      Pile[++i].Adresse = as;
    else if (a < (NbLettres - 1))        /* Sinon, on explore les autres branches */
    {
      if ((Table2[s].NumeroDeParcours < Table2[as].NumeroDeParcours) && (Table2[as].PointDAttache < Table2[s].PointDAttache))  /* Arc de descente de s vers as */
        Table2[s].PointDAttache = Table2[as].PointDAttache;
      else if (((Table2[as].Statut & L_DEJA_CALCULE) == 0) && (Table2[as].NumeroDeParcours <  Table2[s].PointDAttache))  /* Arc de retour de s vers as */
        Table2[s].PointDAttache = Table2[as].NumeroDeParcours;
      Pile[i].Lettre++;
    }
    else                            /* Plus rien a explorer a ce niveau : on met a jour et on remonte */ 
    {
      if ((Table2[s].NumeroDeParcours < Table2[as].NumeroDeParcours) && (Table2[as].PointDAttache < Table2[s].PointDAttache))  /* Arc de descente de s vers as */
        Table2[s].PointDAttache = Table2[as].PointDAttache;
      else if (((Table2[as].Statut & L_DEJA_CALCULE) == 0) && (Table2[as].NumeroDeParcours <  Table2[s].PointDAttache))  /* Arc de retour de s vers as */
        Table2[s].PointDAttache = Table2[as].NumeroDeParcours;
      if (Table2[s].NumeroDeParcours == Table2[s].PointDAttache)  /* Si s est egal a son point d'attache, on a un point d'entree de la L-classe */
      {
        NbLclasses++;
        Dclasse = TableRD[Table2[s].Rclasse];
        if (Dclasse == 0)        /* Nouvelle D-classe */
        {
          NbDclasses++;  
          if (NbDclasses == 1)    /* Ideal minimal */
          {
            NbLclassesMinimales++;
            do
            {
              sk = PileComposanteConnexe[HauteurPile--];
              Table2[sk].Lclasse = NbLclasses;        /* On depile : on parcours la composante connexe issue de s, ce qui donne une L-classe */
              Table2[sk].Statut |= L_DEJA_CALCULE | EST_D_MINIMAL;        /* Le calcul de la L-classe de s est acheve. On peut eliminer ce sommet. */
              Table2[sk].Dclasse = 1;
              TableRD[Table2[sk].Rclasse] = 1;
            }
            while (sk != s);
          }
          else                    /* D-classe autre que l'ideal minimal */
          {
            do
            {
              sk = PileComposanteConnexe[HauteurPile--];
              Table2[sk].Lclasse = NbLclasses;        /* On depile : on parcours la composante connexe issue de s, ce qui donne une L-classe */
              Table2[sk].Statut |= L_DEJA_CALCULE;    /* Le calcul de la L-classe de s est acheve. On peut eliminer ce sommet. */
              Table2[sk].Dclasse = NbDclasses;
              TableRD[Table2[sk].Rclasse] = NbDclasses;
            }
            while (sk != s);
          }
        }
        else
        {
          if (Dclasse == 1)    /* Ideal minimal */
          {
            NbLclassesMinimales++;
            do
            {
              sk = PileComposanteConnexe[HauteurPile--];
              Table2[sk].Lclasse = NbLclasses;        /* On depile : on parcours la composante connexe issue de s, ce qui donne une R-classe */
              Table2[sk].Statut |= L_DEJA_CALCULE | EST_D_MINIMAL;    /* Le calcul de la L-classe de s est acheve. On peut eliminer ce sommet. */
              Table2[sk].Dclasse = 1;
            }
            while (sk != s);
          }
          else                    /* D-classe autre que l'ideal minimal */
          {
            do
            {
              sk = PileComposanteConnexe[HauteurPile--];
              Table2[sk].Lclasse = NbLclasses;        /* On depile : on parcours la composante connexe issue de s, ce qui donne une R-classe */
              Table2[sk].Statut |= L_DEJA_CALCULE;    /* Le calcul de la L-classe de s est acheve. On peut eliminer ce sommet. */
              Table2[sk].Dclasse = Dclasse;
            }
            while (sk != s);
          }
        }
      }
      Pile[i--].Lettre = 0;
    }
  }
  while (i != -1);  
}
/*     if (Table[s].Statut & IDEMPOTENT)        Mise a jour du statut 
      EstElementRegulier = 1;
    if (EstElementRegulier)    
      Table[s].Statut |= REGULIER;
    if (Table2[s].NumeroDeParcours == Table2[s].PointDAttache)   Lorsque l'un sommet est egal a son point d'attache, on commence une nouvelle R-classe
      EstElementRegulier = 0;
  if (((Table[s].Statut & REGULIER) == 0) && ((Table[Table2[s].PointDAttache].Statut & REGULIER) != 0)
    Table[s].Statut |= REGULIER; */ 

/****************************************************
*
* CalculDclasses. OK
*
****************************************************/

void CalculDclasses(void)
{
  if (!(CalculsEffectues & CALCUL_D_CLASSES))
  {
    if (NbLettres == 0)    /* Monoide trivial */
    {
      NbDclasses = NbRclasses = NbLclasses = 1;
      Table2[IDENTITE].Dclasse = Table2[IDENTITE].Rclasse = Table2[IDENTITE].Lclasse = 1;
      Table2[IDENTITE].Statut |= EST_D_MINIMAL;
      CalculsEffectues |= CALCUL_D_CLASSES;    /* On prend note : le calcul des D-classes est termine */
      return;
    }
    CalculRclassesIteratif(); 
    EstElementRegulier = 0;
    TableRD = AlloueMemoireTableau(NbRclasses + 1);   /* +1, car on met une sentinelle en fond de pile */
    ParcoursLclassesIteratif();
    if (!PossedeUnNeutre && (TypeCalcul == Semigroupe))
    {
      NbDclasses--;
      NbLclasses--;
      NbRclasses--;
    }
    if (MemoireAllouee & MEMOIRE_PILE)
    {
      free(Pile);
      MemoireAllouee &= ~MEMOIRE_PILE;    /* On prend note : on a libere la memoire de la pile. */
    }
    free(PileComposanteConnexe);
    free(TableRD);
    CalculsEffectues |= CALCUL_D_CLASSES;    /* On prend note : le calcul des D-classes est termine */
  }
}

/****************************************************
*
* CalculReguliers. OK    Aout 2006
*
****************************************************/

void CalculReguliers(void)
{
  unsigned long i;
  unsigned long *D;
  
  if (!(CalculsEffectues & CALCUL_REGULIERS))
  {
    D = AlloueMemoireTableau(NbDclasses + 1);  
    if (CalculsEffectues & CALCUL_LISTE_IDEMPOTENTS)
      for (i = 0; i < NbIdempotents; i++)
        D[Table2[TableIdempotents[i]].Dclasse] |= 1;
    else  
      for (i = IDENTITE + !PossedeUnNeutre; i <= NbElements; i++)
        if (Table[i].Statut & IDEMPOTENT)
          D[Table2[TableIdempotents[i]].Dclasse] |= 1;
    for (i = IDENTITE + !PossedeUnNeutre; i <= NbElements; i++)  
      if (D[Table2[i].Dclasse] != 0)
        Table[i].Statut |= REGULIER;
    free(D);
  }
  CalculsEffectues |= CALCUL_REGULIERS;    /* On prend note : le calcul des elements reguliers est termine */
}

/****************************************************
*
* CalculHclasses. OK    Aout 2006
*
****************************************************/

void CalculHclasses(void)
{
  unsigned long *TableauTrie, *T, *H, *R;
  unsigned long i, j, r, l;
  
/********************************************************************************* 
* Premiere partie: on trie les elements -- en temps lineaire -- par ordre croissant 
*  de numero de R-classe. Le resultat est mis dans TableauTrie. Pour cela, on decompte
* le nombre d'elements dans chaque R-classe, puis on construit le tableau des cumuls
*  qui donne le rang de l'element courant de la R-classe i dans le tableau trie.
* Voici un exemple.
*
* ---------------------------------------------------
*  Element  | 1  2  3  4  5  6  7  8  9 10 11 12 13 |
* ---------------------------------------------------
*  R-classe | 4  4  3  2  3  1  2  3  1  2  3  1  2 |
* ---------------------------------------------------
*
* Nombre d'elements par R-classe.  Tableau cumule (intialisation)
* ------------------------         ------------
*  R-classe    | 1 2 3 4 |         | 1 2 3  4 |
* ------------------------         ------------
*  Nb Elements | 3 4 4 2 |         | 1 4 8 12 |
* ------------------------         ------------
*
* Le tableau trie est alors
*
* ---------------------------------------------------
*  Element  | 1  2  3  4  5  6  7  8  9 10 11 12 13 |
* ---------------------------------------------------
*  TabTrie  | 6  9 12  4  7 10 13  3  5  8 11  1  2 |
* ---------------------------------------------------
*  R-classe | 1  1  1  2  2  2  2  3  3  3  3  4  4 |
* ---------------------------------------------------
*
*********************************************************************************/

  if (!(CalculsEffectues & CALCUL_H_CLASSES))
  {
    TableauTrie = AlloueMemoireTableau(NbElements+1);
    T = AlloueMemoireTableau(NbRclasses+1);  
    for (i = IDENTITE + !PossedeUnNeutre; i <= NbElements; i++)
      T[Table2[i].Rclasse]++;
    T[0] = 1;
    for (i = 0; i < NbRclasses; i++)
      T[i+1] = T[i+1] + T[i];
    for (i = IDENTITE + !PossedeUnNeutre; i <= NbElements; i++)
    {
      j = Table2[i].Rclasse - 1;   /* On a  0 <= j <= NbRclasses - 1 */
      TableauTrie[T[j]] = i;
      T[j]++;
    }
    free(T);      /* Il est preferable de liberer dans l'ordre inverse d'allocation */

/********************************************************************************* 
* Deuxieme partie: on utilise le tableau des L-classes, que l'on va parcourir
* dans l'ordre de TableauTrie : 6, 9, 12, 4, 7, 10, 13, 3, 5, 8, 11, 1, 2
* On s'assure ainsi que les elements d'une meme R-classe sont parcourus consecutivement.
*
* ---------------------------------------------------
*  Element  | 1  2  3  4  5  6  7  8  9 10 11 12 13 |
* ---------------------------------------------------
*  L-classe | 6  6  5  5  3  4  3  5  2  5  3  1  3 |
* ---------------------------------------------------
*
* et on gere deux tableaux auxilaires, de taille NbLclasses: H et R
* Chaque element 6, 9, 12, 4, a un numero de R-classe r et un numero de L-classe l
*  r est stocke dans R[l] et H[l] contient le numero courant de H-classe.
* Le nombre de H-classes est incremente chaque fois que l'on change de R-classe
* par rapport au dernier passage dans la meme L-classe. Au final, on obtient:
*
* ---------------------------------------------------
*  Element  | 1  2  3  4  5  6  7  8  9 10 11 12 13 |
* ---------------------------------------------------
*  H-classe | 8  8  6  4  7  1  5  6  2  4  7  3  5 |
* ---------------------------------------------------
*********************************************************************************/

    NbHclasses = 0;
    H = AlloueMemoireTableau(NbLclasses + 1);  
    R = AlloueMemoireTableau(NbLclasses + 1);  
    for (i = 1; i <= NbElements; i++)
    {
      j = TableauTrie[i];
      r = Table2[j].Rclasse;
      l = Table2[j].Lclasse;
      if (r != R[l])
      {
        NbHclasses++;
        H[l] = NbHclasses;
        R[l] = r;
      }
      Table2[j].Hclasse = H[l];
/*    printf("Table2[j].Hclasse = %d\n", Table2[j].Hclasse); */
    }
    free(R);  
    free(H);
    free(TableauTrie);
  }
  CalculsEffectues |= CALCUL_H_CLASSES;    /* On prend note : le calcul des H-classes est termine */
}

/****************************************************
*
* CalculDansUnGroupe. OK    Aout 2006
*
****************************************************/

void CalculDansUnGroupe(void)
{
  unsigned long i;
  unsigned long *H;
  
  if (!(CalculsEffectues & CALCUL_DANS_UN_GROUPE))
  {
    H = AlloueMemoireTableau(NbHclasses + 1);  
    if (CalculsEffectues & CALCUL_LISTE_IDEMPOTENTS)
      for (i = 0; i < NbIdempotents; i++)
        H[Table2[TableIdempotents[i]].Hclasse] |= 1;
    else  
      for (i = IDENTITE + !PossedeUnNeutre; i <= NbElements; i++)
        if (Table[i].Statut & IDEMPOTENT)
          H[Table2[TableIdempotents[i]].Hclasse] |= 1;
    for (i = IDENTITE + !PossedeUnNeutre; i <= NbElements; i++)  
      if (H[Table2[i].Hclasse] != 0)
        Table[i].Statut |= DANS_UN_GROUPE;
    free(H);
  }
  CalculsEffectues |= CALCUL_DANS_UN_GROUPE;    /* On prend note : le calcul des elements reguliers est termine */
}