001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.gamma;
018
019/**
020 * Computes the difference between {@link Erf error function values}.
021 */
022public final class ErfDifference {
023    /**
024     * This number solves {@code erf(x) = 0.5} within 1 ulp.
025     * More precisely, the current implementations of
026     * {@link Erf#value(double)} and {@link Erfc#value(double)} satisfy:
027     * <ul>
028     *  <li>{@code Erf.value(X_CRIT) == 0.5},</li>
029     *  <li>{@code Erf.value(Math.nextUp(X_CRIT)) > 0.5},</li>
030     *  <li>{@code Erfc.value(X_CRIT) == 0.5}, and</li>
031     *  <li>{@code Erfc.value(Math.nextUp(X_CRIT)) < 0.5}</li>
032     * </ul>
033     */
034    private static final double X_CRIT = 0.47693627620446993;
035
036    /** Private constructor. */
037    private ErfDifference() {
038        // intentionally empty.
039    }
040
041    /**
042     * The implementation uses either {@link Erf} or {@link Erfc},
043     * depending on which provides the most precise result.
044     *
045     * @param x1 First value.
046     * @param x2 Second value.
047     * @return {@link Erf#value(double) Erf.value(x2) - Erf.value(x1)}.
048     * @throws ArithmeticException if the algorithm fails to converge.
049     */
050    public static double value(double x1,
051                               double x2) {
052        if (x1 > x2) {
053            return -value(x2, x1);
054        }
055        if (x1 < -X_CRIT) {
056            if (x2 < 0) {
057                return Erfc.value(-x2) - Erfc.value(-x1);
058            }
059            return Erf.value(x2) - Erf.value(x1);
060        }
061        if (x2 > X_CRIT &&
062            x1 > 0) {
063            return Erfc.value(x1) - Erfc.value(x2);
064        }
065        return Erf.value(x2) - Erf.value(x1);
066    }
067}