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://en.wikipedia.org/wiki/Trigamma_function">Trigamma function</a>.
021 *
022 * It is the derivative of the {@link Digamma digamma function}:
023 * \( \psi_1(x) = \frac{d^2}{dx^2} (\ln \Gamma(x)) \).
024 */
025public final class Trigamma {
026    /** C limit. */
027    private static final double C_LIMIT = 49;
028
029    /** S limit. */
030    private static final double S_LIMIT = 1e-5;
031    /** Fraction. */
032    private static final double F_1_6 = 1d / 6;
033    /** Fraction. */
034    private static final double F_1_30 = 1d / 30;
035    /** Fraction. */
036    private static final double F_1_42 = 1d / 42;
037
038    /** Private constructor. */
039    private Trigamma() {
040        // intentionally empty.
041    }
042
043    /**
044     * Computes the trigamma function.
045     *
046     * @param x Argument.
047     * @return trigamma(x) to within {@code 1e-8} relative or absolute error whichever is larger.
048     */
049    public static double value(double x) {
050        if (!Double.isFinite(x)) {
051            return x;
052        }
053
054        if (x > 0 && x <= S_LIMIT) {
055            return 1 / (x * x);
056        }
057
058        double trigamma = 0;
059        while (x < C_LIMIT) {
060            trigamma += 1 / (x * x);
061            x += 1;
062        }
063
064        final double inv = 1 / (x * x);
065        //  1    1      1       1       1
066        //  - + ---- + ---- - ----- + -----
067        //  x      2      3       5       7
068        //      2 x    6 x    30 x    42 x
069        trigamma += 1 / x + inv / 2 + inv / x * (F_1_6 - inv * (F_1_30 + F_1_42 * inv));
070
071        return trigamma;
072    }
073}