| Frames | No Frames |
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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: */ 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: }