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