11d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert/*
21d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Copyright (C) 2011 The Guava Authors
31d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert *
41d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Licensed under the Apache License, Version 2.0 (the "License");
51d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * you may not use this file except in compliance with the License.
61d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * You may obtain a copy of the License at
71d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert *
81d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * http://www.apache.org/licenses/LICENSE-2.0
91d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert *
101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * Unless required by applicable law or agreed to in writing, software
111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * distributed under the License is distributed on an "AS IS" BASIS,
121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
131d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * See the License for the specific language governing permissions and
141d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * limitations under the License.
151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */
161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertpackage com.google.common.math;
181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport static com.google.common.base.Preconditions.checkArgument;
207dd252788645e940eada959bdde927426e2531c9Paul Duffinimport static java.lang.Double.POSITIVE_INFINITY;
217dd252788645e940eada959bdde927426e2531c9Paul Duffinimport static java.lang.Double.doubleToRawLongBits;
227dd252788645e940eada959bdde927426e2531c9Paul Duffinimport static java.lang.Double.isNaN;
237dd252788645e940eada959bdde927426e2531c9Paul Duffinimport static java.lang.Double.longBitsToDouble;
241d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
251d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertimport java.math.BigInteger;
261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert/**
287dd252788645e940eada959bdde927426e2531c9Paul Duffin * Utilities for {@code double} primitives.
291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert *
301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert * @author Louis Wasserman
311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert */
321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringertfinal class DoubleUtils {
330888a09821a98ac0680fad765217302858e70fa4Paul Duffin  private DoubleUtils() {
340888a09821a98ac0680fad765217302858e70fa4Paul Duffin  }
35dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin
367dd252788645e940eada959bdde927426e2531c9Paul Duffin  static double nextDown(double d) {
377dd252788645e940eada959bdde927426e2531c9Paul Duffin    return -nextUp(-d);
381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
390888a09821a98ac0680fad765217302858e70fa4Paul Duffin
407dd252788645e940eada959bdde927426e2531c9Paul Duffin  static double nextUp(double d) {
417dd252788645e940eada959bdde927426e2531c9Paul Duffin    if (Double.isNaN(d)) {
427dd252788645e940eada959bdde927426e2531c9Paul Duffin      return d;
437dd252788645e940eada959bdde927426e2531c9Paul Duffin    } else if (d == 0.0) {
447dd252788645e940eada959bdde927426e2531c9Paul Duffin      return Double.MIN_VALUE;
457dd252788645e940eada959bdde927426e2531c9Paul Duffin    } else if (d == Double.POSITIVE_INFINITY) {
467dd252788645e940eada959bdde927426e2531c9Paul Duffin      return d;
47dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin    } else {
487dd252788645e940eada959bdde927426e2531c9Paul Duffin      long bits = Double.doubleToRawLongBits(d);
497dd252788645e940eada959bdde927426e2531c9Paul Duffin      bits += (bits >> 63) | 1;
507dd252788645e940eada959bdde927426e2531c9Paul Duffin      return Double.longBitsToDouble(bits);
51dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin    }
521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  // The mask for the significand, according to the {@link
551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  // Double#doubleToRawLongBits(double)} spec.
561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static final long SIGNIFICAND_MASK = 0x000fffffffffffffL;
571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
581d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  // The mask for the exponent, according to the {@link
591d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  // Double#doubleToRawLongBits(double)} spec.
601d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static final long EXPONENT_MASK = 0x7ff0000000000000L;
611d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
621d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  // The mask for the sign, according to the {@link
631d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  // Double#doubleToRawLongBits(double)} spec.
641d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static final long SIGN_MASK = 0x8000000000000000L;
651d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
661d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static final int SIGNIFICAND_BITS = 52;
671d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
681d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static final int EXPONENT_BIAS = 1023;
690888a09821a98ac0680fad765217302858e70fa4Paul Duffin
707dd252788645e940eada959bdde927426e2531c9Paul Duffin  static final int MAX_EXPONENT = 1023;
710888a09821a98ac0680fad765217302858e70fa4Paul Duffin
727dd252788645e940eada959bdde927426e2531c9Paul Duffin  static final int MIN_EXPONENT = -1022;
73dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin
741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  /**
751d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert   * The implicit 1 bit that is omitted in significands of normal doubles.
761d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert   */
771d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static final long IMPLICIT_BIT = SIGNIFICAND_MASK + 1;
780888a09821a98ac0680fad765217302858e70fa4Paul Duffin
79dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin  static int getExponent(double d) {
807dd252788645e940eada959bdde927426e2531c9Paul Duffin    long bits = Double.doubleToRawLongBits(d) & EXPONENT_MASK;
817dd252788645e940eada959bdde927426e2531c9Paul Duffin    return (int) (bits >>> SIGNIFICAND_BITS) - EXPONENT_BIAS;
82dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin  }
83dbd967a6e5c96cc1a97c5521f88dc1564ba2f81bPaul Duffin
841d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static long getSignificand(double d) {
851d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    checkArgument(isFinite(d), "not a normal value");
861d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    int exponent = getExponent(d);
877dd252788645e940eada959bdde927426e2531c9Paul Duffin    long bits = doubleToRawLongBits(d);
881d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    bits &= SIGNIFICAND_MASK;
890888a09821a98ac0680fad765217302858e70fa4Paul Duffin    return (exponent == MIN_EXPONENT - 1)
900888a09821a98ac0680fad765217302858e70fa4Paul Duffin        ? bits << 1
910888a09821a98ac0680fad765217302858e70fa4Paul Duffin        : bits | IMPLICIT_BIT;
927dd252788645e940eada959bdde927426e2531c9Paul Duffin  }
930888a09821a98ac0680fad765217302858e70fa4Paul Duffin
947dd252788645e940eada959bdde927426e2531c9Paul Duffin  static double copySign(double mag, double sgn) {
957dd252788645e940eada959bdde927426e2531c9Paul Duffin    long bits = Double.doubleToRawLongBits(mag) & ~SIGN_MASK;
967dd252788645e940eada959bdde927426e2531c9Paul Duffin    bits |= Double.doubleToRawLongBits(sgn) & SIGN_MASK;
977dd252788645e940eada959bdde927426e2531c9Paul Duffin    return Double.longBitsToDouble(bits);
981d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
991d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1001d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static boolean isFinite(double d) {
1017dd252788645e940eada959bdde927426e2531c9Paul Duffin    return getExponent(d) <= MAX_EXPONENT;
1021d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
1031d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1041d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static boolean isNormal(double d) {
1057dd252788645e940eada959bdde927426e2531c9Paul Duffin    return getExponent(d) >= MIN_EXPONENT;
1061d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
1071d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1081d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  /*
1091d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert   * Returns x scaled by a power of 2 such that it is in the range [1, 2). Assumes x is positive,
1101d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert   * normal, and finite.
1111d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert   */
1121d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static double scaleNormalize(double x) {
1137dd252788645e940eada959bdde927426e2531c9Paul Duffin    long significand = doubleToRawLongBits(x) & SIGNIFICAND_MASK;
1147dd252788645e940eada959bdde927426e2531c9Paul Duffin    return longBitsToDouble(significand | ONE_BITS);
1151d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
1161d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1171d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  static double bigToDouble(BigInteger x) {
1181d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    // This is an extremely fast implementation of BigInteger.doubleValue().  JDK patch pending.
1191d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    BigInteger absX = x.abs();
1201d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    int exponent = absX.bitLength() - 1;
1211d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    // exponent == floor(log2(abs(x)))
1221d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    if (exponent < Long.SIZE - 1) {
1231d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert      return x.longValue();
1247dd252788645e940eada959bdde927426e2531c9Paul Duffin    } else if (exponent > MAX_EXPONENT) {
1257dd252788645e940eada959bdde927426e2531c9Paul Duffin      return x.signum() * POSITIVE_INFINITY;
1261d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    }
1271d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1281d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    /*
1291d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * We need the top SIGNIFICAND_BITS + 1 bits, including the "implicit" one bit. To make
1301d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * rounding easier, we pick out the top SIGNIFICAND_BITS + 2 bits, so we have one to help us
1311d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * round up or down. twiceSignifFloor will contain the top SIGNIFICAND_BITS + 2 bits, and
1321d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * signifFloor the top SIGNIFICAND_BITS + 1.
1331d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     *
1341d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * It helps to consider the real number signif = absX * 2^(SIGNIFICAND_BITS - exponent).
1351d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     */
1361d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    int shift = exponent - SIGNIFICAND_BITS - 1;
1371d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    long twiceSignifFloor = absX.shiftRight(shift).longValue();
1381d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    long signifFloor = twiceSignifFloor >> 1;
1391d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    signifFloor &= SIGNIFICAND_MASK; // remove the implied bit
1401d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1411d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    /*
1421d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * We round up if either the fractional part of signif is strictly greater than 0.5 (which is
1431d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * true if the 0.5 bit is set and any lower bit is set), or if the fractional part of signif is
1441d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * >= 0.5 and signifFloor is odd (which is true if both the 0.5 bit and the 1 bit are set).
1451d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     */
1461d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    boolean increment = (twiceSignifFloor & 1) != 0
1471d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert        && ((signifFloor & 1) != 0 || absX.getLowestSetBit() < shift);
1481d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    long signifRounded = increment ? signifFloor + 1 : signifFloor;
1491d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    long bits = (long) ((exponent + EXPONENT_BIAS)) << SIGNIFICAND_BITS;
1501d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    bits += signifRounded;
1511d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    /*
1521d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * If signifRounded == 2^53, we'd need to set all of the significand bits to zero and add 1 to
1531d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * the exponent. This is exactly the behavior we get from just adding signifRounded to bits
1541d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * directly.  If the exponent is MAX_DOUBLE_EXPONENT, we round up (correctly) to
1551d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     * Double.POSITIVE_INFINITY.
1561d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert     */
1571d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert    bits |= x.signum() & SIGN_MASK;
1587dd252788645e940eada959bdde927426e2531c9Paul Duffin    return longBitsToDouble(bits);
1597dd252788645e940eada959bdde927426e2531c9Paul Duffin  }
1607dd252788645e940eada959bdde927426e2531c9Paul Duffin
1617dd252788645e940eada959bdde927426e2531c9Paul Duffin  /**
1627dd252788645e940eada959bdde927426e2531c9Paul Duffin   * Returns its argument if it is non-negative, zero if it is negative.
1637dd252788645e940eada959bdde927426e2531c9Paul Duffin   */
1647dd252788645e940eada959bdde927426e2531c9Paul Duffin  static double ensureNonNegative(double value) {
1657dd252788645e940eada959bdde927426e2531c9Paul Duffin    checkArgument(!isNaN(value));
1667dd252788645e940eada959bdde927426e2531c9Paul Duffin    if (value > 0.0) {
1677dd252788645e940eada959bdde927426e2531c9Paul Duffin      return value;
1687dd252788645e940eada959bdde927426e2531c9Paul Duffin    } else {
1697dd252788645e940eada959bdde927426e2531c9Paul Duffin      return 0.0;
1707dd252788645e940eada959bdde927426e2531c9Paul Duffin    }
1711d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert  }
1721d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert
1737dd252788645e940eada959bdde927426e2531c9Paul Duffin  private static final long ONE_BITS = doubleToRawLongBits(1.0);
1741d580d0f6ee4f21eb309ba7b509d2c6d671c4044Bjorn Bringert}
175