This set of Fortran 90 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
coefficients) and routines are provided for calculating them. The
source code is available from gitlab.com, and can be obtained either
by going to
https://gitlab.com/anthonyjs/rrf
and following the download link, or, more directly, by the command
git clone https://gitlab.com/anthonyjs/rrf.git
which will clone the code into a new directory rrf.
There is a program called rrfcalc, included in the distribution, that uses these routines and that can be used interactively to calculate Wigner 3j, 6j and 9j coefficients 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 variable of type(rrf). Such a variable contains a pointer to a type(node) variable, which is the head of a circular list in which the RRF is held in power-of-prime form. It also contains an integer factor which may hold an intermediate number that has not yet been factorized into a product of primes. Such intermediate numbers can get extremely large in complicated manipulations, so it is recommended to use the gfortran compiler, which provides 16-byte integers (up to 38 decimal digits). Other compilers, limited to 8-byte integers (18 decimal digits) will however be adequate for normal use. The program checks for integer overflow and should stop rather than give erroneous results.
All arithmetic and i/o manipulations are carried out by special routines provided in the rrf module, and the list elements are not available directly to a program that uses type(rrf) variables.
Further details of the program, in its original Fortran 77 form, may be found in the paper by A. J. Stone and C. P. Wood, Computer Physics Communications (1980) 21, 213. The current version of the program comprises a Fortran module rrf_module.F90 containing all the basic routines, a module wigner.F90 containing routines to calculate 3j, 6j and 9j coefficients, the calculator program rrfcalc.F90, and programs test3j.f90, test6j.f90 and test9j.f90 which carry out various tests of the wigner routines.
Index of subroutines
Initialization
Housekeeping
Arithmetic and Conversion
- MULT
- MULTI
- SETF
- MULTF
- POWER
- ADD
- ACCUM
- SUBTR
- EQUAL
- ROOT
- INT_TO_RRF
- PP_TO_RRF
- DISPLAY
- RRF_TO_4I
- CHAR4I
- WRITE4I
- RRF_TO_REAL
- RRF_TO_CHAR
- LIST
Angular Momentum Functions
Error handling
Notes
Initialization routine:
- SUBROUTINE INIT_RRF([quiet[,primes]])
- This must be called before any of the other routines. It sets up a table of primes and initializes the table of factorials. QUIET is an optional logical argument, which if true suppresses the initial program banner. PRIMES is an optional integer argument, which if present specifies the number of primes to be generated in the table of primes. Default 200000.
Housekeeping routines:
- SUBROUTINE CLEAR(K)
- Deallocate RRF K. 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.
- SUBROUTINE COPY(K1,K2)
- Set RRF K2 equal to K1. The Fortran assignment statement K2=K1 can also be used; it is overloaded to do the right thing.
- SUBROUTINE RENAME(K1,K2)
- Set RRF K2 equal to K1 and clear K1.
- SUBROUTINE EXCH(K1,K2)
- Exchange the RRFs K1 and K2. Use EXCH or RENAME rather than COPY where possible.
- SUBROUTINE UNIT(K)
- Set K to the RRF value of unity (not the same as K=1).
- FUNCTION ISZERO(K)
- Returns the value .TRUE. if K represents zero and .FALSE. otherwise.
- FUNCTION NONZERO(K)
- Returns the value .TRUE. if K represents a nonzero value and .FALSE. otherwise.
- SUBROUTINE CHSIGN(K)
- Change the sign of K (multiply by -1).
- SUBROUTINE CONJUGATE(K)
- Replace K by its complex conjugate.
Up to index
Arithmetic routines:
- SUBROUTINE MULT(K1,K2,M)
- Multiply K1 by K2**M. K1 and K2 must be different variables. Division is achieved by having M < 0.
- SUBROUTINE MULTI(K1,I,M)
- Multiply RRF K1 by I**M, I, M integers.
- SUBROUTINE SETF(N,K)
- Set RRF K to the value of factorial(N). The factorial table, which initially contains only factorial(0) and factorial(1), is extended if necessary.
- SUBROUTINE MULTF(K,N,M)
- Multiply RRF K by (factorial N)**M. The factorial table is extended if necessary.
- SUBROUTINE POWER(K1,M)
- Replace K1 by K1**M.
- SUBROUTINE 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.
- SUBROUTINE SUBTR(K1,K2)
- Subtract K2 from K1.
- FUNCTION 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.
- SUBROUTINE ROOT(K)
- Replace K by its square root.
Up to index
Conversion routines:
These routines convert RRFs to ordinary Fortran variables, or to values in an A1 character buffer, and vice versa, or print out a character representation of an RRF.
- SUBROUTINE INT_TO_RRF(I,K)
- Express integer I as RRF K.
- SUBROUTINE PP_TO_RRF(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 here.
- SUBROUTINE RRF_TO_4I(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.
- SUBROUTINE 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.
- SUBROUTINE RRF_TO_REAL(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.
- SUBROUTINE RRF_TO_CHAR(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.
Up to index
Angular momentum functions (module wigner):
- SUBROUTINE THREEJ(J1,J2,J3, M1,M2,M3, X, K)
- SUBROUTINE SIXJ(J1,J2,J3, L1,L2,L3, X, K)
- SUBROUTINE 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.
- SUBROUTINE THREE0(J1,J2,J3, K)
- A more efficient equivalent of SUBROUTINE THREEJ(J1,J2,J3, 0,0,0, 1, K)
- SUBROUTINE REGGE3(JP1,JP2,JP3, JM1,JM2,JM3, K)
- An alternative way to get a 3-j coefficient. 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 coefficient required.
Up to index
Error handling
There is a variable hard_errors in module rrfdata which is initially set true. With this setting most errors cause the program to stop after printing an error message. If it is set false, errors are soft, and cause the variable error in module rrflist to be set true after printing the error message. If another error occurs while error is true, the program will treat it as a hard error even if the variable hard_errors is false. If error is set false after an error, the program will continue, but results should be treated with suspicion.
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 coefficients, 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.
- Machine-dependent values are set in rrf_module.F90, and should be checked before compiling the program.
- 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 value for the parameter MAXFCT in rrf_module.F90.
- 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.