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