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)8static __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)18static __inline long int lrint(double x){ 19 return (long)rint(x); 20 } 21 #endif fmax(double x,double y)22static __inline double fmax(double x, double y){ 23 return (x > y) ? x : y; 24 } fmin(double x,double y)25static __inline double fmin(double x, double y){ 26 return (x < y) ? x : y; 27 } log2(double x)28static __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)35double 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