View Javadoc
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   * Function \( \frac{1}{\Gamma(1 + x)} - 1 \).
21   *
22   * Class is immutable.
23   */
24  final class InvGamma1pm1 {
25      /*
26       * Constants copied from DGAM1 in the NSWC library.
27       */
28      /** The constant {@code A0} defined in {@code DGAM1}. */
29      private static final double INV_GAMMA1P_M1_A0 = .611609510448141581788E-08;
30      /** The constant {@code A1} defined in {@code DGAM1}. */
31      private static final double INV_GAMMA1P_M1_A1 = .624730830116465516210E-08;
32      /** The constant {@code B1} defined in {@code DGAM1}. */
33      private static final double INV_GAMMA1P_M1_B1 = .203610414066806987300E+00;
34      /** The constant {@code B2} defined in {@code DGAM1}. */
35      private static final double INV_GAMMA1P_M1_B2 = .266205348428949217746E-01;
36      /** The constant {@code B3} defined in {@code DGAM1}. */
37      private static final double INV_GAMMA1P_M1_B3 = .493944979382446875238E-03;
38      /** The constant {@code B4} defined in {@code DGAM1}. */
39      private static final double INV_GAMMA1P_M1_B4 = -.851419432440314906588E-05;
40      /** The constant {@code B5} defined in {@code DGAM1}. */
41      private static final double INV_GAMMA1P_M1_B5 = -.643045481779353022248E-05;
42      /** The constant {@code B6} defined in {@code DGAM1}. */
43      private static final double INV_GAMMA1P_M1_B6 = .992641840672773722196E-06;
44      /** The constant {@code B7} defined in {@code DGAM1}. */
45      private static final double INV_GAMMA1P_M1_B7 = -.607761895722825260739E-07;
46      /** The constant {@code B8} defined in {@code DGAM1}. */
47      private static final double INV_GAMMA1P_M1_B8 = .195755836614639731882E-09;
48      /** The constant {@code P0} defined in {@code DGAM1}. */
49      private static final double INV_GAMMA1P_M1_P0 = .6116095104481415817861E-08;
50      /** The constant {@code P1} defined in {@code DGAM1}. */
51      private static final double INV_GAMMA1P_M1_P1 = .6871674113067198736152E-08;
52      /** The constant {@code P2} defined in {@code DGAM1}. */
53      private static final double INV_GAMMA1P_M1_P2 = .6820161668496170657918E-09;
54      /** The constant {@code P3} defined in {@code DGAM1}. */
55      private static final double INV_GAMMA1P_M1_P3 = .4686843322948848031080E-10;
56      /** The constant {@code P4} defined in {@code DGAM1}. */
57      private static final double INV_GAMMA1P_M1_P4 = .1572833027710446286995E-11;
58      /** The constant {@code P5} defined in {@code DGAM1}. */
59      private static final double INV_GAMMA1P_M1_P5 = -.1249441572276366213222E-12;
60      /** The constant {@code P6} defined in {@code DGAM1}. */
61      private static final double INV_GAMMA1P_M1_P6 = .4343529937408594255178E-14;
62      /** The constant {@code Q1} defined in {@code DGAM1}. */
63      private static final double INV_GAMMA1P_M1_Q1 = .3056961078365221025009E+00;
64      /** The constant {@code Q2} defined in {@code DGAM1}. */
65      private static final double INV_GAMMA1P_M1_Q2 = .5464213086042296536016E-01;
66      /** The constant {@code Q3} defined in {@code DGAM1}. */
67      private static final double INV_GAMMA1P_M1_Q3 = .4956830093825887312020E-02;
68      /** The constant {@code Q4} defined in {@code DGAM1}. */
69      private static final double INV_GAMMA1P_M1_Q4 = .2692369466186361192876E-03;
70      /** The constant {@code C} defined in {@code DGAM1}. */
71      private static final double INV_GAMMA1P_M1_C = -.422784335098467139393487909917598E+00;
72      /** The constant {@code C0} defined in {@code DGAM1}. */
73      private static final double INV_GAMMA1P_M1_C0 = .577215664901532860606512090082402E+00;
74      /** The constant {@code C1} defined in {@code DGAM1}. */
75      private static final double INV_GAMMA1P_M1_C1 = -.655878071520253881077019515145390E+00;
76      /** The constant {@code C2} defined in {@code DGAM1}. */
77      private static final double INV_GAMMA1P_M1_C2 = -.420026350340952355290039348754298E-01;
78      /** The constant {@code C3} defined in {@code DGAM1}. */
79      private static final double INV_GAMMA1P_M1_C3 = .166538611382291489501700795102105E+00;
80      /** The constant {@code C4} defined in {@code DGAM1}. */
81      private static final double INV_GAMMA1P_M1_C4 = -.421977345555443367482083012891874E-01;
82      /** The constant {@code C5} defined in {@code DGAM1}. */
83      private static final double INV_GAMMA1P_M1_C5 = -.962197152787697356211492167234820E-02;
84      /** The constant {@code C6} defined in {@code DGAM1}. */
85      private static final double INV_GAMMA1P_M1_C6 = .721894324666309954239501034044657E-02;
86      /** The constant {@code C7} defined in {@code DGAM1}. */
87      private static final double INV_GAMMA1P_M1_C7 = -.116516759185906511211397108401839E-02;
88      /** The constant {@code C8} defined in {@code DGAM1}. */
89      private static final double INV_GAMMA1P_M1_C8 = -.215241674114950972815729963053648E-03;
90      /** The constant {@code C9} defined in {@code DGAM1}. */
91      private static final double INV_GAMMA1P_M1_C9 = .128050282388116186153198626328164E-03;
92      /** The constant {@code C10} defined in {@code DGAM1}. */
93      private static final double INV_GAMMA1P_M1_C10 = -.201348547807882386556893914210218E-04;
94      /** The constant {@code C11} defined in {@code DGAM1}. */
95      private static final double INV_GAMMA1P_M1_C11 = -.125049348214267065734535947383309E-05;
96      /** The constant {@code C12} defined in {@code DGAM1}. */
97      private static final double INV_GAMMA1P_M1_C12 = .113302723198169588237412962033074E-05;
98      /** The constant {@code C13} defined in {@code DGAM1}. */
99      private static final double INV_GAMMA1P_M1_C13 = -.205633841697760710345015413002057E-06;
100 
101     /** Private constructor. */
102     private InvGamma1pm1() {
103         // intentionally empty.
104     }
105 
106     /**
107      * Computes the function \( \frac{1}{\Gamma(1 + x)} - 1 \) for {@code -0.5 <= x <= 1.5}.
108      *
109      * This implementation is based on the double precision implementation in
110      * the <em>NSWC Library of Mathematics Subroutines</em>, {@code DGAM1}.
111      *
112      * @param x Argument.
113      * @return \( \frac{1}{\Gamma(1 + x)} - 1 \)
114      * @throws IllegalArgumentException if {@code x < -0.5} or {@code x > 1.5}
115      */
116     public static double value(final double x) {
117         if (x < -0.5 || x > 1.5) {
118             throw new GammaException(GammaException.OUT_OF_RANGE, x, -0.5, 1.5);
119         }
120 
121         final double t = x <= 0.5 ? x : (x - 0.5) - 0.5;
122         if (t < 0) {
123             final double a = INV_GAMMA1P_M1_A0 + t * INV_GAMMA1P_M1_A1;
124             double b = INV_GAMMA1P_M1_B8;
125             b = INV_GAMMA1P_M1_B7 + t * b;
126             b = INV_GAMMA1P_M1_B6 + t * b;
127             b = INV_GAMMA1P_M1_B5 + t * b;
128             b = INV_GAMMA1P_M1_B4 + t * b;
129             b = INV_GAMMA1P_M1_B3 + t * b;
130             b = INV_GAMMA1P_M1_B2 + t * b;
131             b = INV_GAMMA1P_M1_B1 + t * b;
132             b = 1.0 + t * b;
133 
134             double c = INV_GAMMA1P_M1_C13 + t * (a / b);
135             c = INV_GAMMA1P_M1_C12 + t * c;
136             c = INV_GAMMA1P_M1_C11 + t * c;
137             c = INV_GAMMA1P_M1_C10 + t * c;
138             c = INV_GAMMA1P_M1_C9 + t * c;
139             c = INV_GAMMA1P_M1_C8 + t * c;
140             c = INV_GAMMA1P_M1_C7 + t * c;
141             c = INV_GAMMA1P_M1_C6 + t * c;
142             c = INV_GAMMA1P_M1_C5 + t * c;
143             c = INV_GAMMA1P_M1_C4 + t * c;
144             c = INV_GAMMA1P_M1_C3 + t * c;
145             c = INV_GAMMA1P_M1_C2 + t * c;
146             c = INV_GAMMA1P_M1_C1 + t * c;
147             c = INV_GAMMA1P_M1_C + t * c;
148             if (x > 0.5) {
149                 return t * c / x;
150             }
151             return x * ((c + 0.5) + 0.5);
152         }
153         double p = INV_GAMMA1P_M1_P6;
154         p = INV_GAMMA1P_M1_P5 + t * p;
155         p = INV_GAMMA1P_M1_P4 + t * p;
156         p = INV_GAMMA1P_M1_P3 + t * p;
157         p = INV_GAMMA1P_M1_P2 + t * p;
158         p = INV_GAMMA1P_M1_P1 + t * p;
159         p = INV_GAMMA1P_M1_P0 + t * p;
160 
161         double q = INV_GAMMA1P_M1_Q4;
162         q = INV_GAMMA1P_M1_Q3 + t * q;
163         q = INV_GAMMA1P_M1_Q2 + t * q;
164         q = INV_GAMMA1P_M1_Q1 + t * q;
165         q = 1.0 + t * q;
166 
167         double c = INV_GAMMA1P_M1_C13 + (p / q) * t;
168         c = INV_GAMMA1P_M1_C12 + t * c;
169         c = INV_GAMMA1P_M1_C11 + t * c;
170         c = INV_GAMMA1P_M1_C10 + t * c;
171         c = INV_GAMMA1P_M1_C9 + t * c;
172         c = INV_GAMMA1P_M1_C8 + t * c;
173         c = INV_GAMMA1P_M1_C7 + t * c;
174         c = INV_GAMMA1P_M1_C6 + t * c;
175         c = INV_GAMMA1P_M1_C5 + t * c;
176         c = INV_GAMMA1P_M1_C4 + t * c;
177         c = INV_GAMMA1P_M1_C3 + t * c;
178         c = INV_GAMMA1P_M1_C2 + t * c;
179         c = INV_GAMMA1P_M1_C1 + t * c;
180         c = INV_GAMMA1P_M1_C0 + t * c;
181 
182         if (x > 0.5) {
183             return (t / x) * ((c - 0.5) - 0.5);
184         }
185         return x * c;
186     }
187 }