/* This file is part of GNU Pies.
   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
   the Free Software Foundation; either version 3, or (at your option)
   any later version.

   GNU Pies is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with GNU Pies.  If not, see <http://www.gnu.org/licenses/>. */

#include "pies.h"

#define INFINITY INT_MAX

struct pies_depmap
{
  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 */
};

size_t
depmap_dim (struct pies_depmap *dmap)
{
  return dmap->dim;
}

pies_depmap_t
depmap_alloc (size_t count)
{
  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 src)
{
  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;
}

void
depmap_free (pies_depmap_t dmap)
{
  if (dmap)
    {
      free (dmap->dist);
      free (dmap->next);
      free (dmap);
    }
}

void
depmap_set (pies_depmap_t dmap, size_t i, size_t j)
{
  dmap->dist[i][j] = 1;
  dmap->closed = 0;
}

int
depmap_isset (pies_depmap_t dmap, size_t i, size_t j)
{
  return dmap->dist[i][j] != INFINITY;
}

/* Remove nth row and column from the map */
void
depmap_remove (pies_depmap_t dmap, size_t n)
{
  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)
{
  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;
  size_t coord[2];
};

size_t
depmap_next (pies_depmap_t dmap, pies_depmap_pos_t pos)
{
  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->dim ?
	     (size_t) -1 : pos->coord[pos->dir];
}

size_t
depmap_first (pies_depmap_t dmap, enum pies_depmap_direction dir,
	      size_t coord, pies_depmap_pos_t *ppos)
{
  pies_depmap_pos_t pos = grecs_malloc (sizeof *pos);
  *ppos = pos;
  pos->dir = dir;
  pos->coord[!pos->dir] = coord;
  pos->coord[pos->dir] = -1;
  return depmap_next (dmap, pos);
}

void
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;
}

void
depmap_path_free (struct depmap_path *path)
{
  while (path)
    {
      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;
    }
}

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;
	    }
	}
    }
  depmap_free (dp);
  return head;
}