Dec 242017

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

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.

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.

December 24, 2017
Add comments