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.statistics.distribution;
18
19 import org.apache.commons.numbers.core.DD;
20
21 /**
22 * Computes extended precision floating-point operations.
23 *
24 * <p>Extended precision computation is delegated to the {@link DD} class. The methods here
25 * verify the arguments to the computations will not overflow.
26 */
27 final class ExtendedPrecision {
28 /** sqrt(2 pi) as a double-double number.
29 * Divided into two parts from the value sqrt(2 pi) computed to 64 decimal digits. */
30 static final DD SQRT2PI = DD.ofSum(2.5066282746310007, -1.8328579980459167e-16);
31
32 /** Threshold for a big number that may overflow when squared. 2^500. */
33 private static final double BIG = 0x1.0p500;
34 /** Threshold for a small number that may underflow when squared. 2^-500. */
35 private static final double SMALL = 0x1.0p-500;
36 /** Scale up by 2^600. */
37 private static final double SCALE_UP = 0x1.0p600;
38 /** Scale down by 2^600. */
39 private static final double SCALE_DOWN = 0x1.0p-600;
40 /** X squared value where {@code exp(-0.5*x*x)} cannot increase accuracy using the round-off
41 * from x squared. */
42 private static final int EXP_M_HALF_XX_MIN_VALUE = 2;
43 /** Approximate x squared value where {@code exp(-0.5*x*x) == 0}. This is above
44 * {@code -2 * ln(2^-1074)} due to rounding performed within the exp function. */
45 private static final int EXP_M_HALF_XX_MAX_VALUE = 1491;
46
47 /** No instances. */
48 private ExtendedPrecision() {}
49
50 /**
51 * Multiply the term by sqrt(2 pi).
52 *
53 * @param x Value (assumed to be positive)
54 * @return x * sqrt(2 pi)
55 */
56 static double xsqrt2pi(double x) {
57 // Note: Do not convert x to absolute for this use case
58 if (x > BIG) {
59 if (x == Double.POSITIVE_INFINITY) {
60 return Double.POSITIVE_INFINITY;
61 }
62 return computeXsqrt2pi(x * SCALE_DOWN) * SCALE_UP;
63 } else if (x < SMALL) {
64 // Note: Ignore possible zero for this use case
65 return computeXsqrt2pi(x * SCALE_UP) * SCALE_DOWN;
66 } else {
67 return computeXsqrt2pi(x);
68 }
69 }
70
71 /**
72 * Compute {@code a * sqrt(2 * pi)}.
73 *
74 * @param a Value
75 * @return the result
76 */
77 private static double computeXsqrt2pi(double a) {
78 return SQRT2PI.multiply(a).hi();
79 }
80
81 /**
82 * Compute {@code sqrt(2 * x * x)}.
83 *
84 * <p>The result is computed using a high precision computation of
85 * {@code sqrt(2 * x * x)} avoiding underflow or overflow of {@code x}
86 * squared.
87 *
88 * @param x Value (assumed to be positive)
89 * @return {@code sqrt(2 * x * x)}
90 */
91 static double sqrt2xx(double x) {
92 // Note: Do not convert x to absolute for this use case
93 if (x > BIG) {
94 if (x == Double.POSITIVE_INFINITY) {
95 return Double.POSITIVE_INFINITY;
96 }
97 return computeSqrt2aa(x * SCALE_DOWN) * SCALE_UP;
98 } else if (x < SMALL) {
99 // Note: Ignore possible zero for this use case
100 return computeSqrt2aa(x * SCALE_UP) * SCALE_DOWN;
101 } else {
102 return computeSqrt2aa(x);
103 }
104 }
105
106 /**
107 * Compute {@code sqrt(2 * a * a)}.
108 *
109 * @param a Value
110 * @return the result
111 */
112 private static double computeSqrt2aa(double a) {
113 return DD.ofProduct(2 * a, a).sqrt().hi();
114 }
115
116 /**
117 * Compute {@code exp(-0.5*x*x)} with high accuracy. This is performed using information in the
118 * round-off from {@code x*x}.
119 *
120 * <p>This is accurate at large x to 1 ulp until exp(-0.5*x*x) is close to sub-normal. For very
121 * small exp(-0.5*x*x) the adjustment is sub-normal and bits can be lost in the adjustment for a
122 * max observed error of {@code < 2} ulp.
123 *
124 * <p>At small x the accuracy cannot be improved over using exp(-0.5*x*x). This occurs at
125 * {@code x <= sqrt(2)}.
126 *
127 * @param x Value
128 * @return exp(-0.5*x*x)
129 * @see <a href="https://issues.apache.org/jira/browse/STATISTICS-52">STATISTICS-52</a>
130 */
131 static double expmhxx(double x) {
132 final double z = x * x;
133 if (z <= EXP_M_HALF_XX_MIN_VALUE) {
134 return Math.exp(-0.5 * z);
135 } else if (z >= EXP_M_HALF_XX_MAX_VALUE) {
136 // exp(-745.5) == 0
137 return 0;
138 }
139 final DD x2 = DD.ofSquare(x);
140 return expxx(-0.5 * x2.hi(), -0.5 * x2.lo());
141 }
142
143 /**
144 * Compute {@code exp(a+b)} with high accuracy assuming {@code a+b = a}.
145 *
146 * <p>This is accurate at large positive a to 1 ulp. If a is negative and exp(a) is close to
147 * sub-normal a bit of precision may be lost when adjusting result as the adjustment is sub-normal
148 * (max observed error {@code < 2} ulp). For the use case of multiplication of a number less than
149 * 1 by exp(-x*x), a = -x*x, the result will be sub-normal and the rounding error is lost.
150 *
151 * <p>At small |a| the accuracy cannot be improved over using exp(a) as the round-off is too small
152 * to create terms that can adjust the standard result by more than 0.5 ulp. This occurs at
153 * {@code |a| <= 1}.
154 *
155 * @param a High bits of a split number
156 * @param b Low bits of a split number
157 * @return exp(a+b)
158 * @see <a href="https://issues.apache.org/jira/projects/NUMBERS/issues/NUMBERS-177">
159 * Numbers-177: Accurate scaling by exp(z*z)</a>
160 */
161 private static double expxx(double a, double b) {
162 // exp(a+b) = exp(a) * exp(b)
163 // = exp(a) * (exp(b) - 1) + exp(a)
164 // Assuming:
165 // 1. -746 < a < 710 for no under/overflow of exp(a)
166 // 2. a+b = a
167 // As b -> 0 then exp(b) -> 1; expm1(b) -> b
168 // The round-off b is limited to ~ 0.5 * ulp(746) ~ 5.68e-14
169 // and we can use an approximation for expm1 (x/1! + x^2/2! + ...)
170 // The second term is required for the expm1 result but the
171 // bits are not significant to change the following sum with exp(a)
172
173 final double ea = Math.exp(a);
174 // b ~ expm1(b)
175 return ea * b + ea;
176 }
177 }