To simplify the notations for summing equation (12), we can rewrite it using matrices and arrays,

	(42)
where, for the 2nd order R-A approximation (6x6 matrix size)
	(43)

The relationship between the element subscripts and Rehr-Albers expansion indices [l' = (m' n'), l = (m n)] is listed in Tables 5 - 8. Because the F elements converge for large r roughly as (r') ^-(|m'|+2n') (r)^{- (|m|+2n)}, we can easily see, in this notation, there are only one element in each above array or matrix (, G₁ and F₁₁) for zeroth Rehr-Albers expansion order, equivalent to the asymptotic high-energy limit and to the point-scattering approximation. We only need to keep the first three elements in the and G arrays and a (3x3) matrix in F up to first Rehr-Albers expansion order. For the full calculations, Rehr and Albers have demonstrated that in practice it seems sufficient to retain terms only to second order, i.e. the (6x6) matrix, in most cases.

l' = (m', n')	(0,0)	(1,0)	(-1,0)	(0,1)	(2,0)	(-2,0)
Ln=(l,m)	(0,0)	(0,0)	(0,0)	(0,0)	(0,0)	(0,0)

Table 5. The relationship between the subscripts of the factor elements and the Rehr-Albers expansion indices.

From equation (42), we see that the number of array or matrix multiplications dominates the total calculation time. For the (n-1)-th scattering order term, the number of terms in the summation is proportional to N^n-2 times the number of emitters, and times the number of L_f=(l_f,m_f) combinations, where N is the number of atoms in question. This procedure takes much time, even when using pathcut.

We have therefore introduced an iterative summing method, which makes the calculation time be proportional to N³. For the commonly used 6th to 8th scattering order, such a method can cut down the computation time enormously. To that end, we rewrite the equation (42) into a forward summing version as

	(44)
or in a backward summing version as
	(45)

The calculation time will be proportional to N_k·N_mi·(msorder-2)· N_a³ for forward summing and N_k·N_d·(msorder-2)· N_a³ for backward summing, where N_k is the number of photoelectron energies to be simulated, N_mi = (2l_i+1) the number of possible magnetic quantum numbers, N_d the number of photoelectron exit angles, msorder the multiple scattering order, and N_a the number of atoms in the cluster.

The most time-saving contribution which Rehr and Albers made is the derivation of the separable representation of the scattering event, i.e. equation (26). Let's rewrite it as

(46)

Thus, the energy-dependent scattering factor f_ll' (r,r',b) now depends only on the bond lengths and the scattering angle b (which is equivalent to arccos(r,r')),

(47)

This form considerably reduces the number of independent computations of the scattering matrices. Furthermore the only additional angle needed for each element in the calculation is the combinations of a and g at each scattering site.

The most time-consuming computation in each path is that of equation (47). There are N³ choices of r and r', where N is the number of atoms in the cluster. That means we have to compute N³ scattering events f_ll' (r,r',b). However, in practice, in a crystalline surface with its periodic and other symmetries, many events f_ll' have the same r, r' and b values (the absolute orientation of the vectors r, r' does not matter). Therefore, there will be many identical results for f_ll' (r,r',b). So, we can investigate all the scattering events, and calculate only once the identical events. In our program, we further pre-calculate for a preselected series of r and b values, respectively, then use interpolation to produce their values in the computation of f_ll' (r,r',b).