#include <stdio.h>
#include "util.h"
#include "cudd.h"

/* Towers of Hanoi problem using BDDs */

#define SIZE  10   /* number of disks */
/* number of towers = 3 */

DdManager * manager;   /* BDD Manager */
DdNode * v1[SIZE][3];  /* the current state */
DdNode * v2[SIZE][3];  /* the next state */

/* v[i][j] = 1 means that disk i is in tower j */


DdNode * make_a_move(int i,int j,int k);
DdNode * compute_image(DdNode * R, DdNode * T, DdNode * cube);

int main(int argv, char ** argc) {

  DdNode * I;  /* the initial state BDD */
  DdNode * E;  /* the final state BDD */
  DdNode * T;  /* the transition relation */
  DdNode * R;  
  DdNode * image;
  int i,j,found,num_steps;
  DdNode * tmp1, *tmp2;
  DdNode * cube;

  /* Initialize the BDD manager with default options */
  manager = Cudd_Init(0,0,CUDD_UNIQUE_SLOTS,CUDD_CACHE_SLOTS,0);

  /* enable dynamic reordering */
  Cudd_AutodynEnable(manager,CUDD_REORDER_SAME);
  
  /* the cuurent and next state variables are interleaved, usually a 
     good ordering to start with */
  for(i=0;i<SIZE;i++) {
    for(j=0;j<3;j++) {
      v1[i][j] = Cudd_bddNewVar(manager);
      v2[i][j] = Cudd_bddNewVar(manager);
    }
  }
  

  /* In the initial state, all the disks are in tower 0 */
  I = Cudd_ReadOne(manager);
  Cudd_Ref(I);
  for(i=0;i<SIZE;i++) {
    tmp1 = Cudd_bddAnd(manager,I,v1[i][0]);
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,I);
    I = tmp1;
  }
  for(i=0;i<SIZE;i++) {
    tmp1 = Cudd_bddAnd(manager,I,Cudd_Not(v1[i][1]));
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,I);
    I = tmp1;
  }
  for(i=0;i<SIZE;i++) {
    tmp1 = Cudd_bddAnd(manager,I,Cudd_Not(v1[i][2]));
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,I);
    I = tmp1;
  }

  /* In the final state, we want all the disks in tower 2 */
  E = Cudd_ReadOne(manager);
  Cudd_Ref(E);
  for(i=0;i<SIZE;i++) {
    tmp1 = Cudd_bddAnd(manager,E,v1[i][2]);
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,E);
    E = tmp1;
  }
  for(i=0;i<SIZE;i++) {
    tmp1 = Cudd_bddAnd(manager,E,Cudd_Not(v1[i][0]));
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,E);
    E = tmp1;
  }
  for(i=0;i<SIZE;i++) {
    tmp1 = Cudd_bddAnd(manager,E,Cudd_Not(v1[i][1]));
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,E);
    E = tmp1;
  }
  

  T = Cudd_ReadOne(manager);
  Cudd_Ref(T);
  T = Cudd_Not(T);


  /* the transition relation is the disjunction of all possible moves */
  for(i=0;i<SIZE;i++) {
    for(j=0;j<3;j++) {
      tmp1 = make_a_move(i,j,(j+1)%3);
      
      tmp2 = Cudd_bddOr(manager,T,tmp1);
      Cudd_Ref(tmp2);

      Cudd_RecursiveDeref(manager,tmp1);
      Cudd_RecursiveDeref(manager,T);
      
      T = tmp2;
      
      tmp1 = make_a_move(i,j,(j+2)%3);
      
      tmp2 = Cudd_bddOr(manager,T,tmp1);
      Cudd_Ref(tmp2);

      Cudd_RecursiveDeref(manager,tmp1);
      Cudd_RecursiveDeref(manager,T);
      
      T = tmp2;
      
    }
    
  }

  /* Cube of the current state variables, needed by bddAndAbstract */
  cube = Cudd_ReadOne(manager);
  Cudd_Ref(cube);
  for(i=0;i<SIZE;i++)
    for(j=0;j<3;j++) {
      tmp1 = Cudd_bddAnd(manager,cube,v1[i][j]);
      Cudd_Ref(tmp1);
      Cudd_RecursiveDeref(manager,cube);
      cube = tmp1;
    }
  

  num_steps = 0;
  found = 0;

  /* R contains the states reached so far, initialized to I */
  R = I;
  Cudd_Ref(R);

  
  /* Fixed point computation */
  while(1) {

    /* check if we reached the goal state */
    tmp1 = Cudd_bddOr(manager,Cudd_Not(E),R);
    Cudd_Ref(tmp1);
    if (tmp1 == Cudd_ReadOne(manager)) {
      found = 1;   /* goal reached */
      break;
    }
    Cudd_RecursiveDeref(manager,tmp1);

    /* compute the successors of the current state */
    image = compute_image(R,T,cube);
    num_steps ++;

    /* check if we reached a new state */
    /* Note: Fixed point check, easy with BDDs as they are canonical */ 
    tmp1 = Cudd_bddOr(manager,Cudd_Not(image),R);
    Cudd_Ref(tmp1);
    if (tmp1 == Cudd_ReadOne(manager))   
      break;  /* no new state reached */
    Cudd_RecursiveDeref(manager,tmp1);
    
    /* add the new states to the reached set */
    tmp1 = Cudd_bddOr(manager,image,R);
    Cudd_Ref(tmp1);
    
    Cudd_RecursiveDeref(manager,R);
    Cudd_RecursiveDeref(manager,image);
  
    R = tmp1;
  }
  
  Cudd_RecursiveDeref(manager,tmp1);
  Cudd_RecursiveDeref(manager,image);

  if (found) {
    printf("Goal Reached in %d steps\n",num_steps);
  }
  else
    printf("Goal not reached\n");


  /* Free BDDs and the manager */
  Cudd_RecursiveDeref(manager,cube);
  Cudd_RecursiveDeref(manager,I);
  Cudd_RecursiveDeref(manager,T);
  Cudd_RecursiveDeref(manager,E);
  Cudd_Quit(manager);

}



/* Returns a BDD representing move of disk i from tower j to tower k */

DdNode * make_a_move(int i,int j,int k)
{
  DdNode * result;
  int l,m;
  DdNode * tmp1;
  DdNode * tmp2;

  result = v1[i][j];  /* disk i is in tower j */
  Cudd_Ref(result);

  /* this loop enforces that there is no smaller disk in tower j or tower k */
  /* no smaller disk in tower j implies that the disk at hand is at the top */
  /* no smaller disk in tower k implies that we can move the disk to tower k */
  for(l=i-1;l>=0;l--) {
    tmp1 = Cudd_bddAnd(manager,result,Cudd_Not(v1[l][j]));
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,result);
    result = tmp1;
    
    tmp1 = Cudd_bddAnd(manager,result,Cudd_Not(v1[l][k]));
    Cudd_Ref(tmp1);
    Cudd_RecursiveDeref(manager,result);
    result = tmp1;
  }

  /* move the current disk to tower k */
  tmp1 = Cudd_bddAnd(manager,result,v2[i][k]);
  Cudd_Ref(tmp1);
  Cudd_RecursiveDeref(manager,result);
  result = tmp1;

  tmp1 = Cudd_bddAnd(manager,result,Cudd_Not(v2[i][(k+1)%3]));
  Cudd_Ref(tmp1);
  Cudd_RecursiveDeref(manager,result);
  result = tmp1;

  tmp1 = Cudd_bddAnd(manager,result,Cudd_Not(v2[i][(k+2)%3]));
  Cudd_Ref(tmp1);
  Cudd_RecursiveDeref(manager,result);
  result = tmp1;

  /* the other disks stay where they are */
  for(l=0;l<SIZE;l++) {
    if (l!=i) {
      for(m=0;m<3;m++) {
	
	tmp1 = Cudd_bddXnor(manager,v1[l][m],v2[l][m]);
	Cudd_Ref(tmp1);

	tmp2 = Cudd_bddAnd(manager,result,tmp1);
	Cudd_Ref(tmp2);
	Cudd_RecursiveDeref(manager,result);
	Cudd_RecursiveDeref(manager,tmp1);

	result = tmp2;
      }
    }
  }
  
  return result;
}

/* Given a set of states R and a transition relation T, returns the BDD for the
   states that can be reached in one step from R following T */

DdNode * compute_image(DdNode * R, DdNode * T, DdNode * cube)
{
  DdNode * result;
  int i;
  int j;
  DdNode * tmp1;
  
  /* the following Cudd function computes the conjunction of R and T; and quantifies out the variables
     in cube ( the current state variables ) */
  result = Cudd_bddAndAbstract(manager,R,T,cube); 
  Cudd_Ref(result);

  /* rename the next state variables to the current state variables */
  for(i=0;i<SIZE;i++)
    for(j=0;j<3;j++) {
      /* compose replaces v2[i][j] with v1[i][j] */
      tmp1 = Cudd_bddCompose(manager,result,v1[i][j],2*(i*3 + j) + 1); 
      Cudd_Ref(tmp1);
      Cudd_RecursiveDeref(manager,result);
      result = tmp1;
    }
  
  return result;
}



/* 

This example essentially illustrates reachability computation, which
is one of the core operations inside a model checker.  In model
checking, the goal state represents the error state, and the model
checker either:

1. Discovers that the error state is unreachable if the fixed point is
   achieved before hitting the error.

2. Discovers that the error is reachable. In that case it returns a
   counterexample, which is a sequence of steps from the initial state
   to the error state.

The above implementation does not produce a counterexample.

Exercise: Add counterexample generation to the above.

How: 

You will need to store pointers to all the image BDDs.  Starting from
the error state, compute a "reverse image", intersect with the
correspoding forward image, pick a state from the intersection and
repeat till you hit the initial state.

*/