Root-rational-fraction 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 (root-rational-fraction) is referred to via a Fortran integer variable. A zero value denotes a zero RRF; a non-zero value is a pointer to a circular list in which the RRF is held in power-of-prime 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

Housekeeping

• START
• CLEAR
• COPY
• RENAME
• EXCH
• UNIT
• EQUAL

Arithmetic

• MULT
• MULTI
• SETF
• CLEARF
• MULTF
• POWER
• ADD
• ACCUM
• SUBTR
• MINUS
• ROOT

Input-Output

• FROMI
• FROMR
• FROMCH
• TO4I
• CHAR4I
• TOREAL
• TOCHAR
• LIST

Angular Momentum Functions

• THREEJ
• SIXJ
• NINEJ
• THREE0
• REGGE3

Error handling

Notes

Housekeeping routines:

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

Arithmetic routines:

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

Input-output routines:

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 4-integer 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

Angular momentum functions:

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 3-j symbol. JP1=J1+M1, JM1=J1-M1, etc., so all arguments are integers even if the quantum numbers are not. Since the 3-j routine needs J1+M1 etc. anyway, this is actually a more efficient way to define the symbol required.

Up to index

Error handling

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

Notes:

• 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 indeed--much 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.
• Machine-dependent 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.