/**CFile*********************************************************************** FileName [cuddSat.c] PackageName [cudd] Synopsis [Functions for the solution of satisfiability related problems.] Description [External procedures included in this file:
pr is positive
  the first failure is reported to the standard output.]
  SideEffects [None]
  SeeAlso     []
******************************************************************************/
int
Cudd_EqualSupNorm(
  DdManager * dd /* manager */,
  DdNode * f /* first ADD */,
  DdNode * g /* second ADD */,
  CUDD_VALUE_TYPE  tolerance /* maximum allowed difference */,
  int  pr /* verbosity level */)
{
    DdNode *fv, *fvn, *gv, *gvn, *r;
    unsigned int topf, topg;
    statLine(dd);
    /* Check terminal cases. */
    if (f == g) return(1);
    if (Cudd_IsConstant(f) && Cudd_IsConstant(g)) {
	if (ddEqualVal(cuddV(f),cuddV(g),tolerance)) {
	    return(1);
	} else {
	    if (pr>0) {
		(void) fprintf(dd->out,"Offending nodes:\n");
		(void) fprintf(dd->out,
			       "f: address = %p\t value = %40.30f\n",
			       (void *) f, cuddV(f));
		(void) fprintf(dd->out,
			       "g: address = %p\t value = %40.30f\n",
			       (void *) g, cuddV(g));
	    }
	    return(0);
	}
    }
    /* We only insert the result in the cache if the comparison is
    ** successful. Therefore, if we hit we return 1. */
    r = cuddCacheLookup2(dd,(DD_CTFP)Cudd_EqualSupNorm,f,g);
    if (r != NULL) {
	return(1);
    }
    /* Compute the cofactors and solve the recursive subproblems. */
    topf = cuddI(dd,f->index);
    topg = cuddI(dd,g->index);
    if (topf <= topg) {fv = cuddT(f); fvn = cuddE(f);} else {fv = fvn = f;}
    if (topg <= topf) {gv = cuddT(g); gvn = cuddE(g);} else {gv = gvn = g;}
    if (!Cudd_EqualSupNorm(dd,fv,gv,tolerance,pr)) return(0);
    if (!Cudd_EqualSupNorm(dd,fvn,gvn,tolerance,pr)) return(0);
    cuddCacheInsert2(dd,(DD_CTFP)Cudd_EqualSupNorm,f,g,DD_ONE(dd));
    return(1);
} /* end of Cudd_EqualSupNorm */
/**Function********************************************************************
  Synopsis    [Expands cube to a prime implicant of f.]
  Description [Expands cube to a prime implicant of f. Returns the prime
  if successful; NULL otherwise.  In particular, NULL is returned if cube
  is not a real cube or is not an implicant of f.]
  SideEffects [None]
  SeeAlso     [Cudd_bddMaximallyExpand]
******************************************************************************/
DdNode *
Cudd_bddMakePrime(
  DdManager *dd /* manager */,
  DdNode *cube /* cube to be expanded */,
  DdNode *f /* function of which the cube is to be made a prime */)
{
    DdNode *res;
    if (!Cudd_bddLeq(dd,cube,f)) return(NULL);
    do {
	dd->reordered = 0;
	res = cuddBddMakePrime(dd,cube,f);
    } while (dd->reordered == 1);
    return(res);
} /* end of Cudd_bddMakePrime */
/**Function********************************************************************
  Synopsis    [Expands lb to prime implicants of (f and ub).]
  Description [Expands lb to all prime implicants of (f and ub) that contain lb.
  Assumes that lb is contained in ub.  Returns the disjunction of the primes if
  lb is contained in f; returns the zero BDD if lb is not contained in f;
  returns NULL in case of failure.  In particular, NULL is returned if cube is
  not a real cube or is not an implicant of f.  Returning the disjunction of
  all prime implicants works because the resulting function is unate.]
  SideEffects [None]
  SeeAlso     [Cudd_bddMakePrime]
******************************************************************************/
DdNode *
Cudd_bddMaximallyExpand(
  DdManager *dd /* manager */,
  DdNode *lb /* cube to be expanded */,
  DdNode *ub /* upper bound cube */,
  DdNode *f /* function against which to expand */)
{
    DdNode *res;
    if (!Cudd_bddLeq(dd,lb,ub)) return(NULL);
    do {
	dd->reordered = 0;
	res = ddBddMaximallyExpand(dd,lb,ub,f);
    } while (dd->reordered == 1);
    return(res);
} /* end of Cudd_bddMaximallyExpand */
/**Function********************************************************************
  Synopsis    [Find a largest prime of a unate function.]
  Description [Find a largest prime implicant of a unate function.
  Returns the BDD for the prime if succesful; NULL otherwise.  The behavior
  is undefined if f is not unate.  The third argument is used to determine
  whether f is unate positive (increasing) or negative (decreasing)
  in each of the variables in its support.]
  SideEffects [None]
  SeeAlso     [Cudd_bddMaximallyExpand]
******************************************************************************/
DdNode *
Cudd_bddLargestPrimeUnate(
  DdManager *dd /* manager */,
  DdNode *f /* unate function */,
  DdNode *phaseBdd /* cube of the phases */)
{
    DdNode *res;
    int *phases;
    int retval;
    st_table *table;
    /* Extract phase vector for quick access. */
    phases = ALLOC(int, dd->size);
    if (phases == NULL) return(NULL);
    retval = Cudd_BddToCubeArray(dd, phaseBdd, phases);
    if (retval == 0) {
        FREE(phases);
        return(NULL);
    }
    do {
        dd->reordered = 0;
        table = st_init_table(st_ptrcmp,st_ptrhash);
        if (table == NULL) {
            FREE(phases);
            return(NULL);
        }
	(void) ddBddShortestPathUnate(dd, f, phases, table);
        res = ddGetLargestCubeUnate(dd, f, phases, table);
        st_free_table(table);
    } while (dd->reordered == 1);
    FREE(phases);
    return(res);
} /* end of Cudd_bddLargestPrimeUnate */
/*---------------------------------------------------------------------------*/
/* Definition of internal functions                                          */
/*---------------------------------------------------------------------------*/
/**Function********************************************************************
  Synopsis    [Performs the recursive step of Cudd_bddMakePrime.]
  Description [Performs the recursive step of Cudd_bddMakePrime.
  Returns the prime if successful; NULL otherwise.]
  SideEffects [None]
  SeeAlso     []
******************************************************************************/
DdNode *
cuddBddMakePrime(
  DdManager *dd /* manager */,
  DdNode *cube /* cube to be expanded */,
  DdNode *f /* function of which the cube is to be made a prime */)
{
    DdNode *scan;
    DdNode *t, *e;
    DdNode *res = cube;
    DdNode *zero = Cudd_Not(DD_ONE(dd));
    Cudd_Ref(res);
    scan = cube;
    while (!Cudd_IsConstant(scan)) {
	DdNode *reg = Cudd_Regular(scan);
	DdNode *var = dd->vars[reg->index];
	DdNode *expanded = Cudd_bddExistAbstract(dd,res,var);
	if (expanded == NULL) {
            Cudd_RecursiveDeref(dd,res);
	    return(NULL);
	}
	Cudd_Ref(expanded);
	if (Cudd_bddLeq(dd,expanded,f)) {
	    Cudd_RecursiveDeref(dd,res);
	    res = expanded;
	} else {
	    Cudd_RecursiveDeref(dd,expanded);
	}
	cuddGetBranches(scan,&t,&e);
	if (t == zero) {
	    scan = e;
	} else if (e == zero) {
	    scan = t;
	} else {
	    Cudd_RecursiveDeref(dd,res);
	    return(NULL);	/* cube is not a cube */
	}
    }
    if (scan == DD_ONE(dd)) {
	Cudd_Deref(res);
	return(res);
    } else {
	Cudd_RecursiveDeref(dd,res);
	return(NULL);
    }
} /* end of cuddBddMakePrime */
/*---------------------------------------------------------------------------*/
/* Definition of static functions                                            */
/*---------------------------------------------------------------------------*/
/**Function********************************************************************
  Synopsis    [Frees the entries of the visited symbol table.]
  Description [Frees the entries of the visited symbol table. Returns
  ST_CONTINUE.]
  SideEffects [None]
******************************************************************************/
static enum st_retval
freePathPair(
  char * key,
  char * value,
  char * arg)
{
    cuddPathPair *pair;
    pair = (cuddPathPair *) value;
	FREE(pair);
    return(ST_CONTINUE);
} /* end of freePathPair */
/**Function********************************************************************
  Synopsis    [Finds the length of the shortest path(s) in a DD.]
  Description [Finds the length of the shortest path(s) in a DD.
  Uses a local symbol table to store the lengths for each
  node. Only the lengths for the regular nodes are entered in the table,
  because those for the complement nodes are simply obtained by swapping
  the two lenghts.
  Returns a pair of lengths: the length of the shortest path to 1;
  and the length of the shortest path to 0. This is done so as to take
  complement arcs into account.]
  SideEffects [Accumulates the support of the DD in support.]
  SeeAlso     []
******************************************************************************/
static cuddPathPair
getShortest(
  DdNode * root,
  int * cost,
  int * support,
  st_table * visited)
{
    cuddPathPair *my_pair, res_pair, pair_T, pair_E;
    DdNode	*my_root, *T, *E;
    int		weight;
    my_root = Cudd_Regular(root);
    if (st_lookup(visited, my_root, &my_pair)) {
	if (Cudd_IsComplement(root)) {
	    res_pair.pos = my_pair->neg;
	    res_pair.neg = my_pair->pos;
	} else {
	    res_pair.pos = my_pair->pos;
	    res_pair.neg = my_pair->neg;
	}
	return(res_pair);
    }
    /* In the case of a BDD the following test is equivalent to
    ** testing whether the BDD is the constant 1. This formulation,
    ** however, works for ADDs as well, by assuming the usual
    ** dichotomy of 0 and != 0.
    */
    if (cuddIsConstant(my_root)) {
	if (my_root != zero) {
	    res_pair.pos = 0;
	    res_pair.neg = DD_BIGGY;
	} else {
	    res_pair.pos = DD_BIGGY;
	    res_pair.neg = 0;
	}
    } else {
	T = cuddT(my_root);
	E = cuddE(my_root);
	pair_T = getShortest(T, cost, support, visited);
	pair_E = getShortest(E, cost, support, visited);
	weight = WEIGHT(cost, my_root->index);
	res_pair.pos = ddMin(pair_T.pos+weight, pair_E.pos);
	res_pair.neg = ddMin(pair_T.neg+weight, pair_E.neg);
	/* Update support. */
	if (support != NULL) {
	    support[my_root->index] = 1;
	}
    }
    my_pair = ALLOC(cuddPathPair, 1);
    if (my_pair == NULL) {
	if (Cudd_IsComplement(root)) {
	    int tmp = res_pair.pos;
	    res_pair.pos = res_pair.neg;
	    res_pair.neg = tmp;
	}
	return(res_pair);
    }
    my_pair->pos = res_pair.pos;
    my_pair->neg = res_pair.neg;
    st_insert(visited, (char *)my_root, (char *)my_pair);
    if (Cudd_IsComplement(root)) {
	res_pair.pos = my_pair->neg;
	res_pair.neg = my_pair->pos;
    } else {
	res_pair.pos = my_pair->pos;
	res_pair.neg = my_pair->neg;
    }
    return(res_pair);
} /* end of getShortest */
/**Function********************************************************************
  Synopsis    [Build a BDD for a shortest path of f.]
  Description [Build a BDD for a shortest path of f.
  Given the minimum length from the root, and the minimum
  lengths for each node (in visited), apply triangulation at each node.
  Of the two children of each node on a shortest path, at least one is
  on a shortest path. In case of ties the procedure chooses the THEN
  children.
  Returns a pointer to the cube BDD representing the path if
  successful; NULL otherwise.]
  SideEffects [None]
  SeeAlso     []
******************************************************************************/
static DdNode *
getPath(
  DdManager * manager,
  st_table * visited,
  DdNode * f,
  int * weight,
  int  cost)
{
    DdNode	*sol, *tmp;
    DdNode	*my_dd, *T, *E;
    cuddPathPair *T_pair, *E_pair;
    int		Tcost, Ecost;
    int		complement;
    my_dd = Cudd_Regular(f);
    complement = Cudd_IsComplement(f);
    sol = one;
    cuddRef(sol);
    while (!cuddIsConstant(my_dd)) {
	Tcost = cost - WEIGHT(weight, my_dd->index);
	Ecost = cost;
	T = cuddT(my_dd);
	E = cuddE(my_dd);
	if (complement) {T = Cudd_Not(T); E = Cudd_Not(E);}
	st_lookup(visited, Cudd_Regular(T), &T_pair);
	if ((Cudd_IsComplement(T) && T_pair->neg == Tcost) ||
	(!Cudd_IsComplement(T) && T_pair->pos == Tcost)) {
	    tmp = cuddBddAndRecur(manager,manager->vars[my_dd->index],sol);
	    if (tmp == NULL) {
		Cudd_RecursiveDeref(manager,sol);
		return(NULL);
	    }
	    cuddRef(tmp);
	    Cudd_RecursiveDeref(manager,sol);
	    sol = tmp;
	    complement =  Cudd_IsComplement(T);
	    my_dd = Cudd_Regular(T);
	    cost = Tcost;
	    continue;
	}
	st_lookup(visited, Cudd_Regular(E), &E_pair);
	if ((Cudd_IsComplement(E) && E_pair->neg == Ecost) ||
	(!Cudd_IsComplement(E) && E_pair->pos == Ecost)) {
	    tmp = cuddBddAndRecur(manager,Cudd_Not(manager->vars[my_dd->index]),sol);
	    if (tmp == NULL) {
		Cudd_RecursiveDeref(manager,sol);
		return(NULL);
	    }
	    cuddRef(tmp);
	    Cudd_RecursiveDeref(manager,sol);
	    sol = tmp;
	    complement = Cudd_IsComplement(E);
	    my_dd = Cudd_Regular(E);
	    cost = Ecost;
	    continue;
	}
	(void) fprintf(manager->err,"We shouldn't be here!!\n");
	manager->errorCode = CUDD_INTERNAL_ERROR;
	return(NULL);
    }
    cuddDeref(sol);
    return(sol);
} /* end of getPath */
/**Function********************************************************************
  Synopsis    [Finds the size of the largest cube(s) in a DD.]
  Description [Finds the size of the largest cube(s) in a DD.
  This problem is translated into finding the shortest paths from a node
  when both THEN and ELSE arcs have unit lengths.
  Uses a local symbol table to store the lengths for each
  node. Only the lengths for the regular nodes are entered in the table,
  because those for the complement nodes are simply obtained by swapping
  the two lenghts.
  Returns a pair of lengths: the length of the shortest path to 1;
  and the length of the shortest path to 0. This is done so as to take
  complement arcs into account.]
  SideEffects [none]
  SeeAlso     []
******************************************************************************/
static cuddPathPair
getLargest(
  DdNode * root,
  st_table * visited)
{
    cuddPathPair *my_pair, res_pair, pair_T, pair_E;
    DdNode	*my_root, *T, *E;
    my_root = Cudd_Regular(root);
    if (st_lookup(visited, my_root, &my_pair)) {
	if (Cudd_IsComplement(root)) {
	    res_pair.pos = my_pair->neg;
	    res_pair.neg = my_pair->pos;
	} else {
	    res_pair.pos = my_pair->pos;
	    res_pair.neg = my_pair->neg;
	}
	return(res_pair);
    }
    /* In the case of a BDD the following test is equivalent to
    ** testing whether the BDD is the constant 1. This formulation,
    ** however, works for ADDs as well, by assuming the usual
    ** dichotomy of 0 and != 0.
    */
    if (cuddIsConstant(my_root)) {
	if (my_root != zero) {
	    res_pair.pos = 0;
	    res_pair.neg = DD_BIGGY;
	} else {
	    res_pair.pos = DD_BIGGY;
	    res_pair.neg = 0;
	}
    } else {
	T = cuddT(my_root);
	E = cuddE(my_root);
	pair_T = getLargest(T, visited);
	pair_E = getLargest(E, visited);
	res_pair.pos = ddMin(pair_T.pos, pair_E.pos) + 1;
	res_pair.neg = ddMin(pair_T.neg, pair_E.neg) + 1;
    }
    my_pair = ALLOC(cuddPathPair, 1);
    if (my_pair == NULL) {	/* simply do not cache this result */
	if (Cudd_IsComplement(root)) {
	    int tmp = res_pair.pos;
	    res_pair.pos = res_pair.neg;
	    res_pair.neg = tmp;
	}
	return(res_pair);
    }
    my_pair->pos = res_pair.pos;
    my_pair->neg = res_pair.neg;
    /* Caching may fail without affecting correctness. */
    st_insert(visited, (char *)my_root, (char *)my_pair);
    if (Cudd_IsComplement(root)) {
	res_pair.pos = my_pair->neg;
	res_pair.neg = my_pair->pos;
    } else {
	res_pair.pos = my_pair->pos;
	res_pair.neg = my_pair->neg;
    }
    return(res_pair);
} /* end of getLargest */
/**Function********************************************************************
  Synopsis    [Build a BDD for a largest cube of f.]
  Description [Build a BDD for a largest cube of f.
  Given the minimum length from the root, and the minimum
  lengths for each node (in visited), apply triangulation at each node.
  Of the two children of each node on a shortest path, at least one is
  on a shortest path. In case of ties the procedure chooses the THEN
  children.
  Returns a pointer to the cube BDD representing the path if
  successful; NULL otherwise.]
  SideEffects [None]
  SeeAlso     []
******************************************************************************/
static DdNode *
getCube(
  DdManager * manager,
  st_table * visited,
  DdNode * f,
  int  cost)
{
    DdNode	*sol, *tmp;
    DdNode	*my_dd, *T, *E;
    cuddPathPair *T_pair, *E_pair;
    int		Tcost, Ecost;
    int		complement;
    my_dd = Cudd_Regular(f);
    complement = Cudd_IsComplement(f);
    sol = one;
    cuddRef(sol);
    while (!cuddIsConstant(my_dd)) {
	Tcost = cost - 1;
	Ecost = cost - 1;
	T = cuddT(my_dd);
	E = cuddE(my_dd);
	if (complement) {T = Cudd_Not(T); E = Cudd_Not(E);}
	if (!st_lookup(visited, Cudd_Regular(T), &T_pair)) return(NULL);
	if ((Cudd_IsComplement(T) && T_pair->neg == Tcost) ||
	(!Cudd_IsComplement(T) && T_pair->pos == Tcost)) {
	    tmp = cuddBddAndRecur(manager,manager->vars[my_dd->index],sol);
	    if (tmp == NULL) {
		Cudd_RecursiveDeref(manager,sol);
		return(NULL);
	    }
	    cuddRef(tmp);
	    Cudd_RecursiveDeref(manager,sol);
	    sol = tmp;
	    complement =  Cudd_IsComplement(T);
	    my_dd = Cudd_Regular(T);
	    cost = Tcost;
	    continue;
	}
	if (!st_lookup(visited, Cudd_Regular(E), &E_pair)) return(NULL);
	if ((Cudd_IsComplement(E) && E_pair->neg == Ecost) ||
	(!Cudd_IsComplement(E) && E_pair->pos == Ecost)) {
	    tmp = cuddBddAndRecur(manager,Cudd_Not(manager->vars[my_dd->index]),sol);
	    if (tmp == NULL) {
		Cudd_RecursiveDeref(manager,sol);
		return(NULL);
	    }
	    cuddRef(tmp);
	    Cudd_RecursiveDeref(manager,sol);
	    sol = tmp;
	    complement = Cudd_IsComplement(E);
	    my_dd = Cudd_Regular(E);
	    cost = Ecost;
	    continue;
	}
	(void) fprintf(manager->err,"We shouldn't be here!\n");
	manager->errorCode = CUDD_INTERNAL_ERROR;
	return(NULL);
    }
    cuddDeref(sol);
    return(sol);
} /* end of getCube */
/**Function********************************************************************
  Synopsis    [Performs the recursive step of Cudd_bddMaximallyExpand.]
  Description [Performs the recursive step of Cudd_bddMaximallyExpand.
  Returns set of primes or zero BDD if successful; NULL otherwise.  On entry
  to this function, ub and lb should be different from the zero BDD.  The
  function then maintains this invariant.]
  SideEffects [None]
  SeeAlso     []
******************************************************************************/
static DdNode *
ddBddMaximallyExpand(
  DdManager *dd /* manager */,
  DdNode *lb /* cube to be expanded */,
  DdNode *ub /* upper bound cube */,
  DdNode *f /* function against which to expand */)
{
    DdNode *one, *zero, *lbv, *lbvn, *lbnx, *ubv, *ubvn, *fv, *fvn, *res;
    DdNode *F, *UB, *LB, *t, *e;
    unsigned int top, toplb, topub, topf, index;
    statLine(dd);
    /* Terminal cases. */
    one = DD_ONE(dd);
    zero = Cudd_Not(one);
    assert(ub != zero && lb != zero);
    /** There are three major terminal cases in theory:
     **   ub -> f     : return ub
     **   lb == f     : return lb
     **   not(lb -> f): return zero
     ** Only the second case can be checked exactly in constant time.
     ** For the others, we check for sufficient conditions.
     */
    if (ub == f || f == one) return(ub);
    if (lb == f) return(lb);
    if (f == zero || ub == Cudd_Not(f) || lb == one || lb == Cudd_Not(f))
        return(zero);
    if (!Cudd_IsComplement(lb) && Cudd_IsComplement(f)) return(zero);
    /* Here lb and f are not constant. */
    /* Check cache.  Since lb and ub are cubes, their local reference counts
    ** are always 1.  Hence, we only check the reference count of f.
    */
    F = Cudd_Regular(f);
    if (F->ref != 1) {
        DdNode *tmp = cuddCacheLookup(dd, DD_BDD_MAX_EXP_TAG, lb, ub, f);
        if (tmp != NULL) {
            return(tmp);
        }
    }
    /* Compute cofactors.  For lb we use the non-zero one in
    ** both branches of the recursion.
    */
    LB = Cudd_Regular(lb);
    UB = Cudd_Regular(ub);
    topf = dd->perm[F->index];
    toplb = dd->perm[LB->index];
    topub = (ub == one) ? CUDD_CONST_INDEX : dd->perm[UB->index];
    assert(toplb <= topub);
    top = ddMin(topf,toplb);
    if (toplb == top) {
	index = LB->index;
        lbv = cuddT(LB);
        lbvn = cuddE(LB);
        if (lb != LB) {
            lbv = Cudd_Not(lbv);
            lbvn = Cudd_Not(lbvn);
        }
        if (lbv == zero) {
            lbnx = lbvn;
        } else {
            lbnx = lbv;
        }
    } else {
	index = F->index;
        lbnx = lbv = lbvn = lb;
    }
    if (topub == top) {
        ubv = cuddT(UB);
        ubvn = cuddE(UB);
        if (ub != UB) {
            ubv = Cudd_Not(ubv);
            ubvn = Cudd_Not(ubvn);
        }
    } else {
        ubv = ubvn = ub;
    }
    if (topf == top) {
        fv = cuddT(F);
        fvn = cuddE(F);
        if (f != F) {
            fv = Cudd_Not(fv);
            fvn = Cudd_Not(fvn);
        }
    } else {
        fv = fvn = f;
    }
    /* Recursive calls. */
    if (ubv != zero) {
        t = ddBddMaximallyExpand(dd, lbnx, ubv, fv);
        if (t == NULL) return(NULL);
    } else {
        assert(topub == toplb && topub == top && lbv == zero);
        t = zero;
    }
    cuddRef(t);
    /* If the top variable appears only in lb, the positive and negative
    ** cofactors of each operand are the same.  We want to avoid a
    ** needless recursive call, which would force us to give up the
    ** cache optimization trick based on reference counts.
    */
    if (ubv == ubvn && fv == fvn) {
        res = t;
    } else {
        if (ubvn != zero) {
            e = ddBddMaximallyExpand(dd, lbnx, ubvn, fvn);
            if (e == NULL) {
                Cudd_IterDerefBdd(dd,t);
                return(NULL);
            }
        } else {
            assert(topub == toplb && topub == top && lbvn == zero);
            e = zero;
        }
        if (t == e) {
            res = t;
        } else {
            cuddRef(e);
            if (toplb == top) {
                if (lbv == zero) {
                    /* Top variable appears in negative phase. */
                    if (t != one) {
                        DdNode *newT;
                        if (Cudd_IsComplement(t)) {
                            newT = cuddUniqueInter(dd, index, Cudd_Not(t), zero);
                            if (newT == NULL) {
                                Cudd_IterDerefBdd(dd,t);
                                Cudd_IterDerefBdd(dd,e);
                                return(NULL);
                            }
                            newT = Cudd_Not(newT);
                        } else {
                            newT = cuddUniqueInter(dd, index, t, one);
                            if (newT == NULL) {
                                Cudd_IterDerefBdd(dd,t);
                                Cudd_IterDerefBdd(dd,e);
                                return(NULL);
                            }
                        }
                        cuddRef(newT);
                        cuddDeref(t);
                        t = newT;
                    }
                } else if (lbvn == zero) {
                    /* Top variable appears in positive phase. */
                    if (e != one) {
                        DdNode *newE;
                        newE = cuddUniqueInter(dd, index, one, e);
                        if (newE == NULL) {
                            Cudd_IterDerefBdd(dd,t);
                            Cudd_IterDerefBdd(dd,e);
                            return(NULL);
                        }
                        cuddRef(newE);
                        cuddDeref(e);
                        e = newE;
                    }
                } else {
                    /* Not a cube. */
                    Cudd_IterDerefBdd(dd,t);
                    Cudd_IterDerefBdd(dd,e);
                    return(NULL);
                }
            }
            /* Combine results. */
            res = cuddBddAndRecur(dd, t, e);
            if (res == NULL) {
                Cudd_IterDerefBdd(dd,t);
                Cudd_IterDerefBdd(dd,e);
                return(NULL);
            }
            cuddRef(res);
            Cudd_IterDerefBdd(dd,t);
            Cudd_IterDerefBdd(dd,e);
        }
    }
    /* Cache result and return. */
    if (F->ref != 1) {
        cuddCacheInsert(dd, DD_BDD_MAX_EXP_TAG, lb, ub, f, res);
    }
    cuddDeref(res);
    return(res);
} /* end of ddBddMaximallyExpand */
/**Function********************************************************************
  Synopsis    [Performs shortest path computation on a unate function.]
  Description [Performs shortest path computation on a unate function.
  Returns the length of the shortest path to one if successful;
  CUDD_OUT_OF_MEM otherwise.  This function is based on the observation
  that in the BDD of a unate function no node except the constant is
  reachable from the root via paths of different parity.]
  SideEffects [None]
  SeeAlso     [getShortest]
******************************************************************************/
static int
ddBddShortestPathUnate(
  DdManager *dd,
  DdNode *f,
  int *phases,
  st_table *table)
{
    int positive, l, lT, lE;
    DdNode *one = DD_ONE(dd);
    DdNode *zero = Cudd_Not(one);
    DdNode *F, *fv, *fvn;
    if (st_lookup_int(table, f, &l)) {
        return(l);
    }
    if (f == one) {
        l = 0;
    } else if (f == zero) {
        l = DD_BIGGY;
    } else {
        F = Cudd_Regular(f);
        fv = cuddT(F);
        fvn = cuddE(F);
        if (f != F) {
            fv = Cudd_Not(fv);
            fvn = Cudd_Not(fvn);
        }
        lT = ddBddShortestPathUnate(dd, fv, phases, table);
        lE = ddBddShortestPathUnate(dd, fvn, phases, table);
        positive = phases[F->index];
        l = positive ? ddMin(lT+1, lE) : ddMin(lT, lE+1);
    }
    if (st_insert(table, f, (void *)(ptrint) l) == ST_OUT_OF_MEM) {
        return(CUDD_OUT_OF_MEM);
    }
    return(l);
} /* end of ddShortestPathUnate */
/**Function********************************************************************
  Synopsis    [Extracts largest prime of a unate function.]
  Description [Extracts largest prime of a unate function.  Returns the BDD of
  the prime if successful; NULL otherwise.]
  SideEffects [None]
  SeeAlso     [getPath]
******************************************************************************/
static DdNode *
ddGetLargestCubeUnate(
  DdManager *dd,
  DdNode *f,
  int *phases,
  st_table *table)
{
    DdNode *res, *scan;
    DdNode *one = DD_ONE(dd);
    int cost;
    res = one;
    cuddRef(res);
    scan = f;
    st_lookup_int(table, scan, &cost);
    while (!Cudd_IsConstant(scan)) {
        int Pcost, Ncost, Tcost;
        DdNode *tmp, *T, *E;
        DdNode *rscan = Cudd_Regular(scan);
        int index = rscan->index;
        assert(phases[index] == 0 || phases[index] == 1);
        int positive = phases[index] == 1;
        Pcost = positive ? cost - 1 : cost;
        Ncost = positive ? cost : cost - 1;
        T = cuddT(rscan);
        E = cuddE(rscan);
        if (rscan != scan) {
            T = Cudd_Not(T);
            E = Cudd_Not(E);
        }
        tmp = res;
        st_lookup_int(table, T, &Tcost);
        if (Tcost == Pcost) {
            cost = Pcost;
            scan = T;
            if (positive) {
                tmp = cuddBddAndRecur(dd, dd->vars[index], res);
            }
        } else {
            cost = Ncost;
            scan = E;
            if (!positive) {
                tmp = cuddBddAndRecur(dd, Cudd_Not(dd->vars[index]), res);
            }
        }
        if (tmp == NULL) {
            Cudd_IterDerefBdd(dd, res);
            return(NULL);
        }
        cuddRef(tmp);
        Cudd_IterDerefBdd(dd, res);
        res = tmp;
    }
    cuddDeref(res);
    return(res);
} /* end of ddGetLargestCubeUnate */