/*---------------------------------------------------------*/
  /*                                                         */
  /*                       'SYMDIFF'                         */
  /*                                                         */
  /*  Symbolic Differentiation and Algebraic Simplification  */
  /*                                                         */
  /*    The declarative semantics of Prolog are a powerful   */
  /*  tool for problems involving symbolic manipulation.     */
  /*    In this example the rules for differentiation are    */
  /*  expressed as a set of clauses, where, provided that    */
  /*  the head matches the expression in question, the body  */
  /*  furnishes the appropriate result. This result in turn  */
  /*  is passed to the simplifying clauses which operate on  */
  /*  the same principle, but obviously correspond to a      */
  /*  different set of (simplifying) rules.                  */
  /*                                                         */
  /*---------------------------------------------------------*/


/*  This is the main predicate. */

symdiff :-
   nl,
   write('SYMDIFF : Symbolic Differentiation and '),
   write('Algebraic Simplification. Version 1.0'), nl,
   symdiff1.


  /*  symdiff1 controls the main processing; It prompts  */
  /*  the user for data, differentiates the input,       */
  /*  simplifies the result and then displays it on      */
  /*  the terminal.                                      */

symdiff1  :-
   write('Enter a symbolic expression, or "stop" '),
   read(Exp),
   process( Exp ).


process( stop ) :- !.

process( Exp ) :-
   write('Differentiate w.r.t. :'),
   read(Wrt),
   process(Exp,Wrt),
   symdiff1.

process( Exp ,Wrt ):-
   process( Exp ,Wrt, Sdiff),
   nl,
   write('Differential w.r.t '), write(Wrt),
   write(' is '),
   write(Sdiff), nl, nl.


process( Exp ,Wrt, Sdiff):-
   d(Exp, Wrt, Diff ),
   simplify1(Wrt, Diff, Sdiff ),!.


  /*  The following clauses constitute the differentiation   */
  /*  rules, some of which are recursive. Note too the use   */
  /*  of the cut which ensures that once a special case has  */
  /*  been identified, backtracking will not attempt to find */
  /*  an alternative solution.                               */


%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%
% process(x^(x*log(cos(x))), x)
% process(diff(x^(x*log(cos(x))), x), N).
% N = x^(exp(log(1/x))-1)*(exp(log(1/x))+ - (1)/(x*x)/(1/x)*exp(log(1/x))*x*log(x))*cos(x^exp(log(1/x)))
% d(expr, var, deriv)
%%%%%%%%%%%%%%%%%%%%%%%%%%

d(E, _, E) :- var(E),!.

d(d(Expr,Wrt), X, E):- !,d(diff(Expr,Wrt), X, E).
d(d(Expr,Wrt,D), X, E):- !,d(diff(Expr,Wrt,D), X, E).
d(diff(Expr,Wrt), _X, E):- !, d(Expr, Wrt, E).
d(diff(Expr,Wrt,1), X, E):- !, d(diff(Expr,Wrt), X, E).
% d(diff(Expr,Wrt,2), X, E):- !, d(diff(Expr,Wrt), X, E).
% d(f(Wrt), Wrt, 1):- !.
d( U**C, X, O  ) :- !, d( U^C, X, O  ).

% constant
d(E, _, 0) :- numeric_1(E),!.

%d( C, _X, 0 ):- atomic(C).        /* If C is a constant then */
                                   /* d(C)/dX is 0            */
%d( X, X, 1 ):- !.                 /* d(X) w.r.t. X is 1      */

% the same variable
d(E, V, 0) :- atom(E), E \== V, !.

% other variables
d(E, V, 1) :-  atom(E), E == V,!.


d( U+V, X, A+B ):-                 /* d(U+V)/dX = A+B where   */
   d( U, X, A ),                   /* A = d(U)/dX and         */
   d( V, X, B ).                   /* B = d(V)/dX             */

d( U-V, X, A-B ):-                 /* d(U-V)/dX = A-B where   */
   d( U, X, A ),                   /* A = d(U)/dX and         */
   d( V, X, B ).                   /* B = d(V)/dX             */

d( C*U, X, C*A ):-               /* d(C*U)/dX = C*A where     */
   atomic(C),                    /* C is a numeric_1 or variable */
   C \= X,                       /* not equal to X and        */
   d( U, X, A ), !.              /* A = d(U)/dX               */

d( U*V, X, B*U+A*V ):-           /* d(U*V)/dX = B*U+A*V where */
   d( U, X, A ),                 /* A = d(U)/dX and           */
   d( V, X, B ).                 /* B = d(V)/dX               */

d( U/V, X, (A*V-B*U)/(V*V) ):- /* d(U/V)/dX = (A*V-B*U)/(V*V) */
   d( U, X, A),                /* where A = d(U)/dX and       */
   d( V, X, B).                /*       B = d(V)/dX           */

d( U^C, X, C*A*U^(C-1) ):-       /* d(U^C)/dX = C*A*U^(C-1)   */
   atomic(C),                    /* where C is a numeric_1 or    */
   C\=X,                         /* variable not equal to X   */
   d( U, X, A ).                 /* and d(U)/dX = A           */

d( U^C, X, C*A*U^(C-1) ):-       /* d(U^C)/dX = C*A*U^(C-1)   */
   C = -(C1), atomic(C1),        /* where C is a negated numeric_1 or  */
   C1\=X,                        /* variable not equal to X   */
   d( U, X, A ).                 /* and d(U)/dX = A           */

% (f^g)' = f^(g-1)*(gf' + g'f logf)
d(F^G, V, FF) :- !,
    d(F, V, Fd),
    d(G, V, Gd),
    opt_sub(AA, G, 1),
    opt_mul(BB, G,Fd),
    opt_mul(CC, Gd, F),
    opt_mul(DD, CC, log(F)),
    opt_sum(EE, BB, DD),
    opt_mul(FF, F^(AA), EE).

d( sin(W), X, Z*cos(W) ):-       /* d(sin(W))/dX = Z*cos(W)   */
   d( W, X, Z).                  /* where Z = d(W)/dX         */

d( exp(W), X, Z*exp(W) ):-       /* d(exp(W))/dX = Z*exp(W)   */
   d( W, X, Z).                  /* where Z = d(W)/dX         */

d( log(W), X, Z/W ):-            /* d(log(W))/dX = Z/W        */
   d( W, X, Z).                  /* where Z = d(W)/dX         */

d( cos(W), X, -(Z*sin(W)) ):-    /* d(cos(W))/dX = Z*sin(W)   */
   d( W, X, Z).                  /* where Z = d(W)/dX         */


d(_, _, _):-!,fail.

% functions
% (h(g))' = h'(g) * g'
d(sin(E), V, M) :-  d(E,V, Ed), opt_mul(M, Ed, cos(E)).
d(cos(E), V, K) :-  d(E,V, Ed), opt_mul(M, Ed,sin(E)), opt_sub(K, 0, M).

d(tan(E), V, F) :-
    d(E,V, Ed),
    opt_mul(A, Ed, 2),
    opt_mul(B, 2, E),
    opt_mul(C, 2, cos(B)),
    opt_sum(D, C, 1),
    opt_div(F, A, D).

/*
% sum
d(E1 + E2, V, W) :-
    d(E1, V, Ed1),
    d(E2, V, Ed2),
    opt_sum(W, Ed1, Ed2).

% subtraction
d(E1 - E2, V, W) :- !,
    d(E1, V, Ed1),
    d(E2, V, Ed2),
    opt_sub(W, Ed1, Ed2).

% (fg)' = f'g + fg'
d(E1 * E2, V, W) :-
    d(E1, V, Ed1),
    d(E2, V, Ed2),
    opt_mul(L, Ed1, E2),
    opt_mul(P, E1, Ed2),
    opt_sum(W, L, P).


% (f/g)' = (f'g - fg') / (g*g)
d(E1 / E2, V, E) :-
    d(E1, V, Ed1),
    d(E2, V, Ed2),
    opt_mul(A, Ed1,E2),
    opt_mul(B, E1,Ed2),
    opt_mul(C, E2,E2),
    opt_sub(D, A, B),
    opt_div(E, D,C).
*/


d(exp(E), V, F) :-
    d(E,V, Ed),
    opt_mul(F, Ed, exp(E)).

d(log(E), V, F) :-
    d(E,V, Ed),
    opt_div(F, Ed, E).


  /*  The following clauses are rules for the simplification  */
  /*  of the result of the differentiation process. For       */
  /*  example, multiples of zero are dropped, terms gathered  */
  /*  together where possible and factorisation is attempted. */
  /*  (It should be noted that in this program only adjacent  */
  /*  terms are candidates for simplification.)               */

simp(_, X, X ):-          /* an atom or numeric_1 is simplified */
    atomic(X), !.

simp(_, X+0, Y ):-        /* terms of value zero are dropped  */
   simp(_, X, Y ).

simp(_, 0+X, Y ):-        /* terms of value zero are dropped  */
   simp(_, X, Y ).

simp(_, X-0, Y ):-        /* terms of value zero are dropped  */
   simp(_, X, Y ).

simp(_, 0-X, -(Y) ):-     /* terms of value zero are dropped  */
   simp(_, X, Y ).

simp(_, A+B, C ):-        /* sum numbers */
   numeric_1(A),
   numeric_1(B),
   C is A+B.

simp(_, A-A, 0 ).         /* evaluate differences */

simp(_, A-B, C ):-        /* evaluate differences */
   numeric_1(A),
   numeric_1(B),
   C is A-B.

simp(_, _X * 0, 0 ).         /* multiples of zero are zero */

simp(_, 0* _X, 0 ).         /* multiples of zero are zero */

simp(_, 0/_X, 0 ).         /* numerators evaluating to zero are zero */

simp(_, X*1, X ).         /* one is the identity for multiplication */

simp(_, 1*X, X ).         /* one is the identity for multiplication */

simp(_, X/1, X ).         /* divisors of one evaluate to the numerator */

simp(_, X/X, 1 ) :- !.

simp(_, X^1, X ) :- !.

simp(_, _X^0, 1 ) :- !.

simp(_, X*X, X^2 ) :- !.

simp(_, X*X^A, Y ) :-
   simp(_, X^(A+1), Y ), !.

simp(_, X^A*X, Y ) :-
   simp(_, X^(A+1), Y ), !.

simp(_, A*B, X ) :-       /* do evaluation if both are numbers  */
   numeric_1( A ),
   numeric_1( B ),
   X is A*B.

simp(_, A*X+B*X, Z ):-    /* factorisation and recursive */
   A\=X, B\=X,          /* simplification              */
   simp(_, (A+B)*X, Z ).

simp(_, (A+B)*(A-B), X^2-Y^2 ):-   /* difference of two squares */
   simp(_, A, X ),
   simp(_, B, Y ).

simp(_, X^A/X^B, X^C ):-     /* quotient of numeric_1 powers of X */
   numeric_1(A), numeric_1(B),
   C is A-B.

simp(_, A/B, X ) :-       /* do evaluation if both are numbers  */
   numeric_1( A ),
   numeric_1( B ),
   X is A/B.

simp(_, A^B, X ) :-       /* do evaluation if both are numbers  */
   numeric_1( A ),
   numeric_1( B ),
   X is A^B.

simp(_, W+X, Q ):-        /* catchall */
   simp(_, W, Y ),
   simp(_, X, Z ),
   ( W \== Y ; X \== Z ), /* ensure some simplification has taken place */
   simp(_, Y+Z, Q ).

simp(_, W-X, Q ):-        /* catchall */
   simp(_, W, Y ),
   simp(_, X, Z ),
   ( W \== Y ; X \== Z ), /* ensure some simplification has taken place */
   simp(_, Y-Z, Q ).

simp(_, W*X, Q ):-        /* catchall */
   simp(_, W, Y ),
   simp(_, X, Z ),
   ( W \== Y  ; X \== Z ), /* ensure some simplificaion has taken place */
   simp(_, Y*Z, Q ).

simp(_, A/B, C ) :-
   simp(_, A, X ),
   simp(_, B, Y ),
   ( A \== X ; B \== Y ),
   simp(_, X/Y, C ).

simp(_, X^A, C ) :-
   simp(_, A, B ),
   A \== B,
   simp(_, X^B, C).

simp(_, X, X ).           /* if all else fails... */


%numeric_1(A) :- integer(A).

numeric_1(A) :- number(A),!.
numeric_1(A) :- notrace(catch((_ is A),_,fail)).

%% derivatives.pl -- Chapter 12, Section 12.2.3
%% Last Modified: 2/4/14
%% This is a Prolog program which computes derivatives
%% of polynomials.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTING DERIVATIVES of POLYNOMIALS
% For the purpose of this program polynomials are functions
% constructed from a variable x, integers, function symbols
% +, -, * and ^, and parentheses. Exponentiation X^N is
% defined for nonnegative integer N.
% The derivative is only computed once. Multiple calls can lead
% to an infinite loop.

%:-use_module(library(lists)). % not needed for SWI Prolog

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% derivative(F,DF) computes a derivative DF of a polynomial F.
% The result is given in canonical form.
% type F, DF - polynomials.
% mode derivative(+,-)


derivative(F,DF) :-
    XX = x,
    must(d(F,XX,G)),
    must(simplify1(XX,G,DF)),!.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% der(XX,F,DF) computes a derivative DF of a polynomial F.
% type F, DF - polynomials.
% mode der(XX,+,-).

/*
der(_XX,N,0) :-
    numeric_1(N).
der(XX,XX,1).
der(XX,-XX,-1).
der(_XX,_F^0,0).
der(XX,F^N, N*F^M*DF) :- numeric_1(N),
    N>=1,
    M is N-1,
    der(XX,F,DF).
der(XX,E1*E2, E1*DE2 + E2*DE1) :-
    der(XX,E1,DE1),
    der(XX,E2,DE2).
der(XX,E1+E2, DE1+DE2) :-
    der(XX,E1,DE1),
    der(XX,E2,DE2).
der(XX,E1-E2,DE1-DE2) :-
    der(XX,E1,DE1),
    der(XX,E2,DE2).
der(XX,-(E),DE) :-
    der(XX,(-1)*E,DE).
*/
der(XX,I,O):- d(I,XX,O).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%To simplify3a the resulting derivative we represent polynomials
%as lists of the form [...[Ci,Ni]...] where [Ci,Ni] corresponds
%to the term Ci*x^Ni. Polynomials, written in this form will be
%called l-polynomials. The derivative DF of F computed
%by der(XX,F,DF) will be translated into equivalent l-polynomial,
%written in canonical form and transformed back into the term form.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% simplify3a(x,P,SP) iff SP is a canonical form of a polynomial P
% type: P, SP - polynomials
% mode simplify3a(+,+,?)

simplify3a(XX,P,SPO) :-
  simplify1(XX,P,SP),!,
  (P \=@=SP -> simplify3a(XX,SP,SPO) ; SPO=SP).

simplify1(XX,P0,SP) :-
    must(simp(_,P0,P)),
    freeze(XX,atom(XX)), % dmiles
    must(transform(XX,P,TP)),
    must_or_fake(add_similar(TP,S)),
    must_or_fake(remove_zero(S,S1)),
    must_or_fake(poly_sort(S1,S2)),
    (back_to_terms(XX,S2,SP)->true;S2=SP).
%simplify1(_,SP,SP).

must_or_fake(P):- must( P=..[_,A,B]),!,(call(P)->true;A=B).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transform(x,T,[[C1,N1],...,[Ck,Nk]]) implies that
% T = C1*x^N1 +...+ Ck*x^Nk.
% type: T - polynomial, C - integer, N - nonnegative integer.
% mode: transform(XX,+,?).

transform(_XX,C,[[C,0]]) :-
    integer(C),!.
transform(XX,YY,[[1,1]]):- XX=@=YY,!.
transform(XX,-YY,[[-1,1]]):- XX=@=YY,!.
transform(XX,A*B,L) :-
    transform(XX,A,L1),
    transform(XX,B,L2),
    multiply_l(L1,L2,L).
transform(XX,A0^N0,L) :-
     simplify3a(XX,N0,N),
     simplify3a(XX,A0,A),
     transform(XX,A,L1),
     get_degree(L1,N,L),!.
transform(XX,A+B,RT) :-
     transform(XX,A,RA),
     transform(XX,B,RB),
     append(RA,RB,RT).
transform(XX,A-B,RT) :-
     transform(XX,A,RA),
     transform(XX,-B,RB),
     append(RA,RB,RT).
transform(XX,Neg,R) :- unnegative(Neg,A),!,
     transform(XX,(-1)*A,R).
transform(_XX,X,X).

unnegative(Neg,A):- Neg = -(A),!.
% unnegative(Neg,A):- A<0, A is -Neg.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiply_l(L1,L2,L) L is the product of L1 and L2
% type: L1,L2,L - l-polynomials
% mode: multiply_l(+,+,-)

multiply_one(_,[],[]).
multiply_one([C1,N1],[[C2,N2]|R],[[C,N]|T]) :-
    notrace(catch((C is C1*C2, N is N1+N2),E,(wdmsg(E),fail))),
    multiply_one([C1,N1],R,T).

multiply_l([],_,[]).
multiply_l([H1|T1],T2,T) :-
    multiply_one(H1,T2,R1),
    multiply_l(T1,T2,R2),
    append(R1,R2,T).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get_degree(L1,N,L) L is L1^N
% type: L1,L - l-polynomials, N - nonnegative integer
% mode: get_degree(+,+,-)

get_degree(_L1,0,[[1,0]]).
get_degree(L1,1,L1).
get_degree(L1,N,L) :-
    number(N),
    N > 1,
    M is N-1,
    get_degree(L1,M,L2),
    multiply_l(L1,L2,L).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add_similar(L1,L2) implies that L2 is the result of
% adding up similar terms of L2.
% type: L1 and L2 are l-polynomials.
% mode: add_similar(+,?)

add_similar([],[]).
add_similar([[C1,N]|T],S) :-
    append(A,[[C2,N]|B],T),
    C is C1+C2,
    append(A,[[C,N]|B],R),
    add_similar(R,S).

add_similar([[C,N]|T],[[C,N]|ST]) :-
    not_in(N,T),
    add_similar(T,ST).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% not_in(N,L) holds if an l-polynomial L does not contain
% a term with degree N
% type: N - nonnegative integer, L l-polynomial
% mode: not_in(+,+)

not_in(_,[]).
not_in(_,[[]]).
not_in(N,[[_,M]|T]) :-
     is_neq(N,M),
     not_in(N,T).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% remove_zero(L1,L2) iff list L2 is obtained from list L1
% by removing terms of the form [0,N].
% type: L1,L2  l-polynomilas
% mode: remove_zero(+,?)

remove_zero([],[]).
remove_zero([[0,_]|T],T1) :-
    remove_zero(T,T1).
remove_zero([[C,N]|T],[[C,N]|T1]) :-
    is_neq(C,0),
    remove_zero(T,T1).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% poly_sort(L1,L2) sorts l-polynomial L1 in decreasing order
% of exponents and stores the result in L2.
% type: L1,L2  l-polynomilas
% mode: poly_sort(+,-)

poly_sort([X|Xs],Ys) :-
     poly_sort(Xs,Zs),
     insert(X,Zs,Ys).
poly_sort([],[]).

insert(X,[],[X]).
insert(X,[Y|Ys],[Y|Zs]) :-
     greater(Y,X),
     insert(X,Ys,Zs).
insert(X,[Y|Ys],[X,Y|Ys]) :-
     greatereq(X,Y).

greater([_,N1],[_,N2]) :-
     N1 > N2.

greatereq([_,N1],[_,N2]) :-
     N1 >= N2.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% back_to_terms(X,L,T) computes a term-representation T of L.
% L must be an l-polynomial sorted in decreasing order of
% exponents and not containing terms of the form [0,N].
% type: L l-polynomial, T polynomial.
% mode: back_to_terms(+,+,?)

back_to_terms(XX,[[1,1]],XX).
back_to_terms(XX,[[-1,1]],-XX).
back_to_terms(_XX,[[C,0]],C) :-
     integer(C).
back_to_terms(XX,[[1,N]],XX^N):-
      is_neq(N,0),
      is_neq(N,1).
back_to_terms(XX,[[-1,N]],-XX^N):-
      is_neq(N,0),
      is_neq(N,1).
back_to_terms(XX,[[C,1]],C*XX):-
      is_neq(C,1),
      is_neq(C,-1).
back_to_terms(XX,[[C,N]],C*XX^N):-
     is_neq(N,1),
     is_neq(C,1),
     is_neq(C,-1).
back_to_terms(XX,L,F+T) :-
     append(L1,[[C,N]],L),
     is_neq(L1,[]),
     C > 0,
     back_to_terms(XX,L1,F),
     back_to_terms(XX,[[C,N]],T).
back_to_terms(XX,L,F-T) :-
     append(L1,[[C,N]],L),
     is_neq(L1,[]),
     C < 0,
     abs(C,AC),
     back_to_terms(XX,L1,F),
     back_to_terms(XX,[[AC,N]],T).

eq(X,X).
is_neq(X,Y) :-
    \+ eq(X,Y).
abs(X,X) :-
    integer(X),
    X >= 0.
abs(X,Y) :-
    integer(X),
    X < 0,
    Y is (-1)*X.


% dee(+variable, +expression, -derivative)
dee(X, X, 1).
dee(X, Y, 0) :- atomic(Y), X \= Y.

dee(X, U + W, DU + DW) :- aux(X, U, W, DU, DW).
dee(X, U - W, DU - DW) :- aux(X, U, W, DU, DW).

dee(X, U * W, U * DW + W * DU) :- aux(X, U, W, DU, DW).

dee(X, -U, -DU) :- dee(X, U, DU).

aux(X, U, W, DU, DW) :-
    dee(X, U, DU),
    dee(X, W, DW).

% Rather than a larger class of symbolic expressions,
% one might be tempted to implement a simplifier :-)
%
% | ?-
/*
 dee(x, (x + c) * (x - 1), D),
      dee(x, x * (x - 1) + c * (x - 1), E),
      dee(x, x * x - x + c * x - c, F).
*/
%
% D = (x+c)*(1-0)+(x-1)*(1+0)
% E = x*(1-0)+(x-1)*1+(c*(1-0)+(x-1)*0)
% F = x*1+x*1-1+(c*1+x*0)-0 ? ;
%
% no


opt_sum(X, X, 0).
opt_sum(Y, 0, Y).

opt_sum(W, X, Y) :-
    numeric_1(X),
    numeric_1(Y),
    W is X + Y.

opt_sum(2*X, X, Y) :-
    X == Y.

opt_sum(X+Y, X, Y).

%%%%%%%%%%%%%%%%%%%%%%%%%%
opt_sub(X, X, 0).
opt_sub(-Y, 0, Y).

opt_sub(W, X, Y) :-
    numeric_1(X),
    numeric_1(Y),
    W is X - Y.

opt_sub(0, X, Y) :-
    X == Y.

opt_sub(X-Y, X, Y).

%%%%%%%%%%%%%%%%%%%%%%%%%%
opt_mul(0, 0, _).
opt_mul(0, _, 0).

opt_mul(Y, 1, Y).
opt_mul(X, X, 1).

opt_mul(Y, X, Y) :-
    X == 1.

opt_mul(W, X, Y) :-
    numeric_1(X),
    numeric_1(Y),
    W is X * Y.

opt_mul(X*Y, X, Y).

%%%%%%%%%%%%%%%%%%%%%%%%%%
opt_div(_,_,0) :-
    throw(error(evaluation_error(zero_divisor),(is)/2)).

opt_div(0, 0, _N).

opt_div(X, X, 1).

opt_div(1, X, Y) :-
    X == Y.

opt_div(W, X, Y) :-
    numeric_1(X),
    numeric_1(Y),
    W is X / Y.

opt_div(X/Y, X, Y).