1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.statistics.distribution;
18
19 import java.math.BigDecimal;
20 import java.math.MathContext;
21 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.MethodOrderer;
23 import org.junit.jupiter.api.Order;
24 import org.junit.jupiter.api.Test;
25 import org.junit.jupiter.api.TestMethodOrder;
26 import org.junit.jupiter.params.ParameterizedTest;
27 import org.junit.jupiter.params.provider.CsvFileSource;
28 import org.junit.jupiter.params.provider.ValueSource;
29
30
31
32
33 @TestMethodOrder(MethodOrderer.OrderAnnotation.class)
34 class ExtendedPrecisionTest {
35
36 private static final double ROOT2 = Math.sqrt(2.0);
37
38
39
40 private static final double ROOT2PI = Double.parseDouble(SQRT_TWO_PI);
41
42 private static final RMS SQRT2XX_RMS1 = new RMS();
43
44 private static final RMS SQRT2XX_RMS2 = new RMS();
45
46 private static final RMS XSQRT2PI_RMS = new RMS();
47
48 private static final RMS EXPMHXX_RMS1 = new RMS();
49
50 private static final RMS EXPMHXX_RMS2 = new RMS();
51
52
53
54
55
56 private static final class RMS {
57 private double ss;
58 private double max;
59 private int n;
60
61
62
63
64 void add(double x) {
65
66
67 ss += x * x;
68 n++;
69
70 x = Math.abs(x);
71 max = max < x ? x : max;
72 }
73
74
75
76
77
78
79
80
81
82
83
84 double getMax() {
85 return max;
86 }
87
88
89
90
91
92
93
94
95
96 double getRMS() {
97 return Math.sqrt(ss / n);
98 }
99 }
100
101 @Test
102 void testSqrt2PiConstants() {
103 final BigDecimal sqrt2pi = new BigDecimal(SQRT_TWO_PI);
104
105
106
107 final double value = sqrt2pi.doubleValue();
108 final double roundOff = sqrt2pi.subtract(new BigDecimal(value)).doubleValue();
109
110 Assertions.assertEquals(value, value + roundOff, "value + round-off");
111
112 Assertions.assertEquals(value, ExtendedPrecision.SQRT2PI.hi(), "sqrt(2 pi)");
113 Assertions.assertEquals(roundOff, ExtendedPrecision.SQRT2PI.lo(), "sqrt(2 pi) round-off");
114
115 Assertions.assertEquals(value, Math.sqrt(2 * Math.PI), Math.ulp(value), "Math.sqrt(2 pi)");
116 }
117
118 @Test
119 void testSqrt2xxUnderAndOverflow() {
120 final double x = 1.5;
121 final double e = 2.12132034355964257320253308631;
122 Assertions.assertEquals(e, ExtendedPrecision.sqrt2xx(x));
123 for (final int i : new int[] {-1000, -600, -200, 200, 600, 1000}) {
124 final double scale = Math.scalb(1.0, i);
125 final double x1 = x * scale;
126 final double e1 = e * scale;
127 Assertions.assertEquals(e1, ExtendedPrecision.sqrt2xx(x1), () -> Double.toString(x1));
128 }
129 }
130
131 @Test
132 void testSqrt2xxExtremes() {
133
134 Assertions.assertEquals(Double.POSITIVE_INFINITY, ExtendedPrecision.sqrt2xx(Double.MAX_VALUE));
135 Assertions.assertEquals(Double.POSITIVE_INFINITY, ExtendedPrecision.sqrt2xx(Double.POSITIVE_INFINITY));
136 Assertions.assertEquals(0.0, ExtendedPrecision.sqrt2xx(0));
137 Assertions.assertEquals(ROOT2, ExtendedPrecision.sqrt2xx(1));
138 Assertions.assertEquals(Math.sqrt(8), ExtendedPrecision.sqrt2xx(2));
139
140 for (int i = 2; i <= 10; i++) {
141 Assertions.assertEquals(i * Double.MIN_VALUE * Math.sqrt(2), ExtendedPrecision.sqrt2xx(i * Double.MIN_VALUE));
142 }
143
144 Assertions.assertEquals(Double.NaN, ExtendedPrecision.sqrt2xx(Double.NaN));
145
146 Assertions.assertEquals(Double.POSITIVE_INFINITY, ExtendedPrecision.sqrt2xx(-1e300));
147 }
148
149
150
151
152
153
154
155
156
157 @ParameterizedTest
158 @Order(1)
159 @CsvFileSource(resources = "sqrt2xx.csv")
160 void testSqrt2xx(double x, BigDecimal expected) {
161 final double e = expected.doubleValue();
162 Assertions.assertEquals(e, ExtendedPrecision.sqrt2xx(x));
163
164 addError(Math.sqrt(2 * x * x), expected, e, SQRT2XX_RMS1);
165 addError(x * ROOT2, expected, e, SQRT2XX_RMS2);
166 }
167
168 @Test
169 void testSqrt2xxStandardPrecision1() {
170
171 assertPrecision(SQRT2XX_RMS1, 0.9, 0.3);
172 }
173
174 @Test
175 void testSqrt2xxStandardPrecision2() {
176
177 assertPrecision(SQRT2XX_RMS2, 1.3, 0.6);
178 }
179
180 @ParameterizedTest
181 @ValueSource(doubles = {0, 1, Double.MAX_VALUE, Double.POSITIVE_INFINITY, Double.NaN})
182 void testXsqrt2piEdgeCases(double x) {
183 final double expected = x * ROOT2PI;
184 final double actual = ExtendedPrecision.xsqrt2pi(x);
185 Assertions.assertEquals(expected, actual, 1e-15);
186 }
187
188
189
190
191
192
193
194
195
196 @ParameterizedTest
197 @Order(1)
198 @CsvFileSource(resources = "xsqrt2pi.csv")
199 void testXsqrt2pi(double x, BigDecimal expected) {
200 final double e = expected.doubleValue();
201 Assertions.assertEquals(e, ExtendedPrecision.xsqrt2pi(x));
202
203 addError(x * ROOT2PI, expected, e, XSQRT2PI_RMS);
204 }
205
206 @Test
207 void testXsqrt2piPrecision() {
208
209 assertPrecision(XSQRT2PI_RMS, 1.2, 0.6);
210 }
211
212 @ParameterizedTest
213 @ValueSource(doubles = {0, 0.5, 1, 2, 3, 4, 5, 38.5, Double.MAX_VALUE, Double.POSITIVE_INFINITY, Double.NaN})
214 void testExpmhxxEdgeCases(double x) {
215 final double expected = Math.exp(-0.5 * x * x);
216 Assertions.assertEquals(expected, ExtendedPrecision.expmhxx(x));
217 Assertions.assertEquals(expected, ExtendedPrecision.expmhxx(-x));
218 }
219
220
221
222
223
224
225
226
227
228 @ParameterizedTest
229 @Order(1)
230 @CsvFileSource(resources = "expmhxx.csv")
231 void testExpmhxx(double x, BigDecimal expected) {
232 final double e = expected.doubleValue();
233 final double actual = ExtendedPrecision.expmhxx(x);
234 Assertions.assertEquals(e, actual, Math.ulp(e) * 2);
235
236 addError(actual, expected, e, EXPMHXX_RMS1);
237 addError(Math.exp(-0.5 * x * x), expected, e, EXPMHXX_RMS2);
238 }
239
240 @Test
241 void testExpmhxxHighPrecision() {
242
243 assertPrecision(EXPMHXX_RMS1, 1.5, 0.5);
244 }
245
246 @Test
247 void testExpmhxxStandardPrecision() {
248
249 assertPrecision(EXPMHXX_RMS2, 400, 60);
250 }
251
252 private static void assertPrecision(RMS rms, double maxError, double rmsError) {
253 Assertions.assertTrue(rms.getMax() < maxError, () -> "max error: " + rms.getMax());
254 Assertions.assertTrue(rms.getRMS() < rmsError, () -> "rms error: " + rms.getRMS());
255 }
256
257 private static void addError(double z, BigDecimal expected, double e, RMS rms) {
258 double error;
259 if (z == e) {
260 error = 0;
261 } else {
262
263 error = expected.subtract(new BigDecimal(z))
264 .divide(new BigDecimal(Math.ulp(e)), MathContext.DECIMAL64).doubleValue();
265 }
266 rms.add(error);
267 }
268 }