001/** 002 * The MIT License (MIT) 003 * 004 * Copyright (c) 2015-2016 decimal4j (tools4j), Marco Terzer 005 * 006 * Permission is hereby granted, free of charge, to any person obtaining a copy 007 * of this software and associated documentation files (the "Software"), to deal 008 * in the Software without restriction, including without limitation the rights 009 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 010 * copies of the Software, and to permit persons to whom the Software is 011 * furnished to do so, subject to the following conditions: 012 * 013 * The above copyright notice and this permission notice shall be included in all 014 * copies or substantial portions of the Software. 015 * 016 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 017 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 018 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 019 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 020 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 021 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 022 * SOFTWARE. 023 */ 024package org.decimal4j.arithmetic; 025 026import org.decimal4j.api.DecimalArithmetic; 027import org.decimal4j.scale.ScaleMetrics; 028import org.decimal4j.truncate.DecimalRounding; 029import org.decimal4j.truncate.TruncatedPart; 030 031/** 032 * Contains methods to convert from and to double. 033 */ 034final class DoubleConversion { 035 036 private static final long LONG_MASK = 0xffffffffL; 037 038 // The mask for the significand, according to the {@link Double#doubleToRawLongBits(double)} spec. 039 private static final long SIGNIFICAND_MASK = 0x000fffffffffffffL; 040 041 // The mask for the exponent, according to the {@link Double#doubleToRawLongBits(double)} spec. 042 @SuppressWarnings("unused") 043 private static final long EXPONENT_MASK = 0x7ff0000000000000L; 044 045 // The mask for the sign, according to the {@link Double#doubleToRawLongBits(double)} spec. 046 private static final long SIGN_MASK = 0x8000000000000000L; 047 048 private static final int SIGNIFICAND_BITS = 52; 049 050 private static final int EXPONENT_BIAS = 1023; 051 052 /** 053 * The implicit 1 bit that is omitted in significands of normal doubles. 054 */ 055 private static final long IMPLICIT_BIT = SIGNIFICAND_MASK + 1; 056 057 private static final double MIN_LONG_AS_DOUBLE = -0x1p63; 058 /* 059 * We cannot store Long.MAX_VALUE as a double without losing precision. Instead, we store Long.MAX_VALUE + 1 == 060 * -Long.MIN_VALUE, and then offset all comparisons by 1. 061 */ 062 private static final double MAX_LONG_AS_DOUBLE_PLUS_ONE = 0x1p63; 063 064 /** 065 * Converts the specified double value to a long truncating the fractional part if any is present. If the value is 066 * NaN, infinite or outside of the valid long range, an exception is thrown. 067 * 068 * @param value 069 * the value to convert 070 * @return <tt>round<sub>DOWN</sub>(value)</tt> 071 * @throws IllegalArgumentException 072 * if {@code value} is NaN or infinite or if the magnitude is too large for the double to be represented 073 * as a {@code long} 074 */ 075 public static final long doubleToLong(double value) { 076 if (Double.isNaN(value)) { 077 throw new IllegalArgumentException("Cannot convert double to long: " + value); 078 } 079 if (isInLongRange(value)) { 080 return (long) value; 081 } 082 throw new IllegalArgumentException("Overflow for conversion from double to long: " + value); 083 } 084 085 /** 086 * Converts the specified double value to a long rounding the fractional part if necessary using the given 087 * {@code rounding} mode. If the value is NaN, infinite or outside of the valid long range, an exception is thrown. 088 * 089 * @param rounding 090 * the rounding to apply if necessary 091 * @param value 092 * the value to convert 093 * @return <tt>round(value)</tt> 094 * @throws IllegalArgumentException 095 * if {@code value} is NaN or infinite or if the magnitude is too large for the double to be represented 096 * as a {@code long} 097 * @throws ArithmeticException 098 * if {@code roundingMode==UNNECESSARY} and rounding is necessary 099 */ 100 public static final long doubleToLong(DecimalRounding rounding, double value) { 101 if (Double.isNaN(value)) { 102 throw new IllegalArgumentException("Cannot convert double to long: " + value); 103 } 104 if (isInLongRange(value)) { 105 return (long) roundIntermediate(value, rounding); 106 } 107 throw new IllegalArgumentException("Overflow for conversion from double to long: " + value); 108 } 109 110 /* 111 * Copied from guava. This method returns a value y such that rounding y DOWN (towards zero) gives the same result 112 * as rounding x according to the specified mode. PRECONDITION: isFinite(x) 113 */ 114 private static final double roundIntermediate(double x, DecimalRounding mode) { 115 switch (mode) { 116 case UNNECESSARY: 117 if (!isMathematicalInteger(x)) { 118 throw new ArithmeticException("Rounding necessary to convert to an integer value: " + x); 119 } 120 return x; 121 case FLOOR: 122 if (x >= 0.0 || isMathematicalInteger(x)) { 123 return x; 124 } else { 125 return (long)x - 1L; 126 } 127 case CEILING: 128 if (x <= 0.0 || isMathematicalInteger(x)) { 129 return x; 130 } else { 131 return (long)x + 1L; 132 } 133 case DOWN: 134 return x; 135 case UP: 136 if (isMathematicalInteger(x)) { 137 return x; 138 } else { 139 return (long)x + (x > 0 ? 1L : -1L); 140 } 141 case HALF_EVEN: 142 return Math.rint(x); 143 case HALF_UP: { 144 double z = Math.rint(x); 145 if (Math.abs(x - z) == 0.5) { 146 return x + Math.copySign(0.5, x); 147 } else { 148 return z; 149 } 150 } 151 case HALF_DOWN: { 152 double z = Math.rint(x); 153 if (Math.abs(x - z) == 0.5) { 154 return x; 155 } else { 156 return z; 157 } 158 } 159 default: 160 throw new IllegalArgumentException("Unsupported rounding mode: " + mode); 161 } 162 } 163 164 /** 165 * Converts the specified double value to an unscaled decimal truncating extra fractional digits if necessary. If 166 * the value is NaN, infinite or outside of the valid Decimal range, an exception is thrown. 167 * 168 * @param arith 169 * the arithmetic associated with the result value 170 * @param value 171 * the value to convert 172 * @return <tt>round(value)</tt> 173 * @throws IllegalArgumentException 174 * if {@code value} is NaN or infinite or if the magnitude is too large for the double to be represented 175 * as a Decimal of the arithmetic's scale 176 */ 177 public static final long doubleToUnscaled(DecimalArithmetic arith, double value) { 178 return doubleToUnscaled(arith, DecimalRounding.DOWN, value); 179 } 180 181 /** 182 * Converts the specified double value to an unscaled decimal. The specified {@code rounding} mode is used if 183 * rounding is necessary. If the value is NaN, infinite or outside of the valid Decimal range, an exception is 184 * thrown. 185 * 186 * @param arith 187 * the arithmetic associated with the result value 188 * @param rounding 189 * the rounding to apply if necessary 190 * @param value 191 * the value to convert 192 * @return <tt>round(value)</tt> 193 * @throws IllegalArgumentException 194 * if {@code value} is NaN or infinite or if the magnitude is too large for the double to be represented 195 * as a Decimal of the arithmetic's scale 196 * @throws ArithmeticException 197 * if {@code roundingMode==UNNECESSARY} and rounding is necessary 198 */ 199 public static final long doubleToUnscaled(DecimalArithmetic arith, DecimalRounding rounding, double value) { 200 if (value == 0) { 201 return 0; 202 } 203 final int exp = Math.getExponent(value); 204 if (exp >= Long.SIZE) { 205 throw newOverflowException(arith, value); 206 } 207 208 // multiply significand by scale factor into a 128bit integer 209 final ScaleMetrics scaleMetrics = arith.getScaleMetrics(); 210 final long significand = getSignificand(value); 211 212 // HD + Knuth's Algorithm M from [Knu2] section 4.3.1. 213 final int lFactor = (int) (significand & LONG_MASK); 214 final int hFactor = (int) (significand >>> 32); 215 final long w1, w2, w3; 216 long k, t; 217 218 t = scaleMetrics.mulloByScaleFactor(lFactor); 219 w3 = t & LONG_MASK; 220 k = t >>> 32; 221 222 t = scaleMetrics.mulloByScaleFactor(hFactor) + k; 223 w2 = t & LONG_MASK; 224 w1 = t >>> 32; 225 226 t = scaleMetrics.mulhiByScaleFactor(lFactor) + w2; 227 k = t >>> 32; 228 229 final long hScaled = scaleMetrics.mulhiByScaleFactor(hFactor) + w1 + k; 230 final long lScaled = (t << 32) + w3; 231 232 // now multiply or divide by powers of two as instructed by the double 233 // exponent 234 final int shift = exp - SIGNIFICAND_BITS; 235 return doubleToUnscaledShift(arith, rounding, value, hScaled, lScaled, shift); 236 } 237 238 private static final long doubleToUnscaledShift(DecimalArithmetic arith, DecimalRounding rounding, double value, long hScaled, long lScaled, int shift) { 239 if (shift > 0) { 240 // multiply: shift left 241 if (hScaled != 0) { 242 throw newOverflowException(arith, value); 243 } 244 final int zeros = Long.numberOfLeadingZeros(lScaled); 245 if (shift >= zeros) { 246 throw newOverflowException(arith, value); 247 } 248 final long absResult = lScaled << shift; 249 return value >= 0 ? absResult : -absResult; 250 } else if (shift == 0) { 251 if (hScaled != 0 | lScaled < 0) { 252 throw newOverflowException(arith, value); 253 } 254 return value >= 0 ? lScaled : -lScaled; 255 } else {// shift < 0 256 // divide: shift right 257 if (rounding == DecimalRounding.DOWN) { 258 return doubleToUnscaledShiftRight(arith, value, hScaled, lScaled, -shift); 259 } 260 return doubleToUnscaledShiftRight(arith, rounding, value, hScaled, lScaled, -shift); 261 } 262 } 263 264 private static final long doubleToUnscaledShiftRight(DecimalArithmetic arith, double value, long hScaled, long lScaled, int shift) { 265 final long absResult; 266 if (shift < Long.SIZE) { 267 if ((hScaled >>> shift) != 0) { 268 throw newOverflowException(arith, value); 269 } 270 absResult = (hScaled << (Long.SIZE - shift)) | (lScaled >>> shift); 271 } else if (shift < 2 * Long.SIZE) { 272 absResult = (hScaled >>> (shift - Long.SIZE)); 273 } else { 274 return 0;// rounded down 275 } 276 if (absResult < 0) { 277 throw newOverflowException(arith, value); 278 } 279 return value >= 0 ? absResult : -absResult; 280 } 281 282 private static final long doubleToUnscaledShiftRight(DecimalArithmetic arith, DecimalRounding rounding, double value, long hScaled, long lScaled, int shift) { 283 final long absResult; 284 final TruncatedPart truncatedPart; 285 if (shift < Long.SIZE) { 286 if ((hScaled >>> shift) != 0) { 287 throw newOverflowException(arith, value); 288 } 289 absResult = (hScaled << (Long.SIZE - shift)) | (lScaled >>> shift); 290 final long rem = modPow2(lScaled, shift); 291 truncatedPart = Rounding.truncatedPartFor2powN(rem, shift); 292 } else if (shift < 2 * Long.SIZE) { 293 absResult = (hScaled >>> (shift - Long.SIZE)); 294 final long rem = modPow2(hScaled, shift - Long.SIZE); 295 truncatedPart = Rounding.truncatedPartFor2powN(rem, lScaled, shift); 296 } else { 297 absResult = 0;// rounded down 298 truncatedPart = Rounding.truncatedPartFor2powN(hScaled, lScaled, shift); 299 } 300 final int inc = absResult < 0 ? 0 301 : rounding.calculateRoundingIncrement(value >= 0 ? 1 : -1, absResult, truncatedPart); 302 if (absResult < 0 | (value >= 0 & absResult == Long.MAX_VALUE & inc == 1)) { 303 throw newOverflowException(arith, value); 304 } 305 return (value >= 0 ? absResult : -absResult) + inc; 306 } 307 308 /** 309 * Converts the specified long value to a double truncating extra mantissa digits if necessary. 310 * 311 * @param arith 312 * the arithmetic associated with the value 313 * @param value 314 * the long value 315 * @return <tt>round<sub>DOWN</sub>(value)</tt> 316 */ 317 public static final double longToDouble(DecimalArithmetic arith, long value) { 318 return unscaledToDouble(arith, DecimalRounding.DOWN, value); 319 } 320 321 /** 322 * Converts the specified long value to a double rounding extra mantissa digits if necessary. 323 * 324 * @param arith 325 * the arithmetic associated with the value 326 * @param rounding 327 * the rounding to apply if necessary 328 * @param value 329 * the long value 330 * @return <tt>round(value)</tt> 331 * @throws ArithmeticException 332 * if {@code roundingMode==UNNECESSARY} and rounding is necessary 333 */ 334 public static final double longToDouble(DecimalArithmetic arith, DecimalRounding rounding, long value) { 335 if (rounding == DecimalRounding.HALF_EVEN) { 336 return value; 337 } 338 return unscaledToDouble(arith, rounding, value); 339 } 340 341 /** 342 * Converts the specified unscaled decimal value to a double truncating extra precision digits if necessary. 343 * 344 * @param arith 345 * the arithmetic associated with the value 346 * @param unscaled 347 * the unscaled decimal value 348 * @return <tt>round<sub>DOWN</sub>(value)</tt> 349 */ 350 public static final double unscaledToDouble(DecimalArithmetic arith, long unscaled) { 351 return unscaledToDouble(arith, DecimalRounding.DOWN, unscaled); 352 } 353 354 /** 355 * Converts the specified unscaled decimal value to a double rounding extra precision digits if necessary. 356 * 357 * @param arith 358 * the arithmetic associated with the value 359 * @param rounding 360 * the rounding to apply if necessary 361 * @param unscaled 362 * the unscaled decimal value 363 * @return <tt>round(value)</tt> 364 * @throws ArithmeticException 365 * if {@code roundingMode==UNNECESSARY} and rounding is necessary 366 */ 367 public static final double unscaledToDouble(DecimalArithmetic arith, DecimalRounding rounding, long unscaled) { 368 if (unscaled == 0) { 369 return 0; 370 } 371 final ScaleMetrics scaleMetrics = arith.getScaleMetrics(); 372 373 final long absUnscaled = Math.abs(unscaled); 374 if (absUnscaled < (1L << SIGNIFICAND_BITS) & rounding == DecimalRounding.HALF_EVEN) { 375 // Don't have too guard against Math.abs(MIN_VALUE) 376 // because MIN_VALUE also works without loss of precision 377 return (double)unscaled / scaleMetrics.getScaleFactor(); 378 } 379 380 // eliminate sign and trailing power-of-2 zero bits 381 final int pow2 = Long.numberOfTrailingZeros(absUnscaled); 382 final long absVal = absUnscaled >>> pow2; 383 final int nlzAbsVal = Long.numberOfLeadingZeros(absVal); 384 385 /* 386 * NOTE: a) If absVal has no more than 53 bits it can be represented as a double value without loss of precision 387 * (52 mantissa bits plus the implicit leading 1 bit) b) The scale factor has never more than 53 bits if shifted 388 * right by the trailing power-of-2 zero bits ==> For HALF_EVEN rounding mode we can therefore apply the scale 389 * factor via double division without losing information 390 */ 391 if (Long.SIZE - nlzAbsVal <= SIGNIFICAND_BITS + 1 & rounding == DecimalRounding.HALF_EVEN) { 392 return unscaledToDoubleWithDoubleDivisionRoundHalfEven(scaleMetrics, unscaled, pow2, absVal); 393 } 394 395 /* 396 * 1) we align absVal and factor such that: 2*factor > absVal >= factor 397 * then the division absVal/factor == 1.xxxxx, i.e. it is normalized 398 * 2) because we omit the 1 in the mantissa, we calculate 399 * valModFactor = absVal - floor(absVal/factor)*factor = absVal - 1*factor 400 * 3) we shift valModFactor such that the 1 from the division would be on bit 53 401 * 4) we perform the division 402 */ 403 404 // (1) + (2) 405 final int exp; 406 final int mantissaShift; 407 final long valModFactor; 408 final int alignShift = nlzAbsVal - scaleMetrics.getScaleFactorNumberOfLeadingZeros(); 409 if (alignShift >= 0) { 410 final long scaledAbsVal = absVal << alignShift; 411 final long diff = scaledAbsVal - scaleMetrics.getScaleFactor(); 412 exp = -alignShift + (int) (diff >> 63); 413 // if scaledAbsVal < factor we shift left by 1, i.e. we add the absVal 414 valModFactor = diff + ((diff >> 63) & scaledAbsVal); 415 mantissaShift = SIGNIFICAND_BITS; 416 } else { 417 final long scaledFactor = scaleMetrics.getScaleFactor() << -alignShift; 418 if (Unsigned.isLess(absVal, scaledFactor)) { 419 exp = -alignShift - 1; 420 // if absVal < scaledFactor we shift left by 1 (right shift of scaledFactor to avoid overflow) 421 valModFactor = absVal - (scaledFactor >>> 1); 422 mantissaShift = SIGNIFICAND_BITS + alignShift + 1; 423 } else { 424 exp = -alignShift; 425 valModFactor = absVal - scaledFactor; 426 mantissaShift = SIGNIFICAND_BITS + alignShift; 427 } 428 } 429 if (rounding == DecimalRounding.DOWN) { 430 return unscaledToDoubleShiftAndDivideByScaleFactor(scaleMetrics, unscaled, exp + pow2, mantissaShift, 431 valModFactor); 432 } 433 // (3) + (4) 434 return unscaledToDoubleShiftAndDivideByScaleFactor(scaleMetrics, rounding, unscaled, exp + pow2, mantissaShift, 435 valModFactor); 436 } 437 438 private static final double unscaledToDoubleWithDoubleDivisionRoundHalfEven(ScaleMetrics scaleMetrics, long unscaled, int pow2, long absVal) { 439 final int scale = scaleMetrics.getScale(); 440 final double dividend = absVal; 441 final double divisor = scaleMetrics.getScaleFactor() >> scale; 442 final double quotient = dividend / divisor; 443 final int exponent = Math.getExponent(quotient) + pow2 - scale; 444 final long significand = Double.doubleToRawLongBits(quotient) & SIGNIFICAND_MASK; 445 final long raw = (unscaled & SIGN_MASK) | (((long) (exponent + EXPONENT_BIAS)) << SIGNIFICAND_BITS) 446 | significand; 447 return Double.longBitsToDouble(raw); 448 } 449 450 private static final double unscaledToDoubleShiftAndDivideByScaleFactor(ScaleMetrics scaleMetrics, long unscaled, int exp, int mantissaShift, long valModFactor) { 451 final long quot; 452 if (mantissaShift >= 0) { 453 final long hValModFactor = (valModFactor >>> (Long.SIZE - mantissaShift)) & (-mantissaShift >> 63); 454 final long lValModFactor = valModFactor << mantissaShift; 455 if (hValModFactor == 0) { 456 quot = scaleMetrics.divideUnsignedByScaleFactor(lValModFactor); 457 } else { 458 quot = Math.abs(Div.div128by64(DecimalRounding.DOWN, unscaled < 0, hValModFactor, lValModFactor, 459 scaleMetrics.getScaleFactor())); 460 } 461 } else { 462 quot = scaleMetrics.divideByScaleFactor(valModFactor >>> -mantissaShift); 463 } 464 final long raw = (unscaled & SIGN_MASK) | (((long) (exp + EXPONENT_BIAS)) << SIGNIFICAND_BITS) 465 | (quot & SIGNIFICAND_MASK); 466 return Double.longBitsToDouble(raw); 467 } 468 469 private static final double unscaledToDoubleShiftAndDivideByScaleFactor(ScaleMetrics scaleMetrics, DecimalRounding rounding, long unscaled, int exp, int mantissaShift, long valModFactor) { 470 final long quotient; 471 final long scaleFactor = scaleMetrics.getScaleFactor(); 472 if (mantissaShift >= 0) { 473 final long hValModFactor = (valModFactor >>> (Long.SIZE - mantissaShift)) & (-mantissaShift >> 63); 474 final long lValModFactor = valModFactor << mantissaShift; 475 if (hValModFactor == 0) { 476 final long truncated = scaleMetrics.divideUnsignedByScaleFactor(lValModFactor); 477 final long remainder = applySign(unscaled, lValModFactor - scaleMetrics.multiplyByScaleFactor(truncated)); 478 quotient = truncated + Math.abs(Rounding.calculateRoundingIncrementForDivision(rounding, truncated, remainder, scaleFactor)); 479 } else { 480 quotient = Math.abs(Div.div128by64(rounding, unscaled < 0, hValModFactor, lValModFactor, scaleFactor)); 481 // rounding already done by div128by64 482 } 483 } else { 484 final long scaledVal = valModFactor >>> -mantissaShift; 485 final long truncated = scaleMetrics.divideByScaleFactor(scaledVal); 486 final long remainder = applySign(unscaled, ((scaledVal - scaleMetrics.multiplyByScaleFactor(truncated)) << -mantissaShift) 487 | (valModFactor & (-1L >>> (Long.SIZE + mantissaShift)))); 488 // this cannot overflow as min(mantissaShift)=-10 for scale=1, -9 for scale=10, ..., -1 for scale=10^9 489 final long shiftedScaleFactor = scaleFactor << -mantissaShift; 490 quotient = truncated + Math.abs(Rounding.calculateRoundingIncrementForDivision(rounding, truncated, remainder, 491 shiftedScaleFactor)); 492 } 493 final long raw; 494 if (quotient <= SIGNIFICAND_MASK) { 495 raw = (unscaled & SIGN_MASK) | (((long) (exp + EXPONENT_BIAS)) << SIGNIFICAND_BITS) 496 | (quotient & SIGNIFICAND_MASK); 497 } else { 498 // rounding made our value to be 1 instead of smaller than one. 1 + 499 // 1 == 2 i.e. our mantissa is zero due to the implicit 1 and our 500 // exponent increments by 1 501 raw = (unscaled & SIGN_MASK) | (((long) ((exp + 1) + EXPONENT_BIAS)) << SIGNIFICAND_BITS); 502 } 503 return Double.longBitsToDouble(raw); 504 } 505 506 // @return value % (2^n) 507 private static final long modPow2(long value, int n) { 508 // return value & ((1L << n) - 1); 509 return value & (-1L >>> (Long.SIZE - n)) & (-n >> 31);// last bracket is for case n=0 510 } 511 512 private static final long applySign(final long signed, final long value) { 513 return signed >= 0 ? value : -value; 514 } 515 516 private static final boolean isInLongRange(double value) { 517 return MIN_LONG_AS_DOUBLE - value < 1.0 & value < MAX_LONG_AS_DOUBLE_PLUS_ONE; 518 } 519 520 private static final boolean isMathematicalInteger(double x) { 521 return isFinite(x) && (x == 0.0 522 || SIGNIFICAND_BITS - Long.numberOfTrailingZeros(getSignificand(x)) <= Math.getExponent(x)); 523 } 524 525 private static final boolean isFinite(double d) { 526 return Math.abs(d) <= Double.MAX_VALUE; 527 } 528 529 // PRECONDITION: isFinite(d) 530 private static final long getSignificand(double d) { 531 int exponent = Math.getExponent(d); 532 long bits = Double.doubleToRawLongBits(d); 533 bits &= SIGNIFICAND_MASK; 534 return (exponent == Double.MIN_EXPONENT - 1) ? bits << 1 : bits | IMPLICIT_BIT; 535 } 536 537 private static final IllegalArgumentException newOverflowException(final DecimalArithmetic arith, double value) { 538 return new IllegalArgumentException( 539 "Overflow for conversion from double to decimal with scale " + arith.getScale() + ": " + value); 540 } 541 542 // no instances 543 private DoubleConversion() { 544 super(); 545 } 546}