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 >= <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 <= <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 }