Program for the symbolic simplification of algebra



                            ECOMPUTER  ALGEBRA  IN  PROLOGF

                                         by
                                  Sergio Vaghi (*)


                                     April 1987


          (*)  ESTEC,
               P.O. Box 299,
               2200 AG Noordwijk (ZH),
               The Netherlands.

               tel. +31  1719 83453


                                                                     Page 2


               Many   applications   which   are  usually  associated  with
          Artificial Intelligence do not necessarily need an AI language to
          be implemented.  Expert systems,  for  example, can be, and often
          are, written in traditional procedural languages such as FORTRAN,
          Pascal and Basic.

             But are there applications for which the use of  AI  languages
          such  as  LISP  and  Prolog  is  mandatory,  or  at  least highly
          advisable ?  The answer is obviously affirmative, otherwhise such
          languages wouldn't exist, and some nice examples can be provided.
             One of them is  symbolic  mathematical computation or computer
          algebra.  This can be described as the ability of the computer to
          perform  algebraic  manipulations,  as   opposed   to   numerical
          calculations.
          A  computer  performing symbolic differentiation, for example, is
          able to find that the derivative of  y  = xS 2T with respect to x is
          dy/dx = 2x  for any value of x, instead of simply calculating the
          numerical value of the same derivative for a given value of x.

             Computer algebra, at various levels of complexity, is  already
          more  than twenty years old, with some well established programs,
          such as MACSYMA and REDUCE,  and important results at its credit,
          particularly in the fields of  celestial  mechanics  and  quantum
          electrodynamics.
          Algebraic  manipulators  have  been  written, to be sure, both in
          procedural  and  declarative   languages,   but   the  use  of  a
          declarative language has the advantage  of  clarity,  conciseness
          and  elegance,  while  procedural languages were mostly used when
          little else was available.


             In  this   article   I   describe   a   program  for  symbolic
          differentiation in Prolog, which provides a good example  of  the
          power  of  the language and, at the same time, will introduce you
          to the fascinating world of computer algebra.
          Symbolic differentiation is,  in  a  sense,  already  part of the
          Prolog folklore.  An example is given in the book by Clocksin and
          Mellish (ref.1) and another in the collection of Prolog  programs
          by  Coelho,  Cotta  and  Pereira  (ref.2).    But  the  former is
          incomplete (and in the  first  edition  of  the book contained an
          error) and the latter is both incomplete and partially incorrect.
          I have tried to offer something that, I hope, comes a bit  closer
          to   a   reasonably   efficient  symbolic  differentiator.    Its
          application is limited to algebraic functions, but you can easily
          extend it to  include  logarithmic,  trigonometric and hyperbolic
          functions thanks to the incremental nature of the Prolog code.

             While this article is in no way a  comprehensive  introduction
          to  either Prolog or computer algebra, I hope it will give you an
          impression of the many  nice  characteristics of the language and
          perhaps  stimulate  your  interest  in  symbolic  computation  by
          computer.


                                                                     Page 3


           EGetting started.F


             If you are not already familiar with Prolog,  a  good  way  to
          start  is  to read the book 4Programming in Prolog5 by Clocksin and
          Mellish (ref.1).  It describes the basic elements of the language
          and introduces the so-called  Edimburgh  syntax, now used in most
          Prolog implementations.  The book has become a de facto  standard
          of  sorts for Prolog and is highly readable.  The many errors and
          misprints which plagued  the  first  edition  have  now also been
          corrected.
          But you can't learn a new programming language  without  actually
          writing  some program.  For this you need a Prolog interpreter or
          compiler.  Several are available on  the IBM-PC market, some at a
          reasonable price.
          A Prolog interpreter has also been placed in the  public  domain.
          It is called PD Prolog and is a good entry point for newcomers to
          the  language,  who  wish  to get a feeling for the possibilities
          offered by Prolog without having to buy a commercial product.
          I  have  chosen   it   to   develop   the  program  for  symbolic
          differentiation . This has imposed  some  constraints,  since  PD
          Prolog  has a number of limitations with respect to more advanced
          implementations  (for   example,   only   integer  arithmetic  is
          supported), but will give you the possibility to run the  program
          on your IBM-PC or compatible once you have obtained a copy of the
          interpreter   (version   1.9  or  higher)  from  Automata  Design
          Associates,  1570  Arran  Way,  Dresher,   PA  1905,  or  from  a
          distributor of public domain software.
          The program will also run with little or  no  modification  under
          other Prolog implementations following the Edimburgh syntax.


                                                                     Page 4


           EMinimum program.F


             The  rules  of  differentiation  for  algebraic  functions are
          listed in Fig.1.   Before  translating  them  into Prolog code we
          have to define those mathematical operators which are not already
          provided  by  our  Prolog   implementation.      They   are   the
          exponentiation, identified by the symbol '^' (so that x^2 means x
          to  the  power 2), the minus sign as unary operator (as in -x, as
          opposed to the minus sign as  binary operator, as in x-y) denoted
          by ~, and the natural logarithm, ln.
          They are defined as follows
                                       ?-op(11, yfx,  '^').
                                       ?-op( 9,  fx,  '~').
                                       ?-op(11,  fx, 'ln').

          I will not go into the details of these  definitions,  which  are
          explained  in the book by Clocksin and Mellish.  Just assume that
          the above operators are  now  available,  together with the other
          mathematical operators +, -, * and /.

             Writing the program is now fairly simple.   Let's  agree  that
          the Prolog structure
                                d(Y,X,DY).
          means  'DY  is the derivative of Y with respect to X'.  The rules
          of differentiation of Fig.1 can  now  be coded as Prolog clauses.
          The first rule, for example, is coded as
                                                    d(X,X,1).
          meaning that for any X the derivative of X with respect to  X  is
          equal to 1.  The second rule is coded as
                   d(C,X,0):- atomic(C), C \= X.
          which means that the derivative with respect to X of any quantity
          C which is either an atom - that is a non numerical constant - or
          an integer is always 0, provided that C is different from X.
             All  the  rules  of Fig.1 are similarly translated into Prolog
          clauses to form the 4minimum program 5 of Fig.2.
             To find the derivative of  y  =  xS 2T  with respect to x you now
          simply enter
                                                  d(x^2, x, DY).
          after the Prolog prompt '?-'.
          In the Prolog parlance one says that the variables Y and  X  have
          been  4instantiated5  to,  respectively, x^2 and x, and the program
          has been requested to  determine  the  value  of the variable DY.
          Remember that in Prolog any name beginning with a capital  letter
          is  assumed  to  be  a  variable,  which may be instantiated to a
          constant beginning with a lower-case letter.
             Note that  writing  the  program  simply  meant re-writing the
          rules of differentiation in the sintax of  Prolog,  after  having
          arbitrarily defined the meaning of the structure d(Y,X,DY).
          This  is  what is meant when it is said that in logic programming
          the specifications 4are5 the program.
             The  program  is  now  able  to  find  the  derivative  of any
          algebraic function, from the most simple  to  the  most  complex.


                                                                     Page 5


          While  in a procedural language one must tell the computer 4how5 to
          proceed at each step  of  the  differentiation process, in Prolog
          one simply states the basic rules,  and  the  computer  does  the
          rest.
             Too  good  to be true? well, yes and no.  The program of Fig.2
          is  indeed  able  to  handle  the  derivative  of  any  algebraic
          function, but the result it produces  may  not be in the form you
          expect to see it.
          Let's try the very simple example of the derivative  of  y  =  xS 2T
          with  respect  to  x.    The  result  is shown in Fig.3.  It will
          probably take you a  few  seconds  to  realize that the resulting
          expression, DY, is indeed equivalent to 2x.
          The problem is that the program contains only the 4basic5 rules  of
          differentiation.   Other derived rules are not included.  In this
          particular case a human  would  immediately  use the derived rule
          stating that, when c is a constant, the derivative of y = xS cT with
          respect to x is  c xS (c-1)T.  This rule can  be  derived  from  the
          last rule in Fig.1, but it is so much part of the 4heuristics5 used
          by humans that it is advantageous to include it explicitly in the
          program.
             Derived rules of this type are added to the minimum program to
          produce  the program of Fig.4.  Note that the new clauses include
          the symbol '!', called 4cut5  in  Prolog.    Its effect is to avoid
          backtracking, that is the search of an alternative solution  once
          one  has been found.  Indeed the new clauses are special cases of
          those of the minimum program.   If backtracking were allowed, one
          would obtain - in certain cases - two solutions, one from  a  new
          clause and one from an old one.  The former is more specific and,
          therefore,  simpler  and  should be retained.  The latter becomes
          useless and need not be shown.  The cut obtains this effect.

             If you run the new program for    y = xS 2T you will now find the
          expression of Fig.5.  It is already an improvement  with  respect
          to  the  previous result, but still unsatisfactory.  The computer
          doesn't seem to know that 2*1 = 2 and 2-1 = 1. In fact it can not
          simplify an expression.

             Simplification is indeed a  major problem in computer algebra.
          I'll come back to it later.
          First I want to  discuss  two  other  nice  features  of  Prolog.
          Starting from our definition of the structure  d(Y,X,DY). we have
          obtained the derivative of  xS 2T with respect to x by instantiating
          Y  to  x^2  and X to x.  Prolog has then derived the value of the
          variable DY which satisfies the rule, that is  (2*1)*(x^(2-1)).
          Nothing, in principle, forbids us to instantiate X to x and DY to
          (2*1)*(x^(2-1))  and let  Prolog  determine  the value of Y which
          satisfies the rule.  This is shown in Fig.6, where indeed  Prolog
          finds for Y the value x^2.
          This  is  very important because it demonstrates that the program
          of Fig.4 can not only be used for its original pourpose, symbolic
          differentiation of algebraic functions,  but also for the inverse
          operation, that is symbolic integration of the same functions.
          This is the property of logic programming  called  4invertibility5,
          or  the ability to perform inverse computations.  Prolog programs


                                                                     Page 6


          which,  like  the   one   in   Fig.4,   strictly  adhere  to  the
          prescriptions of logic programming allow  invertibility.    Once,
          however,  procedural elements are added to the program, such as a
          numerical  calculation,  the  property   is  lost,  because  such
          elements are not themselves invertible.
          Also note that, for integration, the Prolog variable DY has  been
          instantiated  to    (2*1)*(x^(2-1))    ,  using  exactly the same
          notation used by the  program  when calculating the derivative of
          x^2.   Using the more familiar notation  2*x    would  not  work,
          because  the  computer does not know that the two expressions are
          equivalent.  Even with these restrictions invertibility remains a
          powerful and useful characteristic of the language.

             The second nice feature of Prolog which should be mentioned is
          the incremental nature of its code.
          Suppose  that  you  now  decide  to  extend  the  program  to the
          differentiation of, say, the trigonometric  functions  sinus  and
          cosinus.  You don't have to re-write the entire program for that.
          Simply  add sin and cos to the list of operators already defined,
          and the differentiation rules for  sinus  and cosinus you want to
          use to  the  rules  already  present  in  the  program.  Existing
          definitions and clauses need not be changed.
          Once  sinus  and  cosinus  have  been introduced, you may wish to
          introduce their inverse functions,  then the hyperbolic functions
          and their inverses and so on.  At each step you  simply  add  the
          operators  and  the relevant rules to the program obtained at the
          previous step.    Contrary  to  what  all  too  often  happens in
          procedural languages, adding new features  to  a  Prolog  program
          does not involve the complete re-writing of the code.


                                                                     Page 7


           EThe politics of simplification.F


             I  come  back  now to the problem of simplification.  This has
          always been a major problem  in  computer algebra.  The preferred
          form af an expression may depend on several factors, such as  the
          context  in  which it is derived (in relativity, for example, one
          would prefer to obtain the expression
                                                  E = m cS 2T
          rather than the equivalent expression
                                1 / cS 2T = m / E ),
          the type  of  application  (purely  symbolic  manipulation,  or a
          mixture of symbolic manipulation and numerical  computation)  and
          personal   taste   of  the  user.    Moreover,  when  complicated
          expressions are involved,  it  is  difficult  for the computer to
          recognize that an expression - or a part of it -  is  identically
          equal to zero and can consequently be suppressed.

             Without  entering  into  too many details, it is indicative of
          the complexity of the problem the  fact that already 15 years ago
          Joel Moses of the MIT could write a short dissertation  on  the  4
          politics of simplification5 (ref.3).
          The  criterion  considered  was  the degree to which an algebraic
          manipulation system would force  the  user  to write (or read) an
          expression in a form acceptable to the system.
          A 4radical5 system would insist in completely changing the form  of
          the expression before being able to handle it.  Our program which
          pretends  to  have  2*x written as (2*1)*(x^(2-1))  to be able to
          find its integral can be taken as an example of radical system.
          4Conservative5 systems are those which  would leave to the user the
          burden of simplifying the expressions they found.   Our  program,
          when  used  to  find  the  derivative of xS 2T with respect to x, is
          conservative in this sense,  since  it  leaves to you to simplify
          the expression (2*1)*(x^(2-1)).  Note that,  in  differentiation,
          the program is not a radical one, because you can enter xS 2T in its
          simplest form x^2.
          4Liberal5  systems  will  perform  some  simplifications, but leave
          others to the user.  4New  left5  systems are variations of the old
          radical systems, with some added features to be selected  by  the
          user.
          4Catholic5  systems  are  those  which  take  advantage of the good
          points of the  other  systems,  and  adopt  them depending on the
          application at hand.
          I recommend Moses' article  to  anybody  interested  in  computer
          algebra.

             The  simplification  program I have written is listed in Fig.7
          and  can  be  classified  as  liberal,  in  that  it  imposes  no
          constraint in the form of  the  input expression, and will do its
          best to simplify the resulting expression so that  it  will  look
          sufficiently familiar to you.  It is meant for use in association
          with the differentiation program, but can also be used as a stand
          alone  program  to simplify algebraic expressions.  The fact that


                                                                     Page 8


          it is much longer than the differentiation program shows that, if
          not more complex, it is by far more tedious to write.  In fact it
          essentially contains most of the rules of elementary algebra.

             I will  not  discuss  the  program  in  detail.    It contains
          comments which should help you in understanding  the  code.    It
          also includes some useful features, such as the capability of
          identifying  divisions  by  zero and indefinite forms of the type
          0/0, and gives an example of  the  use of the Prolog recursion to
          calculate the n-th power of a positive or negative integer.
          Since PD  Prolog  uses  integer  arithmetic  only,  divisions  of
          integer  numbers are never actually performed, to avoid round-off
          errors.  Prolog purists will also note the unusually large number
          of cuts in the code.  This  is to force the program to accept the
          first solution it  finds  at  each  step  of  the  simplification
          process,  consistent with the fact that clauses in each group are
          listed in order of increasing generality.
             The simplification is performed by entering
                            simplify(Y,Z).
          where Y is the expression  to  be simplified and Z the simplified
          expression.  If you are not happy with the result you can  always
          try a two pass simplification, that is
                            simplify(Y,T), simplify(T,Z).
          and so on.  In practice I found that two passes are sufficient in
          most cases.  If you enter
                                     s(Y,Z).
          two passes will be automatically performed.

             To  use  the  program  in conjunction with the differentiation
          program applied to the function  y = xS 2T , simply enter
                            Y = x^2 ,
                            d( Y, x, Z),
                            s( Z, DY).
          after the Prolog prompt, as shown  in Fig.8, and you will finally
          obtain the result in the familiar form DY = (2*x).
          This is, of course, the result we were after.   The  price  paid,
          however, is the loss of the invertibility property.  The program,
          or  more  precisely  the combination of the two programs, will no
          longer be able to  integrate  2*x  to  find  x^2.  The procedural
          elements involved in the simplification have  made  invertibility
          impossible.

             How  good  are the differentiation and simplification programs
          when applied to functions less trivial than y = xS 2T ?
          I have tested the programs  using  the problems found in a manual
          of calculus for college level students.  A few examples are found
          in Fig.9.  Both the intermediate, nonsimplified, results and  the
          simplified  ones  are  printed  in  all  cases, to illustrate the
          advantage  of   simplification   in   making   the  results  more
          understandable.  When complicated functions are involved you will
          still have to find your way in  the  resulting  expression,  both
          because  the program includes many brackets to avoid ambiguities,
          and because certain  simplifications  seen  by  the  user are not
          identified by the program.


                                                                     Page 9


          In many applications simbolic  differentiation  is  coupled  with
          numerical calculations, that is a symbolic derivative is fed to a
          numerical  routine  which  calculates  its  values  for the given
          values  of  parameters  and   independent   variable.    In  such
          applications  the  fact  that  the  expression  is  not   totally
          simplified  is  usually  immaterial  since  it  is  the numerical
          routine which performs the  calculations  (but beware of possible
          spurious singularities - writing  x or xS 2T/x  is  4not5  numerically
          equivalent for values of x close to 0).
          Again,  thanks  to the incremental nature of the Prolog code, you
          may also decide to add more  clauses  to the program to take into
          account a larger number of simplification rules.  In view of  the
          limits  of  PD  Prolog the programs of Figs. 4 and 7 go, I think,
          some way towards  a  fairly satisfactory symbolic differentiation
          system.


                                                                    Page 10


           EApplications.F


             A program for symbolic differentiation is a powerful tool  for
          several  applications.    Examples which come immediately to mind
          are the calculation of  the  partial  derivatives of functions of
          several variables and  the  determination  of  the  higher  order
          derivatives of functions of one variable.
          The chain rule
                          dy/dx = dy/du * du/dx
          can  also  be  directly  implemented, and derivatives of implicit
          functions are easily obtained.
          Fig.10 shows with a few examples how to proceed in these cases.

             The most interesting  applications  are,  however, those which
          combine symbolic differentiation with numerical calculation.
          A nice example is the determination of the recursive  formula  of
          the   Newton-Raphson   method  for  solution  of  real  algebraic
          equations of the type f(x) = 0, that is
           xSk+1T = xSkT - (f(xSkT) / f'(xSkT))
          where xSkT and xSk+1T are successive approximations to the solution.
             The  Newton-Raphson  method   has   the   advantage  of  rapid
          convergence (although it may diverge in certain cases),  but  the
          disadvantage  that  the first derivative of the function y = f(x)
          must be  calculated.    If  f(x)  is  a  complicated function its
          derivative might be tedious to  calculate  by  hand,  and  errors
          could  easily  slip  in.    The  alternative  is to calculate the
          derivative of the function or, better, the entire right hand side
          of the recursive  expression  symbolically  by computer, and feed
          the result to a numerical routine.  How this can be done with our
          programs is illustrated in Fig.11 for the simple equation
                  xS 3T - a = 0
          which can be used to calculate iteratively the cubic root of a.
          In practice one would, of  course,  use  the  programs  for  more
          complicated  functions,  when  differentiation by hand takes much
          longer than by computer.

             Other  areas  in   which   the   combination   of  a  symbolic
          differentiator and a numerical routine can be very attractive are
          the solution of  ordinary  differential  equations  using  Taylor
          series   and   the   calculation   of   the  Jacobian  matrix  in
          linearization problems.  If you  are familiar with these subjects
          you should have  no  problem  in  using  the  programs  for  such
          applications.


                                                                    Page 11


           EReferencesF


          1.  W.F.Clocksin  and  C.Mellish,  4Programming  in Prolog5 (second
          edition), 1984

          2. H.Coelo,  J.C.Cotta  and  L.M.Pereira,  4How  to  solve it with
          Prolog5 (3rd edition), 1982

          3. J.Moses, 4Algebraic simplification - a guide for the perplexed5,
          in 4Proceedings of the Second Symposium on Symbolic and  Algebraic
          Manipulation5, 1971


          dx/dx = 1

          dc/dx = 0 , c = const.


          d(-u)/dx = -(du/dx)


          d(u+v)/dx = du/dx + dv/dx

          d(u-v)/dx = du/dx - dv/dx


          d(uv)/dx  = u dv/dx + du/dx v

          d(u/v)/dx = (v du/dx - u dv/dx) / vS 2T

          d(uS vT)/dx  = uS vT ( v/u du/dx + dv/dx ln(u) )


  Fig. 1 - Rules of differentiation for algebraic functions.


  /*  definition of operators */


  ?-op(11, yfx,  '^').                 /*  exponentiation     */
  ?-op( 9,  fx,  '~').                 /*  minus sign         */
  ?-op(11,  fx, 'ln').                 /*  natural logarithm  */


  d(X,X,1).

  d(C,X,0) :- atomic(C), C \= X.

  d(~U,X,~D) :- d(U,X,D).


  d(U+V,X,D1+D2) :- d(U,X,D1), d(V,X,D2).


  d(U-V,X,D1-D2) :- d(U,X,D1), d(V,X,D2).


  d(U*V, X, D2*U+D1*V) :- d(U,X,D1), d(V,X,D2).


  d(U/V, X, (V*D1-U*D2)/(V^2) ) :- d(U,X,D1), d(V,X,D2).


  d(U^V, X, U^V*(V*D1/U+D2*ln(U) )) :- d(U,X,D1), d(V,X,D2).


    Fig. 2 - Symbolic differentiation - minimum program


  ?-/*               derivative of x^2 with respect to x            */


                            d( x^2, x, DY).


  DY = ((x ^ 2) * (((2 * 1) / x) + (0 * ln x )))


     Fig. 3 - Derivative of x^2 with respect to x, using the
              minimum program of Fig.2


   /* -----------------------------------------------------------------

                                DIFFSV.PRO

                               (version 1.0)

                          Copyright 1987 S.Vaghi


                Program for the symbolic differentiation of
                             algebraic functions.

                ...........................................


       Example of how to use the program:

           to find the derivative, DY, of the function
                                                         y = x^2
           with respect to x, one can

                   a) simply enter
                                                  d( x^2, x, DY).
                      after the Prolog prompt,
           or
                   b) enter
                                                   Y = x^2 ,
                                                   d( Y, x, DY).
                      after the prompt,
           or
                   c) enter the complete sequence including
                      simplification, that is
                                                    Y = x^2 ,
                                                    d( Y, x, Z),
                                                    s( Z, DY).

           Method c) is always recommended, in which case the
           program is used in combination with SIMPSV.PRO


     ----------------------------------------------------------------- */


                    /*  definition of operators */


  ?-op(11, yfx,  '^').                 /*  exponentiation     */
  ?-op( 9,  fx,  '~').                 /*  minus sign         */
  ?-op(11,  fx, 'ln').                 /*  natural logarithm  */


  d(X,X,1).


  d(C,X,0) :- atomic(C), C \= X.

  d(~U,X,~D) :- d(U,X,D).


          d(C+U,X,D) :- atomic(C), C \= X, d(U,X,D), ! .
          d(U+C,X,D) :- atomic(C), C \= X, d(U,X,D), ! .

  d(U+V,X,D1+D2) :- d(U,X,D1), d(V,X,D2).


          d(C-U,X,~D) :- atomic(C), C \= X, d(U,X,D), ! .
          d(U-C,X, D) :- atomic(C), C \= X, d(U,X,D), ! .

  d(U-V,X,D1-D2) :- d(U,X,D1), d(V,X,D2).


          d(C*U,X,C*D) :- atomic(C), C \= X, d(U,X,D), ! .
          d(U*C,X,C*D) :- atomic(C), C \= X, d(U,X,D), ! .

  d(U*V, X, D2*U+D1*V) :- d(U,X,D1), d(V,X,D2).


          d(1/U,X, ~D/(U^2) ) :- d(U,X,D), ! .
          d(C/U,X, C*DD) :- atomic(C), C \= X, d(1/U,X,DD), ! .
          d(U/C,X,  D/C) :- atomic(C), C \= X, d(U,X,D), ! .

  d(U/V, X, (V*D1-U*D2)/(V^2) ) :- d(U,X,D1), d(V,X,D2).


          d( U^C, X, C*D*U^(C-1) ) :- atomic(C), C \= X, d(U,X,D), ! .
          d(U^~C, X,            Z) :-  d( 1/(U^C), X, Z), ! .
          d( U^(A/B),X, (A/B)*D*U^(A/B-1) ) :- atomic(A), atomic(B),
                                               A \= X, B \= X,
                                               d(U,X,D), ! .
          d(U^~(A/B),X,                  Z) :- d(1/U^(A/B) , X, Z), ! .

  d(U^V, X, U^V*(V*D1/U+D2*ln(U) )) :- d(U,X,D1), d(V,X,D2).


     Fig. 4 - The program for symbolic differentiation


  ?-/*               derivative of x^2 with respect to x            */


                            d( x^2, x, DY).


  DY = ((2 * 1) * (x ^ (2 - 1)))


     Fig. 5 - Derivative of x^2 with respect to x, using the
              program DIFFSV.PRO


  ?-/*          integral of  (2*1)*(x^(2-1))  with respect to x        */


                        d( Y, x, (2*1)*(x^(2-1)) ).


  Y = (x ^ 2)


     Fig. 6 - Using DIFFSV.PRO to calculate the integral of
              (2*1)*(x^(2-1))


   /* -----------------------------------------------------------------
  --

                                 SIMPSV.PRO

                                (version 1.0)

                            Copyright 1987 S.Vaghi


                  Program for the symbolic simplification of
                              algebraic functions.

                  ..........................................


          Example of how to use the program:

              to simplify the expression
                                               (2*1)*(x^(2-1))
              one can

                  a) simply enter
                                               s( (2*1)*(x^(2-1)), Z).
                     after the Prolog prompt,
              or
                  b) enter
                                                Y = (2*1)*(x^(2-1)),
                                                s( Y, Z).
                     after the prompt.

              In both cases a two pass simplification is performed.


  -----------------------------------------------------------------
  --- */


                   /*  definition of operators */

  ?-op(11, yfx,  '^').                 /*  exponentiation     */
  ?-op( 9,  fx,  '~').                 /*  minus sign         */
  ?-op(11,  fx, 'ln').                 /*  natural logarithm  */


               /*  two pass simplification clause  */


               s(X,Y) :- simplify(X,Z), simplify(Z,Y).


        /*  list processing of the expression to be simplified  */


  simplify(X,X) :- atomic(X), ! .

  simplify(X,Y) :- X =..[Op,Z], simplify(Z,Z1), u(Op,Z1,Y), ! .

  simplify(X,Y) :- X =..[Op,Z,W], simplify(Z,Z1),
                   simplify(W,W1),
                   r(Op,Z1,W1,Y), ! .


        /*  simplification clauses for addition  */


  r('+',~X,~X,Z) :- b('*',2,X,W), u('~',W,Z) , ! .
  r('+', X, X,Z) :- b('*',2,X,Z), ! .

  r('+', X,~Y,Z) :- b('-',X,Y,Z), ! .
  r('+',~X, Y,Z) :- b('-',Y,X,Z), ! .
  r('+',~X,~Y,Z) :- b('+',X,Y,W), u('~',W,Z), ! .

  r('+',   X, Y/Z, W) :- integer(X), integer(Y), integer(Z),
                         T is Z*X+Y,
                         b('/',T,Z,W), ! .
  r('+', X/Z,   Y, W) :- integer(X), integer(Y), integer(Z),
                         T is X+Y*Z,
                         b('/',T,Z,W), ! .

  r('+',   X, Y+Z, W) :- b('+',Y,Z,T), b('+',X,T,W), ! .
  r('+', X+Y,   Z, W) :- b('+',X,Y,T), b('+',T,Z,W), ! .

  r('+', X*Y,Z*Y,W) :- b('+',X,Z,T), b('*',Y,T,W), ! .
  r('+', X*Y,Y*Z,W) :- b('+',X,Z,T), b('*',Y,T,W), ! .
  r('+', Y*X,Z*Y,W) :- b('+',X,Z,T), b('*',Y,T,W), ! .
  r('+', Y*X,Y*Z,W) :- b('+',X,Z,T), b('*',Y,T,W), ! .

  r('+',X,Y,Z) :- integer(Y), b('+',Y,X,Z), ! .


        /*  simplification clauses for subtraction  */


  r('-', X,~X,Z) :- b('*',2,X,Z), ! .
  r('-',~X, X,Z) :- b('*',2,X,W), u('~',W,Z), ! .

  r('-', X,~Y,Z) :- b('+',X,Y,Z), ! .
  r('-',~X, Y,Z) :- b('+',X,Y,W), u('~',W,Z), ! .
  r('-',~X,~Y,Z) :- b('-',Y,X,Z), ! .

  r('-',   X, Y/Z, W) :- integer(X), integer(Y), integer(Z),
                         T is X*Z-Y,


                         b('/',T,Z,W), ! .
  r('-', X/Z,   Y, W) :- integer(X), integer(Y), integer(Z),
                         T is X-Y*Z,
                         b('/',T,Z,W), ! .

  r('-',   X, Y-Z, W) :- b('-',Y,Z,T), b('-',X,T,W), ! .
  r('-', X-Y,   Z, W) :- b('-',X,Y,T), b('-',T,Z,W), ! .

  r('-', X*Y, Z*Y, W) :- b('-',X,Z,T), b('*',Y,T,W), ! .
  r('-', X*Y, Y*Z, W) :- b('-',X,Z,T), b('*',Y,T,W), ! .
  r('-', Y*X, Z*Y, W) :- b('-',X,Z,T), b('*',Y,T,W), ! .
  r('-', Y*X, Y*Z, W) :- b('-',X,Z,T), b('*',Y,T,W), ! .


        /*  simplification clauses for multiplication  */


  r('*', X, X,Z) :- b('^',X,2,Z), ! .
  r('*',~X,~X,Z) :- b('^',X,2,Z), ! .
  r('*',~X, X,Z) :- b('^',X,2,W), u('~',W,Z), ! .
  r('*', X,~X,Z) :- b('^',X,2,W), u('~',W,Z), ! .

  r('*',      X, X^(~1), Z) :- b('/',X,X,Z), ! .
  r('*', X^(~1),      X, Z) :- b('/',X,X,Z), ! .

  r('*',   X, 1/X, Z) :- b('/',X,X,Z), ! .
  r('*', 1/X,   X, Z) :- b('/',X,X,Z), ! .
  r('*',   X, 1/Y, Z) :- b('/',X,Y,Z), ! .
  r('*', 1/X,   Y, Z) :- b('/',Y,X,Z), ! .
  r('*',   M, N/X, Z) :- atomic(M), atomic(N),
                         b('*',M,N,S), b('/',S,X,Z), ! .
  r('*', M/X,   N, Z) :- atomic(M), atomic(N),
                         b('*',M,N,S), b('/',S,X,Z), ! .


  r('*',  X, N/Y, Z) :- atomic(N), b('/',X,Y,S), b('*',N,S,Z), ! .
  r('*',N/Y,   X, Z) :- atomic(N), b('/',X,Y,S), b('*',N,S,Z), ! .

  r('*',     X,Y^(~1),Z) :- b('/',X,Y,Z), ! .
  r('*',     X,  X^~Y,Z) :- b('-',Y,1,S), b('^',X,S,T), b('/',1,T,Z), ! .
  r('*',X^(~1),     Y,Z) :- b('/',Y,X,Z), ! .
  r('*',  X^~Y,     X,Z) :- b('-',Y,1,S), b('^',X,S,T), b('/',1,T,Z), ! .

  r('*',  X,X^Y,Z) :- b('+',1,Y,S), b('^',X,S,Z), ! .
  r('*',X^Y,  X,Z) :- b('+',Y,1,S), b('^',X,S,Z), ! .

  r('*',~X,~Y,Z) :- b('*',X,Y,Z), ! .
  r('*', X,~Y,Z) :- b('*',X,Y,W), u('~',W,Z), ! .
  r('*',~X, Y,Z) :- b('*',X,Y,W), u('~',W,Z), ! .

  r('*',Z^~X,Z^~Y,W) :- b('+',X,Y,S), b('^',Z,S,T), b('/',1,T,W), ! .
  r('*',Z^~X, Z^Y,W) :- b('-',Y,X,S), b('^',Z,S,W), ! .
  r('*', Z^X,Z^~Y,W) :- b('-',X,Y,S), b('^',Z,S,W), ! .


  r('*', Z^X, Z^Y,W) :- b('+',X,Y,T), b('^',Z,T,W), ! .
  r('*',X^~Z,Y^~Z,W) :- b('*',X,Y,S), b('^',S,Z,T), b('/',1,T,W), ! .
  r('*', X^Z,Y^~Z,W) :- b('/',X,Y,S), b('^',S,Z,W), ! .
  r('*',X^~Z, Y^Z,W) :- b('/',Y,X,S), b('^',S,Z,W), ! .
  r('*', X^Z, Y^Z,W) :- b('*',X,Y,T), b('^',T,Z,W), ! .

  r('*', X*Y,   Y, Z) :- b('^',Y,2,S), b('*',X,S,Z), ! .
  r('*', Y*X,   Y, Z) :- b('^',Y,2,S), b('*',X,S,Z), ! .
  r('*',   Y, X*Y, Z) :- b('^',Y,2,S), b('*',X,S,Z), ! .
  r('*',   Y, Y*X, Z) :- b('^',Y,2,S), b('*',X,S,Z), ! .

  r('*', X*Y, X*Z, W) :- b('*',Y,Z,S), b('^',X,2,T), b('*',T,S,W), ! .
  r('*', Y*X, X*Z, W) :- b('*',Y,Z,S), b('^',X,2,T), b('*',T,S,W), ! .
  r('*', X*Y, Z*X, W) :- b('*',Y,Z,S), b('^',X,2,T), b('*',T,S,W), ! .
  r('*', Y*X, Z*X, W) :- b('*',Y,Z,S), b('^',X,2,T), b('*',T,S,W), ! .

  r('*',  M, N*X, W) :- atomic(M), atomic(N),
                        b('*',M,N,P), b('*',P,X,W), ! .
  r('*',M*X,   N, W) :- atomic(M), atomic(N),
                        b('*',M,N,P), b('*',P,X,W), ! .

  r('*',   X, Y*Z, W) :- b('*',Y,Z,T), b('*',X,T,W), ! .
  r('*', X*Y,   Z, W) :- b('*',X,Y,T), b('*',T,Z,W), ! .

  r('*',X,Y,Z) :- integer(Y), b('*',Y,X,Z), ! .


        /*    simplification clauses for division
             (division is never actually performed)  */


  r('/', 1, X/Y, Z) :- b('/',Y,X,Z), ! .
  r('/',~1, X/Y, Z) :- b('/',Y,X,W), u('~',W,Z), ! .

  r('/',~X,~Y,Z) :- b('/',X,Y,Z), ! .
  r('/', X,~Y,Z) :- b('/',X,Y,W), u('~',W,Z), ! .
  r('/',~X, Y,Z) :- b('/',X,Y,W), u('~',W,Z), ! .

  r('/',      X, Y^(~1), Z) :- b('*',X,Y,Z), ! .
  r('/', X^(~1),      Y, Z) :- b('*',X,Y,W), b('/',1,W,Z), ! .

  r('/',   X, Y/Z, W) :- b('*',X,Z,T), b('/',T,Y,W), ! .
  r('/', X/Y,   Z, W) :- b('*',Y,Z,T), b('/',X,T,W), ! .

  r('/',     X,Y^(~Z),W) :- b('^',Y,Z,T), b('*',X,T,W), ! .
  r('/',X^(~Z),     Y,W) :- b('^',X,Z,S), b('*',S,Y,T), b('/',1,T,W), ! .

  r('/',X,X^(~Y),Z) :- b('+',1,Y,S), b('^',X,S,Z), ! .
  r('/',X,   X^Y,Z) :- b('-',Y,1,S), b('^',X,S,T), b('/',1,T,Z), ! .
  r('/',X^(~Y),X,Z) :- b('+',Y,1,S), b('^',X,S,T), b('/',1,T,Z), ! .

  r('/',   X^Y,     X,Z) :- b('-',Y,1,S), b('^',X,S,Z), ! .
  r('/',   X^N,X^(~M),Z) :- b('+',N,M,S), b('^',X,S,Z), ! .


  r('/',X^(~N),   X^M,Z) :- b('+',N,M,S), b('^',X,S,T), b('/',1,T,Z), ! .
  r('/',X^(~N),X^(~M),Z) :- b('-',M,N,S), b('^',X,S,Z), ! .
  r('/',   X^N,   X^M,Z) :- b('-',N,M,W), b('^',X,W,Z), ! .

  r('/',X^(~Z),   Y^Z,W) :- b('*',X,Y,S), b('^',S,Z,T), b('/',1,T,W), ! .
  r('/',   X^Z,Y^(~Z),W) :- b('*',X,Y,S), b('^',S,Z,W), ! .
  r('/',X^(~Z),Y^(~Z),W) :- b('/',Y,X,S), b('^',S,Z,W), ! .
  r('/',   X^Z,   Y^Z,W) :- b('/',X,Y,T), b('^',T,Z,W), ! .


        /*  simplification clauses for exponentiation  */


  r('^',X,~1,Y) :- b('/',1,X,Y), ! .

  r('^',X,~Y,Z) :- b('^',X,Y,W), b('/',1,W,Z), ! .

  r('^',X^Y,Z,W) :- b('*',Y,Z,T), b('^',X,T,W), ! .


        /*  catch all clause to cover cases not covered
            by the previous clauses                      */


  r(X,Y,Z,W) :- b(X,Y,Z,W).


        /*  basic rules for the unary operator '~'  */


  u('~', 0, 0) :- ! .
  u('~',~X, X) :- ! .
  u('~', X,~X) :- ! .
  u('~',X^Y,~(X^Y) ) :- ! .


        /*  basic rules of addition  */


  b('+', X, 0, X) :- ! .
  b('+',~X, 0,~X) :- ! .
  b('+', 0, X, X) :- ! .
  b('+', 0,~X,~X) :- ! .
  b('+', X,~X, 0) :- ! .
  b('+',~X, X, 0) :- ! .

  b('+',X,Y,Z) :- integer(X), integer(Y),
                  Z is X+Y, ! .

  b('+',X,Y,X+Y).


        /*  basic rules of subtraction  */


  b('-', X, 0, X) :- ! .
  b('-',~X, 0,~X) :- ! .
  b('-', 0, X,~X) :- ! .
  b('-', 0,~X, X) :- ! .
  b('-',~X,~X, 0) :- ! .
  b('-', X, X, 0) :- ! .

  b('-',X,Y,Z) :- integer(X), integer(Y),
                  Z is X-Y, ! .

  b('-',X,Y,X-Y).


        /*  basic rules of multiplication  */


  b('*', 0, X,0) :- ! .
  b('*', 0,~X,0) :- ! .
  b('*', X, 0,0) :- ! .
  b('*',~X, 0,0) :- ! .

  b('*', 1, X, X) :- ! .
  b('*', 1,~X,~X) :- ! .
  b('*',~1, X,~X) :- ! .
  b('*',~1,~X, X) :- ! .
  b('*', X, 1, X) :- ! .
  b('*',~X, 1,~X) :- ! .
  b('*', X,~1,~X) :- ! .
  b('*',~X,~1, X) :- ! .

  b('*', X/Y,   Y, X) :- ! .
  b('*',   Y, X/Y, X) :- ! .

  b('*',X,Y,Z) :- integer(X), integer(Y),
                  Z is X*Y, ! .

  b('*',X,Y,X*Y).


        /*  basic rules of division  */


  b('/',0,0,'ERROR - indefinite form 0/0') :- ! .   /* indefinite form */
  b('/',X,0,'ERROR - division by 0      ') :- ! .   /* division by 0   */

  b('/',0, X,0) :- ! .
  b('/',0,~X,0) :- ! .


  b('/', X,1, X) :- ! .
  b('/',~X,1,~X) :- ! .

  b('/', 1,X,    1/X) :- ! .
  b('/',~1,X, ~(1/X)) :- ! .

  b('/', X,~1,~X) :- ! .
  b('/',~X,~1, X) :- ! .

  b('/', 1,  ~X,~(1/X)) :- ! .
  b('/',~1,  ~X,   1/X) :- ! .
  b('/', 1, 1/X,     X) :- ! .
  b('/',~1, 1/X,    ~X) :- ! .

  b('/', X,~X,~1) :- ! .
  b('/',~X, X,~1) :- ! .
  b('/',~X,~X, 1) :- ! .
  b('/', X, X, 1) :- ! .

  b('/',X,Y,X/Y).


        /*  basic rules of exponentiation  */


  b('^',1,X,1) :- ! .

                                                    /* indefinite forms */

  b('^',0, 0,'ERROR - indefinite form 0^0       ') :- ! .
  b('^',0,~1,'ERROR - indefinite form 0^~1 = 1/0') :- ! .
  b('^',0,~K,'ERROR - indefinite form 0^~K = 1/0') :- atomic(K), ! .


  b('^',~X,1,~X) :- ! .
  b('^', X,1, X) :- ! .
  b('^', X,0, 1) :- ! .


                                          /* recursive clauses to
                                             calculate the n-th power
                                             of positive and negative
                                             integers                 */

  b('^', X,N,Y) :- integer(X), integer(N),
                   M is N-1, b('^',X,M,Z),
                   Y is Z*X, ! .
  b('^',~X,N,Y) :- integer(X), integer(N),
                   R is (N mod 2), R \= 0,
                   b('^',X,N,Z), Y = ~Z, ! .
  b('^',~X,N,Y) :- integer(X), integer(Y),
                   R is (N mod 2), R = 0,


                   b('^',X,N,Z), Y =  Z, ! .

  b('^',~X,Y, ~(X^Y)) :- ! .

  b('^',X,Y,X^Y).


     Fig. 7 - The simplification program.


  ?-/*          derivative of x^2 with respect to x with simplification
  */


                                     Y = x^2,

                                     d( Y, x, Z),

                                     s( Z, DY).


  Y = (x ^ 2),

  Z = ((2 * 1) * (x ^ (2 - 1))),

  DY = (2 * x)


     Fig. 8 - Derivative of x^2 with respect to x, using both
              programs DIFFSV.PRO and SIMPSV.PRO


  ?-/*          derivative of  y = (t^2-3)^4  with respect to t          */


                                    Y = (t^2-3)^4,

                                    d( Y, t, Z),

                                    s( Z, DY).


  Y = (((t ^ 2) - 3) ^ 4),

  Z = ((4 * ((2 * 1) * (t ^ (2 - 1)))) * (((t ^ 2) - 3) ^ (4 - 1))),

  DY = ((8 * t) * (((t ^ 2) - 3) ^ 3))


  ?-/*         derivative of  y = x^5 + 5*x^4 - 10*x^2 + 6   w.r.t. x     */


                                Y = x^5 + 5*x^4 - 10*x^2 + 6,

                                d( Y, x, Z),

                                s( Z, DY).


  Y = ((((x ^ 5) + (5 * (x ^ 4))) - (10 * (x ^ 2))) + 6),

  Z = ((((5 * 1) * (x ^ (5 - 1))) + (5 * ((4 * 1) * (x ^ (4 - 1))))) - (10 *
  ((2 *
   1) * (x ^ (2 - 1))))),

  DY = (((5 * (x ^ 4)) + (20 * (x ^ 3))) - (20 * x))


  ?-/*        derivative of  y = (2*x)^(1/2) + 2*x^(1/2)   w.r.t.  x      */


                             Y = (2*x)^(1/2) + 2*x^(1/2),

                             d( Y, x, Z),

                             s( Z, DY).


  Y = (((2 * x) ^ (1 / 2)) + (2 * (x ^ (1 / 2)))),

  Z = ((((1 / 2) * (2 * 1)) * ((2 * x) ^ ((1 / 2) - 1))) + (2 * (((1 / 2) * 1)
  * (
  x ^ ((1 / 2) - 1))))),

  DY = (((2 * x) ^ (-1 / 2)) + (x ^ (-1 / 2)))


     Fig. 9 - Examples of derivatives obtained with the programs.


  ?-/*     partial derivatives of  y = a*u + 1/v   w.r.t. u and v      */


                              Y = a*u + 1/v,

                              d( Y, u, V), s( V, DYDU),

                              d( Y, v, W), s( W, DYDV).


  Y = ((a * u) + (1 / v)),

  V = ((a * 1) + (~ 0  / (v ^ 2))),

  DYDU = a,

  W = ((a * 0) + (~ 1  / (v ^ 2))),

  DYDV = ~ (1 / (v ^ 2))


  ?-/*   1st, 2nd and 3rd derivatives of   y = x^3+3*x^2-8*x+2   w.r.t.  x
  */


                               Y = x^3+3*x^2-8*x+2,

                               d( Y, x, Z), s( Z, DY1),

                               d( DY1, x, V), s( V, DY2),

                               d( DY2, x, W), s( W, DY3).


  Y = ((((x ^ 3) + (3 * (x ^ 2))) - (8 * x)) + 2),

  Z = ((((3 * 1) * (x ^ (3 - 1))) + (3 * ((2 * 1) * (x ^ (2 - 1))))) - (8 *
  1)),
  DY1 = (((3 * (x ^ 2)) + (6 * x)) - 8),

  V = ((3 * ((2 * 1) * (x ^ (2 - 1)))) + (6 * 1)),

  DY2 = (6 + (6 * x)),

  W = (6 * 1),

  DY3 = 6


  ?-/*  y = u^3+4  ,  u = x^2+2*x  ;  derivative of y   w.r.t.  x      */


                                U = x^2+2*x,

                                Y = U^3+4,

                                d( Y, x, Z), s( Z, DY).


  U = ((x ^ 2) + (2 * x)),

  Y = ((((x ^ 2) + (2 * x)) ^ 3) + 4),

  Z = ((3 * (((2 * 1) * (x ^ (2 - 1))) + (2 * 1))) * (((x ^ 2) + (2 * x)) ^ (3
  - 1
  ))),

  DY = ((3 * (2 + (2 * x))) * (((x ^ 2) + (2 * x)) ^ 2))


  ?-/*    x = y * (1-y^2)^(1/2)  ;   derivative of y   w.r.t.  x          */


                             X = y*(1-y^2)^(1/2),

                             d( X, y, Z), s( Z, DX),

                             s( 1/DX, DY).


  X = (y * ((1 - (y ^ 2)) ^ (1 / 2))),

  Z = (((((1 / 2) * ~ ((2 * 1) * (y ^ (2 - 1))) ) * ((1 - (y ^ 2)) ^ ((1 / 2)
  - 1)
  )) * y) + (1 * ((1 - (y ^ 2)) ^ (1 / 2)))),

  DX = (((1 - (y ^ 2)) ^ (1 / 2)) - ((((2 * y) / 2) * ((1 - (y ^ 2)) ^ (-1 /
  2)))
  * y)),

  DY = (1 / (((1 - (y ^ 2)) ^ (1 / 2)) - ((((2 * y) / 2) * ((1 - (y ^ 2)) ^
  (-1 /
  2))) * y)))


     Fig. 10 - Partial derivatives, higher order derivatives, chain rule
               and derivative of implicit function.


  ?-/*   application of the Newton-Raphson method to calculate the real root
  of


                                 x^3 - a = 0


         ( x is the approximation to the solution at the k-th iteration and

           Xnew at the iteration k+1)

    */


                             F = x^3-a,

                             d( F, x, Z), s( Z, F1),

                             s( x-F/F1, Xnew).


  F = ((x ^ 3) - a),

  Z = ((3 * 1) * (x ^ (3 - 1))),

  F1 = (3 * (x ^ 2)),

  Xnew = (x - (((x ^ 3) - a) / (3 * (x ^ 2))))


     Fig. 11 - Example of use of the programs for the Newton-Raphson method.