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