home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power-Programmierung
/
CD1.mdf
/
prolog
/
tutorial
/
algebra
/
algebra.txt
next >
Wrap
Text File
|
1987-04-12
|
59KB
|
2,179 lines
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 Programming in Prolog 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 minimum program 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 instantiated 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 are 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 how 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 basic 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 heuristics 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 cut 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 invertibility,
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
politics of simplification (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 radical 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.
Conservative 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.
Liberal systems will perform some simplifications, but leave
others to the user. New left systems are variations of the old
radical systems, with some added features to be selected by the
user.
Catholic 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 not 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, Programming in Prolog (second
edition), 1984
2. H.Coelo, J.C.Cotta and L.M.Pereira, How to solve it with
Prolog (3rd edition), 1982
3. J.Moses, Algebraic simplification - a guide for the perplexed,
in Proceedings of the Second Symposium on Symbolic and Algebraic
Manipulation, 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.
