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
18 // (C) Copyright John Maddock 2006.
19 // Use, modification and distribution are subject to the
20 // Boost Software License, Version 1.0. (See accompanying file
21 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
22
23 package org.apache.commons.numbers.gamma;
24
25 /**
26 * Math functions used by the Boost functions.
27 *
28 * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
29 * {@code c++} implementations in {@code <boost/math/special_functions/>}.
30 * All work is copyright John Maddock 2006 and subject to the Boost Software License.
31 */
32 final class BoostMath {
33 /** Private constructor. */
34 private BoostMath() {
35 // intentionally empty.
36 }
37
38 /**
39 * Returns {@code pow(x, y) - 1}. This function is accurate when {@code x -> 1} or {@code y} is
40 * small.
41 *
42 * <p>Adapted from {@code boost/math/special_functions/powm1.hpp}. Explicit handling of
43 * edges cases (overflow, domain error) using the policy has been removed.
44 *
45 * @param x the x
46 * @param y the y
47 * @return {@code pow(x, y) - 1}
48 */
49 static double powm1(double x, double y) {
50 if (x > 0) {
51 // Check for small y or x close to 1.
52 // Require term < 0.5
53 // => log(x) * y < 0.5
54 // Assume log(x) ~ (x - 1) [true when x is close to 1]
55 // => |(x-1) * y| < 0.5
56
57 if (Math.abs(y * (x - 1)) < 0.5 || Math.abs(y) < 0.2) {
58 // We don't have any good/quick approximation for log(x) * y
59 // so just try it and see:
60 final double l = y * Math.log(x);
61 if (l < 0.5) {
62 return Math.expm1(l);
63 }
64 // fall through....
65 }
66 } else if (x < 0 &&
67 // y had better be an integer:
68 // x is negative.
69 // pow(x, y) only allowed if y is an integer.
70 // if y is even then we can invert non-zero finite x.
71 Math.rint(y * 0.5) == y * 0.5) {
72 return powm1(-x, y);
73 }
74 return Math.pow(x, y) - 1;
75 }
76 }