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 package org.apache.commons.numbers.examples.jmh.core;
19
20 import java.util.concurrent.TimeUnit;
21 import java.util.function.DoubleBinaryOperator;
22 import java.util.function.DoublePredicate;
23 import java.util.function.DoubleUnaryOperator;
24
25 import org.apache.commons.rng.UniformRandomProvider;
26 import org.apache.commons.rng.simple.RandomSource;
27 import org.openjdk.jmh.annotations.Benchmark;
28 import org.openjdk.jmh.annotations.BenchmarkMode;
29 import org.openjdk.jmh.annotations.Fork;
30 import org.openjdk.jmh.annotations.Measurement;
31 import org.openjdk.jmh.annotations.Mode;
32 import org.openjdk.jmh.annotations.OutputTimeUnit;
33 import org.openjdk.jmh.annotations.Param;
34 import org.openjdk.jmh.annotations.Scope;
35 import org.openjdk.jmh.annotations.Setup;
36 import org.openjdk.jmh.annotations.State;
37 import org.openjdk.jmh.annotations.Warmup;
38 import org.openjdk.jmh.infra.Blackhole;
39
40 /**
41 * Executes a benchmark to measure the speed of operations in the {@link LinearCombination} class.
42 * Benchmarks focus on the split of a double value into high and low parts.
43 */
44 @BenchmarkMode(Mode.AverageTime)
45 @OutputTimeUnit(TimeUnit.NANOSECONDS)
46 @Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
47 @Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
48 @State(Scope.Benchmark)
49 @Fork(value = 1, jvmArgs = {"-server", "-Xms512M", "-Xmx512M"})
50 public class DoubleSplitPerformance {
51 /** The mask for the sign bit and the mantissa. */
52 private static final long SIGN_MATISSA_MASK = 0x800f_ffff_ffff_ffffL;
53
54 /**
55 * The multiplier used to split the double value into high and low parts. From
56 * Dekker (1971): "The constant should be chosen equal to 2^(p - p/2) + 1,
57 * where p is the number of binary digits in the mantissa". Here p is 53
58 * and the multiplier is {@code 2^27 + 1}.
59 */
60 private static final double MULTIPLIER = 1.34217729E8;
61
62 /** The upper limit above which a number may overflow during the split into a high part.
63 * Assuming the multiplier is above 2^27 and the maximum exponent is 1023 then a safe
64 * limit is a value with an exponent of (1023 - 27) = 2^996. */
65 private static final double SAFE_UPPER = 0x1.0p996;
66
67 /** The scale to use when down-scaling during a split into a high part.
68 * This must be smaller than the inverse of the multiplier and a power of 2 for exact scaling. */
69 private static final double DOWN_SCALE = 0x1.0p-30;
70
71 /** The scale to use when re-scaling during a split into a high part.
72 * This is the inverse of {@link #DOWN_SCALE}. */
73 private static final double UP_SCALE = 0x1.0p30;
74
75 /** The mask to zero the lower 27-bits of a long . */
76 private static final long ZERO_LOWER_27_BITS = 0xffff_ffff_f800_0000L;
77
78 /** Constant to no method. */
79 private static final String NONE = "none";
80
81 /**
82 * The numbers to split.
83 */
84 @State(Scope.Benchmark)
85 public static class Numbers {
86 /** The exponent for small numbers. */
87 private static final long EXP_SMALL = Double.doubleToRawLongBits(1.0);
88 /** The exponent for big numbers. */
89 private static final long EXP_BIG = Double.doubleToRawLongBits(SAFE_UPPER);
90
91 /**
92 * The count of numbers.
93 */
94 @Param({"10000"})
95 private int size;
96
97 /**
98 * The fraction of small numbers.
99 *
100 * <p>Note: The split method may employ multiplications.
101 * Big numbers are edge cases that would cause overflow in multiplications.
102 */
103 @Param({"1", "0.999", "0.99", "0.9"})
104 private double edge;
105
106 /** Numbers. */
107 private double[] a;
108
109 /**
110 * Gets the factors.
111 *
112 * @return Factors.
113 */
114 public double[] getNumbers() {
115 return a;
116 }
117
118 /**
119 * Create the factors.
120 */
121 @Setup
122 public void setup() {
123 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create();
124 a = new double[size];
125 for (int i = 0; i < size; i++) {
126 long bits = rng.nextLong() & SIGN_MATISSA_MASK;
127 // The exponent will either be small or big
128 if (rng.nextDouble() < edge) {
129 bits |= EXP_SMALL;
130 } else {
131 bits |= EXP_BIG;
132 }
133 a[i] = Double.longBitsToDouble(bits);
134 }
135 }
136 }
137
138 /**
139 * The factors to multiply.
140 */
141 @State(Scope.Benchmark)
142 public static class BiFactors {
143 /** The exponent for small numbers. */
144 private static final long EXP_SMALL = Double.doubleToRawLongBits(1.0);
145
146 /**
147 * The count of products.
148 */
149 @Param({"5000"})
150 private int size;
151
152 /**
153 * The exponent for big numbers.
154 * Two big numbers multiplied together should overflow.
155 *
156 * <p>Note: If this is set below the safe upper for Dekker's split then a method
157 * that computes the split of the two factors independently and then does the multiply
158 * may create split parts that will overflow during multiplication.
159 */
160 @Param({"600", "1000", "1023"})
161 private int exp;
162
163 /**
164 * The fraction of small numbers.
165 *
166 * <p>Note: The split numbers are used in multiplications.
167 * It is unlikely that many numbers will be larger than the upper limit.
168 * These numbers are edge cases that would cause overflow in multiplications if
169 * the other number is anywhere close to the same magnitude.
170 */
171 @Param({"1", "0.95", "0.9"})
172 private double edge;
173
174 /** Factors. */
175 private double[] a;
176
177 /**
178 * Gets the factors to be multiplied as pairs. The array length will be even.
179 *
180 * @return Factors.
181 */
182 public double[] getFactors() {
183 return a;
184 }
185
186 /**
187 * Create the factors.
188 */
189 @Setup
190 public void setup() {
191 // Validate the big exponent
192 final double d = Math.scalb(1.0, exp);
193 assert Double.isInfinite(d * d) : "Product of big numbers does not overflow";
194 final long expBig = Double.doubleToRawLongBits(d);
195
196 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create();
197 a = new double[size * 2];
198 for (int i = 0; i < a.length; i++) {
199 long bits = rng.nextLong() & SIGN_MATISSA_MASK;
200 // The exponent will either be small or big
201 if (rng.nextDouble() < edge) {
202 bits |= EXP_SMALL;
203 } else {
204 bits |= expBig;
205 }
206 a[i] = Double.longBitsToDouble(bits);
207 }
208 }
209 }
210
211 /**
212 * The numbers to test to determine if they are not normal.
213 */
214 @State(Scope.Benchmark)
215 public static class NonNormalNumbers {
216 /** Non-normal positive numbers. */
217 private static final double[] NON_NORMAL =
218 {Double.POSITIVE_INFINITY, Double.NaN, Double.MIN_NORMAL};
219
220 /**
221 * The count of numbers.
222 */
223 @Param({"10000"})
224 private int size;
225
226 /**
227 * The fraction of non-normal factors.
228 */
229 @Param({"1", "0.999", "0.99", "0.9"})
230 private double edge;
231
232 /** Numbers. */
233 private double[] a;
234
235 /**
236 * Gets the factors.
237 *
238 * @return Factors.
239 */
240 public double[] getFactors() {
241 return a;
242 }
243
244 /**
245 * Create the factors.
246 */
247 @Setup
248 public void setup() {
249 final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_1024_PP.create();
250 a = new double[size];
251 for (int i = 0; i < size; i++) {
252 // Value in (-1, 1)
253 double value = rng.nextDouble() * (rng.nextBoolean() ? -1 : 1);
254 // The number will either be small or non-normal
255 if (rng.nextDouble() < edge) {
256 value *= NON_NORMAL[rng.nextInt(NON_NORMAL.length)];
257 }
258 a[i] = value;
259 }
260 }
261 }
262
263 /**
264 * The split method.
265 */
266 @State(Scope.Benchmark)
267 public static class SplitMethod {
268 /**
269 * The name of the method.
270 */
271 @Param({NONE, "dekker", "dekkerAbs", "dekkerRaw", "bits"})
272 private String name;
273
274 /** The function. */
275 private DoubleUnaryOperator fun;
276
277 /**
278 * Gets the function.
279 *
280 * @return the function
281 */
282 public DoubleUnaryOperator getFunction() {
283 return fun;
284 }
285
286 /**
287 * Create the function.
288 */
289 @Setup
290 public void setup() {
291 if (NONE.equals(name)) {
292 fun = a -> a;
293 } else if ("dekker".equals(name)) {
294 fun = DoubleSplitPerformance::splitDekker;
295 } else if ("dekkerAbs".equals(name)) {
296 fun = DoubleSplitPerformance::splitDekkerAbs;
297 } else if ("dekkerRaw".equals(name)) {
298 fun = DoubleSplitPerformance::splitDekkerRaw;
299 } else if ("bits".equals(name)) {
300 fun = DoubleSplitPerformance::splitBits;
301 } else {
302 throw new IllegalStateException("Unknown split method: " + name);
303 }
304 }
305 }
306
307 /**
308 * The method to test for a non-normal number.
309 */
310 @State(Scope.Benchmark)
311 public static class NonNormalMethod {
312 /**
313 * The name of the method.
314 */
315 @Param({NONE, "if", "exponent", "exponent2"})
316 private String name;
317
318 /** The function. */
319 private DoublePredicate fun;
320
321 /**
322 * Gets the function.
323 *
324 * @return the function
325 */
326 public DoublePredicate getFunction() {
327 return fun;
328 }
329
330 /**
331 * Create the function.
332 */
333 @Setup
334 public void setup() {
335 if (NONE.equals(name)) {
336 fun = a -> false;
337 } else if ("if".equals(name)) {
338 fun = DoubleSplitPerformance::isNotNormalIf;
339 } else if ("exponent".equals(name)) {
340 fun = DoubleSplitPerformance::isNotNormalExponent;
341 } else if ("exponent2".equals(name)) {
342 fun = DoubleSplitPerformance::isNotNormalExponent2;
343 } else {
344 throw new IllegalStateException("Unknown is non-normal method: " + name);
345 }
346 }
347 }
348
349 /**
350 * The method to compute the product round-off.
351 */
352 @State(Scope.Benchmark)
353 public static class RoundoffMethod {
354 /**
355 * The name of the method.
356 */
357 @Param({NONE, "multiply", "multiplyUnscaled",
358 "productLow", "productLow1", "productLow2", "productLow3", "productLowSplit",
359 "productLowUnscaled"})
360 private String name;
361
362 /** The function. */
363 private DoubleBinaryOperator fun;
364
365 /**
366 * Gets the function.
367 *
368 * @return the function
369 */
370 public DoubleBinaryOperator getFunction() {
371 return fun;
372 }
373
374 /**
375 * Create the function.
376 */
377 @Setup
378 public void setup() {
379 if (NONE.equals(name)) {
380 // No actually the round-off but x*y - x*y will be optimised away so this
381 // captures the multiply overhead.
382 fun = (x, y) -> x * y;
383 } else if ("multiply".equals(name)) {
384 final DoublePrecision.Quad result = new DoublePrecision.Quad();
385 fun = (x, y) -> {
386 DoublePrecision.multiply(x, y, result);
387 return result.lo;
388 };
389 } else if ("productLow".equals(name)) {
390 fun = (x, y) -> DoublePrecision.productLow(x, y, x * y);
391 } else if ("productLow1".equals(name)) {
392 fun = (x, y) -> DoublePrecision.productLow1(x, y, x * y);
393 } else if ("productLow2".equals(name)) {
394 fun = (x, y) -> DoublePrecision.productLow2(x, y, x * y);
395 } else if ("productLow3".equals(name)) {
396 fun = (x, y) -> DoublePrecision.productLow3(x, y, x * y);
397 } else if ("productLowSplit".equals(name)) {
398 fun = (x, y) -> DoublePrecision.productLowSplit(x, y, x * y);
399 } else if ("multiplyUnscaled".equals(name)) {
400 final DoublePrecision.Quad result = new DoublePrecision.Quad();
401 fun = (x, y) -> {
402 DoublePrecision.multiplyUnscaled(x, y, result);
403 return result.lo;
404 };
405 } else if ("productLowUnscaled".equals(name)) {
406 fun = (x, y) -> DoublePrecision.productLowUnscaled(x, y, x * y);
407 } else {
408 throw new IllegalStateException("Unknown round-off method: " + name);
409 }
410 }
411 }
412
413 /**
414 * Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) creates
415 * a big value from which to derive the two split parts.
416 * <pre>
417 * c = (2^s + 1) * a
418 * a_big = c - a
419 * a_hi = c - a_big
420 * a_lo = a - a_hi
421 * a = a_hi + a_lo
422 * </pre>
423 *
424 * <p>The multiplicand allows a p-bit value to be split into
425 * (p-s)-bit value {@code a_hi} and a non-overlapping (s-1)-bit value {@code a_lo}.
426 * Combined they have (p-1) bits of significand but the sign bit of {@code a_lo}
427 * contains a bit of information. The constant is chosen so that s is ceil(p/2) where
428 * the precision p for a double is 53-bits (1-bit of the mantissa is assumed to be
429 * 1 for a non sub-normal number) and s is 27.
430 *
431 * @param value Value.
432 * @return the high part of the value.
433 */
434 private static double splitDekker(double value) {
435 // Avoid overflow
436 if (value >= SAFE_UPPER || value <= -SAFE_UPPER) {
437 // Do scaling.
438 final double x = value * DOWN_SCALE;
439 final double c = MULTIPLIER * x;
440 final double hi = (c - (c - x)) * UP_SCALE;
441 if (Double.isInfinite(hi)) {
442 // Number is too large.
443 // This occurs if value is infinite or close to Double.MAX_VALUE.
444 // Note that Dekker's split creates an approximating 26-bit number which may
445 // have an exponent 1 greater than the input value. This will overflow if the
446 // exponent is already +1023. Revert to the raw upper 26 bits of the 53-bit
447 // mantissa (including the assumed leading 1 bit). This conversion will result in
448 // the low part being a 27-bit significand and the potential loss of bits during
449 // addition and multiplication. (Contrast to the Dekker split which creates two
450 // 26-bit numbers with a bit of information moved to the sign of low.)
451 // The conversion will maintain Infinite in the high part where the resulting
452 // low part a_lo = a - a_hi = inf - inf = NaN.
453 return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
454 }
455 return hi;
456 }
457 // normal conversion
458 final double c = MULTIPLIER * value;
459 return c - (c - value);
460 }
461
462 /**
463 * Implement Dekker's method to split a value into two parts.
464 * This is as per {@link #splitDekker(double)} but uses a {@link Math#abs(double)} in the
465 * condition statement.
466 *
467 * @param value Value.
468 * @return the high part of the value.
469 */
470 private static double splitDekkerAbs(double value) {
471 if (Math.abs(value) >= SAFE_UPPER) {
472 final double x = value * DOWN_SCALE;
473 final double c = MULTIPLIER * x;
474 final double hi = (c - (c - x)) * UP_SCALE;
475 if (Double.isInfinite(hi)) {
476 return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
477 }
478 return hi;
479 }
480 final double c = MULTIPLIER * value;
481 return c - (c - value);
482 }
483
484 /**
485 * Implement Dekker's method to split a value into two parts.
486 * This is as per {@link #splitDekker(double)} but has no overflow protection.
487 *
488 * @param value Value.
489 * @return the high part of the value.
490 */
491 private static double splitDekkerRaw(double value) {
492 final double c = MULTIPLIER * value;
493 return c - (c - value);
494 }
495
496 /**
497 * Implement a split using the upper and lower raw bits from the value.
498 *
499 * @param value Value.
500 * @return the high part of the value.
501 */
502 private static double splitBits(double value) {
503 return Double.longBitsToDouble(Double.doubleToRawLongBits(value) & ZERO_LOWER_27_BITS);
504 }
505
506 /**
507 * Checks if the number is not normal. This is functionally equivalent to:
508 * <pre>
509 * </pre>
510 *
511 * @param a The value.
512 * @return true if the value is not normal
513 */
514 private static boolean isNotNormalIf(double a) {
515 final double abs = Math.abs(a);
516 return abs <= Double.MIN_NORMAL || !(abs <= Double.MAX_VALUE);
517 }
518
519 /**
520 * Checks if the number is not normal.
521 *
522 * @param a The value.
523 * @return true if the value is not normal
524 */
525 private static boolean isNotNormalExponent(double a) {
526 // Sub-normal numbers have a biased exponent of 0.
527 // Inf/NaN numbers have a biased exponent of 2047.
528 // Catch both cases by extracting the raw exponent, subtracting 1
529 // and make unsigned. 0 will underflow to a large value.
530 final int baisedExponent = ((int) (Double.doubleToRawLongBits(a) >>> 52)) & 0x7ff;
531 return ((baisedExponent - 1) & 0xffff) >= 2046;
532 }
533
534 /**
535 * Checks if the number is not normal.
536 *
537 * @param a The value.
538 * @return true if the value is not normal
539 */
540 private static boolean isNotNormalExponent2(double a) {
541 // Sub-normal numbers have a biased exponent of 0.
542 // Inf/NaN numbers have a biased exponent of 2047.
543 // Catch both cases by extracting the raw exponent, subtracting 1
544 // and compare unsigned (so 0 underflows to a large value).
545 final int baisedExponent = ((int) (Double.doubleToRawLongBits(a) >>> 52)) & 0x7ff;
546 // Adding int min value is equal to compare unsigned
547 return baisedExponent + Integer.MIN_VALUE - 1 >= 2046 + Integer.MIN_VALUE;
548 }
549
550 // Benchmark methods.
551
552 /**
553 * Benchmark extracting the high part of the split number.
554 *
555 * @param numbers Factors.
556 * @param bh Data sink.
557 * @param method Split method
558 */
559 @Benchmark
560 public void high(Numbers numbers, Blackhole bh, SplitMethod method) {
561 final DoubleUnaryOperator fun = method.getFunction();
562 final double[] a = numbers.getNumbers();
563 for (int i = 0; i < a.length; i++) {
564 bh.consume(fun.applyAsDouble(a[i]));
565 }
566 }
567
568 /**
569 * Benchmark extracting the low part of the split number.
570 *
571 * @param numbers Factors.
572 * @param bh Data sink.
573 * @param method Split method.
574 */
575 @Benchmark
576 public void low(Numbers numbers, Blackhole bh, SplitMethod method) {
577 final DoubleUnaryOperator fun = method.getFunction();
578 final double[] a = numbers.getNumbers();
579 for (int i = 0; i < a.length; i++) {
580 bh.consume(a[i] - fun.applyAsDouble(a[i]));
581 }
582 }
583
584 /**
585 * Benchmark testing if a number is non-normal.
586 *
587 * @param numbers Factors.
588 * @param bh Data sink.
589 * @param method Split method.
590 */
591 @Benchmark
592 public void nonNormal(NonNormalNumbers numbers, Blackhole bh, NonNormalMethod method) {
593 final DoublePredicate fun = method.getFunction();
594 final double[] a = numbers.getFactors();
595 for (int i = 0; i < a.length; i++) {
596 bh.consume(fun.test(a[i]));
597 }
598 }
599
600 /**
601 * Benchmark extracting the round-off from the product of two numbers.
602 *
603 * @param factors Factors.
604 * @param bh Data sink.
605 * @param method Round-off method.
606 */
607 @Benchmark
608 public void productLow(BiFactors factors, Blackhole bh, RoundoffMethod method) {
609 final DoubleBinaryOperator fun = method.getFunction();
610 final double[] a = factors.getFactors();
611 for (int i = 0; i < a.length; i += 2) {
612 bh.consume(fun.applyAsDouble(a[i], a[i + 1]));
613 }
614 }
615 }