1/* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17package org.apache.commons.math.util; 18 19/** 20 * Faster, more accurate, portable alternative to {@link StrictMath}. 21 * <p> 22 * Additionally implements the following methods not found in StrictMath: 23 * <ul> 24 * <li>{@link #asinh(double)}</li> 25 * <li>{@link #acosh(double)}</li> 26 * <li>{@link #atanh(double)}</li> 27 * </ul> 28 * The following methods are found in StrictMath since 1.6 only 29 * <ul> 30 * <li>{@link #copySign(double, double)}</li> 31 * <li>{@link #getExponent(double)}</li> 32 * <li>{@link #nextAfter(double,double)}</li> 33 * <li>{@link #nextUp(double)}</li> 34 * <li>{@link #scalb(double, int)}</li> 35 * <li>{@link #copySign(float, float)}</li> 36 * <li>{@link #getExponent(float)}</li> 37 * <li>{@link #nextAfter(float,double)}</li> 38 * <li>{@link #nextUp(float)}</li> 39 * <li>{@link #scalb(float, int)}</li> 40 * </ul> 41 * @version $Revision: 1074294 $ $Date: 2011-02-24 22:18:59 +0100 (jeu. 24 févr. 2011) $ 42 * @since 2.2 43 */ 44public class FastMath { 45 46 /** Archimede's constant PI, ratio of circle circumference to diameter. */ 47 public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9; 48 49 /** Napier's constant e, base of the natural logarithm. */ 50 public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8; 51 52 /** Exponential evaluated at integer values, 53 * exp(x) = expIntTableA[x + 750] + expIntTableB[x+750]. 54 */ 55 private static final double EXP_INT_TABLE_A[] = new double[1500]; 56 57 /** Exponential evaluated at integer values, 58 * exp(x) = expIntTableA[x + 750] + expIntTableB[x+750] 59 */ 60 private static final double EXP_INT_TABLE_B[] = new double[1500]; 61 62 /** Exponential over the range of 0 - 1 in increments of 2^-10 63 * exp(x/1024) = expFracTableA[x] + expFracTableB[x]. 64 */ 65 private static final double EXP_FRAC_TABLE_A[] = new double[1025]; 66 67 /** Exponential over the range of 0 - 1 in increments of 2^-10 68 * exp(x/1024) = expFracTableA[x] + expFracTableB[x]. 69 */ 70 private static final double EXP_FRAC_TABLE_B[] = new double[1025]; 71 72 /** Factorial table, for Taylor series expansions. */ 73 private static final double FACT[] = new double[20]; 74 75 /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */ 76 private static final double LN_MANT[][] = new double[1024][]; 77 78 /** log(2) (high bits). */ 79 private static final double LN_2_A = 0.693147063255310059; 80 81 /** log(2) (low bits). */ 82 private static final double LN_2_B = 1.17304635250823482e-7; 83 84 /** Coefficients for slowLog. */ 85 private static final double LN_SPLIT_COEF[][] = { 86 {2.0, 0.0}, 87 {0.6666666269302368, 3.9736429850260626E-8}, 88 {0.3999999761581421, 2.3841857910019882E-8}, 89 {0.2857142686843872, 1.7029898543501842E-8}, 90 {0.2222222089767456, 1.3245471311735498E-8}, 91 {0.1818181574344635, 2.4384203044354907E-8}, 92 {0.1538461446762085, 9.140260083262505E-9}, 93 {0.13333332538604736, 9.220590270857665E-9}, 94 {0.11764700710773468, 1.2393345855018391E-8}, 95 {0.10526403784751892, 8.251545029714408E-9}, 96 {0.0952233225107193, 1.2675934823758863E-8}, 97 {0.08713622391223907, 1.1430250008909141E-8}, 98 {0.07842259109020233, 2.404307984052299E-9}, 99 {0.08371849358081818, 1.176342548272881E-8}, 100 {0.030589580535888672, 1.2958646899018938E-9}, 101 {0.14982303977012634, 1.225743062930824E-8}, 102 }; 103 104 /** Coefficients for log, when input 0.99 < x < 1.01. */ 105 private static final double LN_QUICK_COEF[][] = { 106 {1.0, 5.669184079525E-24}, 107 {-0.25, -0.25}, 108 {0.3333333134651184, 1.986821492305628E-8}, 109 {-0.25, -6.663542893624021E-14}, 110 {0.19999998807907104, 1.1921056801463227E-8}, 111 {-0.1666666567325592, -7.800414592973399E-9}, 112 {0.1428571343421936, 5.650007086920087E-9}, 113 {-0.12502530217170715, -7.44321345601866E-11}, 114 {0.11113807559013367, 9.219544613762692E-9}, 115 }; 116 117 /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */ 118 private static final double LN_HI_PREC_COEF[][] = { 119 {1.0, -6.032174644509064E-23}, 120 {-0.25, -0.25}, 121 {0.3333333134651184, 1.9868161777724352E-8}, 122 {-0.2499999701976776, -2.957007209750105E-8}, 123 {0.19999954104423523, 1.5830993332061267E-10}, 124 {-0.16624879837036133, -2.6033824355191673E-8} 125 }; 126 127 /** Sine table (high bits). */ 128 private static final double SINE_TABLE_A[] = new double[14]; 129 130 /** Sine table (low bits). */ 131 private static final double SINE_TABLE_B[] = new double[14]; 132 133 /** Cosine table (high bits). */ 134 private static final double COSINE_TABLE_A[] = new double[14]; 135 136 /** Cosine table (low bits). */ 137 private static final double COSINE_TABLE_B[] = new double[14]; 138 139 /** Tangent table, used by atan() (high bits). */ 140 private static final double TANGENT_TABLE_A[] = new double[14]; 141 142 /** Tangent table, used by atan() (low bits). */ 143 private static final double TANGENT_TABLE_B[] = new double[14]; 144 145 /** Bits of 1/(2*pi), need for reducePayneHanek(). */ 146 private static final long RECIP_2PI[] = new long[] { 147 (0x28be60dbL << 32) | 0x9391054aL, 148 (0x7f09d5f4L << 32) | 0x7d4d3770L, 149 (0x36d8a566L << 32) | 0x4f10e410L, 150 (0x7f9458eaL << 32) | 0xf7aef158L, 151 (0x6dc91b8eL << 32) | 0x909374b8L, 152 (0x01924bbaL << 32) | 0x82746487L, 153 (0x3f877ac7L << 32) | 0x2c4a69cfL, 154 (0xba208d7dL << 32) | 0x4baed121L, 155 (0x3a671c09L << 32) | 0xad17df90L, 156 (0x4e64758eL << 32) | 0x60d4ce7dL, 157 (0x272117e2L << 32) | 0xef7e4a0eL, 158 (0xc7fe25ffL << 32) | 0xf7816603L, 159 (0xfbcbc462L << 32) | 0xd6829b47L, 160 (0xdb4d9fb3L << 32) | 0xc9f2c26dL, 161 (0xd3d18fd9L << 32) | 0xa797fa8bL, 162 (0x5d49eeb1L << 32) | 0xfaf97c5eL, 163 (0xcf41ce7dL << 32) | 0xe294a4baL, 164 0x9afed7ecL << 32 }; 165 166 /** Bits of pi/4, need for reducePayneHanek(). */ 167 private static final long PI_O_4_BITS[] = new long[] { 168 (0xc90fdaa2L << 32) | 0x2168c234L, 169 (0xc4c6628bL << 32) | 0x80dc1cd1L }; 170 171 /** Eighths. 172 * This is used by sinQ, because its faster to do a table lookup than 173 * a multiply in this time-critical routine 174 */ 175 private static final double EIGHTHS[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625}; 176 177 /** Table of 2^((n+2)/3) */ 178 private static final double CBRTTWO[] = { 0.6299605249474366, 179 0.7937005259840998, 180 1.0, 181 1.2599210498948732, 182 1.5874010519681994 }; 183 184 /* 185 * There are 52 bits in the mantissa of a double. 186 * For additional precision, the code splits double numbers into two parts, 187 * by clearing the low order 30 bits if possible, and then performs the arithmetic 188 * on each half separately. 189 */ 190 191 /** 192 * 0x40000000 - used to split a double into two parts, both with the low order bits cleared. 193 * Equivalent to 2^30. 194 */ 195 private static final long HEX_40000000 = 0x40000000L; // 1073741824L 196 197 /** Mask used to clear low order 30 bits */ 198 private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L; 199 200 /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */ 201 private static final double TWO_POWER_52 = 4503599627370496.0; 202 203 // Initialize tables 204 static { 205 int i; 206 207 // Generate an array of factorials 208 FACT[0] = 1.0; 209 for (i = 1; i < FACT.length; i++) { 210 FACT[i] = FACT[i-1] * i; 211 } 212 213 double tmp[] = new double[2]; 214 double recip[] = new double[2]; 215 216 // Populate expIntTable 217 for (i = 0; i < 750; i++) { 218 expint(i, tmp); 219 EXP_INT_TABLE_A[i+750] = tmp[0]; 220 EXP_INT_TABLE_B[i+750] = tmp[1]; 221 222 if (i != 0) { 223 // Negative integer powers 224 splitReciprocal(tmp, recip); 225 EXP_INT_TABLE_A[750-i] = recip[0]; 226 EXP_INT_TABLE_B[750-i] = recip[1]; 227 } 228 } 229 230 // Populate expFracTable 231 for (i = 0; i < EXP_FRAC_TABLE_A.length; i++) { 232 slowexp(i/1024.0, tmp); 233 EXP_FRAC_TABLE_A[i] = tmp[0]; 234 EXP_FRAC_TABLE_B[i] = tmp[1]; 235 } 236 237 // Populate lnMant table 238 for (i = 0; i < LN_MANT.length; i++) { 239 double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L ); 240 LN_MANT[i] = slowLog(d); 241 } 242 243 // Build the sine and cosine tables 244 buildSinCosTables(); 245 } 246 247 /** 248 * Private Constructor 249 */ 250 private FastMath() { 251 } 252 253 // Generic helper methods 254 255 /** 256 * Get the high order bits from the mantissa. 257 * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers 258 * 259 * @param d the value to split 260 * @return the high order part of the mantissa 261 */ 262 private static double doubleHighPart(double d) { 263 if (d > -MathUtils.SAFE_MIN && d < MathUtils.SAFE_MIN){ 264 return d; // These are un-normalised - don't try to convert 265 } 266 long xl = Double.doubleToLongBits(d); 267 xl = xl & MASK_30BITS; // Drop low order bits 268 return Double.longBitsToDouble(xl); 269 } 270 271 /** Compute the square root of a number. 272 * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt} 273 * @param a number on which evaluation is done 274 * @return square root of a 275 */ 276 public static double sqrt(final double a) { 277 return Math.sqrt(a); 278 } 279 280 /** Compute the hyperbolic cosine of a number. 281 * @param x number on which evaluation is done 282 * @return hyperbolic cosine of x 283 */ 284 public static double cosh(double x) { 285 if (x != x) { 286 return x; 287 } 288 289 if (x > 20.0) { 290 return exp(x)/2.0; 291 } 292 293 if (x < -20) { 294 return exp(-x)/2.0; 295 } 296 297 double hiPrec[] = new double[2]; 298 if (x < 0.0) { 299 x = -x; 300 } 301 exp(x, 0.0, hiPrec); 302 303 double ya = hiPrec[0] + hiPrec[1]; 304 double yb = -(ya - hiPrec[0] - hiPrec[1]); 305 306 double temp = ya * HEX_40000000; 307 double yaa = ya + temp - temp; 308 double yab = ya - yaa; 309 310 // recip = 1/y 311 double recip = 1.0/ya; 312 temp = recip * HEX_40000000; 313 double recipa = recip + temp - temp; 314 double recipb = recip - recipa; 315 316 // Correct for rounding in division 317 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip; 318 // Account for yb 319 recipb += -yb * recip * recip; 320 321 // y = y + 1/y 322 temp = ya + recipa; 323 yb += -(temp - ya - recipa); 324 ya = temp; 325 temp = ya + recipb; 326 yb += -(temp - ya - recipb); 327 ya = temp; 328 329 double result = ya + yb; 330 result *= 0.5; 331 return result; 332 } 333 334 /** Compute the hyperbolic sine of a number. 335 * @param x number on which evaluation is done 336 * @return hyperbolic sine of x 337 */ 338 public static double sinh(double x) { 339 boolean negate = false; 340 if (x != x) { 341 return x; 342 } 343 344 if (x > 20.0) { 345 return exp(x)/2.0; 346 } 347 348 if (x < -20) { 349 return -exp(-x)/2.0; 350 } 351 352 if (x == 0) { 353 return x; 354 } 355 356 if (x < 0.0) { 357 x = -x; 358 negate = true; 359 } 360 361 double result; 362 363 if (x > 0.25) { 364 double hiPrec[] = new double[2]; 365 exp(x, 0.0, hiPrec); 366 367 double ya = hiPrec[0] + hiPrec[1]; 368 double yb = -(ya - hiPrec[0] - hiPrec[1]); 369 370 double temp = ya * HEX_40000000; 371 double yaa = ya + temp - temp; 372 double yab = ya - yaa; 373 374 // recip = 1/y 375 double recip = 1.0/ya; 376 temp = recip * HEX_40000000; 377 double recipa = recip + temp - temp; 378 double recipb = recip - recipa; 379 380 // Correct for rounding in division 381 recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip; 382 // Account for yb 383 recipb += -yb * recip * recip; 384 385 recipa = -recipa; 386 recipb = -recipb; 387 388 // y = y + 1/y 389 temp = ya + recipa; 390 yb += -(temp - ya - recipa); 391 ya = temp; 392 temp = ya + recipb; 393 yb += -(temp - ya - recipb); 394 ya = temp; 395 396 result = ya + yb; 397 result *= 0.5; 398 } 399 else { 400 double hiPrec[] = new double[2]; 401 expm1(x, hiPrec); 402 403 double ya = hiPrec[0] + hiPrec[1]; 404 double yb = -(ya - hiPrec[0] - hiPrec[1]); 405 406 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */ 407 double denom = 1.0 + ya; 408 double denomr = 1.0 / denom; 409 double denomb = -(denom - 1.0 - ya) + yb; 410 double ratio = ya * denomr; 411 double temp = ratio * HEX_40000000; 412 double ra = ratio + temp - temp; 413 double rb = ratio - ra; 414 415 temp = denom * HEX_40000000; 416 double za = denom + temp - temp; 417 double zb = denom - za; 418 419 rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr; 420 421 // Adjust for yb 422 rb += yb*denomr; // numerator 423 rb += -ya * denomb * denomr * denomr; // denominator 424 425 // y = y - 1/y 426 temp = ya + ra; 427 yb += -(temp - ya - ra); 428 ya = temp; 429 temp = ya + rb; 430 yb += -(temp - ya - rb); 431 ya = temp; 432 433 result = ya + yb; 434 result *= 0.5; 435 } 436 437 if (negate) { 438 result = -result; 439 } 440 441 return result; 442 } 443 444 /** Compute the hyperbolic tangent of a number. 445 * @param x number on which evaluation is done 446 * @return hyperbolic tangent of x 447 */ 448 public static double tanh(double x) { 449 boolean negate = false; 450 451 if (x != x) { 452 return x; 453 } 454 455 if (x > 20.0) { 456 return 1.0; 457 } 458 459 if (x < -20) { 460 return -1.0; 461 } 462 463 if (x == 0) { 464 return x; 465 } 466 467 if (x < 0.0) { 468 x = -x; 469 negate = true; 470 } 471 472 double result; 473 if (x >= 0.5) { 474 double hiPrec[] = new double[2]; 475 // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1) 476 exp(x*2.0, 0.0, hiPrec); 477 478 double ya = hiPrec[0] + hiPrec[1]; 479 double yb = -(ya - hiPrec[0] - hiPrec[1]); 480 481 /* Numerator */ 482 double na = -1.0 + ya; 483 double nb = -(na + 1.0 - ya); 484 double temp = na + yb; 485 nb += -(temp - na - yb); 486 na = temp; 487 488 /* Denominator */ 489 double da = 1.0 + ya; 490 double db = -(da - 1.0 - ya); 491 temp = da + yb; 492 db += -(temp - da - yb); 493 da = temp; 494 495 temp = da * HEX_40000000; 496 double daa = da + temp - temp; 497 double dab = da - daa; 498 499 // ratio = na/da 500 double ratio = na/da; 501 temp = ratio * HEX_40000000; 502 double ratioa = ratio + temp - temp; 503 double ratiob = ratio - ratioa; 504 505 // Correct for rounding in division 506 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da; 507 508 // Account for nb 509 ratiob += nb / da; 510 // Account for db 511 ratiob += -db * na / da / da; 512 513 result = ratioa + ratiob; 514 } 515 else { 516 double hiPrec[] = new double[2]; 517 // tanh(x) = expm1(2x) / (expm1(2x) + 2) 518 expm1(x*2.0, hiPrec); 519 520 double ya = hiPrec[0] + hiPrec[1]; 521 double yb = -(ya - hiPrec[0] - hiPrec[1]); 522 523 /* Numerator */ 524 double na = ya; 525 double nb = yb; 526 527 /* Denominator */ 528 double da = 2.0 + ya; 529 double db = -(da - 2.0 - ya); 530 double temp = da + yb; 531 db += -(temp - da - yb); 532 da = temp; 533 534 temp = da * HEX_40000000; 535 double daa = da + temp - temp; 536 double dab = da - daa; 537 538 // ratio = na/da 539 double ratio = na/da; 540 temp = ratio * HEX_40000000; 541 double ratioa = ratio + temp - temp; 542 double ratiob = ratio - ratioa; 543 544 // Correct for rounding in division 545 ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da; 546 547 // Account for nb 548 ratiob += nb / da; 549 // Account for db 550 ratiob += -db * na / da / da; 551 552 result = ratioa + ratiob; 553 } 554 555 if (negate) { 556 result = -result; 557 } 558 559 return result; 560 } 561 562 /** Compute the inverse hyperbolic cosine of a number. 563 * @param a number on which evaluation is done 564 * @return inverse hyperbolic cosine of a 565 */ 566 public static double acosh(final double a) { 567 return FastMath.log(a + FastMath.sqrt(a * a - 1)); 568 } 569 570 /** Compute the inverse hyperbolic sine of a number. 571 * @param a number on which evaluation is done 572 * @return inverse hyperbolic sine of a 573 */ 574 public static double asinh(double a) { 575 576 boolean negative = false; 577 if (a < 0) { 578 negative = true; 579 a = -a; 580 } 581 582 double absAsinh; 583 if (a > 0.167) { 584 absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a); 585 } else { 586 final double a2 = a * a; 587 if (a > 0.097) { 588 absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0 - a2 * (1.0 / 15.0 - a2 * (1.0 / 17.0) * 15.0 / 16.0) * 13.0 / 14.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0); 589 } else if (a > 0.036) { 590 absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0 - a2 * (1.0 / 11.0 - a2 * (1.0 / 13.0) * 11.0 / 12.0) * 9.0 / 10.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0); 591 } else if (a > 0.0036) { 592 absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0 - a2 * (1 / 7.0 - a2 * (1 / 9.0) * 7.0 / 8.0) * 5.0 / 6.0) * 3.0 / 4.0) / 2.0); 593 } else { 594 absAsinh = a * (1 - a2 * (1 / 3.0 - a2 * (1 / 5.0) * 3.0 / 4.0) / 2.0); 595 } 596 } 597 598 return negative ? -absAsinh : absAsinh; 599 600 } 601 602 /** Compute the inverse hyperbolic tangent of a number. 603 * @param a number on which evaluation is done 604 * @return inverse hyperbolic tangent of a 605 */ 606 public static double atanh(double a) { 607 608 boolean negative = false; 609 if (a < 0) { 610 negative = true; 611 a = -a; 612 } 613 614 double absAtanh; 615 if (a > 0.15) { 616 absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a)); 617 } else { 618 final double a2 = a * a; 619 if (a > 0.087) { 620 absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0 + a2 * (1.0 / 15.0 + a2 * (1.0 / 17.0))))))))); 621 } else if (a > 0.031) { 622 absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0 + a2 * (1.0 / 11.0 + a2 * (1.0 / 13.0))))))); 623 } else if (a > 0.003) { 624 absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0 + a2 * (1.0 / 7.0 + a2 * (1.0 / 9.0))))); 625 } else { 626 absAtanh = a * (1 + a2 * (1.0 / 3.0 + a2 * (1.0 / 5.0))); 627 } 628 } 629 630 return negative ? -absAtanh : absAtanh; 631 632 } 633 634 /** Compute the signum of a number. 635 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise 636 * @param a number on which evaluation is done 637 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a 638 */ 639 public static double signum(final double a) { 640 return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a 641 } 642 643 /** Compute the signum of a number. 644 * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise 645 * @param a number on which evaluation is done 646 * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a 647 */ 648 public static float signum(final float a) { 649 return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a 650 } 651 652 /** Compute next number towards positive infinity. 653 * @param a number to which neighbor should be computed 654 * @return neighbor of a towards positive infinity 655 */ 656 public static double nextUp(final double a) { 657 return nextAfter(a, Double.POSITIVE_INFINITY); 658 } 659 660 /** Compute next number towards positive infinity. 661 * @param a number to which neighbor should be computed 662 * @return neighbor of a towards positive infinity 663 */ 664 public static float nextUp(final float a) { 665 return nextAfter(a, Float.POSITIVE_INFINITY); 666 } 667 668 /** Returns a pseudo-random number between 0.0 and 1.0. 669 * <p><b>Note:</b> this implementation currently delegates to {@link Math#random} 670 * @return a random number between 0.0 and 1.0 671 */ 672 public static double random() { 673 return Math.random(); 674 } 675 676 /** 677 * Exponential function. 678 * 679 * Computes exp(x), function result is nearly rounded. It will be correctly 680 * rounded to the theoretical value for 99.9% of input values, otherwise it will 681 * have a 1 UPL error. 682 * 683 * Method: 684 * Lookup intVal = exp(int(x)) 685 * Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 ); 686 * Compute z as the exponential of the remaining bits by a polynomial minus one 687 * exp(x) = intVal * fracVal * (1 + z) 688 * 689 * Accuracy: 690 * Calculation is done with 63 bits of precision, so result should be correctly 691 * rounded for 99.9% of input values, with less than 1 ULP error otherwise. 692 * 693 * @param x a double 694 * @return double e<sup>x</sup> 695 */ 696 public static double exp(double x) { 697 return exp(x, 0.0, null); 698 } 699 700 /** 701 * Internal helper method for exponential function. 702 * @param x original argument of the exponential function 703 * @param extra extra bits of precision on input (To Be Confirmed) 704 * @param hiPrec extra bits of precision on output (To Be Confirmed) 705 * @return exp(x) 706 */ 707 private static double exp(double x, double extra, double[] hiPrec) { 708 double intPartA; 709 double intPartB; 710 int intVal; 711 712 /* Lookup exp(floor(x)). 713 * intPartA will have the upper 22 bits, intPartB will have the lower 714 * 52 bits. 715 */ 716 if (x < 0.0) { 717 intVal = (int) -x; 718 719 if (intVal > 746) { 720 if (hiPrec != null) { 721 hiPrec[0] = 0.0; 722 hiPrec[1] = 0.0; 723 } 724 return 0.0; 725 } 726 727 if (intVal > 709) { 728 /* This will produce a subnormal output */ 729 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0; 730 if (hiPrec != null) { 731 hiPrec[0] /= 285040095144011776.0; 732 hiPrec[1] /= 285040095144011776.0; 733 } 734 return result; 735 } 736 737 if (intVal == 709) { 738 /* exp(1.494140625) is nearly a machine number... */ 739 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620; 740 if (hiPrec != null) { 741 hiPrec[0] /= 4.455505956692756620; 742 hiPrec[1] /= 4.455505956692756620; 743 } 744 return result; 745 } 746 747 intVal++; 748 749 intPartA = EXP_INT_TABLE_A[750-intVal]; 750 intPartB = EXP_INT_TABLE_B[750-intVal]; 751 752 intVal = -intVal; 753 } else { 754 intVal = (int) x; 755 756 if (intVal > 709) { 757 if (hiPrec != null) { 758 hiPrec[0] = Double.POSITIVE_INFINITY; 759 hiPrec[1] = 0.0; 760 } 761 return Double.POSITIVE_INFINITY; 762 } 763 764 intPartA = EXP_INT_TABLE_A[750+intVal]; 765 intPartB = EXP_INT_TABLE_B[750+intVal]; 766 } 767 768 /* Get the fractional part of x, find the greatest multiple of 2^-10 less than 769 * x and look up the exp function of it. 770 * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits. 771 */ 772 final int intFrac = (int) ((x - intVal) * 1024.0); 773 final double fracPartA = EXP_FRAC_TABLE_A[intFrac]; 774 final double fracPartB = EXP_FRAC_TABLE_B[intFrac]; 775 776 /* epsilon is the difference in x from the nearest multiple of 2^-10. It 777 * has a value in the range 0 <= epsilon < 2^-10. 778 * Do the subtraction from x as the last step to avoid possible loss of percison. 779 */ 780 final double epsilon = x - (intVal + intFrac / 1024.0); 781 782 /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial. z has 783 full double precision (52 bits). Since z < 2^-10, we will have 784 62 bits of precision when combined with the contant 1. This will be 785 used in the last addition below to get proper rounding. */ 786 787 /* Remez generated polynomial. Converges on the interval [0, 2^-10], error 788 is less than 0.5 ULP */ 789 double z = 0.04168701738764507; 790 z = z * epsilon + 0.1666666505023083; 791 z = z * epsilon + 0.5000000000042687; 792 z = z * epsilon + 1.0; 793 z = z * epsilon + -3.940510424527919E-20; 794 795 /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial 796 expansion. 797 tempA is exact since intPartA and intPartB only have 22 bits each. 798 tempB will have 52 bits of precision. 799 */ 800 double tempA = intPartA * fracPartA; 801 double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB; 802 803 /* Compute the result. (1+z)(tempA+tempB). Order of operations is 804 important. For accuracy add by increasing size. tempA is exact and 805 much larger than the others. If there are extra bits specified from the 806 pow() function, use them. */ 807 final double tempC = tempB + tempA; 808 final double result; 809 if (extra != 0.0) { 810 result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA; 811 } else { 812 result = tempC*z + tempB + tempA; 813 } 814 815 if (hiPrec != null) { 816 // If requesting high precision 817 hiPrec[0] = tempA; 818 hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB; 819 } 820 821 return result; 822 } 823 824 /** Compute exp(x) - 1 825 * @param x number to compute shifted exponential 826 * @return exp(x) - 1 827 */ 828 public static double expm1(double x) { 829 return expm1(x, null); 830 } 831 832 /** Internal helper method for expm1 833 * @param x number to compute shifted exponential 834 * @param hiPrecOut receive high precision result for -1.0 < x < 1.0 835 * @return exp(x) - 1 836 */ 837 private static double expm1(double x, double hiPrecOut[]) { 838 if (x != x || x == 0.0) { // NaN or zero 839 return x; 840 } 841 842 if (x <= -1.0 || x >= 1.0) { 843 // If not between +/- 1.0 844 //return exp(x) - 1.0; 845 double hiPrec[] = new double[2]; 846 exp(x, 0.0, hiPrec); 847 if (x > 0.0) { 848 return -1.0 + hiPrec[0] + hiPrec[1]; 849 } else { 850 final double ra = -1.0 + hiPrec[0]; 851 double rb = -(ra + 1.0 - hiPrec[0]); 852 rb += hiPrec[1]; 853 return ra + rb; 854 } 855 } 856 857 double baseA; 858 double baseB; 859 double epsilon; 860 boolean negative = false; 861 862 if (x < 0.0) { 863 x = -x; 864 negative = true; 865 } 866 867 { 868 int intFrac = (int) (x * 1024.0); 869 double tempA = EXP_FRAC_TABLE_A[intFrac] - 1.0; 870 double tempB = EXP_FRAC_TABLE_B[intFrac]; 871 872 double temp = tempA + tempB; 873 tempB = -(temp - tempA - tempB); 874 tempA = temp; 875 876 temp = tempA * HEX_40000000; 877 baseA = tempA + temp - temp; 878 baseB = tempB + (tempA - baseA); 879 880 epsilon = x - intFrac/1024.0; 881 } 882 883 884 /* Compute expm1(epsilon) */ 885 double zb = 0.008336750013465571; 886 zb = zb * epsilon + 0.041666663879186654; 887 zb = zb * epsilon + 0.16666666666745392; 888 zb = zb * epsilon + 0.49999999999999994; 889 zb = zb * epsilon; 890 zb = zb * epsilon; 891 892 double za = epsilon; 893 double temp = za + zb; 894 zb = -(temp - za - zb); 895 za = temp; 896 897 temp = za * HEX_40000000; 898 temp = za + temp - temp; 899 zb += za - temp; 900 za = temp; 901 902 /* Combine the parts. expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */ 903 double ya = za * baseA; 904 //double yb = za*baseB + zb*baseA + zb*baseB; 905 temp = ya + za * baseB; 906 double yb = -(temp - ya - za * baseB); 907 ya = temp; 908 909 temp = ya + zb * baseA; 910 yb += -(temp - ya - zb * baseA); 911 ya = temp; 912 913 temp = ya + zb * baseB; 914 yb += -(temp - ya - zb*baseB); 915 ya = temp; 916 917 //ya = ya + za + baseA; 918 //yb = yb + zb + baseB; 919 temp = ya + baseA; 920 yb += -(temp - baseA - ya); 921 ya = temp; 922 923 temp = ya + za; 924 //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya); 925 yb += -(temp - ya - za); 926 ya = temp; 927 928 temp = ya + baseB; 929 //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya); 930 yb += -(temp - ya - baseB); 931 ya = temp; 932 933 temp = ya + zb; 934 //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya); 935 yb += -(temp - ya - zb); 936 ya = temp; 937 938 if (negative) { 939 /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */ 940 double denom = 1.0 + ya; 941 double denomr = 1.0 / denom; 942 double denomb = -(denom - 1.0 - ya) + yb; 943 double ratio = ya * denomr; 944 temp = ratio * HEX_40000000; 945 final double ra = ratio + temp - temp; 946 double rb = ratio - ra; 947 948 temp = denom * HEX_40000000; 949 za = denom + temp - temp; 950 zb = denom - za; 951 952 rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr; 953 954 // f(x) = x/1+x 955 // Compute f'(x) 956 // Product rule: d(uv) = du*v + u*dv 957 // Chain rule: d(f(g(x)) = f'(g(x))*f(g'(x)) 958 // d(1/x) = -1/(x*x) 959 // d(1/1+x) = -1/( (1+x)^2) * 1 = -1/((1+x)*(1+x)) 960 // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x)) 961 962 // Adjust for yb 963 rb += yb * denomr; // numerator 964 rb += -ya * denomb * denomr * denomr; // denominator 965 966 // negate 967 ya = -ra; 968 yb = -rb; 969 } 970 971 if (hiPrecOut != null) { 972 hiPrecOut[0] = ya; 973 hiPrecOut[1] = yb; 974 } 975 976 return ya + yb; 977 } 978 979 /** 980 * For x between 0 and 1, returns exp(x), uses extended precision 981 * @param x argument of exponential 982 * @param result placeholder where to place exp(x) split in two terms 983 * for extra precision (i.e. exp(x) = result[0] ° result[1] 984 * @return exp(x) 985 */ 986 private static double slowexp(final double x, final double result[]) { 987 final double xs[] = new double[2]; 988 final double ys[] = new double[2]; 989 final double facts[] = new double[2]; 990 final double as[] = new double[2]; 991 split(x, xs); 992 ys[0] = ys[1] = 0.0; 993 994 for (int i = 19; i >= 0; i--) { 995 splitMult(xs, ys, as); 996 ys[0] = as[0]; 997 ys[1] = as[1]; 998 999 split(FACT[i], as); 1000 splitReciprocal(as, facts); 1001 1002 splitAdd(ys, facts, as); 1003 ys[0] = as[0]; 1004 ys[1] = as[1]; 1005 } 1006 1007 if (result != null) { 1008 result[0] = ys[0]; 1009 result[1] = ys[1]; 1010 } 1011 1012 return ys[0] + ys[1]; 1013 } 1014 1015 /** Compute split[0], split[1] such that their sum is equal to d, 1016 * and split[0] has its 30 least significant bits as zero. 1017 * @param d number to split 1018 * @param split placeholder where to place the result 1019 */ 1020 private static void split(final double d, final double split[]) { 1021 if (d < 8e298 && d > -8e298) { 1022 final double a = d * HEX_40000000; 1023 split[0] = (d + a) - a; 1024 split[1] = d - split[0]; 1025 } else { 1026 final double a = d * 9.31322574615478515625E-10; 1027 split[0] = (d + a - d) * HEX_40000000; 1028 split[1] = d - split[0]; 1029 } 1030 } 1031 1032 /** Recompute a split. 1033 * @param a input/out array containing the split, changed 1034 * on output 1035 */ 1036 private static void resplit(final double a[]) { 1037 final double c = a[0] + a[1]; 1038 final double d = -(c - a[0] - a[1]); 1039 1040 if (c < 8e298 && c > -8e298) { 1041 double z = c * HEX_40000000; 1042 a[0] = (c + z) - z; 1043 a[1] = c - a[0] + d; 1044 } else { 1045 double z = c * 9.31322574615478515625E-10; 1046 a[0] = (c + z - c) * HEX_40000000; 1047 a[1] = c - a[0] + d; 1048 } 1049 } 1050 1051 /** Multiply two numbers in split form. 1052 * @param a first term of multiplication 1053 * @param b second term of multiplication 1054 * @param ans placeholder where to put the result 1055 */ 1056 private static void splitMult(double a[], double b[], double ans[]) { 1057 ans[0] = a[0] * b[0]; 1058 ans[1] = a[0] * b[1] + a[1] * b[0] + a[1] * b[1]; 1059 1060 /* Resplit */ 1061 resplit(ans); 1062 } 1063 1064 /** Add two numbers in split form. 1065 * @param a first term of addition 1066 * @param b second term of addition 1067 * @param ans placeholder where to put the result 1068 */ 1069 private static void splitAdd(final double a[], final double b[], final double ans[]) { 1070 ans[0] = a[0] + b[0]; 1071 ans[1] = a[1] + b[1]; 1072 1073 resplit(ans); 1074 } 1075 1076 /** Compute the reciprocal of in. Use the following algorithm. 1077 * in = c + d. 1078 * want to find x + y such that x+y = 1/(c+d) and x is much 1079 * larger than y and x has several zero bits on the right. 1080 * 1081 * Set b = 1/(2^22), a = 1 - b. Thus (a+b) = 1. 1082 * Use following identity to compute (a+b)/(c+d) 1083 * 1084 * (a+b)/(c+d) = a/c + (bc - ad) / (c^2 + cd) 1085 * set x = a/c and y = (bc - ad) / (c^2 + cd) 1086 * This will be close to the right answer, but there will be 1087 * some rounding in the calculation of X. So by carefully 1088 * computing 1 - (c+d)(x+y) we can compute an error and 1089 * add that back in. This is done carefully so that terms 1090 * of similar size are subtracted first. 1091 * @param in initial number, in split form 1092 * @param result placeholder where to put the result 1093 */ 1094 private static void splitReciprocal(final double in[], final double result[]) { 1095 final double b = 1.0/4194304.0; 1096 final double a = 1.0 - b; 1097 1098 if (in[0] == 0.0) { 1099 in[0] = in[1]; 1100 in[1] = 0.0; 1101 } 1102 1103 result[0] = a / in[0]; 1104 result[1] = (b*in[0]-a*in[1]) / (in[0]*in[0] + in[0]*in[1]); 1105 1106 if (result[1] != result[1]) { // can happen if result[1] is NAN 1107 result[1] = 0.0; 1108 } 1109 1110 /* Resplit */ 1111 resplit(result); 1112 1113 for (int i = 0; i < 2; i++) { 1114 /* this may be overkill, probably once is enough */ 1115 double err = 1.0 - result[0] * in[0] - result[0] * in[1] - 1116 result[1] * in[0] - result[1] * in[1]; 1117 /*err = 1.0 - err; */ 1118 err = err * (result[0] + result[1]); 1119 /*printf("err = %16e\n", err); */ 1120 result[1] += err; 1121 } 1122 } 1123 1124 /** Compute (a[0] + a[1]) * (b[0] + b[1]) in extended precision. 1125 * @param a first term of the multiplication 1126 * @param b second term of the multiplication 1127 * @param result placeholder where to put the result 1128 */ 1129 private static void quadMult(final double a[], final double b[], final double result[]) { 1130 final double xs[] = new double[2]; 1131 final double ys[] = new double[2]; 1132 final double zs[] = new double[2]; 1133 1134 /* a[0] * b[0] */ 1135 split(a[0], xs); 1136 split(b[0], ys); 1137 splitMult(xs, ys, zs); 1138 1139 result[0] = zs[0]; 1140 result[1] = zs[1]; 1141 1142 /* a[0] * b[1] */ 1143 split(b[1], ys); 1144 splitMult(xs, ys, zs); 1145 1146 double tmp = result[0] + zs[0]; 1147 result[1] = result[1] - (tmp - result[0] - zs[0]); 1148 result[0] = tmp; 1149 tmp = result[0] + zs[1]; 1150 result[1] = result[1] - (tmp - result[0] - zs[1]); 1151 result[0] = tmp; 1152 1153 /* a[1] * b[0] */ 1154 split(a[1], xs); 1155 split(b[0], ys); 1156 splitMult(xs, ys, zs); 1157 1158 tmp = result[0] + zs[0]; 1159 result[1] = result[1] - (tmp - result[0] - zs[0]); 1160 result[0] = tmp; 1161 tmp = result[0] + zs[1]; 1162 result[1] = result[1] - (tmp - result[0] - zs[1]); 1163 result[0] = tmp; 1164 1165 /* a[1] * b[0] */ 1166 split(a[1], xs); 1167 split(b[1], ys); 1168 splitMult(xs, ys, zs); 1169 1170 tmp = result[0] + zs[0]; 1171 result[1] = result[1] - (tmp - result[0] - zs[0]); 1172 result[0] = tmp; 1173 tmp = result[0] + zs[1]; 1174 result[1] = result[1] - (tmp - result[0] - zs[1]); 1175 result[0] = tmp; 1176 } 1177 1178 /** Compute exp(p) for a integer p in extended precision. 1179 * @param p integer whose exponential is requested 1180 * @param result placeholder where to put the result in extended precision 1181 * @return exp(p) in standard precision (equal to result[0] + result[1]) 1182 */ 1183 private static double expint(int p, final double result[]) { 1184 //double x = M_E; 1185 final double xs[] = new double[2]; 1186 final double as[] = new double[2]; 1187 final double ys[] = new double[2]; 1188 //split(x, xs); 1189 //xs[1] = (double)(2.7182818284590452353602874713526625L - xs[0]); 1190 //xs[0] = 2.71827697753906250000; 1191 //xs[1] = 4.85091998273542816811e-06; 1192 //xs[0] = Double.longBitsToDouble(0x4005bf0800000000L); 1193 //xs[1] = Double.longBitsToDouble(0x3ed458a2bb4a9b00L); 1194 1195 /* E */ 1196 xs[0] = 2.718281828459045; 1197 xs[1] = 1.4456468917292502E-16; 1198 1199 split(1.0, ys); 1200 1201 while (p > 0) { 1202 if ((p & 1) != 0) { 1203 quadMult(ys, xs, as); 1204 ys[0] = as[0]; ys[1] = as[1]; 1205 } 1206 1207 quadMult(xs, xs, as); 1208 xs[0] = as[0]; xs[1] = as[1]; 1209 1210 p >>= 1; 1211 } 1212 1213 if (result != null) { 1214 result[0] = ys[0]; 1215 result[1] = ys[1]; 1216 1217 resplit(result); 1218 } 1219 1220 return ys[0] + ys[1]; 1221 } 1222 1223 1224 /** 1225 * Natural logarithm. 1226 * 1227 * @param x a double 1228 * @return log(x) 1229 */ 1230 public static double log(final double x) { 1231 return log(x, null); 1232 } 1233 1234 /** 1235 * Internal helper method for natural logarithm function. 1236 * @param x original argument of the natural logarithm function 1237 * @param hiPrec extra bits of precision on output (To Be Confirmed) 1238 * @return log(x) 1239 */ 1240 private static double log(final double x, final double[] hiPrec) { 1241 if (x==0) { // Handle special case of +0/-0 1242 return Double.NEGATIVE_INFINITY; 1243 } 1244 long bits = Double.doubleToLongBits(x); 1245 1246 /* Handle special cases of negative input, and NaN */ 1247 if ((bits & 0x8000000000000000L) != 0 || x != x) { 1248 if (x != 0.0) { 1249 if (hiPrec != null) { 1250 hiPrec[0] = Double.NaN; 1251 } 1252 1253 return Double.NaN; 1254 } 1255 } 1256 1257 /* Handle special cases of Positive infinity. */ 1258 if (x == Double.POSITIVE_INFINITY) { 1259 if (hiPrec != null) { 1260 hiPrec[0] = Double.POSITIVE_INFINITY; 1261 } 1262 1263 return Double.POSITIVE_INFINITY; 1264 } 1265 1266 /* Extract the exponent */ 1267 int exp = (int)(bits >> 52)-1023; 1268 1269 if ((bits & 0x7ff0000000000000L) == 0) { 1270 // Subnormal! 1271 if (x == 0) { 1272 // Zero 1273 if (hiPrec != null) { 1274 hiPrec[0] = Double.NEGATIVE_INFINITY; 1275 } 1276 1277 return Double.NEGATIVE_INFINITY; 1278 } 1279 1280 /* Normalize the subnormal number. */ 1281 bits <<= 1; 1282 while ( (bits & 0x0010000000000000L) == 0) { 1283 exp--; 1284 bits <<= 1; 1285 } 1286 } 1287 1288 1289 if (exp == -1 || exp == 0) { 1290 if (x < 1.01 && x > 0.99 && hiPrec == null) { 1291 /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight 1292 polynomial expansion in higer precision. */ 1293 1294 /* Compute x - 1.0 and split it */ 1295 double xa = x - 1.0; 1296 double xb = xa - x + 1.0; 1297 double tmp = xa * HEX_40000000; 1298 double aa = xa + tmp - tmp; 1299 double ab = xa - aa; 1300 xa = aa; 1301 xb = ab; 1302 1303 double ya = LN_QUICK_COEF[LN_QUICK_COEF.length-1][0]; 1304 double yb = LN_QUICK_COEF[LN_QUICK_COEF.length-1][1]; 1305 1306 for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) { 1307 /* Multiply a = y * x */ 1308 aa = ya * xa; 1309 ab = ya * xb + yb * xa + yb * xb; 1310 /* split, so now y = a */ 1311 tmp = aa * HEX_40000000; 1312 ya = aa + tmp - tmp; 1313 yb = aa - ya + ab; 1314 1315 /* Add a = y + lnQuickCoef */ 1316 aa = ya + LN_QUICK_COEF[i][0]; 1317 ab = yb + LN_QUICK_COEF[i][1]; 1318 /* Split y = a */ 1319 tmp = aa * HEX_40000000; 1320 ya = aa + tmp - tmp; 1321 yb = aa - ya + ab; 1322 } 1323 1324 /* Multiply a = y * x */ 1325 aa = ya * xa; 1326 ab = ya * xb + yb * xa + yb * xb; 1327 /* split, so now y = a */ 1328 tmp = aa * HEX_40000000; 1329 ya = aa + tmp - tmp; 1330 yb = aa - ya + ab; 1331 1332 return ya + yb; 1333 } 1334 } 1335 1336 // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2) 1337 double lnm[] = LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)]; 1338 1339 /* 1340 double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L); 1341 1342 epsilon -= 1.0; 1343 */ 1344 1345 // y is the most significant 10 bits of the mantissa 1346 //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L); 1347 //double epsilon = (x - y) / y; 1348 double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L)); 1349 1350 double lnza = 0.0; 1351 double lnzb = 0.0; 1352 1353 if (hiPrec != null) { 1354 /* split epsilon -> x */ 1355 double tmp = epsilon * HEX_40000000; 1356 double aa = epsilon + tmp - tmp; 1357 double ab = epsilon - aa; 1358 double xa = aa; 1359 double xb = ab; 1360 1361 /* Need a more accurate epsilon, so adjust the division. */ 1362 double numer = bits & 0x3ffffffffffL; 1363 double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L); 1364 aa = numer - xa*denom - xb * denom; 1365 xb += aa / denom; 1366 1367 /* Remez polynomial evaluation */ 1368 double ya = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][0]; 1369 double yb = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][1]; 1370 1371 for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) { 1372 /* Multiply a = y * x */ 1373 aa = ya * xa; 1374 ab = ya * xb + yb * xa + yb * xb; 1375 /* split, so now y = a */ 1376 tmp = aa * HEX_40000000; 1377 ya = aa + tmp - tmp; 1378 yb = aa - ya + ab; 1379 1380 /* Add a = y + lnHiPrecCoef */ 1381 aa = ya + LN_HI_PREC_COEF[i][0]; 1382 ab = yb + LN_HI_PREC_COEF[i][1]; 1383 /* Split y = a */ 1384 tmp = aa * HEX_40000000; 1385 ya = aa + tmp - tmp; 1386 yb = aa - ya + ab; 1387 } 1388 1389 /* Multiply a = y * x */ 1390 aa = ya * xa; 1391 ab = ya * xb + yb * xa + yb * xb; 1392 1393 /* split, so now lnz = a */ 1394 /* 1395 tmp = aa * 1073741824.0; 1396 lnza = aa + tmp - tmp; 1397 lnzb = aa - lnza + ab; 1398 */ 1399 lnza = aa + ab; 1400 lnzb = -(lnza - aa - ab); 1401 } else { 1402 /* High precision not required. Eval Remez polynomial 1403 using standard double precision */ 1404 lnza = -0.16624882440418567; 1405 lnza = lnza * epsilon + 0.19999954120254515; 1406 lnza = lnza * epsilon + -0.2499999997677497; 1407 lnza = lnza * epsilon + 0.3333333333332802; 1408 lnza = lnza * epsilon + -0.5; 1409 lnza = lnza * epsilon + 1.0; 1410 lnza = lnza * epsilon; 1411 } 1412 1413 /* Relative sizes: 1414 * lnzb [0, 2.33E-10] 1415 * lnm[1] [0, 1.17E-7] 1416 * ln2B*exp [0, 1.12E-4] 1417 * lnza [0, 9.7E-4] 1418 * lnm[0] [0, 0.692] 1419 * ln2A*exp [0, 709] 1420 */ 1421 1422 /* Compute the following sum: 1423 * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp; 1424 */ 1425 1426 //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp; 1427 double a = LN_2_A*exp; 1428 double b = 0.0; 1429 double c = a+lnm[0]; 1430 double d = -(c-a-lnm[0]); 1431 a = c; 1432 b = b + d; 1433 1434 c = a + lnza; 1435 d = -(c - a - lnza); 1436 a = c; 1437 b = b + d; 1438 1439 c = a + LN_2_B*exp; 1440 d = -(c - a - LN_2_B*exp); 1441 a = c; 1442 b = b + d; 1443 1444 c = a + lnm[1]; 1445 d = -(c - a - lnm[1]); 1446 a = c; 1447 b = b + d; 1448 1449 c = a + lnzb; 1450 d = -(c - a - lnzb); 1451 a = c; 1452 b = b + d; 1453 1454 if (hiPrec != null) { 1455 hiPrec[0] = a; 1456 hiPrec[1] = b; 1457 } 1458 1459 return a + b; 1460 } 1461 1462 /** Compute log(1 + x). 1463 * @param x a number 1464 * @return log(1 + x) 1465 */ 1466 public static double log1p(final double x) { 1467 double xpa = 1.0 + x; 1468 double xpb = -(xpa - 1.0 - x); 1469 1470 if (x == -1) { 1471 return x/0.0; // -Infinity 1472 } 1473 1474 if (x > 0 && 1/x == 0) { // x = Infinity 1475 return x; 1476 } 1477 1478 if (x>1e-6 || x<-1e-6) { 1479 double hiPrec[] = new double[2]; 1480 1481 final double lores = log(xpa, hiPrec); 1482 if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN 1483 return lores; 1484 } 1485 1486 /* Do a taylor series expansion around xpa */ 1487 /* f(x+y) = f(x) + f'(x)*y + f''(x)/2 y^2 */ 1488 double fx1 = xpb/xpa; 1489 1490 double epsilon = 0.5 * fx1 + 1.0; 1491 epsilon = epsilon * fx1; 1492 1493 return epsilon + hiPrec[1] + hiPrec[0]; 1494 } 1495 1496 /* Value is small |x| < 1e6, do a Taylor series centered on 1.0 */ 1497 double y = x * 0.333333333333333 - 0.5; 1498 y = y * x + 1.0; 1499 y = y * x; 1500 1501 return y; 1502 } 1503 1504 /** Compute the base 10 logarithm. 1505 * @param x a number 1506 * @return log10(x) 1507 */ 1508 public static double log10(final double x) { 1509 final double hiPrec[] = new double[2]; 1510 1511 final double lores = log(x, hiPrec); 1512 if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN 1513 return lores; 1514 } 1515 1516 final double tmp = hiPrec[0] * HEX_40000000; 1517 final double lna = hiPrec[0] + tmp - tmp; 1518 final double lnb = hiPrec[0] - lna + hiPrec[1]; 1519 1520 final double rln10a = 0.4342944622039795; 1521 final double rln10b = 1.9699272335463627E-8; 1522 1523 return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna; 1524 } 1525 1526 /** 1527 * Power function. Compute x^y. 1528 * 1529 * @param x a double 1530 * @param y a double 1531 * @return double 1532 */ 1533 public static double pow(double x, double y) { 1534 final double lns[] = new double[2]; 1535 1536 if (y == 0.0) { 1537 return 1.0; 1538 } 1539 1540 if (x != x) { // X is NaN 1541 return x; 1542 } 1543 1544 1545 if (x == 0) { 1546 long bits = Double.doubleToLongBits(x); 1547 if ((bits & 0x8000000000000000L) != 0) { 1548 // -zero 1549 long yi = (long) y; 1550 1551 if (y < 0 && y == yi && (yi & 1) == 1) { 1552 return Double.NEGATIVE_INFINITY; 1553 } 1554 1555 if (y < 0 && y == yi && (yi & 1) == 1) { 1556 return -0.0; 1557 } 1558 1559 if (y > 0 && y == yi && (yi & 1) == 1) { 1560 return -0.0; 1561 } 1562 } 1563 1564 if (y < 0) { 1565 return Double.POSITIVE_INFINITY; 1566 } 1567 if (y > 0) { 1568 return 0.0; 1569 } 1570 1571 return Double.NaN; 1572 } 1573 1574 if (x == Double.POSITIVE_INFINITY) { 1575 if (y != y) { // y is NaN 1576 return y; 1577 } 1578 if (y < 0.0) { 1579 return 0.0; 1580 } else { 1581 return Double.POSITIVE_INFINITY; 1582 } 1583 } 1584 1585 if (y == Double.POSITIVE_INFINITY) { 1586 if (x * x == 1.0) 1587 return Double.NaN; 1588 1589 if (x * x > 1.0) { 1590 return Double.POSITIVE_INFINITY; 1591 } else { 1592 return 0.0; 1593 } 1594 } 1595 1596 if (x == Double.NEGATIVE_INFINITY) { 1597 if (y != y) { // y is NaN 1598 return y; 1599 } 1600 1601 if (y < 0) { 1602 long yi = (long) y; 1603 if (y == yi && (yi & 1) == 1) { 1604 return -0.0; 1605 } 1606 1607 return 0.0; 1608 } 1609 1610 if (y > 0) { 1611 long yi = (long) y; 1612 if (y == yi && (yi & 1) == 1) { 1613 return Double.NEGATIVE_INFINITY; 1614 } 1615 1616 return Double.POSITIVE_INFINITY; 1617 } 1618 } 1619 1620 if (y == Double.NEGATIVE_INFINITY) { 1621 1622 if (x * x == 1.0) { 1623 return Double.NaN; 1624 } 1625 1626 if (x * x < 1.0) { 1627 return Double.POSITIVE_INFINITY; 1628 } else { 1629 return 0.0; 1630 } 1631 } 1632 1633 /* Handle special case x<0 */ 1634 if (x < 0) { 1635 // y is an even integer in this case 1636 if (y >= TWO_POWER_52 || y <= -TWO_POWER_52) { 1637 return pow(-x, y); 1638 } 1639 1640 if (y == (long) y) { 1641 // If y is an integer 1642 return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y); 1643 } else { 1644 return Double.NaN; 1645 } 1646 } 1647 1648 /* Split y into ya and yb such that y = ya+yb */ 1649 double ya; 1650 double yb; 1651 if (y < 8e298 && y > -8e298) { 1652 double tmp1 = y * HEX_40000000; 1653 ya = y + tmp1 - tmp1; 1654 yb = y - ya; 1655 } else { 1656 double tmp1 = y * 9.31322574615478515625E-10; 1657 double tmp2 = tmp1 * 9.31322574615478515625E-10; 1658 ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000; 1659 yb = y - ya; 1660 } 1661 1662 /* Compute ln(x) */ 1663 final double lores = log(x, lns); 1664 if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN 1665 return lores; 1666 } 1667 1668 double lna = lns[0]; 1669 double lnb = lns[1]; 1670 1671 /* resplit lns */ 1672 double tmp1 = lna * HEX_40000000; 1673 double tmp2 = lna + tmp1 - tmp1; 1674 lnb += lna - tmp2; 1675 lna = tmp2; 1676 1677 // y*ln(x) = (aa+ab) 1678 final double aa = lna * ya; 1679 final double ab = lna * yb + lnb * ya + lnb * yb; 1680 1681 lna = aa+ab; 1682 lnb = -(lna - aa - ab); 1683 1684 double z = 1.0 / 120.0; 1685 z = z * lnb + (1.0 / 24.0); 1686 z = z * lnb + (1.0 / 6.0); 1687 z = z * lnb + 0.5; 1688 z = z * lnb + 1.0; 1689 z = z * lnb; 1690 1691 final double result = exp(lna, z, null); 1692 //result = result + result * z; 1693 return result; 1694 } 1695 1696 /** xi in the range of [1, 2]. 1697 * 3 5 7 1698 * x+1 / x x x \ 1699 * ln ----- = 2 * | x + ---- + ---- + ---- + ... | 1700 * 1-x \ 3 5 7 / 1701 * 1702 * So, compute a Remez approximation of the following function 1703 * 1704 * ln ((sqrt(x)+1)/(1-sqrt(x))) / x 1705 * 1706 * This will be an even function with only positive coefficents. 1707 * x is in the range [0 - 1/3]. 1708 * 1709 * Transform xi for input to the above function by setting 1710 * x = (xi-1)/(xi+1). Input to the polynomial is x^2, then 1711 * the result is multiplied by x. 1712 * @param xi number from which log is requested 1713 * @return log(xi) 1714 */ 1715 private static double[] slowLog(double xi) { 1716 double x[] = new double[2]; 1717 double x2[] = new double[2]; 1718 double y[] = new double[2]; 1719 double a[] = new double[2]; 1720 1721 split(xi, x); 1722 1723 /* Set X = (x-1)/(x+1) */ 1724 x[0] += 1.0; 1725 resplit(x); 1726 splitReciprocal(x, a); 1727 x[0] -= 2.0; 1728 resplit(x); 1729 splitMult(x, a, y); 1730 x[0] = y[0]; 1731 x[1] = y[1]; 1732 1733 /* Square X -> X2*/ 1734 splitMult(x, x, x2); 1735 1736 1737 //x[0] -= 1.0; 1738 //resplit(x); 1739 1740 y[0] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][0]; 1741 y[1] = LN_SPLIT_COEF[LN_SPLIT_COEF.length-1][1]; 1742 1743 for (int i = LN_SPLIT_COEF.length-2; i >= 0; i--) { 1744 splitMult(y, x2, a); 1745 y[0] = a[0]; 1746 y[1] = a[1]; 1747 splitAdd(y, LN_SPLIT_COEF[i], a); 1748 y[0] = a[0]; 1749 y[1] = a[1]; 1750 } 1751 1752 splitMult(y, x, a); 1753 y[0] = a[0]; 1754 y[1] = a[1]; 1755 1756 return y; 1757 } 1758 1759 /** 1760 * For x between 0 and pi/4 compute sine. 1761 * @param x number from which sine is requested 1762 * @param result placeholder where to put the result in extended precision 1763 * @return sin(x) 1764 */ 1765 private static double slowSin(final double x, final double result[]) { 1766 final double xs[] = new double[2]; 1767 final double ys[] = new double[2]; 1768 final double facts[] = new double[2]; 1769 final double as[] = new double[2]; 1770 split(x, xs); 1771 ys[0] = ys[1] = 0.0; 1772 1773 for (int i = 19; i >= 0; i--) { 1774 splitMult(xs, ys, as); 1775 ys[0] = as[0]; ys[1] = as[1]; 1776 1777 if ( (i & 1) == 0) { 1778 continue; 1779 } 1780 1781 split(FACT[i], as); 1782 splitReciprocal(as, facts); 1783 1784 if ( (i & 2) != 0 ) { 1785 facts[0] = -facts[0]; 1786 facts[1] = -facts[1]; 1787 } 1788 1789 splitAdd(ys, facts, as); 1790 ys[0] = as[0]; ys[1] = as[1]; 1791 } 1792 1793 if (result != null) { 1794 result[0] = ys[0]; 1795 result[1] = ys[1]; 1796 } 1797 1798 return ys[0] + ys[1]; 1799 } 1800 1801 /** 1802 * For x between 0 and pi/4 compute cosine 1803 * @param x number from which cosine is requested 1804 * @param result placeholder where to put the result in extended precision 1805 * @return cos(x) 1806 */ 1807 private static double slowCos(final double x, final double result[]) { 1808 1809 final double xs[] = new double[2]; 1810 final double ys[] = new double[2]; 1811 final double facts[] = new double[2]; 1812 final double as[] = new double[2]; 1813 split(x, xs); 1814 ys[0] = ys[1] = 0.0; 1815 1816 for (int i = 19; i >= 0; i--) { 1817 splitMult(xs, ys, as); 1818 ys[0] = as[0]; ys[1] = as[1]; 1819 1820 if ( (i & 1) != 0) { 1821 continue; 1822 } 1823 1824 split(FACT[i], as); 1825 splitReciprocal(as, facts); 1826 1827 if ( (i & 2) != 0 ) { 1828 facts[0] = -facts[0]; 1829 facts[1] = -facts[1]; 1830 } 1831 1832 splitAdd(ys, facts, as); 1833 ys[0] = as[0]; ys[1] = as[1]; 1834 } 1835 1836 if (result != null) { 1837 result[0] = ys[0]; 1838 result[1] = ys[1]; 1839 } 1840 1841 return ys[0] + ys[1]; 1842 } 1843 1844 /** Build the sine and cosine tables. 1845 */ 1846 private static void buildSinCosTables() { 1847 final double result[] = new double[2]; 1848 1849 /* Use taylor series for 0 <= x <= 6/8 */ 1850 for (int i = 0; i < 7; i++) { 1851 double x = i / 8.0; 1852 1853 slowSin(x, result); 1854 SINE_TABLE_A[i] = result[0]; 1855 SINE_TABLE_B[i] = result[1]; 1856 1857 slowCos(x, result); 1858 COSINE_TABLE_A[i] = result[0]; 1859 COSINE_TABLE_B[i] = result[1]; 1860 } 1861 1862 /* Use angle addition formula to complete table to 13/8, just beyond pi/2 */ 1863 for (int i = 7; i < 14; i++) { 1864 double xs[] = new double[2]; 1865 double ys[] = new double[2]; 1866 double as[] = new double[2]; 1867 double bs[] = new double[2]; 1868 double temps[] = new double[2]; 1869 1870 if ( (i & 1) == 0) { 1871 // Even, use double angle 1872 xs[0] = SINE_TABLE_A[i/2]; 1873 xs[1] = SINE_TABLE_B[i/2]; 1874 ys[0] = COSINE_TABLE_A[i/2]; 1875 ys[1] = COSINE_TABLE_B[i/2]; 1876 1877 /* compute sine */ 1878 splitMult(xs, ys, result); 1879 SINE_TABLE_A[i] = result[0] * 2.0; 1880 SINE_TABLE_B[i] = result[1] * 2.0; 1881 1882 /* Compute cosine */ 1883 splitMult(ys, ys, as); 1884 splitMult(xs, xs, temps); 1885 temps[0] = -temps[0]; 1886 temps[1] = -temps[1]; 1887 splitAdd(as, temps, result); 1888 COSINE_TABLE_A[i] = result[0]; 1889 COSINE_TABLE_B[i] = result[1]; 1890 } else { 1891 xs[0] = SINE_TABLE_A[i/2]; 1892 xs[1] = SINE_TABLE_B[i/2]; 1893 ys[0] = COSINE_TABLE_A[i/2]; 1894 ys[1] = COSINE_TABLE_B[i/2]; 1895 as[0] = SINE_TABLE_A[i/2+1]; 1896 as[1] = SINE_TABLE_B[i/2+1]; 1897 bs[0] = COSINE_TABLE_A[i/2+1]; 1898 bs[1] = COSINE_TABLE_B[i/2+1]; 1899 1900 /* compute sine */ 1901 splitMult(xs, bs, temps); 1902 splitMult(ys, as, result); 1903 splitAdd(result, temps, result); 1904 SINE_TABLE_A[i] = result[0]; 1905 SINE_TABLE_B[i] = result[1]; 1906 1907 /* Compute cosine */ 1908 splitMult(ys, bs, result); 1909 splitMult(xs, as, temps); 1910 temps[0] = -temps[0]; 1911 temps[1] = -temps[1]; 1912 splitAdd(result, temps, result); 1913 COSINE_TABLE_A[i] = result[0]; 1914 COSINE_TABLE_B[i] = result[1]; 1915 } 1916 } 1917 1918 /* Compute tangent = sine/cosine */ 1919 for (int i = 0; i < 14; i++) { 1920 double xs[] = new double[2]; 1921 double ys[] = new double[2]; 1922 double as[] = new double[2]; 1923 1924 as[0] = COSINE_TABLE_A[i]; 1925 as[1] = COSINE_TABLE_B[i]; 1926 1927 splitReciprocal(as, ys); 1928 1929 xs[0] = SINE_TABLE_A[i]; 1930 xs[1] = SINE_TABLE_B[i]; 1931 1932 splitMult(xs, ys, as); 1933 1934 TANGENT_TABLE_A[i] = as[0]; 1935 TANGENT_TABLE_B[i] = as[1]; 1936 } 1937 1938 } 1939 1940 /** 1941 * Computes sin(x) - x, where |x| < 1/16. 1942 * Use a Remez polynomial approximation. 1943 * @param x a number smaller than 1/16 1944 * @return sin(x) - x 1945 */ 1946 private static double polySine(final double x) 1947 { 1948 double x2 = x*x; 1949 1950 double p = 2.7553817452272217E-6; 1951 p = p * x2 + -1.9841269659586505E-4; 1952 p = p * x2 + 0.008333333333329196; 1953 p = p * x2 + -0.16666666666666666; 1954 //p *= x2; 1955 //p *= x; 1956 p = p * x2 * x; 1957 1958 return p; 1959 } 1960 1961 /** 1962 * Computes cos(x) - 1, where |x| < 1/16. 1963 * Use a Remez polynomial approximation. 1964 * @param x a number smaller than 1/16 1965 * @return cos(x) - 1 1966 */ 1967 private static double polyCosine(double x) { 1968 double x2 = x*x; 1969 1970 double p = 2.479773539153719E-5; 1971 p = p * x2 + -0.0013888888689039883; 1972 p = p * x2 + 0.041666666666621166; 1973 p = p * x2 + -0.49999999999999994; 1974 p *= x2; 1975 1976 return p; 1977 } 1978 1979 /** 1980 * Compute sine over the first quadrant (0 < x < pi/2). 1981 * Use combination of table lookup and rational polynomial expansion. 1982 * @param xa number from which sine is requested 1983 * @param xb extra bits for x (may be 0.0) 1984 * @return sin(xa + xb) 1985 */ 1986 private static double sinQ(double xa, double xb) { 1987 int idx = (int) ((xa * 8.0) + 0.5); 1988 final double epsilon = xa - EIGHTHS[idx]; //idx*0.125; 1989 1990 // Table lookups 1991 final double sintA = SINE_TABLE_A[idx]; 1992 final double sintB = SINE_TABLE_B[idx]; 1993 final double costA = COSINE_TABLE_A[idx]; 1994 final double costB = COSINE_TABLE_B[idx]; 1995 1996 // Polynomial eval of sin(epsilon), cos(epsilon) 1997 double sinEpsA = epsilon; 1998 double sinEpsB = polySine(epsilon); 1999 final double cosEpsA = 1.0; 2000 final double cosEpsB = polyCosine(epsilon); 2001 2002 // Split epsilon xa + xb = x 2003 final double temp = sinEpsA * HEX_40000000; 2004 double temp2 = (sinEpsA + temp) - temp; 2005 sinEpsB += sinEpsA - temp2; 2006 sinEpsA = temp2; 2007 2008 /* Compute sin(x) by angle addition formula */ 2009 double result; 2010 2011 /* Compute the following sum: 2012 * 2013 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + 2014 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; 2015 * 2016 * Ranges of elements 2017 * 2018 * xxxtA 0 PI/2 2019 * xxxtB -1.5e-9 1.5e-9 2020 * sinEpsA -0.0625 0.0625 2021 * sinEpsB -6e-11 6e-11 2022 * cosEpsA 1.0 2023 * cosEpsB 0 -0.0625 2024 * 2025 */ 2026 2027 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + 2028 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; 2029 2030 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB; 2031 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB; 2032 double a = 0; 2033 double b = 0; 2034 2035 double t = sintA; 2036 double c = a + t; 2037 double d = -(c - a - t); 2038 a = c; 2039 b = b + d; 2040 2041 t = costA * sinEpsA; 2042 c = a + t; 2043 d = -(c - a - t); 2044 a = c; 2045 b = b + d; 2046 2047 b = b + sintA * cosEpsB + costA * sinEpsB; 2048 /* 2049 t = sintA*cosEpsB; 2050 c = a + t; 2051 d = -(c - a - t); 2052 a = c; 2053 b = b + d; 2054 2055 t = costA*sinEpsB; 2056 c = a + t; 2057 d = -(c - a - t); 2058 a = c; 2059 b = b + d; 2060 */ 2061 2062 b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB; 2063 /* 2064 t = sintB; 2065 c = a + t; 2066 d = -(c - a - t); 2067 a = c; 2068 b = b + d; 2069 2070 t = costB*sinEpsA; 2071 c = a + t; 2072 d = -(c - a - t); 2073 a = c; 2074 b = b + d; 2075 2076 t = sintB*cosEpsB; 2077 c = a + t; 2078 d = -(c - a - t); 2079 a = c; 2080 b = b + d; 2081 2082 t = costB*sinEpsB; 2083 c = a + t; 2084 d = -(c - a - t); 2085 a = c; 2086 b = b + d; 2087 */ 2088 2089 if (xb != 0.0) { 2090 t = ((costA + costB) * (cosEpsA + cosEpsB) - 2091 (sintA + sintB) * (sinEpsA + sinEpsB)) * xb; // approximate cosine*xb 2092 c = a + t; 2093 d = -(c - a - t); 2094 a = c; 2095 b = b + d; 2096 } 2097 2098 result = a + b; 2099 2100 return result; 2101 } 2102 2103 /** 2104 * Compute cosine in the first quadrant by subtracting input from PI/2 and 2105 * then calling sinQ. This is more accurate as the input approaches PI/2. 2106 * @param xa number from which cosine is requested 2107 * @param xb extra bits for x (may be 0.0) 2108 * @return cos(xa + xb) 2109 */ 2110 private static double cosQ(double xa, double xb) { 2111 final double pi2a = 1.5707963267948966; 2112 final double pi2b = 6.123233995736766E-17; 2113 2114 final double a = pi2a - xa; 2115 double b = -(a - pi2a + xa); 2116 b += pi2b - xb; 2117 2118 return sinQ(a, b); 2119 } 2120 2121 /** 2122 * Compute tangent (or cotangent) over the first quadrant. 0 < x < pi/2 2123 * Use combination of table lookup and rational polynomial expansion. 2124 * @param xa number from which sine is requested 2125 * @param xb extra bits for x (may be 0.0) 2126 * @param cotanFlag if true, compute the cotangent instead of the tangent 2127 * @return tan(xa+xb) (or cotangent, depending on cotanFlag) 2128 */ 2129 private static double tanQ(double xa, double xb, boolean cotanFlag) { 2130 2131 int idx = (int) ((xa * 8.0) + 0.5); 2132 final double epsilon = xa - EIGHTHS[idx]; //idx*0.125; 2133 2134 // Table lookups 2135 final double sintA = SINE_TABLE_A[idx]; 2136 final double sintB = SINE_TABLE_B[idx]; 2137 final double costA = COSINE_TABLE_A[idx]; 2138 final double costB = COSINE_TABLE_B[idx]; 2139 2140 // Polynomial eval of sin(epsilon), cos(epsilon) 2141 double sinEpsA = epsilon; 2142 double sinEpsB = polySine(epsilon); 2143 final double cosEpsA = 1.0; 2144 final double cosEpsB = polyCosine(epsilon); 2145 2146 // Split epsilon xa + xb = x 2147 double temp = sinEpsA * HEX_40000000; 2148 double temp2 = (sinEpsA + temp) - temp; 2149 sinEpsB += sinEpsA - temp2; 2150 sinEpsA = temp2; 2151 2152 /* Compute sin(x) by angle addition formula */ 2153 2154 /* Compute the following sum: 2155 * 2156 * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + 2157 * sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; 2158 * 2159 * Ranges of elements 2160 * 2161 * xxxtA 0 PI/2 2162 * xxxtB -1.5e-9 1.5e-9 2163 * sinEpsA -0.0625 0.0625 2164 * sinEpsB -6e-11 6e-11 2165 * cosEpsA 1.0 2166 * cosEpsB 0 -0.0625 2167 * 2168 */ 2169 2170 //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB + 2171 // sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; 2172 2173 //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB; 2174 //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB; 2175 double a = 0; 2176 double b = 0; 2177 2178 // Compute sine 2179 double t = sintA; 2180 double c = a + t; 2181 double d = -(c - a - t); 2182 a = c; 2183 b = b + d; 2184 2185 t = costA*sinEpsA; 2186 c = a + t; 2187 d = -(c - a - t); 2188 a = c; 2189 b = b + d; 2190 2191 b = b + sintA*cosEpsB + costA*sinEpsB; 2192 b = b + sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB; 2193 2194 double sina = a + b; 2195 double sinb = -(sina - a - b); 2196 2197 // Compute cosine 2198 2199 a = b = c = d = 0.0; 2200 2201 t = costA*cosEpsA; 2202 c = a + t; 2203 d = -(c - a - t); 2204 a = c; 2205 b = b + d; 2206 2207 t = -sintA*sinEpsA; 2208 c = a + t; 2209 d = -(c - a - t); 2210 a = c; 2211 b = b + d; 2212 2213 b = b + costB*cosEpsA + costA*cosEpsB + costB*cosEpsB; 2214 b = b - (sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB); 2215 2216 double cosa = a + b; 2217 double cosb = -(cosa - a - b); 2218 2219 if (cotanFlag) { 2220 double tmp; 2221 tmp = cosa; cosa = sina; sina = tmp; 2222 tmp = cosb; cosb = sinb; sinb = tmp; 2223 } 2224 2225 2226 /* estimate and correct, compute 1.0/(cosa+cosb) */ 2227 /* 2228 double est = (sina+sinb)/(cosa+cosb); 2229 double err = (sina - cosa*est) + (sinb - cosb*est); 2230 est += err/(cosa+cosb); 2231 err = (sina - cosa*est) + (sinb - cosb*est); 2232 */ 2233 2234 // f(x) = 1/x, f'(x) = -1/x^2 2235 2236 double est = sina/cosa; 2237 2238 /* Split the estimate to get more accurate read on division rounding */ 2239 temp = est * HEX_40000000; 2240 double esta = (est + temp) - temp; 2241 double estb = est - esta; 2242 2243 temp = cosa * HEX_40000000; 2244 double cosaa = (cosa + temp) - temp; 2245 double cosab = cosa - cosaa; 2246 2247 //double err = (sina - est*cosa)/cosa; // Correction for division rounding 2248 double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa; // Correction for division rounding 2249 err += sinb/cosa; // Change in est due to sinb 2250 err += -sina * cosb / cosa / cosa; // Change in est due to cosb 2251 2252 if (xb != 0.0) { 2253 // tan' = 1 + tan^2 cot' = -(1 + cot^2) 2254 // Approximate impact of xb 2255 double xbadj = xb + est*est*xb; 2256 if (cotanFlag) { 2257 xbadj = -xbadj; 2258 } 2259 2260 err += xbadj; 2261 } 2262 2263 return est+err; 2264 } 2265 2266 /** Reduce the input argument using the Payne and Hanek method. 2267 * This is good for all inputs 0.0 < x < inf 2268 * Output is remainder after dividing by PI/2 2269 * The result array should contain 3 numbers. 2270 * result[0] is the integer portion, so mod 4 this gives the quadrant. 2271 * result[1] is the upper bits of the remainder 2272 * result[2] is the lower bits of the remainder 2273 * 2274 * @param x number to reduce 2275 * @param result placeholder where to put the result 2276 */ 2277 private static void reducePayneHanek(double x, double result[]) 2278 { 2279 /* Convert input double to bits */ 2280 long inbits = Double.doubleToLongBits(x); 2281 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023; 2282 2283 /* Convert to fixed point representation */ 2284 inbits &= 0x000fffffffffffffL; 2285 inbits |= 0x0010000000000000L; 2286 2287 /* Normalize input to be between 0.5 and 1.0 */ 2288 exponent++; 2289 inbits <<= 11; 2290 2291 /* Based on the exponent, get a shifted copy of recip2pi */ 2292 long shpi0; 2293 long shpiA; 2294 long shpiB; 2295 int idx = exponent >> 6; 2296 int shift = exponent - (idx << 6); 2297 2298 if (shift != 0) { 2299 shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift); 2300 shpi0 |= RECIP_2PI[idx] >>> (64-shift); 2301 shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift)); 2302 shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift)); 2303 } else { 2304 shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1]; 2305 shpiA = RECIP_2PI[idx]; 2306 shpiB = RECIP_2PI[idx+1]; 2307 } 2308 2309 /* Multiply input by shpiA */ 2310 long a = inbits >>> 32; 2311 long b = inbits & 0xffffffffL; 2312 2313 long c = shpiA >>> 32; 2314 long d = shpiA & 0xffffffffL; 2315 2316 long ac = a * c; 2317 long bd = b * d; 2318 long bc = b * c; 2319 long ad = a * d; 2320 2321 long prodB = bd + (ad << 32); 2322 long prodA = ac + (ad >>> 32); 2323 2324 boolean bita = (bd & 0x8000000000000000L) != 0; 2325 boolean bitb = (ad & 0x80000000L ) != 0; 2326 boolean bitsum = (prodB & 0x8000000000000000L) != 0; 2327 2328 /* Carry */ 2329 if ( (bita && bitb) || 2330 ((bita || bitb) && !bitsum) ) { 2331 prodA++; 2332 } 2333 2334 bita = (prodB & 0x8000000000000000L) != 0; 2335 bitb = (bc & 0x80000000L ) != 0; 2336 2337 prodB = prodB + (bc << 32); 2338 prodA = prodA + (bc >>> 32); 2339 2340 bitsum = (prodB & 0x8000000000000000L) != 0; 2341 2342 /* Carry */ 2343 if ( (bita && bitb) || 2344 ((bita || bitb) && !bitsum) ) { 2345 prodA++; 2346 } 2347 2348 /* Multiply input by shpiB */ 2349 c = shpiB >>> 32; 2350 d = shpiB & 0xffffffffL; 2351 ac = a * c; 2352 bc = b * c; 2353 ad = a * d; 2354 2355 /* Collect terms */ 2356 ac = ac + ((bc + ad) >>> 32); 2357 2358 bita = (prodB & 0x8000000000000000L) != 0; 2359 bitb = (ac & 0x8000000000000000L ) != 0; 2360 prodB += ac; 2361 bitsum = (prodB & 0x8000000000000000L) != 0; 2362 /* Carry */ 2363 if ( (bita && bitb) || 2364 ((bita || bitb) && !bitsum) ) { 2365 prodA++; 2366 } 2367 2368 /* Multiply by shpi0 */ 2369 c = shpi0 >>> 32; 2370 d = shpi0 & 0xffffffffL; 2371 2372 bd = b * d; 2373 bc = b * c; 2374 ad = a * d; 2375 2376 prodA += bd + ((bc + ad) << 32); 2377 2378 /* 2379 * prodA, prodB now contain the remainder as a fraction of PI. We want this as a fraction of 2380 * PI/2, so use the following steps: 2381 * 1.) multiply by 4. 2382 * 2.) do a fixed point muliply by PI/4. 2383 * 3.) Convert to floating point. 2384 * 4.) Multiply by 2 2385 */ 2386 2387 /* This identifies the quadrant */ 2388 int intPart = (int)(prodA >>> 62); 2389 2390 /* Multiply by 4 */ 2391 prodA <<= 2; 2392 prodA |= prodB >>> 62; 2393 prodB <<= 2; 2394 2395 /* Multiply by PI/4 */ 2396 a = prodA >>> 32; 2397 b = prodA & 0xffffffffL; 2398 2399 c = PI_O_4_BITS[0] >>> 32; 2400 d = PI_O_4_BITS[0] & 0xffffffffL; 2401 2402 ac = a * c; 2403 bd = b * d; 2404 bc = b * c; 2405 ad = a * d; 2406 2407 long prod2B = bd + (ad << 32); 2408 long prod2A = ac + (ad >>> 32); 2409 2410 bita = (bd & 0x8000000000000000L) != 0; 2411 bitb = (ad & 0x80000000L ) != 0; 2412 bitsum = (prod2B & 0x8000000000000000L) != 0; 2413 2414 /* Carry */ 2415 if ( (bita && bitb) || 2416 ((bita || bitb) && !bitsum) ) { 2417 prod2A++; 2418 } 2419 2420 bita = (prod2B & 0x8000000000000000L) != 0; 2421 bitb = (bc & 0x80000000L ) != 0; 2422 2423 prod2B = prod2B + (bc << 32); 2424 prod2A = prod2A + (bc >>> 32); 2425 2426 bitsum = (prod2B & 0x8000000000000000L) != 0; 2427 2428 /* Carry */ 2429 if ( (bita && bitb) || 2430 ((bita || bitb) && !bitsum) ) { 2431 prod2A++; 2432 } 2433 2434 /* Multiply input by pio4bits[1] */ 2435 c = PI_O_4_BITS[1] >>> 32; 2436 d = PI_O_4_BITS[1] & 0xffffffffL; 2437 ac = a * c; 2438 bc = b * c; 2439 ad = a * d; 2440 2441 /* Collect terms */ 2442 ac = ac + ((bc + ad) >>> 32); 2443 2444 bita = (prod2B & 0x8000000000000000L) != 0; 2445 bitb = (ac & 0x8000000000000000L ) != 0; 2446 prod2B += ac; 2447 bitsum = (prod2B & 0x8000000000000000L) != 0; 2448 /* Carry */ 2449 if ( (bita && bitb) || 2450 ((bita || bitb) && !bitsum) ) { 2451 prod2A++; 2452 } 2453 2454 /* Multiply inputB by pio4bits[0] */ 2455 a = prodB >>> 32; 2456 b = prodB & 0xffffffffL; 2457 c = PI_O_4_BITS[0] >>> 32; 2458 d = PI_O_4_BITS[0] & 0xffffffffL; 2459 ac = a * c; 2460 bc = b * c; 2461 ad = a * d; 2462 2463 /* Collect terms */ 2464 ac = ac + ((bc + ad) >>> 32); 2465 2466 bita = (prod2B & 0x8000000000000000L) != 0; 2467 bitb = (ac & 0x8000000000000000L ) != 0; 2468 prod2B += ac; 2469 bitsum = (prod2B & 0x8000000000000000L) != 0; 2470 /* Carry */ 2471 if ( (bita && bitb) || 2472 ((bita || bitb) && !bitsum) ) { 2473 prod2A++; 2474 } 2475 2476 /* Convert to double */ 2477 double tmpA = (prod2A >>> 12) / TWO_POWER_52; // High order 52 bits 2478 double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits 2479 2480 double sumA = tmpA + tmpB; 2481 double sumB = -(sumA - tmpA - tmpB); 2482 2483 /* Multiply by PI/2 and return */ 2484 result[0] = intPart; 2485 result[1] = sumA * 2.0; 2486 result[2] = sumB * 2.0; 2487 } 2488 2489 /** 2490 * Sine function. 2491 * @param x a number 2492 * @return sin(x) 2493 */ 2494 public static double sin(double x) { 2495 boolean negative = false; 2496 int quadrant = 0; 2497 double xa; 2498 double xb = 0.0; 2499 2500 /* Take absolute value of the input */ 2501 xa = x; 2502 if (x < 0) { 2503 negative = true; 2504 xa = -xa; 2505 } 2506 2507 /* Check for zero and negative zero */ 2508 if (xa == 0.0) { 2509 long bits = Double.doubleToLongBits(x); 2510 if (bits < 0) { 2511 return -0.0; 2512 } 2513 return 0.0; 2514 } 2515 2516 if (xa != xa || xa == Double.POSITIVE_INFINITY) { 2517 return Double.NaN; 2518 } 2519 2520 /* Perform any argument reduction */ 2521 if (xa > 3294198.0) { 2522 // PI * (2**20) 2523 // Argument too big for CodyWaite reduction. Must use 2524 // PayneHanek. 2525 double reduceResults[] = new double[3]; 2526 reducePayneHanek(xa, reduceResults); 2527 quadrant = ((int) reduceResults[0]) & 3; 2528 xa = reduceResults[1]; 2529 xb = reduceResults[2]; 2530 } else if (xa > 1.5707963267948966) { 2531 /* Inline the Cody/Waite reduction for performance */ 2532 2533 // Estimate k 2534 //k = (int)(xa / 1.5707963267948966); 2535 int k = (int)(xa * 0.6366197723675814); 2536 2537 // Compute remainder 2538 double remA; 2539 double remB; 2540 while (true) { 2541 double a = -k * 1.570796251296997; 2542 remA = xa + a; 2543 remB = -(remA - xa - a); 2544 2545 a = -k * 7.549789948768648E-8; 2546 double b = remA; 2547 remA = a + b; 2548 remB += -(remA - b - a); 2549 2550 a = -k * 6.123233995736766E-17; 2551 b = remA; 2552 remA = a + b; 2553 remB += -(remA - b - a); 2554 2555 if (remA > 0.0) 2556 break; 2557 2558 // Remainder is negative, so decrement k and try again. 2559 // This should only happen if the input is very close 2560 // to an even multiple of pi/2 2561 k--; 2562 } 2563 quadrant = k & 3; 2564 xa = remA; 2565 xb = remB; 2566 } 2567 2568 if (negative) { 2569 quadrant ^= 2; // Flip bit 1 2570 } 2571 2572 switch (quadrant) { 2573 case 0: 2574 return sinQ(xa, xb); 2575 case 1: 2576 return cosQ(xa, xb); 2577 case 2: 2578 return -sinQ(xa, xb); 2579 case 3: 2580 return -cosQ(xa, xb); 2581 default: 2582 return Double.NaN; 2583 } 2584 } 2585 2586 /** 2587 * Cosine function 2588 * @param x a number 2589 * @return cos(x) 2590 */ 2591 public static double cos(double x) { 2592 int quadrant = 0; 2593 2594 /* Take absolute value of the input */ 2595 double xa = x; 2596 if (x < 0) { 2597 xa = -xa; 2598 } 2599 2600 if (xa != xa || xa == Double.POSITIVE_INFINITY) { 2601 return Double.NaN; 2602 } 2603 2604 /* Perform any argument reduction */ 2605 double xb = 0; 2606 if (xa > 3294198.0) { 2607 // PI * (2**20) 2608 // Argument too big for CodyWaite reduction. Must use 2609 // PayneHanek. 2610 double reduceResults[] = new double[3]; 2611 reducePayneHanek(xa, reduceResults); 2612 quadrant = ((int) reduceResults[0]) & 3; 2613 xa = reduceResults[1]; 2614 xb = reduceResults[2]; 2615 } else if (xa > 1.5707963267948966) { 2616 /* Inline the Cody/Waite reduction for performance */ 2617 2618 // Estimate k 2619 //k = (int)(xa / 1.5707963267948966); 2620 int k = (int)(xa * 0.6366197723675814); 2621 2622 // Compute remainder 2623 double remA; 2624 double remB; 2625 while (true) { 2626 double a = -k * 1.570796251296997; 2627 remA = xa + a; 2628 remB = -(remA - xa - a); 2629 2630 a = -k * 7.549789948768648E-8; 2631 double b = remA; 2632 remA = a + b; 2633 remB += -(remA - b - a); 2634 2635 a = -k * 6.123233995736766E-17; 2636 b = remA; 2637 remA = a + b; 2638 remB += -(remA - b - a); 2639 2640 if (remA > 0.0) 2641 break; 2642 2643 // Remainder is negative, so decrement k and try again. 2644 // This should only happen if the input is very close 2645 // to an even multiple of pi/2 2646 k--; 2647 } 2648 quadrant = k & 3; 2649 xa = remA; 2650 xb = remB; 2651 } 2652 2653 //if (negative) 2654 // quadrant = (quadrant + 2) % 4; 2655 2656 switch (quadrant) { 2657 case 0: 2658 return cosQ(xa, xb); 2659 case 1: 2660 return -sinQ(xa, xb); 2661 case 2: 2662 return -cosQ(xa, xb); 2663 case 3: 2664 return sinQ(xa, xb); 2665 default: 2666 return Double.NaN; 2667 } 2668 } 2669 2670 /** 2671 * Tangent function 2672 * @param x a number 2673 * @return tan(x) 2674 */ 2675 public static double tan(double x) { 2676 boolean negative = false; 2677 int quadrant = 0; 2678 2679 /* Take absolute value of the input */ 2680 double xa = x; 2681 if (x < 0) { 2682 negative = true; 2683 xa = -xa; 2684 } 2685 2686 /* Check for zero and negative zero */ 2687 if (xa == 0.0) { 2688 long bits = Double.doubleToLongBits(x); 2689 if (bits < 0) { 2690 return -0.0; 2691 } 2692 return 0.0; 2693 } 2694 2695 if (xa != xa || xa == Double.POSITIVE_INFINITY) { 2696 return Double.NaN; 2697 } 2698 2699 /* Perform any argument reduction */ 2700 double xb = 0; 2701 if (xa > 3294198.0) { 2702 // PI * (2**20) 2703 // Argument too big for CodyWaite reduction. Must use 2704 // PayneHanek. 2705 double reduceResults[] = new double[3]; 2706 reducePayneHanek(xa, reduceResults); 2707 quadrant = ((int) reduceResults[0]) & 3; 2708 xa = reduceResults[1]; 2709 xb = reduceResults[2]; 2710 } else if (xa > 1.5707963267948966) { 2711 /* Inline the Cody/Waite reduction for performance */ 2712 2713 // Estimate k 2714 //k = (int)(xa / 1.5707963267948966); 2715 int k = (int)(xa * 0.6366197723675814); 2716 2717 // Compute remainder 2718 double remA; 2719 double remB; 2720 while (true) { 2721 double a = -k * 1.570796251296997; 2722 remA = xa + a; 2723 remB = -(remA - xa - a); 2724 2725 a = -k * 7.549789948768648E-8; 2726 double b = remA; 2727 remA = a + b; 2728 remB += -(remA - b - a); 2729 2730 a = -k * 6.123233995736766E-17; 2731 b = remA; 2732 remA = a + b; 2733 remB += -(remA - b - a); 2734 2735 if (remA > 0.0) 2736 break; 2737 2738 // Remainder is negative, so decrement k and try again. 2739 // This should only happen if the input is very close 2740 // to an even multiple of pi/2 2741 k--; 2742 } 2743 quadrant = k & 3; 2744 xa = remA; 2745 xb = remB; 2746 } 2747 2748 if (xa > 1.5) { 2749 // Accurracy suffers between 1.5 and PI/2 2750 final double pi2a = 1.5707963267948966; 2751 final double pi2b = 6.123233995736766E-17; 2752 2753 final double a = pi2a - xa; 2754 double b = -(a - pi2a + xa); 2755 b += pi2b - xb; 2756 2757 xa = a + b; 2758 xb = -(xa - a - b); 2759 quadrant ^= 1; 2760 negative ^= true; 2761 } 2762 2763 double result; 2764 if ((quadrant & 1) == 0) { 2765 result = tanQ(xa, xb, false); 2766 } else { 2767 result = -tanQ(xa, xb, true); 2768 } 2769 2770 if (negative) { 2771 result = -result; 2772 } 2773 2774 return result; 2775 } 2776 2777 /** 2778 * Arctangent function 2779 * @param x a number 2780 * @return atan(x) 2781 */ 2782 public static double atan(double x) { 2783 return atan(x, 0.0, false); 2784 } 2785 2786 /** Internal helper function to compute arctangent. 2787 * @param xa number from which arctangent is requested 2788 * @param xb extra bits for x (may be 0.0) 2789 * @param leftPlane if true, result angle must be put in the left half plane 2790 * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true) 2791 */ 2792 private static double atan(double xa, double xb, boolean leftPlane) { 2793 boolean negate = false; 2794 int idx; 2795 2796 if (xa == 0.0) { // Matches +/- 0.0; return correct sign 2797 return leftPlane ? copySign(Math.PI, xa) : xa; 2798 } 2799 2800 if (xa < 0) { 2801 // negative 2802 xa = -xa; 2803 xb = -xb; 2804 negate = true; 2805 } 2806 2807 if (xa > 1.633123935319537E16) { // Very large input 2808 return (negate ^ leftPlane) ? (-Math.PI/2.0) : (Math.PI/2.0); 2809 } 2810 2811 /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */ 2812 if (xa < 1.0) { 2813 idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5); 2814 } else { 2815 double temp = 1.0/xa; 2816 idx = (int) (-((-1.7168146928204136 * temp * temp + 8.0) * temp) + 13.07); 2817 } 2818 double epsA = xa - TANGENT_TABLE_A[idx]; 2819 double epsB = -(epsA - xa + TANGENT_TABLE_A[idx]); 2820 epsB += xb - TANGENT_TABLE_B[idx]; 2821 2822 double temp = epsA + epsB; 2823 epsB = -(temp - epsA - epsB); 2824 epsA = temp; 2825 2826 /* Compute eps = eps / (1.0 + xa*tangent) */ 2827 temp = xa * HEX_40000000; 2828 double ya = xa + temp - temp; 2829 double yb = xb + xa - ya; 2830 xa = ya; 2831 xb += yb; 2832 2833 //if (idx > 8 || idx == 0) 2834 if (idx == 0) { 2835 /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */ 2836 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]); 2837 double denom = 1.0 / (1.0 + (xa + xb) * (TANGENT_TABLE_A[idx] + TANGENT_TABLE_B[idx])); 2838 //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]); 2839 ya = epsA * denom; 2840 yb = epsB * denom; 2841 } else { 2842 double temp2 = xa * TANGENT_TABLE_A[idx]; 2843 double za = 1.0 + temp2; 2844 double zb = -(za - 1.0 - temp2); 2845 temp2 = xb * TANGENT_TABLE_A[idx] + xa * TANGENT_TABLE_B[idx]; 2846 temp = za + temp2; 2847 zb += -(temp - za - temp2); 2848 za = temp; 2849 2850 zb += xb * TANGENT_TABLE_B[idx]; 2851 ya = epsA / za; 2852 2853 temp = ya * HEX_40000000; 2854 final double yaa = (ya + temp) - temp; 2855 final double yab = ya - yaa; 2856 2857 temp = za * HEX_40000000; 2858 final double zaa = (za + temp) - temp; 2859 final double zab = za - zaa; 2860 2861 /* Correct for rounding in division */ 2862 yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za; 2863 2864 yb += -epsA * zb / za / za; 2865 yb += epsB / za; 2866 } 2867 2868 2869 epsA = ya; 2870 epsB = yb; 2871 2872 /* Evaluate polynomial */ 2873 double epsA2 = epsA*epsA; 2874 2875 /* 2876 yb = -0.09001346640161823; 2877 yb = yb * epsA2 + 0.11110718400605211; 2878 yb = yb * epsA2 + -0.1428571349122913; 2879 yb = yb * epsA2 + 0.19999999999273194; 2880 yb = yb * epsA2 + -0.33333333333333093; 2881 yb = yb * epsA2 * epsA; 2882 */ 2883 2884 yb = 0.07490822288864472; 2885 yb = yb * epsA2 + -0.09088450866185192; 2886 yb = yb * epsA2 + 0.11111095942313305; 2887 yb = yb * epsA2 + -0.1428571423679182; 2888 yb = yb * epsA2 + 0.19999999999923582; 2889 yb = yb * epsA2 + -0.33333333333333287; 2890 yb = yb * epsA2 * epsA; 2891 2892 2893 ya = epsA; 2894 2895 temp = ya + yb; 2896 yb = -(temp - ya - yb); 2897 ya = temp; 2898 2899 /* Add in effect of epsB. atan'(x) = 1/(1+x^2) */ 2900 yb += epsB / (1.0 + epsA * epsA); 2901 2902 double result; 2903 double resultb; 2904 2905 //result = yb + eighths[idx] + ya; 2906 double za = EIGHTHS[idx] + ya; 2907 double zb = -(za - EIGHTHS[idx] - ya); 2908 temp = za + yb; 2909 zb += -(temp - za - yb); 2910 za = temp; 2911 2912 result = za + zb; 2913 resultb = -(result - za - zb); 2914 2915 if (leftPlane) { 2916 // Result is in the left plane 2917 final double pia = 1.5707963267948966*2.0; 2918 final double pib = 6.123233995736766E-17*2.0; 2919 2920 za = pia - result; 2921 zb = -(za - pia + result); 2922 zb += pib - resultb; 2923 2924 result = za + zb; 2925 resultb = -(result - za - zb); 2926 } 2927 2928 2929 if (negate ^ leftPlane) { 2930 result = -result; 2931 } 2932 2933 return result; 2934 } 2935 2936 /** 2937 * Two arguments arctangent function 2938 * @param y ordinate 2939 * @param x abscissa 2940 * @return phase angle of point (x,y) between {@code -PI} and {@code PI} 2941 */ 2942 public static double atan2(double y, double x) { 2943 if (x !=x || y != y) { 2944 return Double.NaN; 2945 } 2946 2947 if (y == 0.0) { 2948 double result = x*y; 2949 double invx = 1.0/x; 2950 double invy = 1.0/y; 2951 2952 if (invx == 0.0) { // X is infinite 2953 if (x > 0) { 2954 return y; // return +/- 0.0 2955 } else { 2956 return copySign(Math.PI, y); 2957 } 2958 } 2959 2960 if (x < 0.0 || invx < 0.0) { 2961 if (y < 0.0 || invy < 0.0) { 2962 return -Math.PI; 2963 } else { 2964 return Math.PI; 2965 } 2966 } else { 2967 return result; 2968 } 2969 } 2970 2971 // y cannot now be zero 2972 2973 if (y == Double.POSITIVE_INFINITY) { 2974 if (x == Double.POSITIVE_INFINITY) { 2975 return Math.PI/4.0; 2976 } 2977 2978 if (x == Double.NEGATIVE_INFINITY) { 2979 return Math.PI*3.0/4.0; 2980 } 2981 2982 return Math.PI/2.0; 2983 } 2984 2985 if (y == Double.NEGATIVE_INFINITY) { 2986 if (x == Double.POSITIVE_INFINITY) { 2987 return -Math.PI/4.0; 2988 } 2989 2990 if (x == Double.NEGATIVE_INFINITY) { 2991 return -Math.PI*3.0/4.0; 2992 } 2993 2994 return -Math.PI/2.0; 2995 } 2996 2997 if (x == Double.POSITIVE_INFINITY) { 2998 if (y > 0.0 || 1/y > 0.0) { 2999 return 0.0; 3000 } 3001 3002 if (y < 0.0 || 1/y < 0.0) { 3003 return -0.0; 3004 } 3005 } 3006 3007 if (x == Double.NEGATIVE_INFINITY) 3008 { 3009 if (y > 0.0 || 1/y > 0.0) { 3010 return Math.PI; 3011 } 3012 3013 if (y < 0.0 || 1/y < 0.0) { 3014 return -Math.PI; 3015 } 3016 } 3017 3018 // Neither y nor x can be infinite or NAN here 3019 3020 if (x == 0) { 3021 if (y > 0.0 || 1/y > 0.0) { 3022 return Math.PI/2.0; 3023 } 3024 3025 if (y < 0.0 || 1/y < 0.0) { 3026 return -Math.PI/2.0; 3027 } 3028 } 3029 3030 // Compute ratio r = y/x 3031 final double r = y/x; 3032 if (Double.isInfinite(r)) { // bypass calculations that can create NaN 3033 return atan(r, 0, x < 0); 3034 } 3035 3036 double ra = doubleHighPart(r); 3037 double rb = r - ra; 3038 3039 // Split x 3040 final double xa = doubleHighPart(x); 3041 final double xb = x - xa; 3042 3043 rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x; 3044 3045 double temp = ra + rb; 3046 rb = -(temp - ra - rb); 3047 ra = temp; 3048 3049 if (ra == 0) { // Fix up the sign so atan works correctly 3050 ra = copySign(0.0, y); 3051 } 3052 3053 // Call atan 3054 double result = atan(ra, rb, x < 0); 3055 3056 return result; 3057 } 3058 3059 /** Compute the arc sine of a number. 3060 * @param x number on which evaluation is done 3061 * @return arc sine of x 3062 */ 3063 public static double asin(double x) { 3064 if (x != x) { 3065 return Double.NaN; 3066 } 3067 3068 if (x > 1.0 || x < -1.0) { 3069 return Double.NaN; 3070 } 3071 3072 if (x == 1.0) { 3073 return Math.PI/2.0; 3074 } 3075 3076 if (x == -1.0) { 3077 return -Math.PI/2.0; 3078 } 3079 3080 if (x == 0.0) { // Matches +/- 0.0; return correct sign 3081 return x; 3082 } 3083 3084 /* Compute asin(x) = atan(x/sqrt(1-x*x)) */ 3085 3086 /* Split x */ 3087 double temp = x * HEX_40000000; 3088 final double xa = x + temp - temp; 3089 final double xb = x - xa; 3090 3091 /* Square it */ 3092 double ya = xa*xa; 3093 double yb = xa*xb*2.0 + xb*xb; 3094 3095 /* Subtract from 1 */ 3096 ya = -ya; 3097 yb = -yb; 3098 3099 double za = 1.0 + ya; 3100 double zb = -(za - 1.0 - ya); 3101 3102 temp = za + yb; 3103 zb += -(temp - za - yb); 3104 za = temp; 3105 3106 /* Square root */ 3107 double y; 3108 y = sqrt(za); 3109 temp = y * HEX_40000000; 3110 ya = y + temp - temp; 3111 yb = y - ya; 3112 3113 /* Extend precision of sqrt */ 3114 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y); 3115 3116 /* Contribution of zb to sqrt */ 3117 double dx = zb / (2.0*y); 3118 3119 // Compute ratio r = x/y 3120 double r = x/y; 3121 temp = r * HEX_40000000; 3122 double ra = r + temp - temp; 3123 double rb = r - ra; 3124 3125 rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y; // Correct for rounding in division 3126 rb += -x * dx / y / y; // Add in effect additional bits of sqrt. 3127 3128 temp = ra + rb; 3129 rb = -(temp - ra - rb); 3130 ra = temp; 3131 3132 return atan(ra, rb, false); 3133 } 3134 3135 /** Compute the arc cosine of a number. 3136 * @param x number on which evaluation is done 3137 * @return arc cosine of x 3138 */ 3139 public static double acos(double x) { 3140 if (x != x) { 3141 return Double.NaN; 3142 } 3143 3144 if (x > 1.0 || x < -1.0) { 3145 return Double.NaN; 3146 } 3147 3148 if (x == -1.0) { 3149 return Math.PI; 3150 } 3151 3152 if (x == 1.0) { 3153 return 0.0; 3154 } 3155 3156 if (x == 0) { 3157 return Math.PI/2.0; 3158 } 3159 3160 /* Compute acos(x) = atan(sqrt(1-x*x)/x) */ 3161 3162 /* Split x */ 3163 double temp = x * HEX_40000000; 3164 final double xa = x + temp - temp; 3165 final double xb = x - xa; 3166 3167 /* Square it */ 3168 double ya = xa*xa; 3169 double yb = xa*xb*2.0 + xb*xb; 3170 3171 /* Subtract from 1 */ 3172 ya = -ya; 3173 yb = -yb; 3174 3175 double za = 1.0 + ya; 3176 double zb = -(za - 1.0 - ya); 3177 3178 temp = za + yb; 3179 zb += -(temp - za - yb); 3180 za = temp; 3181 3182 /* Square root */ 3183 double y = sqrt(za); 3184 temp = y * HEX_40000000; 3185 ya = y + temp - temp; 3186 yb = y - ya; 3187 3188 /* Extend precision of sqrt */ 3189 yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y); 3190 3191 /* Contribution of zb to sqrt */ 3192 yb += zb / (2.0*y); 3193 y = ya+yb; 3194 yb = -(y - ya - yb); 3195 3196 // Compute ratio r = y/x 3197 double r = y/x; 3198 3199 // Did r overflow? 3200 if (Double.isInfinite(r)) { // x is effectively zero 3201 return Math.PI/2; // so return the appropriate value 3202 } 3203 3204 double ra = doubleHighPart(r); 3205 double rb = r - ra; 3206 3207 rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x; // Correct for rounding in division 3208 rb += yb / x; // Add in effect additional bits of sqrt. 3209 3210 temp = ra + rb; 3211 rb = -(temp - ra - rb); 3212 ra = temp; 3213 3214 return atan(ra, rb, x<0); 3215 } 3216 3217 /** Compute the cubic root of a number. 3218 * @param x number on which evaluation is done 3219 * @return cubic root of x 3220 */ 3221 public static double cbrt(double x) { 3222 /* Convert input double to bits */ 3223 long inbits = Double.doubleToLongBits(x); 3224 int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023; 3225 boolean subnormal = false; 3226 3227 if (exponent == -1023) { 3228 if (x == 0) { 3229 return x; 3230 } 3231 3232 /* Subnormal, so normalize */ 3233 subnormal = true; 3234 x *= 1.8014398509481984E16; // 2^54 3235 inbits = Double.doubleToLongBits(x); 3236 exponent = (int) ((inbits >> 52) & 0x7ff) - 1023; 3237 } 3238 3239 if (exponent == 1024) { 3240 // Nan or infinity. Don't care which. 3241 return x; 3242 } 3243 3244 /* Divide the exponent by 3 */ 3245 int exp3 = exponent / 3; 3246 3247 /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */ 3248 double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) | 3249 (long)(((exp3 + 1023) & 0x7ff)) << 52); 3250 3251 /* This will be a number between 1 and 2 */ 3252 final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L); 3253 3254 /* Estimate the cube root of mant by polynomial */ 3255 double est = -0.010714690733195933; 3256 est = est * mant + 0.0875862700108075; 3257 est = est * mant + -0.3058015757857271; 3258 est = est * mant + 0.7249995199969751; 3259 est = est * mant + 0.5039018405998233; 3260 3261 est *= CBRTTWO[exponent % 3 + 2]; 3262 3263 // est should now be good to about 15 bits of precision. Do 2 rounds of 3264 // Newton's method to get closer, this should get us full double precision 3265 // Scale down x for the purpose of doing newtons method. This avoids over/under flows. 3266 final double xs = x / (p2*p2*p2); 3267 est += (xs - est*est*est) / (3*est*est); 3268 est += (xs - est*est*est) / (3*est*est); 3269 3270 // Do one round of Newton's method in extended precision to get the last bit right. 3271 double temp = est * HEX_40000000; 3272 double ya = est + temp - temp; 3273 double yb = est - ya; 3274 3275 double za = ya * ya; 3276 double zb = ya * yb * 2.0 + yb * yb; 3277 temp = za * HEX_40000000; 3278 double temp2 = za + temp - temp; 3279 zb += za - temp2; 3280 za = temp2; 3281 3282 zb = za * yb + ya * zb + zb * yb; 3283 za = za * ya; 3284 3285 double na = xs - za; 3286 double nb = -(na - xs + za); 3287 nb -= zb; 3288 3289 est += (na+nb)/(3*est*est); 3290 3291 /* Scale by a power of two, so this is exact. */ 3292 est *= p2; 3293 3294 if (subnormal) { 3295 est *= 3.814697265625E-6; // 2^-18 3296 } 3297 3298 return est; 3299 } 3300 3301 /** 3302 * Convert degrees to radians, with error of less than 0.5 ULP 3303 * @param x angle in degrees 3304 * @return x converted into radians 3305 */ 3306 public static double toRadians(double x) 3307 { 3308 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign 3309 return x; 3310 } 3311 3312 // These are PI/180 split into high and low order bits 3313 final double facta = 0.01745329052209854; 3314 final double factb = 1.997844754509471E-9; 3315 3316 double xa = doubleHighPart(x); 3317 double xb = x - xa; 3318 3319 double result = xb * factb + xb * facta + xa * factb + xa * facta; 3320 if (result == 0) { 3321 result = result * x; // ensure correct sign if calculation underflows 3322 } 3323 return result; 3324 } 3325 3326 /** 3327 * Convert radians to degrees, with error of less than 0.5 ULP 3328 * @param x angle in radians 3329 * @return x converted into degrees 3330 */ 3331 public static double toDegrees(double x) 3332 { 3333 if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign 3334 return x; 3335 } 3336 3337 // These are 180/PI split into high and low order bits 3338 final double facta = 57.2957763671875; 3339 final double factb = 3.145894820876798E-6; 3340 3341 double xa = doubleHighPart(x); 3342 double xb = x - xa; 3343 3344 return xb * factb + xb * facta + xa * factb + xa * facta; 3345 } 3346 3347 /** 3348 * Absolute value. 3349 * @param x number from which absolute value is requested 3350 * @return abs(x) 3351 */ 3352 public static int abs(final int x) { 3353 return (x < 0) ? -x : x; 3354 } 3355 3356 /** 3357 * Absolute value. 3358 * @param x number from which absolute value is requested 3359 * @return abs(x) 3360 */ 3361 public static long abs(final long x) { 3362 return (x < 0l) ? -x : x; 3363 } 3364 3365 /** 3366 * Absolute value. 3367 * @param x number from which absolute value is requested 3368 * @return abs(x) 3369 */ 3370 public static float abs(final float x) { 3371 return (x < 0.0f) ? -x : (x == 0.0f) ? 0.0f : x; // -0.0 => +0.0 3372 } 3373 3374 /** 3375 * Absolute value. 3376 * @param x number from which absolute value is requested 3377 * @return abs(x) 3378 */ 3379 public static double abs(double x) { 3380 return (x < 0.0) ? -x : (x == 0.0) ? 0.0 : x; // -0.0 => +0.0 3381 } 3382 3383 /** 3384 * Compute least significant bit (Unit in Last Position) for a number. 3385 * @param x number from which ulp is requested 3386 * @return ulp(x) 3387 */ 3388 public static double ulp(double x) { 3389 if (Double.isInfinite(x)) { 3390 return Double.POSITIVE_INFINITY; 3391 } 3392 return abs(x - Double.longBitsToDouble(Double.doubleToLongBits(x) ^ 1)); 3393 } 3394 3395 /** 3396 * Compute least significant bit (Unit in Last Position) for a number. 3397 * @param x number from which ulp is requested 3398 * @return ulp(x) 3399 */ 3400 public static float ulp(float x) { 3401 if (Float.isInfinite(x)) { 3402 return Float.POSITIVE_INFINITY; 3403 } 3404 return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1)); 3405 } 3406 3407 /** 3408 * Multiply a double number by a power of 2. 3409 * @param d number to multiply 3410 * @param n power of 2 3411 * @return d × 2<sup>n</sup> 3412 */ 3413 public static double scalb(final double d, final int n) { 3414 3415 // first simple and fast handling when 2^n can be represented using normal numbers 3416 if ((n > -1023) && (n < 1024)) { 3417 return d * Double.longBitsToDouble(((long) (n + 1023)) << 52); 3418 } 3419 3420 // handle special cases 3421 if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) { 3422 return d; 3423 } 3424 if (n < -2098) { 3425 return (d > 0) ? 0.0 : -0.0; 3426 } 3427 if (n > 2097) { 3428 return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY; 3429 } 3430 3431 // decompose d 3432 final long bits = Double.doubleToLongBits(d); 3433 final long sign = bits & 0x8000000000000000L; 3434 int exponent = ((int) (bits >>> 52)) & 0x7ff; 3435 long mantissa = bits & 0x000fffffffffffffL; 3436 3437 // compute scaled exponent 3438 int scaledExponent = exponent + n; 3439 3440 if (n < 0) { 3441 // we are really in the case n <= -1023 3442 if (scaledExponent > 0) { 3443 // both the input and the result are normal numbers, we only adjust the exponent 3444 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa); 3445 } else if (scaledExponent > -53) { 3446 // the input is a normal number and the result is a subnormal number 3447 3448 // recover the hidden mantissa bit 3449 mantissa = mantissa | (1L << 52); 3450 3451 // scales down complete mantissa, hence losing least significant bits 3452 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent)); 3453 mantissa = mantissa >>> (1 - scaledExponent); 3454 if (mostSignificantLostBit != 0) { 3455 // we need to add 1 bit to round up the result 3456 mantissa++; 3457 } 3458 return Double.longBitsToDouble(sign | mantissa); 3459 3460 } else { 3461 // no need to compute the mantissa, the number scales down to 0 3462 return (sign == 0L) ? 0.0 : -0.0; 3463 } 3464 } else { 3465 // we are really in the case n >= 1024 3466 if (exponent == 0) { 3467 3468 // the input number is subnormal, normalize it 3469 while ((mantissa >>> 52) != 1) { 3470 mantissa = mantissa << 1; 3471 --scaledExponent; 3472 } 3473 ++scaledExponent; 3474 mantissa = mantissa & 0x000fffffffffffffL; 3475 3476 if (scaledExponent < 2047) { 3477 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa); 3478 } else { 3479 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY; 3480 } 3481 3482 } else if (scaledExponent < 2047) { 3483 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa); 3484 } else { 3485 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY; 3486 } 3487 } 3488 3489 } 3490 3491 /** 3492 * Multiply a float number by a power of 2. 3493 * @param f number to multiply 3494 * @param n power of 2 3495 * @return f × 2<sup>n</sup> 3496 */ 3497 public static float scalb(final float f, final int n) { 3498 3499 // first simple and fast handling when 2^n can be represented using normal numbers 3500 if ((n > -127) && (n < 128)) { 3501 return f * Float.intBitsToFloat((n + 127) << 23); 3502 } 3503 3504 // handle special cases 3505 if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) { 3506 return f; 3507 } 3508 if (n < -277) { 3509 return (f > 0) ? 0.0f : -0.0f; 3510 } 3511 if (n > 276) { 3512 return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY; 3513 } 3514 3515 // decompose f 3516 final int bits = Float.floatToIntBits(f); 3517 final int sign = bits & 0x80000000; 3518 int exponent = (bits >>> 23) & 0xff; 3519 int mantissa = bits & 0x007fffff; 3520 3521 // compute scaled exponent 3522 int scaledExponent = exponent + n; 3523 3524 if (n < 0) { 3525 // we are really in the case n <= -127 3526 if (scaledExponent > 0) { 3527 // both the input and the result are normal numbers, we only adjust the exponent 3528 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa); 3529 } else if (scaledExponent > -24) { 3530 // the input is a normal number and the result is a subnormal number 3531 3532 // recover the hidden mantissa bit 3533 mantissa = mantissa | (1 << 23); 3534 3535 // scales down complete mantissa, hence losing least significant bits 3536 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent)); 3537 mantissa = mantissa >>> (1 - scaledExponent); 3538 if (mostSignificantLostBit != 0) { 3539 // we need to add 1 bit to round up the result 3540 mantissa++; 3541 } 3542 return Float.intBitsToFloat(sign | mantissa); 3543 3544 } else { 3545 // no need to compute the mantissa, the number scales down to 0 3546 return (sign == 0) ? 0.0f : -0.0f; 3547 } 3548 } else { 3549 // we are really in the case n >= 128 3550 if (exponent == 0) { 3551 3552 // the input number is subnormal, normalize it 3553 while ((mantissa >>> 23) != 1) { 3554 mantissa = mantissa << 1; 3555 --scaledExponent; 3556 } 3557 ++scaledExponent; 3558 mantissa = mantissa & 0x007fffff; 3559 3560 if (scaledExponent < 255) { 3561 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa); 3562 } else { 3563 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY; 3564 } 3565 3566 } else if (scaledExponent < 255) { 3567 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa); 3568 } else { 3569 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY; 3570 } 3571 } 3572 3573 } 3574 3575 /** 3576 * Get the next machine representable number after a number, moving 3577 * in the direction of another number. 3578 * <p> 3579 * The ordering is as follows (increasing): 3580 * <ul> 3581 * <li>-INFINITY</li> 3582 * <li>-MAX_VALUE</li> 3583 * <li>-MIN_VALUE</li> 3584 * <li>-0.0</li> 3585 * <li>+0.0</li> 3586 * <li>+MIN_VALUE</li> 3587 * <li>+MAX_VALUE</li> 3588 * <li>+INFINITY</li> 3589 * <li></li> 3590 * <p> 3591 * If arguments compare equal, then the second argument is returned. 3592 * <p> 3593 * If {@code direction} is greater than {@code d}, 3594 * the smallest machine representable number strictly greater than 3595 * {@code d} is returned; if less, then the largest representable number 3596 * strictly less than {@code d} is returned.</p> 3597 * <p> 3598 * If {@code d} is infinite and direction does not 3599 * bring it back to finite numbers, it is returned unchanged.</p> 3600 * 3601 * @param d base number 3602 * @param direction (the only important thing is whether 3603 * {@code direction} is greater or smaller than {@code d}) 3604 * @return the next machine representable number in the specified direction 3605 */ 3606 public static double nextAfter(double d, double direction) { 3607 3608 // handling of some important special cases 3609 if (Double.isNaN(d) || Double.isNaN(direction)) { 3610 return Double.NaN; 3611 } else if (d == direction) { 3612 return direction; 3613 } else if (Double.isInfinite(d)) { 3614 return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE; 3615 } else if (d == 0) { 3616 return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE; 3617 } 3618 // special cases MAX_VALUE to infinity and MIN_VALUE to 0 3619 // are handled just as normal numbers 3620 3621 final long bits = Double.doubleToLongBits(d); 3622 final long sign = bits & 0x8000000000000000L; 3623 if ((direction < d) ^ (sign == 0L)) { 3624 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1)); 3625 } else { 3626 return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1)); 3627 } 3628 3629 } 3630 3631 /** 3632 * Get the next machine representable number after a number, moving 3633 * in the direction of another number. 3634 * <p> 3635 * The ordering is as follows (increasing): 3636 * <ul> 3637 * <li>-INFINITY</li> 3638 * <li>-MAX_VALUE</li> 3639 * <li>-MIN_VALUE</li> 3640 * <li>-0.0</li> 3641 * <li>+0.0</li> 3642 * <li>+MIN_VALUE</li> 3643 * <li>+MAX_VALUE</li> 3644 * <li>+INFINITY</li> 3645 * <li></li> 3646 * <p> 3647 * If arguments compare equal, then the second argument is returned. 3648 * <p> 3649 * If {@code direction} is greater than {@code f}, 3650 * the smallest machine representable number strictly greater than 3651 * {@code f} is returned; if less, then the largest representable number 3652 * strictly less than {@code f} is returned.</p> 3653 * <p> 3654 * If {@code f} is infinite and direction does not 3655 * bring it back to finite numbers, it is returned unchanged.</p> 3656 * 3657 * @param f base number 3658 * @param direction (the only important thing is whether 3659 * {@code direction} is greater or smaller than {@code f}) 3660 * @return the next machine representable number in the specified direction 3661 */ 3662 public static float nextAfter(final float f, final double direction) { 3663 3664 // handling of some important special cases 3665 if (Double.isNaN(f) || Double.isNaN(direction)) { 3666 return Float.NaN; 3667 } else if (f == direction) { 3668 return (float) direction; 3669 } else if (Float.isInfinite(f)) { 3670 return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE; 3671 } else if (f == 0f) { 3672 return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE; 3673 } 3674 // special cases MAX_VALUE to infinity and MIN_VALUE to 0 3675 // are handled just as normal numbers 3676 3677 final int bits = Float.floatToIntBits(f); 3678 final int sign = bits & 0x80000000; 3679 if ((direction < f) ^ (sign == 0)) { 3680 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1)); 3681 } else { 3682 return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1)); 3683 } 3684 3685 } 3686 3687 /** Get the largest whole number smaller than x. 3688 * @param x number from which floor is requested 3689 * @return a double number f such that f is an integer f <= x < f + 1.0 3690 */ 3691 public static double floor(double x) { 3692 long y; 3693 3694 if (x != x) { // NaN 3695 return x; 3696 } 3697 3698 if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) { 3699 return x; 3700 } 3701 3702 y = (long) x; 3703 if (x < 0 && y != x) { 3704 y--; 3705 } 3706 3707 if (y == 0) { 3708 return x*y; 3709 } 3710 3711 return y; 3712 } 3713 3714 /** Get the smallest whole number larger than x. 3715 * @param x number from which ceil is requested 3716 * @return a double number c such that c is an integer c - 1.0 < x <= c 3717 */ 3718 public static double ceil(double x) { 3719 double y; 3720 3721 if (x != x) { // NaN 3722 return x; 3723 } 3724 3725 y = floor(x); 3726 if (y == x) { 3727 return y; 3728 } 3729 3730 y += 1.0; 3731 3732 if (y == 0) { 3733 return x*y; 3734 } 3735 3736 return y; 3737 } 3738 3739 /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers. 3740 * @param x number from which nearest whole number is requested 3741 * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5 3742 */ 3743 public static double rint(double x) { 3744 double y = floor(x); 3745 double d = x - y; 3746 3747 if (d > 0.5) { 3748 if (y == -1.0) { 3749 return -0.0; // Preserve sign of operand 3750 } 3751 return y+1.0; 3752 } 3753 if (d < 0.5) { 3754 return y; 3755 } 3756 3757 /* half way, round to even */ 3758 long z = (long) y; 3759 return (z & 1) == 0 ? y : y + 1.0; 3760 } 3761 3762 /** Get the closest long to x. 3763 * @param x number from which closest long is requested 3764 * @return closest long to x 3765 */ 3766 public static long round(double x) { 3767 return (long) floor(x + 0.5); 3768 } 3769 3770 /** Get the closest int to x. 3771 * @param x number from which closest int is requested 3772 * @return closest int to x 3773 */ 3774 public static int round(final float x) { 3775 return (int) floor(x + 0.5f); 3776 } 3777 3778 /** Compute the minimum of two values 3779 * @param a first value 3780 * @param b second value 3781 * @return a if a is lesser or equal to b, b otherwise 3782 */ 3783 public static int min(final int a, final int b) { 3784 return (a <= b) ? a : b; 3785 } 3786 3787 /** Compute the minimum of two values 3788 * @param a first value 3789 * @param b second value 3790 * @return a if a is lesser or equal to b, b otherwise 3791 */ 3792 public static long min(final long a, final long b) { 3793 return (a <= b) ? a : b; 3794 } 3795 3796 /** Compute the minimum of two values 3797 * @param a first value 3798 * @param b second value 3799 * @return a if a is lesser or equal to b, b otherwise 3800 */ 3801 public static float min(final float a, final float b) { 3802 if (a > b) { 3803 return b; 3804 } 3805 if (a < b) { 3806 return a; 3807 } 3808 /* if either arg is NaN, return NaN */ 3809 if (a != b) { 3810 return Float.NaN; 3811 } 3812 /* min(+0.0,-0.0) == -0.0 */ 3813 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */ 3814 int bits = Float.floatToRawIntBits(a); 3815 if (bits == 0x80000000) { 3816 return a; 3817 } 3818 return b; 3819 } 3820 3821 /** Compute the minimum of two values 3822 * @param a first value 3823 * @param b second value 3824 * @return a if a is lesser or equal to b, b otherwise 3825 */ 3826 public static double min(final double a, final double b) { 3827 if (a > b) { 3828 return b; 3829 } 3830 if (a < b) { 3831 return a; 3832 } 3833 /* if either arg is NaN, return NaN */ 3834 if (a != b) { 3835 return Double.NaN; 3836 } 3837 /* min(+0.0,-0.0) == -0.0 */ 3838 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */ 3839 long bits = Double.doubleToRawLongBits(a); 3840 if (bits == 0x8000000000000000L) { 3841 return a; 3842 } 3843 return b; 3844 } 3845 3846 /** Compute the maximum of two values 3847 * @param a first value 3848 * @param b second value 3849 * @return b if a is lesser or equal to b, a otherwise 3850 */ 3851 public static int max(final int a, final int b) { 3852 return (a <= b) ? b : a; 3853 } 3854 3855 /** Compute the maximum of two values 3856 * @param a first value 3857 * @param b second value 3858 * @return b if a is lesser or equal to b, a otherwise 3859 */ 3860 public static long max(final long a, final long b) { 3861 return (a <= b) ? b : a; 3862 } 3863 3864 /** Compute the maximum of two values 3865 * @param a first value 3866 * @param b second value 3867 * @return b if a is lesser or equal to b, a otherwise 3868 */ 3869 public static float max(final float a, final float b) { 3870 if (a > b) { 3871 return a; 3872 } 3873 if (a < b) { 3874 return b; 3875 } 3876 /* if either arg is NaN, return NaN */ 3877 if (a != b) { 3878 return Float.NaN; 3879 } 3880 /* min(+0.0,-0.0) == -0.0 */ 3881 /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */ 3882 int bits = Float.floatToRawIntBits(a); 3883 if (bits == 0x80000000) { 3884 return b; 3885 } 3886 return a; 3887 } 3888 3889 /** Compute the maximum of two values 3890 * @param a first value 3891 * @param b second value 3892 * @return b if a is lesser or equal to b, a otherwise 3893 */ 3894 public static double max(final double a, final double b) { 3895 if (a > b) { 3896 return a; 3897 } 3898 if (a < b) { 3899 return b; 3900 } 3901 /* if either arg is NaN, return NaN */ 3902 if (a != b) { 3903 return Double.NaN; 3904 } 3905 /* min(+0.0,-0.0) == -0.0 */ 3906 /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */ 3907 long bits = Double.doubleToRawLongBits(a); 3908 if (bits == 0x8000000000000000L) { 3909 return b; 3910 } 3911 return a; 3912 } 3913 3914 /** 3915 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y} 3916 * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br/> 3917 * avoiding intermediate overflow or underflow. 3918 * 3919 * <ul> 3920 * <li> If either argument is infinite, then the result is positive infinity.</li> 3921 * <li> else, if either argument is NaN then the result is NaN.</li> 3922 * </ul> 3923 * 3924 * @param x a value 3925 * @param y a value 3926 * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>) 3927 */ 3928 public static double hypot(final double x, final double y) { 3929 if (Double.isInfinite(x) || Double.isInfinite(y)) { 3930 return Double.POSITIVE_INFINITY; 3931 } else if (Double.isNaN(x) || Double.isNaN(y)) { 3932 return Double.NaN; 3933 } else { 3934 3935 final int expX = getExponent(x); 3936 final int expY = getExponent(y); 3937 if (expX > expY + 27) { 3938 // y is neglectible with respect to x 3939 return abs(x); 3940 } else if (expY > expX + 27) { 3941 // x is neglectible with respect to y 3942 return abs(y); 3943 } else { 3944 3945 // find an intermediate scale to avoid both overflow and underflow 3946 final int middleExp = (expX + expY) / 2; 3947 3948 // scale parameters without losing precision 3949 final double scaledX = scalb(x, -middleExp); 3950 final double scaledY = scalb(y, -middleExp); 3951 3952 // compute scaled hypotenuse 3953 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY); 3954 3955 // remove scaling 3956 return scalb(scaledH, middleExp); 3957 3958 } 3959 3960 } 3961 } 3962 3963 /** 3964 * Computes the remainder as prescribed by the IEEE 754 standard. 3965 * The remainder value is mathematically equal to {@code x - y*n} 3966 * where {@code n} is the mathematical integer closest to the exact mathematical value 3967 * of the quotient {@code x/y}. 3968 * If two mathematical integers are equally close to {@code x/y} then 3969 * {@code n} is the integer that is even. 3970 * <p> 3971 * <ul> 3972 * <li>If either operand is NaN, the result is NaN.</li> 3973 * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li> 3974 * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li> 3975 * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li> 3976 * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li> 3977 * </ul> 3978 * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder} 3979 * @param dividend the number to be divided 3980 * @param divisor the number by which to divide 3981 * @return the remainder, rounded 3982 */ 3983 public static double IEEEremainder(double dividend, double divisor) { 3984 return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation 3985 } 3986 3987 /** 3988 * Returns the first argument with the sign of the second argument. 3989 * A NaN {@code sign} argument is treated as positive. 3990 * 3991 * @param magnitude the value to return 3992 * @param sign the sign for the returned value 3993 * @return the magnitude with the same sign as the {@code sign} argument 3994 */ 3995 public static double copySign(double magnitude, double sign){ 3996 long m = Double.doubleToLongBits(magnitude); 3997 long s = Double.doubleToLongBits(sign); 3998 if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK 3999 return magnitude; 4000 } 4001 return -magnitude; // flip sign 4002 } 4003 4004 /** 4005 * Returns the first argument with the sign of the second argument. 4006 * A NaN {@code sign} argument is treated as positive. 4007 * 4008 * @param magnitude the value to return 4009 * @param sign the sign for the returned value 4010 * @return the magnitude with the same sign as the {@code sign} argument 4011 */ 4012 public static float copySign(float magnitude, float sign){ 4013 int m = Float.floatToIntBits(magnitude); 4014 int s = Float.floatToIntBits(sign); 4015 if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK 4016 return magnitude; 4017 } 4018 return -magnitude; // flip sign 4019 } 4020 4021 /** 4022 * Return the exponent of a double number, removing the bias. 4023 * <p> 4024 * For double numbers of the form 2<sup>x</sup>, the unbiased 4025 * exponent is exactly x. 4026 * </p> 4027 * @param d number from which exponent is requested 4028 * @return exponent for d in IEEE754 representation, without bias 4029 */ 4030 public static int getExponent(final double d) { 4031 return (int) ((Double.doubleToLongBits(d) >>> 52) & 0x7ff) - 1023; 4032 } 4033 4034 /** 4035 * Return the exponent of a float number, removing the bias. 4036 * <p> 4037 * For float numbers of the form 2<sup>x</sup>, the unbiased 4038 * exponent is exactly x. 4039 * </p> 4040 * @param f number from which exponent is requested 4041 * @return exponent for d in IEEE754 representation, without bias 4042 */ 4043 public static int getExponent(final float f) { 4044 return ((Float.floatToIntBits(f) >>> 23) & 0xff) - 127; 4045 } 4046 4047} 4048