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.descriptive;
18
19 import java.math.BigDecimal;
20 import java.math.BigInteger;
21 import org.apache.commons.numbers.core.DD;
22
23 /**
24 * Support class for integer math.
25 *
26 * @since 1.1
27 */
28 final class IntMath {
29 /** Mask for the lower 32-bits of a long. */
30 private static final long MASK32 = 0xffff_ffffL;
31 /** Mask for the lower 52-bits of a long. */
32 private static final long MASK52 = 0xf_ffff_ffff_ffffL;
33 /** Bias offset for the exponent of a double. */
34 private static final int EXP_BIAS = 1023;
35 /** Shift for the exponent of a double. */
36 private static final int EXP_SHIFT = 52;
37 /** 0.5. */
38 private static final double HALF = 0.5;
39 /** 2^53. */
40 private static final long TWO_POW_53 = 1L << 53;
41
42 /** No instances. */
43 private IntMath() {}
44
45 /**
46 * Square the values as if an unsigned 64-bit long to produce the high 64-bits
47 * of the 128-bit unsigned result.
48 *a
49 * <p>This method computes the equivalent of:
50 * <pre>{@code
51 * Math.multiplyHigh(x, x)
52 * Math.unsignedMultiplyHigh(x, x) - (((x >> 63) & x) << 1)
53 * }</pre>
54 *
55 * <p>Note: The method {@code Math.multiplyHigh} was added in JDK 9
56 * and should be used as above when the source code targets Java 11
57 * to exploit the intrinsic method.
58 *
59 * <p>Note: The method uses the unsigned multiplication. When the input is negative
60 * it can be adjusted to the signed result by subtracting the argument twice from the
61 * result.
62 *
63 * @param x Value
64 * @return the high 64-bits of the 128-bit result
65 */
66 static long squareHigh(long x) {
67 // Computation is based on the following observation about the upper (a and x)
68 // and lower (b and y) bits of unsigned big-endian integers:
69 // ab * xy
70 // = b * y
71 // + b * x0
72 // + a0 * y
73 // + a0 * x0
74 // = b * y
75 // + b * x * 2^32
76 // + a * y * 2^32
77 // + a * x * 2^64
78 //
79 // Summation using a character for each byte:
80 //
81 // byby byby
82 // + bxbx bxbx 0000
83 // + ayay ayay 0000
84 // + axax axax 0000 0000
85 //
86 // The summation can be rearranged to ensure no overflow given
87 // that the result of two unsigned 32-bit integers multiplied together
88 // plus two full 32-bit integers cannot overflow 64 bits:
89 // > long x = (1L << 32) - 1
90 // > x * x + x + x == -1 (all bits set, no overflow)
91 //
92 // The carry is a composed intermediate which will never overflow:
93 //
94 // byby byby
95 // + bxbx 0000
96 // + ayay ayay 0000
97 //
98 // + bxbx 0000 0000
99 // + axax axax 0000 0000
100
101 final long a = x >>> 32;
102 final long b = x & MASK32;
103
104 final long aa = a * a;
105 final long ab = a * b;
106 final long bb = b * b;
107
108 // Cannot overflow
109 final long carry = (bb >>> 32) +
110 (ab & MASK32) +
111 ab;
112 // Note:
113 // low = (carry << 32) | (bb & MASK32)
114 // Benchmarking shows outputting low to a long[] output argument
115 // has no benefit over computing 'low = value * value' separately.
116
117 final long hi = (ab >>> 32) + (carry >>> 32) + aa;
118 // Adjust to the signed result:
119 // if x < 0:
120 // hi - 2 * x
121 return hi - (((x >> 63) & x) << 1);
122 }
123
124 /**
125 * Multiply the two values as if unsigned 64-bit longs to produce the high 64-bits
126 * of the 128-bit unsigned result.
127 *
128 * <p>This method computes the equivalent of:
129 * <pre>{@code
130 * Math.multiplyHigh(a, b) + ((a >> 63) & b) + ((b >> 63) & a)
131 * }</pre>
132 *
133 * <p>Note: The method {@code Math.multiplyHigh} was added in JDK 9
134 * and should be used as above when the source code targets Java 11
135 * to exploit the intrinsic method.
136 *
137 * <p>Note: The method {@code Math.unsignedMultiplyHigh} was added in JDK 18
138 * and should be used when the source code target allows.
139 *
140 * <p>Taken from {@code o.a.c.rng.core.source64.LXMSupport}.
141 *
142 * @param value1 the first value
143 * @param value2 the second value
144 * @return the high 64-bits of the 128-bit result
145 */
146 static long unsignedMultiplyHigh(long value1, long value2) {
147 // Computation is based on the following observation about the upper (a and x)
148 // and lower (b and y) bits of unsigned big-endian integers:
149 // ab * xy
150 // = b * y
151 // + b * x0
152 // + a0 * y
153 // + a0 * x0
154 // = b * y
155 // + b * x * 2^32
156 // + a * y * 2^32
157 // + a * x * 2^64
158 //
159 // Summation using a character for each byte:
160 //
161 // byby byby
162 // + bxbx bxbx 0000
163 // + ayay ayay 0000
164 // + axax axax 0000 0000
165 //
166 // The summation can be rearranged to ensure no overflow given
167 // that the result of two unsigned 32-bit integers multiplied together
168 // plus two full 32-bit integers cannot overflow 64 bits:
169 // > long x = (1L << 32) - 1
170 // > x * x + x + x == -1 (all bits set, no overflow)
171 //
172 // The carry is a composed intermediate which will never overflow:
173 //
174 // byby byby
175 // + bxbx 0000
176 // + ayay ayay 0000
177 //
178 // + bxbx 0000 0000
179 // + axax axax 0000 0000
180
181 final long a = value1 >>> 32;
182 final long b = value1 & MASK32;
183 final long x = value2 >>> 32;
184 final long y = value2 & MASK32;
185
186 final long by = b * y;
187 final long bx = b * x;
188 final long ay = a * y;
189 final long ax = a * x;
190
191 // Cannot overflow
192 final long carry = (by >>> 32) +
193 (bx & MASK32) +
194 ay;
195 // Note:
196 // low = (carry << 32) | (by & INT_TO_UNSIGNED_BYTE_MASK)
197 // Benchmarking shows outputting low to a long[] output argument
198 // has no benefit over computing 'low = value1 * value2' separately.
199
200 return (bx >>> 32) + (carry >>> 32) + ax;
201 }
202
203 /**
204 * Multiply the arguments as if unsigned integers to a {@code double} result.
205 *
206 * @param x Value.
207 * @param y Value.
208 * @return the double
209 */
210 static double unsignedMultiplyToDouble(long x, long y) {
211 final long lo = x * y;
212 // Fast case: check the arguments cannot overflow a long.
213 // This is true if neither has the upper 33-bits set.
214 if (((x | y) >>> 31) == 0) {
215 // Implicit conversion to a double.
216 return lo;
217 }
218 return uint128ToDouble(unsignedMultiplyHigh(x, y), lo);
219 }
220
221 /**
222 * Convert an unsigned 128-bit integer to a {@code double}.
223 *
224 * @param hi High 64-bits.
225 * @param lo Low 64-bits.
226 * @return the double
227 */
228 static double uint128ToDouble(long hi, long lo) {
229 // Require the representation:
230 // 2^exp * mantissa / 2^53
231 // The mantissa has an implied leading 1-bit.
232
233 // We have the mantissa final bit as xxx0 or xxx1.
234 // To perform correct rounding we maintain the 54-th bit (a) and
235 // a check bit (b) of remaining bits.
236 // Cases:
237 // xxx0 00 - round-down [1]
238 // xxx0 0b - round-down [1]
239 // xxx0 a0 - half-even, round-down [4]
240 // xxx0 ab - round-up [2]
241 // xxx1 00 - round-down [1]
242 // xxx1 0b - round-down [1]
243 // xxx1 a0 - half-even, round-up [3]
244 // xxx1 ab - round-up [2]
245 // [1] If the 54-th bit is 0 always round-down.
246 // [2] Otherwise round-up if the check bit is set or
247 // [3] the final bit is odd (half-even rounding up).
248 // [4] half-even rounding down.
249
250 if (hi == 0) {
251 // If lo is a 63-bit result then we are done
252 if (lo >= 0) {
253 return lo;
254 }
255 // Create a 63-bit number with a sticky bit for rounding, rescale the result
256 return 2 * (double) ((lo >>> 1) | (lo & 0x1));
257 }
258
259 // Initially we create the most significant 64-bits.
260 final int shift = Long.numberOfLeadingZeros(hi);
261 // Shift the high bits and add trailing low bits.
262 // The mask is for the bits from low that are *not* used.
263 // Flipping the mask obtains the bits we concatenate
264 // after shifting (64 - shift).
265 final long maskLow = -1L >>> shift;
266 long bits64 = (hi << shift) | ((lo & ~maskLow) >>> -shift);
267 // exponent for 2^exp is the index of the highest bit in the 128 bit integer
268 final int exp = 127 - shift;
269 // Some of the low bits are lost. If non-zero set
270 // a sticky bit for rounding.
271 bits64 |= (lo & maskLow) == 0 ? 0 : 1;
272
273 // We have a 64-bit unsigned fraction magnitude and an exponent.
274 // This must be converted to a IEEE double by mapping the fraction to a base of 2^53.
275
276 // Create the 53-bit mantissa without the implicit 1-bit
277 long bits = (bits64 >>> 11) & MASK52;
278 // Extract 54-th bit and a sticky bit
279 final long a = (bits64 >>> 10) & 0x1;
280 final long b = (bits64 << 54) == 0 ? 0 : 1;
281 // Perform half-even rounding.
282 bits += a & (b | (bits & 0x1));
283
284 // Add the exponent.
285 // No worry about overflow to the sign bit as the max exponent is 127.
286 bits += (long) (exp + EXP_BIAS) << EXP_SHIFT;
287
288 return Double.longBitsToDouble(bits);
289 }
290
291 /**
292 * Return the whole number that is nearest to the {@code double} argument {@code x}
293 * as an {@code int}, with ties rounding towards positive infinity.
294 *
295 * <p>This will raise an {@link ArithmeticException} if the closest
296 * integer result is not within the range {@code [-2^31, 2^31)},
297 * i.e. it overflows an {@code int}; or the argument {@code x}
298 * is not finite.
299 *
300 * <p>Note: This method is equivalent to:
301 * <pre>
302 * Math.toIntExact(Math.round(x))
303 * </pre>
304 *
305 * <p>The behaviour has been re-implemented for consistent error handling
306 * for {@code int}, {@code long} and {@code BigInteger} types.
307 *
308 * @param x Value.
309 * @return rounded value
310 * @throws ArithmeticException if the {@code result} overflows an {@code int},
311 * or {@code x} is not finite
312 * @see Math#round(double)
313 * @see Math#toIntExact(long)
314 */
315 static int toIntExact(double x) {
316 final double r = roundToInteger(x);
317 if (r >= -0x1.0p31 && r < 0x1.0p31) {
318 return (int) r;
319 }
320 throw new ArithmeticException("integer overflow: " + x);
321 }
322
323 /**
324 * Return the whole number that is nearest to the {@code double} argument {@code x}
325 * as an {@code long}, with ties rounding towards positive infinity.
326 *
327 * <p>This will raise an {@link ArithmeticException} if the closest
328 * integer result is not within the range {@code [-2^63, 2^63)},
329 * i.e. it overflows a {@code long}; or the argument {@code x}
330 * is not finite.
331 *
332 * @param x Value.
333 * @return rounded value
334 * @throws ArithmeticException if the {@code result} overflows a {@code long},
335 * or {@code x} is not finite
336 */
337 static long toLongExact(double x) {
338 final double r = roundToInteger(x);
339 if (r >= -0x1.0p63 && r < 0x1.0p63) {
340 return (long) r;
341 }
342 throw new ArithmeticException("long integer overflow: " + x);
343 }
344
345 /**
346 * Return the whole number that is nearest to the {@code double} argument {@code x}
347 * as an {@code int}, with ties rounding towards positive infinity.
348 *
349 * <p>This will raise an {@link ArithmeticException} if the argument {@code x}
350 * is not finite.
351 *
352 * @param x Value.
353 * @return rounded value
354 * @throws ArithmeticException if {@code x} is not finite
355 */
356 static BigInteger toBigIntegerExact(double x) {
357 if (!Double.isFinite(x)) {
358 throw new ArithmeticException("BigInteger overflow: " + x);
359 }
360 final double r = roundToInteger(x);
361 if (r >= -0x1.0p63 && r < 0x1.0p63) {
362 // Representable as a long
363 return BigInteger.valueOf((long) r);
364 }
365 // Large result
366 return new BigDecimal(r).toBigInteger();
367 }
368
369 /**
370 * Get the whole number that is the nearest to x, with ties rounding towards positive infinity.
371 *
372 * <p>This method is intended to perform the equivalent of
373 * {@link Math#round(double)} without converting to a {@code long} primitive type.
374 * This allows the domain of the result to be checked against the range {@code [-2^63, 2^63)}.
375 *
376 * <p>Note: Adapted from {@code o.a.c.math4.AccurateMath.rint} and
377 * modified to perform rounding towards positive infinity.
378 *
379 * @param x Number from which nearest whole number is requested.
380 * @return a double number r such that r is an integer {@code r - 0.5 <= x < r + 0.5}
381 */
382 private static double roundToInteger(double x) {
383 final double y = Math.floor(x);
384 final double d = x - y;
385 if (d >= HALF) {
386 // Here we do not preserve the sign of the operand in the case
387 // of -0.5 < x <= -0.0 since the rounded result is required as an integer.
388 // if y == -1.0:
389 // return -0.0
390 return y + 1.0;
391 }
392 return y;
393 }
394
395 /**
396 * Divide value {@code x} by the count {@code n}.
397 *
398 * @param x Value.
399 * @param n Count.
400 * @return the quotient
401 */
402 static double divide(Int128 x, long n) {
403 final DD a = x.toDD();
404 if (n < TWO_POW_53) {
405 // n is a representable double
406 return a.divide(n).doubleValue();
407 }
408 // Extended precision divide when n > 2^53
409 return a.divide(DD.of(n)).doubleValue();
410 }
411 }