:- module(_, [poly_exp/3], [assertions, regtypes, nativeprops]).

% Ralph Haygood (based on Prolog version by Rick McGeer
% based on Lisp version by R. P. Gabriel)
%
% raise a polynomial (1+x+y+z) to the 10th power (symbolically)
%
% The regtype of polynomials
% E.g.
%   x ---> poly(x,[term(1,1)])
%   y ---> poly(y,[term(1,1)])
%   x+y ----> poly(x,[term(0,poly(y,[term(1,1)])),term(1,1)])
%   1+x+y ----> poly(x,[term(0,poly(y,[term(0,1),term(1,1)])),term(1,1)])

% ----------------------------------------------------------------------

:- regtype poly/1.
poly(A) :-
        num(A).
poly(poly(A,B)) :-
        term(A),
        polyterms(B).

:- regtype polyterms/1.
polyterms([]).
polyterms([X|Xs]) :- polyterm(X), polyterms(Xs).

:- regtype polyterm/1.
polyterm(term(X,Poly)) :-
	term(X),
	poly(Poly).

:- pred poly_add(A,B,C) : poly * poly * term => poly * poly * poly.

poly_add(poly(Var,Terms1), poly(Var,Terms2), poly(Var,Terms)) :- !,
        term_add(Terms1, Terms2, Terms).
poly_add(poly(Var1,Terms1), poly(Var2,Terms2), poly(Var1,Terms)) :-
        Var1 @< Var2, !,
        add_to_order_zero_term(Terms1, poly(Var2,Terms2), Terms).
poly_add(Poly, poly(Var,Terms2), poly(Var,Terms)) :- !,
        add_to_order_zero_term(Terms2, Poly, Terms).
poly_add(poly(Var,Terms1), C, poly(Var,Terms)) :- !,
        add_to_order_zero_term(Terms1, C, Terms).
poly_add(C1, C2, C) :-
        C is C1 + C2.

:- pred term_add(A,B,C) : polyterms * polyterms * term => polyterms * polyterms * polyterms.

term_add([], X, X) :- !.
term_add(X, [], X) :- !.
term_add([term(E,C1)|Terms1], [term(E,C2)|Terms2], [term(E,C)|Terms]) :- !,
        poly_add(C1, C2, C),
        term_add(Terms1, Terms2, Terms).
term_add([term(E1,C1)|Terms1], [term(E2,C2)|Terms2], [term(E1,C1)|Terms]) :-
        E1 < E2, !,
        term_add(Terms1, [term(E2,C2)|Terms2], Terms).
term_add(Terms1, [term(E2,C2)|Terms2], [term(E2,C2)|Terms]) :-
        term_add(Terms1, Terms2, Terms).

add_to_order_zero_term([term(0,C1)|Terms], C2, [term(0,C)|Terms]) :- !,
        poly_add(C1, C2, C).
add_to_order_zero_term(Terms, C, [term(0,C)|Terms]).

:- pred poly_exp(N, Poly, R) : num * poly * term => num * poly * poly.

poly_exp(0, _, 1) :- !.
poly_exp(N, Poly, Result) :-
        N/\1 =:= 0, !,
        M is N>>1,
        poly_exp(M, Poly, Part),
        poly_mul(Part, Part, Result).
poly_exp(N, Poly, Result) :-
        M is N - 1,
        poly_exp(M, Poly, Part),
        poly_mul(Poly, Part, Result).

:- pred poly_mul(A,B,C) : poly * poly * term => poly * poly * poly.

poly_mul(poly(Var,Terms1), poly(Var,Terms2), poly(Var,Terms)) :- !,
        term_mul(Terms1, Terms2, Terms).
poly_mul(poly(Var1,Terms1), poly(Var2,Terms2), poly(Var1,Terms)) :-
        Var1 @< Var2, !,
        mul_through(Terms1, poly(Var2,Terms2), Terms).
poly_mul(P, poly(Var,Terms2), poly(Var,Terms)) :- !,
        mul_through(Terms2, P, Terms).
poly_mul(poly(Var,Terms1), C, poly(Var,Terms)) :- !,
        mul_through(Terms1, C, Terms).
poly_mul(C1, C2, C) :-
        C is C1 * C2.

:- pred term_mul(A,B,C) : polyterms * polyterms * term => polyterms * polyterms * polyterms.

term_mul([], _, []) :- !.
term_mul(_, [], []) :- !.
term_mul([Term|Terms1], Terms2, Terms) :-
        single_term_mul(Terms2, Term, PartA),
        term_mul(Terms1, Terms2, PartB),
        term_add(PartA, PartB, Terms).

:- pred single_term_mul(A,B,C) : polyterms * polyterm * term => polyterms * polyterm * polyterms.

single_term_mul([], _, []).
single_term_mul([term(E1,C1)|Terms1], term(E2,C2), [term(E,C)|Terms]) :-
        E is E1 + E2,
        poly_mul(C1, C2, C),
        single_term_mul(Terms1, term(E2,C2), Terms).

:- pred mul_through(A,B,C) : polyterms * poly * term => polyterms * poly * polyterms.

mul_through([], _, []).
mul_through([term(E,Term)|Terms], Poly, [term(E,NewTerm)|NewTerms]) :-
        poly_mul(Term, Poly, NewTerm),
        mul_through(Terms, Poly, NewTerms).
