% Latex file % New efficient programs to calculate general recoupling coefficients % Part I: Generation of a summation formula % V. Fack, S.N. Pitre, J. Van der Jeugt % Comp. Phys. Commun. 83 (1994), 275-292. % (pictures have been commented out here; see the original for pictures) % \documentclass[12pt]{article} %\usepackage{a4} \newcommand{\sixj}[6]{\left\{ \begin{array}{ccc} #1 & #2 & #3 \\ #4 & #5 & #6 \end{array} \right\}} \newcommand{\ninej}[9]{\left\{ \begin{array}{ccc} #1 & #2 & #3 \\ #4 & #5 & #6 \\ #7 & #8 & #9 \end{array} \right\}} \newcommand{\recoup}[2]{\langle \: #1 \: | \: #2 \: \rangle} \newcommand{\bra}[1]{\langle \: #1} \newcommand{\ket}[1]{| \: #1 \: \rangle} \begin{document} % % TITLE % ----- % \title{New efficient programs to calculate general recoupling coefficients \\ Part I~: Generation of a summation formula} \author{V.~Fack, S.N.~Pitre, J.~Van der Jeugt\thanks{Senior Research Associate N.F.W.O.~(National Fund for Scientific Research, Belgium)} \\ \small Department of Applied Mathematics and Computer Science, \\ \small Universiteit Gent, Krijgslaan 281--S9, B--9000 Gent, Belgium \\ \small E-mail~: Veerle.Fack@rug.ac.be} \date{Published in~: Comput.\ Phys.\ Commun.\ {\bf 83} (1994) 275--292} \maketitle \vfill \begin{abstract} A program to generate a summation formula in terms of 6-$j$ coefficients for a general angular momentum recoupling coefficient is described. This algorithm makes use of binary tree transformations as introduced by Burke~\cite{Bur} in the program {\tt NJSYM}. Due attention is paid to finding an optimal summation formula with a minimal number of summation variables, thereby improving the results of {\tt NJSYM}. The results obtained here are at least as good as (and often better than) the results of the alternative approach {\tt NJGRAF}, introduced by Bar-Shalom and Klapisch~\cite{Bar}, using more advanced graphical methods. \end{abstract} \vfill % % PROGRAM SUMMARY % --------------- % \addtolength{\baselineskip}{0.2\baselineskip} \pagebreak {\large \bf PROGRAM SUMMARY} \begin{description} \item[Title of program~:] {\tt NJFORMULA} \item[Catalogue number~:] \item[Program obtainable from~:] CPC Program Library, Queen's University of Belfast, \\ N.~Ireland (see application form in this issue) \item[Computers used (Operating systems)~:] 386/486-based PCs (MS-DOS, Linux~\cite{Linux}), Sun Sparc (Unix), T800 transputer system (MS-DOS). \item[Programming language used~:] C \\ (Compilers~: Turbo C++~\cite{TC}, GNU CC~\cite{GnuC}, SPARCompiler C, 3L Parallel C~\cite{ParC}). \item[No.\ of lines in source program (module + 2 example main programs)~:] 950 \item[Keywords~:] atomic structure, nuclear structure, scattering, general recoupling coefficient, angular momentum, Racah coefficient, 3$n$-$j$ coefficient, coupling tree, binary trees, recursive search process. \item[Nature of physical problem~:] A general recoupling coefficient for an arbitrary number of (integer or half-integer) angular momenta is expressed as a multiple sum over products of 6-$j$ coefficients, including phase factors and square root factors. This summation formula can then be evaluated for given values of the angular momenta (for this purpose we will provide a program {\tt NJSUMMATION}~\cite{FPV}). \item[Method of solution~:] A summation formula for a general recoupling coefficient is obtained by determining a path between the binary coupling trees corresponding to the state vectors, using so-called flop operations. Such a path consists of a number of `successful flop sequences', each of which is found by a recursive search process. To obtain an optimal summation formula, i.e.\ with a minimal number of summation variables, the total number of flop operations must be minimal. Heuristic techniques are used to keep the length of the path as short as possible. \item[Typical running time~:] For typical examples, obtaining single, double or triple summation formulae, the running time is of the order of 0.01~second on the Sun Sparc and on a 486-based PC. For more complicated examples, such as generating a 7-fold and a 9-fold summation formula, the program takes less than 0.5~second. On the whole we find our running times to be comparable with those of {\tt NJGRAF}~\cite{Bar}. \end{description} % % INTRODUCTION % ------------- % \pagebreak {\large \bf LONG WRITE-UP} \bigskip \section{Introduction} % In this paper we tackle an old problem by means of a new method. The problem is the calculation of a general recoupling coefficient for an arbitrary number of integer or half integer angular momenta. Such computations appear typically in atomic and nuclear structure or scattering calculations. A famous program of Burke~\cite{Bur}, {\tt NJSYM}, is dealing with the same problem. In Burke's approach, the problem is divided in essentially two parts~: first, the general recoupling coefficient is expressed as a (multiple) sum over a product of 6-$j$ coefficients (including phase factors and square root factors); secondly, this expression is evaluated for given data of the angular momenta. Our approaches to both parts are strikingly different from the classical ones, particularly because we make extensive use of recursive programming. With this view, it is appropriate to treat the two parts of the problem as two separate topics. In the present paper we shall be dealing with the first part of the problem, namely the construction of a summation formula for the evaluation of a general recoupling coefficient. In a second paper~\cite{FPV}, we shall present the new approaches, coming from recursive programming and parallelisation, to the remaining part of the problem, i.e.\ that of the actual evaluation of such a summation formula. It follows already from the work of Burke~\cite{Bur} that a general recoupling coefficient can be expressed as a sum over products of 6-$j$ coefficients. In Biedenharn and Louck~\cite{Bie} a special topic (Topic~12) is dealing with this problem and there this fact is formulated as a {\em fundamental theorem}~: the transformation coefficient between each pair of binary coupling schemes is expressible in terms of sums over products of Racah coefficients (ignoring the presence of dimension and phase factors). A binary coupling scheme is simply expressing the order in which $n$ angular momenta are being coupled, e.g.\ ($n=4$)~: \begin{eqnarray} \lefteqn{\ket{((j_1,j_2)j_5,(j_3,j_4)j_6)j_7 \; m}=} & & \nonumber\\ &&\sum_{m_1,\ldots,m_6} C^{j_1 j_2 j_5}_{m_1 m_2 m_5} C^{j_3 j_4 j_6}_{m_3 m_4 m_6} C^{j_5 j_6 j_7}_{m_5 m_6 m} \ket{j_1 m_1} \otimes \ket{j_2 m_2} \otimes \ket{j_3 m_3} \otimes \ket{j_4 m_4} . \end{eqnarray} Herein, $C^{j j' j''}_{m m' m''}$ is a vector-coupling (Wigner or Clebsch-Gordan) coefficient~\cite{Edm,Bie}. The relation between such a coupling scheme and a binary tree is clear~\cite{Bur,Bie}. A general recoupling coefficient is then a transformation coefficient between two such coupling schemes, e.g.~: \begin{equation} \recoup{((j_1,j_2)j_5,(j_3,j_4)j_6)j_7}{(j_1,((j_2,j_3)j_8,j_4)j_9)j_7}. \label{grc} \end{equation} Here, the $m$-dependence is dropped since such coefficients are independent of $m$~\cite{Edm}. To determine a formula for the calculation of an arbitrary general recoupling coefficient (of $n$ angular momenta), Burke~\cite{Bur} determined an algorithm to express this coefficient as a sum over Racah or 6-$j$ coefficients. For example, for the recoupling coefficient~(\ref{grc}) Burke's program {\tt NJSYM} would yield the expression~: \begin{eqnarray} \lefteqn{(-1)^{j_1+j_2-j_5}[(2j_5+1)(2j_6+1)(2j_8+1)(2j_9+1)]^{1/2} } \nonumber\\ &&\times\sum_k (-1)^{j_1+j_8+k}(2k+1) \sixj{j_4}{j_3}{j_6}{j_5}{j_7}{k} \sixj{j_1}{j_2}{j_5}{j_3}{k}{j_8} \sixj{j_1}{j_8}{k}{j_4}{j_7}{j_9}. \label{grcx1} \end{eqnarray} This is a single sum (one summation variable $k$). However, this is by no means the simplest expression for~(\ref{grc}). In fact, using for example the Biedenharn-Elliott sum rule, formula~(\ref{grcx1}) becomes \begin{eqnarray} \lefteqn{(-1)^{j_1+2j_2+j_3+j_4+j_6+j_7+j_9}}\nonumber\\ &&\times[(2j_5+1)(2j_6+1)(2j_8+1)(2j_9+1)]^{1/2} \sixj{j_1}{j_2}{j_5}{j_6}{j_7}{j_9} \sixj{j_2}{j_3}{j_8}{j_4}{j_9}{j_6}. \label{grcx2} \end{eqnarray} From the computational point of view, it is clearly better to use formula (\ref{grcx2}) instead of (\ref{grcx1}) in order to calculate the recoupling coefficient (\ref{grc}), because (\ref{grcx2}) contains no sum and (\ref{grcx1}) is a single sum. As we shall see in the following sections, (\ref{grcx1}) and (\ref{grcx2}) are both related to a path between the binary trees corresponding to the coupling schemes in the bra- and ket-vector of (\ref{grc}); (\ref{grcx1}) corresponds to a path with 3 steps (i.e.\ 3 {\em flops}), whereas (\ref{grcx2}) corresponds to a path with 2 steps. In order to find an optimal expression for a general recoupling coefficient, a shortest path (one with the smallest number of flops) should be found. Herein, a flop is a certain operation on a binary tree; this will be discussed in more detail in Section~2. Thus, the problem of finding an optimal expression for the calculation of a general recoupling coefficient is reduced to a graph theoretical problem~: how to find the shortest path between two (labelled) binary trees using flop-operations. This turns out to be a rather difficult and intriguing problem in graph theory. In Burke's approach~\cite{Bur}, his method is equivalent to finding a certain path between the two binary trees. However, it is usually not a shortest path, as can be seen from the above example. In fact, the path found by means of {\tt NJSYM} is generally rather long, thus yielding expressions which are far from optimal. In order to improve {\tt NJSYM}, Bar-Shalom and Klapisch~\cite{Bar} developed a new program {\tt NJGRAF}. Their way of producing a ``best formula'' (where the number and values of summation variables are minimal) was to use graphical methods as developed by Yutsis, Levinson and Vanagas~\cite{Yut} and explained in many books~\cite{Bie,Bri,Jud,Elb}. This graphical method is very powerful, and for many transformation coefficients it will yield the most optimal expression in terms of 6-$j$ coefficients. However, the algorithms to reduce such graphs are not easy to implement, and for a complicated graph many tests and searches need to be performed in order to find an appropriate reduction, leading to an expression in terms of sums over products of 6-$j$'s. Nevertheless, it should be said that the program using such graphical representation, {\tt NJGRAF}, performs quite well. The purpose of the present paper is to point out that the method of graphical representations is not needed, and that the much simpler method of binary tree transformations can succesfully be used. One advantage is that the algorithm is much easier and rather straightforward~: there is in fact only one kind of search. The data structure is also easier~: this makes the program very transparent. In order to implement such a search it is natural to use a recursive algorithm; this will be analysed and explained in the following sections. It should be mentioned that the algorithm developed here is heuristic in the sense that it produces a short path between the binary trees, however it does not necessarily produce the shortest possible path. Roughly speaking, this is because it uses consecutive local minima and does not try to determine a global minimum. To determine such a global minimum (and thus a shortest path) would be too time and space consuming and one has to find a compromise between the optimality of the generated formula and the search time and space. Nevertheless, our method is so good that for all practical cases it does produce a shortest path and thus the most optimal expression. For general recoupling coefficients that are more complicated our algorithm yields a very short path (although it is not always possible to prove that it is indeed a shortest path), and in fact for all test examples in the literature and for the examples we give ourselves, our algorithm produces an expression which is often better than (and sometimes as good as) the one obtained by means of {\tt NJGRAF}. By a better expression, we mean a (multiple) sum over products of 6-$j$ coefficients with fewer summation variables and fewer 6-$j$'s. These advantages (simpler algorithm, easy data structure, better result) convinced us that our approach is worth following and in this paper we present our study in further detail. In a second paper devoted to this topic~\cite{FPV}, we shall present a number of novelties concerning the actual calculation of a multiple sum over products of 6-$j$'s; this will lead to a new program {\tt NJSUMMATION} as a replacement for the classical subroutine {\tt GENSUM}~\cite{Bur,Bar}. The structure of the paper is as follows. In section~\ref{sec-meth} we describe how a recoupling coefficient can be represented by a pair of binary coupling trees and how a corresponding summation formula can be obtained by a sequence of basic tree operations (so-called flops and exchanges). Section~\ref{sec-algo} describes an algorithm for constructing such a path, performing the transformation of one tree into the other; here it is also explained which ideas are used to keep the path as short as possible. In section~\ref{sec-impl} we give some implementation details about the datastructures used to represent a coupling tree and a summation formula, about the input of general recoupling coefficients, and about the transformation algorithm. Finally in section~\ref{sec-res} some results obtained by our new approach are discussed and compared with those obtained by the program {\tt NJGRAF}~\cite{Bar}. % % METHOD OF SOLUTION % ------------------ % %\pagebreak \section{Method of solution} \label{sec-meth} % In a general recoupling coefficient, such as \begin{equation} \recoup{((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9}{ (((j_1,j_4)j_{10},(j_2,j_3)j_{11})j_{12},j_5)j_9}, \label{eq-exa} \end{equation} each binary coupling scheme or state vector can be represented by a so-called {\em coupling tree\/}, which is a binary tree, as is shown in figure~\ref{fig-tree}. \begin{figure}[htp] \begin{center} %%%%%%\input tree.pic \begin{minipage}{12cm} \caption{Coupling trees for the general recoupling coefficient $ \recoup{((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9}{ (((j_1,j_4)j_{10},(j_2,j_3)j_{11})j_{12},j_5)j_9} $.} \label{fig-tree} \end{minipage} \end{center} \end{figure} Each $j_i$ coefficient corresponds to a {\em node\/} of the coupling tree. Some of the $j_i$ coefficients are the coupling of two other $j_i$ coefficients, e.g.\ $j_6, \ldots, j_{12}$. The corresponding node in the coupling tree has nodes down to the left and to the right, called the left and right {\em children\/} of the node. We call this a {\em coupled node}, to distinguish it from a {\em leaf node}, i.e.\ a node corresponding to a $j_i$ coefficient at the lowest level of the tree, e.g.\ $j_1, \ldots, j_5$. Leaf nodes have no left or right children. The {\em root node\/} is the top node of a coupling tree. Note that every coupled node can be considered as the root node of a {\em subtree\/} of the coupling tree. In a recoupling coefficient the two coupling trees will have a common root node, e.g.\ $j_9$, otherwise the recoupling coefficient is trivially zero. For all practical purposes this means that the root node, and in general every coupled node, is uniquely identified by the set of its descendent leaf nodes, an idea we will use in the algorithm described in section~\ref{sec-algo}. As described by Burke~\cite{Bur}, a general recoupling coefficient can be expressed as a (multiple) sum over products of 6-$j$ coefficients, including phase factors and square root factors. Such a summation formula can be obtained by a transformation process on the coupling trees corresponding to the recoupling coefficient. For that purpose a set of basic operations ({\em exchanges} and {\em flops}) on a coupling tree are defined, which correspond to simple recoupling coefficients, and thus to simple summation expressions. An {\em exchange\/} is a transformation of a state vector of the form $\ket{(a,b)c}$ to $\ket{(b,a)c}$, for which the coupling trees are shown in figure~\ref{fig-exch}, and corresponds to the following recoupling coefficient~\cite{Edm}~: \begin{eqnarray} \recoup{(a,b)c}{(b,a)c} & = & (-1)^{\pm(a+b-c)}. \end{eqnarray} \begin{figure}[htp] \begin{center} %%%%%\input exch.pic \begin{minipage}{12cm} \caption{Coupling trees for an exchange~: $ \recoup{(a,b)c}{(b,a)c} $.} \label{fig-exch} \end{minipage} \end{center} \end{figure} A {\em flop\/} is a transformation of a state vector of the form $\ket{((a,b)d,c)f}$ to $\ket{(a,(b,c)e)f}$ or vice versa, for which the coupling trees are shown in figure~\ref{fig-flop}, and corresponds to the following recoupling coefficient~\cite{Edm}~: \begin{eqnarray} \recoup{((a,b)d,c)f}{(a,(b,c)e)f} & \equiv & \recoup{(a,(b,c)e)f}{((a,b)d,c)f} \nonumber \\ & = & (-1)^{\pm(a+b+c+f)} \sqrt{(2d+1)(2e+1)} \sixj{a}{b}{d}{c}{f}{e}. \label{eq-flop} \end{eqnarray} \begin{figure}[htp] \begin{center} %%%%%\input flop.pic \begin{minipage}{12cm} \caption{Coupling trees for a flop~: $ \recoup{((a,b)d,c)f}{(a,(b,c)e)f} $.} \label{fig-flop} \end{minipage} \end{center} \end{figure} In order to obtain a summation formula for a recoupling coefficient, we start from the coupling tree of the initial state vector and try to transform it into the coupling tree of the final state vector, by applying a sequence of basic operations on subtrees of the original coupling tree. Each transformation step will work on some subtree, which matches the appropriate tree for the basic operation, in the following way~: for an exchange the subtree has to be of the form $(a,b)c$, where the nodes $a$ and $b$ can be coupled nodes or leaf nodes of the complete coupling tree; for a flop the subtree has to be of the form $((a,b)d,c)f$ or $(a,(b,c)e)f$, where the nodes $a$, $b$ and $c$ can be coupled nodes or leaf nodes of the complete coupling tree. Each transformation step will contribute some part in the summation formula~: For an exchange on a subtree $(a,b)c$ the only contribution is a sign factor of the following form~: \[ (-1)^{\pm(a+b-c)}. \] When applying a flop on a subtree $((a,b)d,c)f$, the node $d$ is replaced by a new node $e$ in the transformed subtree $(a,(b,c)e)f$. However, if the nodes $b$ and $c$ are already coupled to a node $j_i$ in the final coupling tree, then the node $e=j_i$ is not truly new. In that case the contribution in the summation formula is only a factor, consisting of a number and a 6-$j$ coefficient~: \[ (-1)^{\pm(a+b+c+f)} \sqrt{(2d+1)(2j_i+1)} \sixj{a}{b}{d}{c}{f}{j_i}. \] In the other case, the node $e$ is truly new, and gives rise to a new summation variable in the summation formula~: \[ \sum_e (-1)^{\pm(a+b+c+f)} \sqrt{(2d+1)(2e+1)} \sixj{a}{b}{d}{c}{f}{e}. \] The complete summation formula is then obtained as the consecutive product of all the contributions of the basic operations needed to transform the initial coupling tree to the final coupling tree. As an example we show how to obtain a summation formula for the following recoupling coefficient, which is known to be proportional to a 9-$j$ coefficient~\cite{Edm}~: \begin{eqnarray} \lefteqn{\recoup{((j_1,j_2)j_{12},(j_3,j_4)j_{34})j}{ ((j_1,j_3)j_{13},(j_2,j_4)j_{24})j}} \qquad & & \nonumber \\ & = &\sqrt{(2j_{12}+1)(2j_{34}+1)(2j_{13}+1)(2j_{24}+1)} \ninej{j_1}{j_2}{j_{12}}{j_3}{j_4}{j_{34}}{j_{13}}{j_{24}}{j}. \label{eq-9j1} \end{eqnarray} Table~\ref{tab-summa} lists a sequence of flops and exchanges which performs the transformation of the coupling tree $((j_1,j_2)j_{12},(j_3,j_4)j_{34})j$ to the coupling tree $((j_1,j_3)j_{13},(j_2,j_4)j_{24})j$, thereby denoting the newly introduced node by $k$. Figure~\ref{fig-summa} shows the complete coupling tree at each stage in the transformation process. \begin{table}[htp] \begin{center} {\footnotesize \begin{tabular}{|ll|l|} \hline \multicolumn{2}{|c|}{Operation on subtree \quad thus creating} & \multicolumn{1}{c|}{Contribution to summation formula} \\ \hline $((j_1,j_2)j_{12},j_{34})j$ & $(j_1,(j_2,j_{34})k)j$ & $ \sum_k (-1)^{j_1+j_2+j_{34}+j} \sqrt{(2j_{12}+1)(2k+1)} \sixj{j_1}{j_2}{j_{12}}{j_{34}}{j}{k} $ \\ $(j_2,j_{34})k$ & $(j_{34},j_2)k$ & $ (-1)^{-j_2-j_{34}+k} $ \\ $((j_3,j_4)j_{34},j_2)k$ & $(j_3,(j_4,j_2)j_{24})k$ & $ (-1)^{j_2+j_3+j_4+k} \sqrt{(2j_{34}+1)(2j_{24}+1)} \sixj{j_3}{j_4}{j_{34}}{j_2}{k}{j_{24}} $ \\ $(j_1,(j_3,j_{24})k)j$ & $((j_1,j_3)j_{13},j_{24})j$ & $ (-1)^{-j_1-j_3-j_{24}-j} \sqrt{(2j_{13}+1)(2k+1)} \sixj{j_1}{j_3}{j_{13}}{j_{24}}{j}{k} $ \\ $ (j_4,j_2)j_{24}$ & $(j_2,j_4)j_{24}$ & $ (-1)^{-j_2-j_4+j_{24}} $ \\ \hline \end{tabular} } \begin{minipage}{14cm} \caption{Sequence of flops and exchanges transforming the coupling tree $((j_1,j_2)j_{12},(j_3,j_4)j_{34})j$ to the coupling tree $((j_1,j_3)j_{13},(j_2,j_4)j_{24})j$.} \label{tab-summa} \end{minipage} \end{center} \end{table} This results in the following summation formula~: \begin{eqnarray} % \lefteqn{\recoup{((j_1,j_2)j_{12},(j_3,j_4)j_{34})j}{ % ((j_1,j_3)j_{13},(j_2,j_4)j_{24})j}} \qquad \nonumber \\ & & \sqrt{(2j_{12}+1)(2j_{34}+1)(2j_{13}+1)(2j_{24}+1)} \\ & & \times \sum_k (-1)^{2k} (2k+1) \sixj{j_1}{j_2}{j_{12}}{j_{34}}{j}{k} \sixj{j_3}{j_4}{j_{34}}{j_2}{k}{j_{24}} \sixj{j_1}{j_3}{j_{13}}{j_{24}}{j}{k}, \nonumber \end{eqnarray} yielding the well known expression for 9-$j$ coefficients in terms of 6-$j$ coefficients. \begin{figure}[htp] \begin{center} %%%%%\input summa.pic \begin{minipage}{12cm} \caption{Sequence of transformations for the recoupling coefficient $ \recoup{((j_1,j_2)j_{12},(j_3,j_4)j_{34})j}{((j_1,j_3)j_{13},(j_2,j_4)j_{24})j} $ } \label{fig-summa} \end{minipage} \end{center} \end{figure} % % ALGORITHM % ---------- % %\pagebreak \section{Algorithm description} \label{sec-algo} % In this section we develop an algorithm to generate a summation formula for a given recoupling coefficient, by constructing a sequence of basic operations, as described in section~\ref{sec-meth}, which transforms the initial state coupling tree into the final state coupling tree. By means of heuristic techniques, we will try to generate an {\em optimal\/} summation formula, which can be evaluated afterwards for actual $j_i$ values in an efficient way. An optimal summation formula will have a minimal number of summation variables. In the transformation process a new summation variable can be introduced by a flop operation, so the number of summation variables will depend on the number of flop operations in the transformation sequence. However not every flop operation introduces a new summation variable~: a flop operation involving a `new' node which is one of the (coupled) nodes of the final coupling tree, only contributes a factor containing a \mbox{6-$j$} coefficient. This means that the number of summation variables equals the number of flop operations in the transformation process, minus the number of coupled nodes in the final state coupling tree that do not appear in the initial coupling tree (i.c.\ the root node). So, generating a summation formula with a minimal number of summation variables corresponds to constructing a transformation sequence for the coupling trees with a {\em minimal number of flop operations}. Note that an arbitrary number of exchange operations can be allowed in the transformation sequence, since an exchange will never contribute a new summation variable. % \subsection{Constructing a transformation sequence} % In the following algorithm the basic idea is to try and obtain the coupled nodes of the final tree as fast as possible. This happens in a number of steps, as follows~: Each step starts from a {\em current\/} tree, which is obtained from the initial tree by a sequence of flop (and exchange) operations and which already contains some of the coupled nodes of the final tree. This current tree is then transformed in such a way that one more coupled node of the final tree is obtained, in a {\em minimal\/} number of flop (and an arbitrary number of exchange) operations. How this can be done, is described below. This transformed tree is then the current tree for the next step. Starting from the initial tree, we thus arrive at the final tree in a finite number of steps, i.e.\ the number of coupled nodes in the final tree not appearing in the initial tree. The number of flop operations needed for the complete transformation is the sum of the number of flop operations of all the algorithm steps. The fact that each `local' step obtains its result in a minimal number of flops, does not {\em guarantee\/} that the final result is obtained in a {\em minimal\/} number of flops, but this `global' minimum will be approached quite close. In each algorithm step we search for a sequence of flop (and exchange) operations in the following way~: First we try out {\em all possible single flop sequences\/} on the current tree, to see if any one of them obtains one of the coupled nodes of the final tree. By `single flop sequences' we mean flop operations of the form $ \recoup{((a,b)d,c)f}{(a,(b,c)e)f} $ (as defined in~(\ref{eq-flop})), but also operations of the form $ \recoup{((a,b)d,c)f}{(b,(a,c)e)f} $, which consist of an exchange followed by a flop operation. By `all possible single flop sequences on a tree' we mean the set of single flop sequences that can be performed on the subtrees of the tree. In the next subsection we describe how this can be obtained by traversing the coupling tree. By `obtaining a node of the final tree' we mean obtaining a node which has the same set of leaf nodes as one of the nodes in the final tree, since a coupled node is uniquely identified by the set of its leaf nodes. We can perform any succesful single flop sequence; this gives rise to a new factor (containing one 6-$j$ coefficient) in the summation formula, but not to a new summation variable. If no single flop sequence obtains a coupled node of the final tree, then we try all possible double flop sequences to see if any one of them results in a coupled node of the final tree. We can perform any successful double flop sequence, giving rise to two 6-$j$ coefficients and one summation variable. If no double flop sequence obtains a coupled node of the final tree, we try all possible sequences with three flops, and so on until at last we obtain a coupled node of the final tree. The next subsection describes in detail how this can be performed by a recursive algorithm. For such a sequence with $i$ flops, the contribution in the summation formula will be a product of $i$ 6-$j$ coefficients and $(i-1)$ summation variables. In each step of the algorithm there may be several minimal flop sequences which obtain a node of the final tree, and we are left with the problem of choosing one of them. At this stage all of them are `best choices', but at a later stage some of them might turn out to be better choices than others, in the sense that some of them might continue in shorter paths than others. The problem is to decide at this stage which sequences will later on give rise to the shortest paths. A rule guaranteeing success can not be given, but the idea is that those sequences which obtain a node at the highest possible level in the tree are the best choices, for the following reasons. When a node of the final tree is obtained, then `fixing' the rest of this subtree (i.e.\ obtaining its nodes in the final tree) involves only the subtree itself, while the rest of the complete tree can be ignored. Thus, in a way, part of the work is split off as an independent subtask of a similar nature. When preferring first to obtain nodes at a high level in the tree, the first subtasks splitting off from the original one, will be the larger ones, while later on smaller subtasks will be split off. This seems to be a more stable approach than splitting off smaller subtasks first. To illustrate the way in which the above algorithm searches for a transformation sequence, we show in detail some of the algorithm steps in generating a summation formula for the recoupling coefficient~(\ref{eq-exa}) shown in figure~\ref{fig-tree}. Starting from the initial tree $((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9$, no single flop sequence obtains a coupled node of the final tree $(((j_1,j_4)j_{10},(j_2,j_3)j_{11})j_{12},j_5)j_9$, however several double flop sequences do, such as~: \begin{eqnarray} \lefteqn{((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9} \qquad & & \nonumber \\ & \to & (((j_1,j_2)j_6,j_3)k_1,(j_4,j_5)j_7)j_9 \to ((((j_1,j_2)j_6,j_3)k_1,j_4)j_{12},j_5)j_9, \label{eq-sel1} \\ {\rm or} \nonumber \\ & \to & (j_1,(j_2,(j_3,(j_4,j_5)j_7)j_8)k_1)j_9 \to (j_1,((j_2,j_3)j_{11},(j_4,j_5)j_7)k_1)j_9, \label{eq-sel2} \end{eqnarray} the latter of which obtains the node $j_{11}$ of the final tree, while the first sequence obtains the node $j_{12}$ (which is the node with $\{j_1,j_2,j_3,j_4\}$ as the set of descendent leaf nodes). Continuing from the sequence~(\ref{eq-sel1}), the final tree can be obtained by two more algorithm steps, each consisting of a single flop sequence~: \begin{eqnarray*} ((((j_1,j_2)j_6,j_3)k_1,j_4)j_{12},j_5)j_9 & \to & (((j_1,(j_2,j_3)j_{11})k_1,j_4)j_{12},j_5)j_9 \\ & \to & (((j_1,j_4)j_{10},(j_2,j_3)j_{11})j_{12},j_5)j_9. \end{eqnarray*} The resulting summation formula is a single sum over a product of four 6-$j$ coefficients~: \begin{eqnarray} & & (-1)^{2j_1+j_2+j_5+3j_6+j_7+2j_9-j_{10}} \nonumber \\ & & \times \sqrt{(2j_6+1)(2j_7+1)(2j_8+1)(2j_{10}+1)(2j_{11}+1)(2j_{12}+1)} \label{eq-exfor1} \\ & & \times \sum_{k_1} (-1)^{k_1+j_1+j_{11}} (2k_1+1) \sixj{j_7}{j_3}{j_8}{j_6}{j_9}{k_1} \sixj{j_5}{j_4}{j_7}{k_1}{j_9}{j_{12}} \nonumber \\ & & \mbox{} \qquad \qquad \times \sixj{j_1}{j_2}{j_6}{j_3}{k_1}{j_{11}} \sixj{j_{11}}{j_1}{k_1}{j_4}{j_{12}}{j_{10}}. \nonumber \end{eqnarray} However, if we continue from the sequence~(\ref{eq-sel2}), the next algorithm step again requires a double flop sequence to obtain a node of the final tree, e.g.~: \begin{eqnarray*} \lefteqn{(j_1,((j_2,j_3)j_{11},(j_4,j_5)j_7)k_1)j_9} \qquad & & \\ & \to & ((j_1,(j_2,j_3)j_{11})k_2,(j_4,j_5)j_7)j_9 \to (((j_1,(j_2,j_3)j_{11})k_2,j_4)j_{12},j_5)j_9. \end{eqnarray*} Continuing from this result, the final tree can be obtained in a single flop sequence. The resulting summation formula is a double sum over a product of five 6-$j$ coefficients~: \begin{eqnarray} & & (-1)^{j_2+j_3+j_5+2j_7+3j_9-j_{10}+2j_{11}} \nonumber \\ & & \times \sqrt{2j_6+1)(2j_7+1)(2j_8+1)(2j_{10}+1)(2j_{11}+1)(2j_{12}+1)} \label{eq-exfor2} \\ & & \times \sum_{k_1} \sum_{k_2} (-1)^{k_1+2k_2+j_2+j_8} (2k_1+1) (2k_2+1) \sixj{j_1}{j_2}{j_6}{j_8}{j_9}{k_1} \sixj{j_7}{j_3}{j_8}{j_2}{k_1}{j_{11}} \nonumber \\ & & \mbox{} \qquad \qquad \times \sixj{j_7}{j_{11}}{k_1}{j_1}{j_9}{k_2} \sixj{j_5}{j_4}{j_7}{k_2}{j_9}{j_{12}} \sixj{j_{11}}{j_1}{k_2}{j_4}{j_{12}}{j_{10}}. \nonumber \end{eqnarray} It is obvious that, in the first step of the algorithm, the sequence~(\ref{eq-sel1}) obtaining the node $j_{12}$ is a better choice than the sequence~(\ref{eq-sel2}) obtaining the node $j_{11}$. Note that the node $j_{12}$ is at a higher level in the final coupling tree than the node $j_{11}$. It is also obvious that the single sum formula is optimal, since in the first algorithm step no single flop sequence is successful, which means that at least one summation variable needs to be introduced. % \subsection{Traversing a coupling tree} % Trying out {\em all possible single flop sequences\/} on a coupling tree is done by traversing the tree in a systematic way, starting from the root node, and examining each node when it is visited. Each coupled node is examined in the following way~: On the subtree with this node as the root node, we can consider four possible single flop sequences, i.e.\ of the following forms~: \begin{quote} $\recoup{((a,b)d,c)f}{(a,(b,c)e)f}$, \qquad $\recoup{(a,(b,c)d)f}{((a,b)e,c)f}$, \\ $\recoup{((a,b)d,c)f}{(b,(a,c)e)f}$, \qquad $\recoup{(b,(a,c)d)f}{((a,b)e,c)f}$, \end{quote} where $f$ denotes the current node. For each of these operations we check if it can be performed on the current subtree, i.e.\ if this subtree `matches' the appropriate form for the operation. If so, we check if the node $e$, that is `created' in the operation, is one of the nodes of the final coupling tree, that has not yet been obtained. If so, a successful single flop sequence is found and the search process can be stopped, otherwise the search process has to be continued. This algorithm can easily be extended to {\em try out all possible $N$ flop sequences}, with a given number $N$ of flops. For that purpose we view an $N$ flop sequence as a single flop sequence followed by an $(N-1)$ flop sequence, a recursive formulation giving rise to the following algorithm. For every coupled node in the tree, we try to find a successful $N$ flop sequence, starting with one of the four possible single flop operations on the current node. This involves checking for each of these four operations if it can be performed on the current subtree, and if so, trying out all $(N-1)$ flop sequences on that subtree (which is done recursively). When $N=1$, we are at the final stage in the recursive process, where we have to test if the current flop sequence obtains a new node of the final coupling tree. Again, if so, a successful single flop sequence is found and the search process can be stopped, otherwise the search process has to be continued. A search process traversing a tree, continues until either the examination of the current node is successful or all the nodes have been examined. If all the nodes have been examined without finding a successful flop sequence of given length, then the search is said to be unsuccessful and we have to continue by trying out longer flop sequences. We are now left with the problem of implementing a search process, which has to examine every node of a tree. Two possible approaches for such a tree traversal are commonly used, i.e.\ so-called {\em depth-first\/} or {\em breadth-first traversal\/}~\cite{DSA}. For the problem considered here we will use a breadth-first traversal, which starts from the root node, then considers every coupled node at the first level of the tree, then the nodes at the second level, etc., until every level in the tree is considered. E.g., in the tree $((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9$ the coupled nodes at the first level are $j_6$ and $j_8$; on the second level the only coupled node is $j_7$. When implementing a breadth-first traversal, a problem is caused by the fact that a node in the tree only contains references to its children, while the information needed is a reference to the next node on the same level in the tree (or to the first node on the next level, in case it is the last node on a level), which is not contained in the node. But this information can be built and stored during the traversal process, in the following way. Assuming that we know which are the nodes on a certain level, we can build a list of nodes which are on the next level, during the process of visiting the nodes on this level; indeed, they are the left and right children of the nodes on the current level. Thus during the traversal process we will keep track of (1) the nodes on the current level which still have to be visited, and (2) the nodes on the next level which are children of visited nodes on the current level and which have to be visited when this level is finished. To store this information we can use the common concept of a {\em queue}, an abstract datastructure simulating real-life waiting queues, where the first customer in the queue is served first and new customers join at the tail of the queue. At each stage in the traversal the first node is taken from the queue and it is examined. If this is not successful, then the left and right (or right and left) children, if they are coupled nodes, are added at the tail of the queue (to be examined at some later stage). The search continues on the (next) first node in the queue, which will be the next node on the same level in the tree or the first node of the next level in the tree. E.g., a breadth-first traversal, traversing each level from left to right, on the tree $((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9$, visiting only the coupled nodes, proceeds as shown in table~\ref{tab-bf}. \begin{table}[htp] \begin{center} \begin{tabular}{|c|c|l|} \hline step & node & queue \\ \hline 0 & & $(j_9)$ \\ 1 & $j_9$ & $(j_6,j_8)$ \\ 2 & $j_6$ & $(j_8)$ \\ 3 & $j_8$ & $(j_7)$ \\ 4 & $j_7$ & $()$ \\ \hline \end{tabular} \begin{minipage}{14cm} \caption{Breadth-first traversal of the tree $((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9$, visiting only the coupled nodes.} \label{tab-bf} \end{minipage} \end{center} \end{table} % % IMPLEMENTATION DETAILS % ---------------------- % %\pagebreak \section{Implementation details} \label{sec-impl} In the program {\tt NJFORMULA} as we present it here, the algorithms for {\em generating\/} a summation formula, as described above, are implemented in a set of routines and these are combined in a module, which is used from a main program. The program is designed such that afterwards it can easily be linked with the program {\tt NJSUMMATION}~\cite{FPV}, which will supply routines to {\em evaluate\/} the formulae generated by the current program. The module implemented here supplies first of all a function {\tt generate\_formula} to generate a summation formula for a given recoupling coefficient. Apart from that it also supplies procedures for the input of a recoupling coefficient in several forms ({\tt read\_expr} and {\tt read\_njsym}), as well as a procedure {\tt write\_formula} for the output of a summation formula. In this section we describe some implementation details as well as some special features concerning the functionality of these routines. We also discuss in detail a some examples of simple main programs, which illustrate how to use the described routines to generate summation formulae for a number of general recoupling coefficients. % \subsection{Datastructures} % Throughout the program the index $i$ of a $j_i$ coefficient will be used to identify the $j_i$ coefficient, while a negative index $-i$ will be used to identify a summation variable $k_i$ of the summation formula. A coupling tree is represented by a datastructure, where each node of the tree is a record structure {\tt NODE}, with a field {\tt index} which contains the index $i$ of the $j_i$ coefficient or $-i$ of the $k_i$ summation variable of the node, fields {\tt left} and {\tt right} containing pointers to the left child and the right child of the node in the tree (leaf nodes, which have no children, contain null pointers), and a field {\tt leafs} containing the set of indices of the $j_i$ coefficients of the leaf nodes which are descendants of this node. E.g., in the coupling tree for the state vector $\ket{((j_1,j_2)j_6, (j_3,(j_4,j_5)j_7)j_8)j_9}$ the record for the node $j_8$ contains its index 8, the set $\{3,4,5\}$, a left pointer to the node $j_3$ and a right pointer to the node $j_7$. A summation formula is represented by a record structure {\tt FORMULA} containing the following data~: The fields {\tt nrjs}, {\tt nrks}, {\tt nrsixjs} contain resp.\ the number of $j_i$ coefficients in the recoupling coefficient, the number of summation variables $k_i$ and the number of 6-$j$ coefficients in the products of the formula. Fields {\tt jsigns}, resp.\ {\tt ksigns}, and {\tt jsqrts}, resp.\ {\tt ksqrts[i]}, are arrays with an entry for every $j_i$ coefficient, resp.\ $k_i$ summation variable, such that \begin{quote} {\tt jsigns[i]}, resp.\ {\tt ksigns[i]} = power of $(-1)^{j_i}$, resp.\ $(-1)^{k_i}$, in the formula, \\ {\tt jsqrts[i]}, resp.\ {\tt ksqrts[i]} = power of $\sqrt{(2 j_i + 1)}$, resp.\ $\sqrt{(2 k_i + 1)}$. \end{quote} The field {\tt sixjs} is an array for which every entry is a representation of one of the 6-$j$ coefficients of the formula; a 6-$j$ coefficient is represented as an array with 6 elements, which are $j_i$ coefficients or $k_i$ summation variables; a $j_i$ coefficient is represented by its index $i$, while a $k_i$ summation variable is represented by its negative index $-i$. E.g., for the double summation~(\ref{eq-exfor2}), obtained above, the structure contains the following data~: \begin{center} {\small \begin{tabular}{rlrl} {\tt nrjs} \quad = & 12 \qquad \qquad {\tt nrks} \quad = \quad 2 & {\tt nrsixjs} \quad = & 5 \\[0.25cm] {\tt jsigns} \quad = & $(0,2,1,0,1,0,2,1,3,-1,2,0)$ & {\tt ksigns} \quad = & $(1,2)$ \\ {\tt jsqrts} \quad = & $(0,0,0,0,0,1,1,1,0,1,1,1)$ & {\tt ksqrts} \quad = & $(2,2)$ \\[0.25cm] {\tt sixjs[1]} \quad = & $(1,2,6,8,9,-1)$ & {\tt sixjs[2]} \quad = & $(7,3,8,2,-1,11)$ \\ {\tt sixjs[3]} \quad = & $(7,11,-1,1,9,-2)$ & {\tt sixjs[4]} \quad = & $(5,4,7,-2,9,12)$ \\ {\tt sixjs[5]} \quad = & $(11,1,-2,4,12,10)$ \end{tabular}} \end{center} % % INPUT OF RECOUPLING COEFFIENTS % \subsection{Input of recoupling coefficients} Several forms of input for recoupling coefficients are possible. For reasons of compatibility with the programs {\tt NJSYM}~\cite{Bur} and {\tt NJGRAF}~\cite{Bar}, we have implemented a procedure {\tt read\_njsym}, which reads a recoupling coefficient in the form used by these programs. It returns the number of $j_i$ coefficients read and the two pointers to the corresponding coupling trees. However, we also propose another input style for a recoupling coefficient, which is well suited for interactive input. The procedure {\tt read\_expr} reads a recoupling coefficient in the form of an expression similar to the mathematical bra-ket notation. It also returns the number of $j_i$ coefficients read and the two pointers to the corresponding coupling trees. Both input routines read their input from the standard input device {\tt stdin}, which is usually the keyboard. Reading input from a file can be done easily by `redirecting' the standard input device when executing the program. % \subsubsection{Input in the style of {\tt NJSYM}} % In the programs {\tt NJSYM} and {\tt NJGRAF} the input of a recoupling coefficient consists of the number of $j_i$ coefficients, the number of couplings, and two lists of triads $(j_{i_1},j_{i_2},j_{i_3})$ describing the coupling in the initial and final state vector, where a triad $(j_{i_1},j_{i_2},j_{i_3})$ signifies that $j_{i_1}$ and $j_{i_2}$ are coupled to $j_{i_3}$. To read a recoupling coefficient given in this form, the following actions have to be performed. In order to build the coupling trees corresponding to this input, the common root node of the trees has to be determined first. We assume that for each state vector the triad describing the root coupling is given as the last triad in the list. Next a coupling tree has to be built for each of the state vectors. This is implemented in a function {\tt build\_tree}, which returns a pointer to the root node of the tree. Building a coupling tree from a list of triads, given the $j_i$ coefficient of the root node, can easily be done by a recursive algorithm, in the following way~: When memory space for the node is allocated, the index $i$ of the root node $j_i$ is filled in. Next the list of triads is searched to find the $j_{i_1}$ and $j_{i_2}$ coefficients that are coupled to $j_i$. If no such coefficients are found, then $j_i$ is a leaf node, pointing nowhere left and right and the set of leaf nodes is simply the singleton $\{j_i\}$. Else $j_i$ is a coupled node, pointing left to a subtree with $j_{i_1}$ as the root node and pointing right to a subtree with $j_{i_2}$ as the root node. The building of these subtrees is done in the same way, i.e.\ by two recursive calls to the function {\tt build\_tree}, returning the pointers to be assigned to the fields {\tt left} and {\tt right} of the node $j_i$. The set of leaf nodes for $j_i$ is obtained as the union of the set of leaf nodes of its left and right subtrees $j_{i_1}$ and $j_{i_2}$. % \subsubsection{Input in the form of expressions} % Another representation of recoupling coefficients, well suited for input and output, is an expression, which closely resembles the bra-ket notation of a recoupling coefficient. The coupling of $j_{i_1}$ and $j_{i_2}$ to $j_{i_3}$ is expressed as $(i_1,i_2)i_3$. In their turn $j_{i_1}$ and $j_{i_2}$ can be the couplings of other $j_i$ coefficients, which is expressed in the same way. This gives rise to a string consisting of nested couplings, describing the coupling tree for a state vector. Two such state vector expressions, enclosed in brackets ({\tt <}) and ({\tt >}) and separated by a vertical bar ({\tt |}), then form the input expression for a general recoupling coefficient. E.g., the input expression for the recoupling coefficient~(\ref{eq-exa}) is~: \begin{quote} {\tt < ((1,2)6,(3,(4,5)7)8)9 | (((1,4)10,(2,3)11)12,5)9 > }. \end{quote} The handling of such an input expression is implemented in the procedure {\tt read\_expr}. A state vector expression can be read and transformed into a coupling tree by means of a recursive algorithm, which we have implemented in a function {\tt read\_tree}. To that purpose, we formally define a state vector as either a positive number, or an expression of the form ( {\sl state vector\/} , {\sl state vector\/} ) {\sl number\/}. The latter is a recursive definition, saying that a state vector can consist of two state vectors, separated by a `,' and enclosed in brackets, followed by a number. The number is the index of the root node in the coupling tree of the state vector, while the two state vectors are its left and right children. From this definition the algorithm to build a coupling tree corresponding to a given state vector expression, which can be viewed as a string of tokens, follows easily. A general step in the algorithm works on an expression corresponding to a subtree, starting from a given node, and builds this subtree in the following way~: Memory space is allocated for the root node of the subtree, a pointer to which will be returned by the algorithm when it finishes. If the next token in the input string is a number, then this subroot node is a leaf node, for which the index is this number, the set of leaf nodes is the singleton containing this number and the left and right pointers are set to null. Else, if the next token in the input string is an opening bracket `(', then this subroot node is a coupled node. This means that the next part of the input string is an expression corresponding to its left child subtree, which can be built by a recursive call to the same algorithm. After that the next token in the input string must be a comma and the next part will be an expression corresponding to the right child subtree of the subroot node, again to be built recursively. Next in the input string will be a closing bracket `)', followed by a number, which is the index of the subroot node. Finally the set of leaf nodes of the subroot node is obtained as the union of the set of leaf nodes of the left and right subtrees. During the recursive process the input is checked for errors, such as forgotten commas or non-matching brackets, in which case we finish the recursive process and report it to be unsuccessful. In this implementation the recursive function {\tt read\_tree} is called twice from the procedure {\tt read\_expr}, to read the bra-vector and the ket-vector of the recoupling coefficient. Also here the input is checked for errors; if errors occur, the procedure allows to retry the input of an expression, until a correct input expression is given. Another feature we have implemented, to allow simplification of the input expression, is automatic numbering of the coupled nodes. In that case only the indices of the leaf nodes have to be given explicitly, while for the coupled nodes an index is generated automatically during the reading process by means of consecutive integers. E.g., for the recoupling coefficient~(\ref{eq-exa}) the simplified input expression is \begin{quote} {\tt < ((1,2),(3,(4,5))) | (((1,4),(2,3)),5) > } \end{quote} and the indices for the coupled nodes generated by the program will be the same as in the input expression given above. This automatic numbering is only possible if the number of leaf nodes is given first. The option to use automatic numbering or not is left to the user and is controlled by a global variable {\tt AUTONUM}, which can be set to {\tt ON} or {\tt OFF} (by default it is set to {\tt ON}). For a better readability of an input expression, we also allow separators (such as spaces, tabulators and newline characters) to be put arbitrarily in between the tokens. % \subsubsection{Remark} When a recoupling coefficient is given, it can happen that some of the coupled nodes in both trees are the same, in the sense that they have the same set of leaf nodes, while the index $i$ of the $j_i$ coefficient for the nodes is different in the two trees. In such cases the value of the recoupling coefficient will be trivially zero, if the corresponding $j_i$ coefficients are not the same. To avoid unnecessary computations (in evaluating the summation formula, or sometimes even in generating a summation formula), we want to report on these cases when reading a recoupling coefficient, during which process they can easily be detected. For that purpose we build a list of delta functions $\delta(j_{i_1},j_{i_2})$, that is returned as an extra parameter by the input routines. In the program presented here we will not yet use this information, as here we are mainly dealing with the algorithm for the generation of a summation formula, independent of the actual $j_i$ values. But the program {\tt NJSUMMATION}, which will deal with the actual $j_i$ values, will indeed need this information to check if any computations are actually needed for a given recoupling coefficient. For that reason we have already implemented this feature in the input routines. % % GENERATION OF A SUMMATION FORMULA % \subsection{Generation of a summation formula} The function {\tt generate\_formula} searches a transformation sequence for a recoupling coefficient, given by the corresponding coupling trees and the number of $j_i$ coefficients, and returns a pointer to the generated summation formula. If for some reason the search was not successful, then a null pointer is returned. It is a rather straightforward implementation of the algorithms described in section~\ref{sec-algo}. The main body is a loop, which searches in each step for a minimal flop sequence obtaining one more coupled node of the final coupling tree, until all the nodes of the final coupling tree are obtained. Each step in the main loop is in its turn a loop trying out all $N$ flop sequences, for $N=1,2,3,\ldots$, until a successful flop sequence has been found. Trying out all $N$ flop sequences, for given $N$, is implemented in a recursive algorithm as described earlier, examining every node of the coupling tree in a breadth-first traversal. For practical purposes we have set an upper limit for the values of $N$ to be tried out, which is kept in a constant {\tt MAXFLEN} and is here set to~5 (in the examples described below no higher value was needed). To keep track of which nodes of the final coupling tree are already obtained, we make use of the field {\tt leafs} in each node, defined to contain the set of descendent leaf nodes of the node. When a node of the final coupling tree is obtained, we set its field {\tt leafs} to be the empty set, thus making clear that in a next step one should not try to obtain this node anymore. When all the nodes of the final coupling tree have their leafset set to empty, then the transformation process is finished. Before starting the transformation process, a sweep of the final coupling tree is performed, checking for each node if it already occurs in the initial coupling tree (by checking for equal leafsets). If so, for such a node in the final tree, the field {\tt leafs} is set to empty. The root node will be such a node. In each step the contribution of the successful flop sequence is recorded in the datastructure that will finally contain the complete summation formula. For every flop operation in the sequence, the appropriate entries in the {\tt jsigns}, {\tt ksigns}, {\tt jsqrts} and {\tt ksqrts} arrays are updated accordingly, and a new 6-$j$ coefficient is added. Every exchange operation in the sequence will influence some entries in the {\tt jsigns}, {\tt ksigns}, {\tt jsqrts} and {\tt ksqrts} arrays. Apart from updating the summation formula in generation, the basic operations in the successful sequence also involve updating of the current coupling tree in each step. In an exchange operation the left and right children of a node are interchanged, which can be done easily by interchanging the {\tt left} and {\tt right} pointer in the node. A flop operation of the form $\recoup{((a,b)d,c)f}{(a,(b,c)e)f}$, working on the current node $f$, changes the structure of the subtree more drastically~: the left child $d$ is replaced by its own left child $a$, while the right child $c$ is replaced by a new node $e$, which has the node $c$ as its right child and the node $b$ as its left child. This transformation involves only updating of the appropriate pointers, no copying of subtrees is necessary. The other flop operations involve similar actions. Finally we would like to make the following remark. Trying out several examples, we noticed that in some cases a different result was obtained when traversing the tree levels from left to right or from right to left. Also trying to transform the final tree into the initial tree instead of vice versa, sometimes produced a better result. For that reason, the function {\tt generate\_formula} tries out several of these possible search processes, and keeps the best result that is obtained in this way. The number of possibilities tried out is controlled by a global variable {\tt SEARCH}, which can be set by the user to either {\tt SHALLOW}, {\tt DEEP} or {\tt VERYDEEP}; a typical choice is {\tt DEEP}, which is also the default value. No matter what the value of {\tt SEARCH} is, the search always stops if a formula with less than two summation variables is obtained, since such a formula is guaranteed to be optimal. For debugging purposes we have also implemented an option to let the program output some intermediate results, which can be used to follow the generation process. This debug information is written to a file, the name of which can be chosen. An overview of each search process is given as a list of coupling trees (in the form of an expression) obtained in each step of a search process, as well as the number of flops in the sequence obtaining the current coupling tree in a step. For each search process the obtained formula is also given. Debugging can be initialised by calling the procedure {\tt init\_debug}, giving the filename as a parameter, and must be closed by calling a procedure {\tt close\_debug}. % % MAIN PROGRAM % \subsection{Main program} The following simple main program can be used to generate summation formulae for an arbitrary number of recoupling coefficients, which are read interactively in the form of expressions. The headerfile {\tt njformul.h} contains a set of declarations, which are needed by the program. The global variable {\tt AUTONUM} is initialised to {\tt ON}, meaning that only leaf nodes of the recoupling coefficients will be given explicitly. The global variable {\tt SEARCH} is initialised to {\tt DEEP}, signifying that both traversing left-to-right and right-to-left will be tried out in a transformation of bra to ket. The {\tt do\ldots\,while} loop reads in each step a new recoupling coefficient, generates the corresponding summation formula and writes this formula. This loop continues until an answer {\tt n} or {\tt N} is given to the question {\tt Continue(y/n)?}. {\small \begin{verbatim} #include "njformul.h" main () { int nrjs; NODE *bra, *ket; DELTAS delta; FORMULA *form; char ch; AUTONUM = ON; SEARCH = DEEP; do { read_expr (&nrjs, &bra, &ket, &delta); form = generate_formula (nrjs, bra, ket); write_formula (form); printf ("Continue (y/n) ? "); scanf ("%c", &ch); } while (ch != 'n' && ch != 'N'); } \end{verbatim}} The following main program reads one recoupling recoupling coefficient in the style of {\tt NJSYM} and writes the generated formula. Debugging is initialised and the debug information will be written to a file {\tt formula.deb}. The variable {\tt SEARCH} is set to {\tt VERYDEEP}, meaning that both left-to-right and right-to-left traversal will be tried out, for transforming bra to ket as well as ket to bra. The variable {\tt AUTONUM} is not significant in this case. {\small \begin{verbatim} #include "njformul.h" main () { int nrjs; NODE *bra, *ket; DELTAS delta; FORMULA *form; init_debug ("formula.deb"); SEARCH = VERYDEEP; read_njsym (&nrjs, &bra, &ket, &delta); form = generate_formula (nrjs, bra, ket); write_formula (form); close_debug (); } \end{verbatim}} % % RESULTS % ------- % %\pagebreak \section{Results and discussion} \label{sec-res} The program {\tt NJFORMULA} is written in C and can be run on any system providing a C (or C++) compiler. We have tested the program on a range of machines (such as 386/486-based PCs running MS-DOS and Linux~\cite{Linux}, and a Sun Sparc running Unix) using several compilers (such as Turbo C++~\cite{TC}, Gnu CC~\cite{GnuC} and SPARCompiler~C). Since this program will be joined later with the program {\tt NJSUMMATION}, which will include a parallel implementation on a T800 transputer system (hosted by PCs under MS-DOS), we have also tested the program {\tt NJFORMULA} on this transputer system, using the 3L Parallel C compiler~\cite{ParC}. Of course the (sequential) program uses none of the parallel features of this compiler, and it is executed on a single transputer. To compare the results obtained by our program {\tt NJFORMULA} with the results obtained by the program {\tt NJGRAF}, we have compiled the program {\tt NJGRAF} (omitting the calls to the subroutine {\tt GENSUM}) on the Sun Sparc using the Sparc Fortran-77 compiler and on the transputer system using the 3L Parallel Fortran compiler~\cite{ParF}. First we want to compare our results for the test cases 1, 2 and 4 discussed by Bar-Shalom and Klapisch~\cite{Bar}, with the results of their program {\tt NJGRAF}. These test cases correspond to the following recoupling coefficients (note that ($G_1$) corresponds to the recoupling coefficient~(\ref{grc}))~: {\small \begin{itemize} \item[($G_1$)] $\recoup{((j_1,j_2)j_5,(j_3,j_4)j_6)j_7}{(j_1,((j_2,j_3)j_8,j_4)j_9)j_7}$ \item[($G_2$)] $\bra{(((j_1,j_2)j_8,((j_3,j_5)j_9,(j_6,j_7)j_{10})j_{11})j_{12},j_4)j_{13}}$ \\ $\ket{(((j_1,j_2)j_{14},((j_3,(j_7,j_5)j_{15})j_{16},j_6)j_{17})j_{18},j_4)j_{13}}$ \item[($G_4$)] $\bra{((j_1,j_2)j_{12},(j_3,(j_4,(j_5,((j_6,j_7)j_{13},(j_8,(j_9,(j_{10},j_{11})j_{14})j_{15})j_{16})j_{17})j_{18})j_{19})j_{20})j_{21}}$ \\ $\ket{(j_5,(j_7,((j_6,(j_{11},(j_9,(j_{10},j_8)j_{22})j_{23})j_{24})j_{25},(j_1,(j_3,(j_2,j_4)j_{26})j_{27})j_{28})j_{29})j_{30})j_{21}}$ \end{itemize} } Table~\ref{tab-njgf} shows for each test case the number of angular momenta involved, the number of summation variables and the number of 6-$j$ coefficients in the generated summation formula, for both programs {\tt NJGRAF} and {\tt NJFORMULA}. To give an idea of the complexity of the search process in the program {\tt NJFORMULA}, we also give the lengths of the minimal flop sequences obtained in each step of the generation algorithm. For each case {\tt NJFORMULA} obtains a summation formula which is as good as the one obtained by {\tt NJGRAF}, i.e.\ a formula with the same number of summation variables and 6-$j$ coefficients. Running times of both programs are comparable and are of the order of 0.01~s on the Sun Sparc and on the 486~PC (under Linux). Apart from these test cases described in~\cite{Bar}, we also discuss the results for some other recoupling coefficients, some of which give rise to quite a complicated search process~: {\small \begin{itemize} \item[($F_0$)] $\recoup{(j_1,j_2)j_5,(j_3,j_4)j_6)j_7}{((j_1,j_3)j_8,(j_2,j_4)j_9)j_7}$ \item[($F_1$)] $\bra{((j_1,j_2)j_6,(j_3,(j_4,j_5)j_7)j_8)j_9}$ $\ket{(((j_1,j_4)j_{10},(j_2,j_3)j_{11})j_{12},j_5)j_9}$ \item[($F_2$)] $\bra{((j_1,(j_2,j_3)j_7)j_8,(j_4,(j_5,j_6)j_9)j_{10})j_{11}}$ $\ket{(((j_1,(j_4,j_5)j_{12})j_{13},j_2)j_{14},(j_3,j_6)j_{15})j_{11}}$ \item[($F_3$)] $\bra{(((j_1,j_2)j_7,(j_3,j_4)j_8)j_9,(j_5,j_6)j_{10})j_{11}}$ $\ket{((j_3,j_6)j_{12},((j_2,j_4)j_{13},(j_1,j_5)j_{14})j_{15})j_{11}}$ \item[($F_4$)] $\bra{(((j_1,j_2)j_7,(j_3,j_4)j_8)j_9,(j_5,j_6)j_{10})j_{11}}$ $\ket{((j_1,j_6)j_{12},((j_3,j_5)j_{13},(j_2,j_4)j_{14})j_{15})j_{11}}$ \item[($F_5$)] $\bra{(((j_1,j_2)j_8,(j_3,j_4)j_9)j_{10},((j_5,j_6)j_{11},j_7)j_{12})j_{13}}$ \\ $\ket{((j_1,j_7)j_{14},(((j_3,j_5)j_{15},(j_2,j_4)j_{16})j_{17},j_6)j_{18})j_{13}}$ \item[($F_6$)] $\bra{(((j_1,j_2)j_8,(j_3,j_4)j_9)j_{10},(j_5,(j_6,j_7)j_{11})j_{12})j_{13}}$ \\ $\ket{(j_5,((j_6,(j_2,j_4)j_{14})j_{15},(j_3,(j_1,j_7)j_{16})j_{17})j_{18})j_{13}}$ \item[($F_7$)] $\bra{((j_1,(j_2,j_3)j_8)j_9,((j_4,j_5)j_{10},(j_6,j_7)j_{11})j_{12})j_{13}}$ \\ $\ket{(((j_1,j_4)j_{14},j_6),j_{15},((j_5,j_2)j_{16},(j_7,j_3)j_{17})j_{18})j_{13}}$ \item[($F_8$)] $\bra{(((j_1,j_2)j_9,j_3)j_{10},(((j_4,j_5)j_{11},j_6)j_{12},(j_7,j_8)j_{13})j_{14})j_{15}}$ \\ $\ket{(((j_1,j_4)j_{16},j_7)j_{17},(((j_2,j_5)j_{18},(j_8,j_3)j_{19})j_{20},j_6)j_{21})j_{15}}$ \item[($F_9$)] $\bra{(((j_1,(j_2,j_3)j_{11})j_{12},((j_4,j_5)j_{13},j_6)j_{14})j_{15},(((j_7,j_8)j_{16},j_9)j_{17},j_{10})j_{18})j_{19}}$ \\ $\ket{(((j_2,j_4)j_{20},j_7)j_{21},((((j_1,j_8)j_{22},(j_9,j_5)j_{23})j_{24},j_{10})j_{25},(j_6,j_3)j_{26})j_{27})j_{19}}$ \end{itemize} } Table~\ref{tab-njgf} also shows the results obtained for these recoupling coefficients, for both programs {\tt NJGRAF} and {\tt NJFORMULA}. For the simplest cases, i.e.~($F_0$)--($F_6$), the running time of both {\tt NJFORMULA} and {\tt NJGRAF} is again of the order of 0.01~s. For some of these cases, i.e.~($F_0$), ($F_2$), ($F_3$) and ($F_6$), both programs obtain a similar formula, with the same number of summation variables. However, for several cases, i.e.~($F_1$), ($F_4$) and ($F_5$), {\tt NJFORMULA} obtains a better formula than {\tt NJGRAF}, with less summation variables. Note that ($F_1$) corresponds to the recoupling coefficient~(\ref{eq-exa}) and that {\tt NJFORMULA} obtains the single sum~(\ref{eq-exfor1}), while {\tt NJGRAF} obtains a double summation similar to formula~(\ref{eq-exfor2}). For ($F_5$) we even get a double summation formula, where {\tt NJGRAF} obtains a 4-fold summation. For the most complicated cases, i.e.~($F_7$)--($F_9$), {\tt NJFORMULA} still generates a summation formula in very little time: ($F_7$) takes about 0.02~s, while ($F_8$) takes 0.15~s and ($F_9$) takes 0.50~s. However for these cases {\tt NJGRAF} gives no result at all~: the program crashes during the search process. \begin{table}[htp] \begin{center} {\small \begin{tabular}{|c|c|cc|ccl|} \hline & & \multicolumn{2}{c|}{\tt NJGRAF} & \multicolumn{3}{c|}{\tt NJFORMULA} \\ Test case & \#$j$'s & \#$k$'s & \#6-$j$'s & \#$k$'s & \#6-$j$'s & length of sequences \\ \hline ($G_1$) & 9 & 0 & 2 & 0 & 2 & 1-1 \\ ($G_2$) & 18 & 0 & 2 & 0 & 2 & 1-1 \\ ($G_4$) & 30 & 3 & 11 & 3 & 11 & 1-2-1-1-2-1-2-1 \\ \hline ($F_0$) & 9 & 1 & 3 & 1 & 3 & 2-1 \\ ($F_1$) & 12 & 2 & 5 & 1 & 4 & 2-1-1 \\ ($F_2$) & 15 & 1 & 5 & 1 & 5 & 1-2-1-1 \\ ($F_3$) & 15 & 2 & 6 & 2 & 6 & 2-2-1-1 \\ ($F_4$) & 15 & 3 & 7 & 2 & 6 & 2-2-1-1 \\ ($F_5$) & 18 & 4 & 9 & 2 & 7 & 2-2-1-1-1 \\ ($F_6$) & 18 & 2 & 7 & 2 & 7 & 1-2-2-1-1 \\ ($F_7$) & 18 & -- & -- & 4 & 9 & 3-2-2-1-1 \\ ($F_8$) & 21 & -- & -- & 7 & 13 & 3-4-1-3-1-1 \\ ($F_9$) & 27 & -- & -- & 9 & 17 & 3-2-4-4-1-1-1-1 \\ \hline \end{tabular} } \begin{minipage}{14cm} \caption{Results obtained by {\tt NJGRAF} and {\tt NJFORMULA} for a range of recoupling coefficients.} \label{tab-njgf} \end{minipage} \end{center} \end{table} From these examples it is clear that our program {\tt NJFORMULA} performs very well. The results are at least as good as the results of the program {\tt NJGRAF} and are obtained in a comparable running time. In many cases we obtain a better formula than {\tt NJGRAF}, and {\tt NJFORMULA} can also handle very complicated cases where {\tt NJGRAF} crashes. It is a small program (the executable code is about 20KB, while the compiled code of {\tt NJGRAF} is 260KB), which runs very well on various systems, also on small systems such as PCs. % % ACKNOWLEDGEMENTS % %\pagebreak \section*{Acknowledgements} It is a pleasure to thank Prof.\ K.\ Srinivasa Rao (IMSC, Madras, India) for stimulating discussions. This work was partly supported by the EEC (contract no.\ CI$1^*$--CT92--0101). % % BIBLIOGRAPHY % \begin{thebibliography}{99} \bibitem{Bur} P.G.\ Burke, Comput.\ Phys.\ Commun.\ {\bf 1} (1970) 241. \bibitem{Bar} A.\ Bar-Shalom and M.\ Klapisch, Comput.\ Phys.\ Commun.\ {\bf 50} (1988) 375. \bibitem{FPV} V.\ Fack, S.N.\ Pitre and J.\ Van der Jeugt, {\em New efficient programs to calculate general recoupling coefficients, part II (in preparation)} \bibitem{Bie} L.C.\ Biedenharn and J.D.\ Louck, The Racah-Wigner algebra in Quantum Theory, Encyclopedia of Mathematics and its Applications, Vol.~9, ed.\ G.-C.\ Rota (Addison Wesley, Reading, MA, 1981). \bibitem{Edm} A.\ Edmonds, Angular momentum in quantum mechanics (Princeton University Press, Princeton, 1957). \bibitem{Yut} A.P.\ Yutsis, I.B.\ Levinson and V.V.\ Vanagas, The Theory of Angular Momentum (Israel Program for Scientific Translation, Jerusalem, 1962). \bibitem{Bri} D.M.\ Brink and G.R.\ Satchler, Theory of Angular Momentum (Clarendon Press, Oxford, 1968). \bibitem{Jud} B.R.\ Judd, Operator Techniques in Atomic Spectroscopy (McGraw-Hill, New York, 1963). \bibitem{Elb} E.\ El Baz and B.\ Castel, Graphical Methods of Spin Algebras (Marcel Dekker, New York, 1972). \bibitem{DSA} E.\ Horowitz and S.\ Sahni, Fundamentals of Data Structures in Pascal (Pitman, London, 1984). \bibitem{Linux} Linux (Unix clone for 386/486-based PCs) version 1.0; publicly available via anonymous ftp to {\tt nic.funet.fi} in directory {\tt /pub/OS/Linux}. \bibitem{TC} Turbo C++ (version 1.01) User Guide, Borland International Inc.\ (1990). \bibitem{GnuC} GNU CC version 2.5, Free Software Foundation, Cambridge, MA, USA; publicly available via anonymous ftp to {\tt prep.ai.mit.edu} in directory {\tt /pub/gnu}. \bibitem{ParC} Parallel C (version 2.2.2) User Guide, 3L Ltd.\, UK (1991). \bibitem{ParF} Parallel Fortran (version 2.1.3) User Guide, 3L Ltd.\, UK (1990). \end{thebibliography} \end{document}