Merge branch 'NUMBERS-120__heinrich'

Closes #63.

commons-numbers-fraction/src/main/java/org/apache/commons/numbers/fraction/BigFraction.java

@@ -116,6 +116,7 @@ private BigFraction(BigInteger num, BigInteger den) {
      * in absolute terms.
      * @param maxDenominator Maximum denominator value allowed.
      * @param maxIterations Maximum number of convergents.
+     * @return a new instance.
      * @throws ArithmeticException if the continued fraction failed to converge.
      */
     private static BigFraction from(final double value,
@@ -660,6 +661,167 @@ public BigFraction divide(final BigFraction fraction) {
         return multiply(fraction.reciprocal());
     }
 
+    /**
+     * Calculates the sign bit, the biased exponent and the significand for a
+     * binary floating-point representation of this {@code BigFraction}
+     * according to the IEEE 754 standard, and encodes these values into a {@code long}
+     * variable. The representative bits are arranged adjacent to each other and
+     * placed at the low-order end of the returned {@code long} value, with the
+     * least significant bits used for the significand, the next more
+     * significant bits for the exponent, and next more significant bit for the
+     * sign.
+     * @param exponentLength the number of bits allowed for the exponent; must be
+     *                       between 1 and 32 (inclusive), and must not be greater
+     *                       than {@code 63 - significandLength}
+     * @param significandLength the number of bits allowed for the significand
+     *                          (excluding the implicit leading 1-bit in
+     *                          normalized numbers, e.g. 52 for a double-precision
+     *                          floating-point number); must be between 1 and
+     *                          {@code 63 - exponentLength} (inclusive)
+     * @return the bits of an IEEE 754 binary floating-point representation of
+     * this fraction encoded in a {@code long}, as described above.
+     * @throws IllegalArgumentException if the arguments do not meet the
+     *         criteria described above
+     */
+    private long toFloatingPointBits(int exponentLength, int significandLength) {
+        if (exponentLength < 1 ||
+            significandLength < 1 ||
+            exponentLength > Math.min(32, 63 - significandLength)) {
+            throw new IllegalArgumentException("exponent length: " + exponentLength +
+                                               "; significand length: " + significandLength);
+        }
+        if (numerator.signum() == 0) {
+            return 0L;
+        }
+
+        final long sign = numerator.signum() == -1 ? 1L : 0L;
+        final BigInteger positiveNumerator = numerator.abs();
+
+        /*
+         * The most significant 1-bit of a non-zero number is not explicitly
+         * stored in the significand of an IEEE 754 normalized binary
+         * floating-point number, so we need to round the value of this fraction
+         * to (significandLength + 1) bits. In order to do this, we calculate
+         * the most significant (significandLength + 2) bits, and then, based on
+         * the least significant of those bits, find out whether we need to
+         * round up or down.
+         *
+         * First, we'll remove all powers of 2 from the denominator because they
+         * are not relevant for the significand of the prospective binary
+         * floating-point value.
+         */
+        final int denRightShift = denominator.getLowestSetBit();
+        final BigInteger divisor = denominator.shiftRight(denRightShift);
+
+        /*
+         * Now, we're going to calculate the (significandLength + 2) most
+         * significant bits of the fraction's value using integer division. To
+         * guarantee that the quotient of the division has at least
+         * (significandLength + 2) bits, the bit length of the dividend must
+         * exceed that of the divisor by at least that amount.
+         *
+         * If the denominator has prime factors other than 2, i.e. if the
+         * divisor was not reduced to 1, an excess of exactly
+         * (significandLength + 2) bits is sufficient, because the knowledge
+         * that the fractional part of the precise quotient's binary
+         * representation does not terminate is enough information to resolve
+         * cases where the most significant (significandLength + 2) bits alone
+         * are not conclusive.
+         *
+         * Otherwise, the quotient must be calculated exactly and the bit length
+         * of the numerator can only be reduced as long as no precision is lost
+         * in the process (meaning it can have powers of 2 removed, like the
+         * denominator).
+         */
+        int numRightShift = positiveNumerator.bitLength() - divisor.bitLength() - (significandLength + 2);
+        if (numRightShift > 0 &&
+            divisor.equals(BigInteger.ONE)) {
+            numRightShift = Math.min(numRightShift, positiveNumerator.getLowestSetBit());
+        }
+        final BigInteger dividend = positiveNumerator.shiftRight(numRightShift);
+
+        final BigInteger quotient = dividend.divide(divisor);
+
+        int quotRightShift = quotient.bitLength() - (significandLength + 1);
+        long significand = roundAndRightShift(
+                quotient,
+                quotRightShift,
+                !divisor.equals(BigInteger.ONE)
+        ).longValue();
+
+        /*
+         * If the significand had to be rounded up, this could have caused the
+         * bit length of the significand to increase by one.
+         */
+        if ((significand & (1L << (significandLength + 1))) != 0) {
+            significand >>= 1;
+            quotRightShift++;
+        }
+
+        /*
+         * Now comes the exponent. The absolute value of this fraction based on
+         * the current local variables is:
+         *
+         * significand * 2^(numRightShift - denRightShift + quotRightShift)
+         *
+         * To get the unbiased exponent for the floating-point value, we need to
+         * add (significandLength) to the above exponent, because all but the
+         * most significant bit of the significand will be treated as a
+         * fractional part. To convert the unbiased exponent to a biased
+         * exponent, we also need to add the exponent bias.
+         */
+        final int exponentBias = (1 << (exponentLength - 1)) - 1;
+        long exponent = numRightShift - denRightShift + quotRightShift + significandLength + exponentBias;
+        final long maxExponent = (1L << exponentLength) - 1L; //special exponent for infinities and NaN
+
+        if (exponent >= maxExponent) { //infinity
+            exponent = maxExponent;
+            significand = 0L;
+        } else if (exponent > 0) { //normalized number
+            significand &= -1L >>> (64 - significandLength); //remove implicit leading 1-bit
+        } else { //smaller than the smallest normalized number
+            /*
+             * We need to round the quotient to fewer than
+             * (significandLength + 1) bits. This must be done with the original
+             * quotient and not with the current significand, because the loss
+             * of precision in the previous rounding might cause a rounding of
+             * the current significand's value to produce a different result
+             * than a rounding of the original quotient.
+             *
+             * So we find out how many high-order bits from the quotient we can
+             * transfer into the significand. The absolute value of the fraction
+             * is:
+             *
+             * quotient * 2^(numRightShift - denRightShift)
+             *
+             * To get the significand, we need to right shift the quotient so
+             * that the above exponent becomes (1 - exponentBias - significandLength)
+             * (the unbiased exponent of a subnormal floating-point number is
+             * defined as equivalent to the minimum unbiased exponent of a
+             * normalized floating-point number, and (- significandLength)
+             * because the significand will be treated as the fractional part).
+             */
+            significand = roundAndRightShift(quotient,
+                                             (1 - exponentBias - significandLength) - (numRightShift - denRightShift),
+                                             !divisor.equals(BigInteger.ONE)).longValue();
+            exponent = 0L;
+
+            /*
+             * Note: It is possible that an otherwise subnormal number will
+             * round up to the smallest normal number. However, this special
+             * case does not need to be treated separately, because the
+             * overflowing highest-order bit of the significand will then simply
+             * become the lowest-order bit of the exponent, increasing the
+             * exponent from 0 to 1 and thus establishing the implicity of the
+             * leading 1-bit.
+             */
+        }
+
+        return (sign << (significandLength + exponentLength)) |
+            (exponent << significandLength) |
+            significand;
+    }
+
     /**
      * <p>
      * Gets the fraction as a {@code double}. This calculates the fraction as
@@ -671,19 +833,37 @@ public BigFraction divide(final BigFraction fraction) {
      */
     @Override
     public double doubleValue() {
-        double doubleNum = numerator.doubleValue();
-        double doubleDen = denominator.doubleValue();
-        double result = doubleNum / doubleDen;
-        if (Double.isInfinite(doubleNum) ||
-            Double.isInfinite(doubleDen) ||
-            Double.isNaN(result)) {
-            // Numerator and/or denominator must be out of range:
-            // Calculate how far to shift them to put them in range.
-            int shift = Math.max(numerator.bitLength(),
-                                 denominator.bitLength()) - Math.getExponent(Double.MAX_VALUE);
-            result = numerator.shiftRight(shift).doubleValue() /
-                denominator.shiftRight(shift).doubleValue();
+        return Double.longBitsToDouble(toFloatingPointBits(11, 52));
+    }
+
+    /**
+     * Rounds an integer to the specified power of two (i.e. the minimum number of
+     * low-order bits that must be zero) and performs a right-shift by this
+     * amount. The rounding mode applied is round to nearest, with ties rounding
+     * to even (meaning the prospective least significant bit must be zero). The
+     * number can optionally be treated as though it contained at
+     * least one 0-bit and one 1-bit in its fractional part, to influence the result in cases
+     * that would otherwise be a tie.
+     * @param value the number to round and right-shift
+     * @param bits the power of two to which to round; must be positive
+     * @param hasFractionalBits whether the number should be treated as though
+     *                          it contained a non-zero fractional part
+     * @return a {@code BigInteger} as described above
+     * @throws IllegalArgumentException if {@code bits <= 0}
+     */
+    private static BigInteger roundAndRightShift(BigInteger value, int bits, boolean hasFractionalBits) {
+        if (bits <= 0) {
+            throw new IllegalArgumentException("bits: " + bits);
         }
+
+        BigInteger result = value.shiftRight(bits);
+        if (value.testBit(bits - 1) &&
+            (hasFractionalBits ||
+             (value.getLowestSetBit() < bits - 1) ||
+             value.testBit(bits))) {
+            result = result.add(BigInteger.ONE); //round up
+        }
+
         return result;
     }
 
@@ -727,20 +907,7 @@ public boolean equals(final Object other) {
      */
     @Override
     public float floatValue() {
-        float floatNum = numerator.floatValue();
-        float floatDen = denominator.floatValue();
-        float result = floatNum / floatDen;
-        if (Float.isInfinite(floatNum) ||
-            Float.isInfinite(floatDen) ||
-            Float.isNaN(result)) {
-            // Numerator and/or denominator must be out of range:
-            // Calculate how far to shift them to put them in range.
-            int shift = Math.max(numerator.bitLength(),
-                                 denominator.bitLength()) - Math.getExponent(Float.MAX_VALUE);
-            result = numerator.shiftRight(shift).floatValue() /
-                denominator.shiftRight(shift).floatValue();
-        }
-        return result;
+        return Float.intBitsToFloat((int) toFloatingPointBits(8, 23));
     }
 
     /**