You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

MSVC_erfc.h 2.5KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687
  1. #ifndef LEPTON_MSVC_ERFC_H_
  2. #define LEPTON_MSVC_ERFC_H_
  3. /*
  4. * Up to version 11 (VC++ 2012), Microsoft does not support the
  5. * standard C99 erf() and erfc() functions so we have to fake them here.
  6. * These were added in version 12 (VC++ 2013), which sets _MSC_VER=1800
  7. * (VC11 has _MSC_VER=1700).
  8. */
  9. #if defined(_MSC_VER)
  10. #define M_PI 3.14159265358979323846264338327950288
  11. #if _MSC_VER <= 1700 // 1700 is VC11, 1800 is VC12
  12. /***************************
  13. * erf.cpp
  14. * author: Steve Strand
  15. * written: 29-Jan-04
  16. ***************************/
  17. #include <cmath>
  18. static const double rel_error= 1E-12; //calculate 12 significant figures
  19. //you can adjust rel_error to trade off between accuracy and speed
  20. //but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)
  21. static double erfc(double x);
  22. static double erf(double x)
  23. //erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
  24. // = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
  25. // = 1-erfc(x)
  26. {
  27. static const double two_sqrtpi= 1.128379167095512574; // 2/sqrt(pi)
  28. if (fabs(x) > 2.2) {
  29. return 1.0 - erfc(x); //use continued fraction when fabs(x) > 2.2
  30. }
  31. double sum= x, term= x, xsqr= x*x;
  32. int j= 1;
  33. do {
  34. term*= xsqr/j;
  35. sum-= term/(2*j+1);
  36. ++j;
  37. term*= xsqr/j;
  38. sum+= term/(2*j+1);
  39. ++j;
  40. } while (fabs(term)/sum > rel_error);
  41. return two_sqrtpi*sum;
  42. }
  43. static double erfc(double x)
  44. //erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
  45. // = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
  46. // = 1-erf(x)
  47. //expression inside [] is a continued fraction so '+' means add to denominator only
  48. {
  49. static const double one_sqrtpi= 0.564189583547756287; // 1/sqrt(pi)
  50. if (fabs(x) < 2.2) {
  51. return 1.0 - erf(x); //use series when fabs(x) < 2.2
  52. }
  53. // Don't look for x==0 here!
  54. if (x < 0) { //continued fraction only valid for x>0
  55. return 2.0 - erfc(-x);
  56. }
  57. double a=1, b=x; //last two convergent numerators
  58. double c=x, d=x*x+0.5; //last two convergent denominators
  59. double q1, q2= b/d; //last two convergents (a/c and b/d)
  60. double n= 1.0, t;
  61. do {
  62. t= a*n+b*x;
  63. a= b;
  64. b= t;
  65. t= c*n+d*x;
  66. c= d;
  67. d= t;
  68. n+= 0.5;
  69. q1= q2;
  70. q2= b/d;
  71. } while (fabs(q1-q2)/q2 > rel_error);
  72. return one_sqrtpi*exp(-x*x)*q2;
  73. }
  74. #endif // _MSC_VER <= 1700
  75. #endif // _MSC_VER
  76. #endif // LEPTON_MSVC_ERFC_H_