1 /* $Id: boost_erf.c 498248 2016-04-14 15:42:26Z boratyng $
2 * ===========================================================================
3 *
4 * This code is ported from the open source C++ library Boost, for which
5 * the following notice applies:
6 *
7 * (C) Copyright John Maddock 2006.
8 * Use, modification and distribution are subject to the Boost Software
9 * License, Version 1.0. (See below or copy at
10 * http://www.boost.org/LICENSE_1_0.txt)
11 *
12 * Boost Software License - Version 1.0 - August 17th, 2003
13 *
14 * Permission is hereby granted, free of charge, to any person or organization
15 * obtaining a copy of the software and accompanying documentation covered by
16 * this license (the "Software") to use, reproduce, display, distribute,
17 * execute, and transmit the Software, and to prepare derivative works of the
18 * Software, and to permit third-parties to whom the Software is furnished to
19 * do so, all subject to the following:
20 *
21 * The copyright notices in the Software and this entire statement, including
22 * the above license grant, this restriction and the following disclaimer,
23 * must be included in all copies of the Software, in whole or in part, and
24 * all derivative works of the Software, unless such copies or derivative
25 * works are solely in the form of machine-executable object code generated by
26 * a source language processor.
27 *
28 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
29 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
30 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
31 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
32 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
33 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
34 * DEALINGS IN THE SOFTWARE.
35 *
36 * ===========================================================================
37 *
38 * Author (editor, really): Greg Boratyn, NCBI
39 *
40 * File Description:
41 * Error function and complementary error function for normal distribution.
42 *
43 * ===========================================================================
44 */
45
46 #include "boost_erf.h"
47 #include <math.h>
48
49 #define FALSE 0
50 #define TRUE 1
51
ErfImpl(double z,int invert)52 static double ErfImpl(double z, int invert)
53 {
54 if (z < 0) {
55 if (!invert) {
56 return -ErfImpl(-z, invert);
57 }
58 else if (z < -0.5) {
59 return 2.0 - ErfImpl(-z, invert);
60 }
61 else {
62 return 1 + ErfImpl(-z, FALSE);
63 }
64 }
65
66 double result;
67
68 /* Big bunch of selection statements now to pick
69 which implementation to use,
70 try to put most likely options first: */
71 if (z < 0.5) {
72
73 /* We're going to calculate erf: */
74 if (z < 1e-10) {
75 if (z == 0) {
76 result = 0.0;
77 }
78 else {
79 static const double c = 0.003379167095512573896158903121545171688;
80 result = z * 1.125 + z * c;
81 }
82 }
83 else {
84
85 /* Maximum Deviation Found: 1.561e-17
86 Expected Error Term: 1.561e-17
87 Maximum Relative Change in Control Points: 1.155e-04
88 Max Error found at double precision = 2.961182e-17 */
89
90 static const double Y = 1.044948577880859375;
91 static const double P[] = {0.0834305892146531832907,
92 -0.338165134459360935041,
93 -0.0509990735146777432841,
94 -0.00772758345802133288487,
95 -0.000322780120964605683831};
96
97 static const double Q[] = {1.0,
98 0.455004033050794024546,
99 0.0875222600142252549554,
100 0.00858571925074406212772,
101 0.000370900071787748000569};
102
103 double zz = z * z;
104
105 double p = (((P[4] * zz + P[3]) * zz + P[2]) * zz + P[1]) * zz + P[0];
106 double q = (((Q[4] * zz + Q[3]) * zz + Q[2]) * zz + Q[1]) * zz + Q[0];
107
108 result = z * (Y + p / q);
109 }
110 }
111 else if (invert ? (z < 28.0) : (z < 5.8)) {
112
113 /* We'll be calculating erfc: */
114 invert = !invert;
115
116 if(z < 1.5) {
117
118 /* Maximum Deviation Found: 3.702e-17
119 Expected Error Term: 3.702e-17
120 Maximum Relative Change in Control Points: 2.845e-04
121 Max Error found at double precision = 4.841816e-17 */
122 static const double Y = 0.405935764312744140625;
123 static const double P[] = {-0.098090592216281240205,
124 0.178114665841120341155,
125 0.191003695796775433986,
126 0.0888900368967884466578,
127 0.0195049001251218801359,
128 0.00180424538297014223957};
129
130 static const double Q[] = {1.0,
131 1.84759070983002217845,
132 1.42628004845511324508,
133 0.578052804889902404909,
134 0.12385097467900864233,
135 0.0113385233577001411017,
136 0.337511472483094676155e-5};
137
138 double x = z - 0.5;
139 double p = ((((P[5] * x + P[4]) * x + P[3]) * x + P[2]) * x + P[1]) * x + P[0];
140 double q = ((((Q[5] * x + Q[4]) * x + Q[3]) * x + Q[2]) * x + Q[1]) * x + Q[0];
141
142 result = Y + p / q;
143 result *= expl(-z * z) / z;
144 }
145 else if(z < 2.5) {
146
147 /* Max Error found at double precision = 6.599585e-18
148 Maximum Deviation Found: 3.909e-18
149 Expected Error Term: 3.909e-18
150 Maximum Relative Change in Control Points: 9.886e-05 */
151 static const double Y = 0.50672817230224609375;
152 static const double P[] = {-0.0243500476207698441272,
153 0.0386540375035707201728,
154 0.04394818964209516296,
155 0.0175679436311802092299,
156 0.00323962406290842133584,
157 0.000235839115596880717416};
158
159 static const double Q[] = {1.0,
160 1.53991494948552447182,
161 0.982403709157920235114,
162 0.325732924782444448493,
163 0.0563921837420478160373,
164 0.00410369723978904575884};
165
166 double x = z - 1.5;
167 double p = ((((P[5] * x + P[4]) * x + P[3]) * x + P[2]) * x + P[1]) * x + P[0];
168 double q = ((((Q[5] * x + Q[4]) * x + Q[3]) * x + Q[2]) * x + Q[1]) * x + Q[0];
169
170 result = Y + p / q;
171 result *= expl(-z * z) / z;
172 }
173 else if(z < 4.5) {
174
175 /* Maximum Deviation Found: 1.512e-17
176 Expected Error Term: 1.512e-17
177 Maximum Relative Change in Control Points: 2.222e-04
178 Max Error found at double precision = 2.062515e-17 */
179 static const double Y = 0.5405750274658203125;
180 static const double P[] = {0.00295276716530971662634,
181 0.0137384425896355332126,
182 0.00840807615555585383007,
183 0.00212825620914618649141,
184 0.000250269961544794627958,
185 0.113212406648847561139e-4};
186
187 static const double Q[] = {1.0,
188 1.04217814166938418171,
189 0.442597659481563127003,
190 0.0958492726301061423444,
191 0.0105982906484876531489,
192 0.000479411269521714493907};
193
194 double x = z - 3.5;
195 double p = ((((P[5] * x + P[4]) * x + P[3]) * x + P[2]) * x + P[1]) * x + P[0];
196 double q = ((((Q[5] * x + Q[4]) * x + Q[3]) * x + Q[2]) * x + Q[1]) * x + Q[0];
197
198 result = Y + p / q;
199 result *= expl(-z * z) / z;
200 }
201 else {
202
203 /* Max Error found at double precision = 2.997958e-17
204 Maximum Deviation Found: 2.860e-17
205 Expected Error Term: 2.859e-17
206 Maximum Relative Change in Control Points: 1.357e-05 */
207 static const double Y = 0.5579090118408203125;
208 static const double P[] = {0.00628057170626964891937,
209 0.0175389834052493308818,
210 -0.212652252872804219852,
211 -0.687717681153649930619,
212 -2.5518551727311523996,
213 -3.22729451764143718517,
214 -2.8175401114513378771};
215
216 static const double Q[] = {1.0,
217 2.79257750980575282228,
218 11.0567237927800161565,
219 15.930646027911794143,
220 22.9367376522880577224,
221 13.5064170191802889145,
222 5.48409182238641741584};
223
224 double x = 1.0 / z;
225 double p = (((((P[6] * x + P[5]) * x + P[4]) * x + P[3]) * x + P[2]) * x + P[1]) * x + P[0];
226 double q = (((((Q[6] * x + Q[5]) * x + Q[4]) * x + Q[3]) * x + Q[2]) * x + Q[1]) * x + Q[0];
227
228 result = Y + p / q;
229 result *= expl(-z * z) / z;
230 }
231 }
232 else {
233
234 /* Any value of z larger than 28 will underflow to zero: */
235 result = 0;
236 invert = !invert;
237 }
238
239 if(invert) {
240 result = 1 - result;
241 }
242
243 return result;
244 }
245
246
Erf(double z)247 double Erf(double z)
248 {
249 return ErfImpl(z, FALSE);
250 }
251
ErfC(double z)252 double ErfC(double z)
253 {
254 return ErfImpl(z, TRUE);
255 }
256
257