Add Ziggurat as a method for standard_normal

[stress-tess]

Right now our implementation of standard_normal does a box-muller transform, but numpy uses the Ziggurat algorithm. Ziggurat algorithm's wikipedia page says:

Nevertheless, the algorithm is computationally much faster than the two most commonly used methods of generating normally distributed random numbers, the Marsaglia polar method and the Box–Muller transform, which require at least one logarithm and one square root calculation for each pair of generated values. However, since the ziggurat algorithm is more complex to implement it is best used when large quantities of random numbers are required.

considering it's used by numpy and we're likely to need to generate large quantities of numbers, I think it would be worthwhile. It's also used in numpy's implementation of other distributions such as exponential which will be added in #3183

Numpy's c based implementation can be found here

This will require saving off some constants . The nature of chpl means we'll want these constants on each locale or else we'll kill performance by going back and forth to the lookup tables on locale 0. This has caused us issues previously (see #2354)

[stress-tess]

hmm it seems like box-muller's wikipedia page argues for it's use:

The Box–Muller transform was developed as a more computationally efficient alternative to the inverse transform sampling method.The ziggurat algorithm gives a more efficient method for scalar processors (e.g. old CPUs), while the Box–Muller transform is superior for processors with vector units (e.g. GPUs or modern CPUs).

So it's not immediately obvious which would be better for our use case. It is worth noting that box-muller can't produce values more than 9.419 standard deviations from the mean

[stress-tess]

I think we're going to have to add the stuff need for ziggurat for exponential anyway, so we could always add a method argument like in standard_exponential. Then do a performance test to decide which should be the default

we also be implementing the inverse cdf method / inverse transform sampling method for standard_exponential. So it seems like we'll end up using them all eventually