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/**
021 * <a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
022 * function</a>.
023 * <p>
024 * The <a href="http://mathworld.wolfram.com/GammaFunction.html">gamma
025 * function</a> can be seen to extend the factorial function to cover real and
026 * complex numbers, but with its argument shifted by {@code -1}. This
027 * implementation supports real numbers.
028 * </p>
029 * <p>
030 * This class is immutable.
031 * </p>
032 */
033public final class Gamma {
034    /** The threshold value for choosing the Lanczos approximation. */
035    private static final double LANCZOS_THRESHOLD = 20;
036
037    /** &radic;(2&pi;). */
038    private static final double SQRT_TWO_PI = 2.506628274631000502;
039
040    /** Private constructor. */
041    private Gamma() {
042        // intentionally empty.
043    }
044
045    /**
046     * Computes the value of \( \Gamma(x) \).
047     * <p>
048     * Based on the <em>NSWC Library of Mathematics Subroutines</em> double
049     * precision implementation, {@code DGAMMA}.
050     *
051     * @param x Argument.
052     * @return \( \Gamma(x) \)
053     */
054    public static double value(final double x) {
055
056        if ((x == Math.rint(x)) && (x <= 0.0)) {
057            return Double.NaN;
058        }
059
060        final double absX = Math.abs(x);
061        if (absX <= LANCZOS_THRESHOLD) {
062            if (x >= 1) {
063                /*
064                 * From the recurrence relation
065                 * Gamma(x) = (x - 1) * ... * (x - n) * Gamma(x - n),
066                 * then
067                 * Gamma(t) = 1 / [1 + InvGamma1pm1.value(t - 1)],
068                 * where t = x - n. This means that t must satisfy
069                 * -0.5 <= t - 1 <= 1.5.
070                 */
071                double prod = 1;
072                double t = x;
073                while (t > 2.5) {
074                    t -= 1;
075                    prod *= t;
076                }
077                return prod / (1 + InvGamma1pm1.value(t - 1));
078            } else {
079                /*
080                 * From the recurrence relation
081                 * Gamma(x) = Gamma(x + n + 1) / [x * (x + 1) * ... * (x + n)]
082                 * then
083                 * Gamma(x + n + 1) = 1 / [1 + InvGamma1pm1.value(x + n)],
084                 * which requires -0.5 <= x + n <= 1.5.
085                 */
086                double prod = x;
087                double t = x;
088                while (t < -0.5) {
089                    t += 1;
090                    prod *= t;
091                }
092                return 1 / (prod * (1 + InvGamma1pm1.value(t)));
093            }
094        } else {
095            final double y = absX + LanczosApproximation.g() + 0.5;
096            final double gammaAbs = SQRT_TWO_PI / absX *
097                                    Math.pow(y, absX + 0.5) *
098                                    Math.exp(-y) * LanczosApproximation.value(absX);
099            if (x > 0) {
100                return gammaAbs;
101            } else {
102                /*
103                 * From the reflection formula
104                 * Gamma(x) * Gamma(1 - x) * sin(pi * x) = pi,
105                 * and the recurrence relation
106                 * Gamma(1 - x) = -x * Gamma(-x),
107                 * it is found
108                 * Gamma(x) = -pi / [x * sin(pi * x) * Gamma(-x)].
109                 */
110                return -Math.PI / (x * Math.sin(Math.PI * x) * gammaAbs);
111            }
112        }
113    }
114}