{'fraction_non_alphanumeric': 0.036097510157308055, 'fraction_numerical': 0.016390596242664167, 'mean_word_length': 4.668799212598425, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 3, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Approximating inverse cumulative distribution functions to produce approximate random variables

For random variables produced through the inverse transform method, approximate random variables are introduced, which are produced by approximations to a distribution's inverse cumulative distribution function. These approximations are designed to be computationally inexpensive, and much cheaper than exact library functions, and thus highly suitable for use in Monte Carlo simulations. Two approximations are presented for the Gaussian distribution: a piecewise constant on equally spaced intervals, and a piecewise linear using geometrically decaying intervals. The error of the approximations are bounded and the convergence demonstrated, and the computational savings measured for C and C++ implementations. Implementations tailored for Intel and Arm hardwares are inspected, alongside hardware agnostic implementations built using OpenMP. The savings are incorporated into a nested multilevel Monte Carlo framework with the Euler-Maruyama scheme to exploit the speed ups without losing accuracy, offering speed ups by a factor of 5--7. These ideas are empirically extended to the Milstein scheme, and the Cox-Ingersoll-Ross process' non central chi-squared distribution, which offer speed ups by a factor of 250 or more.


Introduction
Random number generation is a fundamental building block of a wide range of computational tasks, in particular financial simulations [7,32,43,72]. A frequent performance bottleneck is the generation of random numbers from a specific statistical distribution such as in Monte Carlo simulations. While computers have many excellent and fast implementations of random number generators for the uniform distribution on the interval (0, 1), sampling random variables from a generic distribution is often computationally much more expensive.
We consider random variables produced using the inverse transform method [32], which enables sampling from any uni-variate distribution, and thus is very widely applicable. Additionally, the inverse transform method is crucial for quasi-Monte Carlo simulations [29,47] as no samples are rejected and the low-discrepancy property is preserved [66]; such applications are common in financial simulations [43,72]. Furthermore, we demonstrate how the inverse transform method is particularly well suited to analysis, and that it naturally provides a coupling mechanism for multilevel Monte Carlo applications for a range of stochastic processes, statistical distributions, and numerical schemes.
Our analysis focuses on the Gaussian distribution (a.k.a. the Normal distribution), and the motivations are threefold. Firstly, the distribution is representative of several continuous distributions, and due to the central limit theorem is often the limiting case. Secondly, it is analytically tractable and often admits exact results or is amenable to approximation. Lastly, it is ubiquitous in both academic analysis and scientific computation, with its role cemented within Itô calculus and financial simulations.
To produce Gaussian random variables will require the Gaussian distribution's inverse cumulative distribution function. Constructing approximations accurate to machine precision has long been investigated by the scientific community [8,27,35,48,54,69], where the de facto routine implemented in most libraries is by Wichura [69]. While some applications may require such accurate approximations, for many applications such accuracy is excessive and unnecessarily costly, as is typically the case in Monte Carlo simulations.
To alleviate the cost of exactly sampling from the Gaussian distribution, a popular circumvention is to substitute these samples with random variables with similar statistics. The bulk of such schemes follow a moment matching procedure, where the most well known is to use Rademacher random variables (which take the values ±1 with equal probability [44, page XXXII], giving rise to the weak Euler-Maruyama scheme), matching the mean and variance. Another is to sum twelve uniform random variables and subtract the mean [52, page 500], which also matches the mean and variance, and is still computationally cheap. The most recent work in this direction is by Müller et al. [51], who produce either a three or four point distribution, where the probability mass is positioned so the resulting distribution's moments match the

Approximate Gaussian random variables
The inverse transform method [32, 2.2.1] produce univariate random variable samples from a desired distribution by first generating a uniform random variable U on the unit interval (0, 1), and then computing C −1 (U ), where C −1 is the inverse of the distribution's cumulative distribution function (sometimes called the percentile or quantile function). We will focus on sampling from the Gaussian distribution, whose inverse cumulative distribution function we denote by Φ −1 : (0, 1) → R (some authors use N −1 ), and similarly whose cumulative distribution function and probability density function we denoted by Φ and φ respectively.
Our key proposal is to use the inverse transform method with an approximation to the inverse cumulative distribution function. As the resulting distribution will not exactly match the desired distribution, we call random variables produced in this way approximate random variables, and those without the approximation as exact random variables for added clarity. In general, we will denote exact Gaussian random variables by Z and approximate Gaussian random variables by Z. The key motivation for introducing approximate random variables is that they are computationally cheaper to generate than exact random variables. Consequently, our key motivating criterion in forming approximations will be a simple mathematical construction, hoping this fosters fast implementations. As such there is a trade-off between simplicity and fidelity, where we will primarily be targetting simplicity.
In this section, we present two approximation types: a piecewise constant, and a piecewise linear. For both we will bound the L p error, focusing on their mathematical constructions and analyses. Their implementations and utilisation will be detailed in sections 3 and 4. Both analyses will share and frequently use approximations for moments and tail values of the Gaussian distribution, which we gather together in section 2.1. Thereafter, the piecewise constant and linear approximations are analysed in sections 2.2 and 2.3.

Approximating tail values and high order moments
We require bounds on the approximations of tail values and high order moments. We use the notation of Giles et al. [31] that f (z) ≈ g(z) denotes limz→∞ f (z) g(z) = 1. Our key results are lemmas 2.1 and 2.2 which bound the tail values and high order moments. Lemma 2.1 is an extension of a similar result by Giles et al. [31, lemma 7], extended to give enclosing bounds. Similarly, lemma 2.2 is partly an extension of a result by Giles et al. [31, lemma 9], but extended to arbitrarily high moments rather than just the second. As such, neither of these lemmas are particularly noteworthy in themselves and their proofs resemble work by Giles et al. [31, appendix A]. Similarly, our resulting error bounds for the piecewise constant approximation in section 2.2 will closely resemble a related result by Giles et al. [31, theorem 1]. However, our main result for the piecewise linear approximation in section 2.3 is novel and will require these results, and thus we include them here primarily for completeness.
) and hence we obtain the first inequality.

Piecewise constant approximations on equal intervals
Mathematically it is straightforward to motivate a piecewise constant approximation as the simplest possible approximation to use, especially using equally spaced intervals. As the range of values will go from a continuous to a discrete set, we say the distribution has become quantised and denote our approximation as Q ≈ Φ −1 , where for a uniform random variable U ∼ U (0, 1) we have Z := Φ −1 (U ) and Z := Q(U ). A preview of such an approximation is shown in figure 1(a), where the error will be measured using the L p norm f p := ( 1 0 |f (u)| p du) 1/p . Before presenting theorem 2.3, we can briefly comment on the error seen in figure 1(b). Specifically looking at the root mean squared error (RMSE), corresponding to the L 2 norm, we can see that increasing the number of intervals from 2 to 10 3 gives a drop of 10 2 in the RMSE. For our multilevel Monte Carlo applications in section 4, having approximately 1000 intervals gives a very good fidelity, whereas Rademacher random variables have a very low fidelity. Being able to achieve a reasonable fidelity from our approximation ensures that we achieve the largest portion of the possible computational savings offered by our approximations.
As we have already mentioned, the piecewise constant approximation is closely related to the resulting approximation produced by Giles et al. [31], whose approximation arises from considering uniform random variables truncated to a finite number of bits of precision. Thus our main result from this section, theorem 2.3, closely resembles a related result by Giles et al. [31, theorem 1]. To put our extension into context, we paraphrase the similar result from Giles et al. [31], which is that for a piecewise constant approximation using 2 q intervals, for some integer q 1, then for constant values equal to each interval's midpoint value they have Z − Z 2 2 = O(2 −q q −1 ). Our result from theorem 2.3 extends this to Z − Z p p = O(2 −q q −p/2 ) for p ≥ 2, and numerous other possible constant values other than the midpoint's. Our (a) A piecewise constant approximation using 8 intervals.
10 0 10 1 10 2 10 3 Intervals The L p error and the bound from theorem 2.3. result enables us to increase the order of the error to arbitrarily high norms and is interesting in its own right. It shows that as the intervals become increasing small (corresponding to q → ∞), the dominant term effecting the L p error is the geometric decay O(2 −q/p ) due to the error in the two end intervals, and thus the convergence exists but is slower in higher norms, with the polynomial O(q −1/2 ) term being comparatively negligible (as we can see in figure 1(b)). Additionally, in the related analysis incorporating approximate random variables into a nested multilevel Monte Carlo framework by Giles and Sheridan-Methven [28,63], their bounds on the variance of the multilevel Monte Carlo correction term (discussed more in section 4) rely on the existence of the L p error for p > 2. Hence, while this strengthening of the result may appear only slight, it is crucial for nested multilevel Monte Carlo.
We can now present our main result concerning piecewise constant approximations, namely theorem 2.3. In this we will leave the interval values largely unspecified, and later demonstrate in corollary 2.4 several choices fit within the scope of theorem 2.3. Theorem 2.3. Let a piecewise constant approximation Q ≈ Φ −1 use 2 q equally spaced intervals for some integer q > 1. Denote the intervals On each interval the approximation constant is Q k := Q(u) for u ∈ I k . We assume there exists a constant C independent of q such that:

Then for any even integer
We use condition (1) to reduce our considerations to the domain ( 1 2 , 1), where Φ −1 is strictly positive and convex, and thus obtain Q − Φ −1 p p = 2 K−1 k=2 q−1 A k + 2AK . As Φ −1 is convex, then from the intermediate value theorem there exists a ξ k ∈ I k such that , then from the mean value theorem for any ) . An identical argument follows for any u ∈ [ξ k , u k+1 ) giving the same bound. Using this to bound K−1 k=2 q−1 A k in our expression for Q − Φ −1 p p gives where the last bound comes from considering a translated integral. Changing integration variables the integral becomes z K 0 φ(z) 1−p dz, from which we can use lemma 2.2 to give where the last bound follows from lemma 2.1. Turning our attention to the final interval's contribution AK , then using Jensen's inequality, condition (3), and lemmas 2.1 and 2.2 we obtain Combining our two bounds for K−1 k=2 q−1 A k and AK into our expression for Q − Φ −1 p p we obtain where the coefficients inside the O-notation are only a function of p and not of q.
Corollary 2.4. For approximate Gaussian random variables Z using constants Z k on intervals I k constructed as either then all finite moments of | Z| are uniformly bounded, and conditions (1-3) of theorem 2.3 are satisfied.
Proof. As all three approximations are anti-symmetric about u = 1/2, and therefore all satisfy conditions (1, we consider the domain ( 1 2 , 1) where Z and Z are both positive. Letting the interval I k lie in this domain, then from Jensen's inequality and the law of iterated expectations we obtain for any 1 ≤ n < ∞, where the last inequality is a standard result. Mirroring this result to the full domain (0, 1), we can directly see that the moments of | Z E | are uniformly bounded. Furthermore, as Φ −1 is increasing and convex in [ 1 2 , 1), then the value Z E k is an upper bound on Z C k and Z I k , and so these too have uniformly bounded moments. It is immediately clear that all three choices satisfy condition 2) and the lower bound in 3). As Z E k is an upper bound for the other two choices, it will suffice to show that this satisfies the upper bound in 3). Inspecting the difference between the smallest value in the final IK interval, namely zK , and the average value QK = E(Z | Φ(z) ∈ IK ), we obtain where the first inequality follows from lemma 2.2, the second from lemma 2.1, and the last uses (1−z −2 K ) −1 ≤ 9 8 for q ≥ 10. We bound the final 1 z K term using lemma 2.1 for q ≥ 10 to give 1 z K < 1.1 q −1/2 and therefore QK ≤ zK + 1.25 q −1/2 .
We remark that the Z C k construction from corollary 2.4 is the same as the value used in the analysis of Giles et al. [31, (4)]. In any piecewise constant approximations we use, such as that in figure 1(a), we will use the constants defined by the Z E k construction from corollary 2.4. Furthermore, we can see that our bound from theorem 2.3 appears tight in figure 1(b).

Piecewise linear approximations on geometric intervals
Looking at the piecewise constant approximation in figure 1(a), it is clear there are two immediate improvements that can be made. The first is to use a piecewise linear approximation, which is considerably more appropriate for the central region. Secondly, the intervals should not be of equal sizes, but denser near the singularities at either end. We will make both these modifications in a single step, where we will construct a piecewise linear approximation with geometrically decreasing intervals which are dense near the singularities. For brevity we will denote this just as the piecewise linear approximation. An example piecewise linear approximation using 8 intervals is shown in figure 2(a). The precise nature of the interval widths, and how the linear functions are fitted will be detailed shortly, but by direct comparison against figure 1(a) it is clear that the fidelity of a piecewise linear approximation is much better than the piecewise constant.
The main result from this section will be theorem 2.5, which will bound the L p error of our piecewise linear approximation. The proof will proceed in a similar fashion the proof of theorem 2.3, where we will bound the sum of the central intervals and the end intervals separately. For the central intervals we will use the Peano kernel theorem to bound the point wise error, and in the end intervals several results will be a mixture of exact results and bounds from lemmas 2.1 and 2.2.
Theorem 2.5. For an approximation D ≈ Φ −1 which is anti-symmetric about u = 1/2, with K intervals in (0, 1 2 ), we define the k-th interval I k := [ r k 2 , r k−1 2 ) for k ∈ {1, 2, . . . , K − 1} and IK := (0, r K−1 2 ) for some decay rate r ∈ (0, 1). Each interval uses a piecewise linear approximation D k (u) ≡ D(u) for any u ∈ I k . The gradient and intercept in each interval is set by the L 2 minimisation D k := argmin D ∈P 1 I k |Φ −1 (u) − D (u)| 2 du where P1 is the set of all 1-st order polynomials. Then we have for any 2 ≤ p < ∞ Proof. Considering the domain (0, 1 2 ), we split the contribution into those from the intervals without the singularity, and that from the final interval with the singularity, where where the factors of 2 correct for us only considering the lower half of the domain (0, 1).
(a) A piecewise linear approximation using 8 intervals.  Beginning with the non-singular intervals, we express the pointwise error using the Peano kernel theorem [41,57], which we will later bound. For notational simplicity, we denote the approximated function as f , where f ≡ Φ −1 , and a given interval as To determine the intercept and gradient, we use that they are L 2 optimal, and so the functional derivatives of b a |ε(u)| 2 du with respect to α and β are zero, giving the simultaneous equations It is important to notice that because we chose the L 2 norm, these are a set of linear simultaneous equations, and thus α and β are linear functionals, thus showing that L is linear, (a requirement of the Peano kernel theorem). Evaluating these for the kernel function (f → (· − ξ) + ) gives Thus, the pointwise error is where to achieve the last equality we rescaled our interval [a, b] → [0, 1] and variables η →η whereη := η−a b−a . Taking the absolute value and applying Jensen's inequality immediately gives where for the first inequality we used that d 2 du 2 Φ −1 is maximal at the lower boundary, and for the second inequality we bound this by the maximum with respect to a (which is at a ≈ 0.177). Using this expression for the pointwise error in our summation of the non-singular intervals gives (as r → 1) We now consider the interval [0, b] ≡ IK containing the singularity at 0, which has intercept and gradient α = These integrals can be calculated exactly (with a change of variables), where denoting For the gradient, as the interval becomes ever smaller and b → 0, we can use lemma 2.1 to give where in the first approximation we used φ( √ 2z) ≡ √ 2πφ 2 (z). Interestingly, this means the range of values βb ≈ − 3 z b → 0, and our approximation "flattens" relative to the interval [0, b] as b → 0.
With the intercept and gradient in the singular interval IK known exactly, we define the two points u− and u+ where the error is zero, where 0 < u− < u+ < b, and there are two as Φ −1 is concave in (0, 1 2 ). Corresponding to these we define z− := Φ −1 (u−) and z+ := Φ −1 (u+), where −∞ < z− < z+ < z b < 0. Thus in the singular interval we have Using lemmas 2.1 and 2.2, then for the first of these integrals we obtain and for the second integral we similarly obtain Combining the results for the central intervals and the singular interval we obtain where for the second equality we used lemma 2.1 in the limit r K−1 → 0.
We can see from theorem 2.5 that the O((1 − r) 2p ) term comes from the central regions, and is reduced by taking r → 1, and the o(r K−1 ) term is from the singular interval, and is reduced by taking r K−1 → 0. The key point of interest with this result is that the error from the central regions and the singular region is decoupled. In order to decrease the overall error, it is not sufficient to only increase the number of intervals (K → ∞), which would only improve the error from the singular interval, but the decay rate must also be decreased (r → 1). The independence and interplay of these two errors is important for balancing the fidelity between the central and edge regions.
It is possible to generalise this construction and analysis to piecewise polynomial approximations, constructing the best L 2 approximation on each interval. We could also require the approximations to be continuous over the entire interval (0, 1) by turning our previous minimisation into a constrained minimisation, with a coupled set of linear constraints. While such a continuous approximation may be more aesthetically pleasing, it is of no practical consequence for the inverse transform method, and by definition will have a worse overall L 2 error than the discontinuous approximation from theorem 2.5, and thus we choose not to do this.
The errors from piecewise polynomial approximations for various polynomial orders and interval numbers are shown in figure 2(b), where we have set the decay rate r = 1 2 . As we increase the number of intervals used, then for the piecewise linear function, the L 2 error plateaus to approximately 10 −2 for 16 intervals, which is approximately equal to the error for the piecewise constant approximation using 10 3 intervals. Thus we can appreciate the vast increase in fidelity that the piecewise linear approximation naturally offers over the piecewise constant approximation. Furthermore, this plateau demonstrates that the central regions are limiting the accuracy of the approximation. Inspecting the approximation using 16 intervals in (0, 1 2 ), we see that as we increase the polynomial order from linear to cubic there is a considerable drop in the error approximately equal to a factor of 10 2 , where the piecewise cubic approximation is achieving L 2 errors of around 10 −4 . However, we see that increasing the polynomial order any further does not lead to any significant reduction in the L 2 error, indicating that in this regime it is the error in the singular interval which is dominating the overall error.

High performance implementations
Our motivation for presenting the piecewise constant and linear approximations was to speed up the inverse transform method, where the approximations were inherently simple in their constructions, with the suggestion that any implementations would likely be simple and performant. However, it is not obvious that the implementations from various software libraries, especially heavily optimised commercial libraries such as Intel's Maths Kernel Library (MKL), should be slow. In this section, we will mention how most libraries implement these routines, and why many are inherently ill posed for modern vector hardware. Our two approximations capitalise on the features where most libraries stumble, avoiding division and conditional branching, which will be the key to their success. We give a brief overview of their implementations in C and showcase their superior speed across Intel and Arm hardwares. While we will briefly outline the implementations here in moderate detail, a more comprehensive suite of implementations and experiments, including Intel AVX-512 and Arm SVE specialisations, are hosted in a centralised repository by Sheridan-Methven [62].

The shortcomings in standard library implementations
Several libraries providers, both open source and commercial, offer implementations of the inverse Gaussian cumulative distribution function, (Intel, NAG, Nvidia, MathWorks, Cephes, GNU, Boost, etc.). The focus of these implementations is primarily ensuring near machine precision is achieved for all possible valid input values, and handling errors and edge cases appropriately. While some offer single precision implementations, most offer (or default to) double precision, and to achieve these precisions algorithms have improved over the years [8,27,35,48,54,69]. The most widely used is by Wichura [69], which is used in GNU's scientific library (GSL) [25], the Cephes library [50] (used by Python's SciPy package [68]), and the NAG library [53] (NAG uses a slight modification).
These implementations split the function into two regions: an easy central region, and the difficult tail regions. Inside the central region a rational Padé approximation is used with seventh order polynomials [69]. In the tail regions, the square root of the logarithm is computed, and then a similar rational approximation performed on this. While these are very precise routines, they have several shortcomings with respect to performance.
The first shortcoming is the division operation involved in computing the rational approximation. Division is considerably more expensive than addition or multiplication, and with a much higher latency [24,71]. Thus, for the very highest performance, avoiding division is preferable. (The implementation by Giles [27] already capitalises on this). Similarly, taking logarithms and square roots is expensive, making the tail regions very expensive, and so are best avoided too.
The second short coming is that the routine branches: requesting very expensive calculations for infrequent tail values, and less expensive calculations for more frequent central values. Branching is caused by using "if-else" conditions, but unfortunately this often inhibits vectorisation, resulting in compilers being unable to issue single instruction multiple data (SIMD) operations [67]. Even if a vectorised SIMD implementation is produced, typically the results from both branches are computed, and the correct result selected a posteriori by predication/masking. Thus, irrespective of the input, both the very expensive and less expensive calculations may be performed, resulting in the relatively infrequent tail values producing a disproportionately expensive overall calculation. This effect is commonly known as warp divergence in GPU programs, and in our setting we will call this vector divergence. The wider the vector and the smaller the data type, (and thus the greater the amount of vector parallelisation), the worse the impact of vector divergence. Noting that half precision uses only 16 bits, both Intel's AVX-512 and Fujitsu's A64FX are 512 bits wide, and Arm's SVE vectors can be up to 2048 bits [56,65], the degree of vector parallelisation can be substantial, and the impact of conditional branching becomes evermore crippling as vector parallelisation increases.
In light of these two shortcomings, we will see our approximations can be implemented in a vectorisation friendly manner which avoids conditional branching and are homogenous in their calculations. Furthermore, their piecewise constant and polynomial designs avoid division operations, and only require addition, multiplication, and integer bit manipulations, and thus result in extremely fast executables. We will see that our piecewise constant approximation will rely on the high speed of querying the cache. Furthermore, for the piecewise linear approximation using a decay rate r = 1 2 and 16 intervals, then in single precision all the coefficients can be held in 512 bit wide vector registers, bypassing the need to even query the cache, and thus are extremely fast provided the vector widths are sufficiently large.

Implementing the piecewise linear approximation
The piecewise linear approximation is a bit more involved than the piecewise constant approximation. Not only do we have a polynomial to evaluate (albeit only linear), but the varying widths of the intervals means identifying which interval a given value corresponds to is more involved. Once the appropriate interval's polynomial's coefficients are found, evaluating the polynomial is trivial, where the linear polynomial can be evaluated using a fused multiply and add (FMA) instruction. Higher order polynomials can be similarly computed (e.g. with Horner's rule).
The primary challenge is determining which interval a given value corresponds to based on the construction from theorem 2.5 for a given decay rate r. As floating point numbers are stored in binary using their sign, exponent, and mantissa, the natural choice most amenable for computation is r = 1 2 , which we call the dyadic rate, producing the dyadic intervals shown in table 1. From table 1, we notice the intervals are only dense near the singularity at 0, but not at 1. This is not problematic, as we said in theorem 2.5 that the approximation is anti-symmetric about u = 1/2, and thus we can use the N intervals in (0, 1 2 ), and if our input value is within ( 1 2 , 1), it is straight forward to simply negate the value computed for the input reflected about 1 2 (equivalent to using the inverse complementary cumulative distribution function).
Such an implementation can handle any value in (0, 1 2 )∪( 1 2 , 1), but unfortunately 1 2 is not included within this. While a mathematician may say the individual value 1 2 has zero measure, from a computational perspective such a value is perfectly feasible input, and quite likely to appear in any software tests. Thus our implementation will correctly handle this value. The reader will notice that when we reflect an input value x > 1 2 about 1 2 by using x → 1 − x, then the value 1 2 will not be in any of the intervals indexed between 1 and N in table 1, but remains in the interval indexed by 0. Thus in our later implementations, the arrays of coefficients will always hold as their first entry (index 0) zero values (Φ −1 ( 1 2 ) = 0) so 1 2 is correctly handled. We found during the development of the implementation, being able to correctly handle this value is of considerable practical importance, so we encourage practitioners to also correctly handle this value.
If we are using the dyadic decay rate, then the interval an input x belongs to is −log 2 (x) − 1. However, to compute this we need not take any logarithms, which are expensive to compute. We can obtain the same result by reading off the exponent bits in the floating point representation, treating these as an integer, and then correcting for the exponent bias. This only involves simple bit manipulations, and interpreting the bits in a floating point representation as an integer. In C this is a technique called type punning, and can be achieved either by pointer aliasing, using a union, or using specialised intrinsic functions. Of these, pointer aliasing technically breaks the strict aliasing rule in C (and C++) [40, 6.5.2.3] [64, pages 163-164]. Similarly, type punning with a union in C89 is implementation defined, whereas in C11 the bits are re-interpreted as desired. In the presented implementation we will leave this detail undefined and just use a macro to indicate the type punning. (The implementation approximating the Gaussian distribution uses pointer aliasing, whereas the later approximation to the non central χ 2 distribution in section 5 uses a union, to demonstrate both possibilities [62]).
Input: Floating point uniform random variable U ∈ (0, 1). Output: Floating point approximate Gaussian random variable Z. 1 Form predicate using U > 1 2 . 2 Reflect about 1 2 to obtain U ∈ (0, 1 2 ]. 3 Interpret U as an unsigned integer. 4 Read the exponent bits using bit wise AND. 5 Right shift away the mantissa's bits. 6 Obtain an array index by correcting for exponent bias. 7 Cap the array index to avoid overflow. 8 Read the polynomial's coefficients. 9 Re-interpret U as a float. 10 Form the polynomial approximation Z. 11 Correct sign of approximation based on the predicate. Algorithm 1: Piecewise polynomial approximation using dyadic intervals.
Overall then, the general construction from theorem 2.5 is given in algorithm 1. Furthermore, a single precision C implementation is shown in code 2, where we assume the single precision floats are stored in IEEE 32 bit floating point format [39].
The reason why we decide to implement the approximation in single precision using coefficient arrays of 16 entries is because each requires only 512 bits (16 × 32 bits) to store all the possible values for a given monomial's coefficient. The significance of 512 bits cannot be overstated, as it is the width of an AVX-512 and A64FX vector register. Thus, instead of querying the cache to retrieve the coefficients, they can be held in vector registers, bypassing the cache entirely, and achieving extremely fast speeds. Recognising the coefficients can be stored in a single coalesced vector register is currently largely beyond most compilers' capabilities using OpenMP directives and compiler flags alone. However, in the repository by Sheridan-Methven [62], specialised implementations using Intel vector intrinsics and Arm inline assembly code achieve this, obtaining the ultimate in performance.

Performance of the implementations
Both the piecewise constant and linear implementations in codes 1 and 2 are non branching, vector capable, and use only basic arithmetic and bit wise operations, and thus we anticipate their performance should be exceptionally good. Indeed, their performance, along with several other implementations', is shown in table 2, with experiments performed on an Intel Skylake Xeon Gold CPU and an Arm based Cavium ThunderX2 [62]. typedef unsigned int uint32 ; // A s s u m i n g 32 bit floats . typedef float float32 ; // A s s u m i n g 32 bit i n t e g e r s .

# define TABLE_SIZE 16 # define N_MANTISSA_32 23 // For IEEE 754 # define F L O A T 3 2 _ E X P O N E N T _ B I A S 127 // For IEEE 754 # define F L O A T 3 2 _ E X P O N E N T _ B I A S _ T A B L E _ O F F S E T ( F L O A T 3 2 _ E X P O N E N T _ B I A S -1) # define TA BL E_ M AX _I ND E X ( TABLE_SIZE -1) // Zero i n d e x i n g . # define F L O A T 3 2 _ A S _ U I N T 3 2 ( x ) (...)
// Type punning . Looking at the results from table 2, of all of these, the piecewise linear implementation using Intel vector intrinsics achieves the very fastest speeds, closely approaching the maximum speed of just reading and writing. Unsurprisingly, the freely available implementations from Cephes and GSL are not competitive with the commercial offerings from Intel. Nonetheless, even in single precision, our approximations, on Intel hardware, consistently beat the performance achieved from the Intel high accuracy (HA) or low accuracy (LA) offerings. Comparing the high accuracy Intel offering and the piecewise linear implementation, there stands to be a speed up by a factor of seven by switching to our approximation. These results vindicate our efforts, and that our simple approximations offer considerable speed improvements. It is also needless to say, that compared to the freely available open source offerings, the savings become even more significant.

Multilevel Monte Carlo
One of the core use cases for our high speed approximate random variables is in Monte Carlo applications. Frequently, Monte Carlo is used to estimate expectations of the form E(P ) of functionals P which act on solutions X of stochastic differential equations of the form dXt = a(t, Xt) dt + b(t, Xt) dWt for given drift and diffusion processes a and b. The underlying stochastic process Xt is itself usually approximated by some X using a numerical method, such as the Euler-Maruyama or Milstein schemes [7,44,45]. These approximation schemes simulate the stochastic process from time t = 0 to t = T over N time steps of size ∆t = δ = T N , where the incremental update at the n-th iteration requires a Gaussian random variable Zn to simulate the increment to the underlying Wiener process Wt, where ∆Wn = √ δZn. Such types of Monte Carlo simulations are widespread, with the most famous application being to price financial options.
Approximate random variables come into this picture by substituting the exact random variable samples in the numerical scheme with approximate ones. This facilitates running faster simulations, at the detriment of introducing error. However, using the multilevel Monte Carlo method [26], this error can be compensated for with negligible cost. Thus, the speed improvements offered by switching to approximate random variables can be largely recovered, and the original accuracy can be maintained. A detailed inspection of the error introduced from incorporating approximate Gaussian random variables into the Euler-Maruyama scheme and the associated multilevel Monte Carlo analysis is presented by Giles and Sheridan-Methven [28,63]. As such, we will only briefly review the key points of the setup, and focus on detailing the resultant computational savings that can be expected from using approximate random variables.
For the Euler-Maruyama scheme, the unmodified version using exact Gaussian random variables Z produces an approximation X, whereas the modified scheme using approximate random variables Z produces an approximation X, where tn := nδ. The regular multilevel Monte Carlo construction varies the discretisation between two levels, producing fine and coarse simulations. The functional P would act on each path simulation, producing the fine and coarse approximations P f and P c respectively. If the path simulation uses the exact random variables we denote these as P f and P c , and alternatively if it uses approximate random variables as P f and P c . In general there may be multiple tiers of fine and coarse levels, so we use l to index these, where increasing values of l correspond to finer path simulations. Thus for a given l we have P l ≡ P f and P l−1 ≡ P c , and similarly P l ≡ P f and P l−1 ≡ P c . If we have levels l ∈ {0, 1, 2, . . . , L} and use the convention P−1 := P−1 := 0, then Giles and Sheridan-Methven [28] suggest the nested multilevel Monte Carlo where the first approximation is the regular Monte Carlo procedure [32], the first equality is the usual multilevel Monte Carlo decomposition [26], and the final equality is the nested multilevel Monte Carlo framework [28,63]. Giles and Sheridan-Methven [28,63] show that the two way differences in the regular and nested multilevel Monte Carlo settings behave almost identically, and the main result of their analysis is determining the behaviour of the final four way difference's variance [28, lemmas 4.10 and 4.11] [63, corollaries 6.2.6.2 and 6.2.6.3]. They find that for Lipschitz continuous and differentiable functionals P f − P c − P f + P c p ≤ O(δ 1/2 Z − Z p ) for any p such that 2 ≤ p < p < ∞, and that for Lipschitz continuous but non-differentiable functionals that for any > 0. We can see that in all circumstances covered by their analysis that there is a dependence on the approximation error Z − Z p for some L p norm where p > 2.

Expected time savings
The regular multilevel estimatorθ isθ , and the nested multilevel estimatorθ is where m l , m l , and M l are the number of paths generated, each with a computational time cost ofĉ l ,c l , and C l , and variancev l ,ṽ l , and V l respectively, and all terms with the same suffix (i, l) are computed using the same uniform random numbers. Each of these estimators will have an error due to the finite number of paths used and the approximation scheme employed. The total error arising from these two factors is commonly referred to as the variance bias trade-off, where the mean squared error (MSE) of an estimator θ ∈ {θ,θ} is given by MSE(θ) = V(θ) + Bias 2 (θ) [32, page 16]. Setting the desired MSE to ε 2 and choosing the maximum simulation fidelity such that we just satisfy Bias 2 (θ) ≤ ε 2 2 , then we can derive an expression for the total computational time T . Forcing T to be minimal is achieved by performing a constrained minimisation with an objective function F. Considering the estimatorθ, the corresponding objective function is F := L l=0 m lĉl + µ( L l=0v l m l − ε 2 2 ), where µ is a Lagrange multiplier enforcing the constraint V(θ) = ε 2 2 . Treating the number of paths as a continuous variable this is readily minimised to give where the minimal number of paths required are given by and hence an overall saving of For modest fidelity approximations whereṽ l ≈v l , the termṽ lcl v lĉl ≈c l c l measures the potential time savings, and the term (1 + ( V l C l /ṽ lcl ) 1/2 ) 2 assesses the efficiency of realising these savings. A balance needs to be achieved, and the approximations should be sufficiently fast so there is the potential for large savings, but of a sufficient fidelity so the variance of the expensive four way difference is considerably lower than the variance of the cheaper two way difference.
Although we have estimates for the costs of each estimator from table 2, the variance will depend on the stochastic process being simulated, the numerical method being used, and the approximation employed. To make our estimates more concrete and quantify the possible variance reductions, we consider a geometric Brownian motion where a(t, Xt) := µXt and b(t, Xt) := σXt for two strictly positive constants µ and σ, where we take µ = 0.05, σ = 0.2, X0 = 1, and T = 1 [26, 6.1]. The coarse level's Weiner increments are formed by pair wise summing the fine level's (which uses the fine time increment δ f ). Importantly, the uniform random variable samples producing the exact Gaussian random variables will be the same as those used for producing the approximate Gaussian random variables, ensuring a tight coupling. Thus, for some approximation Φ −1 ≈ Φ −1 , such as those from section 2, for a single uniform sample Un ∼ U(0, 1), then Zn := Φ −1 (Un) and Zn := Φ −1 (Un), where the Un is the same for both of these.
The variance for the various multilevel terms for different time increments for the Euler-Maruyama and Milstein schemes are shown in figure 3 for various approximations. We consider the underlying process itself, corresponding to the functional P (X) = X. The piecewise constant approximation uses 1024 intervals, which from figure 1(b) has an RMSE of approximately 10 −2 . The piecewise linear approximation uses 15 dyadic intervals (16 coefficients stored to correctly handle 1 2 ), which from figure 2(b) has a similar RMSE of approximately 10 −2 . The piecewise cubic approximation uses the same number of intervals, which from figure 2(b) has an RMSE of approximately 10 −4 . We also include the approximation using Rademacher random variables for comparison.
Fine time increment δ f  We can see from figure 3 that the two way differences exhibit the usual strong convergence order 1 2 and 1 of the Euler-Maruyama and Milstein schemes, as expected [44]. Furthermore, as the functional is differentiable and Lipschitz continuous, this strong convergence rate is preserved for the four way differences [28,63]. Giles and Sheridan-Methven [28,63] derive this result for the Euler-Maruyama scheme, but the analysis for the Milstein scheme remains an open problem. Aside from this, the key point to note from figure 3 is the substantial drop in the variance between the two way and four way variances using our approximations. This drop in the variance is large for the piecewise constant and linear approximations, and huge for the piecewise cubic approximation. The approximation using Rademacher random variables has only a very small drop in variance.
Estimating the best cost savings seen compared to Intel (HA) from table 2, the variance reductions from figure 3(a), and using the simplifications C l ≈ĉ l +c l andṽ l ≈v l , then the estimated speed ups and their efficiencies are shown in table 3. For the Rademacher random variables we use the optimistic approximations that this achieves the maximum speed set by reading and writing and that it offers a variance reduction of 2 −1 . We can see from this that the piecewise linear approximation would give the largest savings, although the savings from the piecewise constant and cubic approximations are quite similar. Notice that while the piecewise cubic is achieving near perfect efficiency, its cost savings are not substantial enough to beat the marginally less efficient piecewise linear approximation which offers the larger potential savings. For all of our approximations the four way difference simulations are required very infrequently. Lastly, while the Rademacher random variables may offer the best time savings, the absence of any considerable variance reduction caused by the extremely low fidelity approximation results in a very inefficient multilevel scheme, demonstrating the balance required between speed and fidelity. These cost savings are idealised in the respect that we have attributed the cost entirely to the generation of the random numbers. While this is quite a common assumption, the validity of this assumption will diminish the faster the approximations become, as the basic cost of the other arithmetic operations becomes significant. Thus, while a practitioner should have their ambitions set to achieve these savings, they should set more modest expectations.

The non central χ 2 distribution
A second distribution of considerable practical and theoretical interest is the non-central χ 2 distribution, which regularly arises from the Cox-Ingersoll-Ross (CIR) interest rate model [19] (and the Heston model [36]). The distribution is parametrised as χ 2 ν (λ), where ν > 0 denotes the degrees of freedom and λ ≥ 0 the non-centrality parameter, where we denote the inverse cumulative distribution function as C −1 ν (·; λ). Having a parametrised distribution naturally increases the complexity of any implementation, exact or approximate, and makes the distribution considerably more expensive to compute.
To gauge the procedure used, the associated function in Python's SciPy package (ncx2.ppf) calls the Fortran routine CDFCHN from the CDFLIB library by Brown et al. [13] (C implementations available [15]). This computes the value by root finding [17, algorithm R] on the offset cumulative distribution function Cν (·; λ), where Cν (·; λ) is itself computed by a complicated series expansion [2, (26.4.25)] involving the cumulative distribution function for the central χ 2 distribution. Overall, there are many stages involved, and as remarked by Burkardt [15, cdflib.c]: "Very large values of [λ] can consume immense computer resources". The analogous function ncx2inv in MATLAB, from its statistics and machine learning toolbox, appears to follow a similar approach based on its description [49, page 4301]. To indicate the costs, on an Intel core i7-4870HQ CPU the ratios between sampling from the non-central χ 2 distribution and the Gaussian distribution (norm.ppf and norminv in Python and MATLAB respectively) are shown in table 4, from which it is clear that the non-central χ 2 distribution can be vastly more expensive than the Gaussian distribution.

Approximating the non-central χ 2 distribution
There has been considerable research effort into quickly sampling from the non-central χ 2 distribution by various approximating distributions [1,10,38,42,[58][59][60]70], with the most notable being an approximation using Gaussian random variables by Abdel-Aty [1], (and a similar type of approximation by Sankaran [59]). With the exception of Abdel-Aty [1] and Sankaran [59], the remaining methods do not directly approximate the inverse cumulative distribution function, so do not give rise to an obvious coupling mechanism for multilevel Monte Carlo simulations, whereas our framework does. The schemes by Abdel-Aty [1] and Sankaran [59], both of which require evaluating Φ −1 , are appropriate for comparison to our scheme. Hence, our approximation scheme, while modest in its sophistication, appears novel research in this direction.
For approximating the non-central χ 2 distribution, we simplify our considerations by taking ν to be fixed, and thus the distribution then only has the parameter λ varying. While this may appear a gross simplification, in the financial applications ν is a model constant independent of the numerical method's parameters, and thus this simplification is appropriate.
There are only two difficulties with constructing such an approximation. The first is that the distribution is no longer anti-symmetric about u = 1/2, which is easily remedied by constructing two approximations: one for (0, Time increment δ  The performance of C and C++ implementations are shown in table 5 on Intel Skylake hardware compiled with g++. Neither Intel, NAG, nor GSL offer an exact implementation of the inverse cumulative distribution function for the non-central χ 2 distribution. Consequently, for our C implementation we compare against the cdfchn function from CDFLIB [13,15], and for our C++ implementation the quantile function from the Boost library [3] acting on a non central chi squared object (boost/math/distributions). Like the other implementations, Boost computes the inverse by a numerical inversion of the cumulative distribution function. We can see from table 5 that our approximation is orders of magnitude faster than the exact functions in both languages.
The fidelity of our non central χ 2 approximations C −1 ν (·; λ) across the range of parameter values are quantified by the errors shown in table 6. A piecewise linear approximation achieves a consistently good RMSE, as seen in table 6(a). Similarly, the piecewise cubic approximation demonstrates the same behaviour in table 6(b), and has an RMSE which is typically an order of magnitude less than the piecewise linear approximation's.
The approximations from Abdel-Aty [1] and Sankaran [59] are also shown in tables 6(c) and 6(d) respectively. We can see from table 6(c) that the approximation from Abdel-Aty, while improving for increasing ν, performs relatively poorly for large λ. Conversely, the approximation by Sankaran in table 6(d) is notably better than that by Abdel-Aty and often comparable to our piecewise cubic approximation, although demonstrating a much less consistent performance across the range of parameter values explored.
Having compared the error coming from the approximation by Sankaran [59], we can inspect the computational complexity of the approximation. The approximation can be expressed as where c1 := 1 − 2 3 (k + λ)(k + 3λ)(k + 2λ) −2 , c2 := (k + 2λ)(k + λ) −2 , and c3 := (c1 − 1)(1 − 3c1). Just as in our polynomial approximations, it contains an identical mix of basic arithmetic operations, but a Φ −1 evaluation. This can either be calculated exactly, or for improved speed it can potential be further approximated itself, but this then compounds a second level of approximation on top of the original. Furthermore, the expression contains an exponentiation by a factor of 1/c1, which will not in general be an integer. It can be expected that such a non-integer exponentiation operation will be quite expensive. Consequently, we anticipate this would be measurably slower than our polynomial approximations.
Contrasting the approximations between Sankaran [59] and our piecewise cubic, the piecewise cubic demonstrates a more consistent RMSE alongside a simpler computational composition, makes our piecewise polynomial the more attractive of the two.
Rather than approximating the non-central χ 2 distribution with an exact Gaussian distribution when using the Euler-Maruyama scheme, we propose using approximate non-central χ 2 random variables, such as those from a piecewise linear approximation. This has the benefit of offering vast time savings, whilst introducing far less bias than the Euler-Maruyama scheme. The piecewise linear approximation is on average hundreds of times faster than the exact function, as was seen in table 4, giving vast savings, while still achieving a high fidelity, as was seen from table 6. We can generate path simulations of the CIR process using the exact non-central χ 2 distribution, piecewise linear and cubic approximations, and the truncated Euler-Maruyama scheme from Higham et al. [37] (where √ Xn → |Xn|), where we set κ = 1 2 and θ = σ = T = X0 = 1 (just satisfying the Feller condition 1 ). These respectively give rise to the exact underlying process XT and the three approximations X linear T , X cubic T , and X EM T . The variances of the underlying process and the differences between this and the approximations are shown in figure 4(b). Note that the difference terms are evaluated using the same time increments, and there is no mixing of coarse and fine paths in these two-way differences. The results of figure 4(b) do not change if we instead satisfy the Feller condition with strict inequality by setting e.g. κ = 3 2 . From figure 4(b) we can see the exact underlying process' variance does not vary with the discretisation, which is to be expected as the simulation is from the exact distribution, and thus has no bias. We can see that the Euler-Maruyama approximation to the process exhibits a strong convergence order 1 2 , as expected [33,37]. For the piecewise linear approximation, as the fidelity will not vary with the discretisation, then subsequently the variance also does not vary. Furthermore, given the high fidelity seen in table 6(a), the drop in the variance is approximately 2 −15 . Similarly, the even higher fidelity cubic approximation demonstrates the same behaviour in table 6(b), but with a drop in variance of 2 −25 . These findings demonstrate that the most substantial variance reductions come from using our piecewise polynomial approximations to the non central χ 2 distribution's inverse cumulative distribution function, rather than using the Gaussian distribution with the Euler-Maruyama scheme.
We can repeat these same variance comparisons with a different functional of the terminal value which is path dependent. A common example in financial simulations to consider the arithmetic average of the path, (an Asian option in mathematical finance [32]). Doing so produces an exactly identical variance structure as seen in figure 4(b). From this we conclude that our findings seen and explained thus far appear to carry over immediately to simulations requiring an entire path simulation, not only a terminal value, justifying the use of our approximations in time marching schemes.
Quantifying the savings that can be expected, we conservatively approximate the cost savings for both the linear and cubic approximations as a factor of 300 using table 5. Coupling this with the variance reductions seen in figure 4(b), the anticipated cost savings are shown in table 7. We can see from this that both offer impressive time savings, and the piecewise cubic again achieves a near perfect efficiency. However, the piecewise linear approximation's efficiency, while good, is losing an appreciable fraction, making the piecewise cubic the preferred choice. Ultimately, both offer vast time savings by factors of 250 or higher.

Conclusions
The expense of sampling random variables can be significant. In the work presented, we proposed, developed, analysed, and implemented approximate random variables as a means of circumventing this for a variety of modern hardwares. By incorporating these into a nested multilevel Monte Carlo framework we showcased how the full speed improvements offered can be recovered with near perfect efficiency without losing any accuracy. With a detailed treatment, we showed that even for basic simulations of geometric Brownian motions requiring Gaussian random variables, speedups of a factor by Giles and Sheridan-Methven [28]. Lastly, these approximations were able to achieve very high fidelities, with the piecewise linear approximation using a geometric sequence of intervals providing high resolution of singularities. With the approximations detailed from a mathematical perspective, we highlighted two possible implementations in C [62]. The benefit of these approximations is that they are by design ideally suited for modern vector hardware and achieve the highest possible computational speeds. They can be readily vectorised using OpenMP SIMD directives, have no conditional branching, avoid division and expensive function evaluations, and only require simple additions, multiplications, and bit manipulations. The piecewise constant and linear implementations were orders of magnitude faster than most freely available open source libraries, and typically a factor of 5-7 times faster than the proprietary Intel Maths Kernel Library, achieving close to the maximum speed of reading and writing to memory. This speed comes from the simplicity of their operations, and heavy capitalisation on the fast speeds of querying the cache and vector registers.
Incorporating approximate random variables into the nested multilevel Monte Carlo framework by Giles and Sheridan-Methven [28], the low errors and fast speeds of the approximations can be exploited to obtain their full speed benefits without losing accuracy. Inspecting the appropriate multilevel variance reductions, we demonstrated how practitioners can expect to obtain speed improvements of a factor of 5-7 by using approximate Gaussian random variables, where the fastest approximation was the piecewise linear approximation. This appears to be the case when using either the Euler-Maruyama or Milstein schemes.
Considering the Cox-Ingersoll-Ross process [19], this is known to give rise to the non-central χ 2 distribution, which is an example of a parametrised distribution, and is also extremely expensive to sample from, as we demonstrated. We applied our approximation framework to this to produce a parametrised approximation. The error of using our approximate random variables was orders of magnitude lower than approximating paths using the Euler-Maruyama scheme, showing our approximate random variables are considerably more suitable for generating path simulations than the Euler-Maruyama scheme. This circumvents the problem of the Euler-Maruyama scheme having a very large bias for the Cox-Ingersoll-Ross process [12]. Furthermore, our implementation was substantially quicker than those by CDFLIB and Boost, offering a speed improvement of 250 times or higher.
All of the code to produce these figures, tables, and approximations is freely available and hosted in repositories by Sheridan-Methven [61,62], along with an introductory guide aimed at practitioners [62].