| GNU Classpath (0.95) | |
| 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