Ticket #13744: 13744_rebased.patch

File 13744_rebased.patch, 35.4 KB (added by jdemeyer, 8 years ago)

  • sage/graphs/graph.py 

    # HG changeset patch
# User Jeroen Demeyer <jdemeyer@cage.ugent.be>
# Date 1358416461 -3600
# Node ID cb7d2191b5c56b89609eb5687573d363be2e3397
# Parent  e610e237c6fd5cd0469a8428b99e27d16f1c56f7
Bug in modular_decomposition

diff --git a/sage/graphs/graph.py b/sage/graphs/graph.py
    a b  
    53095309
    53105310        - :meth:`is_prime` -- Tests whether a graph is prime.
    53115311
     5312        TESTS:
     5313
     5314        :trac:`13744`::
     5315
     5316            sage: P = Graph('Dl_')
     5317            sage: P.is_prime()
     5318            False
     5319
    53125320        REFERENCE:
    53135321
    53145322        .. [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
    … …  
    139143  for (s = S[0]->classe; s != NULL; s = s->suiv) {
    140144    printf("[ ");
    141145    for (i = s->debut; i <= s->fin; i++)
    142       printf("%i ", 1 + S[i]->nom);
     146      printf("%i ", S[i]->nom);
    143147    printf("] ");
    144148  }
    145149  printf("\n");
    … …  
    192196  int *ipivot, *imodule, *numclasse, n;
    193197
    194198  if (DEBUG)
    195     printf("Raffinage avec le pivot %i\n", 1 + p->nom);
     199    printf("Raffinage avec le pivot %i\n", p->nom);
    196200  pivot = I->pivot;
    197201  module = I->module;
    198202  ipivot = I->ipivot;
    … …  
    302306          X->inmodule = -1;
    303307          if (DEBUG)
    304308            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);
     309                    S[Xa->debut]->nom, S[Xa->fin]->nom,
     310                    S[X->debut]->nom, S[X->fin]->nom);
    307311        }
    308312        else {
    309313          if (X->inmodule == -1) {
    … …  
    317321            (*imodule)++;
    318322            if (DEBUG)
    319323              printf("module empile:%i-%i\n",
    320                      1 + S[Z->debut]->nom, 1 + S[Z->fin]->nom);
     324                      S[Z->debut]->nom, S[Z->fin]->nom);
    321325          }
    322326        }
    323327      }
    … …  
    333337      Z->inpivot = (*ipivot);
    334338      (*ipivot)++;
    335339      if (DEBUG)
    336         printf("pivot empile: %i-%i\n", 1 + S[Z->debut]->nom,
    337                1 + S[Z->fin]->nom);
     340        printf("pivot empile: %i-%i\n", S[Z->debut]->nom,
     341                S[Z->fin]->nom);
    338342      X->whereXa = 0;
    339343      Xa->whereXa = 0;
    340344    }
    … …  
    366370    sclasse *C1;                /*premiere classe, tete de la chaine */
    367371    sclasse *Y;                        /*classe qui raffine */
    368372    sclasse *X;                        /*classe raffinee */
    369     sclasse *Xa, *Xc;                /* morceaux de X */
     373    sclasse *Xc;                /* morceau utilise de X */
     374    /* sclasse *Xa, *Xc; morceau inutilise de X */
    370375    sommet *x;                        /* x in X */
    371376    sommet *y;                        /* y in Y */
    372377    sommet *centre;                /* le centre du raffinage actuel */
    … …  
    461466                Y->firstpivot = y;        /* le premier!! */
    462467                if (DEBUG)
    463468                    printf("module-pivot %i-%i: sommet %i\n",
    464                            1 + S[Y->debut]->nom, 1 + S[Y->fin]->nom,
    465                            1 + y->nom);
     469                            S[Y->debut]->nom, S[Y->fin]->nom,
     470                            y->nom);
    466471                Raffiner(S, y, centre, &Inf);
    467472            }
    468473        }
    … …  
    498503
    499504        if (DEBUG)
    500505            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);
     506                    S[X->debut]->nom, S[X->fin]->nom, x->nom);
    502507
    503508        centre = x;                /*important! */
    504509        /* astuce: on place {x} en tete de X
    … …  
    537542        if (DEBUG)
    538543            printS(S, n);
    539544    }
     545  free(module);
     546  free(pivot);
    540547}
    541548
    542549/***************************************************************
    … …  
    656663 etablir l'arbre de decomposition modulaire.
    657664 Tous les details sont dans mon memoire de DEA
    658665****************************************************************/
    659 noeud *nouvnoeud(int type, noeud * pere, int sommet, int n)
     666noeud *nouvnoeud(int type, noeud * pere, int sommet, int fv, int lv)
    660667{
    661668    /* cree un nouveau noeud. Noter que l'on est oblige de passer n
    662669       comme parametre car les bords et separateurs droits doivent
    … …  
    669676    nn->type = type;
    670677    nn->pere = pere;
    671678    /* 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 
     679    nn->nom = sommet;
     680    nn->fv=fv;
     681    nn->lv=lv;
    681682    nn->fils = NULL;
    682683    nn->lastfils = NULL;
    683684    nn->id = compteur;
    … …  
    701702    nfils->fpere = nf;
    702703}
    703704
    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 
    763705void printnoeud(noeud * N, int level)
    764706{
    765707    /* imprime recursivement l'arbre par parcours en profondeur */
    … …  
    801743            for (i = 0; i < level; i++)
    802744                printf("  |");
    803745            printf("  +--");
    804             printf("%i\n", 1 + nfils->nom);
     746            printf("%i\n", nfils->nom);
    805747        }
    806748        else {
    807749            printnoeud(nfils, level + 1);
    … …  
    816758    printnoeud(N, 0);
    817759}
    818760
     761
     762void typenoeud(noeud * N, int *L2, int *R2, sommet **S)
     763{
     764  // type recursivement l'arbre
     765  if(N->type == FEUILLE)
     766    return;
     767  else if (N->type != UNKN)
     768    {
     769      printf("Erreur typage!\n");
     770      return;
     771    }
     772  // maintenant N est supposé avoir deux fils
     773  noeud *premier, *second;
     774  premier = N->fils->pointe;
     775  if(premier==NULL)
     776    {
     777      printf("Erreur typage!\n");
     778      return;
     779    }
     780  second = N->fils->suiv->pointe;
     781  if(second==NULL)
     782    {
     783      printf("Erreur typage!\n");
     784      return;
     785    }
     786
     787  sadj *a1;
     788  int i = premier->fv;
     789  int j = second->lv;
     790  if(L2[j]<=i && j<= R2[i])
     791    // ah ah les deux premiers fils forment un module
     792    // on recherche si adjacents
     793    {
     794      N->type = PARALLELE; // sera modifie si adjacence
     795      a1=S[i]->adj;
     796      while(a1!=NULL)
     797        {
     798          if(a1->pointe->place == S[j]->place)
     799            // adjacence trouvee !
     800            {
     801              N->type = SERIE;
     802              break;
     803            }
     804          a1=a1->suiv;
     805        }
     806    }
     807  else
     808    N->type = PREMIER;
     809
     810    fils *ffils;
     811    ffils = N->fils;
     812    do {
     813      typenoeud( ffils -> pointe, L2, R2, S);
     814      ffils = ffils->suiv;
     815    }
     816    while (ffils != NULL);
     817}
     818
     819
    819820noeud *algo2(graphe G, sommet **S)
    820821{
    821822/* algorithme de decomposition modulaire, deuxieme passe
    … …  
    824825*/
    825826    /* debug: S n'est utilise que pour mettre vrainom a jour */
    826827    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 */
    830828
    831     int *fermantes;                /* idem fermantes[i]=2 ssi i)))i+1
    832                                    fermante [n-1]=2 ssi n))) */
    833829    int *ps;                        /* ps[i]=premier separateur de (i,i+1) */
    834830    int *ds;
    835831
    … …  
    837833
    838834    sadj *a1, *a2;                /* parcours de liste d'adjacence */
    839835
    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 */
    845836
    846     int *adjii;                 /* adjii[i]=1 ssi S[i] et S[i+1] sont */
    847                                    /*                             adjacents */
     837    int *R2, *L2; // le generateur des intervalle-modules produit Capelle-like en phase 2
     838    int *R, *L; // le  generateur canonique des intervalle-modules
     839    int *SupR, *SupL; // le  support pour passer de (R2,L2) a (R,L)
     840
     841    int *pile;  // pile utilisee pour calcule generateur
     842    int sommet; // sommet de ladite pile. Indique la derniere place PLEINE
     843
     844    int *tab1,*tab2; // utilise pour le tri lineaire avant algo 6
     845    int *fortg, *fortd; // bornes de mes modules forts
     846    int nb_forts=0; // nombre de modules forts
     847
     848    int *newfortg, *newfortd; // tableau de modules forts utilises temporairement pour bucket sort
     849
    848850    /*PROPHASE : initialisations */
    849851    n=G.n;
    850     ouvrantes = (int *) fabmalloc(n * sizeof(int));
    851     fermantes = (int *) fabmalloc(n * sizeof(int));
    852852    ps = (int *) fabmalloc(n * sizeof(int));
    853853    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
     854    R2 = (int *)fabmalloc(n * sizeof(int));
     855    L2 = (int *)fabmalloc(n * sizeof(int));
     856    SupR = (int *)fabmalloc(n * sizeof(int));
     857    SupL = (int *)fabmalloc(n * sizeof(int));
     858    R = (int *)fabmalloc(n * sizeof(int));
     859    L = (int *)fabmalloc(n * sizeof(int));
     860    pile = (int *)fabmalloc(n * sizeof(int));
     861    tab1 = (int *)fabmalloc(4 * n * sizeof(int));
     862    tab2 = (int *)fabmalloc(2 * n * sizeof(int));
     863    fortg=fabmalloc(2 * n * sizeof(int));
     864    fortd=fabmalloc(2 * n * sizeof(int));
     865    newfortg= (int *)fabmalloc(2 * n * sizeof(int));
     866    newfortd=(int *) fabmalloc(2 * n * sizeof(int));
    875867
    876868    /* PREMIERE PASSE
    877869       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 */
     870       complexite: O(m) */
     871    // utilise ps ds
    884872
    885873    for (i = 0; i < n - 1; i++) {
    886874      /*recherche de ps(i,i+1) */
    … …  
    899887          ||((a2==NULL) && (a1->pointe->place >= i))
    900888          ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place >= i) && (a2->pointe->place >= i)))
    901889        //pas de separateur
    902         ps[i]=i+1;
     890        ps[i]=-1;
    903891      else
    904892        {
    905893          if((a1==NULL) || (a1->pointe->place >= i))
    … …  
    915903              else
    916904                ps[i]=min(a1->pointe->place , a2->pointe->place);
    917905            }
    918           ouvrantes[ps[i]]++;        /* marque la fracture gauche, if any */
    919           fermantes[i]++;
    920906        }
    921907      if (DEBUG)
    922         printf("ps(%i,%i)=%i\n", i , i+1, ps[i]);
    923 
     908        printf("ps(%i,%i)=%i ps(%i,%i)=%i\n", i , i+1, ps[i],
     909               S[i]->nom, S[i + 1]->nom, ps[i]==-1?-1:S[ps[i]]->nom );
    924910        /*recherche de ds(i,i+1)
    925911          plus penible encore!*/
    926912      a1=S[i]->adj;
    … …  
    942928          ||((a2==NULL) && (a1->pointe->place <= i+1))
    943929          ||((a1!=NULL) && (a2!=NULL) && (a1->pointe->place <= i+1) && (a2->pointe->place <= i+1)))
    944930        //pas de separateur
    945         ds[i]=i+1;
     931        ds[i]=-1;
    946932      else
    947933        {
    948934          if((a1==NULL) || (a1->pointe->place <= i+1))
    … …  
    961947
    962948
    963949          //ds[i] = j;
    964         ouvrantes[i + 1]++;        /* marque la fracture gauche, if any */
    965         fermantes[ds[i]]++;        /* attention aux decalages d'indices */
    966950        }
    967951      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);
     952        printf("ds(%i,%i)=%i ds(%i,%i)=%i\n", i,i+1,ds[i],
     953               S[i]->nom, S[i + 1]->nom, ds[i]==-1?-1:S[ds[i]]->nom );
    970954    }
    971955
    972     /*DEUXIEME PASSE: construction du pseudocoarbre */
     956    /* DEUXIEMME PASSE
     957       Calcule le generateur des interval-modules comme algo 9 du papier DAM
     958    */
     959    int v;
     960    // utilise ps ds L2 R2 pile
    973961
    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("(");
     962    // algo 9
     963    // Attention s[v] dans l'algo est  ds[v-1] !!
     964    sommet = -1;
     965
     966    for(v = n-1; v>=0; v--)
     967      if(ds[v-1] != -1)
     968        {
     969          L2[v]=v;
     970          while( pile[sommet] < ds[v-1])
     971            L2[pile[sommet--]]=v;
    985972        }
    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);
     973      else
     974        pile[++sommet]=v;
    997975
    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(")");
     976    while(sommet>=0)
     977      L2[pile[sommet--]]=0;
     978
     979    sommet = -1;
     980    // algo 9 symétrique
     981    //re-attention là c'est ps[v]
     982    for(v = 0 ; v<=n-1; v++)
     983      if(ps[v] != -1)
     984        {
     985          R2[v]=v;
     986          while(  pile[sommet] > ps[v])
     987            R2[pile[sommet--]]=v;
    1025988        }
    1026     }
     989      else
     990        pile[++sommet]=v;
     991
     992    while(sommet>=0)
     993      R2[pile[sommet--]]=n-1;
     994
    1027995    if (DEBUG)
     996      {    printf("L2 = [");
     997        for(v = 0;v<n;v++)
     998          printf("%i ",L2[v]);
     999        printf("]\nR2 = [");
     1000        for(v = 0;v<n;v++)
     1001          printf("%i ",R2[v]);
     1002        printf("]\n");
     1003        for(i=0;i<n;i++)
     1004          for(j=i+1;j<n;j++)
     1005            if(L2[j]<=i && j<= R2[i])
     1006              printf("[%i-%i] module\n",i,j);
     1007
     1008      }
     1009
     1010
     1011
     1012    /* TROISIEME PASSE
     1013       generateur canonique. Algos 3 et 5 du papier de DAM
     1014    */
     1015
     1016    // Algo 3, R
     1017    pile[0]=0;
     1018    sommet =0;
     1019    for(i=1;i<n;i++)
     1020      {
     1021        while(R2[pile[sommet]] < i)
     1022          sommet --;
     1023        SupR[i]=pile[sommet];
     1024        pile[++sommet]=i;
     1025      }
     1026
     1027    // Algo 3, L
     1028
     1029    pile[0]=n-1;
     1030    sommet =0;
     1031    for(i=n-2;i>=0;i--)
     1032      {
     1033        while(L2[pile[sommet]] > i)
     1034          sommet --;
     1035        SupL[i]=pile[sommet];
     1036        pile[++sommet]=i;
     1037      }
     1038
     1039    if (DEBUG)
     1040      {
     1041        printf("SupL = [");
     1042        for(v = 0;v<n;v++)
     1043          printf("%i ",SupL[v]);
     1044        printf("]\nSupR = [");
     1045        for(v = 0;v<n;v++)
     1046          printf("%i ",SupR[v]);
     1047        printf("]\n");
     1048      }
     1049
     1050    // Algo 5, R
     1051    int k;
     1052    R[0]=n-1;
     1053    for(k=1;k<n;k++)
     1054      R[k]=k;
     1055    for(k=n-1;k>=1;k--)
     1056      if(L2[R[k]] <= SupR[k] && R[k] <= R2[SupR[k]])
     1057        R[SupR[k]] = max(R[k], R[SupR[k]]);
     1058
     1059    // Algo 5, L
     1060    L[n-1]=0;
     1061    for(k=0;k<n-1;k++)
     1062      L[k]=k;
     1063    for(k=0;k<n-1;k++)
     1064      if(L2[SupL[k]] <= L[k] && SupL[k] <= R2[L[k]])
     1065        L[SupL[k]] = min(L[k], L[SupL[k]]);
     1066
     1067    if (DEBUG)
     1068      {    printf("L = [");
     1069        for(v = 0;v<n;v++)
     1070          printf("%i ",L[v]);
     1071        printf("]\nR = [");
     1072        for(v = 0;v<n;v++)
     1073          printf("%i ",R[v]);
     1074        printf("]\n");
     1075        for(i=0;i<n;i++)
     1076          for(j=i+1;j<n;j++)
     1077            if(L[j]<=i && j<= R[i])
     1078              printf("[%i-%i] module\n",i,j);
     1079
     1080      }
     1081
     1082    /* QUATRIEME PASSE
     1083       strong modules. Algo 6 du papier de DAM
     1084    */
     1085    //1 remplit le tableau des bornes
     1086    // astuce : on met 2a pour a borne gauche et 2a+1 pour a borne droite. Puis tri lineaire
     1087
     1088    for(i = 0; i<n; i++)
     1089      {
     1090        // intervalle (L[i]..i)
     1091        tab1[4*i] = 2*L[i];
     1092        tab1[4*i+1] = 2*i+1;
     1093        // intervalle (i..R[i])
     1094        tab1[4*i+2] = 2*i;
     1095        tab1[4*i+3] = 2*R[i]+1;
     1096      }
     1097    // tableau des fréquences
     1098
     1099    for( i=0;i<2*n;i++)
     1100      tab2[i]=0;
     1101    for( i=0;i<4*n;i++)
     1102      tab2[tab1[i]]++;
     1103    // tri
     1104    j = 0;
     1105    for( i=0;i<2*n;i++)
     1106      while(tab2[i]-->0)
     1107        {
     1108          tab1[j++] = i;
     1109        }
     1110
     1111    // Algo 6 du papier de DAM
     1112    // j'utilise maintenant tab2 comme stack de sommet : sommet
     1113    sommet = -1;
     1114
     1115
     1116    for( i=0;i<4*n;i++)
     1117      {
     1118        int ai,s; //meme noms de variables que papier
     1119        ai=tab1[i];
     1120        if((ai%2)==0)
     1121          tab2[++sommet] = ai;
     1122        else
     1123          {
     1124            ai/=2;
     1125            s = tab2[sommet--]/2;
     1126
     1127            if(nb_forts == 0 || (fortg[nb_forts - 1] != s) || (fortd[nb_forts - 1] != ai))
     1128              {
     1129                fortg[nb_forts] = s;
     1130                fortd[nb_forts] = ai;
     1131                nb_forts++;
     1132              }
     1133          }
     1134      }
     1135    if(DEBUG)
     1136      {
     1137        printf("Modules forts :\n");
     1138        for( i=0;i<nb_forts;i++)
     1139          printf("[%i %i] ", fortg[i], fortd[i]);
    10281140        printf("\n");
     1141      }
    10291142
    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);
    1037     }
     1143    // TRI des modules forts par bornes droites decroissantes
     1144    for( i=0;i<n;i++)
     1145      tab1[i]=0;
     1146    for( i=0;i<nb_forts;i++)
     1147      tab1[fortd[i]]++;
     1148    for( i=1;i<n;i++)
     1149      tab1[i]+=tab1[i-1];
    10381150
    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 */
     1151    // tab1 = tableau des frequences cumulées
     1152    for( i=0;i<nb_forts;i++)
     1153      {
     1154        int pos = (nb_forts-1)-(tab1[fortd[i]]-1);
     1155        newfortg[ pos ] = fortg[i];
     1156        newfortd[ pos ] = fortd[i];
     1157        tab1[fortd[i]]--;
     1158      }
     1159    if(DEBUG)
     1160      {
     1161        printf("Modules forts :\n");
     1162        for( i=0;i<nb_forts;i++)
     1163          printf("[%i %i] ", newfortg[i], newfortd[i]);
     1164        printf("\n");
     1165      }
    10461166
    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 */
     1167   // now tri par bornes gauches croissantes
     1168    for( i=0;i<n;i++)
     1169      tab1[i]=0;
     1170    for( i=0;i<nb_forts;i++)
     1171      tab1[newfortg[i]]++;
     1172    for( i=1;i<n;i++)
     1173      tab1[i]+=tab1[i-1];
    10611174
    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]);
     1175    // tab1 = tableau des frequences cumulées
     1176    for( i = nb_forts-1 ; i>=0 ; i--)// parcours de droite a gauche car on retourne tout
     1177      {
     1178        int pos = tab1[newfortg[i]]-1;
     1179        // printf("[%i %i] en pos %i\n",newfortg[i],newfortd[i],pos);
     1180        fortg[ pos ] = newfortg[i];
     1181        fortd[ pos ] = newfortd[i];
     1182        tab1[newfortg[i]]--;
     1183      }
     1184    if(DEBUG)
     1185      {
     1186        printf("Modules forts :\n");
     1187        for( i=0;i<nb_forts;i++)
     1188          printf("[%i %i] ", fortg[i], fortd[i]);
     1189        printf("\n");
     1190      }
    10791191
    1080             nextfils = nextfils->suiv;
    1081         }
    1082         while (nextfils != NULL);
     1192    // Finally algo 7 du papier : construit l'arbre
     1193    noeud* racine = nouvnoeud( UNKN, NULL, -1, 0, n-1);
    10831194
     1195    noeud *F = racine;
     1196    noeud *Nouv;
     1197    for(k=1; k< nb_forts; )
     1198      {
     1199        if(F->fv <= fortg[k] && F->lv >= fortd[k] )
     1200          {
     1201            if(fortg[k]==fortd[k]) // on a une feuille
     1202              Nouv = nouvnoeud(FEUILLE, F, S[fortg[k]]->nom, fortg[k], fortd[k]);
     1203            else
     1204              Nouv = nouvnoeud(UNKN, F, -1, fortg[k], fortd[k]);
     1205            ajoutfils(F, Nouv);
     1206            // printf("nouv [%i %i] de pere [%i %i]\n",
     1207            // fortg[k], fortd[k], F->fv, F->lv);
     1208            F = Nouv;
     1209            k++;
    10841210
    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);
     1211          }
     1212        else
     1213          {
     1214            // printf("remonte car [%i %i] pas fils de [%i %i]\n",
     1215            // fortg[k], fortd[k], F->fv, F->lv);
     1216            F = F->pere;
     1217          }
     1218      }
     1219    if(DEBUG)
     1220      {
     1221        printf("arbre avant typage :\n");
     1222        printarbre(racine);
     1223      }
    10881224
    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     }
     1225    /* CINQUIEME PASSE : typage */
     1226    // Methode : on teste grace au generateur la modularite des deux premiers fils.
     1227    // Ensuite si module est serie si adjacence parallele sinon
     1228    typenoeud(racine, L2, R2, S);
    11031229
    1104     if (DEBUG) {
    1105         printf("Arbre apres la troisieme passe:\n");
    1106         printnoeud(racine, 0);
    1107     }
     1230    if(DEBUG)
     1231      {
     1232        printf("\n arbre final :\n");
     1233        printarbre(racine);
     1234      }
     1235    /* nettoie un peu */
     1236    free(ps);
     1237    free(ds);
     1238    free(R2);
     1239    free(L2);
     1240    free(SupR);
     1241    free(SupL);
     1242    free(R);
     1243    free(L);
     1244    free(pile);
     1245    free(tab1);
     1246    free(tab2);
     1247    free(fortg);
     1248    free(fortd);
     1249    free(newfortg);
     1250    free(newfortd);
    11081251
    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!) */
    11121252
    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     }
    1129 
    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 */
    1146 
    1147         scanne = pileinterne[indicepileinterne];
    1148         if (scanne->type != MODULE)
    1149             continue;                /* je traite que les modules */
    1150 
    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     }
    12141253    return racine;
    12151254}
    12161255
    … …  
    12521291     /* affiche la permutation factorisante */
    12531292{
    12541293  int i;
    1255   printf(  "Place (nouvelle num) ");
     1294  printf("Attention dans debug maintenant usage des NOUVEAUX noms\n");
     1295  printf(  "Place (nouveaux noms, nouvelle numerotation) ");
    12561296  for(i=0;i<n;i++)
    12571297    printf("%3i ",S[i]->place);
    1258   printf("\nNom (ancienne num) : ");
     1298   printf("\nNom (ancienne num) :                        ");
    12591299  for(i=0;i<n;i++)
    12601300    printf("%3i ",S[i]->nom);
    12611301  printf("\n");
    12621302}
    12631303
     1304void nettoie(sommet **S, int n)
     1305{
     1306  // Detruit la copie locale du graphe
     1307  // normalement, apres, tous les mallocs sont
     1308  // annules sauf ceux de l'arbre (le resultat)
     1309  int i;
     1310  sommet * scourant;
     1311  sadj *a,*atmp;
     1312  for (i = 0; i < n; i++) {
     1313        /* initialisation des sommets */
     1314        /* notre bebe est le sommet i dans M */
     1315    scourant = S[i];
     1316    a=scourant->adj;
     1317      while(a!=NULL)
     1318        {
     1319          atmp=a->suiv;
     1320          free(a);
     1321          a=atmp;
     1322        }
     1323      free(scourant->classe); // efface toutes les classes car sont singletons
     1324      free(scourant);
     1325  }
     1326  free(S);
     1327}
    12641328
    12651329/* la fonction principale; qui fait pas grand'chose....*/
    12661330noeud *decomposition_modulaire(graphe G)
    … …  
    12821346        PrintS2(S,G.n);
    12831347      }
    12841348    Racine = algo2(G, S);       /* deuxieme partie: calcul de l'arbre */
     1349
     1350    nettoie(S, G.n);
     1351
    12851352    return Racine;
    12861353}

  • 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)