OpenMM
 All Classes Namespaces Functions Variables Enumerations Enumerator Friends
MSVC_erfc.h
1 #ifndef OPENMM_MSVC_ERFC_H_
2 #define OPENMM_MSVC_ERFC_H_
3 
4 /*
5  * At least up to version 8 (VC++ 2005), Microsoft does not support the
6  * standard C99 erf() and erfc() functions. For now we're including these
7  * definitions for an MSVC compilation; if these are added later then
8  * the #ifdef below should change to compare _MSC_VER with a particular
9  * version level.
10  */
11 
12 #ifdef _MSC_VER
13 
14 
15 /***************************
16 * erf.cpp
17 * author: Steve Strand
18 * written: 29-Jan-04
19 ***************************/
20 
21 #include <cmath>
22 
23 static const double rel_error= 1E-12; //calculate 12 significant figures
24 //you can adjust rel_error to trade off between accuracy and speed
25 //but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)
26 
27 static double erfc(double x);
28 
29 static double erf(double x)
30 //erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
31 // = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
32 // = 1-erfc(x)
33 {
34  static const double two_sqrtpi= 1.128379167095512574; // 2/sqrt(pi)
35  if (fabs(x) > 2.2) {
36  return 1.0 - erfc(x); //use continued fraction when fabs(x) > 2.2
37  }
38  double sum= x, term= x, xsqr= x*x;
39  int j= 1;
40  do {
41  term*= xsqr/j;
42  sum-= term/(2*j+1);
43  ++j;
44  term*= xsqr/j;
45  sum+= term/(2*j+1);
46  ++j;
47  } while (fabs(term)/sum > rel_error);
48  return two_sqrtpi*sum;
49 }
50 
51 
52 static double erfc(double x)
53 //erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
54 // = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
55 // = 1-erf(x)
56 //expression inside [] is a continued fraction so '+' means add to denominator only
57 {
58  static const double one_sqrtpi= 0.564189583547756287; // 1/sqrt(pi)
59  if (fabs(x) < 2.2) {
60  return 1.0 - erf(x); //use series when fabs(x) < 2.2
61  }
62  // Don't look for x==0 here!
63  if (x < 0) { //continued fraction only valid for x>0
64  return 2.0 - erfc(-x);
65  }
66  double a=1, b=x; //last two convergent numerators
67  double c=x, d=x*x+0.5; //last two convergent denominators
68  double q1, q2= b/d; //last two convergents (a/c and b/d)
69  double n= 1.0, t;
70  do {
71  t= a*n+b*x;
72  a= b;
73  b= t;
74  t= c*n+d*x;
75  c= d;
76  d= t;
77  n+= 0.5;
78  q1= q2;
79  q2= b/d;
80  } while (fabs(q1-q2)/q2 > rel_error);
81  return one_sqrtpi*exp(-x*x)*q2;
82 }
83 
84 #endif // _MSC_VER
85 
86 #endif // OPENMM_MSVC_ERFC_H_