001/* java.lang.StrictMath -- common mathematical functions, strict Java
002   Copyright (C) 1998, 2001, 2002, 2003, 2006 Free Software Foundation, Inc.
003
004This file is part of GNU Classpath.
005
006GNU Classpath is free software; you can redistribute it and/or modify
007it under the terms of the GNU General Public License as published by
008the Free Software Foundation; either version 2, or (at your option)
009any later version.
010
011GNU Classpath is distributed in the hope that it will be useful, but
012WITHOUT ANY WARRANTY; without even the implied warranty of
013MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
014General Public License for more details.
015
016You should have received a copy of the GNU General Public License
017along with GNU Classpath; see the file COPYING.  If not, write to the
018Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
01902110-1301 USA.
020
021Linking this library statically or dynamically with other modules is
022making a combined work based on this library.  Thus, the terms and
023conditions of the GNU General Public License cover the whole
024combination.
025
026As a special exception, the copyright holders of this library give you
027permission to link this library with independent modules to produce an
028executable, regardless of the license terms of these independent
029modules, and to copy and distribute the resulting executable under
030terms of your choice, provided that you also meet, for each linked
031independent module, the terms and conditions of the license of that
032module.  An independent module is a module which is not derived from
033or based on this library.  If you modify this library, you may extend
034this exception to your version of the library, but you are not
035obligated to do so.  If you do not wish to do so, delete this
036exception statement from your version. */
037
038/*
039 * Some of the algorithms in this class are in the public domain, as part
040 * of fdlibm (freely-distributable math library), available at
041 * http://www.netlib.org/fdlibm/, and carry the following copyright:
042 * ====================================================
043 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
044 *
045 * Developed at SunSoft, a Sun Microsystems, Inc. business.
046 * Permission to use, copy, modify, and distribute this
047 * software is freely granted, provided that this notice
048 * is preserved.
049 * ====================================================
050 */
051
052package java.lang;
053
054import gnu.classpath.Configuration;
055
056import java.util.Random;
057
058/**
059 * Helper class containing useful mathematical functions and constants.
060 * This class mirrors {@link Math}, but is 100% portable, because it uses
061 * no native methods whatsoever.  Also, these algorithms are all accurate
062 * to less than 1 ulp, and execute in <code>strictfp</code> mode, while
063 * Math is allowed to vary in its results for some functions. Unfortunately,
064 * this usually means StrictMath has less efficiency and speed, as Math can
065 * use native methods.
066 *
067 * <p>The source of the various algorithms used is the fdlibm library, at:<br>
068 * <a href="http://www.netlib.org/fdlibm/">http://www.netlib.org/fdlibm/</a>
069 *
070 * Note that angles are specified in radians.  Conversion functions are
071 * provided for your convenience.
072 *
073 * @author Eric Blake (ebb9@email.byu.edu)
074 * @since 1.3
075 */
076public final strictfp class StrictMath
077{
078  /**
079   * StrictMath is non-instantiable.
080   */
081  private StrictMath()
082  {
083  }
084
085  /**
086   * A random number generator, initialized on first use.
087   *
088   * @see #random()
089   */
090  private static Random rand;
091
092  /**
093   * The most accurate approximation to the mathematical constant <em>e</em>:
094   * <code>2.718281828459045</code>. Used in natural log and exp.
095   *
096   * @see #log(double)
097   * @see #exp(double)
098   */
099  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}