1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024.
6 //
7 // This file is part of Force Field X.
8 //
9 // Force Field X is free software; you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License version 3 as published by
11 // the Free Software Foundation.
12 //
13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16 // details.
17 //
18 // You should have received a copy of the GNU General Public License along with
19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20 // Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // Linking this library statically or dynamically with other modules is making a
23 // combined work based on this library. Thus, the terms and conditions of the
24 // GNU General Public License cover the whole combination.
25 //
26 // As a special exception, the copyright holders of this library give you
27 // permission to link this library with independent modules to produce an
28 // executable, regardless of the license terms of these independent modules, and
29 // to copy and distribute the resulting executable under terms of your choice,
30 // provided that you also meet, for each linked independent module, the terms
31 // and conditions of the license of that module. An independent module is a
32 // module which is not derived from or based on this library. If you modify this
33 // library, you may extend this exception to your version of the library, but
34 // you are not obligated to do so. If you do not wish to do so, delete this
35 // exception statement from your version.
36 //
37 // ******************************************************************************
38 package ffx.numerics.special;
39
40 import static org.apache.commons.math3.util.FastMath.PI;
41 import static org.apache.commons.math3.util.FastMath.abs;
42 import static org.apache.commons.math3.util.FastMath.exp;
43 import static org.apache.commons.math3.util.FastMath.floor;
44 import static org.apache.commons.math3.util.FastMath.sqrt;
45
46 /**
47 * Static methods to evaluate erf(x) and erfc(x) for a real argument x. Rational functions are used
48 * that approximate erf(x) and erfc(x) to machine precision (approximately 15 decimal digits).
49 *
50 * <p>Adapted from an original program written by W. J. Cody, Mathematics and Computer Science
51 * Division, Argonne National Laboratory, Argonne, IL 60439
52 *
53 * <p>
54 *
55 * @author Michael J. Schnieders
56 * @since 1.0
57 */
58 public class Erf {
59
60 /**
61 * Mathematical and machine-dependent constants.
62 */
63 private static final double sqrtPI = 1.0 / sqrt(PI);
64
65 private static final double oneSixteenth = 1.0 / 16.0;
66 private static final double thresh = 0.46875;
67 /**
68 * The xSmall argument below which erf(x) may be represented by 2*x/sqrt(pi) and above which x*x won't underflow.
69 * <p>
70 * A conservative value is the largest machine number X such that 1.0 + X = 1.0 to machine precision.
71 */
72 private static final double xSmall = 1.11e-16;
73 /**
74 * xBig is the largest argument acceptable for erfc.
75 * <p>
76 * Solution to the equation:
77 * <p>
78 * W(x) * (1-0.5/x**2) = xMin, where
79 * <p>
80 * W(x) = exp(-x*x)/[x*sqrt(pi)]
81 */
82 private static final double xBig = 26.543;
83
84 private Erf() {
85 }
86
87 /**
88 * Evaluates erf(x) for a real argument x.
89 *
90 * @param arg the value to evaluate erf at.
91 * @return erf of the argument.
92 * @since 1.0
93 */
94 public static double erf(double arg) {
95 return erfCore(arg, false);
96 }
97
98 /**
99 * Evaluate erfc(x) for a real argument x.
100 *
101 * @param arg the value to evaluate erfc at.
102 * @return erfc of the argument.
103 * @since 1.0
104 */
105 public static double erfc(double arg) {
106 return erfCore(arg, true);
107 }
108
109 /**
110 * Evaluates erf(x) or erfc(x) for a real argument x. When called with mode = false, erf is
111 * returned, while with mode = true, erfc is returned.
112 *
113 * @param x the value to evaluate erf or erfc at.
114 * @param mode if mode is true, evaluate erfc, otherwise evaluate erf.
115 * @return if (!mode) erf(arg), else erfc(arg)
116 * @since 1.0
117 */
118 private static double erfCore(double x, boolean mode) {
119
120 // Store the argument and its absolute value.
121 final double y = abs(x);
122 double result = 0.0;
123
124 // Evaluate error function for |x| less than 0.46875.
125 if (y <= thresh) {
126 double ysq = 0.0;
127 if (y > xSmall) {
128 ysq = y * y;
129 }
130 double xNum = 1.85777706184603153e-1 * ysq;
131 double xDen = ysq;
132 xNum = (xNum + 3.16112374387056560e0) * ysq;
133 xDen = (xDen + 2.36012909523441209e1) * ysq;
134 xNum = (xNum + 1.13864154151050156e2) * ysq;
135 xDen = (xDen + 2.44024637934444173e2) * ysq;
136 xNum = (xNum + 3.77485237685302021e2) * ysq;
137 xDen = (xDen + 1.28261652607737228e3) * ysq;
138 result = x * (xNum + 3.20937758913846947e3) / (xDen + 2.84423683343917062e3);
139 if (mode) {
140 result = 1.0 - result;
141 }
142 } else if (y <= 4.0) {
143 // Get complementary error function for 0.46875 <= |x| <= 4.0.
144 double xNum = 2.15311535474403846e-8 * y;
145 double xDen = y;
146 xNum = (xNum + 5.64188496988670089e-1) * y;
147 xDen = (xDen + 1.57449261107098347e1) * y;
148 xNum = (xNum + 8.88314979438837594e0) * y;
149 xDen = (xDen + 1.17693950891312499e2) * y;
150 xNum = (xNum + 6.61191906371416295e1) * y;
151 xDen = (xDen + 5.37181101862009858e2) * y;
152 xNum = (xNum + 2.98635138197400131e2) * y;
153 xDen = (xDen + 1.62138957456669019e3) * y;
154 xNum = (xNum + 8.81952221241769090e2) * y;
155 xDen = (xDen + 3.29079923573345963e3) * y;
156 xNum = (xNum + 1.71204761263407058e3) * y;
157 xDen = (xDen + 4.36261909014324716e3) * y;
158 xNum = (xNum + 2.05107837782607147e3) * y;
159 xDen = (xDen + 3.43936767414372164e3) * y;
160 result = (xNum + 1.23033935479799725e3) / (xDen + 1.23033935480374942e3);
161 double ysq = floor(16.0 * y) * oneSixteenth;
162 double del = (y - ysq) * (y + ysq);
163 result = exp(-ysq * ysq - del) * result;
164 if (!mode) {
165 result = 1.0 - result;
166 if (x < 0.0) {
167 result = -result;
168 }
169 } else if (x < 0.0) {
170 result = 2.0 - result;
171 }
172 } else {
173 // Get complementary error function for |x| greater than 4.0.
174 if (y < xBig) {
175 double iy = 1.0 / y;
176 double ysq = iy * iy;
177 double xNum = 1.63153871373020978e-2 * ysq;
178 double xDen = ysq;
179 xNum = (xNum + 3.05326634961232344e-1) * ysq;
180 xDen = (xDen + 2.56852019228982242e0) * ysq;
181 xNum = (xNum + 3.60344899949804439e-1) * ysq;
182 xDen = (xDen + 1.87295284992346047e0) * ysq;
183 xNum = (xNum + 1.25781726111229246e-1) * ysq;
184 xDen = (xDen + 5.27905102951428412e-1) * ysq;
185 xNum = (xNum + 1.60837851487422766e-2) * ysq;
186 xDen = (xDen + 6.05183413124413191e-2) * ysq;
187 result = ysq * (xNum + 6.58749161529837803e-4) / (xDen + 2.33520497626869185e-3);
188 result = (sqrtPI - result) * iy;
189 ysq = floor(16.0 * y) * oneSixteenth;
190 double del = (y - ysq) * (y + ysq);
191 result = exp(-ysq * ysq - del) * result;
192 }
193 if (!mode) {
194 result = 1.0 - result;
195 if (x < 0.0) {
196 result = -result;
197 }
198 } else {
199 if (x < 0.0) {
200 result = 2.0 - result;
201 }
202 }
203 }
204 return result;
205 }
206 }