Computation of the Multi-Spheres Scattering Coefficient Using the Prime Index Method

Abstract

The analytical-numerical evaluation of the scattering of electromagnetic waves by multiple spheres requires the computation of numerous coefficients. For this purpose, many contributions, available in the literature, have traditionally employed the recursion method. In the present paper, we introduce a novel approach, based on primes and indices, which can be conveniently applied to the computation of the Wigner 3-j symbols, the Wigner D-function, and the Gaunt coefficients. By considering a series-expansion form, our method proves to be easily applicable to a variety of similar problems. We provide examples of coefficient calculations and compare the results with those retrieved from previous publications, demonstrating the advantages of our approach.

1. Introduction

Factorization is an age-old problem arising in various fields of pure and applied mathematics. Due to computational limitations, the evaluation of $n!$, when n is very large, remains a challenging task. In some cases, the complete calculation of the factorial function is not necessary, as illustrated by the evaluation of the quantity ${\displaystyle \frac{500!}{499!}}$. As we know, the value of the ratio is 500. However, any attempt to calculate it directly using factorial functions, available in many programming platforms such as Matlab, is likely not to produce any result. This type of computational problem is akin to what needs to be addressed in multi-sphere scattering. As demonstrated in the following sections, the solutions of scattering problems according to Mie theory are based on numerous coefficients, such as the Wigner 3-j symbols, the Wigner D-functions, and the Gaunt coefficients. These quantities involve the evaluation of ratios of high-valued factorial of natural numbers, numerically hindered by limitations in the finite representation of integers and possible loss of precision. Although several computational strategies for the evaluaton of these quantities have been proposed in the literature [1,2,3,4,5,6,7,8,9,10], the high computational complexity relevant to the retrievement of these quantities is still an issue deserving further investigation. The contribution presented in [1], employing graphical techniques to the computation of Wigner 3-j symbols, suffers from a limited extensibility to the computation of other quantities involved in the solution of multi-sphere scattering problems. Well-established methods for the evaluation of the Wigner 3-j symbols and the Gaunt coefficients, available in the relevant literature, are mainly based on the exploitation of recursive computational schemes [5,6,10], optimized by the identification of the numerous symmetry properties of the quantities to evaluate [6], or are focused on the achievement of high numerical speed by devising efficient storage schemes to allow for fast memorization and retrieval of the large number of values typically required by numerical simulations [3]. While analythical approaches are useful for the definition of asymptotic forms in the case of high-degree coefficients [2], floating-point arithmetic routines must find a trade-off between numerical precision and speed. The well-documented computational strategy devised by Xu [7,8], based on the solution of a linear system of equations, although appreciable under many aspects, is not devoid of some level of inaccuracy and presents quite a high level of complexity, failing to provide easy implementation in a numerical code. In this paper, we present a computational strategy for the evaluation of the aforementioned coefficients, named the Prime Index Method (PIM). The proposed approach leverages on the prime factorization of integers’ factorials in index form, simplifying the computation of ratios of large factorials of natural numbers and allowing for the solution of multi-sphere scattering problems with very high precision and speed.

This paper is organized as follows: in Section 1, the Prime Index Method (PIM) is explained and its main idea is illustrated by means of a simple example. In Section 2, we provide examples of relevant applications for the proposed method and compare the outcomes to the results available in the literature. Conclusions and further developments are discussed in Section 3.

2. The Prime Index Method

The proposed Prime Index Method (PIM) is based on Legendre’s theorem, as described in [11], page 8. The theorem states that the factorial contains the prime factor p exactly $\kappa$ times, the value $\kappa$ being expressed by the following summation:

$$\kappa =\sum _{i=1}^{\infty }⌊{\displaystyle \frac{n}{p^i}}⌋$$

(1)

where the symbol $\lfloor \rfloor$ denotes the floor function. Thus, the quantity $\kappa$ provides the exponent at which the prime factor p must be raised the prime number factorization of $n!$. The infinite summation in Equation (1) may be conveniently truncated, as the quantity $⌊{\displaystyle \frac{n}{p^i}}⌋$ is null if $p^i>n$. By applying Legendre’s theorem, it is possible to achieve the prime-number factorization of $n!$. All prime numbers and their respective exponents of this factorization can be stored in two vectors $p_n$${\kappa }_n$, repectively. Each element $p_n\left(j\right)$ of vector $p_n$ represents a prime factor less than or equal to n and each element ${\kappa }_n\left(j\right)$ represents the exponent of the prime $p_n\left(j\right)$ in the factorization of $n!$. By multiplying all factors $p_n{\left(j\right)}^{{\kappa }_n\left(j\right)}$ toghether, the quantity $n!$ can be retrieved.

For this purpose, we define the vector $p_n$ of length $j_{max}\left(n\right)$, where $p_n\left(j\right)$ ($j=1\dots j_{max}\left(n\right)$) contains all and only the prime numbers greater than 1 and less than or equal to n. The symbol $p_n$ has been chosen to state that the dimension of the vector $p_n\left(j\right)$ is a function of n. The vector $p_n\left(j\right)$ may be easily generated by using the sieve algorithm or library functions available in many commercial software platforms, e.g., the primes function available in Matlab. The value $j_{max}\left(n\right)$ represents the number of primes less than or equal to n: its value (the well-known prime counting function $\pi \left(n\right)$[12]) can be easily retrieved by evaluating the dimension of vector $p_n$.

Secondly, we define an index vector ${\kappa }_n\left(j\right)$ ($j=1\dots j_{max}\left(n\right)$), using Equation (1)

$${\kappa }_n\left(j\right)=\sum _{i=1}^{\infty }⌊{\displaystyle \frac{n}{p_n{\left(j\right)}^i}}⌋.$$

(2)

As stated before, in the numerical implementation of the algorithm, the infinite summation in Equation (2) can be truncated when $i>i_{max}\left(n\right)$, where $i_{max}\left(n\right)$ is the first value for index i, for which the addend reaches zero. Having defined the index vector ${\kappa }_n\left(j\right)$, the factorization of $n!$ can be obtained by raising all $j_{max}\left(n\right)$ prime fators $p_n\left(j\right)$ to their respective exponents ${\kappa }_n\left(j\right)$ and multipying all $p_n{\left(j\right)}^{{\kappa }_n\left(j\right)}$ factors together, i.e.,

$$n!=\prod _{j=1}^{j_{max}\left(n\right)}{p}_n{\left(j\right)}^{{\kappa }_n\left(j\right)}.$$

(3)

Thus, the value of $n!$ can be conveniently expressed by means of two vectors of integers, i.e., $p_n$ and ${\kappa }_n$. Equation (3) represents a more convenient expression of $n!$, which will be used in the following sections. A careful inquiry would confirm that the two expressions of $n!$ do obtain the same results.

Consider the calculation of the following ratio:

$${\displaystyle \frac{n!}{m!}}={\displaystyle \frac{{\prod }_{j=1}^{j_{max}\left(n\right)}{p}_n{\left(j\right)}^{{\kappa }_n\left(j\right)}}{{\prod }_{j^{\prime }=1}^{j_{max}\left(m\right)}{p}_m{\left(j^{\prime }\right)}^{{\kappa }_m\left(j^{\prime }\right)}}}$$

(4)

The fraction in Equation (4) can be reduced to its lowest terms by simplifying the common factors present in both the numerator and denominator. To illustrate the simplification algorithm for Equation (4), we suppose that $n<m$: given the prime vectors $p_n$ and $p_m$, their respective index vectors ${\kappa }_n$ and ${\kappa }_m$ may be set by means of Equation (2). Supposing $n<m$, it is also true that $j_{max}\left(n\right)\le j_{max}\left(m\right)$, i.e., the prime factorization of $m!$ contains $j_{max}\left(m\right)-j_{max}\left(n\right)$ additional prime factors, which are not included in the prime factorization of $n!$. Nevertheless, it is possible to express the prime factorization of $n!$ by means of the prime vector $p_m$, provided that all the $j_{max}\left(m\right)-j_{max}\left(n\right)$ additional prime factors contained in the prime factorization of $m!$ (but not included in the prime factorization of $n!$) are all raised to a null exponent. As a matter of fact, the introduction of further prime factors does not affect the factorization of $n!$, provided all the exponents of the introduced factors are set to zero. In order to acheive this:

• the dimension of the index vector ${\kappa }_n$ must be increased to match the dimension of vector ${\kappa }_m$ by inserting $j_{max}\left(m\right)-j_{max}\left(n\right)$ elements.
• the new elements appended to vector ${\kappa }_n$ are all null.

The modified version of index vector ${\kappa }_n$ will be denoted by ${\kappa }_n^{\prime }$. The new expression for the prime factorization of $n!$ is:

$$n!=\prod _{j=1}^{j_{max}\left(m\right)}{p}_m{\left(j\right)}^{{\kappa }_n^{\prime }\left(j\right)}$$

(5)

By inserting Equation (5) into Equation (4), we obtain:

$${\displaystyle \frac{n!}{m!}}={\displaystyle \frac{{\prod }_{j=1}^{j_{max}\left(m\right)}{p}_m{\left(j\right)}^{{\kappa }_n^{\prime }\left(j\right)}}{{\prod }_{j^{\prime }=1}^{j_{max}\left(m\right)}{p}_m{\left(j^{\prime }\right)}^{{\kappa }_m\left(j^{\prime }\right)}}}=\prod _{j=1}^{j_{max}\left(m\right)}{\displaystyle \frac{p_m{\left(j\right)}^{{\kappa }_n^{\prime }\left(j\right)}}{p_m{\left(j\right)}^{{\kappa }_m\left(j\right)}}}=\prod _{j=1}^{j_{max}\left(m\right)}{p}_m{\left(j\right)}^{\left[{\kappa }_n^{\prime }\left(j\right)-{\kappa }_m\left(j\right)\right]}.$$

(6)

Similar considerations apply if $n>m$. The computational strategy expressed by Equation (6) is applicable even when the value of the factorial functions exceed the computational limits allowed by the computer. We will illustrate the advantages introduced by Equation (6) by means of a simple example: the computation of ${\displaystyle \frac{4!}{5!}}$. In this case, the maximum value between the two arguments of the factorial functions is 5, and the vector $p_m$, listing the the primes under 5, is $[2,3,5]$. According to Equation (2), the index vector for $5!$ is ${\kappa }_m=[3,1,1]$, while the index vector for $4!$ is ${\kappa }_n^{\prime }=[3,1,0]$—the last element, 0, has been added to the elements of ${\kappa }_n=[3,1]$. Therefore, the index difference vector ${\kappa }_n^{\prime }\left(j\right)-{\kappa }_m\left(j\right)$ in Equation (6) is $[0,0,-1]$, thus:

$${\displaystyle \frac{4!}{5!}}=2^0\times 3^0\times 5^{-1}={\displaystyle \frac{1}{5}}.$$

(7)

This is the expected result. The previous example is purely illustrative, as the presented method does not show any advantages when both n and m have small values.

The convenience of the presented method is apparent if n and m are large: the computation of the ratio ${\displaystyle \frac{500!}{499!}}=500$ by prime factorization is straightforward, while any attempt of direct computation retrieves a $NaN$ result for the factorial function in Matlab. Moreover, the computational load in the calculation of ${\displaystyle \frac{500!}{499!}}$ by PIM is negligible, being of the order of milliseconds using an average commercial laptop computer.

It is possible to address more complex situations, e.g., the computation of expressions like ${\displaystyle \frac{n!}{\sqrt{m!}}}$, or, in general, expressions in which the numerator or denominator of a fraction consist of the product of several factorial terms, each one raised to a specific rational exponent, e.g.,

$${\displaystyle \frac{{(n_1!)}^{q_1}{(n_2!)}^{q_2}\dots }{{(m_1!)}^{r_1}{(m_2!)}^{r_2}\dots }},q_i,r_j\in \mathbb{Q}.$$

(8)

The presence of additional factorial factors and their rational exponents may be taken into account, in the prime index representation of Equation (6), by the following algorithm:

• we compute the index vectors ${\kappa }_{n_1}$, ${\kappa }_{m_1}$, ${\kappa }_{n_2}$, and ${\kappa }_{m_2}$…, corresponding to the arguments of the factorial functions $n_1$, $m_1$, $n_2$, and $m_2$ …, respectively;
• supposing that $\gamma =max\{n_1,m_1,n_2,m_2\dots \}$, we increase the dimension of vectors ${\kappa }_{n_i},{\kappa }_{m_i}$ ($n_i,m_i\ne \gamma$) in order to match the dimension of vector ${\kappa }_{\gamma }$, obtaining vectors ${\kappa }_{n_1}^{\prime }$, ${\kappa }_{m_1}^{\prime }$, ${\kappa }_{n_2}^{\prime }$, ${\kappa }_{m_2}^{\prime }$…, ${\kappa }_{\gamma }$,…, as previously described;
• we define a matrix $\kappa (i,j)$, whose rows are formed by the vectors ${\kappa }_{n_1}^{\prime }$, ${\kappa }_{m_1}^{\prime }$, ${\kappa }_{n_2}^{\prime }$, ${\kappa }_{m_2}^{\prime }$…, ${\kappa }_{\gamma }$,…;
• we define a vector $\xi$ whose elements are the exponents $q_1$, $q_2$, …, , $-r_1$, $-r_2$…, thus $\xi =[q_1,,q_2,\dots ,-r_1,-r_2\dots ]$.

Equation (6) becomes:

$${\displaystyle \frac{{(n_1!)}^{q_1}{(n_2!)}^{q_2}\dots }{{(m_1!)}^{r_1}{(m_2!)}^{r_2}\dots }},=\prod _{j=1}^{j_{max}\left(\gamma \right)}{p}_{\gamma }{\left(j\right)}^{{\sum }_i\xi \left(i\right)\kappa (i,j)}.$$

(9)

For convenience, in the numerical implementation of the algorithm, arguments $n_i$, $m_i$ of the factorial functions in Equation (9) may be represented as components of a vector x, i.e., $x=[n_1,n_2,\dots ,m_1,m_2,\dots ]$, as previously done for the relevant exponents represented as the elements of vector $\xi$.

A Matlab Code for the Evaluation of the quantity in Equation (9) is reported in the Supplementary Materials, provided along with the present paper. In the following sections, we will consider relevant application of Equation (9).

3. Results and Discussion

In this section, Equation (9) is applied to the evaluation of coefficients involved in the numerical solution of multi-sphere scattering problems. The accuracy of the computational algorithms is assessed and a comparison with previous results, available in the literature, is provided.

3.1. The Wigner 3-j Symbols

The expression of the Wigner 3-j symbols, as available in [3] is:

$$\begin{array}{c}\left(\begin{array}{ccc}{j}_1& j_2& j_3\\ m_1& m_2& m_3\end{array}\right)=ϵ\left(j_1,j_2,j_3\right)\Delta \left(j_1,j_2,j_3\right){\delta }_{m_1+m_2+m_3,0}{(-1)}^{j_1-j_2-m_3}\times \hfill \\ \hfill \sqrt{\left(j_1+m_1\right)!\left(j_1-m_1\right)!\left(j_2+m_2\right)!\left(j_2-m_2\right)!\left(j_3+m_3\right)!\left(j_3-m_3\right)!}\times \hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\sum _{k=k_{min}}^{k_{max}}{\displaystyle \frac{{(-1)}^k}{k!\left(j_1+j_2-j_3-k\right)!\left(j_1-m_1-k\right)!\left(j_2+m_2-k\right)!}}\times {\displaystyle \frac{1}{\left(j_3-j_2+m_1+k\right)!\left(j_3-j_1-m_2+k\right)!}}.\end{array}$$

(10)

where $k_{max}=min\left(j_1+j_2-j_3,j_1-m_1,j_2+m_2\right)$ and $k_{min}=max\left(-j_3+j_2-m_1,-j_3+j_1\right.$$\left.+m2,0\right)$. The quantities $m_1$, $m_2$, and $m_3$ represent an angular momentum quantum numbers: to avoid ambiguity, we refer to the arguments $m_i$ of the factorial functions in Equation (9) as components of vector x.

The triangle coefficient is defined as follows:

$$\Delta \left(j_1,j_2,j_3\right)=\sqrt{{\displaystyle \frac{\left(j_1+j_2-j_3\right)!\left(j_1-j_2+j_3\right)!\left(-j_1+j_2+j_3\right)!}{\left(j_1+j_2+j_3+1\right)!}}},$$

(11)

$$ϵ\left(j_1,j_2,j_3\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\left(j_1,j_2,j_3\right)\mathrm{form}\mathrm{a}\mathrm{triangle}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.,$$

(12)

$${\delta }_{m_1+m_2+m_3,0}=\left\{\begin{array}{cc}1\hfill & \mathrm{when}{m}_1+m_2+m_3=0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right..$$

(13)

Although the expression may seem rather complex, we will focus on the numerical evaluation of the ratio of factorial functions. We define the x-vector components as follows:

$$x\left(1\right)=j_1+j_2-j_3,x\left(2\right)=j_1-j_2+j_3,x\left(3\right)=-j_1+j_2+j_3,$$

$$x\left(4\right)=j_1+m_1,x\left(5\right)=j_1-m_1,x\left(6\right)=j_2+m_2,x\left(7\right)=j_2-m_2,$$

$$x\left(8\right)=j_3+m_3,x\left(9\right)=j_3-m_3$$

$$x\left(10\right)=j_1+j_2+j_3+1,x\left(11\right)=k,x\left(12\right)=j_1+j_2-j_3-k,x\left(13\right)=j_1-m_1-k,$$

$$x\left(14\right)=j_2+m_2-k,x\left(15\right)=j_3-j_2+m_1+k,x\left(16\right)=j_3-j_1+m_1+k.$$

The correspondent $\xi$ vector is as follows:

$$\xi =[1/2,1/2,1/2,1/2,1/2,1/2,1/2,1/2,1/2,-1/2,-1,-1,-1,-1,-1,-1].$$

The numerical results can be compared to those presented in a previous work by Xu [7]. In [7], the values of the Wigner 3-j symbols are obtained by applying an iterative method, which is much more complex than the one presented here.

To provide an example of the convenience of the proposed algorithm, we consider the following computation:

$$\left(\begin{array}{ccc}260& 280& j_3\\ 228& 268& -496\end{array}\right),$$

(14)

where $j_3$ ranges from 496 to 540, and the results are presented in Figure 1. In Figure 1a, values of the Wigner 3-j symbol in (14), retrieved by means of the PIM (solid line), are plotted versus values of $j_3$ and compared with those retrieved by Xu’s method [7] (circle markers): the two sets of plotted values match with a high precision. Figure 1b reports values of the maximum relative error between the results reported in Figure 1a, retrieved by the two methods, plotted versus values of $j_3$: the relative error keeps very low in the interval of values of $j_3$ under analysis. A Matlab Code for the Evaluation Wigner 3-j symbols is reported in the Supplementary Materials, provided along with the present paper. Regarding the computational speed, Xu’s method is still hindered by the necessity to retrieve the factorial of large numbers. The results of Wigner3-j coefficients, presented in his papers, are retrieved by implementing a very complex algorithm in Fortran code running on a workstation. Xu points out that the computational time, when $j_{max}\left(n\right)$ approaches the value of 80, results in approximately 30 h. When applying the PIM, the computational load is extremely reduced.The typical evaluation times do not exceed a tenth of a second on a commercial PC running interpreted-language Matlab code, even for orders $j_{max}\left(n\right)$ much larger than 80. Such an improvement in numerical performance can be reasonably claimed, even taking into consideration the technological evolution of computers in 25 years.

3.2. Winger D-Function

The expression of the Wigner D-function, available in [13], is as follows:

$$\begin{array}{cc}\hfill d_{mn}^S\left(\theta \right)=& \sqrt{(s+m)!(s-m)!(s+n)!(s-n)!}\hfill \\ \hfill & \times \sum _{k=k_{min}}^{k_{max}}{(-1)}^k{\displaystyle \frac{{\left(cos\frac{1}{2}\theta \right)}^{2s-2k+m-n}{\left(sin\frac{1}{2}\theta \right)}^{2k-m+n}}{k!(s+m-k)!(s-n-k)!(n-m+k)!}}\hfill \end{array}$$

(15)

where $k_{min}=max(0,m-n)$, and $k_{max}=min(s+m,s-n)$.

To apply the proposed algorithm, we define a x vector, as previously obtained:

$$x\left(1\right)=s+m,x\left(2\right)=s-m,x\left(3\right)=s+n,x\left(4\right)=s-n,$$

$$x\left(5\right)=k,x\left(6\right)=s+m-k,x\left(7\right)=s-n-k,x\left(8\right)=n-m+k,$$

and define a correspondent $\xi$ vector:

$$\xi =[1/2,1/2,1/2,1/2,-1,-1,-1,-1].$$

The results of PIM can be compared with the analytical results from [14].

$$d_{0,0}^4\left(\theta \right)={\displaystyle \frac{1}{8}}\left(3-30{cos}^2\theta +35{cos}^4\theta \right)$$

(16)

$$d_{2,2}^4\left(\theta \right)={\displaystyle \frac{1}{4}}{(1+cos\theta )}^2\left(1-7cos\theta +7{cos}^2\theta \right)$$

(17)

$$d_{2,-2}^4\left(\theta \right)={\displaystyle \frac{1}{4}}{(1-cos\theta )}^2\left(1+7cos\theta +7{cos}^2\theta \right)$$

(18)

The results are presented in Figure 2. In Figure 2a, the values of Wigner $d_{0,0}^4$, $d_{2,2}^4$, and $d_{2,-2}^4$ functions retrieved by numerical computations (solid line) are plotted versus values of $\theta$ and compared with the analytical results (circle markers). Figure 2b reports the values of the maximum relative error between the results reported in Figure 2a, retreived by the two methods, plotted versus the values of $\theta$: the relative error remains very low in the interval of values of $\theta$ under analysis.

3.3. Gaunt Coefficient

The expression of the Gaunt coefficient, available in [8], is as follows:

$$\begin{array}{cc}\hfill a(m,n,\mu ,\nu ,p)=& {(-1)}^{m+\mu }(2p+1){\left[{\displaystyle \frac{(n+m)!(\nu +\mu )!(p-m-\mu )!}{(n-m)!(\nu -\mu )!(p+m+\mu )!}}\right]}^{1/2}\hfill \\ \hfill & \times \left(\genfrac{}{}{0ex}{}{n \nu p}{0 0 0}\right)\left(\genfrac{}{}{0ex}{}{n \nu p}{m \mu -m-\mu }\right),\hfill \end{array}$$

(19)

In Equation (19), the evaluation of the two Wigner 3-j symbols and the factor in square brackets may be performed according to the procedure explained previously. The results are presented in Figure 3. In Figure 3a, the values of the Gaunt coefficient retrieved by the PIM are plotted versus the values of $n=\nu$ and compared with those retrieved by Xu’s method. The values of other parameters are $m=1$, $\mu =-1$, and $p=n+\nu$. In Figure 3b, the values of the relative error between the results reported in Figure 3a, retrieved by the two methods, are plotted versus values of $n=\nu$: the relative error remains very low for the interval values of $n=\nu$ under analysis.

3.4. Vector Translation Coefficients

The expression of the Vector Translation Coefficients, available in [7], is as follows:

$$\begin{array}{cc}\hfill A_{-mn\mu \nu }^{lj}=& {(-1)}^m{\displaystyle \frac{(2\nu +1)(n-m)!(\nu -\mu )!}{2n(n+1)(n+m)!(\nu +\mu )!}}exp\left[i(\mu +m){\varphi }_{lj}\right]\hfill \\ & \times \sum _{q=0}^{q_{max}}{i}^p[n(n+1)+\nu (\nu +1)-p(p+1)]a_q\hfill \\ & \times \left[\genfrac{}{}{0ex}{}{h_p^{\left(1\right)}\left(k{d}_{lj}\right)}{j_p\left(k{d}_{lj}\right)}\right]P_p^{\mu +m}\left(cos{\theta }_{lj}\right)\left(\genfrac{}{}{0ex}{}{r\le d_{lj}}{r>d_{lj}}\right),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill B_{-mn\mu \nu }^{lj}=& {(-1)}^{m+1}{\displaystyle \frac{(2\nu +1)(n-m)!(\nu -\mu )!}{2n(n+1)(n+m)!(\nu +\mu )!}}exp\left[i(\mu +m){\varphi }_{lj}\right]\hfill \\ & \times \sum _{q=1}^{Q_{max}}{i}^{p+1}{\left\{\left[{(p+1)}^2-{(n-\nu )}^2\right]\left[{(n+\nu +1)}^2-{(p+1)}^2\right]\right\}}^{1/2}\hfill \\ & \times b(m,n,\mu ,\nu ,p+1,p)\left[\begin{array}{c}{h}_{p+1}^{\left(1\right)}\left(k{d}_{lj}\right)\\ j_{p+1}\left(k{d}_{lj}\right)\end{array}\right]P_{p+1}^{\mu +m}\left(cos{\theta }_{lj}\right)\left(\begin{array}{c}r\le d_{lj}\hfill \\ r>d_{lj}\hfill \end{array}\right),\hfill \end{array}$$

(21)

where

$$\begin{array}{cc}\hfill Q_{max}=& min[n,\nu ,(n+\nu +1-|m+\mu |)/2],\hfill \\ \hfill b(m,n,\mu ,\nu ,p+1,p)=& {(-1)}^{\mu +m}(2p+3){\left[{\displaystyle \frac{(n+m)!(\nu +\mu )!(p-m-\mu +1)!}{(n-m)!(\nu -\mu )!(p+m+\mu +1)!}}\right]}^{1/2}\hfill \\ & \times \left(\begin{array}{ccc}n& \nu & p+1\\ m& \mu & -m-\mu \end{array}\right)\left(\begin{array}{ccc}n\hfill & \nu \hfill & p\hfill \\ 0\hfill & 0\hfill & 0\hfill \end{array}\right).\hfill \end{array}$$

(22)

In Equations (20) and (21), k is the wavenumber and $d_{ij},{\theta }_{ij},{\varphi }_{ij}$ are the spherical coordinates of the origin of the $j^{th}$ coordinate system in the $l^{th}$ coordinate system. Note that Equation (22) involves the use of the Gaunt coefficient. A comparison can be carried out between Xu’s results, presented in [7], where $(k{d}_{ij},{\theta }_{ij},{\varphi }_{ij})=(2,0.5,0.5)$, results are reported in

Table 1. Comparison between the values of the Vector Translation Coefficients retreived by Xu’s method [7] and the PIM for $(k{d}_{ij},{\theta }_{ij},{\varphi }_{ij})=(2,0.5,0.5)$ and different values of m, n, $\mu$, and $\nu$ parameters.

				${\mathit{A}}_{\mathit{m}\mathit{n}\mathit{\mu }\mathit{\nu }}^{\mathit{l}\mathit{j}}({\mathit{kd}}_{\mathit{ij}},{\mathit{\theta }}_{\mathit{ij}},{\mathit{\varphi }}_{\mathit{ij}})$
				XU		PIM
m	n	$\mu$	$\nu$	Real	Imag	Real	Imag
8	10	−9	12	3.663964990 × 1034	−2.762412192 × 1034	3.663965099 × 1034	−2.762412274 × 1034
0	10	0	10	2.969682019 × 10−1	−1.928601440 × 1017	2.969682108 × 10−1	−1.928601498 × 1017
−2	11	3	9	7.726121583 × 1011	1.034255820 × 1012	7.726121813 × 1011	1.034255850 × 1012
10	18	15	22	−6.206840651 × 1035	−8.308775621 × 1035	−6.206840836 × 1035	−8.308775868 × 1035
36	36	−38	38	4.146334728 × 10190	−4.931584782 × 10190	4.146334852 × 10190	−4.931584929 × 10190
				${\mathit{B}}_{\mathit{m}\mathit{n}\mathit{\mu }\mathit{\nu }}^{\mathit{l}\mathit{j}}({\mathit{kd}}_{\mathit{ij}},{\mathit{\theta }}_{\mathit{ij}},{\mathit{\varphi }}_{\mathit{ij}})$
				XU		PIM
m	n	$\mu$	$\nu$	Real	Imag	Real	Imag
8	10	−9	12	−8.370892023 × 1031	−1.110285257 × 1031	−8.370892272 × 1031	−1.110285290 × 1032
0	10	0	10	0.000000000 × 100	0.000000000 × 100	−7.341232630 × 10−4	−1.492552121 × 10−18
−2	11	3	9	1.222239141 × 1010	9.130398908 × 109	1.222239177 × 1010	−9.130399180 × 109
10	18	15	22	−3.610252125 × 1034	2.696938836 × 1034	−3.610252233 × 1034	2.696938917 × 1034
36	36	−38	38	0.000000000 × 100	0.000000000 × 100	0.000000000 × 100	0.000000000 × 100

. In particular, when $m=-\mu$ and $n=\nu$, values $B_{mn\mu \nu }^{lj}(k{d}_{ij},{\theta }_{ij},{\varphi }_{ij})$, obtained by our method, are different from Xu’s. In [9], Xu reported a difference between the results of his method and the results retrieved by other methods, without indicating which results are correct: yet, the author points out that his method evidently failed to retrieve the value of coefficient B for some specific values of the parameters $(m,n,\mu ,\nu )$. Computing the Vector Translation Coefficients B, in cases in which Xu’s method cannot provide the correct answer, e.g., when $(m,n,\mu ,\nu )=(0,10,0,10)$, PIM still retrieves reasonably reliable values, although we could not perform a comparison with the previous results, owing to the small amount of publications on this specific topic. It is worth noting that the incorrect results for A and B, reported in the literature, were due to errors in the calculation methods of the Wigner 3-j symbol. Although the values of the parameters under consideration may be relatively small, the accumulation of small errors may jeopardize the reliability and stability of the numerical solution of complex scattering problems, typically relying on the computation of a very large amount of such coefficient values.

4. Conclusions

In this paper, we present a novel approach for the solution of many numerically challenging computational problems in which the presence of high-valued factorial functions, in the numerator and denominator of an expression, exceed any reasonable computational limit, resulting in a computer-memory overflow or the loss of numerical precision. The proposed PIM approach, based on the evaluation of the series expansions, is applicable to the calculation of ratios of factorial terms, each one raised to a specific rational exponent and, which is proven to be extremely useful for the computation of the Wigner 3-j symbols, the Wigner D-function, and the Gaunt coefficient. The method is fast, reliable, and easily implementable in code, with the key code requiring only around 30 lines. Computational times are extremely reduced even in the evaluation of more complex functions, such as Vector Translation Coefficients $\mathbf{A}$ and $\mathbf{B}$. Advantages are apparent in many other cases, such as the evaluation of the spherical Hankel functions, finding its straightforward application in the numerical solution of electromagnetic scattering problems. In this case, the series-expansion truncation criterion imposes a heavy computational toll when the length of the incident wave is much smaller than the size of the scatterer. By applying the proposed method, computational limitations may be easily overcome, achieving a very high numerical precision [15].