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        // intentional 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.isNaN(x) || Double.isInfinite(x)) {
051            return x;
052        }
053
054        if (x > 0 && x <= S_LIMIT) {
055            return 1 / (x * x);
056        }
057
058        if (x >= C_LIMIT) {
059            final double inv = 1 / (x * x);
060            //  1    1      1       1       1
061            //  - + ---- + ---- - ----- + -----
062            //  x      2      3       5       7
063            //      2 x    6 x    30 x    42 x
064            return 1 / x + inv / 2 + inv / x * (F_1_6 - inv * (F_1_30 + F_1_42 * inv));
065        }
066
067        return value(x + 1) + 1 / (x * x);
068    }
069}