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 * <a href="http://mathworld.wolfram.com/LanczosApproximation.html">
021 * Lanczos approximation</a> to the Gamma function.
022 *
023 * It is related to the Gamma function by the following equation
024 * \[
025 * \Gamma(x) = \sqrt{2\pi} \, \frac{(g + x + \frac{1}{2})^{x + \frac{1}{2}} \, e^{-(g + x + \frac{1}{2})} \, \mathrm{lanczos}(x)}
026 *                                 {x}
027 * \]
028 * where \( g \) is the Lanczos constant.
029 *
030 * See equations (1) through (5), and Paul Godfrey's
031 * <a href="http://my.fit.edu/~gabdo/gamma.txt">Note on the computation
032 * of the convergent Lanczos complex Gamma approximation</a>.
033 */
034public final class LanczosApproximation {
035    /** \( g = \frac{607}{128} \). */
036    private static final double LANCZOS_G = 607d / 128d;
037    /** Lanczos coefficients. */
038    private static final double[] LANCZOS = {
039        0.99999999999999709182,
040        57.156235665862923517,
041        -59.597960355475491248,
042        14.136097974741747174,
043        -0.49191381609762019978,
044        .33994649984811888699e-4,
045        .46523628927048575665e-4,
046        -.98374475304879564677e-4,
047        .15808870322491248884e-3,
048        -.21026444172410488319e-3,
049        .21743961811521264320e-3,
050        -.16431810653676389022e-3,
051        .84418223983852743293e-4,
052        -.26190838401581408670e-4,
053        .36899182659531622704e-5,
054    };
055
056    /** Private constructor. */
057    private LanczosApproximation() {
058        // intentionally empty.
059    }
060
061    /**
062     * Computes the Lanczos approximation.
063     *
064     * @param x Argument.
065     * @return the Lanczos approximation.
066     */
067    public static double value(final double x) {
068        double sum = 0;
069        for (int i = LANCZOS.length - 1; i > 0; i--) {
070            sum += LANCZOS[i] / (x + i);
071        }
072        return sum + LANCZOS[0];
073    }
074
075    /**
076     * @return the Lanczos constant \( g = \frac{607}{128} \).
077     */
078    public static double g() {
079        return LANCZOS_G;
080    }
081}