Source for java.lang.StrictMath

   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 &gt;= <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 &lt;= <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: }