polynomialmani.pl/polymani.pl

%% -*- Mode: Prolog-*-
%% vim: set softtabstop=4 shiftwidth=4 tabstop=4 expandtab:

/**
 *
 * polimani.pl
 *
 * Assignment 1 - Polynomial Manipulator
 * Programming in Logic - DCC-FCUP
 *
 * Diogo Peralta Cordeiro
 * up201705417@fc.up.pt
 *
 * Hugo David Cordeiro Sales
 * up201704178@fc.up.pt
 *
 *********************************************
 *   Follows 'Coding guidelines for Prolog'  *
 * https://doi.org/10.1017/S1471068411000391 *
 *********************************************/

/*
 * Import the Constraint Logic Programming over Finite Domains library
 * Essentially, this library improves the way Prolog deals with integers,
 * allowing more predicates to be reversible.
 * For instance, number(N) is always false, which prevents the
 * reversing of a predicate.
 */
:- use_module(library(clpfd)).
/*
 * Import Constraint Logic Programming for Reals library, which is somewhat
 * similar to clpfd, but for real numbers
 */
:- use_module(library(clpr)).

/*
 * The porter_stem library implements the stemming algorithm described by
 * Porter in Porter, 1980, ``An algorithm for suffix stripping''.
 * The library also includes a functional tokenizer
*/
:- use_module(library(porter_stem)).

                 /*******************************
                 *              NLP             *
                 *******************************/

%% polyplay() is det
%
%  Interactive prompt for the NLP Interface
%
polyplay :-
    %% Replace the prompt and save the old one
    prompt(OldPrompt, '> '),
    %% Read a line as codes, from the 'user_input' stream
    read_line_to_codes(user_input, InCodes),
    %% Restore old prompt
    prompt(_, OldPrompt),
    %% Use the porter_stem library to tokenize the input
    %% This does more than splitting at the spaces, such as
    %% splitting at operators. Aditionally, gives
    %% atoms already and handles 'P1' well
    tokenize_atom(InCodes, LA),
    (
        %% If we read a 'bye', terminate
        LA == [bye],
        write("See ya"),
        nl,
        !
    ;
        (
            %% Parse the input into a tree
            parse_input(TIn, LA, NC),
            (
                %% If the tree is empty, it means nothing was understood
                TIn == void,
                writeln("I didn't understand what you want."),
                writeln(NC)
            ;
                (
                    %% If there is unconsumed input, it means there was a syntatic error
                    NC \== [],
                    write("Syntax error in token: "),
                    %% Take the head of the list
                    cons(H, _, NC),
                    write(H),
                    nl
                ;
                    %% Otherwise, process the parse tree and execute the commands within
                    process_input(TIn)
                )
            )
        ;
            %% Parsing failed
            writeln("Could not parse input, so I didn't understand what you want.")
        ),
	%% Go back to the beginning
	polyplay
    ),
    !.

text2poly(S, P) :-
    split_string(S, " ", "\r\t", LS),
    parse_polynomial(T, LS, []),
    polynomial_tree_to_polynomial(T, P1),
    simpoly(P1, P).

%% process_input(+Tree) is det
%
%  Execute the commands from the parse tree
%
process_input(command(CL, void)) :-
    %% Process only command left
    do_process_input(CL).
process_input(command(CL, TCR)) :-
    %% Process first command
    do_process_input(CL),
    %% If there's a tree on the right
    TCR \== void,
    %% recurse
    process_input(TCR).

%% do_process_input(+Tree) is det
%
%  Process a single command from the input
%  Takes, as the only argument, the tree representing
%  the work to be done
%
do_process_input(help_menu) :-
    writeln("Please allow me to introduce myself"),
    writeln("I'm a man[ual] of wealth and taste"),
    writeln("I've been around for a long, long semester"),
    writeln("Saved many a man's soul and faith"),
    writeln("Pleased to meet you"),
    writeln("Hope you guess my name"),
    writeln("But what's puzzling you"),
    writeln("Is the nature of my game"),
    nl,
    writeln("I'm a Polynomial Manipulator and the following commands are available:"),
    writeln("=> Show - Allows to print a polynomial mathematically"),
    writeln("=> Multiply - Allows to make multiplications"),
    writeln("=> Simplify - Allows to simplify a given polynomial"),
    writeln("=> Add - Allows to make sums"),
    writeln("=> bye - Hello darkness, my old friend"),
    writeln("Use 'tell me about {command}' to learn more about a specific command"),
    nl,
    writeln("Furthermore, I'm capable of memorizing polynomials during runtime. To learn more on that, type: tell me about storage").
do_process_input(help(show)) :-
    writeln("It's almost an echo of what you said, but a mathy one."),
    writeln("Example query:"),
    writeln("> show two plus x squared"),
    writeln("2+x^2").
do_process_input(help(multiply)) :-
    writeln("Multiplies a polynomial represented as an expression by a scalar resulting in a second polynomial. The two first arguments are assumed to be ground. The polynomial resulting from the sum is in simplified form."),
    writeln("Example query:"),
    writeln("> multiply three by two plus x squared"),
    writeln("3*x^2+6").
do_process_input(help(simplify)) :-
    writeln("Simplifies a polynomial represented as an expression as another polynomial as an expression."),
    writeln("Example query:"),
    writeln("> simplify two plus two plus one times y"),
    writeln("y+4").
do_process_input(help(add)) :-
    writeln("Adds two polynomials as expressions resulting in a third one. The two first arguments are assumed to be ground. The polynomial resulting from the sum is in simplified form."),
    writeln("Some example queries:"),
    writeln("> add two times x to four times x"),
    writeln("6*x").
do_process_input(help(storage)) :-
    writeln("Polynomials, pressed between the entries of my storage"),
    writeln("Polynomials, simplified through the ages just like predicates"),
    writeln("Quiet goal come floating down"),
    writeln("And settle softly to the ground"),
    writeln("Like golden atom leaves around my root"),
    writeln("I touched them and they burst apart with sweet polynomials"),
    writeln("Sweet polynomials"),
    nl,
    writeln("Storage manipulation is better illustrated with examples. Some example queries:"),
    writeln("Asking me to memorize something:"),
    writeln("> show two plus x squared as P1"),
    writeln("P1 = 2+x^2"),
    writeln("> multiply three by P1 as P2"),
    writeln("P2 = 3*x^2+6"),
    nl,
    writeln("Asking me to forget something:"),
    writeln("> forget P1 and show stored polynomials"),
    writeln("P2 = 6+3*x^2").
do_process_input(help(bye)) :-
    writeln("There must be some kind of way outta here"),
    writeln("Said the joker to the thief"),
    writeln("There's too much confusion"),
    writeln("I can't get no relief").
do_process_input(show_stored_polynomials) :-
    %% If the command is 'show_stored_polynomials'
    %% Store in D the list of all results of
    %% polynomial_store, which is a dynamic predicate
    findall(nm(X,Y), polynomial_store(X,Y), D),
    %% Print the list of stored variables
    print_all_stored_variables(D).
do_process_input(show(store(P), T)) :-
    %% If the command was to store and show
    %% Delegate the work of storing
    do_process_input(store(P, T)),
    %% Show the variable and the polynomial. Delegate
    do_process_input(show(P, T)).
do_process_input(show(load(P), void)) :-
    %% Make sure there's a variable
    P \== void,
    (
        %% Check if the variable is stored
        polynomial_store(P, T),
        %% Show the variable and the polynomial. Delegate
        do_process_input(show(P, T))
    ;
        writeln("Variable not stored")
    ).
do_process_input(show(void, T)) :-
    %% Make sure there's a tree
    T \== void,
    %% Convert it to a flat polynomial
    polynomial_tree_to_polynomial(T, Pl),
    writeln(Pl).
do_process_input(show(P, T)) :-
    %% Make sure the arguments are in the right format and not one
    %% that should be handled by other clauses
    P \== void,
    P \== store(_),
    P \== load(_),
    T \== void,
    %% Write out the variable and the normal
    %% representation of the polynomial
    write(P),
    write(" = "),
    %% Convert the input tree to a flat polynomial, to print
    polynomial_tree_to_polynomial(T, Pl),
    writeln(Pl).
do_process_input(store_simplified(V,PTNS)) :-
    polynomial_tree_to_polynomial(PTNS, PNS),
    simpoly(PNS,P),
    assertz(polynomial_store(V, P)),
    write(V),
    write(" = "),
    write(P),
    nl.
do_process_input(store(P, T)) :-
    (
        %% If there is a polynomial stored with the same name
        polynomial_store(P, _),
        %% Delete all occurences of it
        retractall(polynomial_store(P, _))
    ;
        %% Otherwise, move forward and don't do anything special
        true
    ),
    !,
    %% Dynamically add to the end of the predicate, effectively storing it
    assertz(polynomial_store(P, T)).
do_process_input(forget(P)) :-
    %% Remove all stored predicates with this variable
    retractall(polynomial_store(P, _)).
do_process_input(simplify(PT)) :-
    %% Convert the polynomial to it's flat representation
    polynomial_tree_to_polynomial(PT, P),
    %% Use simpoly from the UI layer to simplify, so it
    %% gives nice error messages
    simpoly(P, SP),
    writeln(SP).
do_process_input(store_multiplication(TN, PT, V)) :-
    polynomial_tree_to_polynomial(TN, N),
    polynomial_tree_to_polynomial(PT, P),
    %% Left is always number. Enforced by parser
    scalepoly(P, N, P2),
    simpoly(P2, SP),
    assertz(polynomial_store(V, SP)),
    write(V),
    write(" = "),
    write(SP),
    nl.
do_process_input(multiply(TN, PT)) :-
    %% Flatten both
    polynomial_tree_to_polynomial(TN, N),
    polynomial_tree_to_polynomial(PT, P),
    %% Use function from UI layer
    %% Left is always the number. Enforced by parser
    scalepoly(P, N, P2),
    simpoly(P2, SP),
    writeln(SP).
do_process_input(op(+, TN, PT)) :-
    polynomial_tree_to_polynomial(TN, N),
    polynomial_tree_to_polynomial(PT, P),
    addpoly(N, P, P2),
    simpoly(P2, SP),
    writeln(SP).
do_process_input(store_addition(TN, PT, V)) :-
    polynomial_tree_to_polynomial(TN, N),
    polynomial_tree_to_polynomial(PT, P),
    addpoly(N, P, P2),
    simpoly(P2, SP),
    assertz(polynomial_store(V, SP)),
    write(V),
    write(" = "),
    write(SP),
    nl.

%% print_all_stored_variables
%
%  Prints the stored variables, given the list
%
print_all_stored_variables([nm(X,T) | TS]) :-
    %% Print
    write(X),
    write(" = "),
    polynomial_tree_to_polynomial(T, P),
    writeln(P),
    %% Recurse on the right
    print_all_stored_variables(TS).
print_all_stored_variables([]).

%% polynomial_tree_to_polynomial(+tree, -polynomial) is det
%
%  Flatten a polynomail tree into a simple polynomial
%
polynomial_tree_to_polynomial(op(neg, T), P1) :-
    %% If it matched on the version with a node, don't go back
    !,
    %% Delegate
    polynomial_tree_to_polynomial(T, P),
    atom_concat('-', P, P1).
polynomial_tree_to_polynomial(op(Op, TL, TR), P) :-
    %% If it matched on the version with a node, don't go back
    !,
    %% Recurse on the left and right
    polynomial_tree_to_polynomial(TL, PL),
    polynomial_tree_to_polynomial(TR, PR),
    %% Convert the results of the recursion to atoms
    term_to_atom(PL, TermL),
    term_to_atom(PR, TermR),
    %% Concat the parts
    atom_concat(TermL, Op, Temp),
    atom_concat(Temp, TermR, PA),
    %% Convert back to terms
    term_to_atom(P, PA).
polynomial_tree_to_polynomial(load(P), Pl) :-
    %% If there's a load tag, don't go back
    !,
    %% Retrieve the polynomial from the store
    polynomial_store(P, TP),
    %% Recurse as if there was a normal tree here
    polynomial_tree_to_polynomial(TP, Pl).
%% If we get here, it means we didn't match any of the previous
%% clauses, so the tree is already flat
polynomial_tree_to_polynomial(T, T).

%% special_word_number(?W:Atom, ?D:Int) is det
%
%  Definition of a Alphabetical and Numerical relation
%  The third argument refers to different types of precedence:
%   - f means a single unit digit, which doesn't
%       combine on the right, but can on the left
%   - g means a number which doesn't combine left or right
%   - fy means a number which can optionally combine on the right
%   - xfy means a number which requires a number on the left and
%         an optional one on the right
%
special_word_number(zero,      0,       f).
special_word_number(a,         1,       f).
special_word_number(one,       1,       f).
special_word_number(two,       2,       f).
special_word_number(three,     3,       f).
special_word_number(four,      4,       f).
special_word_number(five,      5,       f).
special_word_number(six,       6,       f).
special_word_number(seven,     7,       f).
special_word_number(eight,     8,       f).
special_word_number(nine,      9,       f).
special_word_number(ten,       10,      g).
special_word_number(eleven,    11,      g).
special_word_number(twelve,    12,      g).
special_word_number(thirteen,  13,      g).
special_word_number(fourteen,  14,      g).
special_word_number(fifteen,   15,      g).
special_word_number(sixteen,   16,      g).
special_word_number(seventeen, 17,      g).
special_word_number(eighteen,  18,      g).
special_word_number(nineteen,  19,      g).
special_word_number(twenty,    20,      fy).
special_word_number(thirty,    30,      fy).
special_word_number(forty,     40,      fy).
special_word_number(fifty,     50,      fy).
special_word_number(sixty,     60,      fy).
special_word_number(seventy,   70,      fy).
special_word_number(eighty,    80,      fy).
special_word_number(ninety,    90,      fy).
special_word_number(hundred,   100,     xfy).
special_word_number(thousand,  1000,    xfy).
special_word_number(million,   1000000, xfy).
%% special_word_number(half,      0.5,     xf).
%% special_word_number(third,     0.33333, xf).
%% special_word_number(quarter,   0.25,    xf).

%% parse_number_explicit(+precedence, +tree_left, -tree,
%%                       +stream, -not_consumed) is det
%
%  Use the precedence of the last number and the tree
%  on the left to create the tree on the right, parsing
%  the input stream.
%
%  You probably want to call this with void in the first
%  two arguments. They are for internal use
%
parse_number_explicit(void, void, T, [WN | In], NC) :-
    %% Entry point
    %% Convert a word to a number with precedence f, g, or fy
    special_word_number(WN, N, P),
    member(P, [f, g, fy]),
    !,
    %% Recurse with the previous precedence
    %% and tree as arguments. The same in the other clauses
    parse_number_explicit(P, N, T, In, NC).
parse_number_explicit(fy, NL, T, [WN | In], NC) :-
    special_word_number(WN, N, f),
    !,
    %% Add them on the left tree and recurse
    parse_number_explicit(f, op(+, NL, N), T, In, NC).
parse_number_explicit(xfy, TL, T, [WN | In], NC) :-
    TL \= void,
    special_word_number(WN, N, P),
    member(P, [f, g, fy]),
    !,
    %% Add the parts
    parse_number_explicit(P, op(+, TL, N), T, In, NC).
parse_number_explicit(_, TL, T, [WN | In], NC) :-
    %% Whatever the precedence on the left
    special_word_number(WN, N, xfy),
    %% Enforce there was a left argument
    TL \= void,
    !,
    %% Multiply this and the left
    parse_number_explicit(xfy, op(*, TL, N), T, In, NC).
parse_number_explicit(P, TL, T, [and, WN | In], NC) :-
    %% And is part of a valid number if there's a valid
    %% number on the left and right
    special_word_number(WN, _, _),
    %% Passthrough
    parse_number_explicit(P, TL, T, [WN | In], NC),
    !.
parse_number_explicit(_, T, T, [WN | In], [WN | In]) :-
    %% If there's a tree
    T \= void,
    %% And can't consume more words as
    %% numbers (fail on first), then we're done
    not(special_word_number(WN, _, _)),
    !.
parse_number_explicit(_, T, T, [], []) :-
    %% If consumed all innput and there's a
    %% tree, we're done
    T \= void,
    !.

%% parse_floating_number(-tree, +stream, -not_consumed) is det
%
%  Parse a floating point number
%
parse_floating_number(N) -->
    [N],
    { number(N) }.
parse_floating_number(op('.', TL, TR)) -->
    %% A float is a node with a dot as the operator
    %% If there's a number on the left
    parse_integer_number(TL),
    %% Followed by either point or dot
    [Sep],
    {
        member(Sep, [point, dot]),
        !
    },
    %% Followed by another number
    parse_positive_integer_number(TR).
parse_floating_number(TN) -->
    %% Or any integer
    parse_integer_number(TN).
%% Tests:
%% ?- parse_floating_number(TN, [two], _).
%@ TN = 2.
%% ?- parse_floating_number(TN, [two, dot, 4], _).
%@ TN = op('.', 2, 4).
%% ?- parse_floating_number(TN, [negative, two, dot, 4], _).
%@ TN = op('.', op(neg, 2), 4).
%@ TN = op('.', op(neg, 2), 4).

%% parse_positive_integer_number(-tree, +stream, -not_consumed) is det
%
%  Parse an positive integer number
%
parse_positive_integer_number(N) -->
    [N],
    {
        integer(N),
        %% CLPFD, a number greater than 0
        N #>= 1, !
    }.
parse_positive_integer_number(T) -->
    parse_number_explicit(void, void, T).

%% parse_integer_number(-tree, +stream, -not_consumed) is det
%
%  Parse an integer number
%
parse_integer_number(N) -->
    [N],
    { integer(N), ! }.
parse_integer_number(op(neg, T)) -->    % TODO
    %% A number can start with "negative", to negate it
    [negative],
    parse_integer_number(T).
parse_integer_number(T) -->
    parse_number_explicit(void, void, T).
%% Tests:
%% ?- parse_integer_number(T, [two], _).
%@ T = 2.
%% ?- parse_integer_number(T, [43], _).
%@ T = 43.
%% ?- parse_integer_number(T, [nineteen, two], _).
%@ false.
%% ?- parse_integer_number(T, [twenty], _).
%@ T = 20.
%% ?- parse_integer_number(T, [twenty, point, two], NC).
%@ T = op('.', 20, 2),
%@ NC = [].
%% ?- parse_integer_number(T, [twenty, twenty], _).
%@ false.
%% ?- parse_integer_number(T, [twenty, one], _).
%@ T = op(+, 20, 1).
%% ?- parse_integer_number(T, [negative, one, hundred], _).
%@ T = op(neg, op(*, 1, 100)) ;
%@ false.
%% ?- parse_integer_number(T, [three, hundred], _).
%@ T = op(*, 3, 100).
%% ?- parse_integer_number(T, [twenty, hundred], _).
%@ T = op(*, 20, 100).
%% ?- parse_integer_number(T, [twenty, one, hundred], _).
%@ T = op(*, op(+, 20, 1), 100).
%% ?- parse_integer_number(T, [two, hundred, and, one], _).
%@ T = op(+, op(*, 2, 100), 1).
%% ?- parse_integer_number(T, [twenty, one, hundred, and, twenty, one], _).
%@ T = op(+, op(+, op(*, op(+, 20, 1), 100), 20), 1).
%% ?- parse_integer_number(T, [twenty, one, hundred, and, twenty, one, foo, bar, blah], NC).
%@ T = op(+, op(+, op(*, op(+, 20, 1), 100), 20), 1),
%@ NC = [foo, bar, blah].
%% ?- parse_integer_number(T, [twenty, one, hundred, and, bleg, twenty, quux, one, foo, bar], NC).
%@ T = op(*, op(+, 20, 1), 100),
%@ NC = [and, bleg, twenty, quux, one, foo, bar].
%% ?- parse_integer_number(T, [two, hundred, thousand], _).
%@ T = op(*, op(*, 2, 100), 1000).
%% ?- parse_integer_number(T, [twenty, one, hundred, thousand], _).
%@ T = op(*, op(*, op(+, 20, 1), 100), 1000).
%% ?- parse_integer_number(T, [thirty, five, million, five, hundred, thirty, four], _).
%@ T = op(+, op(+, op(*, op(+, op(*, op(+, 30, 5), 1000000), 5), 100), 30), 4).
%% ?- parse_integer_number(T, [foo, five, million], NC).
%@ false.

%% nlp_parse_power(?List, ?List) is det
%
%  Parse powers
%
parse_power(op(^, TB, 2)) -->
    parse_polynomial_variable(TB),
    [squared].
parse_power(op(^, TB, 3)) -->
    parse_polynomial_variable(TB),
    [cubed].
parse_power(op(^, TB, TN)) -->
    parse_polynomial_variable(TB),
    [raised, to],
    parse_positive_integer_number(TN).
parse_power(TB) -->
    parse_polynomial_variable(TB).


%% parse_operation(+Op, +stream, -not_consumed) is det
%
%  Associate an operator with a word
%
parse_operation(-) --> [minus] | [-].
parse_operation(+) --> [plus]  | [+].
parse_operation(*) --> [times] | [*].

%% parse_polynomial_operand(-tree) is det
%
%  Parse a valid operand for either side of an
%  operation in a polynomial
%
parse_polynomial_operand(T) -->
    parse_floating_number(T).
parse_polynomial_operand(T) -->
    parse_power(T).
parse_polynomial_operand(load(T)) -->
    %% Tag the variable, to be loaded later
    parse_stored_variable(T),
    { ! }.
%% Tests:
%% ?- parse_polynomial_operand(T, [two], NC).
%@ T = 2,
%@ NC = [] ;
%@ false.

%% Declare polynomial_store as a dynamic predicate with two arguments
%% Used to store and retrieve polynomials associated with a variable
%% First being the variable name and the second its content
:- dynamic polynomial_store/2.

%% parse_stored_variable(-var) is det
%
%  Parse a valid variable name, for storage
%  A valid varaible is one that would be valid in Prolog
%
parse_stored_variable(P) -->
    [P],
    { valid_store_variable(P) }.
%% Tests:
%% ?- parse_stored_variable(P, ['P1'], _).
%@ P = 'P1'.
%% ?- parse_stored_variable(P, ['x1'], _).
%@ false.

%% parse_stored_variable(+var) is det
%
%  Check if var is a valid variable name, for storage
%  A valid varaible is one that would be valid in Prolog
%
valid_store_variable(P) :-
    %% Get the list of codes of the atom
    atom_codes(P, L),
    %% Extract the first element
    cons(F, R, L),
    %% Partially builtin facilities:
    %% Check if it's a valid start to a Prolog variable
    code_type(F, prolog_var_start),
    %% Check that each of the following is a valid
    %% continuation to a variable name
    maplist(code_type_swap(prolog_identifier_continue), R).

%% code_type_swap - swap arguments of code_type, so
%% it can be used with maplist. No lambdas...
code_type_swap(X, Y) :- code_type(Y, X).

%% parse_polynomial_variable(-base, +stream, -not_consumed) is det
%
%  Parse a polynomial variable
%
parse_polynomial_variable(B) -->
    [B],
    { polynomial_variable(B) }.

%% parse_polynomial(-tree, +stream, -not_consumed) is det
%
%  Parse a polynomial. Delegates to explicit variant
%
parse_polynomial(T) -->
    %% Ignore "polynomial", if followed by a valid polynomial
    [polynomial],
    { ! },
    %% Delegate
    parse_polynomial(T).
parse_polynomial(T) -->
    %% Delegate
    parse_polynomial_explicit(_-_, T),
    !.
%% Tests:
%% ?- parse_polynomial(T, [], _).
%@ false.
%% ?- parse_polynomial(T, [two], _).
%@ T = 2,
%@ NC = [].
%% ?- parse_polynomial(T, [two, times, three], _).
%@ T = op(*, 2, 3),
%@ NC = [].
%% ?- parse_polynomial(T, [two, times, three, plus, four], _).
%@ T = op(+, op(*, 2, 3), 4).
%% ?- parse_polynomial(T, [two, plus, three, times, four], _).
%@ T = op(+, 2, op(*, 3, 4)).
%% ?- parse_polynomial(T, [two, plus, three, times, four, plus, six, times, five], _).
%@ T = op(+, 2, op(+, op(*, 3, 4), op(*, 6, 5))).
%% ?- parse_polynomial(T, [two, times, times, two], NC).
%@ T = void,
%@ NC = [two, times, times, two].
%% ?- parse_polynomial(T, [two, plus, x, times, four], _).
%@ T = op(+, 2, op(*, x, 4)).
%% ?- parse_polynomial(T, [two, plus, x, times, four, plus, y, raised, to, five], _).
%@ T = op(+, 2, op(+, op(*, x, 4), op(^, y, 5))).
%% ?- parse_polynomial(T, [two, plus, two, plus, one, times, y], _).
%@ T = op(+, op(+, 2, 2), op(*, 1, y)).
%% ?- parse_polynomial(T, [polynomial, 2, plus, 3, plus, 4, y], _).
%@ T = op(+, op(+, 2, 3), op(*, 4, y)).
%% ?- parse_polynomial(T, [2, x], _).
%@ true.

%% parse_polynomial_explicit(+tree-at, -tree) is det
%
%  Explicitly parse a polynomial. You probably don't want
%  to call this directly. The first argument is a difference
%  structure of the tree on the left; the second of the pair
%  represents where to put the rest of the tree, in
%  subsequent passes
%
%  Call with _-_ on the first argument, outputs the
%  tree on the right
%
parse_polynomial_explicit(void-_, T) -->
    %% Entry point. If the three on the left is void
    %% Parse an operand and an op
    parse_polynomial_operand(TL),
    parse_operation(Op),
    %% If consumed an op, double down and don't go back
    !,
    %% Recurse with the tree on the left populated with
    %% what we found and a reference to the place to
    %% place the next tree
    parse_polynomial_explicit(op(Op, TL, TRP)-TRP, T).
parse_polynomial_explicit(void-_, T) -->
    %% Entry point. If the three on the left is void
    parse_polynomial_operand(TL),
    %% There was no op, assume a times
    %% Recurse with the tree on the left populated with
    %% what we found and a reference to the place to
    %% place the next tree
    parse_polynomial_explicit(op(*, TL, TRP)-TRP, T).
parse_polynomial_explicit(TLP-TL, T) -->
    %% Parse an operand on the left and place it in the tree
    %% on the left, through the difference structure,
    parse_polynomial_operand(TL),
    %% parse either a minus or a plus
    { member(COp, [-, +]) },
    parse_operation(COp),
    !,
    %% Recurse on the right with add; position the sub tree on the right
    parse_polynomial_explicit(op(COp, TLP, TRP)-TRP, T).
parse_polynomial_explicit(TLP-T, TLP) -->
    %% Parse an operand on the left, placing the result of the
    %% recusion on the TLP, through the difference structure and return that
    parse_polynomial_operand(TL),
    %% If we find a times,
    parse_operation(*),
    !,
    %% Place the operand on the left, the following on the
    %% right and get a T to place in the tree on the left from above
    parse_polynomial_explicit(op(*, TL, TRP)-TRP, T).
parse_polynomial_explicit(TLP-T, TLP) -->
    %% Same as above, but without an operator,
    parse_polynomial_operand(TL),
    %% Assume times
    parse_polynomial_explicit(op(*, TL, TRP)-TRP, T),
    !.
parse_polynomial_explicit(TLP-TL, TLP) -->
    { TLP \= void },
    %% Base case. Parse a last operand
    parse_polynomial_operand(TL),
    !,
    { TL \= void }.
parse_polynomial_explicit(void-_, T) -->
    %% Base case. Parse a single operand
    parse_polynomial_operand(T),
    !,
    { T \= void }.

%% parse_command(-tree, +stream, -not_consumed) is semidet
%
%  Parse each individual command
%
parse_command(help_menu) -->
    [help].
parse_command(help(C)) -->
    [tell, me, about, C].
parse_command(show_stored_polynomials) -->
    [show, stored, polynomials].
parse_command(forget(P)) -->
    [forget],
    parse_stored_variable(P).
parse_command(show(store(P), T)) -->
    [show],
    parse_polynomial(T),
    [as],
    parse_stored_variable(P).
parse_command(show(P, void)) -->
    [show],
    parse_stored_variable(P).
parse_command(show(void, T)) -->
    [show],
    parse_polynomial(T),
    { nonvar(T) }.
parse_command(store(P, T)) -->
    [let],
    parse_stored_variable(P),
    [be],
    parse_polynomial(T).
parse_command(store(P, T)) -->
    [store],
    parse_polynomial(T),
    [as],
    parse_stored_variable(P).
parse_command(store_simplified(V, P)) -->
    [simplify],
    parse_polynomial(P),
    [as],
    [V].
parse_command(simplify(T)) -->
    [simplify],
    parse_polynomial(T).
parse_command(store_multiplication(TN, TP, V)) -->
    [multiply],
    parse_floating_number(TN),
    [by],
    parse_polynomial(TP),
    [as],
    [V].
parse_command(store_multiplication(TN, TP, V)) -->
    [multiply],
    parse_polynomial(TP),
    [by],
    parse_floating_number(TN),
    [as],
    [V].
parse_command(multiply(TN, TP)) -->
    [multiply],
    parse_floating_number(TN),
    [by],
    parse_polynomial(TP).
parse_command(multiply(TN, TP)) -->
    [multiply],
    parse_polynomial(TP),
    [by],
    parse_floating_number(TN).
parse_command(op(-, TN, TP)) -->
    [sub],
    parse_polynomial(TN),
    [X],
    { member(X, [to, from, with]) },
    parse_polynomial(TP).
parse_command(store_addition(TN, TP, V)) -->
    [add],
    parse_polynomial(TN),
    [X],
    { member(X, [to, with]) },
    parse_polynomial(TP),
    [as],
    [V].
parse_command(op(+, TN, TP)) -->
    [add],
    parse_polynomial(TN),
    [X],
    { member(X, [to, with]) },
    parse_polynomial(TP).
%% Tests:
%% ?- parse_command(T, [show, 3], NC).
%@ T = show(void, 3),
%@ NC = [] .
%% ?- parse_command(T, [add, 3, plus, x, to, 4, plus, x], NC).
%@ false.
%@ T = op(+, op(+, 3, x), op(+, 4, x)),
%@ NC = [] ;
%@ false.

%% parse_input(-tree) is det
%
%  Parse each command and string into a list
%
parse_input(command(TCL, TCR)) -->
    %% Result is a struct with a command on the left
    %% and a struct of commands on the right
    %% parse a single command
    parse_command(TCL),
    %% separator
    [and],
    !,
    %% Recurse
    parse_input(TCR).
parse_input(command(TC, void)) -->
    %% If there's no separator, we're done recursing
    %% parse a single command
    parse_command(TC).
%% ?- parse_input(CT, [show, 3], _).
%@ CT = command(show(void, 3), void).


                 /*******************************
                 *             UI               *
                 *******************************/

/*
    poly2list/2 transforms a list representing a polynomial (second
    argument) into a polynomial represented as an expression (first
    argument) and vice-versa.
*/
poly2list(P, L) :-
    is_polynomial_valid_in_predicate(P, "poly2list"),
    polynomial_to_list(P, L),
    !.

/*
    simpolylist/2 simplifies a polynomial represented as a list into
    another polynomial as a list.
*/
simpoly_list(L, S) :-
    is_polynomial_as_list_valid_in_predicate(L, "simpoly_list"),
    simplify_polynomial_as_list(L, S),
    !.

/*
    simpoly/2 simplifies a polynomial represented as an expression
    as another polynomial as an expression.
*/
simpoly(P, S) :-
    is_polynomial_valid_in_predicate(P, "simpoly"),
    simplify_polynomial(P, S),
    !.
%% Tests:
%% ?- simpoly(2+2+1*y, S).
%@ S = y+4.

/*
    scalepoly/3 multiplies a polynomial represented as an expression by a scalar
    resulting in a second polynomial. The two first arguments are assumed to
    be ground. The polynomial resulting from the sum is in simplified form.
*/
scalepoly(P1, C, S) :-
    is_polynomial_valid_in_predicate(P1, "scalepoly"),
    is_number_valid_in_predicate(C, "scalepoly"),
    scale_polynomial(P1, C, S),
    !.
%% Tests:
%% ?- scalepoly(3*x*z+2*z, 4, S).
%@ S = 12*x*z+8*z.
%% ?- scalepoly(3*x*z+2*z, 2, S).
%@ S = 6*x*z+4*z.

/*
    addpoly/3 adds two polynomials as expressions resulting in a
    third one. The two first arguments are assumed to be ground.
    The polynomial resulting from the sum is in simplified form.
*/
addpoly(P1, P2, S) :-
    is_polynomial_valid_in_predicate(P1, "addpoly"),
    is_polynomial_valid_in_predicate(P2, "addpoly"),
    add_polynomial(P1, P2, S),
    !.
%% Tests:
%% ?- addpoly(3 + x, 3 - x, S).
%@ S = 6.

%% is_polynomial_valid_in_predicate(+T, +F) is det
%
%  Returns true if valid polynomial, fails with UI message otherwise.
%  The failure message reports which polynomial is invalid and in which
%  predicate the problem ocurred.
%
is_polynomial_valid_in_predicate(P, _) :-
    %% If P is a valid polynomial, return true
    polynomial(P),
    !.
is_polynomial_valid_in_predicate(P, F) :-
    %% Otherwise, write the polynomial and fails
    write("Invalid polynomial in "),
    write(F),
    write(": "),
    writeln(P),
    fail.
%% Tests:
%% ?- is_polynomial_valid_in_predicate(1-x, "Test").
%@ true.
%% ?- is_polynomial_valid_in_predicate(a*4-0*x, "Test").
%@ Invalid polynomial in Test: a*4-0*x
%@ false.

%% is_polynomial_as_list_valid_in_predicate(+L, +F) is det
%
%  Returns true if the polynomial represented as list is valid,
%  fails with UI message otherwise.
%  The failure message reports which polynomial is invalid and
%  in which predicate the problem ocurred.
%
is_polynomial_as_list_valid_in_predicate(L, F) :-
    %% If L is a valid polynomial, return true
    list_to_polynomial(L, P),
    is_polynomial_valid_in_predicate(P, F).
%% Tests:
%% ?- is_polynomial_as_list_valid_in_predicate([1], "Test").
%@ true.
%% ?- is_polynomial_as_list_valid_in_predicate([0*x, a*4], "Test").
%@ Invalid polynomial in Test: a*4+0*x
%@ false.

%% is_number_valid_in_predicate(+C:number, +F:string) is det
%
%  Validates that C is a number or prints F and it then it
%
is_number_valid_in_predicate(C, _) :-
    number(C),
    !.
is_number_valid_in_predicate(C, F) :-
    %% Writes the argument and fails
    write("Invalid number in "),
    write(F),
    write(": "),
    writeln(C),
    fail.


                 /*******************************
                 *            BACKEND           *
                 *******************************/

%% polynomial_variable_list(-List) is det
%
%  List of possible polynomial variables
%
polynomial_variable_list([x, y, z]).

%% polynomial_variable(?X:atom) is semidet
%
%  Returns true if X is a polynomial variable, false otherwise.
%
polynomial_variable(X) :-
    polynomial_variable_list(V),
    member(X, V).
%% Tests:
%% ?- polynomial_variable(x).
%@ true .
%% ?- polynomial_variable(a).
%@ false.

%% power(+X:atom) is semidet
%
%  Returns true if X is a power term, false otherwise.
%
power(P^N) :-
    %% CLPFD comparison. Reversible
    N #>= 1,
    polynomial_variable(P).
power(X) :-
    polynomial_variable(X).
%% Tests:
%% ?- power(x).
%@ true .
%% ?- power(a).
%@ false.
%% ?- power(x^1).
%@ true .
%% ?- power(x^3).
%@ true .
%% ?- power(x^(-3)).
%@ false.
%% ?- power(-x).
%@ false.
%% ?- power(X).
%@ X = x^_462546,
%@ _462546 in 1..sup ;
%@ X = y^_462546,
%@ _462546 in 1..sup ;
%@ X = z^_462546,
%@ _462546 in 1..sup ;
%@ X = x ;
%@ X = y ;
%@ X = z.

%% term(+N:atom) is semidet
%
%  Returns true if N is a term, false otherwise.
%
term(N) :-
    % If N is not a free variable
    nonvar(N),
    % Assert it as a number
    number(N).
term(N) :-
    % If N is a free variable and not compound
    not(compound(N)),
    var(N),
    % Assert it must be between negative and positive infinity
    % This uses the CLPR library, which makes this reversible,
    % whereas `number(N)` is always false, since it only succeeds
    % if the argument is bound (to a integer or float)
    {N >= 0; N < 0}.
term(X) :-
    power(X).
term(-X) :-
    power(X).
term(L * In) :-
    term(L),
    term(In).
%% Tests:
%% ?- term(2*x^3).
%@ true .
%% ?- term(x^(-3)).
%@ false.
%% ?- term(a).
%@ false.
%% ?- term(-1*x).
%@ true .
%% ?- term(-x).
%@ true .
%% ?- term((-3)*x^2).
%@ true .
%% ?- term(3.2*x).
%@ true .
%% ?- term(-x*(-z)).
%@ true .
%% ?- term(X).
%@ {X>=0.0} ;
%@ {X<0.0} ;
%@ X = x^_111514,
%@ _111514 in 1..sup ;
%@ X = y^_111514,
%@ _111514 in 1..sup ;
%@ X = z^_111514,
%@ _111514 in 1..sup ;
%@ X = x ;
%@ X = y ;
%@ X = z ;
%@ X = -x^_111522,
%@ _111522 in 1..sup ;
%@ X = -y^_111522,
%@ _111522 in 1..sup ;
%@ X = -z^_111522,
%@ _111522 in 1..sup ;
%@ X = -x ;
%@ X = -y ;
%@ X = -z ;

%% polynomial(+M:atom) is semidet
%
%  Returns true if polynomial, false otherwise.
%
polynomial(M) :-
    %% A polynomial is either a term
    term(M).
polynomial(L + In) :-
    %% Or a sum of terms
    polynomial(L),
    term(In).
polynomial(L - In) :-
    %% Or a subtraction of terms
    polynomial(L),
    term(In).
%% Tests:
%% ?- polynomial(x).
%@ true .
%% ?- polynomial(x^3).
%@ true .
%% ?- polynomial(3*x^7).
%@ true .
%% ?- polynomial(2+3*x+4*x*y^3).
%@ true .
%% ?- polynomial(2 - 3*x + 4*x*y^3).
%@ true .
%% ?- polynomial(a).
%@ false.
%% ?- polynomial(x^(-3)).
%@ false.
%% ?- polynomial(-x + 3).
%@ true .
%% ?- polynomial(-x - -z).
%@ true .

%% power_to_canon(+T:atom, -T^N:atom) is semidet
%
%  Returns a canon power term.
%
power_to_canon(T^N, T^N) :-
    polynomial_variable(T),
    % CLP(FD) operator to ensure N is different from 1,
    % in a reversible way
    N #\= 1.
power_to_canon(T, T^1) :-
    polynomial_variable(T).
%% Tests:
%% ?- power_to_canon(x, X).
%@ X = x^1 .
%% ?- power_to_canon(-x, X).
%@ false.
%@ X = -1*x^1 .
%% ?- power_to_canon(X, x^1).
%@ X = x .
%% ?- power_to_canon(X, x^4).
%@ X = x^4 .
%% ?- power_to_canon(X, a^1).
%@ false.
%% ?- power_to_canon(X, x^(-3)).
%@ X = x^ -3 .
%% ?- power_to_canon(X, -1*x^1).
%@ X = -x .

%% term_to_list(?T, ?List) is semidet
%
%  Converts a term to a list of its monomials and vice versa.
%  Can verify if term and monomials list are compatible.
%
term_to_list(L * N, [N | TS]) :-
    number(N),
    term_to_list(L, TS).
term_to_list(L * P, [P2 | TS]) :-
    power(P),
    power_to_canon(P, P2),
    term_to_list(L, TS).
term_to_list(L * -P, [-P2 | TS]) :-
    power(P),
    power_to_canon(P, P2),
    term_to_list(L, TS).
term_to_list(N, [N]) :-
    number(N).
term_to_list(P, [P2]) :-
    power(P),
    power_to_canon(P, P2).
term_to_list(-P, [-P2]) :-
    power(P),
    power_to_canon(P, P2).
%% Tests:
%% ?- term_to_list(1, X).
%@ X = [1] .
%% ?- term_to_list(-1, X).
%@ X = [-1] .
%% ?- term_to_list(x, X).
%@ X = [x^1] .
%% ?- term_to_list(-x, X).
%@ X = [-x^1] .
%% ?- term_to_list(2 * 3, X).
%@ X = [3, 2] .
%% ?- term_to_list(1*2*y*z*23*x*y*x^3*x, X).
%@ X = [x^1, x^3, y^1, x^1, 23, z^1, y^1, 2, 1] .
%% ?- term_to_list(1*2*y*z*23*x*y*(-1), X).
%@ X = [-1, y^1, x^1, 23, z^1, y^1, 2, 1] .
%% ?- term_to_list(X, [-1]).
%@ X = -1 .
%% ?- term_to_list(X, [x^1, -1]).
%@ X = -1*x .
%% ?- term_to_list(X, [-x^1]).
%@ X = -x .
%% ?- term_to_list(X, [y^1, x^1]).
%@ X = x*y .
%% ?- term_to_list(X, [x^4]).
%@ X = x^4 .
%% ?- term_to_list(X, [y^6, z^2, x^4]).
%@ X = x^4*z^2*y^6 .
%% ?- term_to_list(X, [y^6, z^2, x^4, -2]).
%@ X = -2*x^4*z^2*y^6 .
%% ?- term_to_list(X, [x^1, 0]).
%@ X = 0*x .
%% ?- term_to_list(X, [y^1, -2]).
%@ X = -2*y .

%% simplify_term(+Term_In:term, ?Term_Out:term) is det
%
%  Simplifies a given term.
%  This function can also be be used to verify if
%  a term is simplified.
%
simplify_term(Term_In, Term_Out) :-
    term_to_list(Term_In, L),
    %% Sort the list of numbers and power to group them,
    %% simplifying the job of `join_similar_parts_of_term`
    sort(0, @=<, L, L2),
    (
        %% If there's a 0 in the list, then the whole term is 0
        member(0, L2),
        Term_Out = 0
    ;
        %% Otherwise
        (
            %% If there's only one element, then the term was already simplified
            %% This is done so that the `exclude` following doesn't remove all ones
            length(L2, 1),
            Term_Out = Term_In
        ;
            %% Remove all remaining ones
            exclude(==(1), L2, L3),
            join_similar_parts_of_term(L3, L4),
            %% Reverse the list, since the following call gives the result in the
            %% reverse order otherwise
            reverse(L4, L5),
            term_to_list(Term_Out, L5)
        )
    ),
    % First result is always the most simplified form.
    !.
%% Tests:
%% ?- simplify_term(1, X).
%@ X = 1.
%% ?- simplify_term(x, X).
%@ X = x.
%% ?- simplify_term(2*y*z*x^3*x, X).
%@ X = 2*x^4*y*z.
%% ?- simplify_term(1*y*z*x^3*x, X).
%@ X = x^4*y*z.
%% ?- simplify_term(0*y*z*x^3*x, X).
%@ X = 0.
%% ?- simplify_term(6*y*z*7*x*y*x^3*x, X).
%@ X = 42*x^5*y^2*z.
%% ?- simplify_term(-x, X).
%@ X = -x.
%% ?- simplify_term(-x*y*(-z)*3, X).
%@ X = 3* -x* -z*y.
%% ?- simplify_term(a, X).
%@ false.
%% ?- simplify_term(x^(-3), X).
%@ false.

%% join_similar_parts_of_term(+List, -List) is det
%
%  Combine powers of the same variable in the given list.
%  Requires that the list be sorted.
%
join_similar_parts_of_term([P1, P2 | L], L2) :-
    %% If both symbols are powers
    power(P1),
    power(P2),
    %% Decompose them into their parts
    B^N1 = P1,
    B^N2 = P2,
    %% Sum the exponent
    N is N1 + N2,
    join_similar_parts_of_term([B^N | L], L2),
    % First result is always the most simplified form.
    !.
join_similar_parts_of_term([N1, N2 | L], L2) :-
    %% If they are both numbers
    number(N1),
    number(N2),
    %% Multiply them
    N is N1 * N2,
    join_similar_parts_of_term([N | L], L2),
    % First result is always the most simplified form.
    !.
join_similar_parts_of_term([X | L], [X | L2]) :-
    %% Otherwise consume one element and recurse
    join_similar_parts_of_term(L, L2),
    % First result is always the most simplified form.
    !.
join_similar_parts_of_term([], []).
%% Tests:
%% ?- join_similar_parts_of_term([3], T).
%@ T = [3].
%% ?- join_similar_parts_of_term([x^2], T).
%@ T = [x^2].
%% ?- join_similar_parts_of_term([x^1, x^1, x^1, x^1], T).
%@ T = [x^4].
%% ?- join_similar_parts_of_term([2, 3, x^1, x^2], T).
%@ T = [6, x^3].
%% ?- join_similar_parts_of_term([2, 3, x^1, x^2, y^1, y^6], T).
%@ T = [6, x^3, y^7].
%% ?- join_similar_parts_of_term([2, 3, -x^1, -x^2], T).
%@ T = [6, -x^1, -x^2].

%% simplify_polynomial(+P:atom, -P2:atom) is det
%
%  Simplifies a polynomial.
%
simplify_polynomial(0, 0) :-
    % 0 is already fully simplified. This is an
    % exception to the following algorithm
    !.
simplify_polynomial(P, P2) :-
    polynomial_to_list(P, L),
    simplify_polynomial_as_list(L, L2),
    list_to_polynomial(L2, P2),
    %% The first result is the most simplified one
    !.
%% Tests:
%% ?- simplify_polynomial(1, X).
%@ X = 1.
%% ?- simplify_polynomial(0, X).
%@ X = 0.
%% ?- simplify_polynomial(x, X).
%@ X = x.
%% ?- simplify_polynomial(x*x, X).
%@ X = x^2.
%% ?- simplify_polynomial(2 + 2, X).
%@ X = 2*2.
%% ?- simplify_polynomial(x + x, X).
%@ X = 2*x.
%% ?- simplify_polynomial(0 + x*x, X).
%@ X = x^2.
%% ?- simplify_polynomial(x^2*x + 3*x^3, X).
%@ X = 4*x^3.
%% ?- simplify_polynomial(x^2*x + 3*x^3 + x^3 + x*x*x, X).
%@ X = 6*x^3.
%% ?- simplify_polynomial(x^2*x + 3*x^3 + x^3 + x*x*4 + z, X).
%@ X = 5*x^3+4*x^2+z.
%% ?- simplify_polynomial(x^2*x + 3*x^3 - x^3 - x*x*4 + z, X).
%@ X = 3*x^3-4*x^2+z.
%% ?- simplify_polynomial(x + 1 + x, X).
%@ X = 2*x+1.
%% ?- simplify_polynomial(x + 1 + x + 1 + x + 1 + x, X).
%@ X = 4*x+3.

%% simplify_polynomial_as_list(+L1:List,-L3:List) is det
%
%  Simplifies a polynomial represented as a list.
%
simplify_polynomial_as_list(L, L14) :-
    %% Convert each term to a list
    maplist(term_to_list, L, L2),
    %% Sort each sublist so that the next
    %% sort gives the correct results
    maplist(sort(0, @>=), L2, L3),
    %% Sort the outer list
    sort(0, @>=, L3, L4),
    %% For each of the parts of the terms, join them
    maplist(join_similar_parts_of_term, L4, L5),
    %% Sort each of the sublists
    %% Done so the next call simplifies has less work
    maplist(sort(0, @=<), L5, L6),
    join_similar_terms(L6, L7),
    %% Exclude any sublist that includes a 0 (such as the
    %% equivalent to the term 0*x)
    exclude(member(0), L7, L8),
    %% Reverse each sublist, because the next call
    %% reverses the result
    maplist(reverse, L8, L9),
    maplist(term_to_list, L10, L9),
    %% Delete any 0 from the list
    delete(L10, 0, L11),
    %% Sort list converting back gives the result in the correct order
    sort(0, @=<, L11, L12),
    (
        %% If the list is empty, the result is a list with 0
        L12 = [], L13 = [0]
    ;
        %% Otherwise, this is the result
        L13 = L12
    ),
    %% Further make sure all terms are simplified
    maplist(simplify_term, L13, L14).
%% Tests:
%% ?- simplify_polynomial_as_list([x, 1, x^2, x*y, 3*x^2, 4*x], L).
%@ L = [1, 4*x^2, 5*x, x*y] .
%% ?- simplify_polynomial_as_list([1, x^2, x*y, 3*x^2, -4, -1*x], L).
%@ L = [-3, -1*x, 4*x^2, x*y] .
%% ?- simplify_polynomial_as_list([1, 1*x], L).
%@ L = [1, x] .
%% ?- simplify_polynomial_as_list([0*x, 0], L).
%@ L = [0] .

%% join_similar_terms(+P:List, -P2:List) is det
%
%  Joins similar sublists representing terms by using
%  `add_terms` to check if they can be merged and perform
%  the addition. Requires the list of list be sorted with
%    `maplist(sort(0, @>=), L, L2),
%     sort(0, @>=, L2, L3)`
%  and that the sublists to be sorted with
%  `sort(0, @=<)` since that is inherited from `add_terms`.
%
join_similar_terms([TL, TR | L], L2) :-
    %% Check if terms can be added and add them
    add_terms(TL, TR, T2),
    %% Recurse, accumulation on the first element
    join_similar_terms([T2 | L], L2),
    %% Give only first result. Red cut
    !.
join_similar_terms([X | L], [X | L2]) :-
    %% If a pair of elements can't be added, skip one
    %% and recurse
    join_similar_terms(L, L2),
    %% Give only first result. Red cut
    !.
join_similar_terms([], []).
%% Tests:
%% ?- join_similar_terms([[2, x^3], [3, x^3], [x^3]], L).
%@ L = [[6, x^3]].

%% term_to_canon(+T:List, -T2:List) is det
%
%  Adds the coefficient of the term as the first element of the list
%
%% Special cases to make this predicate reversible
term_to_canon([1], [1]) :-
    !.
term_to_canon(L2, [1 | L]) :-
    nonvar(L),
    L2 = L,
    !.
term_to_canon([-1], [-1]) :-
    !.
term_to_canon([-P | L2], [-1, P | L]) :-
    nonvar(L),
    L2 = L,
    !.
term_to_canon([N2 | L], [N | L]) :-
    number(N),
    N2 = N,
    !.
%% Normal case
term_to_canon(L, [N | L2]) :-
    term_to_canon_with_coefficient(N, L, L2),
    !.
%% Tests:
%% ?- term_to_canon([2], T).
%@ T = [2].
%% ?- term_to_canon([-x], T).
%@ T = [-1, x].
%% ?- term_to_canon([-x^3], T).
%@ T = [-1, x^3].
%% ?- term_to_canon([x^1], T).
%@ T = [1, x^1].
%% ?- term_to_canon([x^3], T).
%@ T = [1, x^3].
%% ?- term_to_canon([x^3, z], T).
%@ T = [1, x^3, z].
%% ?- term_to_canon([2, x^3], T).
%@ T = [2, x^3].
%% ?- term_to_canon([2, -x^3], T).
%@ T = [-2, x^3].
%% ?- term_to_canon([2, -x^3, -z], T).
%@ T = [2, x^3, z].
%% ?- term_to_canon(L, [-1]).
%@ L = [-1].
%% ?- term_to_canon(L, [1]).
%@ L = [1].
%% ?- term_to_canon(L, [-2]).
%@ L = [-2].
%% ?- term_to_canon(L, [-2, x]).
%@ L = [-2, x].
%% ?- term_to_canon(L, [1, x]).
%@ L = [x].
%% ?- term_to_canon(L, [-1, x]).
%@ L = [-x].
%% ?- term_to_canon(L, [1, x, z, y]).
%@ L = [x, z, y].
%% ?- term_to_canon(L, [-1, x, z, y]).
%@ L = [-x, z, y].

%% term_to_canon_with_coefficient(-N:number, +L:List, -L2:List) is semidet
%
%  Calculates the coefficient of the term and removes negations of powers,
%  accumulating the results in N
%
term_to_canon_with_coefficient(N, [N2 | TS], TS2) :-
    number(N2),
    term_to_canon_with_coefficient(N3, TS, TS2),
    N is N2 * N3,
    !.
term_to_canon_with_coefficient(N, [P | TS], [P2 | TS2]) :-
    sign_of_power(P, N2 * P2),
    term_to_canon_with_coefficient(N3, TS, TS2),
    N is N2 * N3,
    !.
term_to_canon_with_coefficient(N, [], []) :-
    nonvar(N);
    N = 1.
%% Tests:
%% ?- term_to_canon_with_coefficient(N, [x], L).
%@ N = 1,
%@ L = [x].
%% ?- term_to_canon_with_coefficient(N, [x, x^2, 2], L).
%@ N = 2,
%@ L = [x^1, x^2].
%% ?- term_to_canon_with_coefficient(N, [x, x^2, 2, 4, z], L).
%@ N = 8,
%@ L = [x, x^2, z].
%% ?- term_to_canon_with_coefficient(N, [x, x^2, 2, 4, -z], L).
%@ N = -8,
%@ L = [x, x^2, z].
%% ?- term_to_canon_with_coefficient(N, [x, -x^2, 2, 4, -z], L).
%@ N = 8,
%@ L = [x, x^2, z].
%% ?- term_to_canon_with_coefficient(N, L, [x]).
%@ N = 1,
%@ L = [x].
%% ?- term_to_canon_with_coefficient(N, L, [1]).
%@ N = 1,
%@ L = [1].
%% ?- term_to_canon_with_coefficient(N, L, [2]).
%@ N = 1,
%@ L = [2].

%% sign_of_power(P:power, P:term) is det
%
%  If there isn't a leading minus, multiplies the power by 1,
%  otherwise by a -1. This way it prefers the positive version.
%  Not idempotent
%
sign_of_power(P, 1*P) :-
    %% If P can't unify with a minus followed by an unnamed variable
    P \= -_,
    !.
sign_of_power(-P, -1*P).
%% Tests:
%% ?- sign_of_power(x, X).
%@ X = 1*x.
%% ?- sign_of_power(-x, X).
%@ X = -1*x.
%% ?- sign_of_power(X, 1*x).
%@ X = x.
%% ?- sign_of_power(X, -1*x).
%@ X = -x.

%% add_terms(+L:List, +In:List, -Result:List) is det
%
%  Adds two terms represented as list by adding
%  the coeficients if the power is the same.
%  Returns false if they can't be added
%  Requires the list of terms to be simplified.
%
add_terms([NL | TL], [NR | TR], [N2 | TL2]) :-
    %% Convert each term to a canon form. This ensures they
    %% have a number in front, so it can be added
    term_to_canon([NL | TL], [NL2 | TL2]),
    term_to_canon([NR | TR], [NR2 | TR2]),
    %% If the rest of the term is the same
    TL2 == TR2,
    %% Add the coeficients
    N2 is NL2 + NR2.
%% Tests
%% ?- add_terms([1], [1], In).
%@ In = [2].
%% ?- add_terms([x], [x], In).
%@ In = [2, x].
%% ?- add_terms([2, x^3], [x^3], In).
%@ In = [3, x^3].
%% ?- add_terms([2, x^3], [3, x^3], In).
%@ In = [5, x^3].
%% ?- add_terms([2, x^3], [3, x^2], In).
%@ false.

%% polynomial_to_list(+P:polynomial, -L:List) is det
%
%  Converts a polynomial in a list.
%
polynomial_to_list(L - T, [T2 | LS]) :-
    term(T),
    negate_term(T, T2),
    polynomial_to_list(L, LS),
    !.
polynomial_to_list(L + T, [T | LS]) :-
    term(T),
    polynomial_to_list(L, LS),
    !.
polynomial_to_list(T, [T]) :-
    term(T),
    !.
%% Tests:
%% ?- polynomial_to_list(2, S).
%@ S = [2].
%% ?- polynomial_to_list(x^2, S).
%@ S = [x^2].
%% ?- polynomial_to_list(x^2 + x^2, S).
%@ S = [x^2, x^2].
%% ?- polynomial_to_list(2*x^2+5+y*2, S).
%@ S = [y*2, 5, 2*x^2].
%% ?- polynomial_to_list(2*x^2+5-y*2, S).
%@ S = [-2*y, 5, 2*x^2].
%% ?- polynomial_to_list(2*x^2-5-y*2, S).
%@ S = [-2*y, -5, 2*x^2].

%% list_to_polynomial(+L:List, -P:Polynomial) is det
%
%  Converts a list in a polynomial.
%  An empty list will return false.
%
list_to_polynomial([T1|T2], P) :-
    % Start recursive calls until we are in the
    % end of the list. We know that the `-` will
    % always be at the left of a term.
    list_to_polynomial(T2, L1),
    (
        % If this is a negative term
        term_string(T1, S1),
        string_chars(S1, [First|_]),
        First = -,
        % Concat them
        term_string(L1, S2),
        string_concat(S2,S1,S3),
        term_string(P, S3)
    ;
        % Otherwise sum them
        P = L1+T1
    ),
    % The others computations are semantically meaningless
    !.
list_to_polynomial([T], T).
%% Tests:
%% ?- list_to_polynomial([1, x, x^2], P).
%@ P = x^2+x+1.
%% ?- list_to_polynomial([-1, -x, -x^2], P).
%@ P = -x^2-x-1.
%% ?- list_to_polynomial([1, -x, x^2], P).
%@ P = x^2-x+1.
%% ?- list_to_polynomial([x^2, x, 1], P).
%@ P = 1+x+x^2.
%% ?- list_to_polynomial([a,-e], P).
%@ P = -e+a.
%% ?- list_to_polynomial([], P).
%@ false.
%% ?- list_to_polynomial([a], P).
%@ P = a.

%% negate_term(T, T2) is det
%
%  Negate the coeficient of a term and return the negated term.
%
negate_term(T, T2) :-
    term_to_list(T, L),
    %% Ensure there is a coeficient
    term_to_canon(L, L2),
    [N | In] = L2,
    %% (-)/1 is an operator, needs to be evaluated, otherwise
    %% it gives a symbolic result, which messes with further processing
    N2 is -N,
    %% Convert the term back from canonic form
    term_to_canon(L3, [N2 | In]),
    %% Reverse the order of the list, because converting
    %% implicitly reverses it
    reverse(L3, L4),
    term_to_list(T2, L4),
    !.
%% Tests:
%% ?- negate_term(1, In).
%@ In = -1.
%% ?- negate_term(x, In).
%@ In = -x.
%% ?- negate_term(-x, In).
%@ In = x.
%% ?- negate_term(x^2, In).
%@ In = -x^2.
%% ?- negate_term(3*x*y^2, In).
%@ In = -3*y^2*x.

%% scale_polynomial(+P:Polynomial,+C:Constant,-S:Polynomial) is det
%
%  Multiplies a polynomial by a scalar.
%
scale_polynomial(P, C, S) :-
    polynomial_to_list(P, L),
    %% Convert each term to a list
    maplist(term_to_list, L, L2),
    %% Canonize terms
    maplist(term_to_canon, L2, L3),
    %% Append C to the start of each sublist
    maplist(cons(C), L3, L4),
    %% Convert to a list of terms
    maplist(term_to_list, L5, L4),
    %% Simplify the resulting polynomial
    simplify_polynomial_as_list(L5, L6),
    %% Return as a simplified polynomial
    list_to_polynomial(L6, S),
    !.
%% Tests:
%% ?- scale_polynomial(3*x^2, 2, S).
%@ S = 6*x^2.

%% cons(+C:atom, +L:List, -L2:List) is det
%
%  Add an atom C to the head of a list L.
%
cons(C, L, [C | L]).
%% Tests:
%% ?- cons(C, L, L2).
%@ L2 = [C|L].

%% add_polynomial(+P1:polynomial,+P2:polynomial,-S:polynomial) is det
%
%  S = P1 + P2.
%
add_polynomial(P1, P2, S) :-
    %% Convert both polynomials to lists
    polynomial_to_list(P1, L1),
    polynomial_to_list(P2, L2),
    %% Join them
    append(L1, L2, L3),
    %% Simplify the resulting polynomial
    simplify_polynomial_as_list(L3, L4),
    %% Convert back
    list_to_polynomial(L4, S),
    !.
%% Tests:
%% ?- add_polynomial(2, 2, S).
%@ S = 4.
%% ?- add_polynomial(x, x, S).
%@ S = 2*x.
%% ?- add_polynomial(2*x+5*z, 2*z+6*x, S).
%@ S = 8*x+7*z.