org.apache.commons.numbers.gamma.Trigamma Maven / Gradle / Ivy
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.numbers.gamma;
/**
* Trigamma function.
*
* It is the derivative of the {@link Digamma digamma function}:
* \( \psi_1(x) = \frac{d^2}{dx^2} (\ln \Gamma(x)) \).
*/
public final class Trigamma {
/** C limit. */
private static final double C_LIMIT = 49;
/** S limit. */
private static final double S_LIMIT = 1e-5;
/** Fraction. */
private static final double F_1_6 = 1d / 6;
/** Fraction. */
private static final double F_1_30 = 1d / 30;
/** Fraction. */
private static final double F_1_42 = 1d / 42;
/** Private constructor. */
private Trigamma() {
// intentional empty.
}
/**
* Computes the trigamma function.
*
* @param x Argument.
* @return trigamma(x) to within {@code 1e-8} relative or absolute error whichever is larger.
*/
public static double value(double x) {
if (Double.isNaN(x) || Double.isInfinite(x)) {
return x;
}
if (x > 0 && x <= S_LIMIT) {
return 1 / (x * x);
}
double trigamma = 0;
while (x < C_LIMIT) {
trigamma += 1 / (x * x);
x += 1;
}
final double inv = 1 / (x * x);
// 1 1 1 1 1
// - + ---- + ---- - ----- + -----
// x 2 3 5 7
// 2 x 6 x 30 x 42 x
trigamma += 1 / x + inv / 2 + inv / x * (F_1_6 - inv * (F_1_30 + F_1_42 * inv));
return trigamma;
}
}