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
019import java.text.MessageFormat;
020
021import org.apache.commons.numbers.fraction.ContinuedFraction;
022
023/**
024 * <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
025 * Regularized Gamma functions</a>.
026 *
027 * Class is immutable.
028 */
029public final class RegularizedGamma {
030    /** Maximum allowed numerical error. */
031    private static final double DEFAULT_EPSILON = 1e-15;
032
033    /** Private constructor. */
034    private RegularizedGamma() {
035        // intentionally empty.
036    }
037
038    /**
039     * \( P(a, x) \) <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
040     * regularized Gamma function</a>.
041     *
042     * Class is immutable.
043     */
044    public static final class P {
045        /** Prevent instantiation. */
046        private P() {}
047
048        /**
049         * Computes the regularized gamma function \( P(a, x) \).
050         *
051         * @param a Argument.
052         * @param x Argument.
053         * @return \( P(a, x) \).
054         * @throws ArithmeticException if the continued fraction fails to converge.
055         */
056        public static double value(double a,
057                                   double x) {
058            return value(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
059        }
060
061        /**
062         * Computes the regularized gamma function \( P(a, x) \).
063         *
064         * The implementation of this method is based on:
065         * <ul>
066         *  <li>
067         *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
068         *   Regularized Gamma Function</a>, equation (1)
069         *  </li>
070         *  <li>
071         *   <a href="http://mathworld.wolfram.com/IncompleteGammaFunction.html">
072         *   Incomplete Gamma Function</a>, equation (4).
073         *  </li>
074         *  <li>
075         *   <a href="http://mathworld.wolfram.com/ConfluentHypergeometricFunctionoftheFirstKind.html">
076         *   Confluent Hypergeometric Function of the First Kind</a>, equation (1).
077         *  </li>
078         * </ul>
079         *
080         * @param a Argument.
081         * @param x Argument.
082         * @param epsilon Tolerance in continued fraction evaluation.
083         * @param maxIterations Maximum number of iterations in continued fraction evaluation.
084         * @return \( P(a, x) \).
085         * @throws ArithmeticException if the continued fraction fails to converge.
086         */
087        public static double value(double a,
088                                   double x,
089                                   double epsilon,
090                                   int maxIterations) {
091            if (Double.isNaN(a) ||
092                Double.isNaN(x) ||
093                a <= 0 ||
094                x < 0) {
095                return Double.NaN;
096            } else if (x == 0) {
097                return 0;
098            } else if (x >= a + 1) {
099                // Q should converge faster in this case.
100                return 1 - RegularizedGamma.Q.value(a, x, epsilon, maxIterations);
101            } else {
102                // Series.
103                double n = 0; // current element index
104                double an = 1 / a; // n-th element in the series
105                double sum = an; // partial sum
106                while (Math.abs(an / sum) > epsilon &&
107                       n < maxIterations &&
108                       sum < Double.POSITIVE_INFINITY) {
109                    // compute next element in the series
110                    n += 1;
111                    an *= x / (a + n);
112
113                    // update partial sum
114                    sum += an;
115                }
116                if (n >= maxIterations) {
117                    throw new ArithmeticException(
118                            MessageFormat.format("Failed to converge within {0} iterations", maxIterations));
119                } else if (Double.isInfinite(sum)) {
120                    return 1;
121                } else {
122                    // Ensure result is in the range [0, 1]
123                    final double result = Math.exp(-x + (a * Math.log(x)) - LogGamma.value(a)) * sum;
124                    return result > 1.0 ? 1.0 : result;
125                }
126            }
127        }
128    }
129
130    /**
131     * Creates the \( Q(a, x) \equiv 1 - P(a, x) \) <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
132     * regularized Gamma function</a>.
133     *
134     * Class is immutable.
135     */
136    public static final class Q {
137        /** Prevent instantiation. */
138        private Q() {}
139
140        /**
141         * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \).
142         *
143         * @param a Argument.
144         * @param x Argument.
145         * @return \( Q(a, x) \).
146         * @throws ArithmeticException if the continued fraction fails to converge.
147         */
148        public static double value(double a,
149                                   double x) {
150            return value(a, x, DEFAULT_EPSILON, Integer.MAX_VALUE);
151        }
152
153        /**
154         * Computes the regularized gamma function \( Q(a, x) = 1 - P(a, x) \).
155         *
156         * The implementation of this method is based on:
157         * <ul>
158         *  <li>
159         *   <a href="http://mathworld.wolfram.com/RegularizedGammaFunction.html">
160         *   Regularized Gamma Function</a>, equation (1).
161         *  </li>
162         *  <li>
163         *   <a href="http://functions.wolfram.com/GammaBetaErf/GammaRegularized/10/0003/">
164         *   Regularized incomplete gamma function: Continued fraction representations
165         *   (formula 06.08.10.0003)</a>
166         *  </li>
167         * </ul>
168         *
169         * @param a Argument.
170         * @param x Argument.
171         * @param epsilon Tolerance in continued fraction evaluation.
172         * @param maxIterations Maximum number of iterations in continued fraction evaluation.
173         * @throws ArithmeticException if the continued fraction fails to converge.
174         * @return \( Q(a, x) \).
175         */
176        public static double value(final double a,
177                                   double x,
178                                   double epsilon,
179                                   int maxIterations) {
180            if (Double.isNaN(a) ||
181                Double.isNaN(x) ||
182                a <= 0 ||
183                x < 0) {
184                return Double.NaN;
185            } else if (x == 0) {
186                return 1;
187            } else if (x < a + 1) {
188                // P should converge faster in this case.
189                return 1 - RegularizedGamma.P.value(a, x, epsilon, maxIterations);
190            } else {
191                final ContinuedFraction cf = new ContinuedFraction() {
192                        /** {@inheritDoc} */
193                        @Override
194                        protected double getA(int n, double x) {
195                            return n * (a - n);
196                        }
197
198                        /** {@inheritDoc} */
199                        @Override
200                        protected double getB(int n, double x) {
201                            return ((2 * n) + 1) - a + x;
202                        }
203                    };
204
205                return Math.exp(-x + (a * Math.log(x)) - LogGamma.value(a)) /
206                    cf.evaluate(x, epsilon, maxIterations);
207            }
208        }
209    }
210}