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 * Function \( \ln \Gamma(x) \).
021 *
022 * Class is immutable.
023 */
024public final class LogGamma {
025    /** Lanczos constant. */
026    private static final double LANCZOS_G = 607d / 128d;
027    /** Performance. */
028    private static final double HALF_LOG_2_PI = 0.5 * Math.log(2.0 * Math.PI);
029
030    /** Private constructor. */
031    private LogGamma() {
032        // intentionally empty.
033    }
034
035    /**
036     * Computes the function \( \ln \Gamma(x) \) for {@code x >= 0}.
037     *
038     * For {@code x <= 8}, the implementation is based on the double precision
039     * implementation in the <em>NSWC Library of Mathematics Subroutines</em>,
040     * {@code DGAMLN}. For {@code x >= 8}, the implementation is based on
041     * <ul>
042     * <li><a href="http://mathworld.wolfram.com/GammaFunction.html">Gamma
043     *     Function</a>, equation (28).</li>
044     * <li><a href="http://mathworld.wolfram.com/LanczosApproximation.html">
045     *     Lanczos Approximation</a>, equations (1) through (5).</li>
046     * <li><a href="http://my.fit.edu/~gabdo/gamma.txt">Paul Godfrey, A note on
047     *     the computation of the convergent Lanczos complex Gamma
048     *     approximation</a></li>
049     * </ul>
050     *
051     * @param x Argument.
052     * @return \( \ln \Gamma(x) \), or {@code NaN} if {@code x <= 0}.
053     */
054    public static double value(double x) {
055        if (Double.isNaN(x) || (x <= 0.0)) {
056            return Double.NaN;
057        } else if (x < 0.5) {
058            return LogGamma1p.value(x) - Math.log(x);
059        } else if (x <= 2.5) {
060            return LogGamma1p.value((x - 0.5) - 0.5);
061        } else if (x <= 8.0) {
062            final int n = (int) Math.floor(x - 1.5);
063            double prod = 1.0;
064            for (int i = 1; i <= n; i++) {
065                prod *= x - i;
066            }
067            return LogGamma1p.value(x - (n + 1)) + Math.log(prod);
068        } else {
069            final double sum = LanczosApproximation.value(x);
070            final double tmp = x + LANCZOS_G + .5;
071            return ((x + .5) * Math.log(tmp)) - tmp +
072                HALF_LOG_2_PI + Math.log(sum / x);
073        }
074    }
075}