The Silicon Graphics Teaching Laboratory has been replaced by the Chemi
cal Information Laboratory and this information is for historical interest only
Rootrationalfraction program package
This set of routines performs arithmetic on numbers whose squares are rational
fractions (possibly negative, so that the numbers may be pure imaginary).
In particular, the theory of angular momentum yields various numbers of
this form (the Wigner 3j, 6j and 9j symbols) and routines are provided for
calculating them. The source code is available as a compressed tarfile (about 25 kbytes).
There is a program, included in the distribution, that uses these routines
and that can be used interactively to calculate Wigner 3j, 6j and 9j symbols
and to perform simple arithmetic on them in the manner of a desk calculator.
Each RRF (rootrationalfraction) is referred to via a Fortran integer variable.
A zero value denotes a zero RRF; a nonzero value is a pointer to a circular
list in which the RRF is held in powerofprime form. The variable must
be set to zero initially and subsequently referred to only via the RRF subroutines,
until it has been cleared by means of the CLEAR subroutine, when it becomes available again for normal use.
Further details of the program may be found in the paper by A. J. Stone
and C. P. Wood, Computer Physics Communications (1980) 21, 213.
Index of subroutines
 CALL START(NL)
 This must be called before any of the other routines. It initialises the
list and sets up a table of primes. NL specifies the number of integer locations
in blank common which can be used by the package, and there may be a declaration
of the type
COMMON I(NL)
in the calling program. If it is omitted, 10,000 locations are reserved
by declarations in the subroutines.
 CALL CLEAR(K)
 Set K to zero and release its list space. This should always be used to
release the space occupied by temporary variables before exit from subroutines
 otherwise the list space will rapidly become full. There is no automatic
garbage collection.
 CALL COPY(K1,K2)
 Set the RRF K2 equal to K1. (Note that this is not the same as the Fortran
statement K2=K1, which must not be used.)
CALL RENAME(K1,K2) Set K2 equal to K1 and clear K1. CALL EXCH(K1,K2) Exchange the RRFs K1 and K2. Use EXCH or RENAME rather than COPY
where possible.
 CALL UNIT(K)
 Set K to the RRF value of unity (not the same as K=1).
 EQUAL(K1,K2)
 The logical function EQUAL has the value .TRUE. if K1 and K2 represent the
same number, .FALSE. otherwise. Note that this is not the same as K1 .EQ.
K2.
Up to index
 CALL MULT(K1,K2,M)
 Multiply K1 by K2**M. K1 and K2 must be different variables. Division is
achieved by having M<0.
 CALL MULTI(K1,I,M)
 Multiply RRF K1 by I**M, I, M integers.
 CALL SETF(N,K)
 Set RRF K to the value of factorial(N). The factorial table, which initially
contains only factorial(1), is extended if necessary.
 CALL CLEARF(N1,N2)
 Release the list space occupied by factorial(N1) to factorial(N2) inclusive.
This should only be necessary if the list space becomes congested.
 CALL MULTF(K,N,M)
 Multiply RRF K by (factorial N)**M. The factorial table is extended if necessary.
 CALL POWER(K1,M)
 Replace K1 by K1**M.
 >CALL ADD(K1,K2)
 Add K2 to K1, leaving K2 unchanged. K1 and K2 must be different variables.
An additional restriction, which is less of a limitation than it might seem,
is that the result of the addition must itself be expressible as an RRF.
Thus, for example, an attempt to add sqrt(3) to sqrt(5) will provoke an
error message. See Note 2.
 CALL ACCUM(K1,A1,K2)
 Sometimes the intermediate result in a long sequence of additions involves
prime factors which are too large to fit into an integer location, even
though the final sum can be factorized satisfactorily. This routine is provided
to meet this eventuality. It adds K2 to the number represented by A1 times
the rrf K1, and returns the result in the same form. If the factorisation
is successful, A1 will contain 1 on return.
 CALL SUBTR(K1,K2)
 Subtract K2 from K1.
 CALL MINUS(K)
 Change the sign of K.
 CALL ROOT(K)
 Replace K by its square root.
Up to index
These routines convert RRFs to ordinary Fortran variables, or to values
in an A1 character buffer, and vice versa.
 CALL FROMI(I,K)
 Express integer I as RRF K.
 CALL FROMR(A,K)
 Express the integer value held in the double precision variable A as RRF
K. This should work if the value is small enough to be held exactly, but
is not guaranteed. If the number in A cannot be completely factorized into
primes small enough to be held in integer locations, the remaining factor
is returned in A, which should therefore be tested on return. If the factorisation
is successful, A is set to 1.
 CALL FROMCH(M,N, NP, K)
 Read an RRF from the first N characters of the A1 character buffer M. The
detailed syntax of the RRF representation is given elsewhere.
 CALL TO4I(K, I1,I2,I3,I4)
 Express RRF K as (I1/I2)sqrt(I3/I4). I1 and I3 may be negative. I2 and I4
are always positive, unless integer overflow occurs, in which case zero
is returned in I2, and I1, I3 and I4 contain rubbish.
 CALL CHAR4I(K, M,M1,M2)
 Convert the rrf K to 4integer form as characters in the buffer M of length
M2, starting at position M1. M1 is updated to point at the next available
position in the buffer. If integer overflow occurs, 12 asterisks are returned.
 CALL TOREAL(K,A,I)
 Express RRF K as A*i**I, where A is double precision and the integer I has
the value 0 or 1.
 CALL TOCHAR(K, M, M1,MAX, NP)
 Express RRF K in A1 character form in the buffer M(MAX), starting at position
M1. The sign term is output as +, , +i or i; the exponents of the first
NP primes follow, and further terms appear as '. prime^exp'. NP may be zero.
A value of zero, +1 or 1 is output as 0, +1 or 1 in the sign position.
An error message is printed if there is not enough space in the buffer.
M1 is updated to point at the next available position in the buffer.
 CALLLIST(K)
 This is a debugging aid. It prints the actual list entries for RRF K on
unit 6. If it is called with a negative argument N, it prints the last
N list entries. Each entry occupies three integer locations; the first contains
a pointer to the next entry, the second contains a prime, and the third
is twice the exponent of that prime. The primes occur in increasing order,
and the last entry points back to the first, for which the 'prime' is always
1. 'Primes' larger than the square of the largest entry in the prime table
may in fact be composite, but this will not affect the working of the program.
Up to index
 CALL THREEJ(J1,J2,J3, M1,M2,M3, X, K)

 CALL SIXJ(J1,J2,J3, L1,L2,L3, X, K)
 CALL NINEJ(A,B,C, D,E,F, G,H,I, X, K)
 These routines set the appropriate vector coupling coefficient in the RRF
K. X is a Fortran integer which should have the value 1 or 2; each argument
(J1 etc.) is interpreted as J1 or J1/2 according as X=1 or 2.
 CALL THREE0(J1,J2,J3, K)
 A more efficient equivalent of CALL THREEJ(J1,J2,J3, 0,0,0, 1, K)
 CALL REGGE3(JP1,JP2,JP3, JM1,JM2,JM3, K)
 An alternative way to get a 3j symbol. JP1=J1+M1, JM1=J1M1, etc., so all
arguments are integers even if the quantum numbers are not. Since the 3j
routine needs J1+M1 etc. anyway, this is actually a more efficient way to
define the symbol required.
Up to index
The first integer location of blank common (called LE in the subroutines)
is used as an error flag. If LE is positive or zero (hard error option)
most errors cause the program to stop after printing an error message. If
LE is negative (soft error option) the action in the event of an error is
to set LE to the error number and continue, first printing the error message
if LE was 1 but not if it was 2. Note that this action implicitly sets
the hard error option for subsequent errors, so the value of LE should be
tested and reset if necessary if the soft error option is to be kept in
operation. LE is set to 0 (hard errors) when subroutine START is called.
Up to index
 The program should not be used unless exact results are essential, as is
it inevitably much slower than ordinary arithmetic. In particular, addition
may be very slow indeedmuch slower than multiplication.
 The embargo on adding numbers if the result cannot be expressed as an RRF
will seem a very severe restriction. In fact, the program has been used
for several fairly elaborate calculations, involving heavy use of 3j, 6j
and 9j symbols, without ever coming up against the restriction. I would
be interested to learn of a serious problem in which it is an important
restriction. It could be removed in principle, but would involve a complete
rewrite of the program.
 The program has been checked with the PFORT verifier for adherence to the
ANSI Fortran standard. The verifier flags a number of 'unsafe references'
which are in fact valid, but there are no violations of the standard.
 Machinedependent values are in COMMON /MCDATA/ and are set in the BLOCK
DATA subprogram.
 The angular momentum routines all use factorials. The largest available
value is factorial(200), which should suffice for most purposes. An error
message is printed if this limit is exceeded, in which case it is necessary
to recompile with a larger array FACT in the common block /TABLES/, and
with MAXFCT in the same common block set (in BLOCK DATA) to the length of
this array.
 The program has been fairly thoroughly tested and has been used successfully
on a number of problems, but some bugs doubtless remain. I should appreciate
being told about them.
