/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  Presprover -- Prove formulas of Presburger arithmetic
  Copyright (C) 2005 Markus Triska triska@gmx.at

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

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

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

  http://stud4.tuwien.ac.at/~e0225855/presprover/presprover.html


  Cf.: Constraint Solving on Terms, H. Comon and C. Kirchner
  published in: "Constraints in Computational Logics", Springer 2001
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


:- use_module(library(bounds)).

:- op(745, xfx, <=).
:- op(750, xfy, /\).
:- op(751, xfy, \/).
:- op(760, xfy, =>).


teq(eq([m(1,x),m(2,y),m(-3,z)],1)).  % x + 2*y - 3z = 1
teq(eq([m(1,x),m(2,y)],3)).
teq(eq([m(2,x),m(5,y)],20)).

tineq(ineq([m(2,x),m(-1,y)],-1)).

append_without([], [], _Without, []).
append_without([], [Y|Ys], Without, As) :-
	( memberchk(Y, Without) ->
		append_without([], Ys, Without, As)
	;
		As = [Y|Rs],
		append_without([], Ys, Without, Rs)
	).
append_without([X|Xs], Ys, Without, [X|As]) :-
	append_without(Xs, Ys, Without, As).

:- append_without([1,2,3], [7,8], [7], [1,2,3,8]).


automaton_empty(aut(_Q,QF,Q0,Delta)) :-
	\+ find_halting(Q0, [], _Visited, QF, Delta).

% for testing
automaton_haltingpath(A, Path) :-
	A = aut(_,QF,Q0,Delta),
	find_halting(Q0, [], Path, QF, Delta).


find_halting(Current, Visited, Visited, QF, _Delta) :-
	memberchk(Current, QF).
find_halting(Current, Visited0, Visited, QF, Delta) :-
	member(delta(Current,_String,Next), Delta),
	append(Visited0, [Current], Visited1),
	\+ memberchk(Next, Visited1),
	find_halting(Next, Visited1, Visited, QF, Delta).

linsum(Ms, S) :- linsum(Ms, 0, S).

linsum([], S, S).
linsum([m(F,V)|Ms], Acc0, S) :-
	Acc1 #= Acc0 + F*V,
	linsum(Ms, Acc1, S).


linsum_withcoeffs(Ms, Coeff, S) :-
	linsum_withcoeffs(Ms, Coeff, 0, S).

linsum_withcoeffs([], [], S, S).
linsum_withcoeffs([m(Factor,_Var)|Ms], [Coeff|Cs], Acc0, S) :-
	Acc1 is Acc0 + (Factor*Coeff),
	linsum_withcoeffs(Ms, Cs, Acc1, S).

mterm_var(m(_C,Var), Var).

states_tuples_transitions([], _Q, [], []).
states_tuples_transitions([S|Ss], Q, [T|Ts], [D|Ds]) :-
	D = delta(Q,T,S),
	states_tuples_transitions(Ss, Q, Ts, Ds).


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
			       Equality
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


factor_mod2(m(C0,_), m(C,Var)) :-
	C is abs(C0 mod 2),
	Var in 0..1.


eq_mod2(eq(Ms,C0), Tuples) :-
	C is abs(C0 mod 2),
	maplist(factor_mod2, Ms, Xs),
	maplist(mterm_var, Xs, Vs),
	linsum(Xs, S),
	findall(Vs, (S mod 2 #= C, label(Vs)), Tuples).

negative(A, B) :- B is (-1)*A.

eq_tuple_newstate(Ms, C, Tuple, Newstate) :-
	maplist(negative, Tuple, Negtuple),
	linsum_withcoeffs(Ms, Negtuple, Sum),
	D is (Sum + C) / 2,
	Newstate = q(D).


saturate_eq([], _Eq, Q, Q, D, D).
saturate_eq([q(C)|QIterRest], Eq, Q0, Q, Delta0, Delta) :-
	Eq = eq(Ms,_),
	eq_mod2(eq(Ms,C), Tuples),
	maplist(eq_tuple_newstate(Ms,C), Tuples, Newstates),
	states_tuples_transitions(Newstates, q(C), Tuples, Newdeltas),
	append(Delta0, Newdeltas, Delta1),
	append_without(QIterRest, Newstates, Q0, QIter1),
	append(Q0, Newstates, Q1),
	saturate_eq(QIter1, Eq, Q1, Q, Delta1, Delta).


eq_automaton(Eq, A) :-
	Eq = eq(_,B),
	Q0 = [q(B)],
	Delta0 = [],
	saturate_eq(Q0, Eq, Q0, Q, Delta0, Delta),
	sort(Q, Q1),
	sort(Delta, Delta1),
	( memberchk(q(0), Q) ->
		QF = [q(0)]
	;
		QF = []
	),
	A = aut(Q1, QF, q(B), Delta1).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
			      Inequality
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

stategeq0(q(N)) :- N >= 0.

ineq_automaton(Ineq, A) :-
	Ineq = ineq(Ms,B),
	Q0 = [q(B)],
	Delta0 = [],
	length(Ms, LMs),
	length(Thetas, LMs),
	Thetas in 0..1,
	findall(Thetas, label(Thetas), Tuples),
	saturate_ineq(Q0, Ineq, Tuples, Q0, Q, Delta0, Delta),
	sort(Q, Q1),
	sort(Delta, Delta1),
	sublist(stategeq0, Q1, QF),
	A = aut(Q1, QF, q(B), Delta1).

ineq_tuple_newstate(Ms, C, Tuple, State) :-
	linsum_withcoeffs(Ms, Tuple, Sum),
	D is floor((C - Sum) / 2),
	State = q(D).

saturate_ineq([], _Ineq, _Tuples, Q, Q, Delta, Delta).
saturate_ineq([q(C)|QIterRest], Ineq, Tuples, Q0, Q, Delta0, Delta) :-
	Ineq = ineq(Ms,_),
	maplist(ineq_tuple_newstate(Ms,C), Tuples, Newstates),
	states_tuples_transitions(Newstates, q(C), Tuples, Newdeltas),
	append(Delta0, Newdeltas, Delta1),
	append_without(QIterRest, Newstates, Q0, QIter1),
	append(Q0, Newstates, Q1),
	saturate_ineq(QIter1, Ineq, Tuples, Q1, Q, Delta1, Delta).


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		     Intersection of two automata
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


uncurry(A, B, (A,B)).

pairs([], _Ys, []).
pairs([E|Es], Ys, Pairs) :-
	maplist(uncurry(E), Ys, EPairs),
	append(EPairs, Restpairs, Pairs),
	pairs(Es, Ys, Restpairs).

:- pairs([1,2], [a,b], [(1,a),(1,b),(2,a),(2,b)]).

pair_state((A,B), q(i(A,B))).

statepairs(As, Bs, Qs) :-
	pairs(As, Bs, Qs1),
	maplist(pair_state, Qs1, Qs).

:- statepairs([q(1),q(0)],[q(2)], [q(i(q(1),q(2))), q(i(q(0),q(2)))]).

intersec_ishalting(QF1, QF2, Q) :-
	Q = q(i(A,B)),
	memberchk(A, QF1),
	memberchk(B, QF2).

:- intersec_ishalting([q(0),q(1)],[q(0)], q(i(q(1),q(0)))).

transform_delta(Q, S-E0, delta(Q,S,E)) :-
	pair_state(E0, E).

intersec_delta([], _Delta1, _Delta2, []).
intersec_delta([Q|Qs], Delta1, Delta2, Ds) :-
	Q = q(i(Q1,Q2)),
	findall(S-(E1,E2), (member(delta(Q1,S,E1),Delta1),member(delta(Q2,S,E2),Delta2)), Ns),
	maplist(transform_delta(Q), Ns, Newdeltas),
	append(Newdeltas, Restdeltas, Ds),
	intersec_delta(Qs, Delta1, Delta2, Restdeltas).


aut_intersection(A1, A2, I) :-
	A1 = aut(Q1,QF1,S1,Delta1),
	A2 = aut(Q2,QF2,S2,Delta2),
	statepairs(Q1, Q2, Qs),
	intersec_delta(Qs, Delta1, Delta2, Delta),
	sublist(intersec_ishalting(QF1,QF2), Qs, QF),
	pair_state((S1,S2), IntStart),
	I = aut(Qs,QF,IntStart,Delta).


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
			Union of two automata
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


state_rename(N, P, n(N,P)).

delta_rename([], _N, []).
delta_rename([D|Ds], N, [R|Rs]) :-
	D = delta(A0,Trans,B0),
	state_rename(N, A0, A),
	state_rename(N, B0, B),
	R = delta(A, Trans, B),
	delta_rename(Ds, N, Rs).


aut_rename(Num, Aut0, Aut) :-
	Aut0 = aut(Qs0, QFs0, Q00, Delta0),
	maplist(state_rename(Num), Qs0, Qs),
	maplist(state_rename(Num), QFs0, QFs),
	delta_rename(Delta0, Num, Delta),
	state_rename(Num, Q00, Q0),
	Aut = aut(Qs, QFs, Q0, Delta).

append_sort(A, B, C) :-
	append(A, B, C1),
	sort(C1, C).

aut_union(A1, A2, U) :-
	aut_rename(1, A1, A11),
	aut_rename(2, A2, A22),
	A11 = aut(Q1,QF1,S1,Delta1),
	A22 = aut(Q2,QF2,S2,Delta2),
	append(Q1, Q2, Allstates),
	length(Allstates, N),
	Q0 = start(N),
	Trans1 = [delta(Q0,epsilon,S1),delta(Q0,epsilon,S2)|Delta2],
	append_sort(Delta1, Trans1, Trans2),
	append_sort([start(N)|Q1], Q2, Qs),
	append_sort(QF1, QF2, QF),
	U = aut(Qs,QF,Q0,Trans2).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
			NDFA -> DFA conversion
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

is_deterministic(Aut) :-
	Aut = aut(_Qs,_QF,_Q0,Delta),
	\+ memberchk(delta(_,epsilon,_), Delta),
	\+ exists_choicepath(Delta).

exists_choicepath(Delta) :-
	member(delta(Q0,Seq,Q1), Delta),
	member(delta(Q0,Seq,Q2), Delta),
	Q1 \== Q2.

transform_table([], _Q, [], []).
transform_table([A|As], Q, [R|Rs], [T|Ts]) :-
	T = t(Q,A,R),
	transform_table(As, Q, Rs, Ts).


state_table(Alphabet, Delta, Q, Table) :-
	maplist(ndfa_allreachable(Q,Delta), Alphabet, Rs),
	transform_table(Alphabet, Q, Rs, Table).

ndfa_table(NDFA, Table) :-
	NDFA = aut(Qs,_QF,_Q0,Delta),
	delta_alphabet(Delta, Alphabet),
	maplist(state_table(Alphabet,Delta), Qs, Tables),
	flatten(Tables, Table).



symbol_union(Table, A, Qs, Us) :-
	phrase(symbol_union(Qs,Table,A), Us1),
	sort(Us1, Us).

symbol_union([], _Table, _A) --> [].
symbol_union([Q|Qs], Table, A) -->
	{ memberchk(t(Q,A,Rs), Table) },
	dlist(Rs),
	symbol_union(Qs, Table, A).

ttable1([t(q0, 0, [q0, q1, q2]), t(q0, 1, [q1, q2]), t(q0, 2, [q2]), t(q1, 0, []),t(q1, 1, [q1, q2]), t(q1, 2, [q2]), t(q2, 0, []), t(q2, 1, []), t(q2, 2, [q2])]).

dlist([]) --> [].
dlist([E|Es]) --> [E], dlist(Es).

state_union([], _Table, _Qs, []).
state_union([A|As], Table, Qs, [NewTrans|Transs]) :-
	symbol_union(Table, A, Qs, Us),
	NewTrans = u(Qs, A, Us),
	state_union(As, Table, Qs, Transs).


:- ttable1(T), symbol_union(T,1,[q1,q2],[q1,q2]).

% final state in DFA: if one of its "sub"-states is final in NDFA
is_dfafinal(Fs, d(Qs)) :-
	sort(Qs, Qs1),
	intersection(Fs, Qs1, Is), % Fs must be sorted
	Is \= [].


dfa_transform_table(u(Q,A,R), delta(d(Q),A,d(R))).

can_reach_final_from_start(Q0, QF, Delta) :-
	ndfa_reach_through_epsilon(Q0, Delta, R),
	memberchk(R, QF).

ndfa_dfa(NDFA, NDFA) :- is_deterministic(NDFA).
ndfa_dfa(NDFA, DFA) :-
	\+ is_deterministic(NDFA),
	NDFA = aut(_,QF,Q0,Delta),
	ndfa_table(NDFA, Table),
	delta_alphabet(Delta, Alphabet),
	state_union(Alphabet, Table, [Q0], Us0),
	saturate_det_table(Alphabet, Table, Us0, Us),
	maplist(dfa_transform_table, Us, DFADelta),
	delta_states(DFADelta, DFAStates1),
	DFAStates2 = [d([Q0])|DFAStates1],
	sort(DFAStates2, DFAStates),
	sort(QF, QF1),
	sublist(is_dfafinal(QF1), DFAStates, DFAFinal),
	( can_reach_final_from_start(Q0, QF, Delta) ->
		DFAFinal1 = [d([Q0])|DFAFinal]
	;
		DFAFinal = DFAFinal1
	),
	sort(DFAFinal1, DFAFinal2),
	DFA = aut(DFAStates,DFAFinal2,d([Q0]),DFADelta).


dettable_notcovered(Us, Q) :-
	member(u(_,_,Q), Us),
	\+ memberchk(u(Q,_,_), Us).

saturate_det_table(Alphabet, Table, Us0, Us) :-
	( dettable_notcovered(Us0, Qs) ->
		state_union(Alphabet, Table, Qs, Us1),
		append(Us0, Us1, Us2),
		saturate_det_table(Alphabet, Table, Us2, Us)
	;
		Us0 = Us
	).



% Rs: all states that are reachable from Q via symbol A
ndfa_allreachable(Q, Delta, A, Rs) :-
	findall(R, ndfa_reachable(Q,A,Delta,R), Rs0),
	sort(Rs0, Rs).


tndfa(aut([q0,q1,q2],[q2],q0,[delta(q0,0,q0),delta(q0,epsilon,q1),delta(q1,1,q1),delta(q1,epsilon,q2),delta(q2,2,q2)])).


ndfa_reach_through_epsilon(Q, Delta, R) :-
	ndfa_reach_through_epsilon(Q, Delta, [], R).

ndfa_reach_through_epsilon(Q, Delta, Visited0, R) :-
	member(delta(Q,epsilon,E), Delta),
	( R = E
	; \+ memberchk(E, Visited0),
	  ndfa_reach_through_epsilon(E, Delta, [E|Visited0], R)
	).

ndfa_reachable(Q, A, Delta, R) :-
	ndfa_reachable(Q, A, Delta, [], R).

ndfa_reachable(Q, A, Delta, _Visited, R) :-
	member(delta(Q,A,S), Delta),
	( R = S
	; ndfa_reach_through_epsilon(S, Delta, R)
	).
ndfa_reachable(Q, A, Delta, Visited0, R) :-
	ndfa_reach_through_epsilon(Q, Delta, E),
	\+ memberchk(E, Visited0),
	ndfa_reachable(E, A, Delta, [E|Visited0], R).


delta_alphabet(Delta, As) :-
	findall(A, (member(delta(_,A,_),Delta), A \= epsilon), As1),
	sort(As1, As).

delta_states(Delta, Qs) :-
	findall(Q, member(delta(Q,_,_),Delta), Qs1),
	findall(R, member(delta(_,_,R),Delta), Qs2),
	append_sort(Qs1, Qs2, Qs).


:- delta_alphabet([delta(a,1,b),delta(c,2,d),delta(f,2,g),delta(2,epsilon,5)],[1,2]).

:- delta_states([delta(a,1,b),delta(c,2,d),delta(f,2,g),delta(2,epsilon,5)],[2,5,a,b,c,d,f,g]).

:- N = aut([q0, q1, q2], [q2], q0, [delta(q0, 0, q0), delta(q0, epsilon, q1), delta(q1, 1, q1), delta(q1, epsilon, q2), delta(q2, 2, q2)]), ndfa_dfa(N, aut([d([]), d([q0]), d([q0, q1, q2]), d([q1, q2]), d([q2])],          [d([q0]), d([q0, q1, q2]), d([q1, q2]), d([q2])],         d([q0]),           [delta(d([q0]), 0, d([q0, q1, q2])), delta(d([q0]), 1, d([q1, q2])), delta(d([q0]), 2, d([q2])), delta(d([q0, q1, q2]), 0, d([q0, q1, q2])), delta(d([q0, q1, q2]), 1, d([q1, q2])), delta(d([q0, q1, q2]), 2, d([q2])), delta(d([q1, q2]), 0, d([])), delta(d([q1, q2]), 1, d([q1,q2])), delta(d([q1, q2]), 2, d([q2])), delta(d([q2]), 0, d([])), delta(d([q2]),1, d([])), delta(d([q2]), 2, d([q2])), delta(d([]), 0, d([])), delta(d([]), 1, d([])), delta(d([]), 2, d([]))])).


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		       Complement of automaton
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


complete_states([], _Alphabet, _Trap, _Delta) --> [].
complete_states([Q|Qs], Alphabet, Trap, Delta) -->
	complete_state(Alphabet, Q, Trap, Delta),
	complete_states(Qs, Alphabet, Trap, Delta).

complete_state([], _Q, _Trap, _Delta) --> [].
complete_state([A|As], Q, Trap, Delta) -->
	( { member(delta(Q,A,_), Delta) } ->
		[]
	;
		[delta(Q,A,Trap)]
	),
	complete_state(As, Q, Trap, Delta).

:- phrase(complete_state([1,2],a,trap,[delta(a,1,b)]), Cs), Cs = [delta(a, 2, trap)].

trapdelta(Trap, A, delta(Trap,A,Trap)).

% complete an automaton with respect to a given alphabet
aut_complete(Alphabet, Aut, Comp) :-
	Aut = aut(Qs,QF,Q0,Delta),
	length(Qs, LQs),
	Trap = trap(LQs),
	phrase(complete_states(Qs,Alphabet,Trap,Delta), Delta1),
	append_sort(Delta, Delta1, Delta2),
	maplist(trapdelta(Trap), Alphabet, TrapDeltas),
	append_sort(Delta2, TrapDeltas, Delta3),
	delta_states(Delta3, CompStates),
	Comp = aut(CompStates, QF, Q0, Delta3).


not_in(List, Elem) :-
	\+ memberchk(Elem, List).

exists_path(Q, Q, _Delta, _Visited0) --> [].
exists_path(Q, R, Delta, Visited0) -->
	{ member(delta(Q,Seq,S), Delta) },
	{ \+ memberchk(S, Visited0) },
	[Seq],
	exists_path(S, R, Delta, [S|Visited0]).

exists_path(Q, R, Delta, Path) :-
	phrase(exists_path(Q,R,Delta,[]), Path).


can_reach_final_through_0(Q, QF, Delta) :-
	can_reach_final_through_0(Q, QF, Delta, []).

can_reach_final_through_0(Q, QF, _Delta, _Visited) :-
	memberchk(Q, QF).
can_reach_final_through_0(Q, QF, Delta, Visited) :-
	member(delta(Q,Seq,R), Delta),
	\+ memberchk(R, Visited),
	maplist(=(0), Seq),
	can_reach_final_through_0(R, QF, Delta, [R|Visited]).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   There is a slight subtlety involved when building the complement,
   resulting from the ambiguity of number representation: an arbitrary
   amount of zeros can be appended - therefore, we can not turn a
   state into a final state if an *actual* final state can be reached
   via zeros. I call this a "pseudo"-final state.  Also, the initial
   state can not turn into a halting state unless there exists a
   proper path to it (that makes sense in terms of the accepted
   number), or the problem is already reduced to reachability.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

pseudo_final_state(QF, Delta, Q) :-
	once(can_reach_final_through_0(Q, QF, Delta)).


aut_complement(Aut, Complement) :-
	ndfa_dfa(Aut, DFA),
	DFA = aut(DQs,DQF,DQ0,Delta),
	( Delta == [] ->
		( memberchk(DQ0, DQF) ->
			Complement = aut(DQs,[],DQ0,[])
		;
			Complement = aut(DQs,[DQ0|DQF],DQ0,[])
		)
	;
		Delta = [delta(_,Seq,_)|_],
		length(Seq, LSeq),
		length(Binary, LSeq),
		Binary in 0..1,
		findall(Binary, label(Binary), Binaries),
		aut_complete(Binaries, DFA, Complete),
		Complete = aut(Qs,QF,Q0,Complete_Delta),
		sublist(not_in(QF), Qs, CFinals1),
		sublist(pseudo_final_state(QF,Complete_Delta), CFinals1, CFinals2),
		list_delete(CFinals2, CFinals1, CFinals3),
		( exists_path(Q0, Q0, Delta, [_|_]) ->
			CFinals3 = CFinals4
		;
			delete(CFinals3, Q0, CFinals4)
		),
		Complement = aut(Qs,CFinals4,Q0,Complete_Delta)
	).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
		     Delete a track of automaton
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

delete_nth([L|Ls0], Curr0, N, Ls) :-
	( Curr0 == N ->
		Ls = Ls0
	;
		Curr is Curr0 + 1,
		Ls = [L|Ls1],
		delete_nth(Ls0, Curr, N, Ls1)
	).

delete_nth(Ls0, N, Ls) :- delete_nth(Ls0, 0, N, Ls).

:- delete_nth([a,b,c,d], 2, [a,b,d]).

delta_remove_track(Track, delta(Q0,Seq0,Q1), delta(Q0,Seq,Q1)) :-
	( Track == 0, Seq0 = [_] ->
		Seq = epsilon
	;
		delete_nth(Seq0, Track, Seq)
	).

aut_delete_track(Aut0, Track, Aut) :-
	Aut0 = aut(Qs,QF,Q0,Delta0),
	maplist(delta_remove_track(Track), Delta0, Delta),
	Aut = aut(Qs,QF,Q0,Delta).



/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
			    Term rewriting
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

coeffs_negate([], []).
coeffs_negate([C0*V|Ss], [C*V|Ns]) :-
	C is (-1)*C0,
	coeffs_negate(Ss, Ns).

rebalance(Ls0-Rs0, Ls-Rs) :-
	Ls0 = [CL|Ls0Rest],
	Rs0 = [CR|Rs0Rest],
	C is CR - CL,
	Rs = [C],
	coeffs_negate(Rs0Rest, Neg),
	append(Ls0Rest, Neg, Ls1),
	sumup(Ls1, Ls).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  normalize input formula to:
    A <= c
    A = c
    A /\ B (and)
    A \/ B (or)
    not(F)
    exists(X, F)
  furthermore, simplify atomic formulas to
      a_1 * x_1 +  a_2 * x_2 + ... + a_n * x_n  (<)=   c
  and represent this using lists.
  A = B could be reduced to A <= B /\ B <= A, but we don't.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


nf(Term0, NF) :-
	Term0 =.. [Op,A0,B0],
	memberchk(Op, [\/,/\]),
	nf(A0, NA),
	nf(B0, NB),
	NF =.. [Op,NA,NB].
nf(Term0, NF) :-
	Term0 =.. [Op,A0,B0],
	memberchk(Op, [<=,=]),
	fix_rewrite(A0, A),
	fix_rewrite(B0, B),
	term_summands(A, As),
	term_summands(B, Bs),
	rebalance(As-Bs, NAs-NBs),
	NF =.. [Op,NAs,NBs].
nf(A0 >= B0, NF) :- nf(B0 <= A0, NF).
nf(A0 < B0, NF) :- nf(A0 + 1 <= B0, NF).
nf(A0 > B0, NF) :- nf(B0 < A0, NF).
nf(forall(X,F), NF) :- nf(not(exists(X,not(F))), NF).
nf(not(F), not(NF)) :- nf(F, NF).
nf(exists(X,F), exists(X,NF)) :- nf(F, NF).
nf(A => B, NF) :- nf(not(A) \/ B, NF).


xyz(x).
xyz(y).
xyz(z).

mvar(V) :-
	atom(V),
	atom_chars(V, Cs),
	Cs = [First|_],
	xyz(First).

coeff_one([], []).
coeff_one([S|Ss], [C|Cs]) :-
	( mvar(S) ->
		C = 1*S
	;
		C = S
	),
	coeff_one(Ss, Cs).

% compute summand list of term
term_summands(S) -->
	{ integer(S) },
	[S].
term_summands(V) -->
	{ mvar(V) },
	[V].
term_summands(A*B) --> [A*B].
term_summands(A+B) -->
	term_summands(A),
	term_summands(B).

term_summands(T, Ss) :-
	phrase(term_summands(T), Ss1),
	coeff_one(Ss1, Ss).


:- term_summands(1+2+3*x+y, [1,2,3*x,1*y]).


% like term_summands, for the other way around (list bound, construct term)
% computes a left-associatve term and can therefore not replace the former
summands(S) -->
	[S],
	{ integer(S) }.
summands(V) -->
	[V],
	{ mvar(V) }.
summands(A*B) -->
	[A*B].
summands(A*B+N) -->
	[A*B],
	summands(N).
summands(S+N) -->
	[S],
	{ integer(S) },
	summands(N).
summands(S+N) -->
	[S],
	{ mvar(S) },
	summands(N).

summands(S, Ss) :- phrase(summands(S), Ss).

:- summands(1+1, [1,1]).

varadd(Vars0, Coeff, Var, Vars) :-
	( member(v(C,Var), Vars0) ->
		C1 is C + Coeff,
		delete(Vars0, v(C,Var), Vars1),
		Vars = [v(C1,Var)|Vars1]
	;
		Vars = [v(Coeff,Var)|Vars0]
	).


vars_sumlist([]) --> [].
vars_sumlist([E|Es]) -->
	vars_sumlist_(E),
	vars_sumlist(Es).

vars_sumlist_(c(C)) --> [C].
vars_sumlist_(v(Coeff, V)) --> [Coeff*V].
vars_sumlist_(unmod(Us)) --> dlist(Us).


sumup(Ss0, Ss) :-
	sumup(Ss0, [unmod([]),c(0)], Vars),
	phrase(vars_sumlist(Vars), Ss).


sumup([], Vars, Vars).
sumup([S|Ss], Vars0, Vars) :-
	( mvar(S) ->
		sumup([1*S|Ss], Vars0, Vars)
	; S = A*B, integer(A), mvar(B) ->
		varadd(Vars0, A, B, Vars1),
		sumup(Ss, Vars1, Vars)
	; integer(S) ->
		memberchk(c(C), Vars0),
		delete(Vars0, c(C), Vars1),
		C1 is C + S,
		Vars2 = [c(C1)|Vars1],
		sumup(Ss, Vars2, Vars)
	; % could not simplify
		memberchk(unmod(U), Vars0),
		delete(Vars0, unmod(U), Vars1),
		Vars2 = [unmod([S|U])|Vars1],
		sumup(Ss, Vars2, Vars)
	).

:- sumup([5,2*x], [v(3,x),c(1)], [v(5,x),c(6)]).

:- sumup([1,2,5*x,1*x], [6*x,3]).



% compute the fix-point
fix_rewrite(A0, A) :-
	once(rewrite(A0, A1)),
	term_summands(A1, Ss1),
	sumup(Ss1, Ss2),
	sort(Ss1, Ss1_NF),
	sort(Ss2, Ss2_NF),
	summands(A2, Ss2_NF),
	(A0 == A1, Ss1_NF == Ss2_NF ->
		A2 = A
	;
		fix_rewrite(A2, A)
	).

rewrite(A, A) :- (integer(A) ; mvar(A) ).
rewrite(-V, (-1)*V) :- mvar(V).
rewrite(A0+B0, A+B) :-
	rewrite(A0, A),
	rewrite(B0, B).
rewrite(A0-B0, A+(-1)*B) :-
	rewrite(A0, A),
	rewrite(B0, B).
rewrite(A0*B0, C) :-
	rewrite(A0, A),
	rewrite(B0, B),
	( integer(A), integer(B) ->
		C is A*B
	; A = (A1 + A2) ->
		C = A1 * B + A2*B
	; B = (B1 + B2) ->
		C = A*B1 + A*B2
	; B = (B1*B2) ->
		C = (A*B1)*B2
	; integer(B) ->
		C = B*A
	;
		C = A*B
	).

:- fix_rewrite(1+3+4*x+5+y+x+3*(7+y),Rs), Rs = 30 + (4*y+5*x).

formula_variables(Term2) -->
	{ Term2 =.. [_Op,X,Y] },
	formula_variables(X),
	formula_variables(Y).
formula_variables(not(Term)) -->
	formula_variables(Term).
formula_variables(S) -->
	{ integer(S) },
	[].
formula_variables(-V) -->
	{ mvar(V) },
	[V].
formula_variables(V) -->
	{ mvar(V) },
	[V].

formula_variables(T, Vs) :-
	phrase(formula_variables(T), Vs0),
	sort(Vs0, Vs). % remove duplicates


nf_quantified_variables(exists(X, F)) -->
	[X],
	nf_quantified_variables(F).
nf_quantified_variables(not(Term)) -->
	nf_quantified_variables(Term).
nf_quantified_variables(Term2) -->
	{ Term2 =.. [Op,X,Y] },
	{ memberchk(Op, [/\,\/]) },
	nf_quantified_variables(X),
	nf_quantified_variables(Y).
nf_quantified_variables(Term2) -->
	{ Term2 =.. [Op|_] },
	{ memberchk(Op, [=,<=]) },
	[].

nf_quantified_variables(T, Vs) :-
	phrase(nf_quantified_variables(T), Vs0),
	sort(Vs0, Vs). % remove duplicates


sub_variables([]) --> [].
sub_variables([T]) -->
	{ integer(T) },
	[].
sub_variables([_Coeff*Var|Ss]) -->
	[Var],
	sub_variables(Ss).

nf_variables(exists(_X,F)) -->
	nf_variables(F).
nf_variables(not(F)) -->
	nf_variables(F).
nf_variables(Term2) -->
	{ Term2 =.. [Op,X,Y] },
	{ memberchk(Op, [/\,\/]) },
	nf_variables(X),
	nf_variables(Y).
nf_variables(Term2) -->
	{ Term2 =.. [Op,X,_Y] },
	{ memberchk(Op, [=,<=]) },
	sub_variables(X).

nf_variables(N, Vs) :-
	phrase(nf_variables(N), Vs).
				% order of variables must match the
				% tracks in the automaton


merge_sub([], _Left, []).
merge_sub([V0|Vs0], Left, [V|Vs]) :-
	( memberchk(Coeff*V0, Left) ->
		V = Coeff*V0
	;
		V = 0*V0
	),
	merge_sub(Vs0, Left, Vs).


merge_variables(Variables, not(Term0), not(Term)) :-
	merge_variables(Variables, Term0, Term).
merge_variables(Variables, exists(Var,Term0), exists(Var,Term)) :-
	merge_variables([Var|Variables], Term0, Term).
merge_variables(Variables, Term0, Term) :-
	Term0 =.. [Op,A0,B0],
	memberchk(Op, [\/,/\]),
	merge_variables(Variables, A0, A),
	merge_variables(Variables, B0, B),
	Term =.. [Op,A,B].
merge_variables(Variables, Term0, Term) :-
	Term0 =.. [Op,A0,B0],  % B0: constant
	memberchk(Op, [<=,=]),
	merge_sub(Variables, A0, A),
	Term =.. [Op,A,B0].

list_delete([], Ls, Ls).
list_delete([D|Ds], Ls0, Ls) :-
	delete(Ls0, D, Ls1),
	list_delete(Ds, Ls1, Ls).

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
   Main entry point for normalization.  Takes a formula, makes sure
   each variable occurs in each (in)equality, (if necessary, enforce
   with coefficient 0), same order everywhere.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


formula_normalized(F, NF) :-
	formula_variables(F, Variables),
	nf(F, NF0),
	nf_quantified_variables(NF0, QVs),
	list_delete(QVs, Variables, Variables1),
	merge_variables(Variables1, NF0, NF).


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% "x = 1":
taut1(aut([q(0), q(1)], [q(0)], q(1), [delta(q(0), [0], q(0)), delta(q(1), [1], q(0))])).

% "x <= 4":
taut2(aut([q(-1), q(0), q(1), q(2), q(4)], [q(0), q(1), q(2), q(4)], q(4), [delta(q(-1), [0], q(-1)), delta(q(-1), [1], q(-1)), delta(q(0), [0], q(0)), delta(q(0), [1], q(-1)), delta(q(1), [0], q(0)), delta(q(1), [1], q(0)), delta(q(2), [0], q(1)), delta(q(2), [1], q(0)), delta(q(4), [0], q(2)), delta(q(4), [1], q(1))])).

% "x >= 5"  (-x <= -5)
taut3(aut([q(-5), q(-3), q(-2), q(-1), q(0)], [q(0)], q(-5), [delta(q(-5), [0], q(-3)), delta(q(-5), [1], q(-2)), delta(q(-3), [0], q(-2)), delta(q(-3), [1], q(-1)), delta(q(-2), [0], q(-1)), delta(q(-2), [1], q(-1)), delta(q(-1), [0], q(-1)), delta(q(-1), [1], q(0)), delta(q(0), [0], q(0)), delta(q(0), [1], q(0))])).


:- taut1(A1), taut2(A2), aut_intersection(A1,A2,Int), \+ automaton_empty(Int).
:- taut1(A1), taut3(A3), aut_intersection(A1,A3,Int), automaton_empty(Int).

eq_sat(Eq) :-
	eq_automaton(Eq,A),
	\+ automaton_empty(A).

:- teq(Eq), eq_sat(Eq).

ineq_sat(Ineq) :-
	ineq_automaton(Ineq,A),
	\+ automaton_empty(A).



summand_coeff(Coeff*_Var, Coeff).

coeff_mterm(C, m(C,_)).


% syntax check

is_term(Var) :-	mvar(Var).
is_term(-Var) :- mvar(Var).
is_term(Num) :- integer(Num).
is_term(Term) :-
	Term =.. [Op,A,B],
	memberchk(Op, [+,-,*]),
	is_term(A),
	is_term(B).

is_atomic_formula(AF) :-
	AF =.. [Op,A,B],
	memberchk(Op, [<,>,=,<=,>=]),
	is_term(A),
	is_term(B).

is_formula(not(F)) :- is_formula(F).
is_formula(AF) :- is_atomic_formula(AF).
is_formula(F) :-
	F =.. [Op,A,B],
	memberchk(Op, [/\,\/,=>]),
	is_formula(A),
	is_formula(B).
is_formula(F) :-
	F =.. [Op,A,B],
	memberchk(Op, [forall,exists]),
	mvar(A),
	is_formula(B).


formula_automaton(X = [Y], A) :-
	maplist(summand_coeff, X, Cs),
	maplist(coeff_mterm, Cs, Ms),
	Eq = eq(Ms,Y),
	eq_automaton(Eq, A).
formula_automaton(X <= [Y], A) :-
	maplist(summand_coeff, X, Cs),
	maplist(coeff_mterm, Cs, Ms),
	Ineq = ineq(Ms,Y),
	ineq_automaton(Ineq, A).
formula_automaton(not(F), A) :-
	formula_automaton(F, A1),
	aut_complement(A1, A).
formula_automaton(X /\ Y, A) :-
	formula_automaton(X, A1),
	formula_automaton(Y, A2),
	aut_intersection(A1, A2, A).
formula_automaton(X \/ Y, A) :-
	formula_automaton(X, A1),
	formula_automaton(Y, A2),
	aut_union(A1, A2, A).
formula_automaton(exists(Var, F), A) :-
	formula_automaton(F, A1),
	nf_variables(F, Vs),
	nf_quantified_variables(F, QVs),
	list_delete(QVs, Vs, Vs1),
	nth0(N, Vs1, Var),
	aut_delete_track(A1, N, A).

syntax_ok(Formula) :-
	( ground(Formula), is_formula(Formula) ->
		true
	;
		format("\n\tInvalid syntax.\n"),
		fail
	).


satisfiable(Formula) :-
	syntax_ok(Formula),
	satisfiable_(Formula).

satisfiable_(Formula) :-
	formula_normalized(Formula, NF),
	formula_automaton(NF, A),
	\+ automaton_empty(A).

valid(Formula) :-
	syntax_ok(Formula),
	\+ satisfiable_(not(Formula)).


do_tests :-
	valid(y > 1 /\ x = 3 /\ x + y < 19  =>  x + 19 > y), % p. 102
	valid(y > 1 /\ x = 3 /\ not(x + y < 19) =>  y + y > y), % p. 102
	valid(x = 3 /\ y = 1 => 3*x + y = 10), % p. 103
	valid(2*y + 3*x = 30 /\ not(not(x=0)) => y = 15), % p. 78 (addendum)
	valid(y = 15 /\ x = 0 => y = 15),
	\+ valid(x+y > 0),
	valid(forall(x,exists(y,x+y > 5))),
	valid(exists(y, x+y > 5)),
	valid(forall(x,exists(y, y > 2))),
	valid(exists(y, x+y > 1)),
	valid(x>0 \/ x = 0),
	\+ valid(x+y>1).


%:- do_tests.