HomeWare 14

home *** CD-ROM | disk | FTP | other *** search

/ HomeWare 14 / HOMEWARE14.bin / educate / algtut.arj / PFSA.DOC next >

Wrap

Text File | 1985-10-03 | 26KB | 667 lines

TABLE OF CONTENTS I Introduction 2 II Usage 2 III Examples 4 A Basic Commands 4 B Procedures 4 C Canonical Transformation 5 D Binomial Expansion 6 E Inversion of Series 6 IV Commands 9 V Error Messages 13 VI Technical Details 14 -2- I. INTRODUCTION The PFSA program was created to do some rather large but simple algebraic computations. PFSA does not have the large variety of commands which more general purpose symbolic algebra programs (muMath and Macsyma) have but it is faster and can work with bigger expressions. PfSA ran 90 times faster than Macsyma on a very large Hamiltonian canonical transformation project, and it runs more than 200 times faster than muMath when doing example E. PFSA is written in Fortran and the arithmetic involving coefficients of terms is done with the Fortran arithmetic. Two versions are provided, one which uses only rational integer*4 arithmetic (PFSAI.EXE) and another (PFSAF.EXE) which uses real*8 arithmetic during calculation of coefficients. The EXE versions on this PC Blue disk need 300K to run. The program can be compiled with different values for parameters (in dimension statements) to get a smaller executable file or to get capability to handle bigger expressions. If you find a bug and want to report it to the author, put a minimal example "algin" file on a floppy and send it to: Don Stevens, 4 Wash. Sq. Village, NY NY 10012. If you want the floppy returned please include a stamped return envelope. Revised versions may be available in the future. The author of PFSA tried to make the usage of PFSA simple. PFSA reads a file "algin", processes the commands therein, and produces a file "algout". The commands in the algin file are mostly close enough to ordinary English so that the algebra to be performed is clear even to an inexperienced reader. PFSA commands allow one to do standard operations on polynomial expressions (addition, subtraction, multiplication, selection of portions of expressions, substitutions, differentiation). Binomial expansions may be generated with non-integer exponents. Numerical evaluation may be done by doing substitutions. Procedures may be defined. The commands are listed and described in Section IV. Section III is intended to be a tutorial for new users of PFSA. Section V is a listing of some error messages and their meaning. Section VI contains technical details of little interest to most users. -3- II. USAGE The user prepares a file "algin" which contains commands that specify the task to be done by PFSA. The first command should be the "order" command and the last command is the "end input" command. The commands of PFSA have been chosen so that the algin file will be approximately a description of the problem in English. A comment line may be inserted by beginning the line with ";" or "*", or a comment may be included on a command line after the command by beginning the comment with a ";". The user runs PFSA, which produces an output file "algout". The user can control the amount of output which PFSA produces during a run by specifying a parameter (named "np"). If np=0 only error messages are produced. After np=1 output from "print" and "fortran" commands is produced. After np=2 the commands are listed as they are encountered. After np=3 time and memory usage information is printed after each command is done. There are three types of objects in PFSA: variables, sums, and labels. A product of variables with a coefficient is called a term, and a sum of terms (perhaps only one) is called a sum. Locations may be designated by using the "label" command, and can be referred to by the "goto" command. The names of variables, sums, and labels should begin with a letter and should not contain any of the following characters: ; + - * / = . Upper case and lower case are distinguished. The default maximum length of names is 15 characters. It is recommended that no names begin with "flu", "lab", or "pro", see the discussion at the start of Section IV. The first command in an algin file should be "order", for example: order var1 var2 var3. This command can be used only once per run. The effects of the above command are: 1. In all terms printed, var1 would precede var2, which would precede var3. 2. All sums printed will be ordered first on the exponent of var1 (smaller exponents precede larger exponents), then on the exponent of var2, etc. 3. If subsequently the command maxe = 2 were given, thereafter whenever a term containing var1 raised to a power greater than 2 is generated during a computation, this term would be discarded. Previously defined sums are not affected. Only var1 is affected. The commands "compute" and "define" have some overlap in their use but the following distinctions should be noted. The "compute" command is always a one line command operating on sums only, while the "define" command may be many lines and may operate on variables also. Only the "compute" command may be used to raise a sum to a non-integer power (binomial expansion) or to divide one sum by another sum. -4- III. EXAMPLES The following examples give the flavor and style of usage of PFSA. Section IV has descriptions of the full set of commands for PFSA. This section is intended to introduce the use of PFSA, not specify it. Refer to Section IV for more exact definitions. The examples are algin files which may be run using PFSAI.EXE. The algin files can be run with PFSAF.EXE if all the fractional numbers are changed to decimal numbers (Fortran real numbers). A. Basic Commands order x y define sum1 = 3/4 x**2 define sum2 1 4 sum1 7 y 1/3 z define sum3 = sum1 * sum2 print sum3 the derivative of y wrt x is dydx differentiate sum4 = d sum3/d x print sum4 substitute 8 for x in sum4 print sum4 fortran sum4 define one = 1 define six = 6 define i = one label lab2 compute i = i + one print i if i.eq.six end goto lab2 end input B. Procedures order x y z * This example contains a procedure which computes factorials. * The lines from "procedure factorial" to "return" define the * procedure. There is no provision for local variables or * sums in procedures. The sum "facarg" is used to * transmit the argument, and the sum "facans" is used to * transmit the answer. -5- procedure factorial np=0 ;After this command only error messages are printed. define one=1 define facans=one if facarg.eq.one goto facend define temp=facarg if temp.lt.one bomb label fact1 compute facans=facans*temp compute temp=temp-one if temp.ne.one goto fact1 label facend np=2 ;after this command the usual printing resumes return define facarg=5 call factorial print facans end input C. Canonical Transformation Here we will find a transformation of coordinates (p,q) to coordinates (pb,qb) which transforms an expression h=p**2 + q**2 +8*e*p*q to a better form in which there is no term having e to the first power. Here the transformation is specified by the generating function f=pb*q + e*pb**2 - e*q**2; the transformation is defined by p=df/dq, qb=df/dpb. One must start from these two equations and work out the transformation explicitly. order e p q pb qb * e is a small parameter (epsilon) define h p^2 q^2 8 e p q define f pb q e pb^2 -e q^2 differentiate dfdq = d f/ d q differentiate dfdpb = d f/ d pb print dfdq print dfdpb * There is no "solve" command in PFSA, one can use various * commands such as "compute" and "extract" to aid in big * computations. Here we see by a simple calculation: * q=qb - 2 e pb and p=pb - 2 e qb + 4 e^2 pb . define qequals qb -2 e pb -6- define pequals pb -2 e qb 4 e^2 pb substitute qequals for q in h substitute pequals for p in h print h * We see that the e term was eliminated. More transformations * can be done to eliminate other unwanted terms. end input D. Binomial Expansion order x define mhalf = -.5 ;this example should be run with PFSAF.EXE maxe=30 nbinomial=30 define aa 1 -x^2 compute sb = aa^mhalf substitute t for x in sb the derivative of y wrt t is sb define arcsin = y taylor order = 31 taylor expand arcsin in t substitute 0 for y in arcsin substitute x for t in arcsin print arcsin substitute .5 for x in arcsin define pi = 6.*arcsin print pi ;accurate to 10 decimal places end input E. Inversion of Series * if y = x + b1 x^2 + b2 x^3 ... is given, this routine finds * x = y + c1 y^2 + c2 y^3 ... in general terms. * This example was contributed by Kaya Imre. PFSAI does this * example in 40 seconds on the XT (Microsoft compiler). PFSAI * does this example in 29 seconds on a Cromemco (8 Mhz 68000). * The equivalent code in muMath runs in 11160 seconds on the XT. * The muMath code (also by Imre) is given after the algin file. -7- order x y np=1 def gen x b1 x^2 b2 x^3 b3 x^4 b4 x^5 def rev y c1 y^2 c2 y^3 c3 y^4 c4 y^5 def one = 1 def two = 2 def zero= 0 def yy = y def dup = rev / yy com dup = dup - one maxe = 6 sub gen for y in rev def cnt = 0 def six = 6 def ex = x def rrev= y def yp = y lab loop def tem = rev sub q for x in tem sub 0 for q in tem com rev = rev - tem def rev = rev / ex com cnt = cnt + one if cnt.le.two goto loop def yp = yp * yy def dup = dup / yy def cc = dup sub 0 for y in cc com dup = dup -cc sub zero for cc in tem def tt = -tem sub tt for cc in rev def rr = -tem * yp com rrev= rrev + rr np = 2 pri rr np = 0 if cnt.lt.six goto loop end input -8- FUNCTION INVRT(F,G,N), GG:G, GH:G, H:(EXPAND(EVSUB(G,Y,F))/X-1)/X, GH:(GH/Y-1)/Y, LOOP WHEN N EQ 0 EXIT, COEF:EVSUB(GH,Y,0), AA:COEF-EVSUB(H,X,0), H:EVSUB(H,COEF,AA)/X, GH:(GH-COEF)/Y, GG:EVSUB(GG,COEF,AA), N:N-1, ENDLOOP, GG ENDFUN$ RDS()$ -9- IV. COMMANDS The input for PFSA is fairly free form. All commands should be in lower case. Commands may be abreviated. The first 3 letters of the command name are sufficient to identify the command, but some commands require additional words in the command line. PFSA reads the algin file in segments (delimited by the "flush" command), and looks for certain commands during this buffer filling. If a line has "flu" or "lab" or "pro" as the first non blank characters it will be interpreted as a command by the preprocessor (generating an error if this were in fact part of a "define" or an "order" block). A line that starts with an asterisk or semi-colon in column 1 will be considered by PFSA to be a blank line. PFSA will ignore blank lines, except those blank lines used to terminate "define" or "order" blocks. In the following descriptions v1, v2, and v3 represent variables, and s1, s2 and s3 represent sums. The listed constructions are the only valid constructions, for example there is no valid "compute" construction involving a variable. Vors1 and vors2 are used to indicate that a variable or a sum may be used. call label1 compute s1=s2 compute s1=s2+s3 compute s1=s2-s3 compute s1=s2*s3 compute s1=s2/s3 The above division is done only if the result is expressible as a sum of terms. So if s3 contained only a single term the division would be done. If s3 has more than one term then if possible s1 will be found such that s1*s3 is equal to s2 up to order maxe in the first variable. This is not implemented for maxe greater than 2. So if for instance maxe=2, s2=1, and s3=1+var1, then s1 will be 1 - var1 + var1**2. compute s1=s2^s3 S1 will be the binomial expansion of s2 to order nb(nbinomial). S2 must have its first term equal to 1, so variable#1 cannot appear in s2 with a negtive exponent. If we write s2=1+z and if s3=a, then s1 will be 1+a*z+a(a-1)z^2/2+..+a(a-1)..(a-nb+1)z^nb/nb!. count s1 = terms in s2 This command causes the number of terms in s2 to be counted and s1 is set equal to this number. define s1=coefficient vors1 vors2 .... In the above form of the define command, "coefficient" represents a rational number (or integer or nothing) and vors1 and vors2 are variables or sums (perhaps with exponent). The above form may not be used if s1 has more than one term, in which case the form below should be used. -10- define s1 coefficient vors1 vors2 .... coefficient vors1 vors2 .... .. (blank line to terminate sequence of terms) delete s1 s2 ... Variables can not be deleted or used as a sum. derivative specification is done with the "the derivative..." command, see below. differentiate s1=d s2 / d v1 differentiate s1=d s2 / d v1 order=number The order of the derivative to be taken can also be specified by a sum s3 which has been set equal to a positive integer, shown below. differentiate s1=d s2 / d v1 order=s3 dump end end input PFSA reads input lines in algin until it reaches "end input". extract s1 = term s2 of s3 In the above form s2 should be a sum which is currently equal to a positive integer. In the form below i and j are integers, and there are no spaces between these integers and the "-". extract s1 = terms i-j of s2 flush The input file algin is stored in a buffer in the running program. If algin is too long to fit completely in the buffer, it should be segmented using the flush command. When this command is executed the previous buffer contents are discarded and the next segment is read in. The call and goto commands are valid only for labels in the current segment. However, all information about variables and sums is kept. fortran s1 S1 will be printed in a form suitable for use in a fortran program. goto label1 if s1.eq.s2 command if s1.ne.s2 command The above constructions test whether s1 and s2 are identically equal; if the relation is true, the command is executed. The command should be a one line command. One may use the other fortran relational operators .lt. .le. .gt. and .ge. but in these cases only the first coefficients of s1 and s2 are compared. -11- integration such as z = integral y dx may be done by: (1) the deriv of q wrt x is y (2) taylor expand q in x (3) define z = q (4) subst 0 for q in z label label1 The next line has label=label1. maxe = i After this command, any terms newly arising which have variable#1 to a higher power than i will be discarded. nbinomial=i Binomial expansions will be done to order i. ndump=0 No dump ndump=1 A dump will occur after fatal errors. np=0 only error messages are printed np=1 the above plus results of "print" and "fortran" np=2 the above plus all commands are printed np=3 the above plus time and memory usage are printed The "first" variable must be specified at the start of an algin file by the order command. Other variables will be ordered according to the order in which they are encountered as the algin file is processed, and they need not be listed in the order command. order v1 v2... The above form of the order command can be used if the list v1 v2... fits on one line. If there are more variables to be ordered than will fit on one line, use the following form: order v1 v2 .... .. ... (blank line to terminate the list of variables) print s1 procedure label1 All the lines between the above command and the "return" command constitute the procedure with the name label1. The procedure can be invoked in the segment (see the flush command) where it occurs by the call command. Procedures should not call other procedures. All names are global and available to to the procedure. return -12- substitute number for v1 in s1 substitute v1 for v2 in s1 substitute s1 for v1 in s2 In the above command, if v1 occurs to a negative power somewhere in s2 and s1 has more than one term, an error message is printed and the substitution is not done. substitute v1 for s2 in s3 substitute s1 for s2 in s3 A limited form of pattern matching is provided by the above command. It is possible to replace any pattern which is constituted of a single term. If a substitute command is given in which the pattern is a sum s2 with more than one term, the first term (of sum s2) is examined and used in searching for matches and the remainder of the terms (in sum s2) are ignored. taylor expand s1 in v1 v2 .. The order of the taylor expansion wanted should have been previously specified by the following command: taylor order = i where i is a positive integer or a sum. The expansion is made about v1=v2=..=0. the derivative of v1 wrt v2 is vors1 The above command is used to tell PFSA the derivatives of variables. Subsequently PFSA will use this information when differentiating expressions containing these variables. If the derivative vors1 is a sum, a "chain rule" will stop at the sum because the substitution of the sum into the place where it will eventually go is not done until all chaining is finished. -13- V. ERROR MESSAGES bad cmnd bad command call an undefined label used in "call" command comput1 bad syntax in "differentiate" command comput2 invalid numerical order of differentiation specified comput3 no "/" found in "differentiate" command comput4 right hand side not a sum in "compute s1=s2" command comput14 order specified incorrectly in "differentiate" command comput15 the order for differentiation should be positive comput16 the order for differentiation should be an integer count1 bad syntax in "count" command count2 final argument in "count" command should be a sum deriv2 the thing to be differentiated should be a sum deriv4 in doing differentiation a chain rule sequence longer than 5 encountered expnd1 the item to be raised to a power is not a sum expnd2-3 the first term of the item must be 1 extrct1 the final argument in "extract" command should be a sum extrct2 a non-reasonable term was requested in "extract" command extrct5 bad syntax in "extract" command extrct6 terms can be extracted only from sums extrct7 j is less than i in "extract s1=terms i-j of s2 gcdout1 an arithmetic mistake has occurred, a denominator is negative getexp2 input line too long getexp4 bad syntax in input line getexp6 input line too long getexp8 bad syntax in input line getexp9 input line too long iftest1-5 bad syntax in "if" command max jp a sum is being raised to a power greater than 30 in a substitution process. Change PFSA or algin to fix this. monitr2-6 bad command monitr7 a label was searched for and not found in current segment monitr8-10 bad command nxt bad syntax in input line nxtch bad syntax in input line nxtnb bad syntax in input line print2 exponents greater than 99 are not handled by print routine print3 a term is requiring too many characters in being printed -14- putfirst the "order" command should be the first command and should occur only once. readin2-10 bad syntax in substitute command readin12 one should not substitute a number for a sum. Use a term as intermediate substitution. return? bad command taylor1 bad value for order in "taylor" command taylor2-5 bad syntax in "taylor" command taylor6 only sums can be taylor expanded by PFSA term1 input line too long thderv1-2 bad syntax in "the derivative" command useorder the first command in the algin file should be "order" -15- VI. TECHNICAL DETAILS The fortran source file PFSAI.FOR is included so that users can change the size or otherwise modify PFSAI.EXE as wanted. The compiled PFSAI.OBJ should be linked with TSUBS.OBJ (obtained by assembling TSUBS.ASM) to produce PFSAI.EXE. The main routine calls subroutine "monitr" which reads in and executes segments of the "algin" file. The segments are delimited by the "flush" command, the last segment is delimited by the "end input" command. Each segment is kept in memory while it is being executed, and so target locations for "goto" and "call procedure" commands should be in the same segment as the command. (For machines which allow large memory usage this can be avoided by using a single segment.) The routine "monitr" examines input lines, decides which command is intended, and calls the appropriate subroutine to execute the command. The sums and variables are referenced through a single set of lists: list "pv" has the print form of the name, list "ittb" is zero for variables and has the beginning term for sums, list "itte" has the number of the final term for sums, list "ios" has the ordinal position in memory for the sums. So for example if ittb(1)=30, itte(1)=30, ittb(2)=3, itte(2)=29, ittb(3)=1, itte(3)=2, the first two terms in the list of terms constitute sum number 3 the next 27 terms constitute sum number 2 and term 30 constitutes sum number 1, and ios(1)=3, ios(2)=2, ios(3)=1. An inverse of ios is kept in list "iosi". A given sum name will keep the same sum number (unless it is deleted and reintroduced) but it may be moved to different positions in the list of terms. The coefficients of terms are kept in the list "tm" and pointers to the list of factors and powers are kept in lists "itb" and "ite". After every command all the sums are put into a canonical sorted order, which is the same order as is displayed when the sum is printed. All sorting is done from start to end and if an out of sequence term is found the correct location is found by a binary search and then the pointers are changed for the affected terms. When terms are created during a multiplication or substitution command the above ordering procedure is done. At the end of generating the terms a packing operation is done (by a routine named "pack") which arranges all terms in memory so the pointers are ordered. Derivatives are handled by keeping a list of derivatives of variables and searching this list (up to 5 deep) for all chains which end with the given independent variable. Of course, the program knows d x**n/ dx, but other derivatives (if any) must have been previously specified by a "the derivative of u wrt x is smthng" command.