1 #include <math.h>
2 #ifdef MS_WIN32
3 #include "malloc.h"
4 typedef int int32_t;
5 typedef long long int64_t;
6 /* Taken from http://siliconandlithium.blogspot.com/2014/05/msvc-c99-mathh-header.html */
7 #define isnormal(x) ((_fpclass(x) == _FPCLASS_NN) || (_fpclass(x) == _FPCLASS_PN))
rint(double x)8 static __inline double rint(double x){
9     const double two_to_52 = 4.5035996273704960e+15;
10     double fa = fabs(x);
11     if(fa >= two_to_52){
12         return x;
13     } else{
14         return copysign(two_to_52 + fa - two_to_52, x);
15     }
16 }
17 #if _MSC_VER < 1928
lrint(double x)18 static __inline long int lrint(double x){
19     return (long)rint(x);
20 }
21 #endif
fmax(double x,double y)22 static __inline double fmax(double x, double y){
23     return (x > y) ? x : y;
24 }
fmin(double x,double y)25 static __inline double fmin(double x, double y){
26     return (x < y) ? x : y;
27 }
log2(double x)28 static __inline double log2(double x) {
29     return log(x) * M_LOG2E;
30 }
31 
32 /* adapted from http://www.johndcook.com/blog/cpp_erf/
33    code is under public domain license */
34 
erf(double x)35 double erf(double x)
36 {
37     /* constants */
38     double a1 =  0.254829592;
39     double a2 = -0.284496736;
40     double a3 =  1.421413741;
41     double a4 = -1.453152027;
42     double a5 =  1.061405429;
43     double p  =  0.3275911;
44     double t;
45     double y;
46 
47     /* Save the sign of x */
48     int sign = 1;
49     if (x < 0)
50         sign = -1;
51     x = fabs(x);
52 
53     /* A&S formula 7.1.26 */
54     t = 1.0/(1.0 + p*x);
55     y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);
56 
57     return sign*y;
58 }
59 
60 #else
61 #include <stdint.h>
62 #if !defined(__FreeBSD__)
63 #  include "alloca.h"
64 #else
65 #  include <stdlib.h>
66 #endif
67 #include <math.h>
68 #endif
69