Master-Slave Algorithm for Highly Accurate and Rapid Computation of the Voigt / Complex Error Function

We obtain a rational approximation of the Voigt / complex error function by Fourier expansion of the exponential function e − ( t − 2 σ ) 2 and present master-slave algorithm for its e ﬃ cient computation. The error analysis shows that at y > 10 − 5 the computed values match with highly accurate references up to the last decimal digits. The common problem that occurs at y → 0 is e ﬀ ectively resolved by main and supplementary approximations running computation ﬂow in a master-slave mode. Since the proposed approximation is rational function, it can be implemented in a rapid algorithm

At positive y the complex probability function W (z) is reduced to (Faddeyeva & Terent'ev, 1961;Armstrong & Nicholls, 1972) where is the complex error function, also known as the Faddeeva function.Consequently, we can write w (x, y) = K (x, y) + iL (x, y) , y > 0.
The complex error function is a solution of the differential equation with initial condition (Schreier, 1992): The real and imaginary parts of the complex error function satisfy the linear system of the first order partial differential equations: where in accordance with Cauchy-Riemann equations: ∂K(x,y) ∂y = ∂L(x,y)   ∂x ∂K(x,y) ∂x = − ∂L(x,y) ∂y .
It is convenient to rewrite Equation (4) in terms of the error and complementary error functions: or erfc (z) ≡ 1 − erf (z) = e −z 2 w (iz) .
The complex error function has very broad practical applications.In particular, the real part of the complex error function K (x, y) = Re w (x, y) , known as the Voigt function (Faddeyeva & Terent'ev, 1961;Armstrong & Nicholls, 1972;Schreier, 1992;Letchworth, 2007), is widely used in many fields of Physics, Astronomy and Chemistry.Mathematically, the Voigt function represents a convolution of Cauchy (or Lorentzian) and Gaussian (or Doppler) distributions and can be used to describe the spectral line broadening effects associated with photon emission and/or absorption, for example, in atmospheric molecules (Berk, 2013;Quine & Abrarov, 2013;Schreier et al., 2014), photo-luminescent materials (Miyauchi et al., 2013), celestial bodies like stars or planets (Emerson, 1996) and so on.The Voigt function can accurately describe the spectral characteristics of a system as it accounts for simultaneous Lorentzian and Doppler broadening effects that occur primarily due to the Heisenberg uncertainty principle, particle collisions and their velocity distribution.The imaginary part of the complex error function L (x, y) = Im w (x, y) is also a useful function that can be applied, for example, to describe accurately spectral behavior in refractive index of various dispersive materials (Balazs & Tobias, 1969;Chan, 1986).
The modern plasma theory involving small-amplitude waves propagating through Maxwellian media utilizes the plasma dispersion function that is directly proportional to complex error function as given by (Fried & Conte, 1961) An algorithm for numerical solving of Equation ( 4) can be used as a versatile tool in many practical applications as the complex error function has a direct relation to some functions that are widely applied in Applied Mathematics and Physics.This includes the error function (see the Equations (5a) and (5b) above), the normal distribution function that is commonly utilized in theory of probability and statistics (Weisstein, 2003): the Fresnel integral (Abramowitz & Stegun, 1972;McKenna, 1984): and the Dawson's integral (Abramowitz & Stegun, 1972;Reichel, 1968;McKenna, 1984): As we can see, the error function erf (z) (see Equations ( 5a) and (5b) above), the plasma dispersion function Z (z), the normal distribution function Φ (z), the Fresnel integral F (z) and the Dawson's integral daw (z) are just reformulations of the complex error function w (z).
In this work we obtain a rational approximation of the Voigt/complex error function by Fourier expansion of the exponential function exp −(t − 2σ) 2 and present a master-slave algorithm for its rapid and highly accurate computation.The algorithm, based on the main and supplementary equations, runs the computation in a master slave mode.Such approach enables us to sustain the high-accuracy even when the parameter y tends to zero practically without increase of the computational time.

Rational Approximation
Applying the Fourier transform, the real and imaginary parts of the complex error function w (z) can be expressed as (Armstrong & Nicholls, 1972;Pagnini & Mainardi, 2010;Srivastava & Chen, 1992) and respectively.Combining the real (6) and imaginary (7) parts together results in Although exponential function e −t 2 /4 is non-periodic, we can use, nevertheless, the Fourier expansion as a useful approximation methodology.In particular, it has been shown previously that the Fourier expansion series of the exponential function (Abrarov et al., 2010a(Abrarov et al., , 2010b): where the Fourier expansion coefficients are and τ m ≥ 12 is the margin value (Abrarov, Quine, & Jagpal, 2010a, 2010b;Geetha, 2010), can be effectively applied to approximate the Voigt/complex error function.The relation ( 9) is limited by −τ m ≤ t ≤ τ m since its right side is periodic with period T = 2τ m .
Figure 1  since the function e −t 2 /4 rapidly decreases with increasing t and damps to zero the integrand.Consequently, this leads to (Abrarov & Quine, 2011, 2012) Although this approximation is rapid and highly accurate, it is not based on rational functions due to presence of the exponential function e iτ m z .
Figure 1.The Fourier expansion series approximation (9) for the exponential function exp −t 2 /4 .The red arrow shows the integration domain τ m = 12.The blue shadowed line is the original exponential function exp −t 2 /4 For sufficiently large values of the parameter y, the multiplication of exponential function e −yt to the both sides of Fourier expansion series (9) eliminates the restriction −τ m ≤ t ≤ τ m since the exponential function e −yt rapidly tends to zero with increasing t: Figure 2 illustrates the graphs of e −t 2 /4 e −yt function approximation computed according to approximation (11).As the right side of exponential function approximation ( 9) is periodic with period 2τ m = 24, we observe the peaks at t = {24, 48, 72, 96, 120 . ..}.These peaks vanish as the parameter y increases.This behavior can be observed by comparison of results corresponding to y = 0.025 (green curve) and y = 0.05 (red curve).When the parameter y is large enough, say around or greater than 1, the exponential multiplier e −yt effectively suppresses to zero the entire expression on the right side of Equation ( 11) as t increases.As a result, all peaks practically disappear and only the initial part of the integrand contributes for integration.We can apply this observation to obtain the complex error function based on rational function approximation.
From approximation (9) it follows that Geometrically the peak of the exponential function e −(t−2σ) 2 /4 is shifted towards right with respect to origin by value equal to 2σ.We, therefore, will refer to the value σ as the shift constant.Now if the shift constant σ is large enough, say around or greater than 1, the approximation (13) remains valid even for the worst case scenario when y → 0. This enables us to retain the upper limit of the integral equal to infinity.Thus, substituting approximation (13) into Equation ( 12) leads to the complex error function approximation where are the coefficients that are independent of x and y parameters.By default, the values in approximation ( 14) can be taken as N = 23, τ m = 12 and σ = 2.In contrast to Equation (10), the complex error function approximation ( 14) is based on superposition of rational functions.Therefore it is a rational (non-polynomial) approximation and its application can be efficient for rapid computation.

Implementation
The properties of the rational approximation ( 14) are quite similar to those of the Chiarella and Reichel approximation given by (see Appendix A) where h is a small parameter and is the Heaviside step function that can be defined in terms of Dirac delta function δ (s).Particularly, at N = 23, τ m = 12, σ = 2 and h ≡ π/τ m the approximations ( 10) and ( 15) result in highly accurate values for the typical range of the input parameters 0 < x < 40, 000 and 10 −4 < y < 10 2 applied in quantitative spectroscopy (Wells, 1999;Quine & Drummond, 2002).However, the Chiarella and Reichel approximation contains exponential functions e −z 2 and e −2πiz/h that requires longer time for computation for input array z = {z 1 , z 2 , z 3 . ..}, which size in quantitative radiative transfer applications may contain hundreds of millions or sometimes more than a billion elements.Therefore, the rational approximation ( 14) can be advantages for more rapid computation.In particular, the rational approximation ( 14) is faster in computation than approximations ( 10) and ( 15) by about 20% and 50%, respectively.The enhanced performance is achieved because in contrast to Equations ( 10) and ( 15) the approximation ( 14) is a rational function that contains no trigonometric or exponential functions dependent on input parameters x or y.
Approximations ( 14) and ( 15) have about same coverage over the complex plane.Specifically, the computational tests reveal that both these approximations provide high-accuracy in the first quadrant of the complex plane except a very narrow band area x > 4 ∩ y < 10 −5 along y axis.It should be noted, however, that the deterioration in accuracy in computation of the Voigt/complex error function w (x, y) at y → 0 is a common problem (see for example recently published Amamou et al. (2013).Although computation with multi-precision can resolve this problem, the rapid programming languages such as C/C++, FORTRAN and Matlab need multi-precision libraries of function files that cause installation and validation complexities and also make computation more time consuming.As a result, the multi-precision computation is quite inconvenient in implementation.
We propose instead the following supplementary approximation (see derivation in Appendix B) This approximation has several advantages.First, it needs no partition of the complex plane in a program code in order to separate the narrow band area x > 4 ∩ y < 10 −5 , where approximation ( 14) has insufficient accuracy.Secondly, it is unnecessary to use any other approximation for this narrow band area as the approximation ( 14) alone performs de facto the main computation (the approximation (16) itself does not determine the values K (x, y min ) and L (x, y min ) -its role is supplementary).Choosing small y min , say y min = 10 −5 , we can find K (x, y min ) and L (x, y min ) through main (or master) approximation ( 14) and then once these two values are found they are just substituted into supplementary (or slave) approximation ( 16).As the approximation ( 16) performs few operations of multiplication, summation and the exponentiation e −x 2 that is nearly instant compared to calculation time for determination of w (x, y min ) = K (x, y min ) + iL (x, y min ), it returns a rapid result practically without deceleration.
Cases where the input parameter y is less than 10 −6 are rarely required in practice.For example, the input parameter range in radiative transfer applications is 10 −4 < y < 10 2 (Wells, 1999;Quine & Drummond, 2002) and at N = 23, τ m = 12 and σ = 2 the accuracy of the approximation (14) in this domain can be better than 10 −14 .Therefore, the application of approximation ( 14) alone is sufficient almost in all practical tasks.
It should be noted that in program coding it is convenient to rewrite the approximation (14) as where is defined in order to deal with simplified function ψ (z).Consequently, the approximation procedure in the masterslave approach can be reproduced as According to this expression, if {x, y} x > 4 ∩ y < 10 −5 the main body of the program passes the complex argument x + i (y + σ) to function file associated to ψ (z).Otherwise, if x > 4 ∩ y < 10 −5 the main body of the program passes the complex argument x + i (y min + σ) to function file associated to the slave equation.This programming procedure simplifies the script and makes the computation flow optimized.

Error Analysis
In order to perform error analysis, it is convenient to define the relative errors for the real and imaginary parts of the complex error function as Re w re f. (z) and Im w re f. (z) , respectively, where w re f. (z) is the reference.To generate reference values we can use, for example, the Matlab code in recently published Algorithm 916 (Zaghloul & Ali, 2011).Alternatively, these values can be generated by latest versions of Mathematica computational software that supports error function of complex argument.Tables 1a and 1b provide some results computed by main approximation ( 14) for the real and imaginary parts, respectively.The fourth and fifth columns show the references and relative errors.As we can see from these tables, the numbers generated by approximation ( 14) can match up to the last significant digits with highly accurate reference values.
It is interesting to note that the Weideman's approximation (see Equation (38-I) and corresponding Matlab code in Weideman, 1994): where L = 2 −1/4 N 1/2 and are the expansion coefficients that can be determined elegantly by FFT method, is the only known so far rational approximation of the complex error function that can provide accuracy better than 10 −9 .Due to rational function representation, the Weideman's approximation ( 17) is rapid in computation and, consequently, widely applied in practice.However, as it has been shown in a recent publication, in order to sustain accuracy 10 −6 at y ≥ 10 −5 and 0 ≤ x ≤ 15 the integer N determining the number of the summation terms in Weideman's approximation must be increased up to 32 (Abrarov & Quine, 2011).Furthermore, according to Schreier et al. (2014) in the range y < 10 −5 and 4 < x < 15 the Weideman's approximation (17) at N = 32 cannot provide accuracy better than 8 × 10 −5 .Despite the fact that in approximation ( 14) the integer N = 23 is smaller, it provides essentially higher accuracy.Specifically, the accuracy of the approximation (14) in the same region y ≥ 10 −5 and 0 ≤ x ≤ 15 is better than 10 −12 and, in general, it can be better than 10 −14 at y > 10 −5 as we can see in the Tables 1a and 1b.This confirms that the convergence rate in the rational approximation ( 14) is more rapid than that of in the Weideman's approximation (17).
It should also be noted that the Weideman's approximation ( 17) is a rational function, dependent on the power n.As a result of exponentiations in polynomial terms, the Weideman's approximation needs more computational resources and, consequently, it is more sensitive to the size of input array z = {z 1 , z 2 , z 3 . ..}.In particular, the computational tests show that array programing software Matlab requires longer computational time when size of input array z exceeds several million elements (Abrarov & Quine, 2011).Therefore application of the approximation ( 14) may be more efficient for extended size of the input array as it is a rational function, independent of power n.
Tables 2a and 2b show some results computed by supplementary approximation ( 16) for the real and imaginary parts in the narrow area x > 4 ∩ y < 10 −5 .As one can see from the last columns, the accuracy in this narrow area is better than 10 −9 .Although further improvement in accuracy can be obtained by higher orders of Taylor series in derivation of the slave approximation ( 16), this is absolutely unnecessary since the accuracy better than 10 −9 is more than enough for practical applications.The computational accuracy of the master-slave algorithm can be seen in the Figures 3a-4b.Due to logarithmic values of the relative errors, the color bars in these figures indicate the order of the accuracy.Figure 3a illustrates the relative error for the real part of the complex error function in the range 0 < x < 15 and 10 −6 < y < 15. Figure 3b shows the relative error for the real part of the complex error function in the range 0 < x < 15 and 0 < y < 10 −6 .Figure 4a depicts the relative error for the imaginary part of the complex error function in the range 0 < x < 15 and 10 −6 < y < 15. Figure 4b shows the relative error for the imaginary part of the complex error function in the range 0 < x < 15 and 0 < y < 10 −6 .As one can see from these figures, at y > 10 −5 the accuracy of the master-slave algorithm can be better than 10 −14 while at y < 10 −6 the accuracy still remain high and better than 10 −9 .Thus, the computational test reveals that almost in all first quadrant of the complex plane except a very narrow band area x > 4 ∩ y < 10 −5 along y-axis, the approximation ( 14) is rapidly convergent.

Application at Large |z|
When a complex argument z is large enough by absolute value, say |x + iy| > ∼ 15, the computation of the complex error function w (z) is not problematic and many rapid rational approximations are available in literature for algorithmic implementation.For example, if the high accuracy is required, it is reasonable to apply the approximation by truncating a continued fraction equation (Gautschi, 1970;Poppe & Wijers, 1990a, 1990b): Otherwise, if the high accuracy is not a concern, a rational approximation reported by Hui et al. (1978) may be used for rapid and reasonably accurate computation of the complex error function w (z).
As we can see from this implementation, the actual range of the parameter x in slave approximation ( 16) cannot be smaller than 4 and nor greater than 15 at y ≤ 10 −5 .Therefore the second term ε y min K (x, y min ) in approximation ( 16) remains considerably larger compared to the first term 1 − ε y min even for the worst case scenario when ε → y min and x → 15.This prevents the uncertainty in the real part of the complex error function due to computer limitation for the smallest number.
At |x + iy| > ∼ 1000, the high-accuracy can be obtained by using very simple rational approximations (Zaghloul & Ali, 2011): for the real and imaginary parts of the complex error function, respectively.Obviously, these equations can be combined together as

Coverage Extension at y < 0
Approximations (10) and ( 14) of the complex error function are valid only when the parameters y = Im [z] is positive.However, these approximations applicability are not limited only in the first x > 0, y > 0 and second x < 0, y > 0 quadrants, but can also be extended in the third x < 0, y < 0 and fourth x > 0, y < 0 quadrants.Specifically, from the property of the complex error function (McKenna, 1984;Zaghloul & Ali, 2011): As we can see from the last identity, application of the rapid master slave algorithm, based on approximation (14), can be extended for coverage over the entire complex plane even at negative y.

Generalization
For future prospects in applications, the described methodology can be generalized to an integral of kind: If f (t) is an elementary function or a combination of elementary functions, the integrations on the right side of approximation (20) may result in a rational function.
In fact, the complex error function approximation ( 14) is a specific case that can be obtained from the generalized approximation (20).In particular, by substituting f (t) = 2 exp 2 (ixt − yt) / √ π into integral (20), we get the complex error function w (x, y) since change of the variable t → t/2 results in (see Equation ( 8 Thus, at f (t) = 2 exp 2 (ixt − yt) / √ π we can obtain approximation ( 14) directly from the generalized approximation (20).

Conclusion
We derived a rational approximation of the Voigt/complex error function by Fourier expansion of the exponential function e −(t−2σ) 2 and present master-slave algorithm for its rapid and high-accuracy computation.The error estimation shows that at y > 10 −5 , the proposed approximation (14) can generate numbers with accuracy better than 10 −14 .As the main approximation is a rational function, it can be implemented in rapid computation.The common problem that occurs at y → 0 is resolved by main and supplementary approximations running in a master-slave mode.Such approach enables the high-accuracy computation even at y << 1 practically without deceleration in computation.
depicts this periodic function at τ m = 12 (black curve) and the original function e −t 2 /4 (blue curve).As we can see from this figure, the function e −t 2 /4 can be approximated by relation (9) only within the range t ∈ [−τ m , τ m ] near the origin.However, despite such a restriction we can write

Table 1a
. The real part values generated by master-slave algorithm at N = 23, τ m = 12 and σ = 2

Table 2a .
The real part values generated by master-slave algorithm at N = 23, τ m = 12 and σ = 2 and small y

Table 2b .
The imaginary part values generated by master-slave algorithm at N = 23, τ m = 12 and σ = 2 and small y Milone et al. (1988)nsion coefficients and t m is the margin value that according toMilone et al. (1988)can be taken equal to 6.As a conversion from function e −t 2 /4 to function e −t 2 can be obtained by change of the variable t/2 → t, the corresponding margin value t m and shift constant ς for the Fourier expansion series of the exponential function e −t 2 are changed by factor of two, i.e.: t m = τ m /2 and ς = 2σ.Similar to approximation (13) we can write now Lastly, applying the series expansion (19) in integral (18) leads to the generalized approximation