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