sysrap/tests/S4MTRandGaussQ.h

0001 #pragma once
0002 /**
0003 S4MTRandGaussQ.h
0004 ==================
0005
0006 Experiments using Geant4 extracts from below classes::
0007
0008    g4-cls G4MTRandGaussQ
0009    g4-cls RandGaussQ
0010
0011 In order to work out how to do "G4RandGauss::shoot(0.0,sigma_alpha)"
0012 GPU side.
0013
0014 The transformQuick/transformSmall looks like a mapping
0015 of InverseErf onto (0,1) rather than (-1, 1)
0016
0017 So maybe can use the CUDA func::
0018
0019     u = curand_uniform(&rng) ;   // u uniform random in 0 to 1
0020     u2 = 2.f*u - 1.f ;           // u2 uniform random in -1 to 1
0021     float v = erfinvf(u2)
0022
0023
0024 * https://mathworld.wolfram.com/InverseErf.html
0025 * https://mathworld.wolfram.com/InverseErfc.html
0026 * https://en.wikipedia.org/wiki/Error_function
0027 * https://en.wikipedia.org/wiki/Q-function
0028
0029
0030 __device__ float erfcf ( float  x )
0031     Calculate the complementary error function of the input argument.
0032
0033 __device__ float erfcinvf ( float  x )
0034     Calculate the inverse complementary error function of the input argument.
0035
0036     * https://docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__SINGLE.html#group__CUDA__MATH__SINGLE_1g31faaaeab2a785191c3e0e66e030ceca
0037
0038     erfinvf(+-0) returns 0
0039
0040     erfinvf(+1)  return +inf
0041
0042     erfinvf(-1) returns -inf
0043
0044
0045
0046 __device__ float erfcxf ( float  x )
0047     Calculate the scaled complementary error function of the input argument.
0048
0049 __device__ float erff ( float  x )
0050     Calculate the error function of the input argument.
0051
0052 __device__ float erfinvf ( float  x )
0053     Calculate the inverse error function of the input argument.
0054
0055
0056 **/
0057
0058 struct S4MTRandGaussQ
0059 {
0060     static double shoot( double mean, double stdDev );
0061     static double shoot();
0062     static double transformQuick (double r);
0063     static double transformSmall (double r);
0064 };
0065
0066 inline double S4MTRandGaussQ::shoot(double mean, double stdDev)
0067 {
0068     return shoot()*stdDev + mean;
0069 }
0070
0071