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 = xS2T 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 = xS2T 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 = xS2T 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 = xScT 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 = xS2T 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 xS2T 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 cS2T rather than the equivalent expression 1 / cS2T = 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 xS2T 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 xS2T 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 = xS2T , 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 = xS2T ? 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 xS2T/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 xS3T - 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) / vS2T d(uSvT)/dx = uSvT ( 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.