1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 // (C) Copyright John Maddock 2006.
19 // Use, modification and distribution are subject to the
20 // Boost Software License, Version 1.0. (See accompanying file
21 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
22
23 package org.apache.commons.numbers.gamma;
24
25 import org.apache.commons.numbers.core.DD;
26
27 /**
28 * Implementation of the <a href="https://mathworld.wolfram.com/Erf.html">error function</a> and
29 * its inverse.
30 *
31 * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
32 * {@code c++} implementation {@code <boost/math/special_functions/erf.hpp>}.
33 * The erf/erfc functions and their inverses are copyright John Maddock 2006 and subject to
34 * the Boost Software License.
35 *
36 * <p>Additions made to support the erfcx function are original work under the Apache software
37 * license.
38 *
39 * @see
40 * <a href="https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_erf/error_function.html">
41 * Boost C++ Error functions</a>
42 */
43 final class BoostErf {
44 /** 1 / sqrt(pi). Used for the scaled complementary error function erfcx. */
45 private static final double ONE_OVER_ROOT_PI = 0.5641895835477562869480794515607725858;
46 /** Threshold for the scaled complementary error function erfcx
47 * where the approximation {@code (1 / sqrt(pi)) / x} can be used. */
48 private static final double ERFCX_APPROX = 6.71e7;
49 /** Threshold for the erf implementation for |x| where the computation
50 * uses {@code erf(x)}; otherwise {@code erfc(x)} is computed. The final result is
51 * achieved by suitable application of symmetry. */
52 private static final double COMPUTE_ERF = 0.5;
53 /** Threshold for the scaled complementary error function erfcx for negative x
54 * where {@code 2 * exp(x*x)} will overflow. Value is 26.62873571375149. */
55 private static final double ERFCX_NEG_X_MAX = Math.sqrt(Math.log(Double.MAX_VALUE / 2));
56 /** Threshold for the scaled complementary error function erfcx for x
57 * where {@code exp(x*x) == 1; x <= t}. Value is (1 + 5/16) * 2^-27 = 9.778887033462524E-9.
58 * <p>Note: This is used for performance. If set to 0 then the result is computed
59 * using expm1(x*x) with the same final result. */
60 private static final double EXP_XX_1 = 0x1.5p-27;
61
62 /** Private constructor. */
63 private BoostErf() {
64 // intentionally empty.
65 }
66
67 // Code ported from Boost 1.77.0
68 //
69 // boost/math/special_functions/erf.hpp
70 // boost/math/special_functions/detail/erf_inv.hpp
71 //
72 // Original code comments, including measured deviations, are preserved.
73 //
74 // Changes to the Boost implementation:
75 // - Update method names to replace underscores with camel case
76 // - Explicitly inline the polynomial function evaluation
77 // using Horner's method (https://en.wikipedia.org/wiki/Horner%27s_method)
78 // - Support odd function for f(0.0) = -f(-0.0)
79 // - Support the scaled complementary error function erfcx
80 // Erf:
81 // - Change extended precision z*z to compute the square round-off
82 // using Dekker's method
83 // - Change extended precision exp(-z*z) to compute using a
84 // round-off addition to the standard exp result (see NUMBERS-177)
85 // - Change the erf threshold for z when erf(z)=1 from
86 // z > 5.8f to z > 5.930664
87 // - Change the erfc threshold for z when erfc(z)=0 from
88 // z < 28 to z < 27.3
89 // - Change rational function approximation for z > 4 to a function
90 // suitable for erfcx (see NUMBERS-177)
91 // Inverse erf:
92 // - Change inverse erf edge case detection to include NaN
93 // - Change edge case detection for integer z
94 //
95 // Note:
96 // Constants using the 'f' suffix are machine
97 // representable as a float, e.g.
98 // assert 0.0891314744949340820313f == 0.0891314744949340820313;
99 // The values are unchanged from the Boost reference.
100
101 /**
102 * Returns the complementary error function.
103 *
104 * @param x the value.
105 * @return the complementary error function.
106 */
107 static double erfc(double x) {
108 return erfImp(x, true, false);
109 }
110
111 /**
112 * Returns the error function.
113 *
114 * @param x the value.
115 * @return the error function.
116 */
117 static double erf(double x) {
118 return erfImp(x, false, false);
119 }
120
121 /**
122 * 53-bit implementation for the error function.
123 *
124 * <p>Note: The {@code scaled} flag only applies when
125 * {@code z >= 0.5} and {@code invert == true}.
126 * This functionality is used to compute erfcx(z) for positive z.
127 *
128 * @param z Point to evaluate
129 * @param invert true to invert the result (for the complementary error function)
130 * @param scaled true to compute the scaled complementary error function
131 * @return the error function result
132 */
133 private static double erfImp(double z, boolean invert, boolean scaled) {
134 if (Double.isNaN(z)) {
135 return Double.NaN;
136 }
137
138 if (z < 0) {
139 // Here the scaled flag is ignored.
140 if (!invert) {
141 return -erfImp(-z, invert, false);
142 } else if (z < -0.5) {
143 return 2 - erfImp(-z, invert, false);
144 } else {
145 return 1 + erfImp(-z, false, false);
146 }
147 }
148
149 double result;
150
151 //
152 // Big bunch of selection statements now to pick
153 // which implementation to use,
154 // try to put most likely options first:
155 //
156 if (z < COMPUTE_ERF) {
157 //
158 // We're going to calculate erf:
159 //
160 // Here the scaled flag is ignored.
161 if (z < 1e-10) {
162 if (z == 0) {
163 result = z;
164 } else {
165 final double c = 0.003379167095512573896158903121545171688;
166 result = z * 1.125f + z * c;
167 }
168 } else {
169 // Maximum Deviation Found: 1.561e-17
170 // Expected Error Term: 1.561e-17
171 // Maximum Relative Change in Control Points: 1.155e-04
172 // Max Error found at double precision = 2.961182e-17
173
174 final double Y = 1.044948577880859375f;
175 final double zz = z * z;
176 double P;
177 P = -0.000322780120964605683831;
178 P = -0.00772758345802133288487 + P * zz;
179 P = -0.0509990735146777432841 + P * zz;
180 P = -0.338165134459360935041 + P * zz;
181 P = 0.0834305892146531832907 + P * zz;
182 double Q;
183 Q = 0.000370900071787748000569;
184 Q = 0.00858571925074406212772 + Q * zz;
185 Q = 0.0875222600142252549554 + Q * zz;
186 Q = 0.455004033050794024546 + Q * zz;
187 Q = 1.0 + Q * zz;
188 result = z * (Y + P / Q);
189 }
190 // Note: Boost threshold of 5.8f has been raised to approximately 5.93 (6073 / 1024);
191 // threshold of 28 has been lowered to approximately 27.3 (6989/256) where exp(-z*z) = 0.
192 } else if (scaled || (invert ? (z < 27.300781f) : (z < 5.9306640625f))) {
193 //
194 // We'll be calculating erfc:
195 //
196 // Here the scaled flag is used.
197 invert = !invert;
198 if (z < 1.5f) {
199 // Maximum Deviation Found: 3.702e-17
200 // Expected Error Term: 3.702e-17
201 // Maximum Relative Change in Control Points: 2.845e-04
202 // Max Error found at double precision = 4.841816e-17
203 final double Y = 0.405935764312744140625f;
204 final double zm = z - 0.5;
205 double P;
206 P = 0.00180424538297014223957;
207 P = 0.0195049001251218801359 + P * zm;
208 P = 0.0888900368967884466578 + P * zm;
209 P = 0.191003695796775433986 + P * zm;
210 P = 0.178114665841120341155 + P * zm;
211 P = -0.098090592216281240205 + P * zm;
212 double Q;
213 Q = 0.337511472483094676155e-5;
214 Q = 0.0113385233577001411017 + Q * zm;
215 Q = 0.12385097467900864233 + Q * zm;
216 Q = 0.578052804889902404909 + Q * zm;
217 Q = 1.42628004845511324508 + Q * zm;
218 Q = 1.84759070983002217845 + Q * zm;
219 Q = 1.0 + Q * zm;
220 result = Y + P / Q;
221 if (scaled) {
222 result /= z;
223 } else {
224 result *= expmxx(z) / z;
225 }
226 } else if (z < 2.5f) {
227 // Max Error found at double precision = 6.599585e-18
228 // Maximum Deviation Found: 3.909e-18
229 // Expected Error Term: 3.909e-18
230 // Maximum Relative Change in Control Points: 9.886e-05
231 final double Y = 0.50672817230224609375f;
232 final double zm = z - 1.5;
233 double P;
234 P = 0.000235839115596880717416;
235 P = 0.00323962406290842133584 + P * zm;
236 P = 0.0175679436311802092299 + P * zm;
237 P = 0.04394818964209516296 + P * zm;
238 P = 0.0386540375035707201728 + P * zm;
239 P = -0.0243500476207698441272 + P * zm;
240 double Q;
241 Q = 0.00410369723978904575884;
242 Q = 0.0563921837420478160373 + Q * zm;
243 Q = 0.325732924782444448493 + Q * zm;
244 Q = 0.982403709157920235114 + Q * zm;
245 Q = 1.53991494948552447182 + Q * zm;
246 Q = 1.0 + Q * zm;
247 result = Y + P / Q;
248 if (scaled) {
249 result /= z;
250 } else {
251 result *= expmxx(z) / z;
252 }
253 // Lowered Boost threshold from 4.5 to 4.0 as this is the limit
254 // for the Cody erfc approximation
255 } else if (z < 4.0f) {
256 // Maximum Deviation Found: 1.512e-17
257 // Expected Error Term: 1.512e-17
258 // Maximum Relative Change in Control Points: 2.222e-04
259 // Max Error found at double precision = 2.062515e-17
260 final double Y = 0.5405750274658203125f;
261 final double zm = z - 3.5;
262 double P;
263 P = 0.113212406648847561139e-4;
264 P = 0.000250269961544794627958 + P * zm;
265 P = 0.00212825620914618649141 + P * zm;
266 P = 0.00840807615555585383007 + P * zm;
267 P = 0.0137384425896355332126 + P * zm;
268 P = 0.00295276716530971662634 + P * zm;
269 double Q;
270 Q = 0.000479411269521714493907;
271 Q = 0.0105982906484876531489 + Q * zm;
272 Q = 0.0958492726301061423444 + Q * zm;
273 Q = 0.442597659481563127003 + Q * zm;
274 Q = 1.04217814166938418171 + Q * zm;
275 Q = 1.0 + Q * zm;
276 result = Y + P / Q;
277 if (scaled) {
278 result /= z;
279 } else {
280 result *= expmxx(z) / z;
281 }
282 } else {
283 // Rational function approximation for erfc(x > 4.0)
284 //
285 // This approximation is not the Boost implementation.
286 // The Boost function is suitable for [4.5 < z < 28].
287 //
288 // This function is suitable for erfcx(z) as it asymptotes
289 // to (1 / sqrt(pi)) / z at large z.
290 //
291 // Taken from "Rational Chebyshev approximations for the error function"
292 // by W. J. Cody, Math. Comp., 1969, PP. 631-638.
293 //
294 // See NUMBERS-177.
295
296 final double izz = 1 / (z * z);
297 double p;
298 p = 1.63153871373020978498e-2;
299 p = 3.05326634961232344035e-1 + p * izz;
300 p = 3.60344899949804439429e-1 + p * izz;
301 p = 1.25781726111229246204e-1 + p * izz;
302 p = 1.60837851487422766278e-2 + p * izz;
303 p = 6.58749161529837803157e-4 + p * izz;
304 double q;
305 q = 1;
306 q = 2.56852019228982242072e00 + q * izz;
307 q = 1.87295284992346047209e00 + q * izz;
308 q = 5.27905102951428412248e-1 + q * izz;
309 q = 6.05183413124413191178e-2 + q * izz;
310 q = 2.33520497626869185443e-3 + q * izz;
311
312 result = izz * p / q;
313 result = (ONE_OVER_ROOT_PI - result) / z;
314
315 if (!scaled) {
316 // In contrast to smaller z: avoid using result *= expmxx(z) / z.
317 // exp(-z*z) can be sub-normal so
318 // multiply by any sub-normal after divide by z
319 result *= expmxx(z);
320 }
321 }
322 } else {
323 //
324 // Any value of z larger than 27.3 will underflow to zero:
325 //
326 result = 0;
327 invert = !invert;
328 }
329
330 if (invert) {
331 // Note: If 0.5 <= z < 28 and the scaled flag is true then
332 // invert will have been flipped to false and the
333 // the result is unchanged as erfcx(z)
334 result = 1 - result;
335 }
336
337 return result;
338 }
339
340 /**
341 * Returns the scaled complementary error function.
342 * <pre>
343 * erfcx(x) = exp(x^2) * erfc(x)
344 * </pre>
345 *
346 * @param x the value.
347 * @return the scaled complementary error function.
348 */
349 static double erfcx(double x) {
350 if (Double.isNaN(x)) {
351 return Double.NaN;
352 }
353
354 // For |z| < 0.5 erfc is computed using erf
355 final double ax = Math.abs(x);
356 if (ax < COMPUTE_ERF) {
357 // Use the erf(x) result.
358 // (1 - erf(x)) * exp(x*x)
359
360 final double erfx = erf(x);
361 if (ax < EXP_XX_1) {
362 // No exponential required
363 return 1 - erfx;
364 }
365
366 // exp(x*x) - exp(x*x) * erf(x)
367 // Avoid use of exp(x*x) with expm1:
368 // exp(x*x) - 1 - (erf(x) * (exp(x*x) - 1)) - erf(x) + 1
369
370 // Sum small to large: |erf(x)| > expm1(x*x)
371 // -erf(x) * expm1(x*x) + expm1(x*x) - erf(x) + 1
372 // Negative x: erf(x) < 0, summed terms are positive, no cancellation occurs.
373 // Positive x: erf(x) > 0 so cancellation can occur.
374 // When terms are ordered by absolute magnitude the magnitude of the next term
375 // is above the round-off from adding the previous term to the sum. Thus
376 // cancellation is negligible compared to errors in the largest computed term (erf(x)).
377
378 final double em1 = Math.expm1(x * x);
379 return -erfx * em1 + em1 - erfx + 1;
380 }
381
382 // Handle negative arguments
383 if (x < 0) {
384 // erfcx(x) = 2*exp(x*x) - erfcx(-x)
385
386 if (x < -ERFCX_NEG_X_MAX) {
387 // Overflow
388 return Double.POSITIVE_INFINITY;
389 }
390
391 final double e = expxx(x);
392 return e - erfImp(-x, true, true) + e;
393 }
394
395 // Approximation for large positive x
396 if (x > ERFCX_APPROX) {
397 return ONE_OVER_ROOT_PI / x;
398 }
399
400 // Compute erfc scaled
401 return erfImp(x, true, true);
402 }
403
404 /**
405 * Returns the inverse complementary error function.
406 *
407 * @param z Value (in {@code [0, 2]}).
408 * @return t such that {@code z = erfc(t)}
409 */
410 static double erfcInv(double z) {
411 //
412 // Begin by testing for domain errors, and other special cases:
413 //
414 if (z < 0 || z > 2 || Double.isNaN(z)) {
415 // Argument outside range [0,2] in inverse erfc function
416 return Double.NaN;
417 }
418 // Domain bounds must be detected as the implementation computes NaN.
419 // (log(q=0) creates infinity and the rational number is
420 // infinity / infinity)
421 if (z == (int) z) {
422 // z return
423 // 2 -inf
424 // 1 0
425 // 0 inf
426 return z == 1 ? 0 : (1 - z) * Double.POSITIVE_INFINITY;
427 }
428
429 //
430 // Normalise the input, so it's in the range [0,1], we will
431 // negate the result if z is outside that range. This is a simple
432 // application of the erfc reflection formula: erfc(-z) = 2 - erfc(z)
433 //
434 final double p;
435 final double q;
436 final double s;
437 if (z > 1) {
438 q = 2 - z;
439 p = 1 - q;
440 s = -1;
441 } else {
442 p = 1 - z;
443 q = z;
444 s = 1;
445 }
446
447 //
448 // And get the result, negating where required:
449 //
450 return s * erfInvImp(p, q);
451 }
452
453 /**
454 * Returns the inverse error function.
455 *
456 * @param z Value (in {@code [-1, 1]}).
457 * @return t such that {@code z = erf(t)}
458 */
459 static double erfInv(double z) {
460 //
461 // Begin by testing for domain errors, and other special cases:
462 //
463 if (z < -1 || z > 1 || Double.isNaN(z)) {
464 // Argument outside range [-1, 1] in inverse erf function
465 return Double.NaN;
466 }
467 // Domain bounds must be detected as the implementation computes NaN.
468 // (log(q=0) creates infinity and the rational number is
469 // infinity / infinity)
470 if (z == (int) z) {
471 // z return
472 // -1 -inf
473 // -0 -0
474 // 0 0
475 // 1 inf
476 return z == 0 ? z : z * Double.POSITIVE_INFINITY;
477 }
478
479 //
480 // Normalise the input, so it's in the range [0,1], we will
481 // negate the result if z is outside that range. This is a simple
482 // application of the erf reflection formula: erf(-z) = -erf(z)
483 //
484 final double p;
485 final double q;
486 final double s;
487 if (z < 0) {
488 p = -z;
489 q = 1 - p;
490 s = -1;
491 } else {
492 p = z;
493 q = 1 - z;
494 s = 1;
495 }
496 //
497 // And get the result, negating where required:
498 //
499 return s * erfInvImp(p, q);
500 }
501
502 /**
503 * Common implementation for inverse erf and erfc functions.
504 *
505 * @param p P-value
506 * @param q Q-value (1-p)
507 * @return the inverse
508 */
509 private static double erfInvImp(double p, double q) {
510 final double result;
511
512 if (p <= 0.5) {
513 //
514 // Evaluate inverse erf using the rational approximation:
515 //
516 // x = p(p+10)(Y+R(p))
517 //
518 // Where Y is a constant, and R(p) is optimised for a low
519 // absolute error compared to |Y|.
520 //
521 // double: Max error found: 2.001849e-18
522 // long double: Max error found: 1.017064e-20
523 // Maximum Deviation Found (actual error term at infinite precision) 8.030e-21
524 //
525 final float Y = 0.0891314744949340820313f;
526 double P;
527 P = -0.00538772965071242932965;
528 P = 0.00822687874676915743155 + P * p;
529 P = 0.0219878681111168899165 + P * p;
530 P = -0.0365637971411762664006 + P * p;
531 P = -0.0126926147662974029034 + P * p;
532 P = 0.0334806625409744615033 + P * p;
533 P = -0.00836874819741736770379 + P * p;
534 P = -0.000508781949658280665617 + P * p;
535 double Q;
536 Q = 0.000886216390456424707504;
537 Q = -0.00233393759374190016776 + Q * p;
538 Q = 0.0795283687341571680018 + Q * p;
539 Q = -0.0527396382340099713954 + Q * p;
540 Q = -0.71228902341542847553 + Q * p;
541 Q = 0.662328840472002992063 + Q * p;
542 Q = 1.56221558398423026363 + Q * p;
543 Q = -1.56574558234175846809 + Q * p;
544 Q = -0.970005043303290640362 + Q * p;
545 Q = 1.0 + Q * p;
546 final double g = p * (p + 10);
547 final double r = P / Q;
548 result = g * Y + g * r;
549 } else if (q >= 0.25) {
550 //
551 // Rational approximation for 0.5 > q >= 0.25
552 //
553 // x = sqrt(-2*log(q)) / (Y + R(q))
554 //
555 // Where Y is a constant, and R(q) is optimised for a low
556 // absolute error compared to Y.
557 //
558 // double : Max error found: 7.403372e-17
559 // long double : Max error found: 6.084616e-20
560 // Maximum Deviation Found (error term) 4.811e-20
561 //
562 final float Y = 2.249481201171875f;
563 final double xs = q - 0.25f;
564 double P;
565 P = -3.67192254707729348546;
566 P = 21.1294655448340526258 + P * xs;
567 P = 17.445385985570866523 + P * xs;
568 P = -44.6382324441786960818 + P * xs;
569 P = -18.8510648058714251895 + P * xs;
570 P = 17.6447298408374015486 + P * xs;
571 P = 8.37050328343119927838 + P * xs;
572 P = 0.105264680699391713268 + P * xs;
573 P = -0.202433508355938759655 + P * xs;
574 double Q;
575 Q = 1.72114765761200282724;
576 Q = -22.6436933413139721736 + Q * xs;
577 Q = 10.8268667355460159008 + Q * xs;
578 Q = 48.5609213108739935468 + Q * xs;
579 Q = -20.1432634680485188801 + Q * xs;
580 Q = -28.6608180499800029974 + Q * xs;
581 Q = 3.9713437953343869095 + Q * xs;
582 Q = 6.24264124854247537712 + Q * xs;
583 Q = 1.0 + Q * xs;
584 final double g = Math.sqrt(-2 * Math.log(q));
585 final double r = P / Q;
586 result = g / (Y + r);
587 } else {
588 //
589 // For q < 0.25 we have a series of rational approximations all
590 // of the general form:
591 //
592 // let: x = sqrt(-log(q))
593 //
594 // Then the result is given by:
595 //
596 // x(Y+R(x-B))
597 //
598 // where Y is a constant, B is the lowest value of x for which
599 // the approximation is valid, and R(x-B) is optimised for a low
600 // absolute error compared to Y.
601 //
602 // Note that almost all code will really go through the first
603 // or maybe second approximation. After than we're dealing with very
604 // small input values indeed.
605 //
606 // Limit for a double: Math.sqrt(-Math.log(Double.MIN_VALUE)) = 27.28...
607 // Branches for x >= 44 (supporting 80 and 128 bit long double) have been removed.
608 final double x = Math.sqrt(-Math.log(q));
609 if (x < 3) {
610 // Max error found: 1.089051e-20
611 final float Y = 0.807220458984375f;
612 final double xs = x - 1.125f;
613 double P;
614 P = -0.681149956853776992068e-9;
615 P = 0.285225331782217055858e-7 + P * xs;
616 P = -0.679465575181126350155e-6 + P * xs;
617 P = 0.00214558995388805277169 + P * xs;
618 P = 0.0290157910005329060432 + P * xs;
619 P = 0.142869534408157156766 + P * xs;
620 P = 0.337785538912035898924 + P * xs;
621 P = 0.387079738972604337464 + P * xs;
622 P = 0.117030156341995252019 + P * xs;
623 P = -0.163794047193317060787 + P * xs;
624 P = -0.131102781679951906451 + P * xs;
625 double Q;
626 Q = 0.01105924229346489121;
627 Q = 0.152264338295331783612 + Q * xs;
628 Q = 0.848854343457902036425 + Q * xs;
629 Q = 2.59301921623620271374 + Q * xs;
630 Q = 4.77846592945843778382 + Q * xs;
631 Q = 5.38168345707006855425 + Q * xs;
632 Q = 3.46625407242567245975 + Q * xs;
633 Q = 1.0 + Q * xs;
634 final double R = P / Q;
635 result = Y * x + R * x;
636 } else if (x < 6) {
637 // Max error found: 8.389174e-21
638 final float Y = 0.93995571136474609375f;
639 final double xs = x - 3;
640 double P;
641 P = 0.266339227425782031962e-11;
642 P = -0.230404776911882601748e-9 + P * xs;
643 P = 0.460469890584317994083e-5 + P * xs;
644 P = 0.000157544617424960554631 + P * xs;
645 P = 0.00187123492819559223345 + P * xs;
646 P = 0.00950804701325919603619 + P * xs;
647 P = 0.0185573306514231072324 + P * xs;
648 P = -0.00222426529213447927281 + P * xs;
649 P = -0.0350353787183177984712 + P * xs;
650 double Q;
651 Q = 0.764675292302794483503e-4;
652 Q = 0.00263861676657015992959 + Q * xs;
653 Q = 0.0341589143670947727934 + Q * xs;
654 Q = 0.220091105764131249824 + Q * xs;
655 Q = 0.762059164553623404043 + Q * xs;
656 Q = 1.3653349817554063097 + Q * xs;
657 Q = 1.0 + Q * xs;
658 final double R = P / Q;
659 result = Y * x + R * x;
660 } else if (x < 18) {
661 // Max error found: 1.481312e-19
662 final float Y = 0.98362827301025390625f;
663 final double xs = x - 6;
664 double P;
665 P = 0.99055709973310326855e-16;
666 P = -0.281128735628831791805e-13 + P * xs;
667 P = 0.462596163522878599135e-8 + P * xs;
668 P = 0.449696789927706453732e-6 + P * xs;
669 P = 0.149624783758342370182e-4 + P * xs;
670 P = 0.000209386317487588078668 + P * xs;
671 P = 0.00105628862152492910091 + P * xs;
672 P = -0.00112951438745580278863 + P * xs;
673 P = -0.0167431005076633737133 + P * xs;
674 double Q;
675 Q = 0.282243172016108031869e-6;
676 Q = 0.275335474764726041141e-4 + Q * xs;
677 Q = 0.000964011807005165528527 + Q * xs;
678 Q = 0.0160746087093676504695 + Q * xs;
679 Q = 0.138151865749083321638 + Q * xs;
680 Q = 0.591429344886417493481 + Q * xs;
681 Q = 1.0 + Q * xs;
682 final double R = P / Q;
683 result = Y * x + R * x;
684 } else {
685 // x < 44
686 // Max error found: 5.697761e-20
687 final float Y = 0.99714565277099609375f;
688 final double xs = x - 18;
689 double P;
690 P = -0.116765012397184275695e-17;
691 P = 0.145596286718675035587e-11 + P * xs;
692 P = 0.411632831190944208473e-9 + P * xs;
693 P = 0.396341011304801168516e-7 + P * xs;
694 P = 0.162397777342510920873e-5 + P * xs;
695 P = 0.254723037413027451751e-4 + P * xs;
696 P = -0.779190719229053954292e-5 + P * xs;
697 P = -0.0024978212791898131227 + P * xs;
698 double Q;
699 Q = 0.509761276599778486139e-9;
700 Q = 0.144437756628144157666e-6 + Q * xs;
701 Q = 0.145007359818232637924e-4 + Q * xs;
702 Q = 0.000690538265622684595676 + Q * xs;
703 Q = 0.0169410838120975906478 + Q * xs;
704 Q = 0.207123112214422517181 + Q * xs;
705 Q = 1.0 + Q * xs;
706 final double R = P / Q;
707 result = Y * x + R * x;
708 }
709 }
710 return result;
711 }
712
713 /**
714 * Compute {@code exp(x*x)} with high accuracy. This is performed using
715 * information in the round-off from {@code x*x}.
716 *
717 * <p>This is accurate at large x to 1 ulp.
718 *
719 * <p>At small x the accuracy cannot be improved over using exp(x*x).
720 * This occurs at {@code x <= 1}.
721 *
722 * <p>Warning: This has no checks for overflow. The method is never called
723 * when {@code x*x > log(MAX_VALUE/2)}.
724 *
725 * @param x Value
726 * @return exp(x*x)
727 */
728 static double expxx(double x) {
729 // Note: If exp(a) overflows this can create NaN if the
730 // round-off b is negative or zero:
731 // exp(a) * exp1m(b) + exp(a)
732 // inf * 0 + inf or inf * -b + inf
733 final DD x2 = DD.ofSquare(x);
734 return expxx(x2.hi(), x2.lo());
735 }
736
737 /**
738 * Compute {@code exp(-x*x)} with high accuracy. This is performed using
739 * information in the round-off from {@code x*x}.
740 *
741 * <p>This is accurate at large x to 1 ulp until exp(-x*x) is close to
742 * sub-normal. For very small exp(-x*x) the adjustment is sub-normal and
743 * bits can be lost in the adjustment for a max observed error of {@code < 2} ulp.
744 *
745 * <p>At small x the accuracy cannot be improved over using exp(-x*x).
746 * This occurs at {@code x <= 1}.
747 *
748 * @param x Value
749 * @return exp(-x*x)
750 */
751 static double expmxx(double x) {
752 final DD x2 = DD.ofSquare(x);
753 return expxx(-x2.hi(), -x2.lo());
754 }
755
756 /**
757 * Compute {@code exp(a+b)} with high accuracy assuming {@code a+b = a}.
758 *
759 * <p>This is accurate at large positive a to 1 ulp. If a is negative and exp(a) is
760 * close to sub-normal a bit of precision may be lost when adjusting result
761 * as the adjustment is sub-normal (max observed error {@code < 2} ulp).
762 * For the use case of multiplication of a number less than 1 by exp(-x*x), a = -x*x,
763 * the result will be sub-normal and the rounding error is lost.
764 *
765 * <p>At small |a| the accuracy cannot be improved over using exp(a) as the
766 * round-off is too small to create terms that can adjust the standard result by
767 * more than 0.5 ulp. This occurs at {@code |a| <= 1}.
768 *
769 * @param a High bits of a split number
770 * @param b Low bits of a split number
771 * @return exp(a+b)
772 */
773 private static double expxx(double a, double b) {
774 // exp(a+b) = exp(a) * exp(b)
775 // = exp(a) * (exp(b) - 1) + exp(a)
776 // Assuming:
777 // 1. -746 < a < 710 for no under/overflow of exp(a)
778 // 2. a+b = a
779 // As b -> 0 then exp(b) -> 1; expm1(b) -> b
780 // The round-off b is limited to ~ 0.5 * ulp(746) ~ 5.68e-14
781 // and we can use an approximation for expm1 (x/1! + x^2/2! + ...)
782 // The second term is required for the expm1 result but the
783 // bits are not significant to change the following sum with exp(a)
784
785 final double ea = Math.exp(a);
786 // b ~ expm1(b)
787 return ea * b + ea;
788 }
789 }