Dec 242017
Symbolic math and equation solver (like MuMath).
File PFSA.ZIP from The Programmer’s Corner in
Category Science and Education
Symbolic math and equation solver (like MuMath).
File Name File Size Zip Size Zip Type
-README 2852 1129 deflated
PFSA.DOC 26359 8532 deflated
PFSAF.EXE 95836 42489 deflated
PFSAI.EXE 96750 41744 deflated
PFSAI.FOR 70089 15770 deflated
PICOM.FOR 1023 417 deflated
TSUBS.ASM 663 212 deflated

Download File PFSA.ZIP Here

Contents of the PFSA.DOC file


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



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.


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
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.


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
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.


procedure factorial
np=0 ;After this command only error messages are printed.
define one=1
define facans=one
if goto facend
define temp=facarg
if bomb
label fact1
compute facans=facans*temp
compute temp=temp-one
if goto fact1
label facend
np=2 ;after this command the usual printing resumes

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
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
-2 e pb


define pequals
-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
define aa

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.

order x y
def gen
b1 x^2
b2 x^3
b3 x^4
b4 x^5

def rev
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 goto loop
end input




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

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.

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


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

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 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.

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.

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
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:
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.

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.


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

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"


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.

 December 24, 2017  Add comments

Leave a Reply