libjava/classpath/java/lang/StrictMath.java

   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 cosine 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 positive
 645    * infinity.</li>
 646    * <li>If the argument is zero, the result is one.</li>
 647    * </ul>
 648    *
 649    * @param x the argument to <em>cosh</em>
 650    * @return the hyperbolic cosine of <code>x</code>
 651    *
 652    * @since 1.5
 653    */
 654   public static double cosh(double x)
 655   {
 656     // Method :
 657     // mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
 658     // 1. Replace x by |x| (cosh(x) = cosh(-x)).
 659     // 2.
 660     //                                             [ exp(x) - 1 ]^2
 661     //  0        <= x <= ln2/2  :  cosh(x) := 1 + -------------------
 662     //                                                 2*exp(x)
 663     //
 664     //                                        exp(x) +  1/exp(x)
 665     //  ln2/2    <= x <= 22     :  cosh(x) := ------------------
 666     //                                               2
 667     //  22       <= x <= lnovft :  cosh(x) := exp(x)/2
 668     //  lnovft   <= x <= ln2ovft:  cosh(x) := exp(x/2)/2 * exp(x/2)
 669     //  ln2ovft  <  x           :  cosh(x) := +inf  (overflow)
 670
 671     double t, w;
 672     long bits;
 673     int hx;
 674     int lx;
 675
 676     // handle special cases
 677     if (x != x)
 678       return Double.NaN;
 679     if (x == Double.POSITIVE_INFINITY)
 680       return Double.POSITIVE_INFINITY;
 681     if (x == Double.NEGATIVE_INFINITY)
 682       return Double.POSITIVE_INFINITY;
 683
 684     bits = Double.doubleToLongBits(x);
 685     hx = getHighDWord(bits) & 0x7fffffff;  // ignore sign
 686     lx = getLowDWord(bits);
 687
 688     // |x| in [0, 0.5 * ln(2)], return 1 + expm1(|x|)^2 / (2 * exp(|x|))
 689     if (hx < 0x3fd62e43)
 690       {
 691         t = expm1(abs(x));
 692         w = 1.0 + t;
 693
 694         // for tiny arguments return 1.
 695         if (hx < 0x3c800000)
 696           return w;
 697
 698         return 1.0 + (t * t) / (w + w);
 699       }
 700
 701     // |x| in [0.5 * ln(2), 22], return exp(|x|)/2 + 1 / (2 * exp(|x|))
 702     if (hx < 0x40360000)
 703       {
 704         t = exp(abs(x));
 705
 706         return 0.5 * t + 0.5 / t;
 707       }
 708
 709     // |x| in [22, log(Double.MAX_VALUE)], return 0.5 * exp(|x|)
 710     if (hx < 0x40862e42)
 711       return 0.5 * exp(abs(x));
 712
 713     // |x| in [log(Double.MAX_VALUE), overflowthreshold],
 714     // return exp(x/2)/2 * exp(x/2)
 715
 716     // we need to force an unsigned <= compare, thus can not use lx.
 717     if ((hx < 0x408633ce)
 718         || ((hx == 0x408633ce)
 719             && ((bits & 0x00000000ffffffffL) <= 0x8fb9f87dL)))
 720       {
 721         w = exp(0.5 * abs(x));
 722         t = 0.5 * w;
 723
 724         return t * w;
 725       }
 726
 727     // |x| > overflowthreshold
 728     return Double.POSITIVE_INFINITY;
 729   }
 730
 731   /**
 732    * Returns the lower two words of a long. This is intended to be
 733    * used like this:
 734    * <code>getLowDWord(Double.doubleToLongBits(x))</code>.
 735    */
 736   private static int getLowDWord(long x)
 737   {
 738     return (int) (x & 0x00000000ffffffffL);
 739   }
 740
 741   /**
 742    * Returns the higher two words of a long. This is intended to be
 743    * used like this:
 744    * <code>getHighDWord(Double.doubleToLongBits(x))</code>.
 745    */
 746   private static int getHighDWord(long x)
 747   {
 748     return (int) ((x & 0xffffffff00000000L) >> 32);
 749   }
 750
 751   /**
 752    * Returns a double with the IEEE754 bit pattern given in the lower
 753    * and higher two words <code>lowDWord</code> and <code>highDWord</code>.
 754    */
 755   private static double buildDouble(int lowDWord, int highDWord)
 756   {
 757     return Double.longBitsToDouble((((long) highDWord & 0xffffffffL) << 32)
 758                                    | ((long) lowDWord & 0xffffffffL));
 759   }
 760
 761   /**
 762    * Returns the cube root of <code>x</code>. The sign of the cube root
 763    * is equal to the sign of <code>x</code>.
 764    *
 765    * Special cases:
 766    * <ul>
 767    * <li>If the argument is NaN, the result is NaN</li>
 768    * <li>If the argument is positive infinity, the result is positive
 769    * infinity.</li>
 770    * <li>If the argument is negative infinity, the result is negative
 771    * infinity.</li>
 772    * <li>If the argument is zero, the result is zero with the same
 773    * sign as the argument.</li>
 774    * </ul>
 775    *
 776    * @param x the number to take the cube root of
 777    * @return the cube root of <code>x</code>
 778    * @see #sqrt(double)
 779    *
 780    * @since 1.5
 781    */
 782   public static double cbrt(double x)
 783   {
 784     boolean negative = (x < 0);
 785     double r;
 786     double s;
 787     double t;
 788     double w;
 789
 790     long bits;
 791     int l;
 792     int h;
 793
 794     // handle the special cases
 795     if (x != x)
 796       return Double.NaN;
 797     if (x == Double.POSITIVE_INFINITY)
 798       return Double.POSITIVE_INFINITY;
 799     if (x == Double.NEGATIVE_INFINITY)
 800       return Double.NEGATIVE_INFINITY;
 801     if (x == 0)
 802       return x;
 803
 804     x = abs(x);
 805     bits = Double.doubleToLongBits(x);
 806
 807     if (bits < 0x0010000000000000L)   // subnormal number
 808       {
 809         t = TWO_54;
 810         t *= x;
 811
 812         // __HI(t)=__HI(t)/3+B2;
 813         bits = Double.doubleToLongBits(t);
 814         h = getHighDWord(bits);
 815         l = getLowDWord(bits);
 816
 817         h = h / 3 + CBRT_B2;
 818
 819         t = buildDouble(l, h);
 820       }
 821     else
 822       {
 823         // __HI(t)=__HI(x)/3+B1;
 824         h = getHighDWord(bits);
 825         l = 0;
 826
 827         h = h / 3 + CBRT_B1;
 828         t = buildDouble(l, h);
 829       }
 830
 831     // new cbrt to 23 bits
 832     r =  t * t / x;
 833     s =  CBRT_C + r * t;
 834     t *= CBRT_G + CBRT_F / (s + CBRT_E + CBRT_D / s);
 835
 836     // chopped to 20 bits and make it larger than cbrt(x)
 837     bits = Double.doubleToLongBits(t);
 838     h = getHighDWord(bits);
 839
 840     // __LO(t)=0;
 841     // __HI(t)+=0x00000001;
 842     l = 0;
 843     h += 1;
 844     t = buildDouble(l, h);
 845
 846     // one step newton iteration to 53 bits with error less than 0.667 ulps
 847     s = t * t;              // t * t is exact
 848     r = x / s;
 849     w = t + t;
 850     r = (r - t) / (w + r);  // r - s is exact
 851     t = t + t * r;
 852
 853     return negative ? -t : t;
 854   }
 855
 856   /**
 857    * Take <em>e</em><sup>a</sup>.  The opposite of <code>log()</code>. If the
 858    * argument is NaN, the result is NaN; if the argument is positive infinity,
 859    * the result is positive infinity; and if the argument is negative
 860    * infinity, the result is positive zero.
 861    *
 862    * @param x the number to raise to the power
 863    * @return the number raised to the power of <em>e</em>
 864    * @see #log(double)
 865    * @see #pow(double, double)
 866    */
 867   public static double exp(double x)
 868   {
 869     if (x != x)
 870       return x;
 871     if (x > EXP_LIMIT_H)
 872       return Double.POSITIVE_INFINITY;
 873     if (x < EXP_LIMIT_L)
 874       return 0;
 875
 876     // Argument reduction.
 877     double hi;
 878     double lo;
 879     int k;
 880     double t = abs(x);
 881     if (t > 0.5 * LN2)
 882       {
 883         if (t < 1.5 * LN2)
 884           {
 885             hi = t - LN2_H;
 886             lo = LN2_L;
 887             k = 1;
 888           }
 889         else
 890           {
 891             k = (int) (INV_LN2 * t + 0.5);
 892             hi = t - k * LN2_H;
 893             lo = k * LN2_L;
 894           }
 895         if (x < 0)
 896           {
 897             hi = -hi;
 898             lo = -lo;
 899             k = -k;
 900           }
 901         x = hi - lo;
 902       }
 903     else if (t < 1 / TWO_28)
 904       return 1;
 905     else
 906       lo = hi = k = 0;
 907
 908     // Now x is in primary range.
 909     t = x * x;
 910     double c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
 911     if (k == 0)
 912       return 1 - (x * c / (c - 2) - x);
 913     double y = 1 - (lo - x * c / (2 - c) - hi);
 914     return scale(y, k);
 915   }
 916
 917   /**
 918    * Returns <em>e</em><sup>x</sup> - 1.
 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 -1.</li>
 925    * <li>If the argument is zero, the result is zero.</li>
 926    * </ul>
 927    *
 928    * @param x the argument to <em>e</em><sup>x</sup> - 1.
 929    * @return <em>e</em> raised to the power <code>x</code> minus one.
 930    * @see #exp(double)
 931    */
 932   public static double expm1(double x)
 933   {
 934     // Method
 935     //   1. Argument reduction:
 936     //  Given x, find r and integer k such that
 937     //
 938     //            x = k * ln(2) + r,  |r| <= 0.5 * ln(2)
 939     //
 940     //  Here a correction term c will be computed to compensate
 941     //  the error in r when rounded to a floating-point number.
 942     //
 943     //   2. Approximating expm1(r) by a special rational function on
 944     //  the interval [0, 0.5 * ln(2)]:
 945     //  Since
 946     //      r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 - r^4/360 + ...
 947     //  we define R1(r*r) by
 948     //      r*(exp(r)+1)/(exp(r)-1) = 2 + r^2/6 * R1(r*r)
 949     //  That is,
 950     //      R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
 951     //               = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
 952     //               = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
 953     //  We use a special Remes algorithm on [0, 0.347] to generate
 954     //  a polynomial of degree 5 in r*r to approximate R1. The
 955     //  maximum error of this polynomial approximation is bounded
 956     //  by 2**-61. In other words,
 957     //      R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
 958     //  where   Q1  =  -1.6666666666666567384E-2,
 959     //          Q2  =   3.9682539681370365873E-4,
 960     //          Q3  =  -9.9206344733435987357E-6,
 961     //          Q4  =   2.5051361420808517002E-7,
 962     //          Q5  =  -6.2843505682382617102E-9;
 963     //          (where z=r*r, and Q1 to Q5 are called EXPM1_Qx in the source)
 964     //  with error bounded by
 965     //      |                  5           |     -61
 966     //      | 1.0+Q1*z+...+Q5*z   -  R1(z) | <= 2
 967     //      |                              |
 968     //
 969     //  expm1(r) = exp(r)-1 is then computed by the following
 970     //  specific way which minimize the accumulation rounding error:
 971     //                         2     3
 972     //                        r     r    [ 3 - (R1 + R1*r/2)  ]
 973     //        expm1(r) = r + --- + --- * [--------------------]
 974     //                        2     2    [ 6 - r*(3 - R1*r/2) ]
 975     //
 976     //  To compensate the error in the argument reduction, we use
 977     //          expm1(r+c) = expm1(r) + c + expm1(r)*c
 978     //                     ~ expm1(r) + c + r*c
 979     //  Thus c+r*c will be added in as the correction terms for
 980     //  expm1(r+c). Now rearrange the term to avoid optimization
 981     //  screw up:
 982     //                  (      2                                    2 )
 983     //                  ({  ( r    [ R1 -  (3 - R1*r/2) ]  )  }    r  )
 984     //   expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
 985     //                  ({  ( 2    [ 6 - r*(3 - R1*r/2) ]  )  }    2  )
 986     //                      (                                             )
 987     //
 988     //             = r - E
 989     //   3. Scale back to obtain expm1(x):
 990     //  From step 1, we have
 991     //     expm1(x) = either 2^k*[expm1(r)+1] - 1
 992     //              = or     2^k*[expm1(r) + (1-2^-k)]
 993     //   4. Implementation notes:
 994     //  (A). To save one multiplication, we scale the coefficient Qi
 995     //       to Qi*2^i, and replace z by (x^2)/2.
 996     //  (B). To achieve maximum accuracy, we compute expm1(x) by
 997     //    (i)   if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
 998     //    (ii)  if k=0, return r-E
 999     //    (iii) if k=-1, return 0.5*(r-E)-0.5
1000     //        (iv)      if k=1 if r < -0.25, return 2*((r+0.5)- E)
1001     //                 else          return  1.0+2.0*(r-E);
1002     //    (v)   if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
1003     //    (vi)  if k <= 20, return 2^k((1-2^-k)-(E-r)), else
1004     //    (vii) return 2^k(1-((E+2^-k)-r))
1005
1006     boolean negative = (x < 0);
1007     double y, hi, lo, c, t, e, hxs, hfx, r1;
1008     int k;
1009
1010     long bits;
1011     int h_bits;
1012     int l_bits;
1013
1014     c = 0.0;
1015     y = abs(x);
1016
1017     bits = Double.doubleToLongBits(y);
1018     h_bits = getHighDWord(bits);
1019     l_bits = getLowDWord(bits);
1020
1021     // handle special cases and large arguments
1022     if (h_bits >= 0x4043687a)        // if |x| >= 56 * ln(2)
1023       {
1024         if (h_bits >= 0x40862e42)    // if |x| >= EXP_LIMIT_H
1025           {
1026             if (h_bits >= 0x7ff00000)
1027               {
1028                 if (((h_bits & 0x000fffff) | (l_bits & 0xffffffff)) != 0)
1029                   return Double.NaN;               // exp(NaN) = NaN
1030                 else
1031                   return negative ? -1.0 : x;      // exp({+-inf}) = {+inf, -1}
1032               }
1033
1034             if (x > EXP_LIMIT_H)
1035               return Double.POSITIVE_INFINITY;     // overflow
1036           }
1037
1038         if (negative)                // x <= -56 * ln(2)
1039           return -1.0;
1040       }
1041
1042     // argument reduction
1043     if (h_bits > 0x3fd62e42)         // |x| > 0.5 * ln(2)
1044       {
1045         if (h_bits < 0x3ff0a2b2)     // |x| < 1.5 * ln(2)
1046           {
1047             if (negative)
1048               {
1049                 hi = x + LN2_H;
1050                 lo = -LN2_L;
1051                 k = -1;
1052               }
1053             else
1054               {
1055                 hi = x - LN2_H;
1056                 lo = LN2_L;
1057                 k  = 1;
1058               }
1059           }
1060         else
1061           {
1062             k = (int) (INV_LN2 * x + (negative ? - 0.5 : 0.5));
1063             t = k;
1064             hi = x - t * LN2_H;
1065             lo = t * LN2_L;
1066           }
1067
1068         x = hi - lo;
1069         c = (hi - x) - lo;
1070
1071       }
1072     else if (h_bits < 0x3c900000)    // |x| < 2^-54 return x
1073       return x;
1074     else
1075       k = 0;
1076
1077     // x is now in primary range
1078     hfx = 0.5 * x;
1079     hxs = x * hfx;
1080     r1 = 1.0 + hxs * (EXPM1_Q1
1081              + hxs * (EXPM1_Q2
1082              + hxs * (EXPM1_Q3
1083              + hxs * (EXPM1_Q4
1084              + hxs *  EXPM1_Q5))));
1085     t = 3.0 - r1 * hfx;
1086     e = hxs * ((r1 - t) / (6.0 - x * t));
1087
1088     if (k == 0)
1089       {
1090         return x - (x * e - hxs);    // c == 0
1091       }
1092     else
1093       {
1094         e = x * (e - c) - c;
1095         e -= hxs;
1096
1097         if (k == -1)
1098           return 0.5 * (x - e) - 0.5;
1099
1100         if (k == 1)
1101           {
1102             if (x < - 0.25)
1103               return -2.0 * (e - (x + 0.5));
1104             else
1105               return 1.0 + 2.0 * (x - e);
1106           }
1107
1108         if (k <= -2 || k > 56)       // sufficient to return exp(x) - 1
1109           {
1110             y = 1.0 - (e - x);
1111
1112             bits = Double.doubleToLongBits(y);
1113             h_bits = getHighDWord(bits);
1114             l_bits = getLowDWord(bits);
1115
1116             h_bits += (k << 20);     // add k to y's exponent
1117
1118             y = buildDouble(l_bits, h_bits);
1119
1120             return y - 1.0;
1121           }
1122
1123         t = 1.0;
1124         if (k < 20)
1125           {
1126             bits = Double.doubleToLongBits(t);
1127             h_bits = 0x3ff00000 - (0x00200000 >> k);
1128             l_bits = getLowDWord(bits);
1129
1130             t = buildDouble(l_bits, h_bits);      // t = 1 - 2^(-k)
1131             y = t - (e - x);
1132
1133             bits = Double.doubleToLongBits(y);
1134             h_bits = getHighDWord(bits);
1135             l_bits = getLowDWord(bits);
1136
1137             h_bits += (k << 20);     // add k to y's exponent
1138
1139             y = buildDouble(l_bits, h_bits);
1140           }
1141         else
1142           {
1143             bits = Double.doubleToLongBits(t);
1144             h_bits = (0x000003ff - k) << 20;
1145             l_bits = getLowDWord(bits);
1146
1147             t = buildDouble(l_bits, h_bits);      // t = 2^(-k)
1148
1149             y = x - (e + t);
1150             y += 1.0;
1151
1152             bits = Double.doubleToLongBits(y);
1153             h_bits = getHighDWord(bits);
1154             l_bits = getLowDWord(bits);
1155
1156             h_bits += (k << 20);     // add k to y's exponent
1157
1158             y = buildDouble(l_bits, h_bits);
1159           }
1160       }
1161
1162     return y;
1163   }
1164
1165   /**
1166    * Take ln(a) (the natural log).  The opposite of <code>exp()</code>. If the
1167    * argument is NaN or negative, the result is NaN; if the argument is
1168    * positive infinity, the result is positive infinity; and if the argument
1169    * is either zero, the result is negative infinity.
1170    *
1171    * <p>Note that the way to get log<sub>b</sub>(a) is to do this:
1172    * <code>ln(a) / ln(b)</code>.
1173    *
1174    * @param x the number to take the natural log of
1175    * @return the natural log of <code>a</code>
1176    * @see #exp(double)
1177    */
1178   public static double log(double x)
1179   {
1180     if (x == 0)
1181       return Double.NEGATIVE_INFINITY;
1182     if (x < 0)
1183       return Double.NaN;
1184     if (! (x < Double.POSITIVE_INFINITY))
1185       return x;
1186
1187     // Normalize x.
1188     long bits = Double.doubleToLongBits(x);
1189     int exp = (int) (bits >> 52);
1190     if (exp == 0) // Subnormal x.
1191       {
1192         x *= TWO_54;
1193         bits = Double.doubleToLongBits(x);
1194         exp = (int) (bits >> 52) - 54;
1195       }
1196     exp -= 1023; // Unbias exponent.
1197     bits = (bits & 0x000fffffffffffffL) | 0x3ff0000000000000L;
1198     x = Double.longBitsToDouble(bits);
1199     if (x >= SQRT_2)
1200       {
1201         x *= 0.5;
1202         exp++;
1203       }
1204     x--;
1205     if (abs(x) < 1 / TWO_20)
1206       {
1207         if (x == 0)
1208           return exp * LN2_H + exp * LN2_L;
1209         double r = x * x * (0.5 - 1 / 3.0 * x);
1210         if (exp == 0)
1211           return x - r;
1212         return exp * LN2_H - ((r - exp * LN2_L) - x);
1213       }
1214     double s = x / (2 + x);
1215     double z = s * s;
1216     double w = z * z;
1217     double t1 = w * (LG2 + w * (LG4 + w * LG6));
1218     double t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
1219     double r = t2 + t1;
1220     if (bits >= 0x3ff6174a00000000L && bits < 0x3ff6b85200000000L)
1221       {
1222         double h = 0.5 * x * x; // Need more accuracy for x near sqrt(2).
1223         if (exp == 0)
1224           return x - (h - s * (h + r));
1225         return exp * LN2_H - ((h - (s * (h + r) + exp * LN2_L)) - x);
1226       }
1227     if (exp == 0)
1228       return x - s * (x - r);
1229     return exp * LN2_H - ((s * (x - r) - exp * LN2_L) - x);
1230   }
1231
1232   /**
1233    * Take a square root. If the argument is NaN or negative, the result is
1234    * NaN; if the argument is positive infinity, the result is positive
1235    * infinity; and if the result is either zero, the result is the same.
1236    *
1237    * <p>For other roots, use pow(x, 1/rootNumber).
1238    *
1239    * @param x the numeric argument
1240    * @return the square root of the argument
1241    * @see #pow(double, double)
1242    */
1243   public static double sqrt(double x)
1244   {
1245     if (x < 0)
1246       return Double.NaN;
1247     if (x == 0 || ! (x < Double.POSITIVE_INFINITY))
1248       return x;
1249
1250     // Normalize x.
1251     long bits = Double.doubleToLongBits(x);
1252     int exp = (int) (bits >> 52);
1253     if (exp == 0) // Subnormal x.
1254       {
1255         x *= TWO_54;
1256         bits = Double.doubleToLongBits(x);
1257         exp = (int) (bits >> 52) - 54;
1258       }
1259     exp -= 1023; // Unbias exponent.
1260     bits = (bits & 0x000fffffffffffffL) | 0x0010000000000000L;
1261     if ((exp & 1) == 1) // Odd exp, double x to make it even.
1262       bits <<= 1;
1263     exp >>= 1;
1264
1265     // Generate sqrt(x) bit by bit.
1266     bits <<= 1;
1267     long q = 0;
1268     long s = 0;
1269     long r = 0x0020000000000000L; // Move r right to left.
1270     while (r != 0)
1271       {
1272         long t = s + r;
1273         if (t <= bits)
1274           {
1275             s = t + r;
1276             bits -= t;
1277             q += r;
1278           }
1279         bits <<= 1;
1280         r >>= 1;
1281       }
1282
1283     // Use floating add to round correctly.
1284     if (bits != 0)
1285       q += q & 1;
1286     return Double.longBitsToDouble((q >> 1) + ((exp + 1022L) << 52));
1287   }
1288
1289   /**
1290    * Raise a number to a power. Special cases:<ul>
1291    * <li>If the second argument is positive or negative zero, then the result
1292    * is 1.0.</li>
1293    * <li>If the second argument is 1.0, then the result is the same as the
1294    * first argument.</li>
1295    * <li>If the second argument is NaN, then the result is NaN.</li>
1296    * <li>If the first argument is NaN and the second argument is nonzero,
1297    * then the result is NaN.</li>
1298    * <li>If the absolute value of the first argument is greater than 1 and
1299    * the second argument is positive infinity, or the absolute value of the
1300    * first argument is less than 1 and the second argument is negative
1301    * infinity, then the result is positive infinity.</li>
1302    * <li>If the absolute value of the first argument is greater than 1 and
1303    * the second argument is negative infinity, or the absolute value of the
1304    * first argument is less than 1 and the second argument is positive
1305    * infinity, then the result is positive zero.</li>
1306    * <li>If the absolute value of the first argument equals 1 and the second
1307    * argument is infinite, then the result is NaN.</li>
1308    * <li>If the first argument is positive zero and the second argument is
1309    * greater than zero, or the first argument is positive infinity and the
1310    * second argument is less than zero, then the result is positive zero.</li>
1311    * <li>If the first argument is positive zero and the second argument is
1312    * less than zero, or the first argument is positive infinity and the
1313    * second argument is greater than zero, then the result is positive
1314    * infinity.</li>
1315    * <li>If the first argument is negative zero and the second argument is
1316    * greater than zero but not a finite odd integer, or the first argument is
1317    * negative infinity and the second argument is less than zero but not a
1318    * finite odd integer, then the result is positive zero.</li>
1319    * <li>If the first argument is negative zero and the second argument is a
1320    * positive finite odd integer, or the first argument is negative infinity
1321    * and the second argument is a negative finite odd integer, then the result
1322    * is negative zero.</li>
1323    * <li>If the first argument is negative zero and the second argument is
1324    * less than zero but not a finite odd integer, or the first argument is
1325    * negative infinity and the second argument is greater than zero but not a
1326    * finite odd integer, then the result is positive infinity.</li>
1327    * <li>If the first argument is negative zero and the second argument is a
1328    * negative finite odd integer, or the first argument is negative infinity
1329    * and the second argument is a positive finite odd integer, then the result
1330    * is negative infinity.</li>
1331    * <li>If the first argument is less than zero and the second argument is a
1332    * finite even integer, then the result is equal to the result of raising
1333    * the absolute value of the first argument to the power of the second
1334    * argument.</li>
1335    * <li>If the first argument is less than zero and the second argument is a
1336    * finite odd integer, then the result is equal to the negative of the
1337    * result of raising the absolute value of the first argument to the power
1338    * of the second argument.</li>
1339    * <li>If the first argument is finite and less than zero and the second
1340    * argument is finite and not an integer, then the result is NaN.</li>
1341    * <li>If both arguments are integers, then the result is exactly equal to
1342    * the mathematical result of raising the first argument to the power of
1343    * the second argument if that result can in fact be represented exactly as
1344    * a double value.</li>
1345    *
1346    * </ul><p>(In the foregoing descriptions, a floating-point value is
1347    * considered to be an integer if and only if it is a fixed point of the
1348    * method {@link #ceil(double)} or, equivalently, a fixed point of the
1349    * method {@link #floor(double)}. A value is a fixed point of a one-argument
1350    * method if and only if the result of applying the method to the value is
1351    * equal to the value.)
1352    *
1353    * @param x the number to raise
1354    * @param y the power to raise it to
1355    * @return x<sup>y</sup>
1356    */
1357   public static double pow(double x, double y)
1358   {
1359     // Special cases first.
1360     if (y == 0)
1361       return 1;
1362     if (y == 1)
1363       return x;
1364     if (y == -1)
1365       return 1 / x;
1366     if (x != x || y != y)
1367       return Double.NaN;
1368
1369     // When x < 0, yisint tells if y is not an integer (0), even(1),
1370     // or odd (2).
1371     int yisint = 0;
1372     if (x < 0 && floor(y) == y)
1373       yisint = (y % 2 == 0) ? 2 : 1;
1374     double ax = abs(x);
1375     double ay = abs(y);
1376
1377     // More special cases, of y.
1378     if (ay == Double.POSITIVE_INFINITY)
1379       {
1380         if (ax == 1)
1381           return Double.NaN;
1382         if (ax > 1)
1383           return y > 0 ? y : 0;
1384         return y < 0 ? -y : 0;
1385       }
1386     if (y == 2)
1387       return x * x;
1388     if (y == 0.5)
1389       return sqrt(x);
1390
1391     // More special cases, of x.
1392     if (x == 0 || ax == Double.POSITIVE_INFINITY || ax == 1)
1393       {
1394         if (y < 0)
1395           ax = 1 / ax;
1396         if (x < 0)
1397           {
1398             if (x == -1 && yisint == 0)
1399               ax = Double.NaN;
1400             else if (yisint == 1)
1401               ax = -ax;
1402           }
1403         return ax;
1404       }
1405     if (x < 0 && yisint == 0)
1406       return Double.NaN;
1407
1408     // Now we can start!
1409     double t;
1410     double t1;
1411     double t2;
1412     double u;
1413     double v;
1414     double w;
1415     if (ay > TWO_31)
1416       {
1417         if (ay > TWO_64) // Automatic over/underflow.
1418           return ((ax < 1) ? y < 0 : y > 0) ? Double.POSITIVE_INFINITY : 0;
1419         // Over/underflow if x is not close to one.
1420         if (ax < 0.9999995231628418)
1421           return y < 0 ? Double.POSITIVE_INFINITY : 0;
1422         if (ax >= 1.0000009536743164)
1423           return y > 0 ? Double.POSITIVE_INFINITY : 0;
1424         // Now |1-x| is <= 2**-20, sufficient to compute
1425         // log(x) by x-x^2/2+x^3/3-x^4/4.
1426         t = x - 1;
1427         w = t * t * (0.5 - t * (1 / 3.0 - t * 0.25));
1428         u = INV_LN2_H * t;
1429         v = t * INV_LN2_L - w * INV_LN2;
1430         t1 = (float) (u + v);
1431         t2 = v - (t1 - u);
1432       }
1433     else
1434     {
1435       long bits = Double.doubleToLongBits(ax);
1436       int exp = (int) (bits >> 52);
1437       if (exp == 0) // Subnormal x.
1438         {
1439           ax *= TWO_54;
1440           bits = Double.doubleToLongBits(ax);
1441           exp = (int) (bits >> 52) - 54;
1442         }
1443       exp -= 1023; // Unbias exponent.
1444       ax = Double.longBitsToDouble((bits & 0x000fffffffffffffL)
1445                                    | 0x3ff0000000000000L);
1446       boolean k;
1447       if (ax < SQRT_1_5)  // |x|<sqrt(3/2).
1448         k = false;
1449       else if (ax < SQRT_3) // |x|<sqrt(3).
1450         k = true;
1451       else
1452         {
1453           k = false;
1454           ax *= 0.5;
1455           exp++;
1456         }
1457
1458       // Compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5).
1459       u = ax - (k ? 1.5 : 1);
1460       v = 1 / (ax + (k ? 1.5 : 1));
1461       double s = u * v;
1462       double s_h = (float) s;
1463       double t_h = (float) (ax + (k ? 1.5 : 1));
1464       double t_l = ax - (t_h - (k ? 1.5 : 1));
1465       double s_l = v * ((u - s_h * t_h) - s_h * t_l);
1466       // Compute log(ax).
1467       double s2 = s * s;
1468       double r = s_l * (s_h + s) + s2 * s2
1469         * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6)))));
1470       s2 = s_h * s_h;
1471       t_h = (float) (3.0 + s2 + r);
1472       t_l = r - (t_h - 3.0 - s2);
1473       // u+v = s*(1+...).
1474       u = s_h * t_h;
1475       v = s_l * t_h + t_l * s;
1476       // 2/(3log2)*(s+...).
1477       double p_h = (float) (u + v);
1478       double p_l = v - (p_h - u);
1479       double z_h = CP_H * p_h;
1480       double z_l = CP_L * p_h + p_l * CP + (k ? DP_L : 0);
1481       // log2(ax) = (s+..)*2/(3*log2) = exp + dp_h + z_h + z_l.
1482       t = exp;
1483       t1 = (float) (z_h + z_l + (k ? DP_H : 0) + t);
1484       t2 = z_l - (t1 - t - (k ? DP_H : 0) - z_h);
1485     }
1486
1487     // Split up y into y1+y2 and compute (y1+y2)*(t1+t2).
1488     boolean negative = x < 0 && yisint == 1;
1489     double y1 = (float) y;
1490     double p_l = (y - y1) * t1 + y * t2;
1491     double p_h = y1 * t1;
1492     double z = p_l + p_h;
1493     if (z >= 1024) // Detect overflow.
1494       {
1495         if (z > 1024 || p_l + OVT > z - p_h)
1496           return negative ? Double.NEGATIVE_INFINITY
1497             : Double.POSITIVE_INFINITY;
1498       }
1499     else if (z <= -1075) // Detect underflow.
1500       {
1501         if (z < -1075 || p_l <= z - p_h)
1502           return negative ? -0.0 : 0;
1503       }
1504
1505     // Compute 2**(p_h+p_l).
1506     int n = round((float) z);
1507     p_h -= n;
1508     t = (float) (p_l + p_h);
1509     u = t * LN2_H;
1510     v = (p_l - (t - p_h)) * LN2 + t * LN2_L;
1511     z = u + v;
1512     w = v - (z - u);
1513     t = z * z;
1514     t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
1515     double r = (z * t1) / (t1 - 2) - (w + z * w);
1516     z = scale(1 - (r - z), n);
1517     return negative ? -z : z;
1518   }
1519
1520   /**
1521    * Get the IEEE 754 floating point remainder on two numbers. This is the
1522    * value of <code>x - y * <em>n</em></code>, where <em>n</em> is the closest
1523    * double to <code>x / y</code> (ties go to the even n); for a zero
1524    * remainder, the sign is that of <code>x</code>. If either argument is NaN,
1525    * the first argument is infinite, or the second argument is zero, the result
1526    * is NaN; if x is finite but y is infinite, the result is x.
1527    *
1528    * @param x the dividend (the top half)
1529    * @param y the divisor (the bottom half)
1530    * @return the IEEE 754-defined floating point remainder of x/y
1531    * @see #rint(double)
1532    */
1533   public static double IEEEremainder(double x, double y)
1534   {
1535     // Purge off exception values.
1536     if (x == Double.NEGATIVE_INFINITY || ! (x < Double.POSITIVE_INFINITY)
1537         || y == 0 || y != y)
1538       return Double.NaN;
1539
1540     boolean negative = x < 0;
1541     x = abs(x);
1542     y = abs(y);
1543     if (x == y || x == 0)
1544       return 0 * x; // Get correct sign.
1545
1546     // Achieve x < 2y, then take first shot at remainder.
1547     if (y < TWO_1023)
1548       x %= y + y;
1549
1550     // Now adjust x to get correct precision.
1551     if (y < 4 / TWO_1023)
1552       {
1553         if (x + x > y)
1554           {
1555             x -= y;
1556             if (x + x >= y)
1557               x -= y;
1558           }
1559       }
1560     else
1561       {
1562         y *= 0.5;
1563         if (x > y)
1564           {
1565             x -= y;
1566             if (x >= y)
1567               x -= y;
1568           }
1569       }
1570     return negative ? -x : x;
1571   }
1572
1573   /**
1574    * Take the nearest integer that is that is greater than or equal to the
1575    * argument. If the argument is NaN, infinite, or zero, the result is the
1576    * same; if the argument is between -1 and 0, the result is negative zero.
1577    * Note that <code>Math.ceil(x) == -Math.floor(-x)</code>.
1578    *
1579    * @param a the value to act upon
1580    * @return the nearest integer &gt;= <code>a</code>
1581    */
1582   public static double ceil(double a)
1583   {
1584     return -floor(-a);
1585   }
1586
1587   /**
1588    * Take the nearest integer that is that is less than or equal to the
1589    * argument. If the argument is NaN, infinite, or zero, the result is the
1590    * same. Note that <code>Math.ceil(x) == -Math.floor(-x)</code>.
1591    *
1592    * @param a the value to act upon
1593    * @return the nearest integer &lt;= <code>a</code>
1594    */
1595   public static double floor(double a)
1596   {
1597     double x = abs(a);
1598     if (! (x < TWO_52) || (long) a == a)
1599       return a; // No fraction bits; includes NaN and infinity.
1600     if (x < 1)
1601       return a >= 0 ? 0 * a : -1; // Worry about signed zero.
1602     return a < 0 ? (long) a - 1.0 : (long) a; // Cast to long truncates.
1603   }
1604
1605   /**
1606    * Take the nearest integer to the argument.  If it is exactly between
1607    * two integers, the even integer is taken. If the argument is NaN,
1608    * infinite, or zero, the result is the same.
1609    *
1610    * @param a the value to act upon
1611    * @return the nearest integer to <code>a</code>
1612    */
1613   public static double rint(double a)
1614   {
1615     double x = abs(a);
1616     if (! (x < TWO_52))
1617       return a; // No fraction bits; includes NaN and infinity.
1618     if (x <= 0.5)
1619       return 0 * a; // Worry about signed zero.
1620     if (x % 2 <= 0.5)
1621       return (long) a; // Catch round down to even.
1622     return (long) (a + (a < 0 ? -0.5 : 0.5)); // Cast to long truncates.
1623   }
1624
1625   /**
1626    * Take the nearest integer to the argument.  This is equivalent to
1627    * <code>(int) Math.floor(f + 0.5f)</code>. If the argument is NaN, the
1628    * result is 0; otherwise if the argument is outside the range of int, the
1629    * result will be Integer.MIN_VALUE or Integer.MAX_VALUE, as appropriate.
1630    *
1631    * @param f the argument to round
1632    * @return the nearest integer to the argument
1633    * @see Integer#MIN_VALUE
1634    * @see Integer#MAX_VALUE
1635    */
1636   public static int round(float f)
1637   {
1638     return (int) floor(f + 0.5f);
1639   }
1640
1641   /**
1642    * Take the nearest long to the argument.  This is equivalent to
1643    * <code>(long) Math.floor(d + 0.5)</code>. If the argument is NaN, the
1644    * result is 0; otherwise if the argument is outside the range of long, the
1645    * result will be Long.MIN_VALUE or Long.MAX_VALUE, as appropriate.
1646    *
1647    * @param d the argument to round
1648    * @return the nearest long to the argument
1649    * @see Long#MIN_VALUE
1650    * @see Long#MAX_VALUE
1651    */
1652   public static long round(double d)
1653   {
1654     return (long) floor(d + 0.5);
1655   }
1656
1657   /**
1658    * Get a random number.  This behaves like Random.nextDouble(), seeded by
1659    * System.currentTimeMillis() when first called. In other words, the number
1660    * is from a pseudorandom sequence, and lies in the range [+0.0, 1.0).
1661    * This random sequence is only used by this method, and is threadsafe,
1662    * although you may want your own random number generator if it is shared
1663    * among threads.
1664    *
1665    * @return a random number
1666    * @see Random#nextDouble()
1667    * @see System#currentTimeMillis()
1668    */
1669   public static synchronized double random()
1670   {
1671     if (rand == null)
1672       rand = new Random();
1673     return rand.nextDouble();
1674   }
1675
1676   /**
1677    * Convert from degrees to radians. The formula for this is
1678    * radians = degrees * (pi/180); however it is not always exact given the
1679    * limitations of floating point numbers.
1680    *
1681    * @param degrees an angle in degrees
1682    * @return the angle in radians
1683    */
1684   public static double toRadians(double degrees)
1685   {
1686     return (degrees * PI) / 180;
1687   }
1688
1689   /**
1690    * Convert from radians to degrees. The formula for this is
1691    * degrees = radians * (180/pi); however it is not always exact given the
1692    * limitations of floating point numbers.
1693    *
1694    * @param rads an angle in radians
1695    * @return the angle in degrees
1696    */
1697   public static double toDegrees(double rads)
1698   {
1699     return (rads * 180) / PI;
1700   }
1701
1702   /**
1703    * Constants for scaling and comparing doubles by powers of 2. The compiler
1704    * must automatically inline constructs like (1/TWO_54), so we don't list
1705    * negative powers of two here.
1706    */
1707   private static final double
1708     TWO_16 = 0x10000, // Long bits 0x40f0000000000000L.
1709     TWO_20 = 0x100000, // Long bits 0x4130000000000000L.
1710     TWO_24 = 0x1000000, // Long bits 0x4170000000000000L.
1711     TWO_27 = 0x8000000, // Long bits 0x41a0000000000000L.
1712     TWO_28 = 0x10000000, // Long bits 0x41b0000000000000L.
1713     TWO_29 = 0x20000000, // Long bits 0x41c0000000000000L.
1714     TWO_31 = 0x80000000L, // Long bits 0x41e0000000000000L.
1715     TWO_49 = 0x2000000000000L, // Long bits 0x4300000000000000L.
1716     TWO_52 = 0x10000000000000L, // Long bits 0x4330000000000000L.
1717     TWO_54 = 0x40000000000000L, // Long bits 0x4350000000000000L.
1718     TWO_57 = 0x200000000000000L, // Long bits 0x4380000000000000L.
1719     TWO_60 = 0x1000000000000000L, // Long bits 0x43b0000000000000L.
1720     TWO_64 = 1.8446744073709552e19, // Long bits 0x43f0000000000000L.
1721     TWO_66 = 7.378697629483821e19, // Long bits 0x4410000000000000L.
1722     TWO_1023 = 8.98846567431158e307; // Long bits 0x7fe0000000000000L.
1723
1724   /**
1725    * Super precision for 2/pi in 24-bit chunks, for use in
1726    * {@link #remPiOver2(double, double[])}.
1727    */
1728   private static final int TWO_OVER_PI[] = {
1729     0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
1730     0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
1731     0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
1732     0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
1733     0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
1734     0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
1735     0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
1736     0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
1737     0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
1738     0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
1739     0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
1740   };
1741
1742   /**
1743    * Super precision for pi/2 in 24-bit chunks, for use in
1744    * {@link #remPiOver2(double, double[])}.
1745    */
1746   private static final double PI_OVER_TWO[] = {
1747     1.570796251296997, // Long bits 0x3ff921fb40000000L.
1748     7.549789415861596e-8, // Long bits 0x3e74442d00000000L.
1749     5.390302529957765e-15, // Long bits 0x3cf8469880000000L.
1750     3.282003415807913e-22, // Long bits 0x3b78cc5160000000L.
1751     1.270655753080676e-29, // Long bits 0x39f01b8380000000L.
1752     1.2293330898111133e-36, // Long bits 0x387a252040000000L.
1753     2.7337005381646456e-44, // Long bits 0x36e3822280000000L.
1754     2.1674168387780482e-51, // Long bits 0x3569f31d00000000L.
1755   };
1756
1757   /**
1758    * More constants related to pi, used in
1759    * {@link #remPiOver2(double, double[])} and elsewhere.
1760    */
1761   private static final double
1762     PI_L = 1.2246467991473532e-16, // Long bits 0x3ca1a62633145c07L.
1763     PIO2_1 = 1.5707963267341256, // Long bits 0x3ff921fb54400000L.
1764     PIO2_1L = 6.077100506506192e-11, // Long bits 0x3dd0b4611a626331L.
1765     PIO2_2 = 6.077100506303966e-11, // Long bits 0x3dd0b4611a600000L.
1766     PIO2_2L = 2.0222662487959506e-21, // Long bits 0x3ba3198a2e037073L.
1767     PIO2_3 = 2.0222662487111665e-21, // Long bits 0x3ba3198a2e000000L.
1768     PIO2_3L = 8.4784276603689e-32; // Long bits 0x397b839a252049c1L.
1769
1770   /**
1771    * Natural log and square root constants, for calculation of
1772    * {@link #exp(double)}, {@link #log(double)} and
1773    * {@link #pow(double, double)}. CP is 2/(3*ln(2)).
1774    */
1775   private static final double
1776     SQRT_1_5 = 1.224744871391589, // Long bits 0x3ff3988e1409212eL.
1777     SQRT_2 = 1.4142135623730951, // Long bits 0x3ff6a09e667f3bcdL.
1778     SQRT_3 = 1.7320508075688772, // Long bits 0x3ffbb67ae8584caaL.
1779     EXP_LIMIT_H = 709.782712893384, // Long bits 0x40862e42fefa39efL.
1780     EXP_LIMIT_L = -745.1332191019411, // Long bits 0xc0874910d52d3051L.
1781     CP = 0.9617966939259756, // Long bits 0x3feec709dc3a03fdL.
1782     CP_H = 0.9617967009544373, // Long bits 0x3feec709e0000000L.
1783     CP_L = -7.028461650952758e-9, // Long bits 0xbe3e2fe0145b01f5L.
1784     LN2 = 0.6931471805599453, // Long bits 0x3fe62e42fefa39efL.
1785     LN2_H = 0.6931471803691238, // Long bits 0x3fe62e42fee00000L.
1786     LN2_L = 1.9082149292705877e-10, // Long bits 0x3dea39ef35793c76L.
1787     INV_LN2 = 1.4426950408889634, // Long bits 0x3ff71547652b82feL.
1788     INV_LN2_H = 1.4426950216293335, // Long bits 0x3ff7154760000000L.
1789     INV_LN2_L = 1.9259629911266175e-8; // Long bits 0x3e54ae0bf85ddf44L.
1790
1791   /**
1792    * Constants for computing {@link #log(double)}.
1793    */
1794   private static final double
1795     LG1 = 0.6666666666666735, // Long bits 0x3fe5555555555593L.
1796     LG2 = 0.3999999999940942, // Long bits 0x3fd999999997fa04L.
1797     LG3 = 0.2857142874366239, // Long bits 0x3fd2492494229359L.
1798     LG4 = 0.22222198432149784, // Long bits 0x3fcc71c51d8e78afL.
1799     LG5 = 0.1818357216161805, // Long bits 0x3fc7466496cb03deL.
1800     LG6 = 0.15313837699209373, // Long bits 0x3fc39a09d078c69fL.
1801     LG7 = 0.14798198605116586; // Long bits 0x3fc2f112df3e5244L.
1802
1803   /**
1804    * Constants for computing {@link #pow(double, double)}. L and P are
1805    * coefficients for series; OVT is -(1024-log2(ovfl+.5ulp)); and DP is ???.
1806    * The P coefficients also calculate {@link #exp(double)}.
1807    */
1808   private static final double
1809     L1 = 0.5999999999999946, // Long bits 0x3fe3333333333303L.
1810     L2 = 0.4285714285785502, // Long bits 0x3fdb6db6db6fabffL.
1811     L3 = 0.33333332981837743, // Long bits 0x3fd55555518f264dL.
1812     L4 = 0.272728123808534, // Long bits 0x3fd17460a91d4101L.
1813     L5 = 0.23066074577556175, // Long bits 0x3fcd864a93c9db65L.
1814     L6 = 0.20697501780033842, // Long bits 0x3fca7e284a454eefL.
1815     P1 = 0.16666666666666602, // Long bits 0x3fc555555555553eL.
1816     P2 = -2.7777777777015593e-3, // Long bits 0xbf66c16c16bebd93L.
1817     P3 = 6.613756321437934e-5, // Long bits 0x3f11566aaf25de2cL.
1818     P4 = -1.6533902205465252e-6, // Long bits 0xbebbbd41c5d26bf1L.
1819     P5 = 4.1381367970572385e-8, // Long bits 0x3e66376972bea4d0L.
1820     DP_H = 0.5849624872207642, // Long bits 0x3fe2b80340000000L.
1821     DP_L = 1.350039202129749e-8, // Long bits 0x3e4cfdeb43cfd006L.
1822     OVT = 8.008566259537294e-17; // Long bits 0x3c971547652b82feL.
1823
1824   /**
1825    * Coefficients for computing {@link #sin(double)}.
1826    */
1827   private static final double
1828     S1 = -0.16666666666666632, // Long bits 0xbfc5555555555549L.
1829     S2 = 8.33333333332249e-3, // Long bits 0x3f8111111110f8a6L.
1830     S3 = -1.984126982985795e-4, // Long bits 0xbf2a01a019c161d5L.
1831     S4 = 2.7557313707070068e-6, // Long bits 0x3ec71de357b1fe7dL.
1832     S5 = -2.5050760253406863e-8, // Long bits 0xbe5ae5e68a2b9cebL.
1833     S6 = 1.58969099521155e-10; // Long bits 0x3de5d93a5acfd57cL.
1834
1835   /**
1836    * Coefficients for computing {@link #cos(double)}.
1837    */
1838   private static final double
1839     C1 = 0.0416666666666666, // Long bits 0x3fa555555555554cL.
1840     C2 = -1.388888888887411e-3, // Long bits 0xbf56c16c16c15177L.
1841     C3 = 2.480158728947673e-5, // Long bits 0x3efa01a019cb1590L.
1842     C4 = -2.7557314351390663e-7, // Long bits 0xbe927e4f809c52adL.
1843     C5 = 2.087572321298175e-9, // Long bits 0x3e21ee9ebdb4b1c4L.
1844     C6 = -1.1359647557788195e-11; // Long bits 0xbda8fae9be8838d4L.
1845
1846   /**
1847    * Coefficients for computing {@link #tan(double)}.
1848    */
1849   private static final double
1850     T0 = 0.3333333333333341, // Long bits 0x3fd5555555555563L.
1851     T1 = 0.13333333333320124, // Long bits 0x3fc111111110fe7aL.
1852     T2 = 0.05396825397622605, // Long bits 0x3faba1ba1bb341feL.
1853     T3 = 0.021869488294859542, // Long bits 0x3f9664f48406d637L.
1854     T4 = 8.8632398235993e-3, // Long bits 0x3f8226e3e96e8493L.
1855     T5 = 3.5920791075913124e-3, // Long bits 0x3f6d6d22c9560328L.
1856     T6 = 1.4562094543252903e-3, // Long bits 0x3f57dbc8fee08315L.
1857     T7 = 5.880412408202641e-4, // Long bits 0x3f4344d8f2f26501L.
1858     T8 = 2.464631348184699e-4, // Long bits 0x3f3026f71a8d1068L.
1859     T9 = 7.817944429395571e-5, // Long bits 0x3f147e88a03792a6L.
1860     T10 = 7.140724913826082e-5, // Long bits 0x3f12b80f32f0a7e9L.
1861     T11 = -1.8558637485527546e-5, // Long bits 0xbef375cbdb605373L.
1862     T12 = 2.590730518636337e-5; // Long bits 0x3efb2a7074bf7ad4L.
1863
1864   /**
1865    * Coefficients for computing {@link #asin(double)} and
1866    * {@link #acos(double)}.
1867    */
1868   private static final double
1869     PS0 = 0.16666666666666666, // Long bits 0x3fc5555555555555L.
1870     PS1 = -0.3255658186224009, // Long bits 0xbfd4d61203eb6f7dL.
1871     PS2 = 0.20121253213486293, // Long bits 0x3fc9c1550e884455L.
1872     PS3 = -0.04005553450067941, // Long bits 0xbfa48228b5688f3bL.
1873     PS4 = 7.915349942898145e-4, // Long bits 0x3f49efe07501b288L.
1874     PS5 = 3.479331075960212e-5, // Long bits 0x3f023de10dfdf709L.
1875     QS1 = -2.403394911734414, // Long bits 0xc0033a271c8a2d4bL.
1876     QS2 = 2.0209457602335057, // Long bits 0x40002ae59c598ac8L.
1877     QS3 = -0.6882839716054533, // Long bits 0xbfe6066c1b8d0159L.
1878     QS4 = 0.07703815055590194; // Long bits 0x3fb3b8c5b12e9282L.
1879
1880   /**
1881    * Coefficients for computing {@link #atan(double)}.
1882    */
1883   private static final double
1884     ATAN_0_5H = 0.4636476090008061, // Long bits 0x3fddac670561bb4fL.
1885     ATAN_0_5L = 2.2698777452961687e-17, // Long bits 0x3c7a2b7f222f65e2L.
1886     ATAN_1_5H = 0.982793723247329, // Long bits 0x3fef730bd281f69bL.
1887     ATAN_1_5L = 1.3903311031230998e-17, // Long bits 0x3c7007887af0cbbdL.
1888     AT0 = 0.3333333333333293, // Long bits 0x3fd555555555550dL.
1889     AT1 = -0.19999999999876483, // Long bits 0xbfc999999998ebc4L.
1890     AT2 = 0.14285714272503466, // Long bits 0x3fc24924920083ffL.
1891     AT3 = -0.11111110405462356, // Long bits 0xbfbc71c6fe231671L.
1892     AT4 = 0.09090887133436507, // Long bits 0x3fb745cdc54c206eL.
1893     AT5 = -0.0769187620504483, // Long bits 0xbfb3b0f2af749a6dL.
1894     AT6 = 0.06661073137387531, // Long bits 0x3fb10d66a0d03d51L.
1895     AT7 = -0.058335701337905735, // Long bits 0xbfadde2d52defd9aL.
1896     AT8 = 0.049768779946159324, // Long bits 0x3fa97b4b24760debL.
1897     AT9 = -0.036531572744216916, // Long bits 0xbfa2b4442c6a6c2fL.
1898     AT10 = 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L.
1899
1900   /**
1901    * Constants for computing {@link #cbrt(double)}.
1902    */
1903   private static final int
1904     CBRT_B1 = 715094163, // B1 = (682-0.03306235651)*2**20
1905     CBRT_B2 = 696219795; // B2 = (664-0.03306235651)*2**20
1906
1907   /**
1908    * Constants for computing {@link #cbrt(double)}.
1909    */
1910   private static final double
1911     CBRT_C =  5.42857142857142815906e-01, // Long bits  0x3fe15f15f15f15f1L
1912     CBRT_D = -7.05306122448979611050e-01, // Long bits  0xbfe691de2532c834L
1913     CBRT_E =  1.41428571428571436819e+00, // Long bits  0x3ff6a0ea0ea0ea0fL
1914     CBRT_F =  1.60714285714285720630e+00, // Long bits  0x3ff9b6db6db6db6eL
1915     CBRT_G =  3.57142857142857150787e-01; // Long bits  0x3fd6db6db6db6db7L
1916
1917   /**
1918    * Constants for computing {@link #expm1(double)}
1919    */
1920   private static final double
1921     EXPM1_Q1 = -3.33333333333331316428e-02, // Long bits  0xbfa11111111110f4L
1922     EXPM1_Q2 =  1.58730158725481460165e-03, // Long bits  0x3f5a01a019fe5585L
1923     EXPM1_Q3 = -7.93650757867487942473e-05, // Long bits  0xbf14ce199eaadbb7L
1924     EXPM1_Q4 =  4.00821782732936239552e-06, // Long bits  0x3ed0cfca86e65239L
1925     EXPM1_Q5 = -2.01099218183624371326e-07; // Long bits  0xbe8afdb76e09c32dL
1926
1927   /**
1928    * Helper function for reducing an angle to a multiple of pi/2 within
1929    * [-pi/4, pi/4].
1930    *
1931    * @param x the angle; not infinity or NaN, and outside pi/4
1932    * @param y an array of 2 doubles modified to hold the remander x % pi/2
1933    * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4],
1934    *         1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4]
1935    */
1936   private static int remPiOver2(double x, double[] y)
1937   {
1938     boolean negative = x < 0;
1939     x = abs(x);
1940     double z;
1941     int n;
1942     if (Configuration.DEBUG && (x <= PI / 4 || x != x
1943                                 || x == Double.POSITIVE_INFINITY))
1944       throw new InternalError("Assertion failure");
1945     if (x < 3 * PI / 4) // If |x| is small.
1946       {
1947         z = x - PIO2_1;
1948         if ((float) x != (float) (PI / 2)) // 33+53 bit pi is good enough.
1949           {
1950             y[0] = z - PIO2_1L;
1951             y[1] = z - y[0] - PIO2_1L;
1952           }
1953         else // Near pi/2, use 33+33+53 bit pi.
1954           {
1955             z -= PIO2_2;
1956             y[0] = z - PIO2_2L;
1957             y[1] = z - y[0] - PIO2_2L;
1958           }
1959         n = 1;
1960       }
1961     else if (x <= TWO_20 * PI / 2) // Medium size.
1962       {
1963         n = (int) (2 / PI * x + 0.5);
1964         z = x - n * PIO2_1;
1965         double w = n * PIO2_1L; // First round good to 85 bits.
1966         y[0] = z - w;
1967         if (n >= 32 || (float) x == (float) (w))
1968           {
1969             if (x / y[0] >= TWO_16) // Second iteration, good to 118 bits.
1970               {
1971                 double t = z;
1972                 w = n * PIO2_2;
1973                 z = t - w;
1974                 w = n * PIO2_2L - (t - z - w);
1975                 y[0] = z - w;
1976                 if (x / y[0] >= TWO_49) // Third iteration, 151 bits accuracy.
1977                   {
1978                     t = z;
1979                     w = n * PIO2_3;
1980                     z = t - w;
1981                     w = n * PIO2_3L - (t - z - w);
1982                     y[0] = z - w;
1983                   }
1984               }
1985           }
1986         y[1] = z - y[0] - w;
1987       }
1988     else
1989       {
1990         // All other (large) arguments.
1991         int e0 = (int) (Double.doubleToLongBits(x) >> 52) - 1046;
1992         z = scale(x, -e0); // e0 = ilogb(z) - 23.
1993         double[] tx = new double[3];
1994         for (int i = 0; i < 2; i++)
1995           {
1996             tx[i] = (int) z;
1997             z = (z - tx[i]) * TWO_24;
1998           }
1999         tx[2] = z;
2000         int nx = 2;
2001         while (tx[nx] == 0)
2002           nx--;
2003         n = remPiOver2(tx, y, e0, nx);
2004       }
2005     if (negative)
2006       {
2007         y[0] = -y[0];
2008         y[1] = -y[1];
2009         return -n;
2010       }
2011     return n;
2012   }
2013
2014   /**
2015    * Helper function for reducing an angle to a multiple of pi/2 within
2016    * [-pi/4, pi/4].
2017    *
2018    * @param x the positive angle, broken into 24-bit chunks
2019    * @param y an array of 2 doubles modified to hold the remander x % pi/2
2020    * @param e0 the exponent of x[0]
2021    * @param nx the last index used in x
2022    * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4],
2023    *         1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4]
2024    */
2025   private static int remPiOver2(double[] x, double[] y, int e0, int nx)
2026   {
2027     int i;
2028     int ih;
2029     int n;
2030     double fw;
2031     double z;
2032     int[] iq = new int[20];
2033     double[] f = new double[20];
2034     double[] q = new double[20];
2035     boolean recompute = false;
2036
2037     // Initialize jk, jz, jv, q0; note that 3>q0.
2038     int jk = 4;
2039     int jz = jk;
2040     int jv = max((e0 - 3) / 24, 0);
2041     int q0 = e0 - 24 * (jv + 1);
2042
2043     // Set up f[0] to f[nx+jk] where f[nx+jk] = TWO_OVER_PI[jv+jk].
2044     int j = jv - nx;
2045     int m = nx + jk;
2046     for (i = 0; i <= m; i++, j++)
2047       f[i] = (j < 0) ? 0 : TWO_OVER_PI[j];
2048
2049     // Compute q[0],q[1],...q[jk].
2050     for (i = 0; i <= jk; i++)
2051       {
2052         for (j = 0, fw = 0; j <= nx; j++)
2053           fw += x[j] * f[nx + i - j];
2054         q[i] = fw;
2055       }
2056
2057     do
2058       {
2059         // Distill q[] into iq[] reversingly.
2060         for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--)
2061           {
2062             fw = (int) (1 / TWO_24 * z);
2063             iq[i] = (int) (z - TWO_24 * fw);
2064             z = q[j - 1] + fw;
2065           }
2066
2067         // Compute n.
2068         z = scale(z, q0);
2069         z -= 8 * floor(z * 0.125); // Trim off integer >= 8.
2070         n = (int) z;
2071         z -= n;
2072         ih = 0;
2073         if (q0 > 0) // Need iq[jz-1] to determine n.
2074           {
2075             i = iq[jz - 1] >> (24 - q0);
2076             n += i;
2077             iq[jz - 1] -= i << (24 - q0);
2078             ih = iq[jz - 1] >> (23 - q0);
2079           }
2080         else if (q0 == 0)
2081           ih = iq[jz - 1] >> 23;
2082         else if (z >= 0.5)
2083           ih = 2;
2084
2085         if (ih > 0) // If q > 0.5.
2086           {
2087             n += 1;
2088             int carry = 0;
2089             for (i = 0; i < jz; i++) // Compute 1-q.
2090               {
2091                 j = iq[i];
2092                 if (carry == 0)
2093                   {
2094                     if (j != 0)
2095                       {
2096                         carry = 1;
2097                         iq[i] = 0x1000000 - j;
2098                       }
2099                   }
2100                 else
2101                   iq[i] = 0xffffff - j;
2102               }
2103             switch (q0)
2104               {
2105               case 1: // Rare case: chance is 1 in 12 for non-default.
2106                 iq[jz - 1] &= 0x7fffff;
2107                 break;
2108               case 2:
2109                 iq[jz - 1] &= 0x3fffff;
2110               }
2111             if (ih == 2)
2112               {
2113                 z = 1 - z;
2114                 if (carry != 0)
2115                   z -= scale(1, q0);
2116               }
2117           }
2118
2119         // Check if recomputation is needed.
2120         if (z == 0)
2121           {
2122             j = 0;
2123             for (i = jz - 1; i >= jk; i--)
2124               j |= iq[i];
2125             if (j == 0) // Need recomputation.
2126               {
2127                 int k;
2128                 for (k = 1; iq[jk - k] == 0; k++); // k = no. of terms needed.
2129
2130                 for (i = jz + 1; i <= jz + k; i++) // Add q[jz+1] to q[jz+k].
2131                   {
2132                     f[nx + i] = TWO_OVER_PI[jv + i];
2133                     for (j = 0, fw = 0; j <= nx; j++)
2134                       fw += x[j] * f[nx + i - j];
2135                     q[i] = fw;
2136                   }
2137                 jz += k;
2138                 recompute = true;
2139               }
2140           }
2141       }
2142     while (recompute);
2143
2144     // Chop off zero terms.
2145     if (z == 0)
2146       {
2147         jz--;
2148         q0 -= 24;
2149         while (iq[jz] == 0)
2150           {
2151             jz--;
2152             q0 -= 24;
2153           }
2154       }
2155     else // Break z into 24-bit if necessary.
2156       {
2157         z = scale(z, -q0);
2158         if (z >= TWO_24)
2159           {
2160             fw = (int) (1 / TWO_24 * z);
2161             iq[jz] = (int) (z - TWO_24 * fw);
2162             jz++;
2163             q0 += 24;
2164             iq[jz] = (int) fw;
2165           }
2166         else
2167           iq[jz] = (int) z;
2168       }
2169
2170     // Convert integer "bit" chunk to floating-point value.
2171     fw = scale(1, q0);
2172     for (i = jz; i >= 0; i--)
2173       {
2174         q[i] = fw * iq[i];
2175         fw *= 1 / TWO_24;
2176       }
2177
2178     // Compute PI_OVER_TWO[0,...,jk]*q[jz,...,0].
2179     double[] fq = new double[20];
2180     for (i = jz; i >= 0; i--)
2181       {
2182         fw = 0;
2183         for (int k = 0; k <= jk && k <= jz - i; k++)
2184           fw += PI_OVER_TWO[k] * q[i + k];
2185         fq[jz - i] = fw;
2186       }
2187
2188     // Compress fq[] into y[].
2189     fw = 0;
2190     for (i = jz; i >= 0; i--)
2191       fw += fq[i];
2192     y[0] = (ih == 0) ? fw : -fw;
2193     fw = fq[0] - fw;
2194     for (i = 1; i <= jz; i++)
2195       fw += fq[i];
2196     y[1] = (ih == 0) ? fw : -fw;
2197     return n;
2198   }
2199
2200   /**
2201    * Helper method for scaling a double by a power of 2.
2202    *
2203    * @param x the double
2204    * @param n the scale; |n| < 2048
2205    * @return x * 2**n
2206    */
2207   private static double scale(double x, int n)
2208   {
2209     if (Configuration.DEBUG && abs(n) >= 2048)
2210       throw new InternalError("Assertion failure");
2211     if (x == 0 || x == Double.NEGATIVE_INFINITY
2212         || ! (x < Double.POSITIVE_INFINITY) || n == 0)
2213       return x;
2214     long bits = Double.doubleToLongBits(x);
2215     int exp = (int) (bits >> 52) & 0x7ff;
2216     if (exp == 0) // Subnormal x.
2217       {
2218         x *= TWO_54;
2219         exp = ((int) (Double.doubleToLongBits(x) >> 52) & 0x7ff) - 54;
2220       }
2221     exp += n;
2222     if (exp > 0x7fe) // Overflow.
2223       return Double.POSITIVE_INFINITY * x;
2224     if (exp > 0) // Normal.
2225       return Double.longBitsToDouble((bits & 0x800fffffffffffffL)
2226                                      | ((long) exp << 52));
2227     if (exp <= -54)
2228       return 0 * x; // Underflow.
2229     exp += 54; // Subnormal result.
2230     x = Double.longBitsToDouble((bits & 0x800fffffffffffffL)
2231                                 | ((long) exp << 52));
2232     return x * (1 / TWO_54);
2233   }
2234
2235   /**
2236    * Helper trig function; computes sin in range [-pi/4, pi/4].
2237    *
2238    * @param x angle within about pi/4
2239    * @param y tail of x, created by remPiOver2
2240    * @return sin(x+y)
2241    */
2242   private static double sin(double x, double y)
2243   {
2244     if (Configuration.DEBUG && abs(x + y) > 0.7854)
2245       throw new InternalError("Assertion failure");
2246     if (abs(x) < 1 / TWO_27)
2247       return x;  // If |x| ~< 2**-27, already know answer.
2248
2249     double z = x * x;
2250     double v = z * x;
2251     double r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
2252     if (y == 0)
2253       return x + v * (S1 + z * r);
2254     return x - ((z * (0.5 * y - v * r) - y) - v * S1);
2255   }
2256
2257   /**
2258    * Helper trig function; computes cos in range [-pi/4, pi/4].
2259    *
2260    * @param x angle within about pi/4
2261    * @param y tail of x, created by remPiOver2
2262    * @return cos(x+y)
2263    */
2264   private static double cos(double x, double y)
2265   {
2266     if (Configuration.DEBUG && abs(x + y) > 0.7854)
2267       throw new InternalError("Assertion failure");
2268     x = abs(x);
2269     if (x < 1 / TWO_27)
2270       return 1;  // If |x| ~< 2**-27, already know answer.
2271
2272     double z = x * x;
2273     double r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
2274
2275     if (x < 0.3)
2276       return 1 - (0.5 * z - (z * r - x * y));
2277
2278     double qx = (x > 0.78125) ? 0.28125 : (x * 0.25);
2279     return 1 - qx - ((0.5 * z - qx) - (z * r - x * y));
2280   }
2281
2282   /**
2283    * Helper trig function; computes tan in range [-pi/4, pi/4].
2284    *
2285    * @param x angle within about pi/4
2286    * @param y tail of x, created by remPiOver2
2287    * @param invert true iff -1/tan should be returned instead
2288    * @return tan(x+y)
2289    */
2290   private static double tan(double x, double y, boolean invert)
2291   {
2292     // PI/2 is irrational, so no double is a perfect multiple of it.
2293     if (Configuration.DEBUG && (abs(x + y) > 0.7854 || (x == 0 && invert)))
2294       throw new InternalError("Assertion failure");
2295     boolean negative = x < 0;
2296     if (negative)
2297       {
2298         x = -x;
2299         y = -y;
2300       }
2301     if (x < 1 / TWO_28) // If |x| ~< 2**-28, already know answer.
2302       return (negative ? -1 : 1) * (invert ? -1 / x : x);
2303
2304     double z;
2305     double w;
2306     boolean large = x >= 0.6744;
2307     if (large)
2308       {
2309         z = PI / 4 - x;
2310         w = PI_L / 4 - y;
2311         x = z + w;
2312         y = 0;
2313       }
2314     z = x * x;
2315     w = z * z;
2316     // Break x**5*(T1+x**2*T2+...) into
2317     //   x**5(T1+x**4*T3+...+x**20*T11)
2318     // + x**5(x**2*(T2+x**4*T4+...+x**22*T12)).
2319     double r = T1 + w * (T3 + w * (T5 + w * (T7 + w * (T9 + w * T11))));
2320     double v = z * (T2 + w * (T4 + w * (T6 + w * (T8 + w * (T10 + w * T12)))));
2321     double s = z * x;
2322     r = y + z * (s * (r + v) + y);
2323     r += T0 * s;
2324     w = x + r;
2325     if (large)
2326       {
2327         v = invert ? -1 : 1;
2328         return (negative ? -1 : 1) * (v - 2 * (x - (w * w / (w + v) - r)));
2329       }
2330     if (! invert)
2331       return w;
2332
2333     // Compute -1.0/(x+r) accurately.
2334     z = (float) w;
2335     v = r - (z - x);
2336     double a = -1 / w;
2337     double t = (float) a;
2338     return t + a * (1 + t * z + t * v);
2339   }
2340
2341   /**
2342    * <p>
2343    * Returns the sign of the argument as follows:
2344    * </p>
2345    * <ul>
2346    * <li>If <code>a</code> is greater than zero, the result is 1.0.</li>
2347    * <li>If <code>a</code> is less than zero, the result is -1.0.</li>
2348    * <li>If <code>a</code> is <code>NaN</code>, the result is <code>NaN</code>.
2349    * <li>If <code>a</code> is positive or negative zero, the result is the
2350    * same.</li>
2351    * </ul>
2352    *
2353    * @param a the numeric argument.
2354    * @return the sign of the argument.
2355    * @since 1.5.
2356    */
2357   public static double signum(double a)
2358   {
2359     // There's no difference.
2360     return Math.signum(a);
2361   }
2362
2363   /**
2364    * <p>
2365    * Returns the sign of the argument as follows:
2366    * </p>
2367    * <ul>
2368    * <li>If <code>a</code> is greater than zero, the result is 1.0f.</li>
2369    * <li>If <code>a</code> is less than zero, the result is -1.0f.</li>
2370    * <li>If <code>a</code> is <code>NaN</code>, the result is <code>NaN</code>.
2371    * <li>If <code>a</code> is positive or negative zero, the result is the
2372    * same.</li>
2373    * </ul>
2374    *
2375    * @param a the numeric argument.
2376    * @return the sign of the argument.
2377    * @since 1.5.
2378    */
2379   public static float signum(float a)
2380   {
2381     // There's no difference.
2382     return Math.signum(a);
2383   }
2384
2385   /**
2386    * Return the ulp for the given double argument.  The ulp is the
2387    * difference between the argument and the next larger double.  Note
2388    * that the sign of the double argument is ignored, that is,
2389    * ulp(x) == ulp(-x).  If the argument is a NaN, then NaN is returned.
2390    * If the argument is an infinity, then +Inf is returned.  If the
2391    * argument is zero (either positive or negative), then
2392    * {@link Double#MIN_VALUE} is returned.
2393    * @param d the double whose ulp should be returned
2394    * @return the difference between the argument and the next larger double
2395    * @since 1.5
2396    */
2397   public static double ulp(double d)
2398   {
2399     // There's no difference.
2400     return Math.ulp(d);
2401   }
2402
2403   /**
2404    * Return the ulp for the given float argument.  The ulp is the
2405    * difference between the argument and the next larger float.  Note
2406    * that the sign of the float argument is ignored, that is,
2407    * ulp(x) == ulp(-x).  If the argument is a NaN, then NaN is returned.
2408    * If the argument is an infinity, then +Inf is returned.  If the
2409    * argument is zero (either positive or negative), then
2410    * {@link Float#MIN_VALUE} is returned.
2411    * @param f the float whose ulp should be returned
2412    * @return the difference between the argument and the next larger float
2413    * @since 1.5
2414    */
2415   public static float ulp(float f)
2416   {
2417     // There's no difference.
2418     return Math.ulp(f);
2419   }
2420 }