*** 23,32 **** --- 23,33 ---- * questions. */ package java.lang; + import java.math.BigDecimal; import java.util.Random; import jdk.internal.math.FloatConsts; import jdk.internal.math.DoubleConsts; import jdk.internal.HotSpotIntrinsicCandidate;

*** 1448,1457 **** --- 1449,1621 ---- } return (a <= b) ? a : b; } /** + * Returns the fused multiply-accumulate of the three arguments; + * that is, returns the exact product of the first two arguments + * summed with the third argument and then rounded once to the + * nearest {@code double}. + * + * The rounding is done using the {@linkplain + * java.math.RoundingMode#HALF_EVEN round to nearest even + * rounding mode}. + * + * In contrast, if {@code a * b + c} is evaluated as a regular + * floating-point expression, two rounding errors are involved, + * the first for the multiply operation, the second for the + * addition operation. + * + * <p>Special cases: + * <ul> + * <li> If any argument is NaN, the result is NaN. + * + * <li> If one of the first two arguments is infinite and the + * other is zero, the result is NaN. + * + * <li> If the exact product of the first two arguments is infinite + * (in other words, at least one of the arguments is infinite and + * the other is neither zero nor NaN) and the third argument is an + * infinity of the opposite sign, the result is NaN. + * + * </ul> + * + * <p>Note that {@code fusedMac(a, 1.0, c)} returns the same + * result as ({@code a + c}). However, + * {@code fusedMac(a, b, +0.0)} does <em>not</em> always return the + * same result as ({@code a * b}) since + * {@code fusedMac(-0.0, +0.0, +0.0)} is {@code +0.0} while + * ({@code 0.0 * +0.0}) is {@code -0.0}; {@code fusedMac(a, b, -0.0)} is + * equivalent to ({@code a * b}) however. + * + * @param a a value + * @param b a value + * @param c a value + * + * @return (<i>a</i> × <i>b</i> + <i>c</i>) + * computed, as if with unlimited range and precision, and rounded + * once to the nearest {@code double} value + */ + // @HotSpotIntrinsicCandidate + public static double fusedMac(double a, double b, double c) { + if (!Double.isFinite(a) || + !Double.isFinite(b) || + !Double.isFinite(c)) { + // Infinite and NaN arithmetic is not quite the same with + // two roundings as opposed to just one. With two + // roundings, the product and overflow and if the addend + // is infinite, a spurious NaN can be produced. + + if (Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(c)) { + return Double.NaN; + } else { // One or more infinities in the input + if ((Double.isInfinite(a) && b == 0.0) || + (Double.isInfinite(b) && a == 0.0)) { + return Double.NaN; + } else { + double product = a * b; + double result = product + c; + if (Double.isInfinite(product) && + (!Double.isInfinite(a) && !Double.isInfinite(b))) { + // Intermediate overflow; might cause a + // spurious NaN with added to infinite c. + return c; + } + return result; + } + } + } if (a == 1.0 || // Could just return b + c + b == 1.0 || // Could just return a + c + a == 0.0 || // Handle any signed zero issues + b == 0.0 || // Handle any signed zero issues + c == 0.0) { // Could be +0.0 or -0.0) + // The special cases of a == 1.0, b == 1.0, and c == 0.0 + // are handled as double operations rather than falling + // back on BigDecimal. The case where a == 1.0 or b == 1.0 + // would be handled correctly by the BigDecimal code + // below, unlike the signed zero case where the exact + // product is -0.0 and c is -0.0. + + // As noted above, if c is zero, it must still be added in + // since a negative zero can change the sign of the final + // result if the product is zero. At the cost of + // additional checks on a and b, the addition of a zero c + // could be avoided. + + return a * b + c; + } else { + // The set of values representable in BigDecimal includes + // neither signed zero nor infinities nor NaN. The logic + // above handles those cases. All finite double values are + // representable in BigDecimal. The code below exactly + // converts a, b, and c to their equivalent BigDecimal + // numbers, exactly computes the product and sum as + // appropriate, and then performs a single rounding to + // double using the round to nearest rounding mode. + return (((new BigDecimal(a)).multiply(new BigDecimal(b))). + add(new BigDecimal(c))).doubleValue(); + } + } + + /** + * Returns the fused multiply-accumulate of the three arguments; + * that is, returns the exact product of the first two arguments + * summed with the third argument and then rounded once to the + * nearest {@code float}. + * + * The rounding is done using the {@linkplain + * java.math.RoundingMode#HALF_EVEN round to nearest even + * rounding mode}. + * + * In contrast, if {@code a * b + c} is evaluated as a regular + * floating-point expression, two rounding errors are involved, + * the first for the multiply operation, the second for the + * addition operation. + * + * <p>Special cases: + * <ul> + * <li> If any argument is NaN, the result is NaN. + * + * <li> If one of the first two arguments is infinite and the + * other is zero, the result is NaN. + * + * <li> If the exact product of the first two arguments is infinite + * (in other words, at least one of the arguments is infinite and + * the other is neither zero nor NaN) and the third argument is an + * infinity of the opposite sign, the result is NaN. + * + * </ul> + * + * <p>Note that {@code fusedMac(a, 1.0f, c)} returns the same + * result as ({@code a + c}). However, + * {@code fusedMac(a, b, +0.0f)} does <em>not</em> always return the + * same result as ({@code a * b}) since + * {@code fusedMac(-0.0f, +0.0f, +0.0f)} is {@code +0.0f} while + * ({@code 0.0f * +0.0f}) is {@code -0.0f}; {@code fusedMac(a, b, -0.0f)} is + * equivalent to ({@code a * b}) however. + * + * @param a a value + * @param b a value + * @param c a value + * + * @return (<i>a</i> × <i>b</i> + <i>c</i>) + * computed, as if with unlimited range and precision, and rounded + * once to the nearest {@code float} value + */ + // @HotSpotIntrinsicCandidate + public static float fusedMac(float a, float b, float c) { + // Since double has more than twice the precision of float, + // the multiply of a * b is exact in double. The add of c to + // the product then has one rounded error. Since double + // moreover has more than 2p + 2 precision compared to float, + // the double rounding of (a*b + c) from the double format to + // the float format is equivalent to rounding the result + // directly to float precision. + return (float)(((double) a * (double) b ) + (double) c); + } + + /** * Returns the size of an ulp of the argument. An ulp, unit in * the last place, of a {@code double} value is the positive * distance between this floating-point value and the {@code * double} value next larger in magnitude. Note that for non-NaN * <i>x</i>, <code>ulp(-<i>x</i>) == ulp(<i>x</i>)</code>.