001    /* java.lang.StrictMath -- common mathematical functions, strict Java
002       Copyright (C) 1998, 2001, 2002, 2003, 2006 Free Software Foundation, Inc.
003    
004    This file is part of GNU Classpath.
005    
006    GNU Classpath is free software; you can redistribute it and/or modify
007    it under the terms of the GNU General Public License as published by
008    the Free Software Foundation; either version 2, or (at your option)
009    any later version.
010    
011    GNU Classpath is distributed in the hope that it will be useful, but
012    WITHOUT ANY WARRANTY; without even the implied warranty of
013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
014    General Public License for more details.
015    
016    You should have received a copy of the GNU General Public License
017    along with GNU Classpath; see the file COPYING.  If not, write to the
018    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
019    02110-1301 USA.
020    
021    Linking this library statically or dynamically with other modules is
022    making a combined work based on this library.  Thus, the terms and
023    conditions of the GNU General Public License cover the whole
024    combination.
025    
026    As a special exception, the copyright holders of this library give you
027    permission to link this library with independent modules to produce an
028    executable, regardless of the license terms of these independent
029    modules, and to copy and distribute the resulting executable under
030    terms of your choice, provided that you also meet, for each linked
031    independent module, the terms and conditions of the license of that
032    module.  An independent module is a module which is not derived from
033    or based on this library.  If you modify this library, you may extend
034    this exception to your version of the library, but you are not
035    obligated to do so.  If you do not wish to do so, delete this
036    exception 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    
052    package java.lang;
053    
054    import gnu.classpath.Configuration;
055    
056    import 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     */
076    public 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    }