File indexing completed on 2024-06-16 04:17:22
0001 /** 0002 * SPDX-FileCopyrightText: 1982-1989 Donald H. House <x@unknown.com> 0003 * SPDX-FileCopyrightText: 1989 Robert Allen <x@unknown.com> 0004 * SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 #include "gauss.h" 0008 #include <math.h> 0009 #include <stdlib.h> 0010 #include <QtGlobal> 0011 0012 #ifdef Q_OS_WIN 0013 // quoting DRAND48(3) man-page: 0014 // These functions are declared obsolete by SVID 3, 0015 // which states that rand(3) should be used instead. 0016 #define srand48 srand 0017 #define drand48() (static_cast<double>(qrand()) / static_cast<double>(RAND_MAX)) 0018 #endif 0019 0020 double gauss::gaussian(double mean, double std, int seed) 0021 { 0022 0023 #define itblmax 20 /* length - 1 of table describing F inverse */ 0024 #define didu 40.0 /* delta table position / delta ind. variable */ 0025 /* itblmax / 0.5 */ 0026 0027 /* interpolation table for F inverse */ 0028 static double tbl[] = {0.00000E+00, 6.27500E-02, 1.25641E-01, 1.89000E-01, 0029 2.53333E-01, 3.18684E-01, 3.85405E-01, 4.53889E-01, 0030 5.24412E-01, 5.97647E-01, 6.74375E-01, 7.55333E-01, 0031 8.41482E-01, 9.34615E-01, 1.03652E+00, 1.15048E+00, 0032 1.28167E+00, 1.43933E+00, 1.64500E+00, 1.96000E+00, 0033 3.87000E+00 0034 }; 0035 0036 static int first_time = 1; 0037 0038 double u; 0039 double di; 0040 int index, minus; 0041 double delta, gaussian_random_value; 0042 0043 if (first_time) { 0044 srand48(seed); 0045 first_time = 0; 0046 } 0047 /* 0048 compute uniform random number between 0.0 - 0.5, and a sign with 0049 probability 1/2 of being either + or - 0050 */ 0051 u = drand48(); 0052 if (u >= 0.5) { 0053 minus = 0; 0054 u = u - 0.5; 0055 } 0056 else { 0057 minus = 1; 0058 } 0059 /* interpolate gaussian random number using table */ 0060 0061 index = (int)(di = (didu * u)); 0062 if (index == itblmax) { 0063 delta = tbl[index]; 0064 } 0065 else { 0066 di -= index; 0067 delta = tbl[index] + (tbl[index + 1] - tbl[index]) * di; 0068 } 0069 gaussian_random_value = mean + std * (minus ? -delta : delta); 0070 0071 return(gaussian_random_value); 0072 }