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 * Special math functions.
21 *
22 * @since 1.1
23 */
24 final class SpecialMath {
25 /** Minimum x for log1pmx(x). */
26 private static final double X_MIN = -1;
27 /** Low threshold to use log1p(x) - x. */
28 private static final double X_LOW = -0.79149064;
29 /** High threshold to use log1p(x) - x. */
30 private static final double X_HIGH = 1;
31 /** 2^-6. */
32 private static final double TWO_POW_M6 = 0x1.0p-6;
33 /** 2^-12. */
34 private static final double TWO_POW_M12 = 0x1.0p-12;
35 /** 2^-20. */
36 private static final double TWO_POW_M20 = 0x1.0p-20;
37 /** 2^-53. */
38 private static final double TWO_POW_M53 = 0x1.0p-53;
39
40 /** Private constructor. */
41 private SpecialMath() {
42 // intentionally empty.
43 }
44
45 /**
46 * Returns {@code log(1 + x) - x}. This function is accurate when {@code x -> 0}.
47 *
48 * <p>This function uses a Taylor series expansion when x is small ({@code |x| < 0.01}):
49 *
50 * <pre>
51 * ln(1 + x) - x = -x^2/2 + x^3/3 - x^4/4 + ...
52 * </pre>
53 *
54 * <p>or around 0 ({@code -0.791 <= x <= 1}):
55 *
56 * <pre>
57 * ln(x + a) = ln(a) + 2 [z + z^3/3 + z^5/5 + z^7/7 + ... ]
58 *
59 * z = x / (2a + x)
60 * </pre>
61 *
62 * <p>For a = 1:
63 *
64 * <pre>
65 * ln(x + 1) - x = -x + 2 [z + z^3/3 + z^5/5 + z^7/7 + ... ]
66 * = z * (-x + 2z^2 [ 1/3 + z^2/5 + z^4/7 + ... ])
67 * </pre>
68 *
69 * <p>The code is based on the {@code log1pmx} documentation for the <a
70 * href="https://rdrr.io/rforge/DPQ/man/log1pmx.html">R DPQ package</a> with addition of the
71 * direct Taylor series for tiny x.
72 *
73 * <p>See Abramowitz, M. and Stegun, I. A. (1972) Handbook of Mathematical Functions. New York:
74 * Dover. Formulas 4.1.24 and 4.2.29, p.68. <a
75 * href="https://en.wikipedia.org/wiki/Abramowitz_and_Stegun">Wikipedia: Abramowitz_and_Stegun</a>
76 * provides links to the full text which is in public domain.
77 *
78 * @param x Value x
79 * @return {@code log(1 + x) - x}
80 */
81 static double log1pmx(double x) {
82 // -1 is the minimum supported value
83 if (x <= X_MIN) {
84 return x == X_MIN ? Double.NEGATIVE_INFINITY : Double.NaN;
85 }
86 // Use the thresholds documented in the R implementation
87 if (x < X_LOW || x > X_HIGH) {
88 return Math.log1p(x) - x;
89 }
90 final double a = Math.abs(x);
91
92 // Addition to the R version for small x.
93 // Use a direct Taylor series:
94 // ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + ...
95 if (a < TWO_POW_M6) {
96 return log1pmxSmall(x, a);
97 }
98
99 // The use of the following series is fast converging:
100 // ln(x + 1) - x = -x + 2 [z + z^3/3 + z^5/5 + z^7/7 + ... ]
101 // = z * (-x + 2z^2 [ 1/3 + z^2/5 + z^4/7 + ... ])
102 // z = x / (2 + x)
103 //
104 // Tests show this is more accurate when |x| > 1e-4 than the direct Taylor series.
105 // The direct series can be modified to sum multiple terms together for a small
106 // increase in precision to a closer match to this variation but the direct series
107 // takes approximately 3x longer to converge.
108
109 final double z = x / (2 + x);
110 final double zz = z * z;
111
112 // Series sum
113 // sum(k=0,...,Inf; zz^k/(3+k*2)) = 1/3 + zz/5 + zz^2/7 + zz^3/9 + ... )
114
115 double sum = 1.0 / 3;
116 double numerator = 1;
117 int denominator = 3;
118 for (;;) {
119 numerator *= zz;
120 denominator += 2;
121 final double sum2 = sum + numerator / denominator;
122 // Since |x| <= 1 the additional terms will reduce in magnitude.
123 // Iterate until convergence. Expected iterations:
124 // x iterations
125 // -0.79 38
126 // -0.5 15
127 // -0.1 5
128 // 0.1 5
129 // 0.5 10
130 // 1.0 15
131 if (sum2 == sum) {
132 break;
133 }
134 sum = sum2;
135 }
136 return z * (2 * zz * sum - x);
137 }
138
139 /**
140 * Returns {@code log(1 + x) - x}. This function is accurate when
141 * {@code x -> 0}.
142 *
143 * <p>This function uses a Taylor series expansion when x is small
144 * ({@code |x| < 0.01}):
145 *
146 * <pre>
147 * ln(1 + x) - x = -x^2/2 + x^3/3 - x^4/4 + ...
148 * </pre>
149 *
150 * <p>No loop iterations are used as the series is directly expanded
151 * for a set number of terms based on the absolute value of x.
152 *
153 * @param x Value x (assumed to be small)
154 * @param a Absolute value of x
155 * @return {@code log(1 + x) - x}
156 */
157 private static double log1pmxSmall(double x, double a) {
158 // Use a direct Taylor series:
159 // ln(1 + x) = x - x^2/2 + x^3/3 - x^4/4 + ...
160 // Reverse the summation (small to large) for a marginal increase in precision.
161 // To stop the Taylor series the next term must be less than 1 ulp from the
162 // answer.
163 // x^n/n < |log(1+x)-x| * eps
164 // eps = machine epsilon = 2^-53
165 // x^n < |log(1+x)-x| * eps
166 // n < (log(|log(1+x)-x|) + log(eps)) / log(x)
167 // In practice this is a conservative limit.
168
169 final double x2 = x * x;
170
171 if (a < TWO_POW_M53) {
172 // Below machine epsilon. Addition of x^3/3 is not possible.
173 // Subtract from zero to prevent creating -0.0 for x=0.
174 return 0 - x2 / 2;
175 }
176
177 final double x4 = x2 * x2;
178
179 // +/-9.5367431640625e-07: log1pmx = -4.547470617660916e-13 :
180 // -4.5474764000725028e-13
181 // n = 4.69
182 if (a < TWO_POW_M20) {
183 // n=5
184 return x * x4 / 5 -
185 x4 / 4 +
186 x * x2 / 3 -
187 x2 / 2;
188 }
189
190 // +/-2.44140625E-4: log1pmx = -2.9797472637290841e-08 : -2.9807173914456693e-08
191 // n = 6.49
192 if (a < TWO_POW_M12) {
193 // n=7
194 return x * x2 * x4 / 7 -
195 x2 * x4 / 6 +
196 x * x4 / 5 -
197 x4 / 4 +
198 x * x2 / 3 -
199 x2 / 2;
200 }
201
202 // Assume |x| < 2^-6
203 // +/-0.015625: log1pmx = -0.00012081346403474586 : -0.00012335696813916864
204 // n = 10.9974
205
206 // n=11
207 final double x8 = x4 * x4;
208 return x * x2 * x8 / 11 -
209 x2 * x8 / 10 +
210 x * x8 / 9 -
211 x8 / 8 +
212 x * x2 * x4 / 7 -
213 x2 * x4 / 6 +
214 x * x4 / 5 -
215 x4 / 4 +
216 x * x2 / 3 -
217 x2 / 2;
218 }
219 }