get_submodules.sh: take into account imath submodule of isl

[ppn.git] / eqv3.cc

#include <isl/constraint.h>
#include <isl/union_map.h>
#include <isl/union_set.h>
#include <isa/pdg.h>
#include "dependence_graph.h"
#include "eqv3_options.h"
#include "msa.h"

/* Check the equivalence of two programs, given as dependence graphs,
 * by computing what is essentially the cross product of the two
 * dependence graph and checking that we can only reach states corresponding
 * to equal inputs from the state corresponding to equal output.
 * We consider various ways of performing the reachability analysis:
 * applying the transitive closure to the initial state, FAST, Lever and Aspic.
 * In all but the first case, we simply print a model that can be used
 * as input to the corresponding tool.
 */

static unsigned update_out_dim(pdg::PDG *pdg, unsigned out_dim)
{
        for (int i = 0; i < pdg->arrays.size(); ++i) {
                pdg::array *array = pdg->arrays[i];
                if (array->type != pdg::array::output)
                        continue;
                if (array->dims.size() + 1 > out_dim)
                        out_dim = array->dims.size() + 1;
        }

        return out_dim;
}

/* The input arrays of the two programs are supposed to be the same,
 * so they should at least have the same dimension.  Make sure
 * this is true, because we depend on it later on.
 */
static int check_input_arrays(dependence_graph *dg1, dependence_graph *dg2)
{
        for (int i = 0; i < dg1->input_computations.size(); ++i)
                for (int j = 0; j < dg2->input_computations.size(); ++j) {
                        if (strcmp(dg1->input_computations[i]->operation,
                                   dg2->input_computations[j]->operation))
                                continue;
                        if (dg1->input_computations[i]->dim ==
                            dg2->input_computations[j]->dim)
                                continue;
                        fprintf(stderr,
                            "input arrays \"%s\" do not have the same dimension\n",
                            dg1->input_computations[i]->operation);
                        return -1;
                }

        return 0;
}

/* Compute a map between the domains of the output nodes of
 * the dependence graphs that expresses that array1 in dg1
 * and array2 in dg2 are equal.
 * This map is then returned as a set.
 */
static __isl_give isl_set *compute_equal_array(
        dependence_graph *dg1, dependence_graph *dg2, int array1, int array2)
{
        const char *array_name = dg1->output_arrays[array1];
        isl_set *out1 = isl_set_copy(dg1->out->domain);
        isl_set *out2 = isl_set_copy(dg2->out->domain);
        unsigned dim1 = isl_set_n_dim(out1);
        unsigned dim2 = isl_set_n_dim(out2);

        out1 = isl_set_fix_dim_si(out1, dim1 - 1, array1);
        out2 = isl_set_fix_dim_si(out2, dim2 - 1, array2);
        out1 = isl_set_remove_dims(out1, isl_dim_set, dim1 - 1, 1);
        out2 = isl_set_remove_dims(out2, isl_dim_set, dim2 - 1, 1);
        int equal = isl_set_is_equal(out1, out2);
        isl_set_free(out1);
        isl_set_free(out2);
        assert(equal >= 0);
        if (!equal) {
                fprintf(stderr, "different output domains for array %s\n",
                                array_name);
                return NULL;
        }

        out1 = isl_set_copy(dg1->out->domain);
        out2 = isl_set_copy(dg2->out->domain);
        out1 = isl_set_fix_dim_si(out1, dim1 - 1, array1);
        out2 = isl_set_fix_dim_si(out2, dim2 - 1, array2);
        out1 = isl_set_coalesce(out1);
        out2 = isl_set_coalesce(out2);
        isl_space *dim = isl_space_map_from_set(isl_set_get_space(out2));
        isl_map *id = isl_map_identity(dim);
        id = isl_map_remove_dims(id, isl_dim_in, isl_map_n_in(id) - 1, 1);
        id = isl_map_apply_domain(id, isl_map_copy(id));
        id = isl_map_intersect_domain(id, out1);
        id = isl_map_intersect_range(id, out2);

        return isl_map_wrap(id);
}

/* Compute a map between the domains of the output nodes of
 * the dependence graphs that expresses that arrays with
 * the same names in both graphs are equal.
 * The map is returned as a set.
 */
static __isl_give isl_union_set *compute_equal_output(
        dependence_graph *dg1, dependence_graph *dg2)
{
        int i1 = 0, i2 = 0;
        isl_space *dim;
        isl_map *map;
        isl_set *equal_output;
        unsigned d;

        dim = isl_set_get_space(dg1->out->domain);
        d = isl_space_dim(dim, isl_dim_set);
        map = isl_map_empty(isl_space_map_from_set(dim));
        equal_output = isl_map_wrap(map);

        while (i1 < dg1->output_arrays.size() || i2 < dg2->output_arrays.size()) {
                int cmp;
                cmp = i1 == dg1->output_arrays.size() ? 1 :
                      i2 == dg2->output_arrays.size() ? -1 :
                      strcmp(dg1->output_arrays[i1], dg2->output_arrays[i2]);
                if (cmp < 0) {
                        fprintf(stderr,
                                "Array \"%s\" only in program 1; not checked\n",
                                dg1->output_arrays[i1]);
                        ++i1;
                } else if (cmp > 0) {
                        fprintf(stderr,
                                "Array \"%s\" only in program 2; not checked\n",
                                dg2->output_arrays[i2]);
                        ++i2;
                } else {
                        isl_set *eq_array;
                        eq_array = compute_equal_array(dg1, dg2, i1, i2);
                        equal_output = isl_set_union(equal_output, eq_array);
                        ++i1;
                        ++i2;
                }
        }

        equal_output = isl_set_set_tuple_name(equal_output, "s0");

        return isl_union_set_from_set(equal_output);
}

/* For a given state "s" of the MSA, return a set describing all
 * instances of that state.
 */
static __isl_give isl_union_set *state_domain(MSA &msa, int s)
{
        char name[50];
        isl_space *dim1, *dim2;
        isl_set *set;

        snprintf(name, sizeof(name), "s%d", s);
        dim1 = isl_set_get_space(msa.state[s].first->domain);
        dim2 = isl_set_get_space(msa.state[s].second->domain);
        dim1 = isl_space_from_domain(dim1);
        dim2 = isl_space_from_range(dim2);
        dim1 = isl_space_join(dim1, dim2);
        dim1 = isl_space_wrap(dim1);
        dim1 = isl_space_set_tuple_name(dim1, isl_dim_set, name);
        set = isl_set_universe(dim1);

        return isl_union_set_from_set(set);
}

/* For a given state corresponding to a pair of equal input arrays,
 * return a set describing pairs of equal instances of those arrays.
 */
static __isl_give isl_union_set *equal_input_domain(MSA &msa, int s)
{
        char name[50];
        isl_space *dim;
        isl_set *set;

        snprintf(name, sizeof(name), "s%d", s);
        dim = isl_set_get_space(msa.state[s].first->domain);
        set = isl_map_wrap(isl_map_identity(isl_space_map_from_set(dim)));
        set = isl_set_set_tuple_name(set, name);

        return isl_union_set_from_set(set);
}

/* Return all instances of final states, i.e., pairs of equal input
 * arrays, pairs of equal constant functions and all bad combinations.
 */
static __isl_give isl_union_set *compute_final_states(isl_ctx *ctx, MSA &msa)
{
        isl_space *dim;
        isl_union_set *final;

        dim = isl_space_set_alloc(ctx, 0, 0);
        final = isl_union_set_empty(dim);

        for (int i = 0; i < msa.input_state.size(); ++i) {
                int s = msa.input_state[i];
                final = isl_union_set_union(final, state_domain(msa, s));
        }

        for (int i = 0; i < msa.fail_state.size(); ++i) {
                int s = msa.fail_state[i];
                final = isl_union_set_union(final, state_domain(msa, s));
        }

        for (int i = 0; i < msa.cst_state.size(); ++i) {
                int s = msa.cst_state[i];
                final = isl_union_set_union(final, state_domain(msa, s));
        }

        return final;
}

/* Return all instances of good final states, i.e., equal instances
 * of equal input arrays and pairs of equal constant functions.
 */
static __isl_give isl_union_set *compute_good_states(isl_ctx *ctx, MSA &msa)
{
        isl_space *dim;
        isl_union_set *final;

        dim = isl_space_set_alloc(ctx, 0, 0);
        final = isl_union_set_empty(dim);

        for (int i = 0; i < msa.input_state.size(); ++i) {
                int s = msa.input_state[i];
                final = isl_union_set_union(final, equal_input_domain(msa, s));
        }

        for (int i = 0; i < msa.cst_state.size(); ++i) {
                int s = msa.cst_state[i];
                final = isl_union_set_union(final, state_domain(msa, s));
        }

        return final;
}

/* Perform the reachability analysis by first computing the
 * transitive closure and applying it to a set representing
 * equal output.  We then check if all reachable final states
 * are good final states.
 * If we find a bad state that may be reachable, we also
 * apply the inverse of the transitive closure to that state
 * to find out where it may have come from.
 */
static int check_equivalence_tc(isl_ctx *ctx,
        dependence_graph *dg1, dependence_graph *dg2)
{
        int root;
        MSA msa;
        isl_space *dim;
        isl_union_map *transition;
        isl_union_set *equal_output;
        isl_union_set *test;
        isl_union_set *final;
        isl_union_set *good;
        isl_printer *p;
        int ok = 1;

        if (check_input_arrays(dg1, dg2))
                return -1;

        p = isl_printer_to_file(ctx, stdout);

        root = msa.add(dg1->out, dg2->out);

        dim = isl_space_alloc(ctx, 0, 0, 0);
        transition = isl_union_map_empty(dim);
        for (int i = 0; i < msa.firing.size(); ++i) {
                isl_map *map = isl_map_copy(msa.firing[i]);
                transition = isl_union_map_union(transition,
                                                 isl_union_map_from_map(map));
        }

        transition = isl_union_map_transitive_closure(transition, NULL);

        equal_output = compute_equal_output(dg1, dg2);
        test = isl_union_set_apply(isl_union_set_copy(equal_output),
                                    isl_union_map_copy(transition));

        final = compute_final_states(ctx, msa);
        test = isl_union_set_intersect(test, final);

        good = compute_good_states(ctx, msa);
        test = isl_union_set_subtract(test, good);

        if (!isl_union_set_is_empty(test)) {
                ok = 0;
                transition = isl_union_map_reverse(transition);
                test = isl_union_set_apply(test, transition);
                test = isl_union_set_coalesce(test);
                test = isl_union_set_intersect(test, equal_output);

                isl_printer_print_str(p, "NOT equivalent");
                isl_printer_end_line(p);
                isl_printer_print_union_set(p, test);
                isl_printer_end_line(p);
        } else {
                isl_union_set_free(equal_output);
                isl_union_map_free(transition);
        }

        isl_union_set_free(test);

        if (!ok)
                isl_printer_print_str(p, "NOT ");
        isl_printer_print_str(p, "OK");
        isl_printer_end_line(p);

        isl_printer_free(p);

        return 0;
}

struct transition_info {
        FILE *out;
        int nvar;               /* Total number of variables */
        int state;              /* Next extra state/total number of states */
        int transition;         /* Current transition */
        int from;               /* From state of current firing relation */
        int to;                 /* To state of current firing relation */

        int first;              /* Are we printing the first expression? */
        int last_out;           /* Last variable that was explicitly assigned */
        int translation;        /* Current firing relation is a translation */

        int lever;              /* Are we printing Lever output? */
};

static int is_translation(__isl_keep isl_basic_map *bmap)
{
        unsigned n_out = isl_basic_map_dim(bmap, isl_dim_out);
        unsigned n_div = isl_basic_map_dim(bmap, isl_dim_div);

        if (n_div > 0)
                return 0;

        for (int i = 0; i < n_out; ++i) {
                isl_constraint *c = NULL;
                if (!isl_basic_map_has_defining_equality(bmap,
                                                         isl_dim_out, i, &c))
                        return 0;
                isl_constraint_free(c);
        }

        return 1;
}

/* Update the maximal number of variables and the total number of states
 * based on the given firing relation.
 * If the firing relation is a translation, then we don't need anything extra.
 * Otherwise, we need an extra state and
 * to encode the firing relation we need at most as many variables
 * as are needed to describe the relation.  Note that this does
 * not include the parameters as they are treated separately.
 */
static isl_stat check_transition(__isl_take isl_basic_map *firing, void *user)
{
        transition_info *ti = (transition_info *)user;

        if (!is_translation(firing)) {
                unsigned n_in = isl_basic_map_dim(firing, isl_dim_in);
                unsigned n_out = isl_basic_map_dim(firing, isl_dim_out);
                unsigned n_div = isl_basic_map_dim(firing, isl_dim_div);
                unsigned total = n_in + n_out + n_div;

                if (total > ti->nvar)
                        ti->nvar = total;
                ti->state++;
        }

        isl_basic_map_free(firing);
        return isl_stat_ok;
}

/* Compute the maximal number of variables needed as well as the
 * total number of states.
 * The number of variables needs to be at least the maximal
 * number of iterators needed to describe a state in the MSA.
 * For transitions that are not translations, we may need more.
 */
static void compute_max_var(MSA &msa, struct transition_info *ti)
{
        for (int i = 0; i < msa.state.size(); ++i) {
                computation *c1 = msa.state[i].first;
                computation *c2 = msa.state[i].second;
                unsigned d1 = isl_set_dim(c1->domain, isl_dim_set);
                unsigned d2 = isl_set_dim(c2->domain, isl_dim_set);

                if (d1 + d2 > ti->nvar)
                        ti->nvar = d1 + d2;
        }
        for (int i = 0; i < msa.transition.size(); ++i) {
                int f = msa.transition[i].second;
                isl_map_foreach_basic_map(msa.firing[f],
                        &check_transition, ti);
        }
}

/* Print the constraint as a linear combination.
 * Because Lever doesn't like linear expressions that start with
 * "-", we always start off with a non-negative constant (possibly
 * zero).
 * "s" is the sign with which to multiply the linear combination.
 * For constraints that end up in the action, we shouldn't print
 * the coefficient of the output variable being defined.  In this
 * case "no_output" is set to 1.
 */
static __isl_give isl_printer *print_linear(__isl_take isl_printer *p,
        __isl_keep isl_constraint *c, int s, int no_output)
{
        isl_val *v;
        int nparam;
        int n_in, n_out, n_div;

        nparam = isl_constraint_dim(c, isl_dim_param);
        n_in = isl_constraint_dim(c, isl_dim_in);
        n_out = isl_constraint_dim(c, isl_dim_out);
        n_div = isl_constraint_dim(c, isl_dim_div);

        v = isl_constraint_get_constant_val(c);
        if (s < 0)
                v = isl_val_neg(v);
        if (isl_val_is_neg(v))
                p = isl_printer_print_str(p, "0");
        p = isl_printer_print_val(p, v);
        isl_val_free(v);
        for (int i = 0; i < nparam; ++i) {
                const char *name;
                v = isl_constraint_get_coefficient_val(c, isl_dim_param, i);
                if (isl_val_is_zero(v)) {
                        isl_val_free(v);
                        continue;
                }
                if (s < 0)
                        v = isl_val_neg(v);
                if (isl_val_is_pos(v))
                        p = isl_printer_print_str(p, "+");
                p = isl_printer_print_val(p, v);
                name = isl_constraint_get_dim_name(c, isl_dim_param, i);
                p = isl_printer_print_str(p, "*");
                p = isl_printer_print_str(p, name);
                isl_val_free(v);
        }
        for (int i = 0; i < n_in; ++i) {
                v = isl_constraint_get_coefficient_val(c, isl_dim_in, i);
                if (isl_val_is_zero(v)) {
                        isl_val_free(v);
                        continue;
                }
                if (s < 0)
                        v = isl_val_neg(v);
                if (isl_val_is_pos(v))
                        p = isl_printer_print_str(p, "+");
                p = isl_printer_print_val(p, v);
                p = isl_printer_print_str(p, "*x");
                p = isl_printer_print_int(p, i);
                isl_val_free(v);
        }
        if (!no_output) {
                for (int i = 0; i < n_out; ++i) {
                        v = isl_constraint_get_coefficient_val(c,
                                                                isl_dim_out, i);
                        if (isl_val_is_zero(v)) {
                                isl_val_free(v);
                                continue;
                        }
                        if (s < 0)
                                v = isl_val_neg(v);
                        if (isl_val_is_pos(v))
                                p = isl_printer_print_str(p, "+");
                        p = isl_printer_print_val(p, v);
                        p = isl_printer_print_str(p, "*x");
                        p = isl_printer_print_int(p, n_in + i);
                        isl_val_free(v);
                }
                for (int i = 0; i < n_div; ++i) {
                        v = isl_constraint_get_coefficient_val(c,
                                                                isl_dim_div, i);
                        if (isl_val_is_zero(v)) {
                                isl_val_free(v);
                                continue;
                        }
                        if (s < 0)
                                v = isl_val_neg(v);
                        if (isl_val_is_pos(v))
                                p = isl_printer_print_str(p, "+");
                        p = isl_printer_print_val(p, v);
                        p = isl_printer_print_str(p, "*x");
                        p = isl_printer_print_int(p, n_in + n_out + i);
                        isl_val_free(v);
                }
        }

        return p;
}

/* Print the constraints that go into the guard.
 * If we are dealing with a translation, then we shouldn't print any
 * constraints involving the output variables.
 * Those constraint should be equalities defining the translation
 * and will be printed in the action.
 */
static isl_stat print_guard_constraint(__isl_take isl_constraint *c, void *user)
{
        transition_info *ti = (transition_info *)user;
        int n_out;
        int n_in;
        isl_ctx *ctx;
        isl_val *v;
        isl_printer *p;

        n_out = isl_constraint_dim(c, isl_dim_out);

        if (ti->translation) {
                for (int i = 0; i < n_out; ++i) {
                        v = isl_constraint_get_coefficient_val(c,
                                                            isl_dim_out, i);
                        if (!isl_val_is_zero(v)) {
                                assert(isl_constraint_is_equality(c));
                                isl_constraint_free(c);
                                isl_val_free(v);
                                return isl_stat_ok;
                        }
                }
        }

        ctx = isl_constraint_get_ctx(c);
        p = isl_printer_to_file(ctx, ti->out);
        if (!ti->first)
                p = isl_printer_print_str(p, " && ");
        p = print_linear(p, c, 1, 0);
        if (isl_constraint_is_equality(c))
                p = isl_printer_print_str(p, " = 0");
        else
                p = isl_printer_print_str(p, " >= 0");

        ti->first = 0;

        isl_printer_free(p);
        isl_constraint_free(c);
        isl_val_free(v);

        return isl_stat_ok;
}

/* Print the constraints that go into the action.
 * That is, equality constraints that define one of the output variables.
 */
static isl_stat print_action_constraint(__isl_take isl_constraint *c,
        void *user)
{
        transition_info *ti = (transition_info *)user;
        int n_in;
        int n_out;
        isl_ctx *ctx;
        isl_printer *p;
        int out;
        int s;

        if (!isl_constraint_is_equality(c)) {
                isl_constraint_free(c);
                return isl_stat_ok;
        }

        n_in = isl_constraint_dim(c, isl_dim_in);
        n_out = isl_constraint_dim(c, isl_dim_out);

        for (out = 0; out < n_out; ++out) {
                if (isl_constraint_involves_dims(c, isl_dim_out, out, 1))
                        break;
        }

        if (out >= n_out) {
                isl_constraint_free(c);
                return isl_stat_ok;
        }

        if (out > ti->last_out)
                ti->last_out = out;

        if (isl_constraint_is_lower_bound(c, isl_dim_out, out))
                s = -1;
        else
                s = 1;

        ctx = isl_constraint_get_ctx(c);
        p = isl_printer_to_file(ctx, ti->out);
        if (!ti->first)
                p = isl_printer_print_str(p, ", ");

        p = isl_printer_print_str(p, "x");
        p = isl_printer_print_int(p, out);
        p = isl_printer_print_str(p, "' = ");
        p = print_linear(p, c, s, 1);

        ti->first = 0;

        isl_constraint_free(c);
        return isl_stat_ok;
}

/* A firing relation that is not a translation is encoded using
 * several transitions.  We first allow an unconditional transition
 * to an extra state.  In this state, we allow arbitrary increments
 * on the extra variables that will correspond to the output variables
 * and the existentially quantified variables of the relation.  This allows
 * those variables to attain any value, since they are initially set
 * to zero and FAST and Lever only support non-negative variables.
 * Finally, we add a transition from the extra
 * state to the original target state with a guard that corresponds
 * to the firing relation and an action that assigns the output
 * variables to the main variables.
 */
static isl_stat print_non_translation(__isl_take isl_basic_map *firing,
        transition_info *ti)
{
        unsigned n_in = isl_basic_map_dim(firing, isl_dim_in);
        unsigned n_out = isl_basic_map_dim(firing, isl_dim_out);
        unsigned n_div = isl_basic_map_dim(firing, isl_dim_div);

        fprintf(ti->out, "transition t%d := {\n", ti->transition);
        if (!ti->lever) {
                fprintf(ti->out, "\tfrom := s%d;\n", ti->from);
                fprintf(ti->out, "\tto := s%d;\n", ti->state);
                fprintf(ti->out, "\tguard := true;\n");
                fprintf(ti->out, "\taction := x0' = x0");
        } else {
                fprintf(ti->out, "\tguard := state = %d;\n", ti->from);
                fprintf(ti->out, "\taction := state' = %d", ti->state);
        }
        for (int i = 0; i < n_out + n_div; ++i)
                fprintf(ti->out, ", x%d' = 0", n_in + i);
        fprintf(ti->out, ";\n");
        fprintf(ti->out, "};\n");

        ti->transition++;

        for (int i = 0; i < n_out + n_div; ++i) {
                fprintf(ti->out, "transition t%d := {\n", ti->transition);
                if (!ti->lever) {
                        fprintf(ti->out, "\tfrom := s%d;\n", ti->state);
                        fprintf(ti->out, "\tto := s%d;\n", ti->state);
                        fprintf(ti->out, "\tguard := true;\n");
                } else
                        fprintf(ti->out, "\tguard := state = %d;\n", ti->state);
                fprintf(ti->out, "\taction := x%d' = x%d + 1;\n",
                                    n_in + i, n_in + i);
                fprintf(ti->out, "};\n");
                ti->transition++;
        }

        fprintf(ti->out, "transition t%d := {\n", ti->transition);
        if (!ti->lever) {
                fprintf(ti->out, "\tfrom := s%d;\n", ti->state);
                fprintf(ti->out, "\tto := s%d;\n", ti->to);
        }
        fprintf(ti->out, "\tguard := ");
        if (!ti->lever)
                ti->first = 1;
        else {
                fprintf(ti->out, "state = %d", ti->state);
                ti->first = 0;
        }
        isl_basic_map_foreach_constraint(firing, print_guard_constraint, ti);
        if (ti->first)
                fprintf(ti->out, "true");
        fprintf(ti->out, ";\n");
        fprintf(ti->out, "\taction := ");
        for (int i = 0; i < n_out; ++i)
                fprintf(ti->out, "%sx%d' = x%d", i ? ", " : "", i, n_in + i);
        for (int i = n_out; i < ti->nvar; ++i)
                fprintf(ti->out, ", x%d' = 0", i);
        if (ti->lever)
                fprintf(ti->out, ", state' = %d", ti->to);
        fprintf(ti->out, ";\n");
        fprintf(ti->out, "};\n");

        ti->transition++;

        ti->state++;

        isl_basic_map_free(firing);

        return isl_stat_ok;
}

/* A firing relation that is a translation can be encoded directly
 * into a transition.  We simply put the translation in the action
 * and the remaining constraints in the guard.
 */
static isl_stat print_translation(__isl_take isl_basic_map *firing,
        transition_info *ti)
{
        fprintf(ti->out, "transition t%d := {\n", ti->transition);
        if (!ti->lever) {
                fprintf(ti->out, "\tfrom := s%d;\n", ti->from);
                fprintf(ti->out, "\tto := s%d;\n", ti->to);
        }

        fprintf(ti->out, "\tguard := ");
        if (!ti->lever)
                ti->first = 1;
        else {
                fprintf(ti->out, "state = %d", ti->from);
                ti->first = 0;
        }
        isl_basic_map_foreach_constraint(firing, print_guard_constraint, ti);
        if (ti->first)
                fprintf(ti->out, "true");
        fprintf(ti->out, ";\n");

        fprintf(ti->out, "\taction := ");
        if (!ti->lever)
                ti->first = 1;
        else {
                fprintf(ti->out, "state' = %d", ti->to);
                ti->first = 0;
        }
        ti->last_out = -1;
        isl_basic_map_foreach_constraint(firing, print_action_constraint, ti);
        for (int i = ti->last_out + 1; i < ti->nvar; ++i)
                fprintf(ti->out, "%sx%d' = 0",
                                    (i || !ti->first) ? ", " : "", i);
        fprintf(ti->out, ";\n");

        fprintf(ti->out, "};\n");

        ti->transition++;

        isl_basic_map_free(firing);

        return isl_stat_ok;
}

/* Print one or more transitions based on the firing relation.
 * The way a firing relation is encoded depends on whether it is
 * a translation or not.
 */
static isl_stat print_transition(__isl_take isl_basic_map *firing, void *user)
{
        transition_info *ti = (transition_info *)user;

        ti->translation = is_translation(firing);
        if (!ti->translation)
                return print_non_translation(firing, ti);
        return print_translation(firing, ti);
}

struct region_info {
        FILE *out;
        int name_offset;

        int first;
};

static isl_stat print_region_constraint(__isl_take isl_constraint *c,
        void *user)
{
        region_info *ri = (region_info *)user;
        isl_ctx *ctx;
        isl_printer *p;

        ctx = isl_constraint_get_ctx(c);
        p = isl_printer_to_file(ctx, ri->out);
        p = isl_printer_print_str(p, " && ");

        p = print_linear(p, c, 1, 0);
        if (isl_constraint_is_equality(c))
                p = isl_printer_print_str(p, " = 0");
        else
                p = isl_printer_print_str(p, " >= 0");

        isl_printer_free(p);
        isl_constraint_free(c);
        return isl_stat_ok;
}

/* Print a convex part of a region.
 * In FAST, states are designated by labels, while in Lever, they
 * are encoded as integer, so we have to drop the first letter ("s")
 * of the state name.
 */
static isl_stat print_region_basic_set(__isl_take isl_basic_set *bset,
        void *user)
{
        region_info *ri = (region_info *)user;

        if (!ri->first)
                fprintf(ri->out, " || ");

        fprintf(ri->out, "(");
        fprintf(ri->out, "state = %s",
                isl_basic_set_get_tuple_name(bset) + ri->name_offset);

        isl_basic_set_foreach_constraint(bset, &print_region_constraint, user);

        fprintf(ri->out, ")");

        ri->first = 0;

        isl_basic_set_free(bset);
        return isl_stat_ok;
}

static isl_stat print_region_set(__isl_take isl_set *set, void *user)
{
        isl_set_foreach_basic_set(set, &print_region_basic_set, user);

        isl_set_free(set);
        return isl_stat_ok;
}

static void print_region(FILE *out, const char *name,
        __isl_keep isl_union_set *region)
{
        region_info ri;

        ri.out = out;
        ri.first = 1;
        ri.name_offset = 0;

        fprintf(out, "Region %s := { ", name);
        isl_union_set_foreach_set(region, &print_region_set, &ri);
        fprintf(out, " };\n");
}

/* Construct a FAST model that tests whether all reachable
 * final states are a subset of the good final states.
 */
static int check_equivalence_fast(isl_ctx *ctx,
        dependence_graph *dg1, dependence_graph *dg2)
{
        int root;
        MSA msa;
        FILE *out = stdout;
        unsigned nparam;
        transition_info ti;
        isl_space *dim;
        isl_union_set *init;
        isl_union_set *final;
        isl_union_set *good;

        if (check_input_arrays(dg1, dg2))
                return -1;

        root = msa.add(dg1->out, dg2->out);

        ti.nvar = 0;
        ti.state = msa.state.size();
        ti.lever = 0;

        compute_max_var(msa, &ti);
        init = compute_equal_output(dg1, dg2);
        final = compute_final_states(ctx, msa);
        good = compute_good_states(ctx, msa);
        dim = isl_union_set_get_space(init);
        nparam = isl_space_dim(dim, isl_dim_param);

        fprintf(out, "model m1 {\n");

        fprintf(out, "var");
        for (int i = 0; i < ti.nvar; ++i)
                fprintf(out, "%s x%d", i ? "," : "", i);
        for (int i = 0; i < nparam; ++i)
                fprintf(out, ", %s", isl_space_get_dim_name(dim, isl_dim_param, i));
        fprintf(out, ";\n");

        fprintf(out, "states");
        for (int i = 0; i < ti.state; ++i)
                fprintf(out, "%s s%d", i ? "," : "", i);
        fprintf(out, ";\n");

        ti.out = out;
        ti.transition = 0;
        ti.state = msa.state.size();
        for (int i = 0; i < msa.transition.size(); ++i) {
                int f = msa.transition[i].second;
                ti.from = msa.transition[i].first.first;
                ti.to = msa.transition[i].first.second;
                isl_map_foreach_basic_map(msa.firing[f],
                        &print_transition, &ti);
        }

        fprintf(out, "}\n");

        fprintf(out, "strategy s1 {\n");
        print_region(out, "init", init);
        print_region(out, "final", final);
        print_region(out, "good", good);
        fprintf(out, "Transitions t := {");
        for (int i = 0; i < ti.transition; ++i)
                fprintf(out, "%s t%d", i ? "," : "", i);
        fprintf(out, " };\n");
        fprintf(out, "Region reach := post*(init, t);\n");
        fprintf(out, "if (subSet(reach && final, good))\n");
        fprintf(out, "\tthen print(\"OK\\n\");\n");
        fprintf(out, "\telse print(\"NOT OK\\n\");\n");
        fprintf(out, "endif\n");
        fprintf(out, "}\n");

        isl_space_free(dim);
        isl_union_set_free(init);
        isl_union_set_free(final);
        isl_union_set_free(good);

        return 0;
}

static void print_label(FILE *out, const char *name,
        __isl_keep isl_union_set *region)
{
        region_info ri;

        ri.out = out;
        ri.first = 1;
        ri.name_offset = 1;

        fprintf(out, "label %s := { ", name);
        isl_union_set_foreach_set(region, &print_region_set, &ri);
        fprintf(out, " };\n");
}

/* Construct a Lever model that tests whether the bad final states
 * are unreachable.
 */
static int check_equivalence_lever(isl_ctx *ctx,
        dependence_graph *dg1, dependence_graph *dg2)
{
        int root;
        MSA msa;
        FILE *out = stdout;
        unsigned nparam;
        transition_info ti;
        isl_space *dim;
        isl_union_set *init;
        isl_union_set *bad;

        if (check_input_arrays(dg1, dg2))
                return -1;

        root = msa.add(dg1->out, dg2->out);

        ti.nvar = 0;
        ti.state = msa.state.size();
        ti.lever = 1;

        compute_max_var(msa, &ti);
        init = compute_equal_output(dg1, dg2);
        bad = isl_union_set_subtract(compute_final_states(ctx, msa),
                                     compute_good_states(ctx, msa));
        dim = isl_union_set_get_space(init);
        nparam = isl_space_dim(dim, isl_dim_param);

        fprintf(out, "var state");
        for (int i = 0; i < ti.nvar; ++i)
                fprintf(out, ", x%d", i);
        for (int i = 0; i < nparam; ++i)
                fprintf(out, ", %s", isl_space_get_dim_name(dim, isl_dim_param, i));
        fprintf(out, ";\n");

        fprintf(out, "ctl_spec {\n");
        print_label(out, "bad", bad);
        fprintf(out, "}\n");

        fprintf(out, "system {\n");
        print_label(out, "init", init);

        ti.out = out;
        ti.transition = 0;
        ti.state = msa.state.size();
        for (int i = 0; i < msa.transition.size(); ++i) {
                int f = msa.transition[i].second;
                ti.from = msa.transition[i].first.first;
                ti.to = msa.transition[i].first.second;
                isl_map_foreach_basic_map(msa.firing[f],
                        &print_transition, &ti);
        }

        fprintf(out, "}\n");

        isl_space_free(dim);
        isl_union_set_free(init);
        isl_union_set_free(bad);

        return 0;
}

/* A firing relation that is not a translation is encoded in Aspic using
 * using two transitions passing through an extra state.
 * The first transition assigns arbitrary values to the extra variables
 * that will correspond to the output variables and the existentially
 * quantified variables of the relation.
 * The second transition, from the extra
 * state to the original target state, has a guard that corresponds
 * to the firing relation and an action that assigns the output
 * variables to the main variables.
 */
static isl_stat print_aspic_non_translation(__isl_take isl_basic_map *firing,
        transition_info *ti)
{
        unsigned n_in = isl_basic_map_dim(firing, isl_dim_in);
        unsigned n_out = isl_basic_map_dim(firing, isl_dim_out);
        unsigned n_div = isl_basic_map_dim(firing, isl_dim_div);

        fprintf(ti->out, "transition t%d := {\n", ti->transition);
        fprintf(ti->out, "\tfrom := s%d;\n", ti->from);
        fprintf(ti->out, "\tto := s%d;\n", ti->state);
        fprintf(ti->out, "\tguard := true;\n");
        fprintf(ti->out, "\taction := ");
        for (int i = 0; i < n_out + n_div; ++i)
                fprintf(ti->out, "%sx%d' = ?", i ? ", " : "", n_in + i);
        fprintf(ti->out, ";\n");

        fprintf(ti->out, "};\n");

        ti->transition++;

        fprintf(ti->out, "transition t%d := {\n", ti->transition);
        if (!ti->lever) {
                fprintf(ti->out, "\tfrom := s%d;\n", ti->state);
                fprintf(ti->out, "\tto := s%d;\n", ti->to);
        }
        fprintf(ti->out, "\tguard := ");
        if (!ti->lever)
                ti->first = 1;
        else {
                fprintf(ti->out, "state = %d", ti->state);
                ti->first = 0;
        }
        isl_basic_map_foreach_constraint(firing, print_guard_constraint, ti);
        if (ti->first)
                fprintf(ti->out, "true");
        fprintf(ti->out, ";\n");
        fprintf(ti->out, "\taction := ");
        for (int i = 0; i < n_out; ++i)
                fprintf(ti->out, "%sx%d' = x%d", i ? ", " : "", i, n_in + i);
        for (int i = n_out; i < ti->nvar; ++i)
                fprintf(ti->out, ", x%d' = 0", i);
        if (ti->lever)
                fprintf(ti->out, ", state' = %d", ti->to);
        fprintf(ti->out, ";\n");
        fprintf(ti->out, "};\n");

        ti->transition++;

        ti->state++;

        isl_basic_map_free(firing);

        return isl_stat_ok;
}

/* Print one or more transitions based on the firing relation.
 * The way a firing relation is encoded depends on whether it is
 * a translation or not.
 */
static isl_stat print_aspic_transition(__isl_take isl_basic_map *firing,
        void *user)
{
        transition_info *ti = (transition_info *)user;

        ti->translation = is_translation(firing);
        if (!ti->translation)
                return print_aspic_non_translation(firing, ti);
        return print_translation(firing, ti);
}

/* Print a transition from a bad final state to the canonical bad state.
 */
static isl_stat print_bad_aspic_transition(__isl_take isl_map *firing,
        void *user)
{
        transition_info *ti = (transition_info *)user;
        const char *name;

        name = isl_map_get_tuple_name(firing, isl_dim_in);
        ti->from = atoi(name + 1);

        isl_map_foreach_basic_map(firing, &print_aspic_transition, user);

        isl_map_free(firing);

        return isl_stat_ok;
}

/* Construct an Aspic model that tests whether the bad final states
 * are unreachable.  Since Aspic only supports convex "bad" regions,
 * we add an extra state with transitions from all instances of states
 * that are bad.  The bad region is then simply any instance of that
 * extra state.
 */
static int check_equivalence_aspic(isl_ctx *ctx,
        dependence_graph *dg1, dependence_graph *dg2)
{
        int root;
        MSA msa;
        FILE *out = stdout;
        unsigned nparam;
        transition_info ti;
        isl_space *dim;
        isl_union_set *init;
        isl_union_set *bad;
        isl_union_map *bad_map;

        if (check_input_arrays(dg1, dg2))
                return -1;

        root = msa.add(dg1->out, dg2->out);

        ti.nvar = 0;
        ti.state = msa.state.size();
        ti.lever = 0;

        compute_max_var(msa, &ti);
        init = compute_equal_output(dg1, dg2);
        bad = isl_union_set_subtract(compute_final_states(ctx, msa),
                                     compute_good_states(ctx, msa));
        dim = isl_union_set_get_space(init);
        nparam = isl_space_dim(dim, isl_dim_param);

        fprintf(out, "model m1 {\n");

        fprintf(out, "var");
        for (int i = 0; i < ti.nvar; ++i)
                fprintf(out, "%s x%d", i ? "," : "", i);
        for (int i = 0; i < nparam; ++i)
                fprintf(out, ", %s", isl_space_get_dim_name(dim, isl_dim_param, i));
        fprintf(out, ";\n");

        fprintf(out, "states");
        for (int i = 0; i < ti.state + 1; ++i)
                fprintf(out, "%s s%d", i ? "," : "", i);
        fprintf(out, ";\n");

        ti.out = out;
        ti.transition = 0;
        ti.state = msa.state.size();
        for (int i = 0; i < msa.transition.size(); ++i) {
                int f = msa.transition[i].second;
                ti.from = msa.transition[i].first.first;
                ti.to = msa.transition[i].first.second;
                isl_map_foreach_basic_map(msa.firing[f],
                        &print_aspic_transition, &ti);
        }

        bad_map = isl_union_map_from_domain_and_range(bad,
                        isl_union_set_from_set(isl_set_universe(dim)));
        ti.to = ti.state;
        isl_union_map_foreach_map(bad_map, &print_bad_aspic_transition, &ti);

        fprintf(out, "}\n");

        fprintf(out, "strategy s1 {\n");
        print_region(out, "init", init);
        fprintf(out, "Region bad := { state = s%d };\n", ti.state);
        fprintf(out, "}\n");

        isl_union_set_free(init);
        isl_union_map_free(bad_map);

        return 0;
}

int main(int argc, char *argv[])
{
        isl_ctx *ctx;
        struct options *options = options_new_with_defaults();
        isl_set *context = NULL;
        dependence_graph *dg[2];
        int res = 1;

        argc = options_parse(options, argc, argv, ISL_ARG_ALL);

        ctx = isl_ctx_alloc_with_options(&options_args, options);
        if (!ctx) {
                fprintf(stderr, "Unable to allocate ctx\n");
                return -1;
        }

        if (options->context)
                context = isl_set_read_from_str(ctx, options->context);

        pdg::PDG *pdg[2];
        unsigned out_dim = 0;
        for (int i = 0; i < 2; ++i) {
                FILE *in;
                in = fopen(options->program[i], "r");
                if (!in) {
                        fprintf(stderr, "Unable to open %s\n", options->program[i]);
                        return -1;
                }
                pdg[i] = yaml::Load<pdg::PDG>(in, ctx);
                fclose(in);
                if (!pdg[i]) {
                        fprintf(stderr, "Unable to read %s\n", options->program[i]);
                        return -1;
                }
                out_dim = update_out_dim(pdg[i], out_dim);

                if (context &&
                    pdg[i]->params.size() != isl_set_dim(context, isl_dim_param)) {
                        fprintf(stderr,
                            "Parameter dimension mismatch; context ignored\n");
                        isl_set_free(context);
                        context = NULL;
                }
        }

        for (int i = 0; i < 2; ++i) {
                dg[i] = pdg_to_dg(pdg[i], out_dim, isl_set_copy(context));
                pdg[i]->free();
                delete pdg[i];
        }

        if (options->reachability == REACH_TC)
                res = check_equivalence_tc(ctx, dg[0], dg[1]);
        else if (options->reachability == REACH_FAST)
                res = check_equivalence_fast(ctx, dg[0], dg[1]);
        else if (options->reachability == REACH_LEVER)
                res = check_equivalence_lever(ctx, dg[0], dg[1]);
        else if (options->reachability == REACH_ASPIC)
                res = check_equivalence_aspic(ctx, dg[0], dg[1]);

        for (int i = 0; i < 2; ++i)
                delete dg[i];
        isl_set_free(context);
        isl_ctx_free(ctx);

        return res;
}