1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.numbers.gamma;
18
19 /**
20 * <a href="http://en.wikipedia.org/wiki/Trigamma_function">Trigamma function</a>.
21 *
22 * It is the derivative of the {@link Digamma digamma function}:
23 * \( \psi_1(x) = \frac{d^2}{dx^2} (\ln \Gamma(x)) \).
24 */
25 public final class Trigamma {
26 /** C limit. */
27 private static final double C_LIMIT = 49;
28
29 /** S limit. */
30 private static final double S_LIMIT = 1e-5;
31 /** Fraction. */
32 private static final double F_1_6 = 1d / 6;
33 /** Fraction. */
34 private static final double F_1_30 = 1d / 30;
35 /** Fraction. */
36 private static final double F_1_42 = 1d / 42;
37
38 /** Private constructor. */
39 private Trigamma() {
40 // intentionally empty.
41 }
42
43 /**
44 * Computes the trigamma function.
45 *
46 * @param x Argument.
47 * @return trigamma(x) to within {@code 1e-8} relative or absolute error whichever is larger.
48 */
49 public static double value(double x) {
50 if (!Double.isFinite(x)) {
51 return x;
52 }
53
54 if (x > 0 && x <= S_LIMIT) {
55 return 1 / (x * x);
56 }
57
58 double trigamma = 0;
59 while (x < C_LIMIT) {
60 trigamma += 1 / (x * x);
61 x += 1;
62 }
63
64 final double inv = 1 / (x * x);
65 // 1 1 1 1 1
66 // - + ---- + ---- - ----- + -----
67 // x 2 3 5 7
68 // 2 x 6 x 30 x 42 x
69 trigamma += 1 / x + inv / 2 + inv / x * (F_1_6 - inv * (F_1_30 + F_1_42 * inv));
70
71 return trigamma;
72 }
73 }