#include <stdio.h>
#include <assert.h>

#include "caml/memory.h"
#include "caml/mlvalues.h"
#include "caml/alloc.h"

#include "alpha.h"

#ifdef __STDC__
#define inline
#endif

#ifdef __GNUC__
#undef inline
#define inline __inline__
#endif

/***************************************************************************/

void failwith(char *);

static inline void sanity_check(value x)
{
    assert(!(-1024 < x && x < 1024));
}

static inline int word_aligned_p(void *p)
{
    unsigned x = (unsigned)p;
#ifdef ARCH_SIXTYFOUR
    return x % 8 == 0;
#else
    return x % 4 == 0;
#endif
}

static inline int safe_Int_val(value v)
{
    assert(Is_long(v));
    return Int_val(v);
}

static inline value safe_ptr(value v)
{
    assert(Is_block(v));
    return v;
}

/************************************************************************
 * Arenas
 ************************************************************************/

static int num_of_arenas = 0;

typedef struct {
    int num;
    arena_info *info;
    int count;
} ma_arena_obj;

static ma_arena_obj *ma_create_arena()
{
    ma_arena_obj *ao = malloc(sizeof(ma_arena_obj));
    assert(word_aligned_p(ao));
    ao->num = num_of_arenas++;
    ao->info = a_arena_make();
    assert(ao->info);
    ao->count = 0;
    return ao;
}

typedef struct {
    final_fun f; 
    ma_arena_obj *arena;
    void *val;
} ma_ptr_obj;

#define Ptr_arena(x) (((ma_ptr_obj *)(x))->arena)
#define Ptr_val(x) (((ma_ptr_obj *)(x))->val)

static void ma_free_ptr(value v);

#define WOSIZE(T) (sizeof(T) / sizeof(value))

value Val_alpha_ptr(ma_arena_obj *ao, void *p)
{
    value result = Val_unit;
    Begin_roots1(result);
    result = alloc_shr(WOSIZE(ma_ptr_obj), Final_tag);
    Final_fun(result)= ma_free_ptr;
    assert(ao->info != 0);
    Ptr_arena(result) = ao; ao->count++;
    Ptr_val(result) = p;
    End_roots();
    return result;
}

static inline void validate_arena_ptr(value v)
{
    (void)v;
    assert(Ptr_arena(v) != 0);
    assert(Ptr_arena(v)->info != 0);
}

static void ma_free_ptr(value v)
{
    ma_arena_obj *ao;
    validate_arena_ptr(v);
    ao = Ptr_arena(v);
    Ptr_arena(v) = 0;
    assert(ao->count > 0);
    ao->count--;
    if (ao->count == 0) {
	a_arena_free(ao->info);
	ao->info = 0;
#ifdef NDEBUG
	free(ao);
#endif
    }
}

value ma_make_arena(value unit)
{
    ma_arena_obj *ao = ma_create_arena();
    return Val_alpha_ptr(ao,ao);
}

ma_arena_obj *Arena_val(value x)
{
    ma_arena_obj *ao;
    safe_ptr(x);
    validate_arena_ptr(x);
    ao = (ma_arena_obj*)Ptr_val(x);
    assert(ao == Ptr_arena(x));
    assert(ao->info != 0);
    return ao;
}

/************************************************************************
 * Integers
 ************************************************************************/

value ma_gcd(value zx, value zy)
{
    int x = safe_Int_val(zx);
    int y = safe_Int_val(zy);
    int res = a_gcd(x, y);
    return Val_int(res);
}

value ma_egcd(value za, value zb)
{
    int a, b, u, v, d;
    value result = Val_unit;
    Begin_roots3(result, za, zb);
    a = safe_Int_val(za);
    b = safe_Int_val(zb);
    d = a_egcd(a,b,&u,&v);
    result = alloc(3, 0);
    sanity_check(result);
    Field(result, 0) = Val_int(d);
    Field(result, 1) = Val_int(u);
    Field(result, 2) = Val_int(v);
    End_roots();
    return result;
}

value ma_lcm(value zx, value zy)
{
    int x = safe_Int_val(zx);
    int y = safe_Int_val(zy);
    int res = a_lcm(x, y);
    return Val_int(res);
}

/************************************************************************
 * Vectors
 ************************************************************************/

static inline AlphaVector AlphaVector_val(value x)
{
    validate_arena_ptr(x);
    return ((AlphaVector)(Ptr_val(safe_ptr(x))));
}

value ma_make_vector(value zao, value n)
{
    ma_arena_obj *ao = Arena_val(zao);
    return Val_alpha_ptr(ao, a_make_vector(ao->info, safe_Int_val(n)));
}

value ma_length_vec(value zv)
{
    AlphaVector v = AlphaVector_val(zv);
    return Val_int(a_length_vec(v));
}

value ma_elt_vec_ref(value zv, value zi)
{
    AlphaVector v = AlphaVector_val(zv);
    int i = safe_Int_val(zi);
    int x = a_elt_vec(v,i);
    return Val_int(x);
}

value ma_elt_vec_set(value zv, value zi, value zx)
{
    AlphaVector v = AlphaVector_val(zv);
    int i = safe_Int_val(zi);
    int x = safe_Int_val(zx);
    a_elt_vec(v,i) = x;
    return Val_unit;
}

/************************************************************************
 * Matrixs
 ************************************************************************/

static inline AlphaMatrix AlphaMatrix_val(value x)
{
    validate_arena_ptr(x);
    return ((AlphaMatrix)(Ptr_val(safe_ptr(x))));
}

value ma_make_matrix(value zao, value zrs, value zcs)
{
    ma_arena_obj *ao = Arena_val(zao);
    int rs = safe_Int_val(zrs);
    int cs = safe_Int_val(zcs);
    return Val_alpha_ptr(ao, a_make_matrix(ao->info, rs, cs));
}

value ma_rows_mx(value m)
{
    return Val_int(a_rows_mx(AlphaMatrix_val(m)));
}

value ma_cols_mx(value m)
{
    return Val_int(a_cols_mx(AlphaMatrix_val(m)));
}

value ma_elt_mx_ref(value zm, value zi, value zj)
{
    AlphaMatrix m = AlphaMatrix_val(zm);
    int i = safe_Int_val(zi);
    int j = safe_Int_val(zj);
    int x = a_elt_mx(m,i,j);
    return Val_int(x);
}


value ma_elt_mx_set(value zm, value zi, value zj, value zx)
{
    AlphaMatrix m = AlphaMatrix_val(zm);
    int i = safe_Int_val(zi);
    int j = safe_Int_val(zj);
    int x = safe_Int_val(zx);
    a_elt_mx(m,i,j) = x;
    return Val_unit;
}

/************************************************************************
 * Routines
 ************************************************************************/

value ma_set_vec(value zv, value zx)
{
    AlphaVector v = AlphaVector_val(zv);
    int x = safe_Int_val(zx);
    a_set_vec(v,x);
    return Val_unit;
}

value ma_neg_vec(value zv)
{
    AlphaVector v = AlphaVector_val(zv);
    a_neg_vec(v);
    return Val_unit;
}

value ma_saxpy(value za, value zx, value zb, value zy, value zz)
{
    int a = safe_Int_val(za);
    AlphaVector x = AlphaVector_val(zx);
    int b = safe_Int_val(zb);
    AlphaVector y = AlphaVector_val(zy);
    AlphaVector z = AlphaVector_val(zz);
    a_saxpy(a,x,b,y,z);
    return Val_unit;
}

value ma_copy_vec(value zdst, value zsrc)
{
    AlphaVector dst = AlphaVector_val(zdst);
    AlphaVector src = AlphaVector_val(zsrc);
    a_copy_vec(dst,src);
    return Val_unit;
}

value ma_resize_vec(value zao, value zx, value znewn)
{
    ma_arena_obj *ao = Arena_val(zao);
    AlphaVector x = AlphaVector_val(zx);
    int newn = safe_Int_val(znewn);
    a_resize_vec(ao->info,&x,newn);
    return Val_alpha_ptr(ao, x);
}

value ma_is_zero_vec(value zx)
{
    AlphaVector x = AlphaVector_val(zx);
    int res = a_is_zero_vec(x);
    return res ? Val_true : Val_false;
}

value ma_gcd_vec(value zv)
{
    int res;
    AlphaVector v = AlphaVector_val(zv);
    res = a_gcd_vec(v);
    return Val_int(res);
}

value ma_extract_vec(value zx, value zi1, value zi2, value zy)
{
    AlphaVector x = AlphaVector_val(zx);
    int i1 = safe_Int_val(zi1);
    int i2 = safe_Int_val(zi2);
    AlphaVector y = AlphaVector_val(zy);
    a_extract_vec(x,i1,i2,y);
    return Val_unit;
}

value ma_copy_mx(value zdst, value zsrc)
{
    AlphaMatrix dst = AlphaMatrix_val(zdst);
    AlphaMatrix src = AlphaMatrix_val(zsrc);
    a_copy_mx(dst,src);
    return Val_unit;
}
    
value ma_swap_col(value zA, value zi, value zj)
{
    AlphaMatrix A = AlphaMatrix_val(zA);
    int i = safe_Int_val(zi);
    int j = safe_Int_val(zj);
    a_swap_col(A, i, j);
    return Val_unit;
}

value ma_swap_row(value zA, value zi, value zj)
{
    AlphaMatrix A = AlphaMatrix_val(zA);
    int i = safe_Int_val(zi);
    int j = safe_Int_val(zj);
    a_swap_row(A, i, j);
    return Val_unit;
}

value ma_set2ident(value zA)
{
    AlphaMatrix A = AlphaMatrix_val(zA);
    a_set2ident(A);
    return Val_unit;
}

value ma_is_zero_matrix (value zA)
{
    AlphaMatrix A = AlphaMatrix_val(zA);
    int b = a_is_zero_matrix(A);
    return b ? Val_true : Val_false;
}

value ma_is_ident_matrix (value zA)
{
    AlphaMatrix A = AlphaMatrix_val(zA);
    int b = a_is_ident_matrix(A);
    return b ? Val_true : Val_false;
}

value ma_matmult(value zao, value zA, value zB, value zC)
{
    ma_arena_obj *ao = Arena_val(zao);
    AlphaMatrix A = AlphaMatrix_val(zA);
    AlphaMatrix B = AlphaMatrix_val(zB);
    AlphaMatrix C = AlphaMatrix_val(zC);
    a_matmult(ao->info, A, B, C);
    return Val_unit;
}

value ma_mvm(value zao, value zA, value zx, value zy)
{
   ma_arena_obj *ao = Arena_val(zao);
   AlphaMatrix A = AlphaMatrix_val(zA);
   AlphaVector x = AlphaVector_val(zx);
   AlphaVector y = AlphaVector_val(zy);
   a_mvm(ao->info, A, x, y);
   return Val_unit;
}

value ma_transpose(value zA, value zAt)
{
    AlphaMatrix A = AlphaMatrix_val(zA);
    AlphaMatrix At = AlphaMatrix_val(zAt);
    a_transpose(A, At);
    return Val_unit;
}

value ma_extract_native(value zA, value zi1, value zi2, value zj1, value zj2,
			value zB)
{
    AlphaMatrix A;
    int i1, i2, j1, j2;
    AlphaMatrix B;

    /* A i1 i2 j1 j2 B */
    A = AlphaMatrix_val(zA);
    i1 = safe_Int_val(zi1);
    i2 = safe_Int_val(zi2);
    j1 = safe_Int_val(zj1);
    j2 = safe_Int_val(zj2);
    B = AlphaMatrix_val(zB);

    a_extract(A,i1,i2,j1,j2,B);
    return Val_unit;
}

value ma_extract(value *argv, int argn)
{
    assert(argn == 6);
    return ma_extract_native(argv[0],argv[1],argv[2],argv[3],argv[4],
			     argv[5]);
}

value ma_replace_native(value zA, value zi1, value zi2, value zj1, value zj2,
			value zB)
{
    AlphaMatrix A;
    int i1, i2, j1, j2;
    AlphaMatrix B;

    /* A i1 i2 j1 j2 B */
    A = AlphaMatrix_val(zA);
    i1 = safe_Int_val(zi1);
    i2 = safe_Int_val(zi2);
    j1 = safe_Int_val(zj1);
    j2 = safe_Int_val(zj2);
    B = AlphaMatrix_val(zB);

    a_replace(A,i1,i2,j1,j2,B);
    return Val_unit;
}

value ma_replace(value *argv, int argn)
{
    assert(argn == 6);
    return ma_replace_native(argv[0],argv[1],argv[2],argv[3],argv[4],
			     argv[5]);
}


value ma_add_cols(value zao, value zA, value zc)
{
    ma_arena_obj *ao = Arena_val(zao);
    AlphaMatrix A = AlphaMatrix_val(zA);
    int c = safe_Int_val(zc);
    a_add_cols(ao->info,&A,c);
    return Val_alpha_ptr(ao, A);
}

value ma_trunc_cols(value zao, value zA, value zc)
{
    ma_arena_obj *ao = Arena_val(zao);
    AlphaMatrix A = AlphaMatrix_val(zA);
    int c = safe_Int_val(zc);
    a_trunc_cols(ao->info,&A,c);
    return Val_alpha_ptr(ao, A);
}

value ma_lower_echelon(value zao, value zA, value zU)
{
    ma_arena_obj *ao = Arena_val(zao);
    AlphaMatrix A = AlphaMatrix_val(zA);
    AlphaMatrix U = AlphaMatrix_val(zU);
    int prank = 0;
    a_lower_echelon(ao->info, A, U, &prank);
    return Val_int(prank);
}

value ma_snf(value zao, value zA, value zU, value zV)
{
    ma_arena_obj *ao = Arena_val(zao);
    AlphaMatrix A = AlphaMatrix_val(zA);
    AlphaMatrix U = AlphaMatrix_val(zU);
    AlphaMatrix V = AlphaMatrix_val(zV);
    int prank = 0;
    a_snf(ao->info, A, U, V, &prank);
    return Val_int(prank);
}

value ma_solve_linear(value zao, value zA, value zB, value zX0, value zH)
{
    /* returns -1 is there are no solutions */

    ma_arena_obj *ao = Arena_val(zao);
    AlphaMatrix A = AlphaMatrix_val(zA);
    AlphaVector B = AlphaVector_val(zB);
    AlphaVector X0 = AlphaVector_val(zX0);
    AlphaMatrix H = AlphaMatrix_val(zH);
    int rank = 0;
    int has_solutions;
    has_solutions = a_solve_linear(ao->info, A, B, X0, H, &rank);
    rank = has_solutions ? rank : -1;
    return Val_int(rank);
}

value ma_solve_echelon(value zL, value zB, value zX)
{
    AlphaMatrix L = AlphaMatrix_val(zL);
    AlphaVector B = AlphaVector_val(zB);
    AlphaVector X = AlphaVector_val(zX);
    int succ = a_solve_echelon(L, B, X);
    return succ ? Val_true : Val_false;
}

value ma_unit_backsolve (value zU, value zb, value zx)
{
    AlphaMatrix U = AlphaMatrix_val (zU);
    AlphaVector b = AlphaVector_val (zb);
    AlphaVector x = AlphaVector_val (zx);
    a_unit_backsolve (U, b, x);
    return Val_unit;
}

value ma_block_unit_backsolve (value zU, value zB, value zX)
{
    AlphaMatrix U = AlphaMatrix_val (zU);
    AlphaMatrix B = AlphaMatrix_val (zB);
    AlphaMatrix X = AlphaMatrix_val (zX);
    a_block_unit_backsolve (U, B, X);
    return Val_unit;
}


/************************************************************************/
/************************************************************************/
/************************************************************************/
/***			      Fourier-Motzkin                         ***/
/************************************************************************/
/************************************************************************/
/************************************************************************/

static xa_calc_set ml_set_to_alpha_set(arena_info *arena, value set_v);
static value alpha_set_to_ml_set(xa_calc_set set);

value ma_simplify_set(value in_set_v)
{
    arena_info *arena = a_arena_make();
    xa_calc_set in_set;
    xa_calc_set out_set;
    value out_set_v;

    in_set = ml_set_to_alpha_set(arena, in_set_v);
    if (in_set != 0)
	out_set = a_simplify_set(arena, in_set);
    else
	out_set = 0;
    out_set_v = alpha_set_to_ml_set(out_set);
    a_arena_free(arena);

    return out_set_v;
}

value ma_satisfiable_p(value in_set_v)
{
    arena_info *arena = a_arena_make();
    xa_calc_set in_set = ml_set_to_alpha_set(arena, in_set_v);
    int result;
    if (in_set != 0)
	result = a_satisfiable_p(arena, in_set);
    else
	result = 0;
    a_arena_free(arena);
    return result ? Val_true : Val_false;
}

/************************************************************************/
/*				 Accessors                              */
/************************************************************************/

static inline value *field_ref(value x, int tag, int size, int i)
{
    assert(Is_block(x));
    assert(Tag_val(x) == tag);
    assert(Wosize_val(x) == size);
    assert(0 <= i && i < size);
    return &Field(x, i);
}

/* List accessors */
static inline int nullp(value x) { return x == Val_int(0); }
static inline int consp(value x) { return Is_block(x); }
static inline value *car(value x) { return field_ref(x, 0, 2, 0); }
static inline value *cdr(value x) { return field_ref(x, 0, 2, 1); }
static inline value null() { return Val_int(0); }
static inline value cons(value CAR, value CDR)
{
    value result = Val_unit;
    Begin_roots3(result, CAR, CDR);
    result = alloc(2,0);
    sanity_check(result);
    *car(result) = CAR;
    *cdr(result) = CDR;
    End_roots();
    return result;
}

/* fmve_set_type accessors */
static inline int no_solution_p(value x) { return x == Val_int(0); }
static inline value *set_outs(value x)   { return field_ref(x, 0, 3, 0); }
static inline value *set_exists(value x) { return field_ref(x, 0, 3, 1); }
static inline value *set_pred(value x)   { return field_ref(x, 0, 3, 2); }
static inline value make_no_solution()
{
    return Val_int(0);
}
static inline value make_set(value outs, value exists, value pred)
{
    value result = Val_unit;
    Begin_roots4(result,outs,exists,pred);
    result = alloc(3,0);
    sanity_check(result);
    *set_outs(result) = outs;
    *set_exists(result) = exists;
    *set_pred(result) = pred;
    End_roots();
    return result;
}

/* fmve_constraint_type accessors */
static inline value *constr_kind(value x)  { return field_ref(x, 0, 3, 0); }
static inline value *constr_left(value x)  { return field_ref(x, 0, 3, 1); }
static inline value *constr_right(value x) { return field_ref(x, 0, 3, 2); }
static inline value make_constr(xa_calc_op k, value left, value right)
{
    value result = Val_unit, kind = Val_unit;
    Begin_roots4(result,kind,left,right);

    switch (k) {
    case XA_CALC_EQ: kind = Val_int(0); break;
    case XA_CALC_LE: kind = Val_int(1); break;
    case XA_CALC_LT: kind = Val_int(2); break;
    case XA_CALC_GE: kind = Val_int(3); break;
    case XA_CALC_GT: kind = Val_int(4); break;
    default: assert(0);
    }

    result = alloc(3,0);
    sanity_check(result);
    *constr_kind(result) = kind;
    *constr_left(result) = left;
    *constr_right(result) = right;

    End_roots();
    return result;
}

static inline xa_calc_op xlate_constr_kind(value k)
{
    /* These must occur in the same order as in alpha.mli */
    switch (Int_val(k)) {
    case 0: return XA_CALC_EQ;
    case 1: return XA_CALC_LE;
    case 2: return XA_CALC_LT;
    case 3: return XA_CALC_GE;
    case 4: return XA_CALC_GT;
    default: assert(0); return XA_CALC_UNKNOWN;
    }
}

static inline int term_const_p(value x) { return Tag_val(x) == 0; }
static inline value *term_const_val(value x) { return field_ref(x, 0, 1, 0); }
static inline value *term_var_coef(value x) { return field_ref(x, 1, 2, 0); }
static inline value *term_var_name(value x)  { return field_ref(x, 1, 2, 1); }
static inline value make_term_const(int val)
{
    value result = Val_unit;
    Begin_roots1(result);
    result = alloc(1, 0);
    sanity_check(result);
    *term_const_val(result) = Val_int(val);
    End_roots();
    return result;
}
static inline value make_term_var(int coef, char *var)
{
    value result = Val_unit, var_val = Val_unit;
    Begin_roots2(result, var_val);
    var_val = copy_string(var);
    result = alloc(2, 1);
    sanity_check(result);
    *term_var_coef(result) = Val_int(coef);
    *term_var_name(result) = var_val;
    End_roots();
    return result;
}

/************************************************************************/
/*			    ML Set -> Alpha Set                         */
/************************************************************************/

#define DEFINE_ML_LIST_WALKER(Name, Type, Subname)		\
static A_CONS_CELL(Type) Name(arena_info *arena, value lst)	\
{								\
    if (nullp(lst))						\
	return 0;						\
    else {							\
	A_CONS_CELL(Type) p;					\
	a_make_cons(arena, Type, p, Subname(arena, *car(lst)),	\
		    Name(arena, *cdr(lst)));			\
	return p;						\
    }								\
}								\
static A_CONS_CELL(Type) Name(arena_info *arena, value lst)


static A_CONS_CELL(xa_calc_var) walk_ml_vars(arena_info *arena, value lst);
static xa_calc_var walk_ml_var(arena_info *arena, value v);
static A_CONS_CELL(xa_calc_constraint) walk_ml_constrs(arena_info *arena, 
						       value lst);
static xa_calc_constraint walk_ml_constr(arena_info *arena, value v);
static A_CONS_CELL(xa_calc_term) walk_ml_terms(arena_info *arena, value lst);
static xa_calc_term walk_ml_term(arena_info *arena, value v);


static xa_calc_set ml_set_to_alpha_set(arena_info *arena, value set_v)
{
    if (no_solution_p(set_v))
	return 0;
    else
	return a_make_calc_set(
	    arena, walk_ml_vars(arena, *set_outs(set_v)),
	    walk_ml_vars(arena, *set_exists(set_v)),
	    walk_ml_constrs(arena, *set_pred(set_v)));
}

DEFINE_ML_LIST_WALKER(walk_ml_vars, xa_calc_var, walk_ml_var);

static xa_calc_var walk_ml_var(arena_info *arena, value v)
{
    return a_make_calc_var(arena, String_val(v));
}

DEFINE_ML_LIST_WALKER(walk_ml_constrs, xa_calc_constraint, walk_ml_constr);

static xa_calc_constraint walk_ml_constr(arena_info *arena, value v)
{
    return a_make_calc_constraint(arena,
				  xlate_constr_kind(*constr_kind(v)),
				  walk_ml_terms(arena, *constr_left(v)),
				  walk_ml_terms(arena, *constr_right(v)));
}

DEFINE_ML_LIST_WALKER(walk_ml_terms, xa_calc_term, walk_ml_term);

static xa_calc_term walk_ml_term(arena_info *arena, value v)
{
    if (term_const_p(v))
	return a_make_calc_term(arena, 0, safe_Int_val(*term_const_val(v)));
    else
	return a_make_calc_term(
	    arena, walk_ml_var(arena, *term_var_name(v)),
	    safe_Int_val(*term_var_coef(v)));
}


/************************************************************************/
/*			    ML Set -> Alpha Set                         */
/************************************************************************/

#define DEFINE_ALPHA_LIST_WALKER(Name, Type, Subname)		\
static value Name(A_CONS_CELL(Type) lst)			\
{								\
    if (lst == 0)						\
	return null();						\
    else {							\
        value result = Val_unit, CAR = Val_unit, CDR = Val_unit; \
	Begin_roots3(result,CAR,CDR);				\
        CAR = Subname(a_car(lst));				\
        CDR = Name(a_cdr(lst));					\
	result = cons(CAR,CDR);					\
	End_roots();						\
	return result;						\
    }								\
}								\
static value Name(A_CONS_CELL(Type) lst)

static value walk_alpha_vars(A_CONS_CELL(xa_calc_var) lst);
static value walk_alpha_var(xa_calc_var v);
static value walk_alpha_constrs(A_CONS_CELL(xa_calc_constraint) lst);
static value walk_alpha_constr(xa_calc_constraint v);
static value walk_alpha_terms(A_CONS_CELL(xa_calc_term) lst);
static value walk_alpha_term(xa_calc_term v);

static value alpha_set_to_ml_set(xa_calc_set set)
{
    if (set == 0)
	return make_no_solution();
    else {
	value result = Val_unit, set_vars = Val_unit, 
	    exist_vars = Val_unit, pred = Val_unit;
	Begin_roots4(result,set_vars,exist_vars,pred);
	set_vars = walk_alpha_vars(xa_set_vars_calc_set(set));
	exist_vars = walk_alpha_vars(xa_exist_vars_calc_set(set));
	pred = walk_alpha_constrs(xa_pred_calc_set(set));
	result = make_set(set_vars, exist_vars, pred);
	End_roots();
	return result;
    }
}

DEFINE_ALPHA_LIST_WALKER(walk_alpha_vars, xa_calc_var, walk_alpha_var);

static value walk_alpha_var(xa_calc_var v)
{
    return copy_string(xa_name_calc_var(v));
}

DEFINE_ALPHA_LIST_WALKER(walk_alpha_constrs, xa_calc_constraint,
			 walk_alpha_constr);

static value walk_alpha_constr(xa_calc_constraint v)
{
    value result = Val_unit, t1 = Val_unit, t2 = Val_unit;
    Begin_roots3(result, t1, t2);
    t1 = walk_alpha_terms(xa_t1_calc_constraint(v));
    t2 = walk_alpha_terms(xa_t2_calc_constraint(v));
    result = make_constr(xa_op_calc_constraint(v), t1, t2);
    End_roots();
    return result;
}

DEFINE_ALPHA_LIST_WALKER(walk_alpha_terms, xa_calc_term, walk_alpha_term);

static value walk_alpha_term(xa_calc_term v)
{
    if (xa_is_const_calc_term(v))
	return make_term_const(xa_coeff_calc_term(v));
    else
	return make_term_var(
	    xa_coeff_calc_term(v),
	    xa_name_calc_var(xa_var_calc_term(v)));
}

