Ticket #13744: trac_13744.patch

File trac_13744.patch, 31.5 KB (added by ncohen, 8 years ago)
  • sage/graphs/graph.py

    # HG changeset patch
    # User Nathann Cohen <nathann.cohen@gmail.com>
    # Date 1353726318 -3600
    # Node ID 0c71918867d5e77373fa10420c4fb2eff2d571a1
    # Parent  189737bfb42d7c399ba217d91170b5e93cc6bef5
    Bug in modular_decomposition
    
    diff --git a/sage/graphs/graph.py b/sage/graphs/graph.py
    a b  
    46664666
    46674667        - :meth:`is_prime` -- Tests whether a graph is prime.
    46684668
     4669        TESTS:
     4670
     4671        :trac:`13744`::
     4672
     4673            sage: P = Graph('Dl_')
     4674            sage: P.is_prime()
     4675            False
     4676
    46694677        REFERENCE:
    46704678
    46714679        .. [FMDec] Fabien de Montgolfier
  • sage/graphs/modular_decomposition/src/dm.c

    diff --git a/sage/graphs/modular_decomposition/src/dm.c b/sage/graphs/modular_decomposition/src/dm.c
    a b  
    3030(pour cette notion, cf these de Christian Capelle)
    3131grace a une technique d'affinage de partitions (cf Habib, Paul & Viennot)
    3232Le second construit l'arbre de decomposition modulaire a partir de cette
    33 permutation, cf mon memoire de DEA
    34 Montpellier, decembre 2000
     33permutation, cf
     34[Anne Bergeron, Cedric Chauve, Fabien de Montgolfier, Mathieu Raffinot,
     35 Computing Common Intervals of K Permutations, with Applications to Modular
     36 Decomposition of Graphs. SIAM J. Discrete Math. 22(3): 1022-1039 (2008)]
     37
     38Montpellier, decembre 2000 et Paris, novembre 2010
    3539********************************************************/
    3640
    3741//#include "dm_english.h"
    3842#include <stdio.h>
    3943#include <stdlib.h>
    4044
    41 #define DEBUG 0                 /* si 0 aucune sortie graphique!! */
     45#define DEBUG 0         /* si 0 aucune sortie graphique!! */
    4246
    4347/* dm.h definit les constantes FEUILLE, MODULE, etc...
    4448ainsi que les structures noeud et fils. Les autres
     
    5458  /* On a donc sigma(nom)=place */
    5559  struct Sadj *adj;
    5660  struct SClasse *classe;
     61 
     62
    5763} sommet;
    5864
    5965/* liste d'adjacence d'un sommet, DOUBLEMENT chainee*/
     
    139145  for (s = S[0]->classe; s != NULL; s = s->suiv) {
    140146    printf("[ ");
    141147    for (i = s->debut; i <= s->fin; i++)
    142       printf("%i ", 1 + S[i]->nom);
     148      printf("%i ", S[i]->nom);
    143149    printf("] ");
    144150  }
    145151  printf("\n");
     
    192198  int *ipivot, *imodule, *numclasse, n;
    193199
    194200  if (DEBUG)
    195     printf("Raffinage avec le pivot %i\n", 1 + p->nom);
     201    printf("Raffinage avec le pivot %i\n", p->nom);
    196202  pivot = I->pivot;
    197203  module = I->module;
    198204  ipivot = I->ipivot;
     
    302308          X->inmodule = -1;
    303309          if (DEBUG)
    304310            printf("Dans module %i-%i ecrase %i-%i\n",
    305                    1 + S[Xa->debut]->nom, 1 + S[Xa->fin]->nom,
    306                    1 + S[X->debut]->nom, 1 + S[X->fin]->nom);
     311                    S[Xa->debut]->nom, S[Xa->fin]->nom,
     312                    S[X->debut]->nom, S[X->fin]->nom);
    307313        }
    308314        else {
    309315          if (X->inmodule == -1) {
     
    317323            (*imodule)++;
    318324            if (DEBUG)
    319325              printf("module empile:%i-%i\n",
    320                      1 + S[Z->debut]->nom, 1 + S[Z->fin]->nom);
     326                      S[Z->debut]->nom, S[Z->fin]->nom);
    321327          }
    322328        }
    323329      }
     
    333339      Z->inpivot = (*ipivot);
    334340      (*ipivot)++;
    335341      if (DEBUG)
    336         printf("pivot empile: %i-%i\n", 1 + S[Z->debut]->nom,
    337                1 + S[Z->fin]->nom);
     342        printf("pivot empile: %i-%i\n", S[Z->debut]->nom,
     343                S[Z->fin]->nom);
    338344      X->whereXa = 0;
    339345      Xa->whereXa = 0;
    340346    }
     
    366372    sclasse *C1;                /*premiere classe, tete de la chaine */
    367373    sclasse *Y;                 /*classe qui raffine */
    368374    sclasse *X;                 /*classe raffinee */
    369     sclasse *Xa, *Xc;           /* morceaux de X */
     375    sclasse *Xc;                /* morceau utilise de X */
     376    /* sclasse *Xa, *Xc; morceau inutilise de X */
    370377    sommet *x;                  /* x in X */
    371378    sommet *y;                  /* y in Y */
    372379    sommet *centre;             /* le centre du raffinage actuel */
     
    461468                Y->firstpivot = y;      /* le premier!! */
    462469                if (DEBUG)
    463470                    printf("module-pivot %i-%i: sommet %i\n",
    464                            1 + S[Y->debut]->nom, 1 + S[Y->fin]->nom,
    465                            1 + y->nom);
     471                            S[Y->debut]->nom, S[Y->fin]->nom,
     472                            y->nom);
    466473                Raffiner(S, y, centre, &Inf);
    467474            }
    468475        }
     
    498505
    499506        if (DEBUG)
    500507            printf("Relance dans le module %i-%i avec le sommet %i\n",
    501                    1 + S[X->debut]->nom, 1 + S[X->fin]->nom, 1 + x->nom);
     508                    S[X->debut]->nom,  S[X->fin]->nom, x->nom);
    502509
    503510        centre = x;             /*important! */
    504511        /* astuce: on place {x} en tete de X
     
    537544        if (DEBUG)
    538545            printS(S, n);
    539546    }
     547  free(module);
     548  free(pivot);
    540549}
    541550
    542551/***************************************************************
     
    656665 etablir l'arbre de decomposition modulaire.
    657666 Tous les details sont dans mon memoire de DEA
    658667****************************************************************/
    659 noeud *nouvnoeud(int type, noeud * pere, int sommet, int n)
     668noeud *nouvnoeud(int type, noeud * pere, int sommet, int fv, int lv)
    660669{
    661670    /* cree un nouveau noeud. Noter que l'on est oblige de passer n
    662671       comme parametre car les bords et separateurs droits doivent
     
    669678    nn->type = type;
    670679    nn->pere = pere;
    671680    /* nn->fpere ne peut etre deja mis a jour... */
    672     nn->sommet = sommet;
    673     nn->ps = n + 2;
    674     nn->ds = -2;
    675     /*ces valeurs pour distinguer "non calcule" (-2) */
    676     /*   de "abscence de separateur" (-1). De plus, on fera des min et des */
    677     /*   max sur les bords */
    678     nn->bg = n + 2;
    679     nn->bd = -2;                /* idem */
    680 
     681    nn->nom = sommet;
     682    nn->fv=fv;
     683    nn->lv=lv;
    681684    nn->fils = NULL;
    682685    nn->lastfils = NULL;
    683686    nn->id = compteur;
     
    701704    nfils->fpere = nf;
    702705}
    703706
    704 void fusionne(noeud * pere, noeud * artefact)
    705 {
    706     /*fusionne un artefact a son pere.
    707        utilise le champ fpere qui permet de savoir ou se greffer
    708        une structure fils sera detruite dans l'operation: artefact->fils */
    709     fils *greffe;
    710     fils *f;
    711     /* met a jour la liste des peres */
    712     f = artefact->fils;
    713     while (f != NULL) {
    714         f->pointe->pere = pere; /*avant c'etait ancien... */
    715         /* f->pointe->fpere est inchange */
    716         f = f->suiv;
    717     }
    718     /* greffe la liste */
    719     greffe = artefact->fpere;
    720     artefact->lastfils->suiv = greffe->suiv;
    721     greffe->pointe = artefact->fils->pointe;
    722     greffe->suiv = artefact->fils->suiv;
    723     artefact->fils->pointe->fpere = greffe;     /*artefact->fils a disparu */
    724     if (pere->lastfils == greffe)
    725         pere->lastfils = artefact->lastfils;
    726 }
    727 
    728 void
    729 extraire(noeud * ancien, noeud * nouveau, fils * premier, fils * dernier)
    730 {
    731     /* extrait la liste [premier...dernier] des fils de l'ancien noeud,
    732        et en fait la liste des fils du nouveau noeud */
    733     fils *nf;                   /* il faut une structure fils de plus */
    734     fils *f;                    /*indice de mise a jour */
    735     nf = (fils *) fabmalloc(sizeof(fils));
    736     nf->pointe = premier->pointe;
    737     nf->suiv = premier->suiv;
    738     premier->pointe->fpere = nf;
    739     nouveau->pere = ancien;
    740     nouveau->fils = nf;
    741     nouveau->lastfils = dernier;
    742     nouveau->bg = premier->pointe->bg;
    743     nouveau->bd = dernier->pointe->bd;
    744     nouveau->ps = premier->pointe->bg;  /* nouveau est suppose etre un */
    745     nouveau->ds = dernier->pointe->bd;  /* module, donc bords=separateurs! */
    746     if (ancien->lastfils == dernier)
    747         ancien->lastfils = premier;
    748     /* ecrase l'ancier premier */
    749     nouveau->fpere = premier;
    750     premier->pointe = nouveau;
    751     premier->suiv = dernier->suiv;
    752     /* met a jour dernier */
    753     dernier->suiv = NULL;
    754     /* met a jour la liste des peres */
    755     f = nf;
    756     while (f != dernier->suiv) {
    757         f->pointe->pere = nouveau;      /*avant c'etait ancien... */
    758         f->pointe->fpere = premier;
    759         f = f->suiv;
    760     }
    761 }
    762 
    763707void printnoeud(noeud * N, int level)
    764708{
    765709    /* imprime recursivement l'arbre par parcours en profondeur */
     
    801745            for (i = 0; i < level; i++)
    802746                printf("  |");
    803747            printf("  +--");
    804             printf("%i\n", 1 + nfils->nom);
     748            printf("%i\n", nfils->nom);
    805749        }
    806750        else {
    807751            printnoeud(nfils, level + 1);
     
    816760    printnoeud(N, 0);
    817761}
    818762
     763
     764void typenoeud(noeud * N, int *L2, int *R2, sommet **S)
     765{
     766  // type recursivement l'arbre
     767  if(N->type == FEUILLE)
     768    return;
     769  else if (N->type != UNKN)
     770    {
     771      printf("Erreur typage!\n");
     772      return;
     773    }
     774  // maintenant N est supposé avoir deux fils
     775  noeud *premier, *second;
     776  premier = N->fils->pointe;
     777  if(premier==NULL)
     778    {
     779      printf("Erreur typage!\n");
     780      return;
     781    }
     782  second = N->fils->suiv->pointe;
     783  if(second==NULL)
     784    {
     785      printf("Erreur typage!\n");
     786      return;
     787    }
     788
     789  sadj *a1;
     790  int i = premier->fv;
     791  int j = second->lv;
     792  if(L2[j]<=i && j<= R2[i])
     793    // ah ah les deux premiers fils forment un module
     794    // on recherche si adjacents
     795    {
     796      N->type = PARALLELE; // sera modifie si adjacence
     797      a1=S[i]->adj;
     798      while(a1!=NULL)
     799        {
     800          if(a1->pointe->place == S[j]->place)
     801            // adjacence trouvee !
     802            {
     803              N->type = SERIE;
     804              break;
     805            }
     806          a1=a1->suiv;
     807        }
     808    }
     809  else
     810    N->type = PREMIER;
     811 
     812    fils *ffils;
     813    ffils = N->fils;
     814    do {
     815      typenoeud( ffils -> pointe, L2, R2, S);
     816      ffils = ffils->suiv;
     817    }
     818    while (ffils != NULL);
     819}
     820
     821
    819822noeud *algo2(graphe G, sommet **S)
    820823{
    821824/* algorithme de decomposition modulaire, deuxieme passe
     
    824827*/
    825828    /* debug: S n'est utilise que pour mettre vrainom a jour */
    826829    int n; //nombre de sommets du graphe
    827     int *ouvrantes;             /* tableau du nombre de parentheses ouvrantes */
    828     /* ouvrante[i]=3 ssi i-1(((i  */
    829     /* ouvrante [0]=3: (((0 */
    830 
    831     int *fermantes;             /* idem fermantes[i]=2 ssi i)))i+1
    832                                    fermante [n-1]=2 ssi n))) */
     830                       
    833831    int *ps;                    /* ps[i]=premier separateur de (i,i+1) */
    834832    int *ds;
    835833
     
    837835
    838836    sadj *a1, *a2;                /* parcours de liste d'adjacence */
    839837
    840     noeud *racine;              /*racine du pseudocoardre */
    841     noeud *courant, *nouveau;   /* noeud courant du pseudocoarbre */
    842     noeud **pileinterne;        /* pile des modules pour les passes 3,5,5 */
    843     int indicepileinterne = 0;  /*pointeur dans cette pile */
    844     int taillepileinterne;      /* taille de la pile apres la 2eme passe */
     838   
     839    int *R2, *L2; // le generateur des intervalle-modules produit Capelle-like en phase 2
     840    int *R, *L; // le  generateur canonique des intervalle-modules
     841    int *SupR, *SupL; // le  support pour passer de (R2,L2) a (R,L)
    845842
    846     int *adjii;                 /* adjii[i]=1 ssi S[i] et S[i+1] sont */
    847                                    /*                      adjacents */
     843    int *pile;  // pile utilisee pour calcule generateur
     844    int sommet; // sommet de ladite pile. Indique la derniere place PLEINE
     845
     846    int *tab1,*tab2; // utilise pour le tri lineaire avant algo 6
     847    int *fortg, *fortd; // bornes de mes modules forts
     848    int nb_forts=0; // nombre de modules forts
     849 
     850    int *newfortg, *newfortd; // tableau de modules forts utilises temporairement pour bucket sort
     851 
    848852    /*PROPHASE : initialisations */
    849853    n=G.n;
    850     ouvrantes = (int *) fabmalloc(n * sizeof(int));
    851     fermantes = (int *) fabmalloc(n * sizeof(int));
    852854    ps = (int *) fabmalloc(n * sizeof(int));
    853855    ds = (int *) fabmalloc(n * sizeof(int));
    854     pileinterne = (noeud **) fabmalloc((2 * n + 4) * sizeof(noeud *));
    855     adjii= (int *) fabmalloc(n*sizeof(int));
    856     /*pas plus de 2n+4 noeuds internes dans le pseudocoarbre */
    857     for (i = 0; i < n; i++) {
    858       ouvrantes[i] = 0;
    859       fermantes[i] = 0;
    860       adjii[i]=0;
    861     }
    862 
    863     /* remplit adjii qui dit si S[i] adjacent a S[i+1] */
    864     for(i=0; i<n-1; i++)
    865       {
    866          a1=S[i]->adj;
    867          while((a1!=NULL)&&(a1->pointe->place != i+1))
    868            a1=a1->suiv;
    869          if( a1 == NULL)
    870            adjii[i]=0;
    871          else // a1->pointe->place==i+1, donc i adj i+1
    872            adjii[i]=1;
    873       }
    874     adjii[n-1]=0; //perfectionnisme
     856    R2 = (int *)fabmalloc(n * sizeof(int));
     857    L2 = (int *)fabmalloc(n * sizeof(int));
     858    SupR = (int *)fabmalloc(n * sizeof(int));
     859    SupL = (int *)fabmalloc(n * sizeof(int));
     860    R = (int *)fabmalloc(n * sizeof(int));
     861    L = (int *)fabmalloc(n * sizeof(int));
     862    pile = (int *)fabmalloc(n * sizeof(int));
     863    tab1 = (int *)fabmalloc(4 * n * sizeof(int));
     864    tab2 = (int *)fabmalloc(2 * n * sizeof(int));
     865    fortg=fabmalloc(2 * n * sizeof(int));
     866    fortd=fabmalloc(2 * n * sizeof(int));
     867    newfortg= (int *)fabmalloc(2 * n * sizeof(int));
     868    newfortd=(int *) fabmalloc(2 * n * sizeof(int));
    875869
    876870    /* PREMIERE PASSE
    877871       on va parentheser la permutation factorisante.
    878        tout bonnement, on lit S et on cherche les separateurs;
    879        apres quoi ouvrantes et fermantes sont ecrites
    880        complexite: O(n^2) */
    881 
    882     ouvrantes[0] = 1;
    883     fermantes[n - 1] = 1;       /* parentheses des bords */
     872       complexite: O(m) */
     873    // utilise ps ds
    884874
    885875    for (i = 0; i < n - 1; i++) {
    886876      /*recherche de ps(i,i+1) */
     
    899889          ||((a2==NULL) && (a1->pointe->place >= i))
    900890          ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place >= i) && (a2->pointe->place >= i)))
    901891        //pas de separateur
    902         ps[i]=i+1;
     892        ps[i]=-1;
    903893      else
    904894        {
    905895          if((a1==NULL) || (a1->pointe->place >= i))
     
    915905              else
    916906                ps[i]=min(a1->pointe->place , a2->pointe->place);
    917907            }
    918           ouvrantes[ps[i]]++;   /* marque la fracture gauche, if any */
    919           fermantes[i]++;
    920908        }
    921909      if (DEBUG)
    922         printf("ps(%i,%i)=%i\n", i , i+1, ps[i]);
    923 
     910        printf("ps(%i,%i)=%i ps(%i,%i)=%i\n", i , i+1, ps[i],
     911               S[i]->nom,    S[i + 1]->nom, ps[i]==-1?-1:S[ps[i]]->nom );
    924912        /*recherche de ds(i,i+1)
    925913          plus penible encore!*/
    926914      a1=S[i]->adj;
     
    942930          ||((a2==NULL) && (a1->pointe->place <= i+1))
    943931          ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place <= i+1) && (a2->pointe->place <= i+1)))
    944932        //pas de separateur
    945         ds[i]=i+1;
     933        ds[i]=-1;
    946934      else
    947935        {
    948936          if((a1==NULL) || (a1->pointe->place <= i+1))
     
    961949       
    962950       
    963951          //ds[i] = j;
    964         ouvrantes[i + 1]++;     /* marque la fracture gauche, if any */
    965         fermantes[ds[i]]++;     /* attention aux decalages d'indices */
    966952        }
    967953      if (DEBUG)
    968         printf("ds(%i,%i)=%i\n", i,i+1,ds[i]);
    969       //S[i]->nom + 1,         S[i + 1]->nom + 1, S[ds[i]]->nom + 1);
    970     }
    971  
    972     /*DEUXIEME PASSE: construction du pseudocoarbre */
    973 
    974     racine = nouvnoeud(UNKN, NULL, -1, n);
    975     courant = racine;
    976     for (i = 0; i < n; i++) {
    977         /*1: on lit des parentheses ouvrantes: descentes */
    978         for (j = 0; j < ouvrantes[i]; j++) {
    979             /*Descente vers un nouveau noeud */
    980             nouveau = nouvnoeud(UNKN, courant, -1, n);
    981             ajoutfils(courant, nouveau);        /*on l'ajoute... */
    982             courant = nouveau;  /* et on descent */
    983             if (DEBUG)
    984                 printf("(");
    985         }
    986         /* 2: on lit le sommet: feuille */
    987         nouveau = nouvnoeud(FEUILLE, courant, i, n);
    988         ajoutfils(courant, nouveau);    /*on ajoute le bebe... */
    989         (*nouveau).ps = i;
    990         (*nouveau).ds = i;
    991         (*nouveau).bg = i;
    992         (*nouveau).bd = i;
    993         nouveau->nom = S[i]->nom;
    994         /* et pourquoi pas i ? Je m'embrouille... */
    995         if (DEBUG)
    996             printf(" %i ", S[i]->nom + 1);
    997 
    998         /*3: on lit des parentheses fermantes: remontees */
    999         for (j = 0; j < fermantes[i]; j++) {
    1000             /*ASTUCE: ici on va en profiter pour supprimer les
    1001                noeuds a un fils, afin d'economiser une passe */
    1002             if (courant->fils == courant->lastfils) {   /*just one son */
    1003                 courant->pere->lastfils->pointe = courant->fils->pointe;
    1004                 courant->fils->pointe->pere = courant->pere;
    1005                 courant->fils->pointe->fpere = courant->pere->lastfils;
    1006                 /* explication: le dernier fils de courant.pere est
    1007                    actuellement courant himself. Il est donc pointe par
    1008                    courant.pere.lastfils.pointe. Il suffit de changer ce
    1009                    pointeur pour qu'il pointe maintenant non plus sur courant,
    1010                    mais sur l'unique fils de courant: courant.fils.pointe.
    1011                    Ouf! */
    1012                 /* NB: courant est maintenant deconnecte de l'arbre.
    1013                    on pourrait faire un free() mais bon... */
    1014             }
    1015             else {
    1016                 /*on empile ce noeud interne.
    1017                    L'ordre est celui de la postvisite d'un DFS */
    1018                 pileinterne[indicepileinterne] = courant;
    1019                 indicepileinterne++;
    1020             }
    1021             /* et dans tous les cas: on remonte! */
    1022             courant = courant->pere;
    1023             if (DEBUG)
    1024                 printf(")");
    1025         }
    1026     }
    1027     if (DEBUG)
    1028         printf("\n");
    1029 
    1030     /* on enleve un ptit defaut */
    1031     racine = racine->fils->pointe;
    1032     racine->pere = NULL;
    1033     racine->fpere = NULL;
    1034     if (DEBUG) {
    1035         printf("Arbre apres la deuxieme passe:\n");
    1036         printnoeud(racine, 0);
     954        printf("ds(%i,%i)=%i ds(%i,%i)=%i\n", i,i+1,ds[i],
     955               S[i]->nom,   S[i + 1]->nom, ds[i]==-1?-1:S[ds[i]]->nom );
    1037956    }
    1038957
    1039     /*TROISIEME PASSE */
    1040     /* A ce stade, on a un pseudocoarbre de racine racine,
    1041        sans noeuds a un fils, et dont les etiquettes sont
    1042        FEUIILLE ou UNKN. Il y a une pile des noeuds UNKN, stockes
    1043        dans l'ordre de la postvisite d'un parcours en profondeur.
    1044        On va s'en servir pour faire remonter les bords et les
    1045        separateurs de bas en haut */
     958    /* DEUXIEMME PASSE
     959       Calcule le generateur des interval-modules comme algo 9 du papier DAM
     960    */       
     961    int v;
     962    // utilise ps ds L2 R2 pile
    1046963
    1047     taillepileinterne = indicepileinterne;
    1048     for (indicepileinterne = 0; indicepileinterne < taillepileinterne;
    1049          indicepileinterne++) {
    1050         noeud *scanne;
    1051         fils *nextfils;
    1052         noeud *succourant;
    1053         /* scanne est le noeud (pere) que l'on regarde;
    1054            nextfils parcours sa liste de fils;
    1055            courant est le fils actuellement examine et
    1056            succourant=succ(courant) */
    1057         noeud *debutnoeud;
    1058         fils *debutfils;
    1059         /*deux variables utilise pour la recherche de jumeaux:
    1060            debut du bloc max */
     964    // algo 9
     965    // Attention s[v] dans l'algo est  ds[v-1] !!
     966    sommet = -1;
     967   
     968    for(v = n-1; v>=0; v--)
     969      if(ds[v-1] != -1)
     970        {
     971          L2[v]=v;
     972          while( pile[sommet] < ds[v-1])
     973            L2[pile[sommet--]]=v;
     974        }       
     975      else
     976        pile[++sommet]=v;
    1061977
    1062         scanne = pileinterne[indicepileinterne];        /*he oui, la pile */
    1063         nextfils = scanne->fils;        /*on commence au premier fils */
    1064         do {
    1065             /*la boucle chiante... cf mon memoire de DEA */
    1066             courant = nextfils->pointe;
    1067             /* bords */
    1068             scanne->bg = min(scanne->bg, courant->bg);
    1069             scanne->bd = max(scanne->bd, courant->bd);
    1070             /*separateurs */
    1071             scanne->ps = min(scanne->ps, courant->ps);
    1072             if (scanne->fils->pointe != courant)
    1073                 /*ce n'est pas le premier fils */
    1074                 scanne->ps = min(scanne->ps, ps[(courant->bg) - 1]);
    1075             scanne->ds = max(scanne->ds, courant->ds);
    1076             if (scanne->lastfils->pointe != courant)
    1077                 /*ce n'est pas le dernier fils */
    1078                 scanne->ds = max(scanne->ds, ds[courant->bd]);
     978    while(sommet>=0)
     979      L2[pile[sommet--]]=0;
    1079980
    1080             nextfils = nextfils->suiv;
    1081         }
    1082         while (nextfils != NULL);
     981    sommet = -1;
     982    // algo 9 symétrique
     983    //re-attention là c'est ps[v]
     984    for(v = 0 ; v<=n-1; v++)
     985      if(ps[v] != -1)
     986        {
     987          R2[v]=v;
     988          while(  pile[sommet] > ps[v])
     989            R2[pile[sommet--]]=v;
     990        }       
     991      else
     992        pile[++sommet]=v;
     993   
     994    while(sommet>=0)
     995      R2[pile[sommet--]]=n-1;
    1083996
     997    if (DEBUG)
     998      {    printf("L2 = [");
     999        for(v = 0;v<n;v++)
     1000          printf("%i ",L2[v]);
     1001        printf("]\nR2 = [");
     1002        for(v = 0;v<n;v++)
     1003          printf("%i ",R2[v]);
     1004        printf("]\n");
     1005        for(i=0;i<n;i++)
     1006          for(j=i+1;j<n;j++)
     1007            if(L2[j]<=i && j<= R2[i])
     1008              printf("[%i-%i] module\n",i,j);
    10841009
    1085         if (DEBUG)
    1086             printf("noeud %i-%i: ps=%i ds=%i", 1 + scanne->bg,
    1087                    1 + scanne->bd, 1 + scanne->ps, 1 + scanne->ds);
     1010      }
    10881011
    1089         /* maintenant le test tout simple pour savoir si on a un module: */
    1090         if (((scanne->ps) == (scanne->bg))
    1091             && ((scanne->ds) == (scanne->bd))) {
    1092             /*on a un module */
    1093             scanne->type = MODULE;
    1094             if (DEBUG)
    1095                 printf(" Module.\n");
    1096         }
    1097         else {
    1098             scanne->type = ARTEFACT;
    1099             if (DEBUG)
    1100                 printf(" artefact.\n");
    1101         }
    1102     }
    11031012
    1104     if (DEBUG) {
    1105         printf("Arbre apres la troisieme passe:\n");
    1106         printnoeud(racine, 0);
    1107     }
     1013   
     1014    /* TROISIEME PASSE
     1015       generateur canonique. Algos 3 et 5 du papier de DAM
     1016    */
     1017 
     1018    // Algo 3, R
     1019    pile[0]=0;
     1020    sommet =0;
     1021    for(i=1;i<n;i++)
     1022      {
     1023        while(R2[pile[sommet]] < i)
     1024          sommet --;
     1025        SupR[i]=pile[sommet];
     1026        pile[++sommet]=i;
     1027      }
     1028   
     1029    // Algo 3, L
     1030 
     1031    pile[0]=n-1;
     1032    sommet =0;
     1033    for(i=n-2;i>=0;i--)
     1034      {
     1035        while(L2[pile[sommet]] > i)
     1036          sommet --;
     1037        SupL[i]=pile[sommet];
     1038        pile[++sommet]=i;
     1039      }
     1040   
     1041    if (DEBUG)
     1042      {
     1043        printf("SupL = [");
     1044        for(v = 0;v<n;v++)
     1045          printf("%i ",SupL[v]);
     1046        printf("]\nSupR = [");
     1047        for(v = 0;v<n;v++)
     1048          printf("%i ",SupR[v]);
     1049        printf("]\n");
     1050      }
     1051 
     1052    // Algo 5, R
     1053    int k;
     1054    R[0]=n-1;
     1055    for(k=1;k<n;k++)
     1056      R[k]=k;
     1057    for(k=n-1;k>=1;k--)
     1058      if(L2[R[k]] <= SupR[k] && R[k] <= R2[SupR[k]])
     1059        R[SupR[k]] = max(R[k], R[SupR[k]]);
     1060 
     1061    // Algo 5, L
     1062    L[n-1]=0;
     1063    for(k=0;k<n-1;k++)
     1064      L[k]=k;
     1065    for(k=0;k<n-1;k++)
     1066      if(L2[SupL[k]] <= L[k] && SupL[k] <= R2[L[k]])
     1067        L[SupL[k]] = min(L[k], L[SupL[k]]);
     1068 
     1069    if (DEBUG)
     1070      {    printf("L = [");
     1071        for(v = 0;v<n;v++)
     1072          printf("%i ",L[v]);
     1073        printf("]\nR = [");
     1074        for(v = 0;v<n;v++)
     1075          printf("%i ",R[v]);
     1076        printf("]\n");
     1077        for(i=0;i<n;i++)
     1078          for(j=i+1;j<n;j++)
     1079            if(L[j]<=i && j<= R[i])
     1080              printf("[%i-%i] module\n",i,j);
    11081081
    1109     /* QUATRIEME PASSE */
    1110     /* technique:on se contente de fusionner les artefacts a leurs parents
    1111        ca se fait de bas en haut grace a pileinterne (toujours elle!) */
     1082      }
     1083   
     1084    /* QUATRIEME PASSE
     1085       strong modules. Algo 6 du papier de DAM
     1086    */
     1087    //1 remplit le tableau des bornes
     1088    // astuce : on met 2a pour a borne gauche et 2a+1 pour a borne droite. Puis tri lineaire
    11121089
    1113     for (indicepileinterne = 0; indicepileinterne < taillepileinterne;
    1114          indicepileinterne++) {
    1115         noeud *scanne;
    1116         scanne = pileinterne[indicepileinterne];        /*he oui, la pile */
    1117         if (scanne->type == ARTEFACT) {
    1118             /*attention! La pile peut contenir des noeuds deconnectes */
    1119             fusionne(scanne->pere, scanne);
    1120             if (DEBUG)
    1121                 printf("Artefact elimine: %i-%i\n", 1 + scanne->bg,
    1122                        1 + scanne->bd);
    1123         }
    1124     }
    1125     if (DEBUG) {
    1126         printf("Arbre apres la quatrieme passe:\n");
    1127         printnoeud(racine, 0);
    1128     }
     1090    for(i = 0; i<n; i++)
     1091      {
     1092        // intervalle (L[i]..i)
     1093        tab1[4*i] = 2*L[i];
     1094        tab1[4*i+1] = 2*i+1;
     1095        // intervalle (i..R[i])
     1096        tab1[4*i+2] = 2*i;
     1097        tab1[4*i+3] = 2*R[i]+1;
     1098      }
     1099    // tableau des fréquences
     1100 
     1101    for( i=0;i<2*n;i++)
     1102      tab2[i]=0;
     1103    for( i=0;i<4*n;i++)
     1104      tab2[tab1[i]]++;
     1105    // tri
     1106    j = 0;
     1107    for( i=0;i<2*n;i++)
     1108      while(tab2[i]-->0)
     1109        {
     1110          tab1[j++] = i;
     1111        }
    11291112
    1130     /* CINQIEME ET DERNIERE PASSE */
    1131     /* on va typer les noeuds et extraire les fusions.
    1132        comment on fait? Ben, la pile.... */
    1133     for (indicepileinterne = 0; indicepileinterne < taillepileinterne;
    1134          indicepileinterne++) {
    1135         noeud *scanne;
    1136         fils *nextfils;
    1137         noeud *succourant;
    1138         /* scanne est le noeud (pere) que l'on regarde;
    1139            nextfils parcours sa liste de fils;
    1140            courant est le fils actuellement examine et
    1141            succourant=succ(courant) */
    1142         noeud *debutnoeud;
    1143         fils *debutfils;
    1144         /*deux variables utilise pour la recherche de jumeaux:
    1145            debut du bloc max */
     1113    // Algo 6 du papier de DAM
     1114    // j'utilise maintenant tab2 comme stack de sommet : sommet
     1115    sommet = -1;
     1116 
     1117   
     1118    for( i=0;i<4*n;i++)
     1119      {
     1120        int ai,s; //meme noms de variables que papier
     1121        ai=tab1[i];
     1122        if((ai%2)==0)
     1123          tab2[++sommet] = ai;
     1124        else
     1125          {
     1126            ai/=2;
     1127            s = tab2[sommet--]/2;
     1128         
     1129            if(nb_forts == 0 || (fortg[nb_forts - 1] != s) || (fortd[nb_forts - 1] != ai))
     1130              {
     1131                fortg[nb_forts] = s;
     1132                fortd[nb_forts] = ai;
     1133                nb_forts++;
     1134              }
     1135          }
     1136      }
     1137    if(DEBUG)   
     1138      {
     1139        printf("Modules forts :\n");
     1140        for( i=0;i<nb_forts;i++)
     1141          printf("[%i %i] ", fortg[i], fortd[i]);
     1142        printf("\n");
     1143      }
     1144   
     1145    // TRI des modules forts par bornes droites decroissantes
     1146    for( i=0;i<n;i++)
     1147      tab1[i]=0;
     1148    for( i=0;i<nb_forts;i++)
     1149      tab1[fortd[i]]++;
     1150    for( i=1;i<n;i++)
     1151      tab1[i]+=tab1[i-1];
    11461152
    1147         scanne = pileinterne[indicepileinterne];
    1148         if (scanne->type != MODULE)
    1149             continue;           /* je traite que les modules */
     1153    // tab1 = tableau des frequences cumulées
     1154    for( i=0;i<nb_forts;i++)
     1155      {
     1156        int pos = (nb_forts-1)-(tab1[fortd[i]]-1);
     1157        newfortg[ pos ] = fortg[i];
     1158        newfortd[ pos ] = fortd[i];
     1159        tab1[fortd[i]]--;
     1160      }
     1161    if(DEBUG)   
     1162      {
     1163        printf("Modules forts :\n");
     1164        for( i=0;i<nb_forts;i++)
     1165          printf("[%i %i] ", newfortg[i], newfortd[i]);
     1166        printf("\n");
     1167      }
    11501168
    1151         nextfils = scanne->fils;        /*on commence au premier fils */
    1152         while (1) {
    1153             courant = nextfils->pointe;
    1154             succourant = nextfils->suiv->pointe;
    1155             if (ps[courant->bd] >= courant->bg
    1156                 && ds[courant->bd] <= succourant->bd) {
    1157                 /*Des jumeaux!! ->serie ou parallele! */
    1158                 /* on va determiner le bloc max de jumeaux consecutifs */
    1159                 debutnoeud = courant;
    1160                 debutfils = nextfils;
    1161                 while (ps[courant->bd] >= courant->bg &&
    1162                        ds[courant->bd] <= succourant->bd &&
    1163                        nextfils->suiv != NULL) {
    1164                     nextfils = nextfils->suiv;
    1165                     courant = nextfils->pointe;
    1166                     if (nextfils->suiv == NULL)
    1167                         break;
    1168                     succourant = nextfils->suiv->pointe;
    1169                 }
    1170                 /*maintenant on connait la taille du bloc: il va de
    1171                    debutnoeud a courant inclus,
    1172                    et dans la liste des fils de scanne,
    1173                    il s'etend de debutfils a nextfils inclus.
    1174                    On extrait cette liste pour en faire les fils d'un
    1175                    nouveau noeud... sauf si pas la peine! */
    1176                 if (debutfils == scanne->fils
    1177                     && nextfils == scanne->lastfils) {
    1178                     /* le noeud scanne etait serie ou parallele */
    1179                   if ( adjii[debutnoeud->bd] !=0)
    1180                         scanne->type = SERIE;
    1181                     else
    1182                         scanne->type = PARALLELE;
    1183                 }
    1184                 else {
    1185                   if ( adjii[debutnoeud->bd]!=0)
    1186                     /*seule cette ligne distingue G de G' !! */
    1187                     {
    1188                         nouveau = nouvnoeud(SERIE, scanne, -1, n);
    1189                         if (DEBUG)
    1190                             printf("Module serie extrait: %i-%i\n",
    1191                                    1 + debutnoeud->bg, 1 + courant->bd);
    1192                     }
    1193                     else {
    1194                         nouveau = nouvnoeud(PARALLELE, scanne, -1, n);
    1195                         if (DEBUG)
    1196                             printf("Module parallele extrait: %i-%i\n",
    1197                                    1 + debutnoeud->bg, 1 + courant->bd);
    1198                     }
    1199                     extraire(scanne, nouveau, debutfils, nextfils);
    1200                 }
    1201             }
    1202             if (scanne->type == MODULE)
    1203                 scanne->type = PREMIER;
    1204             if (nextfils->suiv == NULL || nextfils->suiv->suiv == NULL
    1205                 || nextfils->suiv->suiv->suiv == NULL)
    1206                 break;
    1207             nextfils = nextfils->suiv;
    1208         }
    1209     }
    1210     if (DEBUG) {
    1211         printf("Arbre final:\n");
    1212         printnoeud(racine, 0);
    1213     }
     1169   // now tri par bornes gauches croissantes
     1170    for( i=0;i<n;i++)
     1171      tab1[i]=0;
     1172    for( i=0;i<nb_forts;i++)
     1173      tab1[newfortg[i]]++;
     1174    for( i=1;i<n;i++)
     1175      tab1[i]+=tab1[i-1];
     1176   
     1177    // tab1 = tableau des frequences cumulées
     1178    for( i = nb_forts-1 ; i>=0 ; i--)// parcours de droite a gauche car on retourne tout
     1179      {
     1180        int pos = tab1[newfortg[i]]-1;
     1181        // printf("[%i %i] en pos %i\n",newfortg[i],newfortd[i],pos);
     1182        fortg[ pos ] = newfortg[i];
     1183        fortd[ pos ] = newfortd[i];
     1184        tab1[newfortg[i]]--;
     1185      }
     1186    if(DEBUG)   
     1187      {
     1188        printf("Modules forts :\n");
     1189        for( i=0;i<nb_forts;i++)
     1190          printf("[%i %i] ", fortg[i], fortd[i]);
     1191        printf("\n");
     1192      }
     1193
     1194    // Finally algo 7 du papier : construit l'arbre
     1195    noeud* racine = nouvnoeud( UNKN, NULL, -1, 0, n-1);
     1196
     1197    noeud *F = racine;
     1198    noeud *Nouv;
     1199    for(k=1; k< nb_forts; )
     1200      {
     1201        if(F->fv <= fortg[k] && F->lv >= fortd[k] )
     1202          {     
     1203            if(fortg[k]==fortd[k]) // on a une feuille
     1204              Nouv = nouvnoeud(FEUILLE, F, S[fortg[k]]->nom, fortg[k], fortd[k]);
     1205            else
     1206              Nouv = nouvnoeud(UNKN, F, -1,  fortg[k], fortd[k]);
     1207            ajoutfils(F, Nouv);
     1208            // printf("nouv [%i %i] de pere [%i %i]\n",
     1209            // fortg[k], fortd[k], F->fv, F->lv);
     1210            F = Nouv;
     1211            k++;
     1212           
     1213          }
     1214        else
     1215          {
     1216            // printf("remonte car [%i %i] pas fils de [%i %i]\n",
     1217            // fortg[k], fortd[k], F->fv, F->lv);
     1218            F = F->pere;
     1219          }
     1220      }
     1221    if(DEBUG)
     1222      {
     1223        printf("arbre avant typage :\n");
     1224        printarbre(racine);
     1225      }
     1226
     1227    /* CINQUIEME PASSE : typage */
     1228    // Methode : on teste grace au generateur la modularite des deux premiers fils.
     1229    // Ensuite si module est serie si adjacence parallele sinon
     1230    typenoeud(racine, L2, R2, S);
     1231
     1232    if(DEBUG)
     1233      {
     1234        printf("\n arbre final :\n");
     1235        printarbre(racine);
     1236      }
     1237    /* nettoie un peu */
     1238    free(ps);
     1239    free(ds);
     1240    free(R2);
     1241    free(L2);
     1242    free(SupR);
     1243    free(SupL);
     1244    free(R);
     1245    free(L);
     1246    free(pile);
     1247    free(tab1);
     1248    free(tab2);
     1249    free(fortg);
     1250    free(fortd);
     1251    free(newfortg);
     1252    free(newfortd);
     1253
     1254
    12141255    return racine;
    12151256}
    12161257
     
    12521293     /* affiche la permutation factorisante */
    12531294{
    12541295  int i;
    1255   printf(  "Place (nouvelle num) ");
     1296  printf("Attention dans debug maintenant usage des NOUVEAUX noms\n");
     1297  printf(  "Place (nouveaux noms, nouvelle numerotation) ");
    12561298  for(i=0;i<n;i++)
    12571299    printf("%3i ",S[i]->place);
    1258   printf("\nNom (ancienne num) : ");
     1300   printf("\nNom (ancienne num) :                        ");
    12591301  for(i=0;i<n;i++)
    12601302    printf("%3i ",S[i]->nom);
    12611303  printf("\n");
    12621304}
    12631305
     1306void nettoie(sommet **S, int n)
     1307{
     1308  // Detruit la copie locale du graphe
     1309  // normalement, apres, tous les mallocs sont
     1310  // annules sauf ceux de l'arbre (le resultat)
     1311  int i;
     1312  sommet * scourant;
     1313  sadj *a,*atmp;
     1314  for (i = 0; i < n; i++) {
     1315        /* initialisation des sommets */
     1316        /* notre bebe est le sommet i dans M */
     1317    scourant = S[i];
     1318    a=scourant->adj;
     1319      while(a!=NULL)
     1320        {
     1321          atmp=a->suiv;
     1322          free(a);
     1323          a=atmp;
     1324        }
     1325      free(scourant->classe); // efface toutes les classes car sont singletons
     1326      free(scourant);
     1327  }
     1328  free(S);
     1329}
    12641330
    12651331/* la fonction principale; qui fait pas grand'chose....*/
    12661332noeud *decomposition_modulaire(graphe G)
     
    12821348        PrintS2(S,G.n);
    12831349      }
    12841350    Racine = algo2(G, S);       /* deuxieme partie: calcul de l'arbre */
     1351
     1352    nettoie(S, G.n);
     1353
    12851354    return Racine;
    12861355}
    12871356
  • sage/graphs/modular_decomposition/src/dm_english.h

    diff --git a/sage/graphs/modular_decomposition/src/dm_english.h b/sage/graphs/modular_decomposition/src/dm_english.h
    a b  
    6666Each internal node is labelled SERIE (for series), PARALLELE (for parallel) or PREMIER (for prime) depending of the quotient's type.
    6767Each leaf is labelled FEUILLE and also contains the vertex number of the leaf.
    6868As the tree is an inclusion tree, the vertex-set corresponding to an internal node correspond to the vertices numbers of the leaves that descend from that tree. The function decomposition_modulaire() return a pointer to the root of the tree.
    69 
     69*/
    7070
    7171
    7272/* define the type of nodes. UNKN,MODULE,ARTEFACT are for internal use*/
     
    8585  int type;     // is FEUILLE, SERIE, PARALLELE or PREMIER             
    8686  struct Noeud *pere;   // adress of parent node, NULL if root
    8787  struct Fils *fpere;   // points the head of the linked list of sons (if type is not FEUILLE, else is NULL)
    88   int ps;       // internal use
    89   int bg;       // internal use
    90   int ds;       // internal use
    91   int bd;       // internal use
    92   int sommet;   // internal use
    9388  int nom;      // if type=FEUILLE, number of the corresponding vertex of the graph
    9489  struct Fils *fils; // points the head of the linked list of sons
    9590  struct Fils *lastfils;  // internal use (points the last item in the listed list of sons)
    9691  int id;       // internal use (node unique ID)
     92  int fv;  // internal use (first vertex in factorizing permutation)
     93  int lv;  // internal use (last vertex in factorizing permutation)
     94
    9795} noeud;
    9896
    9997/* linked list that strore the sons of an internal node (in any order) */
  • sage/graphs/modular_decomposition/src/random.c

    diff --git a/sage/graphs/modular_decomposition/src/random.c b/sage/graphs/modular_decomposition/src/random.c
    a b  
    6767    {
    6868      if(VERBEUX)
    6969        {
    70           printf("Ajacence de %i: ",i+1);
     70          printf("Ajacence de %i: ",i);
    7171          for(a=G.G[i]; a!=NULL; a=a->suiv)
    72             printf("%i ",1+a->s);
     72            printf("%i ",a->s);
    7373        }
    7474      for(j=i+1;j<n;j++)
    7575        {
     
    8686            a->suiv=G.G[j];
    8787            G.G[j]=a;
    8888            if(VERBEUX)
    89               printf("%i ",j+1);
     89              printf("%i ",j);
    9090          }
    9191        }
    9292      if(VERBEUX)