1 files changed, 292 insertions, 114 deletions

diff --git a/src/depmap.c b/src/depmap.c
index fd5521f..a74a7f0 100644
--- a/src/depmap.c
+++ b/src/depmap.c
@@ -1,5 +1,5 @@
 /* This file is part of GNU Pies.
-   Copyright (C) 2008-2020 Sergey Poznyakoff
+   Copyright (C) 2008-2023 Sergey Poznyakoff
 
    GNU Pies is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -16,135 +16,179 @@
 
 #include "pies.h"
 
-#ifndef CHAR_BIT
-# define CHAR_BIT 8
-#endif
-#define BITS_PER_WORD   (sizeof(unsigned)*CHAR_BIT)
-
-#define WORDSIZE(n)     (((n) + BITS_PER_WORD - 1) / BITS_PER_WORD)
-#define SETBIT(x, i)    ((x)[(i)/BITS_PER_WORD] |= (1<<((i) % BITS_PER_WORD)))
-#define CLRBIT(x, i)    ((x)[(i)/BITS_PER_WORD] &= ~(1<<((i) % BITS_PER_WORD)))
-#define BITISSET(x, i)  (((x)[(i)/BITS_PER_WORD] & (1<<((i) % BITS_PER_WORD))) != 0)
-
-void
-TC (unsigned *R, int n)
-{
-  register int rowsize;
-  register unsigned mask;
-  register unsigned *rowj;
-  register unsigned *rp;
-  register unsigned *rend;
-  register unsigned *ccol;
-
-  unsigned *relend;
-  unsigned *cword;
-  unsigned *rowi;
-
-  rowsize = WORDSIZE (n) * sizeof (unsigned);
-  relend = (unsigned *) ((char *) R + (n * rowsize));
-
-  cword = R;
-  mask = 1;
-  rowi = R;
-  while (rowi < relend)
-    {
-      ccol = cword;
-      rowj = R;
-
-      while (rowj < relend)
-	{
-	  if (*ccol & mask)
-	    {
-	      rp = rowi;
-	      rend = (unsigned *) ((char *) rowj + rowsize);
-
-	      while (rowj < rend)
-		*rowj++ |= *rp++;
-	    }
-	  else
-	    {
-	      rowj = (unsigned *) ((char *) rowj + rowsize);
-	    }
-
-	  ccol = (unsigned *) ((char *) ccol + rowsize);
-	}
-
-      mask <<= 1;
-      if (mask == 0)
-	{
-	  mask = 1;
-	  cword++;
-	}
-      rowi = (unsigned *) ((char *) rowi + rowsize);
-    }
-}
+#define INFINITY INT_MAX
 
 struct pies_depmap
 {
-  size_t nrows;         /* Number of rows (== number of mapped elements) */
-  size_t rowlen;        /* Length of a row in words */
-  unsigned r[1];        /* Data rows */
+  size_t dim;     /* Number of vertices */
+  int **dist;     /* Distance between two vertices */
+  int **next;     /* Next vertex on the shortest path from i to j */
+  int closed;     /* True if the transitive closure was computed */
+  int map[1];     /* Storage array for dist and next */
 };
 
-pies_depmap_t
-depmap_alloc (size_t count)
+size_t
+depmap_dim (struct pies_depmap *dmap)
 {
-  size_t size = (count + BITS_PER_WORD - 1) / BITS_PER_WORD;
-  pies_depmap_t dmap = grecs_zalloc (sizeof (*dmap) - 1
-				     + count * size * sizeof (unsigned));
-  dmap->nrows  = count;
-  dmap->rowlen = size;
-  return dmap;
+  return dmap->dim;
 }
 
-/* Return size of a depmap row in bytes */
-static inline size_t
-depmap_row_size (pies_depmap_t dpm)
+pies_depmap_t
+depmap_alloc (size_t count)
 {
-  return dpm->rowlen * sizeof (unsigned);
-}
+  struct pies_depmap *dmap;
+  size_t size;
+  size_t i, j;
+  int *ip;
+
+  if (count < 0)
+    abort ();
+  /* size = sizeof (*dmap) + 2 * count * count * sizeof (int) */
+  if (count > 0)
+    {
+      if (SIZE_MAX / count < count)
+	grecs_alloc_die ();
+      size = count * count;
+      if (SIZE_MAX / 2 < size)
+	grecs_alloc_die ();
+      size += size;
+      if (SIZE_MAX / sizeof (int) < size)
+	grecs_alloc_die ();
+      size *= sizeof (int);
+      if (SIZE_MAX - size < sizeof (*dmap))
+	grecs_alloc_die ();
+      size += sizeof (*dmap);
+    }
+  else
+    size = sizeof (*dmap);
+  
+  dmap = grecs_zalloc (size);
+  dmap->dim = count;
+  dmap->closed = 0;
+  dmap->dist = grecs_calloc (sizeof (dmap->dist[0]), count);
+  dmap->next = grecs_calloc (sizeof (dmap->next[0]), count);
+
+  ip = dmap->map;
+  for (i = 0; i < count; i++)
+    {
+      dmap->dist[i] = ip;
+      for (j = 0; j < count; j++)
+	*ip++ = INFINITY;
+    }
+	
+  for (i = 0; i < count; i++)
+    {
+      dmap->next[i] = ip;
+      for (j = 0; j < count; j++)
+	*ip++ = -1;
+    }
+  
+  return dmap;
+}      
 
 pies_depmap_t
-depmap_copy (pies_depmap_t dpm)
+depmap_copy (pies_depmap_t src)
 {
-  pies_depmap_t copy = depmap_alloc (dpm->nrows);
-  memcpy (copy->r, dpm->r, dpm->nrows * depmap_row_size (dpm));
-  return copy;
+  pies_depmap_t dst = depmap_alloc (src->dim);
+  memcpy (dst->map, src->map, 2 * src->dim * src->dim * sizeof (src->map[0]));
+  dst->closed = src->closed;
+  return dst;
 }
 
-static unsigned *
-depmap_rowptr (pies_depmap_t dmap, size_t row)
+void
+depmap_free (pies_depmap_t dmap)
 {
-  return dmap->r + dmap->rowlen * row;
+  if (dmap)
+    {
+      free (dmap->dist);
+      free (dmap->next);
+      free (dmap);
+    }
 }
 
 void
-depmap_set (pies_depmap_t dmap, size_t row, size_t col)
+depmap_set (pies_depmap_t dmap, size_t i, size_t j)
 {
-  unsigned *rptr = depmap_rowptr (dmap, row);
-  SETBIT (rptr, col);
+  dmap->dist[i][j] = 1;
+  dmap->closed = 0;
 }
 
-void
-depmap_clear (pies_depmap_t dmap, size_t row, size_t col)
+int
+depmap_isset (pies_depmap_t dmap, size_t i, size_t j)
 {
-  unsigned *rptr = depmap_rowptr (dmap, row);
-  CLRBIT (rptr, col);
+  return dmap->dist[i][j] != INFINITY;
 }
 
-int
-depmap_isset (pies_depmap_t dmap, size_t row, size_t col)
+/* Remove nth row and column from the map */
+void
+depmap_remove (pies_depmap_t dmap, size_t n)
 {
-  unsigned *rptr = depmap_rowptr (dmap, row);
-  return BITISSET (rptr, col);
+  if (n < dmap->dim - 1)
+    {
+      size_t i;
+      
+      /* Remove nth row */
+      memmove (&dmap->dist[n], &dmap->dist[n + 1],
+	       (dmap->dim - n - 1) * sizeof (dmap->dist[0]));
+      memmove (&dmap->next[n], &dmap->next[n + 1],
+	       (dmap->dim - n - 1) * sizeof (dmap->next[0]));
+
+      /* Remove nth column */
+      for (i = 0; i < dmap->dim; i++)
+	{
+	  memmove (&dmap->dist[i][n], &dmap->dist[i][n+1],
+		   (dmap->dim - n - 1) * sizeof (dmap->dist[0][0]));
+	  memmove (&dmap->next[i][n], &dmap->next[i][n + 1],
+		   (dmap->dim - n - 1) * sizeof (dmap->next[0][0]));
+	}
+    }
+  dmap->dim--;
 }
 
 void
 depmap_tc (pies_depmap_t dmap)
 {
-  TC (dmap->r, dmap->nrows);
-}
+  size_t i, j, k;
+  
+  if (dmap->closed) return;
+  /*
+   * Initialize next vertex map and reset all distances between
+   * adjacent vertices to 1.
+   */
+  for (i = 0; i < dmap->dim; i++)
+    for (j = 0; j < dmap->dim; j++)
+      {
+	if (depmap_isset (dmap, i, j))
+	  {
+	    dmap->dist[i][j] = 1;
+	    dmap->next[i][j] = j;
+	  }
+	else
+	  {
+	    dmap->next[i][j] = -1;
+	  }
+      }
 
+  /* Modified Floyd-Warshall algorithm */
+  for (k = 0; k < dmap->dim; k++)
+    for (i = 0; i < dmap->dim; i++)
+      for (j = 0; j < dmap->dim; j++)
+	{
+	  if (depmap_isset (dmap, i, k) && depmap_isset (dmap, k, j))
+	    {
+	      int t = dmap->dist[i][k] + dmap->dist[k][j];
+	      if (dmap->dist[i][j] > t)
+		{
+		  dmap->dist[i][j] = t;
+		  dmap->next[i][j] = dmap->next[i][k];
+		}
+	    }
+	}
+
+  /* Mark completion */
+  dmap->closed = 1;
+}
+  
 struct pies_depmap_pos
 {
   enum pies_depmap_direction dir;
@@ -154,12 +198,12 @@ struct pies_depmap_pos
 size_t
 depmap_next (pies_depmap_t dmap, pies_depmap_pos_t pos)
 {
-  for (pos->coord[pos->dir]++; pos->coord[pos->dir] < dmap->nrows;
+  for (pos->coord[pos->dir]++; pos->coord[pos->dir] < dmap->dim;
        pos->coord[pos->dir]++)
     if (depmap_isset (dmap, pos->coord[0], pos->coord[1]))
       break;
 
-  return pos->coord[pos->dir] == dmap->nrows ?
+  return pos->coord[pos->dir] == dmap->dim ?
 	     (size_t) -1 : pos->coord[pos->dir];
 }
 
@@ -180,23 +224,157 @@ depmap_end (pies_depmap_pos_t pos)
 {
   grecs_free (pos);
 }
+
+static struct depmap_path *
+depmap_path_alloc (void)
+{
+  struct depmap_path *path = grecs_malloc (sizeof *path);
+  path->len = 0;
+  path->head = path->tail = NULL;
+  path->next = NULL;
+  return path;
+}
 
-/* Remove nth row and column from the map */
 void
-depmap_remove (pies_depmap_t dmap, size_t n)
+depmap_path_free (struct depmap_path *path)
 {
-  if (n < dmap->nrows - 1)
+  while (path)
     {
-      size_t i, j;
+      struct depmap_path *next_path = path->next;
+      struct depmap_path_elem *p = path->head;
+      while (p)
+	{
+	  struct depmap_path_elem *next = p->next;
+	  free (p);
+	  p = next;
+	}
+      free (path);
+      path = next_path;
+    }
+}
 
-      /* Remove nth row */
-      memmove (depmap_rowptr (dmap, n), depmap_rowptr (dmap, n + 1),
-	       (dmap->nrows - n - 1) * depmap_row_size (dmap));
-      /* Renumber all columns j >= n: j = j - 1 */
-      for (i = 0; i < dmap->nrows; i++)
-	for (j = n; j < dmap->nrows - 1; j++)
-	  (depmap_isset (dmap, i, j + 1) ? depmap_set : depmap_clear)
-	    (dmap, i, j);
+static void
+depmap_path_append (struct depmap_path *path, int idx)
+{
+  struct depmap_path_elem *p = grecs_malloc (sizeof (*p));
+  p->idx = idx;
+  p->next = NULL;
+  if (path->tail)
+    path->tail->next = p;
+  else
+    path->head = p;
+  path->tail = p;
+  path->len++;
+}
+  
+/* Returns the path from vertex i to j. */
+static struct depmap_path *
+depmap_path (pies_depmap_t dmap, int i, int j)
+{
+  struct depmap_path *dp;
+  if (dmap->closed && depmap_isset (dmap, i, j))
+    {
+      dp = depmap_path_alloc ();
+      depmap_path_append (dp, i);
+      while (i != j)
+	{
+	  i = dmap->next[i][j];
+	  depmap_path_append (dp, i);
+	}
+    }
+  else
+    dp = NULL;
+  return dp;
+}
+
+/*
+ * If the vertex i is a part of a cyclic path, return the path adjusted
+ * so that it starts at the vertex with the lowest number among the vertices
+ * traversed by that path.
+ */
+static struct depmap_path *
+depmap_normalized_cyclic_path (pies_depmap_t dmap, int i)
+{
+  struct depmap_path *dp;
+  if (dmap->closed && depmap_isset (dmap, i, i))
+    {
+      struct depmap_path_elem *start, *p;
+      
+      dp = depmap_path (dmap, dmap->next[i][i], i);
+      start = dp->head;
+      for (p = start->next; p; p = p->next)
+	{
+	  if (p->idx < start->idx)
+	    start = p;
+	}
+      if (start != dp->head)
+	{
+	  dp->tail->next = dp->head;
+	  for (p = dp->head; p->next != start; p = p->next)
+	    ;
+	  dp->tail = p;
+	  dp->tail->next = NULL;
+	  dp->head = start;
+	}
+    }
+  else
+    dp = NULL;
+  return dp;
+}
+
+/* Return true (1) if normalized paths A and B are equal. */
+static int
+depmap_path_eq (struct depmap_path *a, struct depmap_path *b)
+{
+  struct depmap_path_elem *ap, *bp;
+  ap = a->head;
+  bp = b->head;
+  while (ap)
+    {
+      if (bp == NULL || ap->idx != bp->idx)
+	return 0;
+      ap = ap->next;
+      bp = bp->next;
+    }
+  return bp == NULL;
+}
+
+static int
+depmap_path_find (struct depmap_path *a, struct depmap_path *b)
+{
+  while (a)
+    {
+      if (depmap_path_eq (a, b))
+	return 1;
+      a = a->next;
+    }
+  return 0;
+}
+
+struct depmap_path *
+depmap_cycle_detect (pies_depmap_t dmap)
+{
+  pies_depmap_t dp = depmap_copy (dmap);
+  size_t i;
+  struct depmap_path *head = NULL, *tail = NULL;
+  
+  depmap_tc (dp);
+  for (i = 0; i < dmap->dim; i++)
+    {
+      if (depmap_isset (dp, i, i))
+	{
+	  struct depmap_path *path = depmap_normalized_cyclic_path (dp, i);
+	  if (!head)
+	    head = tail = path;
+	  else if (depmap_path_find (head, path))
+	    depmap_path_free (path);
+	  else
+	    {
+	      tail->next = path;
+	      tail = path;
+	    }
+	}
     }
-  dmap->nrows--;
+  depmap_free (dp);
+  return head;
 }