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 package org.apache.commons.numbers.gamma;
18
19 /**
20 * Computes the difference between {@link Erf error function values}.
21 */
22 public final class ErfDifference {
23 /**
24 * This number solves {@code erf(x) = 0.5} within 1 ulp.
25 * More precisely, the current implementations of
26 * {@link Erf#value(double)} and {@link Erfc#value(double)} satisfy:
27 * <ul>
28 * <li>{@code Erf.value(X_CRIT) == 0.5},</li>
29 * <li>{@code Erf.value(Math.nextUp(X_CRIT)) > 0.5},</li>
30 * <li>{@code Erfc.value(X_CRIT) == 0.5}, and</li>
31 * <li>{@code Erfc.value(Math.nextUp(X_CRIT)) < 0.5}</li>
32 * </ul>
33 */
34 private static final double X_CRIT = 0.47693627620446993;
35
36 /** Private constructor. */
37 private ErfDifference() {
38 // intentionally empty.
39 }
40
41 /**
42 * The implementation uses either {@link Erf} or {@link Erfc},
43 * depending on which provides the most precise result.
44 *
45 * @param x1 First value.
46 * @param x2 Second value.
47 * @return {@link Erf#value(double) Erf.value(x2) - Erf.value(x1)}.
48 * @throws ArithmeticException if the algorithm fails to converge.
49 */
50 public static double value(double x1,
51 double x2) {
52 if (x1 > x2) {
53 return -value(x2, x1);
54 }
55 if (x1 < -X_CRIT) {
56 if (x2 < 0) {
57 return Erfc.value(-x2) - Erfc.value(-x1);
58 }
59 return Erf.value(x2) - Erf.value(x1);
60 }
61 if (x2 > X_CRIT &&
62 x1 > 0) {
63 return Erfc.value(x1) - Erfc.value(x2);
64 }
65 return Erf.value(x2) - Erf.value(x1);
66 }
67 }