1 /* java.lang.StrictMath -- common mathematical functions, strict Java 2 Copyright (C) 1998, 2001, 2002, 2003, 2006 Free Software Foundation, Inc. 3 4 This file is part of GNU Classpath. 5 6 GNU Classpath is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation; either version 2, or (at your option) 9 any later version. 10 11 GNU Classpath is distributed in the hope that it will be useful, but 12 WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14 General Public License for more details. 15 16 You should have received a copy of the GNU General Public License 17 along with GNU Classpath; see the file COPYING. If not, write to the 18 Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 19 02110-1301 USA. 20 21 Linking this library statically or dynamically with other modules is 22 making a combined work based on this library. Thus, the terms and 23 conditions of the GNU General Public License cover the whole 24 combination. 25 26 As a special exception, the copyright holders of this library give you 27 permission to link this library with independent modules to produce an 28 executable, regardless of the license terms of these independent 29 modules, and to copy and distribute the resulting executable under 30 terms of your choice, provided that you also meet, for each linked 31 independent module, the terms and conditions of the license of that 32 module. An independent module is a module which is not derived from 33 or based on this library. If you modify this library, you may extend 34 this exception to your version of the library, but you are not 35 obligated to do so. If you do not wish to do so, delete this 36 exception statement from your version. */ 37 38 /* 39 * Some of the algorithms in this class are in the public domain, as part 40 * of fdlibm (freely-distributable math library), available at 41 * http://www.netlib.org/fdlibm/, and carry the following copyright: 42 * ==================================================== 43 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 44 * 45 * Developed at SunSoft, a Sun Microsystems, Inc. business. 46 * Permission to use, copy, modify, and distribute this 47 * software is freely granted, provided that this notice 48 * is preserved. 49 * ==================================================== 50 */ 51 52 package java.lang; 53 54 import gnu.classpath.Configuration; 55 56 import java.util.Random; 57 58 /** 59 * Helper class containing useful mathematical functions and constants. 60 * This class mirrors {@link Math}, but is 100% portable, because it uses 61 * no native methods whatsoever. Also, these algorithms are all accurate 62 * to less than 1 ulp, and execute in <code>strictfp</code> mode, while 63 * Math is allowed to vary in its results for some functions. Unfortunately, 64 * this usually means StrictMath has less efficiency and speed, as Math can 65 * use native methods. 66 * 67 * <p>The source of the various algorithms used is the fdlibm library, at:<br> 68 * <a href="http://www.netlib.org/fdlibm/">http://www.netlib.org/fdlibm/</a> 69 * 70 * Note that angles are specified in radians. Conversion functions are 71 * provided for your convenience. 72 * 73 * @author Eric Blake (ebb9@email.byu.edu) 74 * @since 1.3 75 */ 76 public final strictfp class StrictMath 77 { 78 /** 79 * StrictMath is non-instantiable. 80 */ StrictMath()81 private StrictMath() 82 { 83 } 84 85 /** 86 * A random number generator, initialized on first use. 87 * 88 * @see #random() 89 */ 90 private static Random rand; 91 92 /** 93 * The most accurate approximation to the mathematical constant <em>e</em>: 94 * <code>2.718281828459045</code>. Used in natural log and exp. 95 * 96 * @see #log(double) 97 * @see #exp(double) 98 */ 99 public static final double E 100 = 2.718281828459045; // Long bits 0x4005bf0z8b145769L. 101 102 /** 103 * The most accurate approximation to the mathematical constant <em>pi</em>: 104 * <code>3.141592653589793</code>. This is the ratio of a circle's diameter 105 * to its circumference. 106 */ 107 public static final double PI 108 = 3.141592653589793; // Long bits 0x400921fb54442d18L. 109 110 /** 111 * Take the absolute value of the argument. (Absolute value means make 112 * it positive.) 113 * 114 * <p>Note that the the largest negative value (Integer.MIN_VALUE) cannot 115 * be made positive. In this case, because of the rules of negation in 116 * a computer, MIN_VALUE is what will be returned. 117 * This is a <em>negative</em> value. You have been warned. 118 * 119 * @param i the number to take the absolute value of 120 * @return the absolute value 121 * @see Integer#MIN_VALUE 122 */ abs(int i)123 public static int abs(int i) 124 { 125 return (i < 0) ? -i : i; 126 } 127 128 /** 129 * Take the absolute value of the argument. (Absolute value means make 130 * it positive.) 131 * 132 * <p>Note that the the largest negative value (Long.MIN_VALUE) cannot 133 * be made positive. In this case, because of the rules of negation in 134 * a computer, MIN_VALUE is what will be returned. 135 * This is a <em>negative</em> value. You have been warned. 136 * 137 * @param l the number to take the absolute value of 138 * @return the absolute value 139 * @see Long#MIN_VALUE 140 */ abs(long l)141 public static long abs(long l) 142 { 143 return (l < 0) ? -l : l; 144 } 145 146 /** 147 * Take the absolute value of the argument. (Absolute value means make 148 * it positive.) 149 * 150 * @param f the number to take the absolute value of 151 * @return the absolute value 152 */ abs(float f)153 public static float abs(float f) 154 { 155 return (f <= 0) ? 0 - f : f; 156 } 157 158 /** 159 * Take the absolute value of the argument. (Absolute value means make 160 * it positive.) 161 * 162 * @param d the number to take the absolute value of 163 * @return the absolute value 164 */ abs(double d)165 public static double abs(double d) 166 { 167 return (d <= 0) ? 0 - d : d; 168 } 169 170 /** 171 * Return whichever argument is smaller. 172 * 173 * @param a the first number 174 * @param b a second number 175 * @return the smaller of the two numbers 176 */ min(int a, int b)177 public static int min(int a, int b) 178 { 179 return (a < b) ? a : b; 180 } 181 182 /** 183 * Return whichever argument is smaller. 184 * 185 * @param a the first number 186 * @param b a second number 187 * @return the smaller of the two numbers 188 */ min(long a, long b)189 public static long min(long a, long b) 190 { 191 return (a < b) ? a : b; 192 } 193 194 /** 195 * Return whichever argument is smaller. If either argument is NaN, the 196 * result is NaN, and when comparing 0 and -0, -0 is always smaller. 197 * 198 * @param a the first number 199 * @param b a second number 200 * @return the smaller of the two numbers 201 */ min(float a, float b)202 public static float min(float a, float b) 203 { 204 // this check for NaN, from JLS 15.21.1, saves a method call 205 if (a != a) 206 return a; 207 // no need to check if b is NaN; < will work correctly 208 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special 209 if (a == 0 && b == 0) 210 return -(-a - b); 211 return (a < b) ? a : b; 212 } 213 214 /** 215 * Return whichever argument is smaller. If either argument is NaN, the 216 * result is NaN, and when comparing 0 and -0, -0 is always smaller. 217 * 218 * @param a the first number 219 * @param b a second number 220 * @return the smaller of the two numbers 221 */ min(double a, double b)222 public static double min(double a, double b) 223 { 224 // this check for NaN, from JLS 15.21.1, saves a method call 225 if (a != a) 226 return a; 227 // no need to check if b is NaN; < will work correctly 228 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special 229 if (a == 0 && b == 0) 230 return -(-a - b); 231 return (a < b) ? a : b; 232 } 233 234 /** 235 * Return whichever argument is larger. 236 * 237 * @param a the first number 238 * @param b a second number 239 * @return the larger of the two numbers 240 */ max(int a, int b)241 public static int max(int a, int b) 242 { 243 return (a > b) ? a : b; 244 } 245 246 /** 247 * Return whichever argument is larger. 248 * 249 * @param a the first number 250 * @param b a second number 251 * @return the larger of the two numbers 252 */ max(long a, long b)253 public static long max(long a, long b) 254 { 255 return (a > b) ? a : b; 256 } 257 258 /** 259 * Return whichever argument is larger. If either argument is NaN, the 260 * result is NaN, and when comparing 0 and -0, 0 is always larger. 261 * 262 * @param a the first number 263 * @param b a second number 264 * @return the larger of the two numbers 265 */ max(float a, float b)266 public static float max(float a, float b) 267 { 268 // this check for NaN, from JLS 15.21.1, saves a method call 269 if (a != a) 270 return a; 271 // no need to check if b is NaN; > will work correctly 272 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special 273 if (a == 0 && b == 0) 274 return a - -b; 275 return (a > b) ? a : b; 276 } 277 278 /** 279 * Return whichever argument is larger. If either argument is NaN, the 280 * result is NaN, and when comparing 0 and -0, 0 is always larger. 281 * 282 * @param a the first number 283 * @param b a second number 284 * @return the larger of the two numbers 285 */ max(double a, double b)286 public static double max(double a, double b) 287 { 288 // this check for NaN, from JLS 15.21.1, saves a method call 289 if (a != a) 290 return a; 291 // no need to check if b is NaN; > will work correctly 292 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special 293 if (a == 0 && b == 0) 294 return a - -b; 295 return (a > b) ? a : b; 296 } 297 298 /** 299 * The trigonometric function <em>sin</em>. The sine of NaN or infinity is 300 * NaN, and the sine of 0 retains its sign. 301 * 302 * @param a the angle (in radians) 303 * @return sin(a) 304 */ sin(double a)305 public static double sin(double a) 306 { 307 if (a == Double.NEGATIVE_INFINITY || ! (a < Double.POSITIVE_INFINITY)) 308 return Double.NaN; 309 310 if (abs(a) <= PI / 4) 311 return sin(a, 0); 312 313 // Argument reduction needed. 314 double[] y = new double[2]; 315 int n = remPiOver2(a, y); 316 switch (n & 3) 317 { 318 case 0: 319 return sin(y[0], y[1]); 320 case 1: 321 return cos(y[0], y[1]); 322 case 2: 323 return -sin(y[0], y[1]); 324 default: 325 return -cos(y[0], y[1]); 326 } 327 } 328 329 /** 330 * The trigonometric function <em>cos</em>. The cosine of NaN or infinity is 331 * NaN. 332 * 333 * @param a the angle (in radians). 334 * @return cos(a). 335 */ cos(double a)336 public static double cos(double a) 337 { 338 if (a == Double.NEGATIVE_INFINITY || ! (a < Double.POSITIVE_INFINITY)) 339 return Double.NaN; 340 341 if (abs(a) <= PI / 4) 342 return cos(a, 0); 343 344 // Argument reduction needed. 345 double[] y = new double[2]; 346 int n = remPiOver2(a, y); 347 switch (n & 3) 348 { 349 case 0: 350 return cos(y[0], y[1]); 351 case 1: 352 return -sin(y[0], y[1]); 353 case 2: 354 return -cos(y[0], y[1]); 355 default: 356 return sin(y[0], y[1]); 357 } 358 } 359 360 /** 361 * The trigonometric function <em>tan</em>. The tangent of NaN or infinity 362 * is NaN, and the tangent of 0 retains its sign. 363 * 364 * @param a the angle (in radians) 365 * @return tan(a) 366 */ tan(double a)367 public static double tan(double a) 368 { 369 if (a == Double.NEGATIVE_INFINITY || ! (a < Double.POSITIVE_INFINITY)) 370 return Double.NaN; 371 372 if (abs(a) <= PI / 4) 373 return tan(a, 0, false); 374 375 // Argument reduction needed. 376 double[] y = new double[2]; 377 int n = remPiOver2(a, y); 378 return tan(y[0], y[1], (n & 1) == 1); 379 } 380 381 /** 382 * The trigonometric function <em>arcsin</em>. The range of angles returned 383 * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN or 384 * its absolute value is beyond 1, the result is NaN; and the arcsine of 385 * 0 retains its sign. 386 * 387 * @param x the sin to turn back into an angle 388 * @return arcsin(x) 389 */ asin(double x)390 public static double asin(double x) 391 { 392 boolean negative = x < 0; 393 if (negative) 394 x = -x; 395 if (! (x <= 1)) 396 return Double.NaN; 397 if (x == 1) 398 return negative ? -PI / 2 : PI / 2; 399 if (x < 0.5) 400 { 401 if (x < 1 / TWO_27) 402 return negative ? -x : x; 403 double t = x * x; 404 double p = t * (PS0 + t * (PS1 + t * (PS2 + t * (PS3 + t 405 * (PS4 + t * PS5))))); 406 double q = 1 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4))); 407 return negative ? -x - x * (p / q) : x + x * (p / q); 408 } 409 double w = 1 - x; // 1>|x|>=0.5. 410 double t = w * 0.5; 411 double p = t * (PS0 + t * (PS1 + t * (PS2 + t * (PS3 + t 412 * (PS4 + t * PS5))))); 413 double q = 1 + t * (QS1 + t * (QS2 + t * (QS3 + t * QS4))); 414 double s = sqrt(t); 415 if (x >= 0.975) 416 { 417 w = p / q; 418 t = PI / 2 - (2 * (s + s * w) - PI_L / 2); 419 } 420 else 421 { 422 w = (float) s; 423 double c = (t - w * w) / (s + w); 424 p = 2 * s * (p / q) - (PI_L / 2 - 2 * c); 425 q = PI / 4 - 2 * w; 426 t = PI / 4 - (p - q); 427 } 428 return negative ? -t : t; 429 } 430 431 /** 432 * The trigonometric function <em>arccos</em>. The range of angles returned 433 * is 0 to pi radians (0 to 180 degrees). If the argument is NaN or 434 * its absolute value is beyond 1, the result is NaN. 435 * 436 * @param x the cos to turn back into an angle 437 * @return arccos(x) 438 */ 439 public static double acos(double x) 440 { 441 boolean negative = x < 0; 442 if (negative) 443 x = -x; 444 if (! (x <= 1)) 445 return Double.NaN; 446 if (x == 1) 447 return negative ? PI : 0; 448 if (x < 0.5) 449 { 450 if (x < 1 / TWO_57) 451 return PI / 2; 452 double z = x * x; 453 double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z 454 * (PS4 + z * PS5))))); 455 double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); 456 double r = x - (PI_L / 2 - x * (p / q)); 457 return negative ? PI / 2 + r : PI / 2 - r; 458 } 459 if (negative) // x<=-0.5. 460 { 461 double z = (1 + x) * 0.5; 462 double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z 463 * (PS4 + z * PS5))))); 464 double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); 465 double s = sqrt(z); 466 double w = p / q * s - PI_L / 2; 467 return PI - 2 * (s + w); 468 } 469 double z = (1 - x) * 0.5; // x>0.5. 470 double s = sqrt(z); 471 double df = (float) s; 472 double c = (z - df * df) / (s + df); 473 double p = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z 474 * (PS4 + z * PS5))))); 475 double q = 1 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4))); 476 double w = p / q * s + c; 477 return 2 * (df + w); 478 } 479 480 /** 481 * The trigonometric function <em>arcsin</em>. The range of angles returned 482 * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN, the 483 * result is NaN; and the arctangent of 0 retains its sign. 484 * 485 * @param x the tan to turn back into an angle 486 * @return arcsin(x) 487 * @see #atan2(double, double) 488 */ 489 public static double atan(double x) 490 { 491 double lo; 492 double hi; 493 boolean negative = x < 0; 494 if (negative) 495 x = -x; 496 if (x >= TWO_66) 497 return negative ? -PI / 2 : PI / 2; 498 if (! (x >= 0.4375)) // |x|<7/16, or NaN. 499 { 500 if (! (x >= 1 / TWO_29)) // Small, or NaN. 501 return negative ? -x : x; 502 lo = hi = 0; 503 } 504 else if (x < 1.1875) 505 { 506 if (x < 0.6875) // 7/16<=|x|<11/16. 507 { 508 x = (2 * x - 1) / (2 + x); 509 hi = ATAN_0_5H; 510 lo = ATAN_0_5L; 511 } 512 else // 11/16<=|x|<19/16. 513 { 514 x = (x - 1) / (x + 1); 515 hi = PI / 4; 516 lo = PI_L / 4; 517 } 518 } 519 else if (x < 2.4375) // 19/16<=|x|<39/16. 520 { 521 x = (x - 1.5) / (1 + 1.5 * x); 522 hi = ATAN_1_5H; 523 lo = ATAN_1_5L; 524 } 525 else // 39/16<=|x|<2**66. 526 { 527 x = -1 / x; 528 hi = PI / 2; 529 lo = PI_L / 2; 530 } 531 532 // Break sum from i=0 to 10 ATi*z**(i+1) into odd and even poly. 533 double z = x * x; 534 double w = z * z; 535 double s1 = z * (AT0 + w * (AT2 + w * (AT4 + w * (AT6 + w 536 * (AT8 + w * AT10))))); 537 double s2 = w * (AT1 + w * (AT3 + w * (AT5 + w * (AT7 + w * AT9)))); 538 if (hi == 0) 539 return negative ? x * (s1 + s2) - x : x - x * (s1 + s2); 540 z = hi - ((x * (s1 + s2) - lo) - x); 541 return negative ? -z : z; 542 } 543 544 /** 545 * A special version of the trigonometric function <em>arctan</em>, for 546 * converting rectangular coordinates <em>(x, y)</em> to polar 547 * <em>(r, theta)</em>. This computes the arctangent of x/y in the range 548 * of -pi to pi radians (-180 to 180 degrees). Special cases:<ul> 549 * <li>If either argument is NaN, the result is NaN.</li> 550 * <li>If the first argument is positive zero and the second argument is 551 * positive, or the first argument is positive and finite and the second 552 * argument is positive infinity, then the result is positive zero.</li> 553 * <li>If the first argument is negative zero and the second argument is 554 * positive, or the first argument is negative and finite and the second 555 * argument is positive infinity, then the result is negative zero.</li> 556 * <li>If the first argument is positive zero and the second argument is 557 * negative, or the first argument is positive and finite and the second 558 * argument is negative infinity, then the result is the double value 559 * closest to pi.</li> 560 * <li>If the first argument is negative zero and the second argument is 561 * negative, or the first argument is negative and finite and the second 562 * argument is negative infinity, then the result is the double value 563 * closest to -pi.</li> 564 * <li>If the first argument is positive and the second argument is 565 * positive zero or negative zero, or the first argument is positive 566 * infinity and the second argument is finite, then the result is the 567 * double value closest to pi/2.</li> 568 * <li>If the first argument is negative and the second argument is 569 * positive zero or negative zero, or the first argument is negative 570 * infinity and the second argument is finite, then the result is the 571 * double value closest to -pi/2.</li> 572 * <li>If both arguments are positive infinity, then the result is the 573 * double value closest to pi/4.</li> 574 * <li>If the first argument is positive infinity and the second argument 575 * is negative infinity, then the result is the double value closest to 576 * 3*pi/4.</li> 577 * <li>If the first argument is negative infinity and the second argument 578 * is positive infinity, then the result is the double value closest to 579 * -pi/4.</li> 580 * <li>If both arguments are negative infinity, then the result is the 581 * double value closest to -3*pi/4.</li> 582 * 583 * </ul><p>This returns theta, the angle of the point. To get r, albeit 584 * slightly inaccurately, use sqrt(x*x+y*y). 585 * 586 * @param y the y position 587 * @param x the x position 588 * @return <em>theta</em> in the conversion of (x, y) to (r, theta) 589 * @see #atan(double) 590 */ 591 public static double atan2(double y, double x) 592 { 593 if (x != x || y != y) 594 return Double.NaN; 595 if (x == 1) 596 return atan(y); 597 if (x == Double.POSITIVE_INFINITY) 598 { 599 if (y == Double.POSITIVE_INFINITY) 600 return PI / 4; 601 if (y == Double.NEGATIVE_INFINITY) 602 return -PI / 4; 603 return 0 * y; 604 } 605 if (x == Double.NEGATIVE_INFINITY) 606 { 607 if (y == Double.POSITIVE_INFINITY) 608 return 3 * PI / 4; 609 if (y == Double.NEGATIVE_INFINITY) 610 return -3 * PI / 4; 611 return (1 / (0 * y) == Double.POSITIVE_INFINITY) ? PI : -PI; 612 } 613 if (y == 0) 614 { 615 if (1 / (0 * x) == Double.POSITIVE_INFINITY) 616 return y; 617 return (1 / y == Double.POSITIVE_INFINITY) ? PI : -PI; 618 } 619 if (y == Double.POSITIVE_INFINITY || y == Double.NEGATIVE_INFINITY 620 || x == 0) 621 return y < 0 ? -PI / 2 : PI / 2; 622 623 double z = abs(y / x); // Safe to do y/x. 624 if (z > TWO_60) 625 z = PI / 2 + 0.5 * PI_L; 626 else if (x < 0 && z < 1 / TWO_60) 627 z = 0; 628 else 629 z = atan(z); 630 if (x > 0) 631 return y > 0 ? z : -z; 632 return y > 0 ? PI - (z - PI_L) : z - PI_L - PI; 633 } 634 635 /** 636 * Returns the hyperbolic sine of <code>x</code> which is defined as 637 * (exp(x) - exp(-x)) / 2. 638 * 639 * Special cases: 640 * <ul> 641 * <li>If the argument is NaN, the result is NaN</li> 642 * <li>If the argument is positive infinity, the result is positive 643 * infinity.</li> 644 * <li>If the argument is negative infinity, the result is negative 645 * infinity.</li> 646 * <li>If the argument is zero, the result is zero.</li> 647 * </ul> 648 * 649 * @param x the argument to <em>sinh</em> 650 * @return the hyperbolic sine of <code>x</code> 651 * 652 * @since 1.5 653 */ 654 public static double sinh(double x) 655 { 656 // Method : 657 // mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 658 // 1. Replace x by |x| (sinh(-x) = -sinh(x)). 659 // 2. 660 // E + E/(E+1) 661 // 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) 662 // 2 663 // 664 // 22 <= x <= lnovft : sinh(x) := exp(x)/2 665 // lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) 666 // ln2ovft < x : sinh(x) := +inf (overflow) 667 668 double t, w, h; 669 670 long bits; 671 long h_bits; 672 long l_bits; 673 674 // handle special cases 675 if (x != x) 676 return x; 677 if (x == Double.POSITIVE_INFINITY) 678 return Double.POSITIVE_INFINITY; 679 if (x == Double.NEGATIVE_INFINITY) 680 return Double.NEGATIVE_INFINITY; 681 682 if (x < 0) 683 h = - 0.5; 684 else 685 h = 0.5; 686 687 bits = Double.doubleToLongBits(x); 688 h_bits = getHighDWord(bits) & 0x7fffffffL; // ignore sign 689 l_bits = getLowDWord(bits); 690 691 // |x| in [0, 22], return sign(x) * 0.5 * (E+E/(E+1)) 692 if (h_bits < 0x40360000L) // |x| < 22 693 { 694 if (h_bits < 0x3e300000L) // |x| < 2^-28 695 return x; // for tiny arguments return x 696 697 t = expm1(abs(x)); 698 699 if (h_bits < 0x3ff00000L) 700 return h * (2.0 * t - t * t / (t + 1.0)); 701 702 return h * (t + t / (t + 1.0)); 703 } 704 705 // |x| in [22, log(Double.MAX_VALUE)], return 0.5 * exp(|x|) 706 if (h_bits < 0x40862e42L) 707 return h * exp(abs(x)); 708 709 // |x| in [log(Double.MAX_VALUE), overflowthreshold] 710 if ((h_bits < 0x408633ceL) 711 || ((h_bits == 0x408633ceL) && (l_bits <= 0x8fb9f87dL))) 712 { 713 w = exp(0.5 * abs(x)); 714 t = h * w; 715 716 return t * w; 717 } 718 719 // |x| > overflowthershold 720 return h * Double.POSITIVE_INFINITY; 721 } 722 723 /** 724 * Returns the hyperbolic cosine of <code>x</code>, which is defined as 725 * (exp(x) + exp(-x)) / 2. 726 * 727 * Special cases: 728 * <ul> 729 * <li>If the argument is NaN, the result is NaN</li> 730 * <li>If the argument is positive infinity, the result is positive 731 * infinity.</li> 732 * <li>If the argument is negative infinity, the result is positive 733 * infinity.</li> 734 * <li>If the argument is zero, the result is one.</li> 735 * </ul> 736 * 737 * @param x the argument to <em>cosh</em> 738 * @return the hyperbolic cosine of <code>x</code> 739 * 740 * @since 1.5 741 */ 742 public static double cosh(double x) 743 { 744 // Method : 745 // mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 746 // 1. Replace x by |x| (cosh(x) = cosh(-x)). 747 // 2. 748 // [ exp(x) - 1 ]^2 749 // 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- 750 // 2*exp(x) 751 // 752 // exp(x) + 1/exp(x) 753 // ln2/2 <= x <= 22 : cosh(x) := ------------------ 754 // 2 755 // 22 <= x <= lnovft : cosh(x) := exp(x)/2 756 // lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) 757 // ln2ovft < x : cosh(x) := +inf (overflow) 758 759 double t, w; 760 long bits; 761 long hx; 762 long lx; 763 764 // handle special cases 765 if (x != x) 766 return x; 767 if (x == Double.POSITIVE_INFINITY) 768 return Double.POSITIVE_INFINITY; 769 if (x == Double.NEGATIVE_INFINITY) 770 return Double.POSITIVE_INFINITY; 771 772 bits = Double.doubleToLongBits(x); 773 hx = getHighDWord(bits) & 0x7fffffffL; // ignore sign 774 lx = getLowDWord(bits); 775 776 // |x| in [0, 0.5 * ln(2)], return 1 + expm1(|x|)^2 / (2 * exp(|x|)) 777 if (hx < 0x3fd62e43L) 778 { 779 t = expm1(abs(x)); 780 w = 1.0 + t; 781 782 // for tiny arguments return 1. 783 if (hx < 0x3c800000L) 784 return w; 785 786 return 1.0 + (t * t) / (w + w); 787 } 788 789 // |x| in [0.5 * ln(2), 22], return exp(|x|)/2 + 1 / (2 * exp(|x|)) 790 if (hx < 0x40360000L) 791 { 792 t = exp(abs(x)); 793 794 return 0.5 * t + 0.5 / t; 795 } 796 797 // |x| in [22, log(Double.MAX_VALUE)], return 0.5 * exp(|x|) 798 if (hx < 0x40862e42L) 799 return 0.5 * exp(abs(x)); 800 801 // |x| in [log(Double.MAX_VALUE), overflowthreshold], 802 // return exp(x/2)/2 * exp(x/2) 803 if ((hx < 0x408633ceL) 804 || ((hx == 0x408633ceL) && (lx <= 0x8fb9f87dL))) 805 { 806 w = exp(0.5 * abs(x)); 807 t = 0.5 * w; 808 809 return t * w; 810 } 811 812 // |x| > overflowthreshold 813 return Double.POSITIVE_INFINITY; 814 } 815 816 /** 817 * Returns the hyperbolic tangent of <code>x</code>, which is defined as 818 * (exp(x) - exp(-x)) / (exp(x) + exp(-x)), i.e. sinh(x) / cosh(x). 819 * 820 Special cases: 821 * <ul> 822 * <li>If the argument is NaN, the result is NaN</li> 823 * <li>If the argument is positive infinity, the result is 1.</li> 824 * <li>If the argument is negative infinity, the result is -1.</li> 825 * <li>If the argument is zero, the result is zero.</li> 826 * </ul> 827 * 828 * @param x the argument to <em>tanh</em> 829 * @return the hyperbolic tagent of <code>x</code> 830 * 831 * @since 1.5 832 */ 833 public static double tanh(double x) 834 { 835 // Method : 836 // 0. tanh(x) is defined to be (exp(x) - exp(-x)) / (exp(x) + exp(-x)) 837 // 1. reduce x to non-negative by tanh(-x) = -tanh(x). 838 // 2. 0 <= x <= 2^-55 : tanh(x) := x * (1.0 + x) 839 // -t 840 // 2^-55 < x <= 1 : tanh(x) := -----; t = expm1(-2x) 841 // t + 2 842 // 2 843 // 1 <= x <= 22.0 : tanh(x) := 1 - ----- ; t=expm1(2x) 844 // t + 2 845 // 22.0 < x <= INF : tanh(x) := 1. 846 847 double t, z; 848 849 long bits; 850 long h_bits; 851 852 // handle special cases 853 if (x != x) 854 return x; 855 if (x == Double.POSITIVE_INFINITY) 856 return 1.0; 857 if (x == Double.NEGATIVE_INFINITY) 858 return -1.0; 859 860 bits = Double.doubleToLongBits(x); 861 h_bits = getHighDWord(bits) & 0x7fffffffL; // ingnore sign 862 863 if (h_bits < 0x40360000L) // |x| < 22 864 { 865 if (h_bits < 0x3c800000L) // |x| < 2^-55 866 return x * (1.0 + x); 867 868 if (h_bits >= 0x3ff00000L) // |x| >= 1 869 { 870 t = expm1(2.0 * abs(x)); 871 z = 1.0 - 2.0 / (t + 2.0); 872 } 873 else // |x| < 1 874 { 875 t = expm1(-2.0 * abs(x)); 876 z = -t / (t + 2.0); 877 } 878 } 879 else // |x| >= 22 880 z = 1.0; 881 882 return (x >= 0) ? z : -z; 883 } 884 885 /** 886 * Returns the lower two words of a long. This is intended to be 887 * used like this: 888 * <code>getLowDWord(Double.doubleToLongBits(x))</code>. 889 */ 890 private static long getLowDWord(long x) 891 { 892 return x & 0x00000000ffffffffL; 893 } 894 895 /** 896 * Returns the higher two words of a long. This is intended to be 897 * used like this: 898 * <code>getHighDWord(Double.doubleToLongBits(x))</code>. 899 */ 900 private static long getHighDWord(long x) 901 { 902 return (x & 0xffffffff00000000L) >> 32; 903 } 904 905 /** 906 * Returns a double with the IEEE754 bit pattern given in the lower 907 * and higher two words <code>lowDWord</code> and <code>highDWord</code>. 908 */ 909 private static double buildDouble(long lowDWord, long highDWord) 910 { 911 return Double.longBitsToDouble(((highDWord & 0xffffffffL) << 32) 912 | (lowDWord & 0xffffffffL)); 913 } 914 915 /** 916 * Returns the cube root of <code>x</code>. The sign of the cube root 917 * is equal to the sign of <code>x</code>. 918 * 919 * Special cases: 920 * <ul> 921 * <li>If the argument is NaN, the result is NaN</li> 922 * <li>If the argument is positive infinity, the result is positive 923 * infinity.</li> 924 * <li>If the argument is negative infinity, the result is negative 925 * infinity.</li> 926 * <li>If the argument is zero, the result is zero with the same 927 * sign as the argument.</li> 928 * </ul> 929 * 930 * @param x the number to take the cube root of 931 * @return the cube root of <code>x</code> 932 * @see #sqrt(double) 933 * 934 * @since 1.5 935 */ 936 public static double cbrt(double x) 937 { 938 boolean negative = (x < 0); 939 double r; 940 double s; 941 double t; 942 double w; 943 944 long bits; 945 long l; 946 long h; 947 948 // handle the special cases 949 if (x != x) 950 return x; 951 if (x == Double.POSITIVE_INFINITY) 952 return Double.POSITIVE_INFINITY; 953 if (x == Double.NEGATIVE_INFINITY) 954 return Double.NEGATIVE_INFINITY; 955 if (x == 0) 956 return x; 957 958 x = abs(x); 959 bits = Double.doubleToLongBits(x); 960 961 if (bits < 0x0010000000000000L) // subnormal number 962 { 963 t = TWO_54; 964 t *= x; 965 966 // __HI(t)=__HI(t)/3+B2; 967 bits = Double.doubleToLongBits(t); 968 h = getHighDWord(bits); 969 l = getLowDWord(bits); 970 971 h = h / 3 + CBRT_B2; 972 973 t = buildDouble(l, h); 974 } 975 else 976 { 977 // __HI(t)=__HI(x)/3+B1; 978 h = getHighDWord(bits); 979 l = 0; 980 981 h = h / 3 + CBRT_B1; 982 t = buildDouble(l, h); 983 } 984 985 // new cbrt to 23 bits 986 r = t * t / x; 987 s = CBRT_C + r * t; 988 t *= CBRT_G + CBRT_F / (s + CBRT_E + CBRT_D / s); 989 990 // chopped to 20 bits and make it larger than cbrt(x) 991 bits = Double.doubleToLongBits(t); 992 h = getHighDWord(bits); 993 994 // __LO(t)=0; 995 // __HI(t)+=0x00000001; 996 l = 0; 997 h += 1; 998 t = buildDouble(l, h); 999 1000 // one step newton iteration to 53 bits with error less than 0.667 ulps 1001 s = t * t; // t * t is exact 1002 r = x / s; 1003 w = t + t; 1004 r = (r - t) / (w + r); // r - t is exact 1005 t = t + t * r; 1006 1007 return negative ? -t : t; 1008 } 1009 1010 /** 1011 * Take <em>e</em><sup>a</sup>. The opposite of <code>log()</code>. If the 1012 * argument is NaN, the result is NaN; if the argument is positive infinity, 1013 * the result is positive infinity; and if the argument is negative 1014 * infinity, the result is positive zero. 1015 * 1016 * @param x the number to raise to the power 1017 * @return the number raised to the power of <em>e</em> 1018 * @see #log(double) 1019 * @see #pow(double, double) 1020 */ 1021 public static double exp(double x) 1022 { 1023 if (x != x) 1024 return x; 1025 if (x > EXP_LIMIT_H) 1026 return Double.POSITIVE_INFINITY; 1027 if (x < EXP_LIMIT_L) 1028 return 0; 1029 1030 // Argument reduction. 1031 double hi; 1032 double lo; 1033 int k; 1034 double t = abs(x); 1035 if (t > 0.5 * LN2) 1036 { 1037 if (t < 1.5 * LN2) 1038 { 1039 hi = t - LN2_H; 1040 lo = LN2_L; 1041 k = 1; 1042 } 1043 else 1044 { 1045 k = (int) (INV_LN2 * t + 0.5); 1046 hi = t - k * LN2_H; 1047 lo = k * LN2_L; 1048 } 1049 if (x < 0) 1050 { 1051 hi = -hi; 1052 lo = -lo; 1053 k = -k; 1054 } 1055 x = hi - lo; 1056 } 1057 else if (t < 1 / TWO_28) 1058 return 1; 1059 else 1060 lo = hi = k = 0; 1061 1062 // Now x is in primary range. 1063 t = x * x; 1064 double c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); 1065 if (k == 0) 1066 return 1 - (x * c / (c - 2) - x); 1067 double y = 1 - (lo - x * c / (2 - c) - hi); 1068 return scale(y, k); 1069 } 1070 1071 /** 1072 * Returns <em>e</em><sup>x</sup> - 1. 1073 * Special cases: 1074 * <ul> 1075 * <li>If the argument is NaN, the result is NaN.</li> 1076 * <li>If the argument is positive infinity, the result is positive 1077 * infinity</li> 1078 * <li>If the argument is negative infinity, the result is -1.</li> 1079 * <li>If the argument is zero, the result is zero.</li> 1080 * </ul> 1081 * 1082 * @param x the argument to <em>e</em><sup>x</sup> - 1. 1083 * @return <em>e</em> raised to the power <code>x</code> minus one. 1084 * @see #exp(double) 1085 */ 1086 public static double expm1(double x) 1087 { 1088 // Method 1089 // 1. Argument reduction: 1090 // Given x, find r and integer k such that 1091 // 1092 // x = k * ln(2) + r, |r| <= 0.5 * ln(2) 1093 // 1094 // Here a correction term c will be computed to compensate 1095 // the error in r when rounded to a floating-point number. 1096 // 1097 // 2. Approximating expm1(r) by a special rational function on 1098 // the interval [0, 0.5 * ln(2)]: 1099 // Since 1100 // r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 - r^4/360 + ... 1101 // we define R1(r*r) by 1102 // r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 * R1(r*r) 1103 // That is, 1104 // R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) 1105 // = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) 1106 // = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... 1107 // We use a special Remes algorithm on [0, 0.347] to generate 1108 // a polynomial of degree 5 in r*r to approximate R1. The 1109 // maximum error of this polynomial approximation is bounded 1110 // by 2**-61. In other words, 1111 // R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 1112 // where Q1 = -1.6666666666666567384E-2, 1113 // Q2 = 3.9682539681370365873E-4, 1114 // Q3 = -9.9206344733435987357E-6, 1115 // Q4 = 2.5051361420808517002E-7, 1116 // Q5 = -6.2843505682382617102E-9; 1117 // (where z=r*r, and Q1 to Q5 are called EXPM1_Qx in the source) 1118 // with error bounded by 1119 // | 5 | -61 1120 // | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 1121 // | | 1122 // 1123 // expm1(r) = exp(r)-1 is then computed by the following 1124 // specific way which minimize the accumulation rounding error: 1125 // 2 3 1126 // r r [ 3 - (R1 + R1*r/2) ] 1127 // expm1(r) = r + --- + --- * [--------------------] 1128 // 2 2 [ 6 - r*(3 - R1*r/2) ] 1129 // 1130 // To compensate the error in the argument reduction, we use 1131 // expm1(r+c) = expm1(r) + c + expm1(r)*c 1132 // ~ expm1(r) + c + r*c 1133 // Thus c+r*c will be added in as the correction terms for 1134 // expm1(r+c). Now rearrange the term to avoid optimization 1135 // screw up: 1136 // ( 2 2 ) 1137 // ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) 1138 // expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) 1139 // ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) 1140 // ( ) 1141 // 1142 // = r - E 1143 // 3. Scale back to obtain expm1(x): 1144 // From step 1, we have 1145 // expm1(x) = either 2^k*[expm1(r)+1] - 1 1146 // = or 2^k*[expm1(r) + (1-2^-k)] 1147 // 4. Implementation notes: 1148 // (A). To save one multiplication, we scale the coefficient Qi 1149 // to Qi*2^i, and replace z by (x^2)/2. 1150 // (B). To achieve maximum accuracy, we compute expm1(x) by 1151 // (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf) 1152 // (ii) if k=0, return r-E 1153 // (iii) if k=-1, return 0.5*(r-E)-0.5 1154 // (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) 1155 // else return 1.0+2.0*(r-E); 1156 // (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) 1157 // (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else 1158 // (vii) return 2^k(1-((E+2^-k)-r)) 1159 1160 boolean negative = (x < 0); 1161 double y, hi, lo, c, t, e, hxs, hfx, r1; 1162 int k; 1163 1164 long bits; 1165 long h_bits; 1166 long l_bits; 1167 1168 c = 0.0; 1169 y = abs(x); 1170 1171 bits = Double.doubleToLongBits(y); 1172 h_bits = getHighDWord(bits); 1173 l_bits = getLowDWord(bits); 1174 1175 // handle special cases and large arguments 1176 if (h_bits >= 0x4043687aL) // if |x| >= 56 * ln(2) 1177 { 1178 if (h_bits >= 0x40862e42L) // if |x| >= EXP_LIMIT_H 1179 { 1180 if (h_bits >= 0x7ff00000L) 1181 { 1182 if (((h_bits & 0x000fffffL) | (l_bits & 0xffffffffL)) != 0) 1183 return x; // exp(NaN) = NaN 1184 else 1185 return negative ? -1.0 : x; // exp({+-inf}) = {+inf, -1} 1186 } 1187 1188 if (x > EXP_LIMIT_H) 1189 return Double.POSITIVE_INFINITY; // overflow 1190 } 1191 1192 if (negative) // x <= -56 * ln(2) 1193 return -1.0; 1194 } 1195 1196 // argument reduction 1197 if (h_bits > 0x3fd62e42L) // |x| > 0.5 * ln(2) 1198 { 1199 if (h_bits < 0x3ff0a2b2L) // |x| < 1.5 * ln(2) 1200 { 1201 if (negative) 1202 { 1203 hi = x + LN2_H; 1204 lo = -LN2_L; 1205 k = -1; 1206 } 1207 else 1208 { 1209 hi = x - LN2_H; 1210 lo = LN2_L; 1211 k = 1; 1212 } 1213 } 1214 else 1215 { 1216 k = (int) (INV_LN2 * x + (negative ? - 0.5 : 0.5)); 1217 t = k; 1218 hi = x - t * LN2_H; 1219 lo = t * LN2_L; 1220 } 1221 1222 x = hi - lo; 1223 c = (hi - x) - lo; 1224 1225 } 1226 else if (h_bits < 0x3c900000L) // |x| < 2^-54 return x 1227 return x; 1228 else 1229 k = 0; 1230 1231 // x is now in primary range 1232 hfx = 0.5 * x; 1233 hxs = x * hfx; 1234 r1 = 1.0 + hxs * (EXPM1_Q1 1235 + hxs * (EXPM1_Q2 1236 + hxs * (EXPM1_Q3 1237 + hxs * (EXPM1_Q4 1238 + hxs * EXPM1_Q5)))); 1239 t = 3.0 - r1 * hfx; 1240 e = hxs * ((r1 - t) / (6.0 - x * t)); 1241 1242 if (k == 0) 1243 { 1244 return x - (x * e - hxs); // c == 0 1245 } 1246 else 1247 { 1248 e = x * (e - c) - c; 1249 e -= hxs; 1250 1251 if (k == -1) 1252 return 0.5 * (x - e) - 0.5; 1253 1254 if (k == 1) 1255 { 1256 if (x < - 0.25) 1257 return -2.0 * (e - (x + 0.5)); 1258 else 1259 return 1.0 + 2.0 * (x - e); 1260 } 1261 1262 if (k <= -2 || k > 56) // sufficient to return exp(x) - 1 1263 { 1264 y = 1.0 - (e - x); 1265 1266 bits = Double.doubleToLongBits(y); 1267 h_bits = getHighDWord(bits); 1268 l_bits = getLowDWord(bits); 1269 1270 h_bits += (k << 20); // add k to y's exponent 1271 1272 y = buildDouble(l_bits, h_bits); 1273 1274 return y - 1.0; 1275 } 1276 1277 t = 1.0; 1278 if (k < 20) 1279 { 1280 bits = Double.doubleToLongBits(t); 1281 h_bits = 0x3ff00000L - (0x00200000L >> k); 1282 l_bits = getLowDWord(bits); 1283 1284 t = buildDouble(l_bits, h_bits); // t = 1 - 2^(-k) 1285 y = t - (e - x); 1286 1287 bits = Double.doubleToLongBits(y); 1288 h_bits = getHighDWord(bits); 1289 l_bits = getLowDWord(bits); 1290 1291 h_bits += (k << 20); // add k to y's exponent 1292 1293 y = buildDouble(l_bits, h_bits); 1294 } 1295 else 1296 { 1297 bits = Double.doubleToLongBits(t); 1298 h_bits = (0x000003ffL - k) << 20; 1299 l_bits = getLowDWord(bits); 1300 1301 t = buildDouble(l_bits, h_bits); // t = 2^(-k) 1302 1303 y = x - (e + t); 1304 y += 1.0; 1305 1306 bits = Double.doubleToLongBits(y); 1307 h_bits = getHighDWord(bits); 1308 l_bits = getLowDWord(bits); 1309 1310 h_bits += (k << 20); // add k to y's exponent 1311 1312 y = buildDouble(l_bits, h_bits); 1313 } 1314 } 1315 1316 return y; 1317 } 1318 1319 /** 1320 * Take ln(a) (the natural log). The opposite of <code>exp()</code>. If the 1321 * argument is NaN or negative, the result is NaN; if the argument is 1322 * positive infinity, the result is positive infinity; and if the argument 1323 * is either zero, the result is negative infinity. 1324 * 1325 * <p>Note that the way to get log<sub>b</sub>(a) is to do this: 1326 * <code>ln(a) / ln(b)</code>. 1327 * 1328 * @param x the number to take the natural log of 1329 * @return the natural log of <code>a</code> 1330 * @see #exp(double) 1331 */ 1332 public static double log(double x) 1333 { 1334 if (x == 0) 1335 return Double.NEGATIVE_INFINITY; 1336 if (x < 0) 1337 return Double.NaN; 1338 if (! (x < Double.POSITIVE_INFINITY)) 1339 return x; 1340 1341 // Normalize x. 1342 long bits = Double.doubleToLongBits(x); 1343 int exp = (int) (bits >> 52); 1344 if (exp == 0) // Subnormal x. 1345 { 1346 x *= TWO_54; 1347 bits = Double.doubleToLongBits(x); 1348 exp = (int) (bits >> 52) - 54; 1349 } 1350 exp -= 1023; // Unbias exponent. 1351 bits = (bits & 0x000fffffffffffffL) | 0x3ff0000000000000L; 1352 x = Double.longBitsToDouble(bits); 1353 if (x >= SQRT_2) 1354 { 1355 x *= 0.5; 1356 exp++; 1357 } 1358 x--; 1359 if (abs(x) < 1 / TWO_20) 1360 { 1361 if (x == 0) 1362 return exp * LN2_H + exp * LN2_L; 1363 double r = x * x * (0.5 - 1 / 3.0 * x); 1364 if (exp == 0) 1365 return x - r; 1366 return exp * LN2_H - ((r - exp * LN2_L) - x); 1367 } 1368 double s = x / (2 + x); 1369 double z = s * s; 1370 double w = z * z; 1371 double t1 = w * (LG2 + w * (LG4 + w * LG6)); 1372 double t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7))); 1373 double r = t2 + t1; 1374 if (bits >= 0x3ff6174a00000000L && bits < 0x3ff6b85200000000L) 1375 { 1376 double h = 0.5 * x * x; // Need more accuracy for x near sqrt(2). 1377 if (exp == 0) 1378 return x - (h - s * (h + r)); 1379 return exp * LN2_H - ((h - (s * (h + r) + exp * LN2_L)) - x); 1380 } 1381 if (exp == 0) 1382 return x - s * (x - r); 1383 return exp * LN2_H - ((s * (x - r) - exp * LN2_L) - x); 1384 } 1385 1386 /** 1387 * Take a square root. If the argument is NaN or negative, the result is 1388 * NaN; if the argument is positive infinity, the result is positive 1389 * infinity; and if the result is either zero, the result is the same. 1390 * 1391 * <p>For other roots, use pow(x, 1/rootNumber). 1392 * 1393 * @param x the numeric argument 1394 * @return the square root of the argument 1395 * @see #pow(double, double) 1396 */ 1397 public static double sqrt(double x) 1398 { 1399 if (x < 0) 1400 return Double.NaN; 1401 if (x == 0 || ! (x < Double.POSITIVE_INFINITY)) 1402 return x; 1403 1404 // Normalize x. 1405 long bits = Double.doubleToLongBits(x); 1406 int exp = (int) (bits >> 52); 1407 if (exp == 0) // Subnormal x. 1408 { 1409 x *= TWO_54; 1410 bits = Double.doubleToLongBits(x); 1411 exp = (int) (bits >> 52) - 54; 1412 } 1413 exp -= 1023; // Unbias exponent. 1414 bits = (bits & 0x000fffffffffffffL) | 0x0010000000000000L; 1415 if ((exp & 1) == 1) // Odd exp, double x to make it even. 1416 bits <<= 1; 1417 exp >>= 1; 1418 1419 // Generate sqrt(x) bit by bit. 1420 bits <<= 1; 1421 long q = 0; 1422 long s = 0; 1423 long r = 0x0020000000000000L; // Move r right to left. 1424 while (r != 0) 1425 { 1426 long t = s + r; 1427 if (t <= bits) 1428 { 1429 s = t + r; 1430 bits -= t; 1431 q += r; 1432 } 1433 bits <<= 1; 1434 r >>= 1; 1435 } 1436 1437 // Use floating add to round correctly. 1438 if (bits != 0) 1439 q += q & 1; 1440 return Double.longBitsToDouble((q >> 1) + ((exp + 1022L) << 52)); 1441 } 1442 1443 /** 1444 * Raise a number to a power. Special cases:<ul> 1445 * <li>If the second argument is positive or negative zero, then the result 1446 * is 1.0.</li> 1447 * <li>If the second argument is 1.0, then the result is the same as the 1448 * first argument.</li> 1449 * <li>If the second argument is NaN, then the result is NaN.</li> 1450 * <li>If the first argument is NaN and the second argument is nonzero, 1451 * then the result is NaN.</li> 1452 * <li>If the absolute value of the first argument is greater than 1 and 1453 * the second argument is positive infinity, or the absolute value of the 1454 * first argument is less than 1 and the second argument is negative 1455 * infinity, then the result is positive infinity.</li> 1456 * <li>If the absolute value of the first argument is greater than 1 and 1457 * the second argument is negative infinity, or the absolute value of the 1458 * first argument is less than 1 and the second argument is positive 1459 * infinity, then the result is positive zero.</li> 1460 * <li>If the absolute value of the first argument equals 1 and the second 1461 * argument is infinite, then the result is NaN.</li> 1462 * <li>If the first argument is positive zero and the second argument is 1463 * greater than zero, or the first argument is positive infinity and the 1464 * second argument is less than zero, then the result is positive zero.</li> 1465 * <li>If the first argument is positive zero and the second argument is 1466 * less than zero, or the first argument is positive infinity and the 1467 * second argument is greater than zero, then the result is positive 1468 * infinity.</li> 1469 * <li>If the first argument is negative zero and the second argument is 1470 * greater than zero but not a finite odd integer, or the first argument is 1471 * negative infinity and the second argument is less than zero but not a 1472 * finite odd integer, then the result is positive zero.</li> 1473 * <li>If the first argument is negative zero and the second argument is a 1474 * positive finite odd integer, or the first argument is negative infinity 1475 * and the second argument is a negative finite odd integer, then the result 1476 * is negative zero.</li> 1477 * <li>If the first argument is negative zero and the second argument is 1478 * less than zero but not a finite odd integer, or the first argument is 1479 * negative infinity and the second argument is greater than zero but not a 1480 * finite odd integer, then the result is positive infinity.</li> 1481 * <li>If the first argument is negative zero and the second argument is a 1482 * negative finite odd integer, or the first argument is negative infinity 1483 * and the second argument is a positive finite odd integer, then the result 1484 * is negative infinity.</li> 1485 * <li>If the first argument is less than zero and the second argument is a 1486 * finite even integer, then the result is equal to the result of raising 1487 * the absolute value of the first argument to the power of the second 1488 * argument.</li> 1489 * <li>If the first argument is less than zero and the second argument is a 1490 * finite odd integer, then the result is equal to the negative of the 1491 * result of raising the absolute value of the first argument to the power 1492 * of the second argument.</li> 1493 * <li>If the first argument is finite and less than zero and the second 1494 * argument is finite and not an integer, then the result is NaN.</li> 1495 * <li>If both arguments are integers, then the result is exactly equal to 1496 * the mathematical result of raising the first argument to the power of 1497 * the second argument if that result can in fact be represented exactly as 1498 * a double value.</li> 1499 * 1500 * </ul><p>(In the foregoing descriptions, a floating-point value is 1501 * considered to be an integer if and only if it is a fixed point of the 1502 * method {@link #ceil(double)} or, equivalently, a fixed point of the 1503 * method {@link #floor(double)}. A value is a fixed point of a one-argument 1504 * method if and only if the result of applying the method to the value is 1505 * equal to the value.) 1506 * 1507 * @param x the number to raise 1508 * @param y the power to raise it to 1509 * @return x<sup>y</sup> 1510 */ 1511 public static double pow(double x, double y) 1512 { 1513 // Special cases first. 1514 if (y == 0) 1515 return 1; 1516 if (y == 1) 1517 return x; 1518 if (y == -1) 1519 return 1 / x; 1520 if (x != x || y != y) 1521 return Double.NaN; 1522 1523 // When x < 0, yisint tells if y is not an integer (0), even(1), 1524 // or odd (2). 1525 int yisint = 0; 1526 if (x < 0 && floor(y) == y) 1527 yisint = (y % 2 == 0) ? 2 : 1; 1528 double ax = abs(x); 1529 double ay = abs(y); 1530 1531 // More special cases, of y. 1532 if (ay == Double.POSITIVE_INFINITY) 1533 { 1534 if (ax == 1) 1535 return Double.NaN; 1536 if (ax > 1) 1537 return y > 0 ? y : 0; 1538 return y < 0 ? -y : 0; 1539 } 1540 if (y == 2) 1541 return x * x; 1542 if (y == 0.5) 1543 return sqrt(x); 1544 1545 // More special cases, of x. 1546 if (x == 0 || ax == Double.POSITIVE_INFINITY || ax == 1) 1547 { 1548 if (y < 0) 1549 ax = 1 / ax; 1550 if (x < 0) 1551 { 1552 if (x == -1 && yisint == 0) 1553 ax = Double.NaN; 1554 else if (yisint == 1) 1555 ax = -ax; 1556 } 1557 return ax; 1558 } 1559 if (x < 0 && yisint == 0) 1560 return Double.NaN; 1561 1562 // Now we can start! 1563 double t; 1564 double t1; 1565 double t2; 1566 double u; 1567 double v; 1568 double w; 1569 if (ay > TWO_31) 1570 { 1571 if (ay > TWO_64) // Automatic over/underflow. 1572 return ((ax < 1) ? y < 0 : y > 0) ? Double.POSITIVE_INFINITY : 0; 1573 // Over/underflow if x is not close to one. 1574 if (ax < 0.9999995231628418) 1575 return y < 0 ? Double.POSITIVE_INFINITY : 0; 1576 if (ax >= 1.0000009536743164) 1577 return y > 0 ? Double.POSITIVE_INFINITY : 0; 1578 // Now |1-x| is <= 2**-20, sufficient to compute 1579 // log(x) by x-x^2/2+x^3/3-x^4/4. 1580 t = x - 1; 1581 w = t * t * (0.5 - t * (1 / 3.0 - t * 0.25)); 1582 u = INV_LN2_H * t; 1583 v = t * INV_LN2_L - w * INV_LN2; 1584 t1 = (float) (u + v); 1585 t2 = v - (t1 - u); 1586 } 1587 else 1588 { 1589 long bits = Double.doubleToLongBits(ax); 1590 int exp = (int) (bits >> 52); 1591 if (exp == 0) // Subnormal x. 1592 { 1593 ax *= TWO_54; 1594 bits = Double.doubleToLongBits(ax); 1595 exp = (int) (bits >> 52) - 54; 1596 } 1597 exp -= 1023; // Unbias exponent. 1598 ax = Double.longBitsToDouble((bits & 0x000fffffffffffffL) 1599 | 0x3ff0000000000000L); 1600 boolean k; 1601 if (ax < SQRT_1_5) // |x|<sqrt(3/2). 1602 k = false; 1603 else if (ax < SQRT_3) // |x|<sqrt(3). 1604 k = true; 1605 else 1606 { 1607 k = false; 1608 ax *= 0.5; 1609 exp++; 1610 } 1611 1612 // Compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5). 1613 u = ax - (k ? 1.5 : 1); 1614 v = 1 / (ax + (k ? 1.5 : 1)); 1615 double s = u * v; 1616 double s_h = (float) s; 1617 double t_h = (float) (ax + (k ? 1.5 : 1)); 1618 double t_l = ax - (t_h - (k ? 1.5 : 1)); 1619 double s_l = v * ((u - s_h * t_h) - s_h * t_l); 1620 // Compute log(ax). 1621 double s2 = s * s; 1622 double r = s_l * (s_h + s) + s2 * s2 1623 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); 1624 s2 = s_h * s_h; 1625 t_h = (float) (3.0 + s2 + r); 1626 t_l = r - (t_h - 3.0 - s2); 1627 // u+v = s*(1+...). 1628 u = s_h * t_h; 1629 v = s_l * t_h + t_l * s; 1630 // 2/(3log2)*(s+...). 1631 double p_h = (float) (u + v); 1632 double p_l = v - (p_h - u); 1633 double z_h = CP_H * p_h; 1634 double z_l = CP_L * p_h + p_l * CP + (k ? DP_L : 0); 1635 // log2(ax) = (s+..)*2/(3*log2) = exp + dp_h + z_h + z_l. 1636 t = exp; 1637 t1 = (float) (z_h + z_l + (k ? DP_H : 0) + t); 1638 t2 = z_l - (t1 - t - (k ? DP_H : 0) - z_h); 1639 } 1640 1641 // Split up y into y1+y2 and compute (y1+y2)*(t1+t2). 1642 boolean negative = x < 0 && yisint == 1; 1643 double y1 = (float) y; 1644 double p_l = (y - y1) * t1 + y * t2; 1645 double p_h = y1 * t1; 1646 double z = p_l + p_h; 1647 if (z >= 1024) // Detect overflow. 1648 { 1649 if (z > 1024 || p_l + OVT > z - p_h) 1650 return negative ? Double.NEGATIVE_INFINITY 1651 : Double.POSITIVE_INFINITY; 1652 } 1653 else if (z <= -1075) // Detect underflow. 1654 { 1655 if (z < -1075 || p_l <= z - p_h) 1656 return negative ? -0.0 : 0; 1657 } 1658 1659 // Compute 2**(p_h+p_l). 1660 int n = round((float) z); 1661 p_h -= n; 1662 t = (float) (p_l + p_h); 1663 u = t * LN2_H; 1664 v = (p_l - (t - p_h)) * LN2 + t * LN2_L; 1665 z = u + v; 1666 w = v - (z - u); 1667 t = z * z; 1668 t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); 1669 double r = (z * t1) / (t1 - 2) - (w + z * w); 1670 z = scale(1 - (r - z), n); 1671 return negative ? -z : z; 1672 } 1673 1674 /** 1675 * Get the IEEE 754 floating point remainder on two numbers. This is the 1676 * value of <code>x - y * <em>n</em></code>, where <em>n</em> is the closest 1677 * double to <code>x / y</code> (ties go to the even n); for a zero 1678 * remainder, the sign is that of <code>x</code>. If either argument is NaN, 1679 * the first argument is infinite, or the second argument is zero, the result 1680 * is NaN; if x is finite but y is infinite, the result is x. 1681 * 1682 * @param x the dividend (the top half) 1683 * @param y the divisor (the bottom half) 1684 * @return the IEEE 754-defined floating point remainder of x/y 1685 * @see #rint(double) 1686 */ 1687 public static double IEEEremainder(double x, double y) 1688 { 1689 // Purge off exception values. 1690 if (x == Double.NEGATIVE_INFINITY || ! (x < Double.POSITIVE_INFINITY) 1691 || y == 0 || y != y) 1692 return Double.NaN; 1693 1694 boolean negative = x < 0; 1695 x = abs(x); 1696 y = abs(y); 1697 if (x == y || x == 0) 1698 return 0 * x; // Get correct sign. 1699 1700 // Achieve x < 2y, then take first shot at remainder. 1701 if (y < TWO_1023) 1702 x %= y + y; 1703 1704 // Now adjust x to get correct precision. 1705 if (y < 4 / TWO_1023) 1706 { 1707 if (x + x > y) 1708 { 1709 x -= y; 1710 if (x + x >= y) 1711 x -= y; 1712 } 1713 } 1714 else 1715 { 1716 y *= 0.5; 1717 if (x > y) 1718 { 1719 x -= y; 1720 if (x >= y) 1721 x -= y; 1722 } 1723 } 1724 return negative ? -x : x; 1725 } 1726 1727 /** 1728 * Take the nearest integer that is that is greater than or equal to the 1729 * argument. If the argument is NaN, infinite, or zero, the result is the 1730 * same; if the argument is between -1 and 0, the result is negative zero. 1731 * Note that <code>Math.ceil(x) == -Math.floor(-x)</code>. 1732 * 1733 * @param a the value to act upon 1734 * @return the nearest integer >= <code>a</code> 1735 */ 1736 public static double ceil(double a) 1737 { 1738 return -floor(-a); 1739 } 1740 1741 /** 1742 * Take the nearest integer that is that is less than or equal to the 1743 * argument. If the argument is NaN, infinite, or zero, the result is the 1744 * same. Note that <code>Math.ceil(x) == -Math.floor(-x)</code>. 1745 * 1746 * @param a the value to act upon 1747 * @return the nearest integer <= <code>a</code> 1748 */ 1749 public static double floor(double a) 1750 { 1751 double x = abs(a); 1752 if (! (x < TWO_52) || (long) a == a) 1753 return a; // No fraction bits; includes NaN and infinity. 1754 if (x < 1) 1755 return a >= 0 ? 0 * a : -1; // Worry about signed zero. 1756 return a < 0 ? (long) a - 1.0 : (long) a; // Cast to long truncates. 1757 } 1758 1759 /** 1760 * Take the nearest integer to the argument. If it is exactly between 1761 * two integers, the even integer is taken. If the argument is NaN, 1762 * infinite, or zero, the result is the same. 1763 * 1764 * @param a the value to act upon 1765 * @return the nearest integer to <code>a</code> 1766 */ 1767 public static double rint(double a) 1768 { 1769 double x = abs(a); 1770 if (! (x < TWO_52)) 1771 return a; // No fraction bits; includes NaN and infinity. 1772 if (x <= 0.5) 1773 return 0 * a; // Worry about signed zero. 1774 if (x % 2 <= 0.5) 1775 return (long) a; // Catch round down to even. 1776 return (long) (a + (a < 0 ? -0.5 : 0.5)); // Cast to long truncates. 1777 } 1778 1779 /** 1780 * Take the nearest integer to the argument. This is equivalent to 1781 * <code>(int) Math.floor(f + 0.5f)</code>. If the argument is NaN, the 1782 * result is 0; otherwise if the argument is outside the range of int, the 1783 * result will be Integer.MIN_VALUE or Integer.MAX_VALUE, as appropriate. 1784 * 1785 * @param f the argument to round 1786 * @return the nearest integer to the argument 1787 * @see Integer#MIN_VALUE 1788 * @see Integer#MAX_VALUE 1789 */ 1790 public static int round(float f) 1791 { 1792 return (int) floor(f + 0.5f); 1793 } 1794 1795 /** 1796 * Take the nearest long to the argument. This is equivalent to 1797 * <code>(long) Math.floor(d + 0.5)</code>. If the argument is NaN, the 1798 * result is 0; otherwise if the argument is outside the range of long, the 1799 * result will be Long.MIN_VALUE or Long.MAX_VALUE, as appropriate. 1800 * 1801 * @param d the argument to round 1802 * @return the nearest long to the argument 1803 * @see Long#MIN_VALUE 1804 * @see Long#MAX_VALUE 1805 */ 1806 public static long round(double d) 1807 { 1808 return (long) floor(d + 0.5); 1809 } 1810 1811 /** 1812 * Get a random number. This behaves like Random.nextDouble(), seeded by 1813 * System.currentTimeMillis() when first called. In other words, the number 1814 * is from a pseudorandom sequence, and lies in the range [+0.0, 1.0). 1815 * This random sequence is only used by this method, and is threadsafe, 1816 * although you may want your own random number generator if it is shared 1817 * among threads. 1818 * 1819 * @return a random number 1820 * @see Random#nextDouble() 1821 * @see System#currentTimeMillis() 1822 */ 1823 public static synchronized double random() 1824 { 1825 if (rand == null) 1826 rand = new Random(); 1827 return rand.nextDouble(); 1828 } 1829 1830 /** 1831 * Convert from degrees to radians. The formula for this is 1832 * radians = degrees * (pi/180); however it is not always exact given the 1833 * limitations of floating point numbers. 1834 * 1835 * @param degrees an angle in degrees 1836 * @return the angle in radians 1837 */ 1838 public static double toRadians(double degrees) 1839 { 1840 return (degrees * PI) / 180; 1841 } 1842 1843 /** 1844 * Convert from radians to degrees. The formula for this is 1845 * degrees = radians * (180/pi); however it is not always exact given the 1846 * limitations of floating point numbers. 1847 * 1848 * @param rads an angle in radians 1849 * @return the angle in degrees 1850 */ 1851 public static double toDegrees(double rads) 1852 { 1853 return (rads * 180) / PI; 1854 } 1855 1856 /** 1857 * Constants for scaling and comparing doubles by powers of 2. The compiler 1858 * must automatically inline constructs like (1/TWO_54), so we don't list 1859 * negative powers of two here. 1860 */ 1861 private static final double 1862 TWO_16 = 0x10000, // Long bits 0x40f0000000000000L. 1863 TWO_20 = 0x100000, // Long bits 0x4130000000000000L. 1864 TWO_24 = 0x1000000, // Long bits 0x4170000000000000L. 1865 TWO_27 = 0x8000000, // Long bits 0x41a0000000000000L. 1866 TWO_28 = 0x10000000, // Long bits 0x41b0000000000000L. 1867 TWO_29 = 0x20000000, // Long bits 0x41c0000000000000L. 1868 TWO_31 = 0x80000000L, // Long bits 0x41e0000000000000L. 1869 TWO_49 = 0x2000000000000L, // Long bits 0x4300000000000000L. 1870 TWO_52 = 0x10000000000000L, // Long bits 0x4330000000000000L. 1871 TWO_54 = 0x40000000000000L, // Long bits 0x4350000000000000L. 1872 TWO_57 = 0x200000000000000L, // Long bits 0x4380000000000000L. 1873 TWO_60 = 0x1000000000000000L, // Long bits 0x43b0000000000000L. 1874 TWO_64 = 1.8446744073709552e19, // Long bits 0x43f0000000000000L. 1875 TWO_66 = 7.378697629483821e19, // Long bits 0x4410000000000000L. 1876 TWO_1023 = 8.98846567431158e307; // Long bits 0x7fe0000000000000L. 1877 1878 /** 1879 * Super precision for 2/pi in 24-bit chunks, for use in 1880 * {@link #remPiOver2(double, double[])}. 1881 */ 1882 private static final int TWO_OVER_PI[] = { 1883 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62, 1884 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a, 1885 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129, 1886 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41, 1887 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8, 1888 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf, 1889 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5, 1890 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08, 1891 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3, 1892 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880, 1893 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b, 1894 }; 1895 1896 /** 1897 * Super precision for pi/2 in 24-bit chunks, for use in 1898 * {@link #remPiOver2(double, double[])}. 1899 */ 1900 private static final double PI_OVER_TWO[] = { 1901 1.570796251296997, // Long bits 0x3ff921fb40000000L. 1902 7.549789415861596e-8, // Long bits 0x3e74442d00000000L. 1903 5.390302529957765e-15, // Long bits 0x3cf8469880000000L. 1904 3.282003415807913e-22, // Long bits 0x3b78cc5160000000L. 1905 1.270655753080676e-29, // Long bits 0x39f01b8380000000L. 1906 1.2293330898111133e-36, // Long bits 0x387a252040000000L. 1907 2.7337005381646456e-44, // Long bits 0x36e3822280000000L. 1908 2.1674168387780482e-51, // Long bits 0x3569f31d00000000L. 1909 }; 1910 1911 /** 1912 * More constants related to pi, used in 1913 * {@link #remPiOver2(double, double[])} and elsewhere. 1914 */ 1915 private static final double 1916 PI_L = 1.2246467991473532e-16, // Long bits 0x3ca1a62633145c07L. 1917 PIO2_1 = 1.5707963267341256, // Long bits 0x3ff921fb54400000L. 1918 PIO2_1L = 6.077100506506192e-11, // Long bits 0x3dd0b4611a626331L. 1919 PIO2_2 = 6.077100506303966e-11, // Long bits 0x3dd0b4611a600000L. 1920 PIO2_2L = 2.0222662487959506e-21, // Long bits 0x3ba3198a2e037073L. 1921 PIO2_3 = 2.0222662487111665e-21, // Long bits 0x3ba3198a2e000000L. 1922 PIO2_3L = 8.4784276603689e-32; // Long bits 0x397b839a252049c1L. 1923 1924 /** 1925 * Natural log and square root constants, for calculation of 1926 * {@link #exp(double)}, {@link #log(double)} and 1927 * {@link #pow(double, double)}. CP is 2/(3*ln(2)). 1928 */ 1929 private static final double 1930 SQRT_1_5 = 1.224744871391589, // Long bits 0x3ff3988e1409212eL. 1931 SQRT_2 = 1.4142135623730951, // Long bits 0x3ff6a09e667f3bcdL. 1932 SQRT_3 = 1.7320508075688772, // Long bits 0x3ffbb67ae8584caaL. 1933 EXP_LIMIT_H = 709.782712893384, // Long bits 0x40862e42fefa39efL. 1934 EXP_LIMIT_L = -745.1332191019411, // Long bits 0xc0874910d52d3051L. 1935 CP = 0.9617966939259756, // Long bits 0x3feec709dc3a03fdL. 1936 CP_H = 0.9617967009544373, // Long bits 0x3feec709e0000000L. 1937 CP_L = -7.028461650952758e-9, // Long bits 0xbe3e2fe0145b01f5L. 1938 LN2 = 0.6931471805599453, // Long bits 0x3fe62e42fefa39efL. 1939 LN2_H = 0.6931471803691238, // Long bits 0x3fe62e42fee00000L. 1940 LN2_L = 1.9082149292705877e-10, // Long bits 0x3dea39ef35793c76L. 1941 INV_LN2 = 1.4426950408889634, // Long bits 0x3ff71547652b82feL. 1942 INV_LN2_H = 1.4426950216293335, // Long bits 0x3ff7154760000000L. 1943 INV_LN2_L = 1.9259629911266175e-8; // Long bits 0x3e54ae0bf85ddf44L. 1944 1945 /** 1946 * Constants for computing {@link #log(double)}. 1947 */ 1948 private static final double 1949 LG1 = 0.6666666666666735, // Long bits 0x3fe5555555555593L. 1950 LG2 = 0.3999999999940942, // Long bits 0x3fd999999997fa04L. 1951 LG3 = 0.2857142874366239, // Long bits 0x3fd2492494229359L. 1952 LG4 = 0.22222198432149784, // Long bits 0x3fcc71c51d8e78afL. 1953 LG5 = 0.1818357216161805, // Long bits 0x3fc7466496cb03deL. 1954 LG6 = 0.15313837699209373, // Long bits 0x3fc39a09d078c69fL. 1955 LG7 = 0.14798198605116586; // Long bits 0x3fc2f112df3e5244L. 1956 1957 /** 1958 * Constants for computing {@link #pow(double, double)}. L and P are 1959 * coefficients for series; OVT is -(1024-log2(ovfl+.5ulp)); and DP is ???. 1960 * The P coefficients also calculate {@link #exp(double)}. 1961 */ 1962 private static final double 1963 L1 = 0.5999999999999946, // Long bits 0x3fe3333333333303L. 1964 L2 = 0.4285714285785502, // Long bits 0x3fdb6db6db6fabffL. 1965 L3 = 0.33333332981837743, // Long bits 0x3fd55555518f264dL. 1966 L4 = 0.272728123808534, // Long bits 0x3fd17460a91d4101L. 1967 L5 = 0.23066074577556175, // Long bits 0x3fcd864a93c9db65L. 1968 L6 = 0.20697501780033842, // Long bits 0x3fca7e284a454eefL. 1969 P1 = 0.16666666666666602, // Long bits 0x3fc555555555553eL. 1970 P2 = -2.7777777777015593e-3, // Long bits 0xbf66c16c16bebd93L. 1971 P3 = 6.613756321437934e-5, // Long bits 0x3f11566aaf25de2cL. 1972 P4 = -1.6533902205465252e-6, // Long bits 0xbebbbd41c5d26bf1L. 1973 P5 = 4.1381367970572385e-8, // Long bits 0x3e66376972bea4d0L. 1974 DP_H = 0.5849624872207642, // Long bits 0x3fe2b80340000000L. 1975 DP_L = 1.350039202129749e-8, // Long bits 0x3e4cfdeb43cfd006L. 1976 OVT = 8.008566259537294e-17; // Long bits 0x3c971547652b82feL. 1977 1978 /** 1979 * Coefficients for computing {@link #sin(double)}. 1980 */ 1981 private static final double 1982 S1 = -0.16666666666666632, // Long bits 0xbfc5555555555549L. 1983 S2 = 8.33333333332249e-3, // Long bits 0x3f8111111110f8a6L. 1984 S3 = -1.984126982985795e-4, // Long bits 0xbf2a01a019c161d5L. 1985 S4 = 2.7557313707070068e-6, // Long bits 0x3ec71de357b1fe7dL. 1986 S5 = -2.5050760253406863e-8, // Long bits 0xbe5ae5e68a2b9cebL. 1987 S6 = 1.58969099521155e-10; // Long bits 0x3de5d93a5acfd57cL. 1988 1989 /** 1990 * Coefficients for computing {@link #cos(double)}. 1991 */ 1992 private static final double 1993 C1 = 0.0416666666666666, // Long bits 0x3fa555555555554cL. 1994 C2 = -1.388888888887411e-3, // Long bits 0xbf56c16c16c15177L. 1995 C3 = 2.480158728947673e-5, // Long bits 0x3efa01a019cb1590L. 1996 C4 = -2.7557314351390663e-7, // Long bits 0xbe927e4f809c52adL. 1997 C5 = 2.087572321298175e-9, // Long bits 0x3e21ee9ebdb4b1c4L. 1998 C6 = -1.1359647557788195e-11; // Long bits 0xbda8fae9be8838d4L. 1999 2000 /** 2001 * Coefficients for computing {@link #tan(double)}. 2002 */ 2003 private static final double 2004 T0 = 0.3333333333333341, // Long bits 0x3fd5555555555563L. 2005 T1 = 0.13333333333320124, // Long bits 0x3fc111111110fe7aL. 2006 T2 = 0.05396825397622605, // Long bits 0x3faba1ba1bb341feL. 2007 T3 = 0.021869488294859542, // Long bits 0x3f9664f48406d637L. 2008 T4 = 8.8632398235993e-3, // Long bits 0x3f8226e3e96e8493L. 2009 T5 = 3.5920791075913124e-3, // Long bits 0x3f6d6d22c9560328L. 2010 T6 = 1.4562094543252903e-3, // Long bits 0x3f57dbc8fee08315L. 2011 T7 = 5.880412408202641e-4, // Long bits 0x3f4344d8f2f26501L. 2012 T8 = 2.464631348184699e-4, // Long bits 0x3f3026f71a8d1068L. 2013 T9 = 7.817944429395571e-5, // Long bits 0x3f147e88a03792a6L. 2014 T10 = 7.140724913826082e-5, // Long bits 0x3f12b80f32f0a7e9L. 2015 T11 = -1.8558637485527546e-5, // Long bits 0xbef375cbdb605373L. 2016 T12 = 2.590730518636337e-5; // Long bits 0x3efb2a7074bf7ad4L. 2017 2018 /** 2019 * Coefficients for computing {@link #asin(double)} and 2020 * {@link #acos(double)}. 2021 */ 2022 private static final double 2023 PS0 = 0.16666666666666666, // Long bits 0x3fc5555555555555L. 2024 PS1 = -0.3255658186224009, // Long bits 0xbfd4d61203eb6f7dL. 2025 PS2 = 0.20121253213486293, // Long bits 0x3fc9c1550e884455L. 2026 PS3 = -0.04005553450067941, // Long bits 0xbfa48228b5688f3bL. 2027 PS4 = 7.915349942898145e-4, // Long bits 0x3f49efe07501b288L. 2028 PS5 = 3.479331075960212e-5, // Long bits 0x3f023de10dfdf709L. 2029 QS1 = -2.403394911734414, // Long bits 0xc0033a271c8a2d4bL. 2030 QS2 = 2.0209457602335057, // Long bits 0x40002ae59c598ac8L. 2031 QS3 = -0.6882839716054533, // Long bits 0xbfe6066c1b8d0159L. 2032 QS4 = 0.07703815055590194; // Long bits 0x3fb3b8c5b12e9282L. 2033 2034 /** 2035 * Coefficients for computing {@link #atan(double)}. 2036 */ 2037 private static final double 2038 ATAN_0_5H = 0.4636476090008061, // Long bits 0x3fddac670561bb4fL. 2039 ATAN_0_5L = 2.2698777452961687e-17, // Long bits 0x3c7a2b7f222f65e2L. 2040 ATAN_1_5H = 0.982793723247329, // Long bits 0x3fef730bd281f69bL. 2041 ATAN_1_5L = 1.3903311031230998e-17, // Long bits 0x3c7007887af0cbbdL. 2042 AT0 = 0.3333333333333293, // Long bits 0x3fd555555555550dL. 2043 AT1 = -0.19999999999876483, // Long bits 0xbfc999999998ebc4L. 2044 AT2 = 0.14285714272503466, // Long bits 0x3fc24924920083ffL. 2045 AT3 = -0.11111110405462356, // Long bits 0xbfbc71c6fe231671L. 2046 AT4 = 0.09090887133436507, // Long bits 0x3fb745cdc54c206eL. 2047 AT5 = -0.0769187620504483, // Long bits 0xbfb3b0f2af749a6dL. 2048 AT6 = 0.06661073137387531, // Long bits 0x3fb10d66a0d03d51L. 2049 AT7 = -0.058335701337905735, // Long bits 0xbfadde2d52defd9aL. 2050 AT8 = 0.049768779946159324, // Long bits 0x3fa97b4b24760debL. 2051 AT9 = -0.036531572744216916, // Long bits 0xbfa2b4442c6a6c2fL. 2052 AT10 = 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L. 2053 2054 /** 2055 * Constants for computing {@link #cbrt(double)}. 2056 */ 2057 private static final int 2058 CBRT_B1 = 715094163, // B1 = (682-0.03306235651)*2**20 2059 CBRT_B2 = 696219795; // B2 = (664-0.03306235651)*2**20 2060 2061 /** 2062 * Constants for computing {@link #cbrt(double)}. 2063 */ 2064 private static final double 2065 CBRT_C = 5.42857142857142815906e-01, // Long bits 0x3fe15f15f15f15f1L 2066 CBRT_D = -7.05306122448979611050e-01, // Long bits 0xbfe691de2532c834L 2067 CBRT_E = 1.41428571428571436819e+00, // Long bits 0x3ff6a0ea0ea0ea0fL 2068 CBRT_F = 1.60714285714285720630e+00, // Long bits 0x3ff9b6db6db6db6eL 2069 CBRT_G = 3.57142857142857150787e-01; // Long bits 0x3fd6db6db6db6db7L 2070 2071 /** 2072 * Constants for computing {@link #expm1(double)} 2073 */ 2074 private static final double 2075 EXPM1_Q1 = -3.33333333333331316428e-02, // Long bits 0xbfa11111111110f4L 2076 EXPM1_Q2 = 1.58730158725481460165e-03, // Long bits 0x3f5a01a019fe5585L 2077 EXPM1_Q3 = -7.93650757867487942473e-05, // Long bits 0xbf14ce199eaadbb7L 2078 EXPM1_Q4 = 4.00821782732936239552e-06, // Long bits 0x3ed0cfca86e65239L 2079 EXPM1_Q5 = -2.01099218183624371326e-07; // Long bits 0xbe8afdb76e09c32dL 2080 2081 /** 2082 * Helper function for reducing an angle to a multiple of pi/2 within 2083 * [-pi/4, pi/4]. 2084 * 2085 * @param x the angle; not infinity or NaN, and outside pi/4 2086 * @param y an array of 2 doubles modified to hold the remander x % pi/2 2087 * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4], 2088 * 1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4] 2089 */ 2090 private static int remPiOver2(double x, double[] y) 2091 { 2092 boolean negative = x < 0; 2093 x = abs(x); 2094 double z; 2095 int n; 2096 if (Configuration.DEBUG && (x <= PI / 4 || x != x 2097 || x == Double.POSITIVE_INFINITY)) 2098 throw new InternalError("Assertion failure"); 2099 if (x < 3 * PI / 4) // If |x| is small. 2100 { 2101 z = x - PIO2_1; 2102 if ((float) x != (float) (PI / 2)) // 33+53 bit pi is good enough. 2103 { 2104 y[0] = z - PIO2_1L; 2105 y[1] = z - y[0] - PIO2_1L; 2106 } 2107 else // Near pi/2, use 33+33+53 bit pi. 2108 { 2109 z -= PIO2_2; 2110 y[0] = z - PIO2_2L; 2111 y[1] = z - y[0] - PIO2_2L; 2112 } 2113 n = 1; 2114 } 2115 else if (x <= TWO_20 * PI / 2) // Medium size. 2116 { 2117 n = (int) (2 / PI * x + 0.5); 2118 z = x - n * PIO2_1; 2119 double w = n * PIO2_1L; // First round good to 85 bits. 2120 y[0] = z - w; 2121 if (n >= 32 || (float) x == (float) (w)) 2122 { 2123 if (x / y[0] >= TWO_16) // Second iteration, good to 118 bits. 2124 { 2125 double t = z; 2126 w = n * PIO2_2; 2127 z = t - w; 2128 w = n * PIO2_2L - (t - z - w); 2129 y[0] = z - w; 2130 if (x / y[0] >= TWO_49) // Third iteration, 151 bits accuracy. 2131 { 2132 t = z; 2133 w = n * PIO2_3; 2134 z = t - w; 2135 w = n * PIO2_3L - (t - z - w); 2136 y[0] = z - w; 2137 } 2138 } 2139 } 2140 y[1] = z - y[0] - w; 2141 } 2142 else 2143 { 2144 // All other (large) arguments. 2145 int e0 = (int) (Double.doubleToLongBits(x) >> 52) - 1046; 2146 z = scale(x, -e0); // e0 = ilogb(z) - 23. 2147 double[] tx = new double[3]; 2148 for (int i = 0; i < 2; i++) 2149 { 2150 tx[i] = (int) z; 2151 z = (z - tx[i]) * TWO_24; 2152 } 2153 tx[2] = z; 2154 int nx = 2; 2155 while (tx[nx] == 0) 2156 nx--; 2157 n = remPiOver2(tx, y, e0, nx); 2158 } 2159 if (negative) 2160 { 2161 y[0] = -y[0]; 2162 y[1] = -y[1]; 2163 return -n; 2164 } 2165 return n; 2166 } 2167 2168 /** 2169 * Helper function for reducing an angle to a multiple of pi/2 within 2170 * [-pi/4, pi/4]. 2171 * 2172 * @param x the positive angle, broken into 24-bit chunks 2173 * @param y an array of 2 doubles modified to hold the remander x % pi/2 2174 * @param e0 the exponent of x[0] 2175 * @param nx the last index used in x 2176 * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4], 2177 * 1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4] 2178 */ 2179 private static int remPiOver2(double[] x, double[] y, int e0, int nx) 2180 { 2181 int i; 2182 int ih; 2183 int n; 2184 double fw; 2185 double z; 2186 int[] iq = new int[20]; 2187 double[] f = new double[20]; 2188 double[] q = new double[20]; 2189 boolean recompute = false; 2190 2191 // Initialize jk, jz, jv, q0; note that 3>q0. 2192 int jk = 4; 2193 int jz = jk; 2194 int jv = max((e0 - 3) / 24, 0); 2195 int q0 = e0 - 24 * (jv + 1); 2196 2197 // Set up f[0] to f[nx+jk] where f[nx+jk] = TWO_OVER_PI[jv+jk]. 2198 int j = jv - nx; 2199 int m = nx + jk; 2200 for (i = 0; i <= m; i++, j++) 2201 f[i] = (j < 0) ? 0 : TWO_OVER_PI[j]; 2202 2203 // Compute q[0],q[1],...q[jk]. 2204 for (i = 0; i <= jk; i++) 2205 { 2206 for (j = 0, fw = 0; j <= nx; j++) 2207 fw += x[j] * f[nx + i - j]; 2208 q[i] = fw; 2209 } 2210 2211 do 2212 { 2213 // Distill q[] into iq[] reversingly. 2214 for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) 2215 { 2216 fw = (int) (1 / TWO_24 * z); 2217 iq[i] = (int) (z - TWO_24 * fw); 2218 z = q[j - 1] + fw; 2219 } 2220 2221 // Compute n. 2222 z = scale(z, q0); 2223 z -= 8 * floor(z * 0.125); // Trim off integer >= 8. 2224 n = (int) z; 2225 z -= n; 2226 ih = 0; 2227 if (q0 > 0) // Need iq[jz-1] to determine n. 2228 { 2229 i = iq[jz - 1] >> (24 - q0); 2230 n += i; 2231 iq[jz - 1] -= i << (24 - q0); 2232 ih = iq[jz - 1] >> (23 - q0); 2233 } 2234 else if (q0 == 0) 2235 ih = iq[jz - 1] >> 23; 2236 else if (z >= 0.5) 2237 ih = 2; 2238 2239 if (ih > 0) // If q > 0.5. 2240 { 2241 n += 1; 2242 int carry = 0; 2243 for (i = 0; i < jz; i++) // Compute 1-q. 2244 { 2245 j = iq[i]; 2246 if (carry == 0) 2247 { 2248 if (j != 0) 2249 { 2250 carry = 1; 2251 iq[i] = 0x1000000 - j; 2252 } 2253 } 2254 else 2255 iq[i] = 0xffffff - j; 2256 } 2257 switch (q0) 2258 { 2259 case 1: // Rare case: chance is 1 in 12 for non-default. 2260 iq[jz - 1] &= 0x7fffff; 2261 break; 2262 case 2: 2263 iq[jz - 1] &= 0x3fffff; 2264 } 2265 if (ih == 2) 2266 { 2267 z = 1 - z; 2268 if (carry != 0) 2269 z -= scale(1, q0); 2270 } 2271 } 2272 2273 // Check if recomputation is needed. 2274 if (z == 0) 2275 { 2276 j = 0; 2277 for (i = jz - 1; i >= jk; i--) 2278 j |= iq[i]; 2279 if (j == 0) // Need recomputation. 2280 { 2281 int k; // k = no. of terms needed. 2282 for (k = 1; iq[jk - k] == 0; k++) 2283 ; 2284 2285 for (i = jz + 1; i <= jz + k; i++) // Add q[jz+1] to q[jz+k]. 2286 { 2287 f[nx + i] = TWO_OVER_PI[jv + i]; 2288 for (j = 0, fw = 0; j <= nx; j++) 2289 fw += x[j] * f[nx + i - j]; 2290 q[i] = fw; 2291 } 2292 jz += k; 2293 recompute = true; 2294 } 2295 } 2296 } 2297 while (recompute); 2298 2299 // Chop off zero terms. 2300 if (z == 0) 2301 { 2302 jz--; 2303 q0 -= 24; 2304 while (iq[jz] == 0) 2305 { 2306 jz--; 2307 q0 -= 24; 2308 } 2309 } 2310 else // Break z into 24-bit if necessary. 2311 { 2312 z = scale(z, -q0); 2313 if (z >= TWO_24) 2314 { 2315 fw = (int) (1 / TWO_24 * z); 2316 iq[jz] = (int) (z - TWO_24 * fw); 2317 jz++; 2318 q0 += 24; 2319 iq[jz] = (int) fw; 2320 } 2321 else 2322 iq[jz] = (int) z; 2323 } 2324 2325 // Convert integer "bit" chunk to floating-point value. 2326 fw = scale(1, q0); 2327 for (i = jz; i >= 0; i--) 2328 { 2329 q[i] = fw * iq[i]; 2330 fw *= 1 / TWO_24; 2331 } 2332 2333 // Compute PI_OVER_TWO[0,...,jk]*q[jz,...,0]. 2334 double[] fq = new double[20]; 2335 for (i = jz; i >= 0; i--) 2336 { 2337 fw = 0; 2338 for (int k = 0; k <= jk && k <= jz - i; k++) 2339 fw += PI_OVER_TWO[k] * q[i + k]; 2340 fq[jz - i] = fw; 2341 } 2342 2343 // Compress fq[] into y[]. 2344 fw = 0; 2345 for (i = jz; i >= 0; i--) 2346 fw += fq[i]; 2347 y[0] = (ih == 0) ? fw : -fw; 2348 fw = fq[0] - fw; 2349 for (i = 1; i <= jz; i++) 2350 fw += fq[i]; 2351 y[1] = (ih == 0) ? fw : -fw; 2352 return n; 2353 } 2354 2355 /** 2356 * Helper method for scaling a double by a power of 2. 2357 * 2358 * @param x the double 2359 * @param n the scale; |n| < 2048 2360 * @return x * 2**n 2361 */ 2362 private static double scale(double x, int n) 2363 { 2364 if (Configuration.DEBUG && abs(n) >= 2048) 2365 throw new InternalError("Assertion failure"); 2366 if (x == 0 || x == Double.NEGATIVE_INFINITY 2367 || ! (x < Double.POSITIVE_INFINITY) || n == 0) 2368 return x; 2369 long bits = Double.doubleToLongBits(x); 2370 int exp = (int) (bits >> 52) & 0x7ff; 2371 if (exp == 0) // Subnormal x. 2372 { 2373 x *= TWO_54; 2374 exp = ((int) (Double.doubleToLongBits(x) >> 52) & 0x7ff) - 54; 2375 } 2376 exp += n; 2377 if (exp > 0x7fe) // Overflow. 2378 return Double.POSITIVE_INFINITY * x; 2379 if (exp > 0) // Normal. 2380 return Double.longBitsToDouble((bits & 0x800fffffffffffffL) 2381 | ((long) exp << 52)); 2382 if (exp <= -54) 2383 return 0 * x; // Underflow. 2384 exp += 54; // Subnormal result. 2385 x = Double.longBitsToDouble((bits & 0x800fffffffffffffL) 2386 | ((long) exp << 52)); 2387 return x * (1 / TWO_54); 2388 } 2389 2390 /** 2391 * Helper trig function; computes sin in range [-pi/4, pi/4]. 2392 * 2393 * @param x angle within about pi/4 2394 * @param y tail of x, created by remPiOver2 2395 * @return sin(x+y) 2396 */ 2397 private static double sin(double x, double y) 2398 { 2399 if (Configuration.DEBUG && abs(x + y) > 0.7854) 2400 throw new InternalError("Assertion failure"); 2401 if (abs(x) < 1 / TWO_27) 2402 return x; // If |x| ~< 2**-27, already know answer. 2403 2404 double z = x * x; 2405 double v = z * x; 2406 double r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6))); 2407 if (y == 0) 2408 return x + v * (S1 + z * r); 2409 return x - ((z * (0.5 * y - v * r) - y) - v * S1); 2410 } 2411 2412 /** 2413 * Helper trig function; computes cos in range [-pi/4, pi/4]. 2414 * 2415 * @param x angle within about pi/4 2416 * @param y tail of x, created by remPiOver2 2417 * @return cos(x+y) 2418 */ 2419 private static double cos(double x, double y) 2420 { 2421 if (Configuration.DEBUG && abs(x + y) > 0.7854) 2422 throw new InternalError("Assertion failure"); 2423 x = abs(x); 2424 if (x < 1 / TWO_27) 2425 return 1; // If |x| ~< 2**-27, already know answer. 2426 2427 double z = x * x; 2428 double r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6))))); 2429 2430 if (x < 0.3) 2431 return 1 - (0.5 * z - (z * r - x * y)); 2432 2433 double qx = (x > 0.78125) ? 0.28125 : (x * 0.25); 2434 return 1 - qx - ((0.5 * z - qx) - (z * r - x * y)); 2435 } 2436 2437 /** 2438 * Helper trig function; computes tan in range [-pi/4, pi/4]. 2439 * 2440 * @param x angle within about pi/4 2441 * @param y tail of x, created by remPiOver2 2442 * @param invert true iff -1/tan should be returned instead 2443 * @return tan(x+y) 2444 */ 2445 private static double tan(double x, double y, boolean invert) 2446 { 2447 // PI/2 is irrational, so no double is a perfect multiple of it. 2448 if (Configuration.DEBUG && (abs(x + y) > 0.7854 || (x == 0 && invert))) 2449 throw new InternalError("Assertion failure"); 2450 boolean negative = x < 0; 2451 if (negative) 2452 { 2453 x = -x; 2454 y = -y; 2455 } 2456 if (x < 1 / TWO_28) // If |x| ~< 2**-28, already know answer. 2457 return (negative ? -1 : 1) * (invert ? -1 / x : x); 2458 2459 double z; 2460 double w; 2461 boolean large = x >= 0.6744; 2462 if (large) 2463 { 2464 z = PI / 4 - x; 2465 w = PI_L / 4 - y; 2466 x = z + w; 2467 y = 0; 2468 } 2469 z = x * x; 2470 w = z * z; 2471 // Break x**5*(T1+x**2*T2+...) into 2472 // x**5(T1+x**4*T3+...+x**20*T11) 2473 // + x**5(x**2*(T2+x**4*T4+...+x**22*T12)). 2474 double r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11)))); 2475 double v = z * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12))))); 2476 double s = z * x; 2477 r = y + z * (s * (r + v) + y); 2478 r += T0 * s; 2479 w = x + r; 2480 if (large) 2481 { 2482 v = invert ? -1 : 1; 2483 return (negative ? -1 : 1) * (v - 2 * (x - (w * w / (w + v) - r))); 2484 } 2485 if (! invert) 2486 return w; 2487 2488 // Compute -1.0/(x+r) accurately. 2489 z = (float) w; 2490 v = r - (z - x); 2491 double a = -1 / w; 2492 double t = (float) a; 2493 return t + a * (1 + t * z + t * v); 2494 } 2495 2496 /** 2497 * <p> 2498 * Returns the sign of the argument as follows: 2499 * </p> 2500 * <ul> 2501 * <li>If <code>a</code> is greater than zero, the result is 1.0.</li> 2502 * <li>If <code>a</code> is less than zero, the result is -1.0.</li> 2503 * <li>If <code>a</code> is <code>NaN</code>, the result is <code>NaN</code>. 2504 * <li>If <code>a</code> is positive or negative zero, the result is the 2505 * same.</li> 2506 * </ul> 2507 * 2508 * @param a the numeric argument. 2509 * @return the sign of the argument. 2510 * @since 1.5. 2511 */ 2512 public static double signum(double a) 2513 { 2514 // There's no difference. 2515 return Math.signum(a); 2516 } 2517 2518 /** 2519 * <p> 2520 * Returns the sign of the argument as follows: 2521 * </p> 2522 * <ul> 2523 * <li>If <code>a</code> is greater than zero, the result is 1.0f.</li> 2524 * <li>If <code>a</code> is less than zero, the result is -1.0f.</li> 2525 * <li>If <code>a</code> is <code>NaN</code>, the result is <code>NaN</code>. 2526 * <li>If <code>a</code> is positive or negative zero, the result is the 2527 * same.</li> 2528 * </ul> 2529 * 2530 * @param a the numeric argument. 2531 * @return the sign of the argument. 2532 * @since 1.5. 2533 */ 2534 public static float signum(float a) 2535 { 2536 // There's no difference. 2537 return Math.signum(a); 2538 } 2539 2540 /** 2541 * Return the ulp for the given double argument. The ulp is the 2542 * difference between the argument and the next larger double. Note 2543 * that the sign of the double argument is ignored, that is, 2544 * ulp(x) == ulp(-x). If the argument is a NaN, then NaN is returned. 2545 * If the argument is an infinity, then +Inf is returned. If the 2546 * argument is zero (either positive or negative), then 2547 * {@link Double#MIN_VALUE} is returned. 2548 * @param d the double whose ulp should be returned 2549 * @return the difference between the argument and the next larger double 2550 * @since 1.5 2551 */ 2552 public static double ulp(double d) 2553 { 2554 // There's no difference. 2555 return Math.ulp(d); 2556 } 2557 2558 /** 2559 * Return the ulp for the given float argument. The ulp is the 2560 * difference between the argument and the next larger float. Note 2561 * that the sign of the float argument is ignored, that is, 2562 * ulp(x) == ulp(-x). If the argument is a NaN, then NaN is returned. 2563 * If the argument is an infinity, then +Inf is returned. If the 2564 * argument is zero (either positive or negative), then 2565 * {@link Float#MIN_VALUE} is returned. 2566 * @param f the float whose ulp should be returned 2567 * @return the difference between the argument and the next larger float 2568 * @since 1.5 2569 */ 2570 public static float ulp(float f) 2571 { 2572 // There's no difference. 2573 return Math.ulp(f); 2574 } 2575 } 2576