1 /* java.lang.StrictMath -- common mathematical functions, strict Java
2 Copyright (C) 1998, 2001, 2002 Free Software Foundation, Inc.
4 This file is part of GNU Classpath.
6 GNU Classpath is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2, or (at your option)
11 GNU Classpath is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with GNU Classpath; see the file COPYING. If not, write to the
18 Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
21 Linking this library statically or dynamically with other modules is
22 making a combined work based on this library. Thus, the terms and
23 conditions of the GNU General Public License cover the whole
26 As a special exception, the copyright holders of this library give you
27 permission to link this library with independent modules to produce an
28 executable, regardless of the license terms of these independent
29 modules, and to copy and distribute the resulting executable under
30 terms of your choice, provided that you also meet, for each linked
31 independent module, the terms and conditions of the license of that
32 module. An independent module is a module which is not derived from
33 or based on this library. If you modify this library, you may extend
34 this exception to your version of the library, but you are not
35 obligated to do so. If you do not wish to do so, delete this
36 exception statement from your version. */
39 * Some of the algorithms in this class are in the public domain, as part
40 * of fdlibm (freely-distributable math library), available at
41 * http://www.netlib.org/fdlibm/, and carry the following copyright:
42 * ====================================================
43 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
45 * Developed at SunSoft, a Sun Microsystems, Inc. business.
46 * Permission to use, copy, modify, and distribute this
47 * software is freely granted, provided that this notice
49 * ====================================================
54 import java
.util
.Random
;
55 import gnu
.classpath
.Configuration
;
58 * Helper class containing useful mathematical functions and constants.
59 * This class mirrors {@link Math}, but is 100% portable, because it uses
60 * no native methods whatsoever. Also, these algorithms are all accurate
61 * to less than 1 ulp, and execute in <code>strictfp</code> mode, while
62 * Math is allowed to vary in its results for some functions. Unfortunately,
63 * this usually means StrictMath has less efficiency and speed, as Math can
66 * <p>The source of the various algorithms used is the fdlibm library, at:<br>
67 * <a href="http://www.netlib.org/fdlibm/">http://www.netlib.org/fdlibm/</a>
69 * Note that angles are specified in radians. Conversion functions are
70 * provided for your convenience.
72 * @author Eric Blake <ebb9@email.byu.edu>
75 public final strictfp class StrictMath
78 * StrictMath is non-instantiable.
85 * A random number generator, initialized on first use.
89 private static Random rand
;
92 * The most accurate approximation to the mathematical constant <em>e</em>:
93 * <code>2.718281828459045</code>. Used in natural log and exp.
98 public static final double E
99 = 2.718281828459045; // Long bits 0x4005bf0z8b145769L.
102 * The most accurate approximation to the mathematical constant <em>pi</em>:
103 * <code>3.141592653589793</code>. This is the ratio of a circle's diameter
104 * to its circumference.
106 public static final double PI
107 = 3.141592653589793; // Long bits 0x400921fb54442d18L.
110 * Take the absolute value of the argument. (Absolute value means make
113 * <p>Note that the the largest negative value (Integer.MIN_VALUE) cannot
114 * be made positive. In this case, because of the rules of negation in
115 * a computer, MIN_VALUE is what will be returned.
116 * This is a <em>negative</em> value. You have been warned.
118 * @param i the number to take the absolute value of
119 * @return the absolute value
120 * @see Integer#MIN_VALUE
122 public static int abs(int i
)
124 return (i
< 0) ?
-i
: i
;
128 * Take the absolute value of the argument. (Absolute value means make
131 * <p>Note that the the largest negative value (Long.MIN_VALUE) cannot
132 * be made positive. In this case, because of the rules of negation in
133 * a computer, MIN_VALUE is what will be returned.
134 * This is a <em>negative</em> value. You have been warned.
136 * @param l the number to take the absolute value of
137 * @return the absolute value
138 * @see Long#MIN_VALUE
140 public static long abs(long l
)
142 return (l
< 0) ?
-l
: l
;
146 * Take the absolute value of the argument. (Absolute value means make
149 * @param f the number to take the absolute value of
150 * @return the absolute value
152 public static float abs(float f
)
154 return (f
<= 0) ?
0 - f
: f
;
158 * Take the absolute value of the argument. (Absolute value means make
161 * @param d the number to take the absolute value of
162 * @return the absolute value
164 public static double abs(double d
)
166 return (d
<= 0) ?
0 - d
: d
;
170 * Return whichever argument is smaller.
172 * @param a the first number
173 * @param b a second number
174 * @return the smaller of the two numbers
176 public static int min(int a
, int b
)
178 return (a
< b
) ? a
: b
;
182 * Return whichever argument is smaller.
184 * @param a the first number
185 * @param b a second number
186 * @return the smaller of the two numbers
188 public static long min(long a
, long b
)
190 return (a
< b
) ? a
: b
;
194 * Return whichever argument is smaller. If either argument is NaN, the
195 * result is NaN, and when comparing 0 and -0, -0 is always smaller.
197 * @param a the first number
198 * @param b a second number
199 * @return the smaller of the two numbers
201 public static float min(float a
, float b
)
203 // this check for NaN, from JLS 15.21.1, saves a method call
206 // no need to check if b is NaN; < will work correctly
207 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special
208 if (a
== 0 && b
== 0)
210 return (a
< b
) ? a
: b
;
214 * Return whichever argument is smaller. If either argument is NaN, the
215 * result is NaN, and when comparing 0 and -0, -0 is always smaller.
217 * @param a the first number
218 * @param b a second number
219 * @return the smaller of the two numbers
221 public static double min(double a
, double b
)
223 // this check for NaN, from JLS 15.21.1, saves a method call
226 // no need to check if b is NaN; < will work correctly
227 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special
228 if (a
== 0 && b
== 0)
230 return (a
< b
) ? a
: b
;
234 * Return whichever argument is larger.
236 * @param a the first number
237 * @param b a second number
238 * @return the larger of the two numbers
240 public static int max(int a
, int b
)
242 return (a
> b
) ? a
: b
;
246 * Return whichever argument is larger.
248 * @param a the first number
249 * @param b a second number
250 * @return the larger of the two numbers
252 public static long max(long a
, long b
)
254 return (a
> b
) ? a
: b
;
258 * Return whichever argument is larger. If either argument is NaN, the
259 * result is NaN, and when comparing 0 and -0, 0 is always larger.
261 * @param a the first number
262 * @param b a second number
263 * @return the larger of the two numbers
265 public static float max(float a
, float b
)
267 // this check for NaN, from JLS 15.21.1, saves a method call
270 // no need to check if b is NaN; > will work correctly
271 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special
272 if (a
== 0 && b
== 0)
274 return (a
> b
) ? a
: b
;
278 * Return whichever argument is larger. If either argument is NaN, the
279 * result is NaN, and when comparing 0 and -0, 0 is always larger.
281 * @param a the first number
282 * @param b a second number
283 * @return the larger of the two numbers
285 public static double max(double a
, double b
)
287 // this check for NaN, from JLS 15.21.1, saves a method call
290 // no need to check if b is NaN; > will work correctly
291 // recall that -0.0 == 0.0, but [+-]0.0 - [+-]0.0 behaves special
292 if (a
== 0 && b
== 0)
294 return (a
> b
) ? a
: b
;
298 * The trigonometric function <em>sin</em>. The sine of NaN or infinity is
299 * NaN, and the sine of 0 retains its sign.
301 * @param a the angle (in radians)
304 public static double sin(double a
)
306 if (a
== Double
.NEGATIVE_INFINITY
|| ! (a
< Double
.POSITIVE_INFINITY
))
309 if (abs(a
) <= PI
/ 4)
312 // Argument reduction needed.
313 double[] y
= new double[2];
314 int n
= remPiOver2(a
, y
);
318 return sin(y
[0], y
[1]);
320 return cos(y
[0], y
[1]);
322 return -sin(y
[0], y
[1]);
324 return -cos(y
[0], y
[1]);
329 * The trigonometric function <em>cos</em>. The cosine of NaN or infinity is
332 * @param a the angle (in radians).
335 public static double cos(double a
)
337 if (a
== Double
.NEGATIVE_INFINITY
|| ! (a
< Double
.POSITIVE_INFINITY
))
340 if (abs(a
) <= PI
/ 4)
343 // Argument reduction needed.
344 double[] y
= new double[2];
345 int n
= remPiOver2(a
, y
);
349 return cos(y
[0], y
[1]);
351 return -sin(y
[0], y
[1]);
353 return -cos(y
[0], y
[1]);
355 return sin(y
[0], y
[1]);
360 * The trigonometric function <em>tan</em>. The tangent of NaN or infinity
361 * is NaN, and the tangent of 0 retains its sign.
363 * @param a the angle (in radians)
366 public static double tan(double a
)
368 if (a
== Double
.NEGATIVE_INFINITY
|| ! (a
< Double
.POSITIVE_INFINITY
))
371 if (abs(a
) <= PI
/ 4)
372 return tan(a
, 0, false);
374 // Argument reduction needed.
375 double[] y
= new double[2];
376 int n
= remPiOver2(a
, y
);
377 return tan(y
[0], y
[1], (n
& 1) == 1);
381 * The trigonometric function <em>arcsin</em>. The range of angles returned
382 * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN or
383 * its absolute value is beyond 1, the result is NaN; and the arcsine of
384 * 0 retains its sign.
386 * @param x the sin to turn back into an angle
389 public static double asin(double x
)
391 boolean negative
= x
< 0;
397 return negative ?
-PI
/ 2 : PI
/ 2;
401 return negative ?
-x
: x
;
403 double p
= t
* (PS0
+ t
* (PS1
+ t
* (PS2
+ t
* (PS3
+ t
404 * (PS4
+ t
* PS5
)))));
405 double q
= 1 + t
* (QS1
+ t
* (QS2
+ t
* (QS3
+ t
* QS4
)));
406 return negative ?
-x
- x
* (p
/ q
) : x
+ x
* (p
/ q
);
408 double w
= 1 - x
; // 1>|x|>=0.5.
410 double p
= t
* (PS0
+ t
* (PS1
+ t
* (PS2
+ t
* (PS3
+ t
411 * (PS4
+ t
* PS5
)))));
412 double q
= 1 + t
* (QS1
+ t
* (QS2
+ t
* (QS3
+ t
* QS4
)));
417 t
= PI
/ 2 - (2 * (s
+ s
* w
) - PI_L
/ 2);
422 double c
= (t
- w
* w
) / (s
+ w
);
423 p
= 2 * s
* (p
/ q
) - (PI_L
/ 2 - 2 * c
);
425 t
= PI
/ 4 - (p
- q
);
427 return negative ?
-t
: t
;
431 * The trigonometric function <em>arccos</em>. The range of angles returned
432 * is 0 to pi radians (0 to 180 degrees). If the argument is NaN or
433 * its absolute value is beyond 1, the result is NaN.
435 * @param x the cos to turn back into an angle
438 public static double acos(double x
)
440 boolean negative
= x
< 0;
446 return negative ? PI
: 0;
452 double p
= z
* (PS0
+ z
* (PS1
+ z
* (PS2
+ z
* (PS3
+ z
453 * (PS4
+ z
* PS5
)))));
454 double q
= 1 + z
* (QS1
+ z
* (QS2
+ z
* (QS3
+ z
* QS4
)));
455 double r
= x
- (PI_L
/ 2 - x
* (p
/ q
));
456 return negative ? PI
/ 2 + r
: PI
/ 2 - r
;
458 if (negative
) // x<=-0.5.
460 double z
= (1 + x
) * 0.5;
461 double p
= z
* (PS0
+ z
* (PS1
+ z
* (PS2
+ z
* (PS3
+ z
462 * (PS4
+ z
* PS5
)))));
463 double q
= 1 + z
* (QS1
+ z
* (QS2
+ z
* (QS3
+ z
* QS4
)));
465 double w
= p
/ q
* s
- PI_L
/ 2;
466 return PI
- 2 * (s
+ w
);
468 double z
= (1 - x
) * 0.5; // x>0.5.
470 double df
= (float) s
;
471 double c
= (z
- df
* df
) / (s
+ df
);
472 double p
= z
* (PS0
+ z
* (PS1
+ z
* (PS2
+ z
* (PS3
+ z
473 * (PS4
+ z
* PS5
)))));
474 double q
= 1 + z
* (QS1
+ z
* (QS2
+ z
* (QS3
+ z
* QS4
)));
475 double w
= p
/ q
* s
+ c
;
480 * The trigonometric function <em>arcsin</em>. The range of angles returned
481 * is -pi/2 to pi/2 radians (-90 to 90 degrees). If the argument is NaN, the
482 * result is NaN; and the arctangent of 0 retains its sign.
484 * @param x the tan to turn back into an angle
486 * @see #atan2(double, double)
488 public static double atan(double x
)
492 boolean negative
= x
< 0;
496 return negative ?
-PI
/ 2 : PI
/ 2;
497 if (! (x
>= 0.4375)) // |x|<7/16, or NaN.
499 if (! (x
>= 1 / TWO_29
)) // Small, or NaN.
500 return negative ?
-x
: x
;
505 if (x
< 0.6875) // 7/16<=|x|<11/16.
507 x
= (2 * x
- 1) / (2 + x
);
511 else // 11/16<=|x|<19/16.
513 x
= (x
- 1) / (x
+ 1);
518 else if (x
< 2.4375) // 19/16<=|x|<39/16.
520 x
= (x
- 1.5) / (1 + 1.5 * x
);
524 else // 39/16<=|x|<2**66.
531 // Break sum from i=0 to 10 ATi*z**(i+1) into odd and even poly.
534 double s1
= z
* (AT0
+ w
* (AT2
+ w
* (AT4
+ w
* (AT6
+ w
535 * (AT8
+ w
* AT10
)))));
536 double s2
= w
* (AT1
+ w
* (AT3
+ w
* (AT5
+ w
* (AT7
+ w
* AT9
))));
538 return negative ? x
* (s1
+ s2
) - x
: x
- x
* (s1
+ s2
);
539 z
= hi
- ((x
* (s1
+ s2
) - lo
) - x
);
540 return negative ?
-z
: z
;
544 * A special version of the trigonometric function <em>arctan</em>, for
545 * converting rectangular coordinates <em>(x, y)</em> to polar
546 * <em>(r, theta)</em>. This computes the arctangent of x/y in the range
547 * of -pi to pi radians (-180 to 180 degrees). Special cases:<ul>
548 * <li>If either argument is NaN, the result is NaN.</li>
549 * <li>If the first argument is positive zero and the second argument is
550 * positive, or the first argument is positive and finite and the second
551 * argument is positive infinity, then the result is positive zero.</li>
552 * <li>If the first argument is negative zero and the second argument is
553 * positive, or the first argument is negative and finite and the second
554 * argument is positive infinity, then the result is negative zero.</li>
555 * <li>If the first argument is positive zero and the second argument is
556 * negative, or the first argument is positive and finite and the second
557 * argument is negative infinity, then the result is the double value
558 * closest to pi.</li>
559 * <li>If the first argument is negative zero and the second argument is
560 * negative, or the first argument is negative and finite and the second
561 * argument is negative infinity, then the result is the double value
562 * closest to -pi.</li>
563 * <li>If the first argument is positive and the second argument is
564 * positive zero or negative zero, or the first argument is positive
565 * infinity and the second argument is finite, then the result is the
566 * double value closest to pi/2.</li>
567 * <li>If the first argument is negative and the second argument is
568 * positive zero or negative zero, or the first argument is negative
569 * infinity and the second argument is finite, then the result is the
570 * double value closest to -pi/2.</li>
571 * <li>If both arguments are positive infinity, then the result is the
572 * double value closest to pi/4.</li>
573 * <li>If the first argument is positive infinity and the second argument
574 * is negative infinity, then the result is the double value closest to
576 * <li>If the first argument is negative infinity and the second argument
577 * is positive infinity, then the result is the double value closest to
579 * <li>If both arguments are negative infinity, then the result is the
580 * double value closest to -3*pi/4.</li>
582 * </ul><p>This returns theta, the angle of the point. To get r, albeit
583 * slightly inaccurately, use sqrt(x*x+y*y).
585 * @param y the y position
586 * @param x the x position
587 * @return <em>theta</em> in the conversion of (x, y) to (r, theta)
590 public static double atan2(double y
, double x
)
592 if (x
!= x
|| y
!= y
)
596 if (x
== Double
.POSITIVE_INFINITY
)
598 if (y
== Double
.POSITIVE_INFINITY
)
600 if (y
== Double
.NEGATIVE_INFINITY
)
604 if (x
== Double
.NEGATIVE_INFINITY
)
606 if (y
== Double
.POSITIVE_INFINITY
)
608 if (y
== Double
.NEGATIVE_INFINITY
)
610 return (1 / (0 * y
) == Double
.POSITIVE_INFINITY
) ? PI
: -PI
;
614 if (1 / (0 * x
) == Double
.POSITIVE_INFINITY
)
616 return (1 / y
== Double
.POSITIVE_INFINITY
) ? PI
: -PI
;
618 if (y
== Double
.POSITIVE_INFINITY
|| y
== Double
.NEGATIVE_INFINITY
620 return y
< 0 ?
-PI
/ 2 : PI
/ 2;
622 double z
= abs(y
/ x
); // Safe to do y/x.
624 z
= PI
/ 2 + 0.5 * PI_L
;
625 else if (x
< 0 && z
< 1 / TWO_60
)
630 return y
> 0 ? z
: -z
;
631 return y
> 0 ? PI
- (z
- PI_L
) : z
- PI_L
- PI
;
635 * Take <em>e</em><sup>a</sup>. The opposite of <code>log()</code>. If the
636 * argument is NaN, the result is NaN; if the argument is positive infinity,
637 * the result is positive infinity; and if the argument is negative
638 * infinity, the result is positive zero.
640 * @param x the number to raise to the power
641 * @return the number raised to the power of <em>e</em>
643 * @see #pow(double, double)
645 public static double exp(double x
)
650 return Double
.POSITIVE_INFINITY
;
654 // Argument reduction.
669 k
= (int) (INV_LN2
* t
+ 0.5);
681 else if (t
< 1 / TWO_28
)
686 // Now x is in primary range.
688 double c
= x
- t
* (P1
+ t
* (P2
+ t
* (P3
+ t
* (P4
+ t
* P5
))));
690 return 1 - (x
* c
/ (c
- 2) - x
);
691 double y
= 1 - (lo
- x
* c
/ (2 - c
) - hi
);
696 * Take ln(a) (the natural log). The opposite of <code>exp()</code>. If the
697 * argument is NaN or negative, the result is NaN; if the argument is
698 * positive infinity, the result is positive infinity; and if the argument
699 * is either zero, the result is negative infinity.
701 * <p>Note that the way to get log<sub>b</sub>(a) is to do this:
702 * <code>ln(a) / ln(b)</code>.
704 * @param x the number to take the natural log of
705 * @return the natural log of <code>a</code>
708 public static double log(double x
)
711 return Double
.NEGATIVE_INFINITY
;
714 if (! (x
< Double
.POSITIVE_INFINITY
))
718 long bits
= Double
.doubleToLongBits(x
);
719 int exp
= (int) (bits
>> 52);
720 if (exp
== 0) // Subnormal x.
723 bits
= Double
.doubleToLongBits(x
);
724 exp
= (int) (bits
>> 52) - 54;
726 exp
-= 1023; // Unbias exponent.
727 bits
= (bits
& 0x000fffffffffffffL
) | 0x3ff0000000000000L
;
728 x
= Double
.longBitsToDouble(bits
);
735 if (abs(x
) < 1 / TWO_20
)
738 return exp
* LN2_H
+ exp
* LN2_L
;
739 double r
= x
* x
* (0.5 - 1 / 3.0 * x
);
742 return exp
* LN2_H
- ((r
- exp
* LN2_L
) - x
);
744 double s
= x
/ (2 + x
);
747 double t1
= w
* (LG2
+ w
* (LG4
+ w
* LG6
));
748 double t2
= z
* (LG1
+ w
* (LG3
+ w
* (LG5
+ w
* LG7
)));
750 if (bits
>= 0x3ff6174a00000000L
&& bits
< 0x3ff6b85200000000L
)
752 double h
= 0.5 * x
* x
; // Need more accuracy for x near sqrt(2).
754 return x
- (h
- s
* (h
+ r
));
755 return exp
* LN2_H
- ((h
- (s
* (h
+ r
) + exp
* LN2_L
)) - x
);
758 return x
- s
* (x
- r
);
759 return exp
* LN2_H
- ((s
* (x
- r
) - exp
* LN2_L
) - x
);
763 * Take a square root. If the argument is NaN or negative, the result is
764 * NaN; if the argument is positive infinity, the result is positive
765 * infinity; and if the result is either zero, the result is the same.
767 * <p>For other roots, use pow(x, 1/rootNumber).
769 * @param x the numeric argument
770 * @return the square root of the argument
771 * @see #pow(double, double)
773 public static double sqrt(double x
)
777 if (x
== 0 || ! (x
< Double
.POSITIVE_INFINITY
))
781 long bits
= Double
.doubleToLongBits(x
);
782 int exp
= (int) (bits
>> 52);
783 if (exp
== 0) // Subnormal x.
786 bits
= Double
.doubleToLongBits(x
);
787 exp
= (int) (bits
>> 52) - 54;
789 exp
-= 1023; // Unbias exponent.
790 bits
= (bits
& 0x000fffffffffffffL
) | 0x0010000000000000L
;
791 if ((exp
& 1) == 1) // Odd exp, double x to make it even.
795 // Generate sqrt(x) bit by bit.
799 long r
= 0x0020000000000000L
; // Move r right to left.
813 // Use floating add to round correctly.
816 return Double
.longBitsToDouble((q
>> 1) + ((exp
+ 1022L) << 52));
820 * Raise a number to a power. Special cases:<ul>
821 * <li>If the second argument is positive or negative zero, then the result
823 * <li>If the second argument is 1.0, then the result is the same as the
824 * first argument.</li>
825 * <li>If the second argument is NaN, then the result is NaN.</li>
826 * <li>If the first argument is NaN and the second argument is nonzero,
827 * then the result is NaN.</li>
828 * <li>If the absolute value of the first argument is greater than 1 and
829 * the second argument is positive infinity, or the absolute value of the
830 * first argument is less than 1 and the second argument is negative
831 * infinity, then the result is positive infinity.</li>
832 * <li>If the absolute value of the first argument is greater than 1 and
833 * the second argument is negative infinity, or the absolute value of the
834 * first argument is less than 1 and the second argument is positive
835 * infinity, then the result is positive zero.</li>
836 * <li>If the absolute value of the first argument equals 1 and the second
837 * argument is infinite, then the result is NaN.</li>
838 * <li>If the first argument is positive zero and the second argument is
839 * greater than zero, or the first argument is positive infinity and the
840 * second argument is less than zero, then the result is positive zero.</li>
841 * <li>If the first argument is positive zero and the second argument is
842 * less than zero, or the first argument is positive infinity and the
843 * second argument is greater than zero, then the result is positive
845 * <li>If the first argument is negative zero and the second argument is
846 * greater than zero but not a finite odd integer, or the first argument is
847 * negative infinity and the second argument is less than zero but not a
848 * finite odd integer, then the result is positive zero.</li>
849 * <li>If the first argument is negative zero and the second argument is a
850 * positive finite odd integer, or the first argument is negative infinity
851 * and the second argument is a negative finite odd integer, then the result
852 * is negative zero.</li>
853 * <li>If the first argument is negative zero and the second argument is
854 * less than zero but not a finite odd integer, or the first argument is
855 * negative infinity and the second argument is greater than zero but not a
856 * finite odd integer, then the result is positive infinity.</li>
857 * <li>If the first argument is negative zero and the second argument is a
858 * negative finite odd integer, or the first argument is negative infinity
859 * and the second argument is a positive finite odd integer, then the result
860 * is negative infinity.</li>
861 * <li>If the first argument is less than zero and the second argument is a
862 * finite even integer, then the result is equal to the result of raising
863 * the absolute value of the first argument to the power of the second
865 * <li>If the first argument is less than zero and the second argument is a
866 * finite odd integer, then the result is equal to the negative of the
867 * result of raising the absolute value of the first argument to the power
868 * of the second argument.</li>
869 * <li>If the first argument is finite and less than zero and the second
870 * argument is finite and not an integer, then the result is NaN.</li>
871 * <li>If both arguments are integers, then the result is exactly equal to
872 * the mathematical result of raising the first argument to the power of
873 * the second argument if that result can in fact be represented exactly as
874 * a double value.</li>
876 * </ul><p>(In the foregoing descriptions, a floating-point value is
877 * considered to be an integer if and only if it is a fixed point of the
878 * method {@link #ceil(double)} or, equivalently, a fixed point of the
879 * method {@link #floor(double)}. A value is a fixed point of a one-argument
880 * method if and only if the result of applying the method to the value is
881 * equal to the value.)
883 * @param x the number to raise
884 * @param y the power to raise it to
885 * @return x<sup>y</sup>
887 public static double pow(double x
, double y
)
889 // Special cases first.
896 if (x
!= x
|| y
!= y
)
899 // When x < 0, yisint tells if y is not an integer (0), even(1),
902 if (x
< 0 && floor(y
) == y
)
903 yisint
= (y
% 2 == 0) ?
2 : 1;
907 // More special cases, of y.
908 if (ay
== Double
.POSITIVE_INFINITY
)
913 return y
> 0 ? y
: 0;
914 return y
< 0 ?
-y
: 0;
921 // More special cases, of x.
922 if (x
== 0 || ax
== Double
.POSITIVE_INFINITY
|| ax
== 1)
928 if (x
== -1 && yisint
== 0)
930 else if (yisint
== 1)
935 if (x
< 0 && yisint
== 0)
947 if (ay
> TWO_64
) // Automatic over/underflow.
948 return ((ax
< 1) ? y
< 0 : y
> 0) ? Double
.POSITIVE_INFINITY
: 0;
949 // Over/underflow if x is not close to one.
950 if (ax
< 0.9999995231628418)
951 return y
< 0 ? Double
.POSITIVE_INFINITY
: 0;
952 if (ax
>= 1.0000009536743164)
953 return y
> 0 ? Double
.POSITIVE_INFINITY
: 0;
954 // Now |1-x| is <= 2**-20, sufficient to compute
955 // log(x) by x-x^2/2+x^3/3-x^4/4.
957 w
= t
* t
* (0.5 - t
* (1 / 3.0 - t
* 0.25));
959 v
= t
* INV_LN2_L
- w
* INV_LN2
;
960 t1
= (float) (u
+ v
);
965 long bits
= Double
.doubleToLongBits(ax
);
966 int exp
= (int) (bits
>> 52);
967 if (exp
== 0) // Subnormal x.
970 bits
= Double
.doubleToLongBits(ax
);
971 exp
= (int) (bits
>> 52) - 54;
973 exp
-= 1023; // Unbias exponent.
974 ax
= Double
.longBitsToDouble((bits
& 0x000fffffffffffffL
)
975 | 0x3ff0000000000000L
);
977 if (ax
< SQRT_1_5
) // |x|<sqrt(3/2).
979 else if (ax
< SQRT_3
) // |x|<sqrt(3).
988 // Compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5).
989 u
= ax
- (k ?
1.5 : 1);
990 v
= 1 / (ax
+ (k ?
1.5 : 1));
992 double s_h
= (float) s
;
993 double t_h
= (float) (ax
+ (k ?
1.5 : 1));
994 double t_l
= ax
- (t_h
- (k ?
1.5 : 1));
995 double s_l
= v
* ((u
- s_h
* t_h
) - s_h
* t_l
);
998 double r
= s_l
* (s_h
+ s
) + s2
* s2
999 * (L1
+ s2
* (L2
+ s2
* (L3
+ s2
* (L4
+ s2
* (L5
+ s2
* L6
)))));
1001 t_h
= (float) (3.0 + s2
+ r
);
1002 t_l
= r
- (t_h
- 3.0 - s2
);
1005 v
= s_l
* t_h
+ t_l
* s
;
1006 // 2/(3log2)*(s+...).
1007 double p_h
= (float) (u
+ v
);
1008 double p_l
= v
- (p_h
- u
);
1009 double z_h
= CP_H
* p_h
;
1010 double z_l
= CP_L
* p_h
+ p_l
* CP
+ (k ? DP_L
: 0);
1011 // log2(ax) = (s+..)*2/(3*log2) = exp + dp_h + z_h + z_l.
1013 t1
= (float) (z_h
+ z_l
+ (k ? DP_H
: 0) + t
);
1014 t2
= z_l
- (t1
- t
- (k ? DP_H
: 0) - z_h
);
1017 // Split up y into y1+y2 and compute (y1+y2)*(t1+t2).
1018 boolean negative
= x
< 0 && yisint
== 1;
1019 double y1
= (float) y
;
1020 double p_l
= (y
- y1
) * t1
+ y
* t2
;
1021 double p_h
= y1
* t1
;
1022 double z
= p_l
+ p_h
;
1023 if (z
>= 1024) // Detect overflow.
1025 if (z
> 1024 || p_l
+ OVT
> z
- p_h
)
1026 return negative ? Double
.NEGATIVE_INFINITY
1027 : Double
.POSITIVE_INFINITY
;
1029 else if (z
<= -1075) // Detect underflow.
1031 if (z
< -1075 || p_l
<= z
- p_h
)
1032 return negative ?
-0.0 : 0;
1035 // Compute 2**(p_h+p_l).
1036 int n
= round((float) z
);
1038 t
= (float) (p_l
+ p_h
);
1040 v
= (p_l
- (t
- p_h
)) * LN2
+ t
* LN2_L
;
1044 t1
= z
- t
* (P1
+ t
* (P2
+ t
* (P3
+ t
* (P4
+ t
* P5
))));
1045 double r
= (z
* t1
) / (t1
- 2) - (w
+ z
* w
);
1046 z
= scale(1 - (r
- z
), n
);
1047 return negative ?
-z
: z
;
1051 * Get the IEEE 754 floating point remainder on two numbers. This is the
1052 * value of <code>x - y * <em>n</em></code>, where <em>n</em> is the closest
1053 * double to <code>x / y</code> (ties go to the even n); for a zero
1054 * remainder, the sign is that of <code>x</code>. If either argument is NaN,
1055 * the first argument is infinite, or the second argument is zero, the result
1056 * is NaN; if x is finite but y is infinte, the result is x.
1058 * @param x the dividend (the top half)
1059 * @param y the divisor (the bottom half)
1060 * @return the IEEE 754-defined floating point remainder of x/y
1061 * @see #rint(double)
1063 public static double IEEEremainder(double x
, double y
)
1065 // Purge off exception values.
1066 if (x
== Double
.NEGATIVE_INFINITY
|| ! (x
< Double
.POSITIVE_INFINITY
)
1067 || y
== 0 || y
!= y
)
1070 boolean negative
= x
< 0;
1073 if (x
== y
|| x
== 0)
1074 return 0 * x
; // Get correct sign.
1076 // Achieve x < 2y, then take first shot at remainder.
1080 // Now adjust x to get correct precision.
1081 if (y
< 4 / TWO_1023
)
1100 return negative ?
-x
: x
;
1104 * Take the nearest integer that is that is greater than or equal to the
1105 * argument. If the argument is NaN, infinite, or zero, the result is the
1106 * same; if the argument is between -1 and 0, the result is negative zero.
1107 * Note that <code>Math.ceil(x) == -Math.floor(-x)</code>.
1109 * @param a the value to act upon
1110 * @return the nearest integer >= <code>a</code>
1112 public static double ceil(double a
)
1118 * Take the nearest integer that is that is less than or equal to the
1119 * argument. If the argument is NaN, infinite, or zero, the result is the
1120 * same. Note that <code>Math.ceil(x) == -Math.floor(-x)</code>.
1122 * @param a the value to act upon
1123 * @return the nearest integer <= <code>a</code>
1125 public static double floor(double a
)
1128 if (! (x
< TWO_52
) || (long) a
== a
)
1129 return a
; // No fraction bits; includes NaN and infinity.
1131 return a
>= 0 ?
0 * a
: -1; // Worry about signed zero.
1132 return a
< 0 ?
(long) a
- 1.0 : (long) a
; // Cast to long truncates.
1136 * Take the nearest integer to the argument. If it is exactly between
1137 * two integers, the even integer is taken. If the argument is NaN,
1138 * infinite, or zero, the result is the same.
1140 * @param a the value to act upon
1141 * @return the nearest integer to <code>a</code>
1143 public static double rint(double a
)
1147 return a
; // No fraction bits; includes NaN and infinity.
1149 return 0 * a
; // Worry about signed zero.
1151 return (long) a
; // Catch round down to even.
1152 return (long) (a
+ (a
< 0 ?
-0.5 : 0.5)); // Cast to long truncates.
1156 * Take the nearest integer to the argument. This is equivalent to
1157 * <code>(int) Math.floor(f + 0.5f)</code>. If the argument is NaN, the
1158 * result is 0; otherwise if the argument is outside the range of int, the
1159 * result will be Integer.MIN_VALUE or Integer.MAX_VALUE, as appropriate.
1161 * @param f the argument to round
1162 * @return the nearest integer to the argument
1163 * @see Integer#MIN_VALUE
1164 * @see Integer#MAX_VALUE
1166 public static int round(float f
)
1168 return (int) floor(f
+ 0.5f
);
1172 * Take the nearest long to the argument. This is equivalent to
1173 * <code>(long) Math.floor(d + 0.5)</code>. If the argument is NaN, the
1174 * result is 0; otherwise if the argument is outside the range of long, the
1175 * result will be Long.MIN_VALUE or Long.MAX_VALUE, as appropriate.
1177 * @param d the argument to round
1178 * @return the nearest long to the argument
1179 * @see Long#MIN_VALUE
1180 * @see Long#MAX_VALUE
1182 public static long round(double d
)
1184 return (long) floor(d
+ 0.5);
1188 * Get a random number. This behaves like Random.nextDouble(), seeded by
1189 * System.currentTimeMillis() when first called. In other words, the number
1190 * is from a pseudorandom sequence, and lies in the range [+0.0, 1.0).
1191 * This random sequence is only used by this method, and is threadsafe,
1192 * although you may want your own random number generator if it is shared
1195 * @return a random number
1196 * @see Random#nextDouble()
1197 * @see System#currentTimeMillis()
1199 public static synchronized double random()
1202 rand
= new Random();
1203 return rand
.nextDouble();
1207 * Convert from degrees to radians. The formula for this is
1208 * radians = degrees * (pi/180); however it is not always exact given the
1209 * limitations of floating point numbers.
1211 * @param degrees an angle in degrees
1212 * @return the angle in radians
1214 public static double toRadians(double degrees
)
1216 return degrees
* (PI
/ 180);
1220 * Convert from radians to degrees. The formula for this is
1221 * degrees = radians * (180/pi); however it is not always exact given the
1222 * limitations of floating point numbers.
1224 * @param rads an angle in radians
1225 * @return the angle in degrees
1227 public static double toDegrees(double rads
)
1229 return rads
* (180 / PI
);
1233 * Constants for scaling and comparing doubles by powers of 2. The compiler
1234 * must automatically inline constructs like (1/TWO_54), so we don't list
1235 * negative powers of two here.
1237 private static final double
1238 TWO_16
= 0x10000, // Long bits 0x40f0000000000000L.
1239 TWO_20
= 0x100000, // Long bits 0x4130000000000000L.
1240 TWO_24
= 0x1000000, // Long bits 0x4170000000000000L.
1241 TWO_27
= 0x8000000, // Long bits 0x41a0000000000000L.
1242 TWO_28
= 0x10000000, // Long bits 0x41b0000000000000L.
1243 TWO_29
= 0x20000000, // Long bits 0x41c0000000000000L.
1244 TWO_31
= 0x80000000L
, // Long bits 0x41e0000000000000L.
1245 TWO_49
= 0x2000000000000L
, // Long bits 0x4300000000000000L.
1246 TWO_52
= 0x10000000000000L
, // Long bits 0x4330000000000000L.
1247 TWO_54
= 0x40000000000000L
, // Long bits 0x4350000000000000L.
1248 TWO_57
= 0x200000000000000L
, // Long bits 0x4380000000000000L.
1249 TWO_60
= 0x1000000000000000L
, // Long bits 0x43b0000000000000L.
1250 TWO_64
= 1.8446744073709552e19
, // Long bits 0x43f0000000000000L.
1251 TWO_66
= 7.378697629483821e19
, // Long bits 0x4410000000000000L.
1252 TWO_1023
= 8.98846567431158e307
; // Long bits 0x7fe0000000000000L.
1255 * Super precision for 2/pi in 24-bit chunks, for use in
1256 * {@link #remPiOver2()}.
1258 private static final int TWO_OVER_PI
[] = {
1259 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
1260 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
1261 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
1262 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
1263 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
1264 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
1265 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
1266 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
1267 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
1268 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
1269 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
1273 * Super precision for pi/2 in 24-bit chunks, for use in
1274 * {@link #remPiOver2()}.
1276 private static final double PI_OVER_TWO
[] = {
1277 1.570796251296997, // Long bits 0x3ff921fb40000000L.
1278 7.549789415861596e-8, // Long bits 0x3e74442d00000000L.
1279 5.390302529957765e-15, // Long bits 0x3cf8469880000000L.
1280 3.282003415807913e-22, // Long bits 0x3b78cc5160000000L.
1281 1.270655753080676e-29, // Long bits 0x39f01b8380000000L.
1282 1.2293330898111133e-36, // Long bits 0x387a252040000000L.
1283 2.7337005381646456e-44, // Long bits 0x36e3822280000000L.
1284 2.1674168387780482e-51, // Long bits 0x3569f31d00000000L.
1288 * More constants related to pi, used in {@link #remPiOver2()} and
1291 private static final double
1292 PI_L
= 1.2246467991473532e-16, // Long bits 0x3ca1a62633145c07L.
1293 PIO2_1
= 1.5707963267341256, // Long bits 0x3ff921fb54400000L.
1294 PIO2_1L
= 6.077100506506192e-11, // Long bits 0x3dd0b4611a626331L.
1295 PIO2_2
= 6.077100506303966e-11, // Long bits 0x3dd0b4611a600000L.
1296 PIO2_2L
= 2.0222662487959506e-21, // Long bits 0x3ba3198a2e037073L.
1297 PIO2_3
= 2.0222662487111665e-21, // Long bits 0x3ba3198a2e000000L.
1298 PIO2_3L
= 8.4784276603689e-32; // Long bits 0x397b839a252049c1L.
1301 * Natural log and square root constants, for calculation of
1302 * {@link #exp(double)}, {@link #log(double)} and
1303 * {@link #power(double, double)}. CP is 2/(3*ln(2)).
1305 private static final double
1306 SQRT_1_5
= 1.224744871391589, // Long bits 0x3ff3988e1409212eL.
1307 SQRT_2
= 1.4142135623730951, // Long bits 0x3ff6a09e667f3bcdL.
1308 SQRT_3
= 1.7320508075688772, // Long bits 0x3ffbb67ae8584caaL.
1309 EXP_LIMIT_H
= 709.782712893384, // Long bits 0x40862e42fefa39efL.
1310 EXP_LIMIT_L
= -745.1332191019411, // Long bits 0xc0874910d52d3051L.
1311 CP
= 0.9617966939259756, // Long bits 0x3feec709dc3a03fdL.
1312 CP_H
= 0.9617967009544373, // Long bits 0x3feec709e0000000L.
1313 CP_L
= -7.028461650952758e-9, // Long bits 0xbe3e2fe0145b01f5L.
1314 LN2
= 0.6931471805599453, // Long bits 0x3fe62e42fefa39efL.
1315 LN2_H
= 0.6931471803691238, // Long bits 0x3fe62e42fee00000L.
1316 LN2_L
= 1.9082149292705877e-10, // Long bits 0x3dea39ef35793c76L.
1317 INV_LN2
= 1.4426950408889634, // Long bits 0x3ff71547652b82feL.
1318 INV_LN2_H
= 1.4426950216293335, // Long bits 0x3ff7154760000000L.
1319 INV_LN2_L
= 1.9259629911266175e-8; // Long bits 0x3e54ae0bf85ddf44L.
1322 * Constants for computing {@link #log(double)}.
1324 private static final double
1325 LG1
= 0.6666666666666735, // Long bits 0x3fe5555555555593L.
1326 LG2
= 0.3999999999940942, // Long bits 0x3fd999999997fa04L.
1327 LG3
= 0.2857142874366239, // Long bits 0x3fd2492494229359L.
1328 LG4
= 0.22222198432149784, // Long bits 0x3fcc71c51d8e78afL.
1329 LG5
= 0.1818357216161805, // Long bits 0x3fc7466496cb03deL.
1330 LG6
= 0.15313837699209373, // Long bits 0x3fc39a09d078c69fL.
1331 LG7
= 0.14798198605116586; // Long bits 0x3fc2f112df3e5244L.
1334 * Constants for computing {@link #pow(double, double)}. L and P are
1335 * coefficients for series; OVT is -(1024-log2(ovfl+.5ulp)); and DP is ???.
1336 * The P coefficients also calculate {@link #exp(double)}.
1338 private static final double
1339 L1
= 0.5999999999999946, // Long bits 0x3fe3333333333303L.
1340 L2
= 0.4285714285785502, // Long bits 0x3fdb6db6db6fabffL.
1341 L3
= 0.33333332981837743, // Long bits 0x3fd55555518f264dL.
1342 L4
= 0.272728123808534, // Long bits 0x3fd17460a91d4101L.
1343 L5
= 0.23066074577556175, // Long bits 0x3fcd864a93c9db65L.
1344 L6
= 0.20697501780033842, // Long bits 0x3fca7e284a454eefL.
1345 P1
= 0.16666666666666602, // Long bits 0x3fc555555555553eL.
1346 P2
= -2.7777777777015593e-3, // Long bits 0xbf66c16c16bebd93L.
1347 P3
= 6.613756321437934e-5, // Long bits 0x3f11566aaf25de2cL.
1348 P4
= -1.6533902205465252e-6, // Long bits 0xbebbbd41c5d26bf1L.
1349 P5
= 4.1381367970572385e-8, // Long bits 0x3e66376972bea4d0L.
1350 DP_H
= 0.5849624872207642, // Long bits 0x3fe2b80340000000L.
1351 DP_L
= 1.350039202129749e-8, // Long bits 0x3e4cfdeb43cfd006L.
1352 OVT
= 8.008566259537294e-17; // Long bits 0x3c971547652b82feL.
1355 * Coefficients for computing {@link #sin(double)}.
1357 private static final double
1358 S1
= -0.16666666666666632, // Long bits 0xbfc5555555555549L.
1359 S2
= 8.33333333332249e-3, // Long bits 0x3f8111111110f8a6L.
1360 S3
= -1.984126982985795e-4, // Long bits 0xbf2a01a019c161d5L.
1361 S4
= 2.7557313707070068e-6, // Long bits 0x3ec71de357b1fe7dL.
1362 S5
= -2.5050760253406863e-8, // Long bits 0xbe5ae5e68a2b9cebL.
1363 S6
= 1.58969099521155e-10; // Long bits 0x3de5d93a5acfd57cL.
1366 * Coefficients for computing {@link #cos(double)}.
1368 private static final double
1369 C1
= 0.0416666666666666, // Long bits 0x3fa555555555554cL.
1370 C2
= -1.388888888887411e-3, // Long bits 0xbf56c16c16c15177L.
1371 C3
= 2.480158728947673e-5, // Long bits 0x3efa01a019cb1590L.
1372 C4
= -2.7557314351390663e-7, // Long bits 0xbe927e4f809c52adL.
1373 C5
= 2.087572321298175e-9, // Long bits 0x3e21ee9ebdb4b1c4L.
1374 C6
= -1.1359647557788195e-11; // Long bits 0xbda8fae9be8838d4L.
1377 * Coefficients for computing {@link #tan(double)}.
1379 private static final double
1380 T0
= 0.3333333333333341, // Long bits 0x3fd5555555555563L.
1381 T1
= 0.13333333333320124, // Long bits 0x3fc111111110fe7aL.
1382 T2
= 0.05396825397622605, // Long bits 0x3faba1ba1bb341feL.
1383 T3
= 0.021869488294859542, // Long bits 0x3f9664f48406d637L.
1384 T4
= 8.8632398235993e-3, // Long bits 0x3f8226e3e96e8493L.
1385 T5
= 3.5920791075913124e-3, // Long bits 0x3f6d6d22c9560328L.
1386 T6
= 1.4562094543252903e-3, // Long bits 0x3f57dbc8fee08315L.
1387 T7
= 5.880412408202641e-4, // Long bits 0x3f4344d8f2f26501L.
1388 T8
= 2.464631348184699e-4, // Long bits 0x3f3026f71a8d1068L.
1389 T9
= 7.817944429395571e-5, // Long bits 0x3f147e88a03792a6L.
1390 T10
= 7.140724913826082e-5, // Long bits 0x3f12b80f32f0a7e9L.
1391 T11
= -1.8558637485527546e-5, // Long bits 0xbef375cbdb605373L.
1392 T12
= 2.590730518636337e-5; // Long bits 0x3efb2a7074bf7ad4L.
1395 * Coefficients for computing {@link #asin(double)} and
1396 * {@link #acos(double)}.
1398 private static final double
1399 PS0
= 0.16666666666666666, // Long bits 0x3fc5555555555555L.
1400 PS1
= -0.3255658186224009, // Long bits 0xbfd4d61203eb6f7dL.
1401 PS2
= 0.20121253213486293, // Long bits 0x3fc9c1550e884455L.
1402 PS3
= -0.04005553450067941, // Long bits 0xbfa48228b5688f3bL.
1403 PS4
= 7.915349942898145e-4, // Long bits 0x3f49efe07501b288L.
1404 PS5
= 3.479331075960212e-5, // Long bits 0x3f023de10dfdf709L.
1405 QS1
= -2.403394911734414, // Long bits 0xc0033a271c8a2d4bL.
1406 QS2
= 2.0209457602335057, // Long bits 0x40002ae59c598ac8L.
1407 QS3
= -0.6882839716054533, // Long bits 0xbfe6066c1b8d0159L.
1408 QS4
= 0.07703815055590194; // Long bits 0x3fb3b8c5b12e9282L.
1411 * Coefficients for computing {@link #atan(double)}.
1413 private static final double
1414 ATAN_0_5H
= 0.4636476090008061, // Long bits 0x3fddac670561bb4fL.
1415 ATAN_0_5L
= 2.2698777452961687e-17, // Long bits 0x3c7a2b7f222f65e2L.
1416 ATAN_1_5H
= 0.982793723247329, // Long bits 0x3fef730bd281f69bL.
1417 ATAN_1_5L
= 1.3903311031230998e-17, // Long bits 0x3c7007887af0cbbdL.
1418 AT0
= 0.3333333333333293, // Long bits 0x3fd555555555550dL.
1419 AT1
= -0.19999999999876483, // Long bits 0xbfc999999998ebc4L.
1420 AT2
= 0.14285714272503466, // Long bits 0x3fc24924920083ffL.
1421 AT3
= -0.11111110405462356, // Long bits 0xbfbc71c6fe231671L.
1422 AT4
= 0.09090887133436507, // Long bits 0x3fb745cdc54c206eL.
1423 AT5
= -0.0769187620504483, // Long bits 0xbfb3b0f2af749a6dL.
1424 AT6
= 0.06661073137387531, // Long bits 0x3fb10d66a0d03d51L.
1425 AT7
= -0.058335701337905735, // Long bits 0xbfadde2d52defd9aL.
1426 AT8
= 0.049768779946159324, // Long bits 0x3fa97b4b24760debL.
1427 AT9
= -0.036531572744216916, // Long bits 0xbfa2b4442c6a6c2fL.
1428 AT10
= 0.016285820115365782; // Long bits 0x3f90ad3ae322da11L.
1431 * Helper function for reducing an angle to a multiple of pi/2 within
1434 * @param x the angle; not infinity or NaN, and outside pi/4
1435 * @param y an array of 2 doubles modified to hold the remander x % pi/2
1436 * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4],
1437 * 1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4]
1439 private static int remPiOver2(double x
, double[] y
)
1441 boolean negative
= x
< 0;
1445 if (Configuration
.DEBUG
&& (x
<= PI
/ 4 || x
!= x
1446 || x
== Double
.POSITIVE_INFINITY
))
1447 throw new InternalError("Assertion failure");
1448 if (x
< 3 * PI
/ 4) // If |x| is small.
1451 if ((float) x
!= (float) (PI
/ 2)) // 33+53 bit pi is good enough.
1454 y
[1] = z
- y
[0] - PIO2_1L
;
1456 else // Near pi/2, use 33+33+53 bit pi.
1460 y
[1] = z
- y
[0] - PIO2_2L
;
1464 else if (x
<= TWO_20
* PI
/ 2) // Medium size.
1466 n
= (int) (2 / PI
* x
+ 0.5);
1468 double w
= n
* PIO2_1L
; // First round good to 85 bits.
1470 if (n
>= 32 || (float) x
== (float) (w
))
1472 if (x
/ y
[0] >= TWO_16
) // Second iteration, good to 118 bits.
1477 w
= n
* PIO2_2L
- (t
- z
- w
);
1479 if (x
/ y
[0] >= TWO_49
) // Third iteration, 151 bits accuracy.
1484 w
= n
* PIO2_3L
- (t
- z
- w
);
1489 y
[1] = z
- y
[0] - w
;
1493 // All other (large) arguments.
1494 int e0
= (int) (Double
.doubleToLongBits(x
) >> 52) - 1046;
1495 z
= scale(x
, -e0
); // e0 = ilogb(z) - 23.
1496 double[] tx
= new double[3];
1497 for (int i
= 0; i
< 2; i
++)
1500 z
= (z
- tx
[i
]) * TWO_24
;
1506 n
= remPiOver2(tx
, y
, e0
, nx
);
1518 * Helper function for reducing an angle to a multiple of pi/2 within
1521 * @param x the positive angle, broken into 24-bit chunks
1522 * @param y an array of 2 doubles modified to hold the remander x % pi/2
1523 * @param e0 the exponent of x[0]
1524 * @param nx the last index used in x
1525 * @return the quadrant of the result, mod 4: 0: [-pi/4, pi/4],
1526 * 1: [pi/4, 3*pi/4], 2: [3*pi/4, 5*pi/4], 3: [-3*pi/4, -pi/4]
1528 private static int remPiOver2(double[] x
, double[] y
, int e0
, int nx
)
1535 int[] iq
= new int[20];
1536 double[] f
= new double[20];
1537 double[] q
= new double[20];
1538 boolean recompute
= false;
1540 // Initialize jk, jz, jv, q0; note that 3>q0.
1543 int jv
= max((e0
- 3) / 24, 0);
1544 int q0
= e0
- 24 * (jv
+ 1);
1546 // Set up f[0] to f[nx+jk] where f[nx+jk] = TWO_OVER_PI[jv+jk].
1549 for (i
= 0; i
<= m
; i
++, j
++)
1550 f
[i
] = (j
< 0) ?
0 : TWO_OVER_PI
[j
];
1552 // Compute q[0],q[1],...q[jk].
1553 for (i
= 0; i
<= jk
; i
++)
1555 for (j
= 0, fw
= 0; j
<= nx
; j
++)
1556 fw
+= x
[j
] * f
[nx
+ i
- j
];
1562 // Distill q[] into iq[] reversingly.
1563 for (i
= 0, j
= jz
, z
= q
[jz
]; j
> 0; i
++, j
--)
1565 fw
= (int) (1 / TWO_24
* z
);
1566 iq
[i
] = (int) (z
- TWO_24
* fw
);
1572 z
-= 8 * floor(z
* 0.125); // Trim off integer >= 8.
1576 if (q0
> 0) // Need iq[jz-1] to determine n.
1578 i
= iq
[jz
- 1] >> (24 - q0
);
1580 iq
[jz
- 1] -= i
<< (24 - q0
);
1581 ih
= iq
[jz
- 1] >> (23 - q0
);
1584 ih
= iq
[jz
- 1] >> 23;
1588 if (ih
> 0) // If q > 0.5.
1592 for (i
= 0; i
< jz
; i
++) // Compute 1-q.
1600 iq
[i
] = 0x1000000 - j
;
1604 iq
[i
] = 0xffffff - j
;
1608 case 1: // Rare case: chance is 1 in 12 for non-default.
1609 iq
[jz
- 1] &= 0x7fffff;
1612 iq
[jz
- 1] &= 0x3fffff;
1622 // Check if recomputation is needed.
1626 for (i
= jz
- 1; i
>= jk
; i
--)
1628 if (j
== 0) // Need recomputation.
1631 for (k
= 1; iq
[jk
- k
] == 0; k
++); // k = no. of terms needed.
1633 for (i
= jz
+ 1; i
<= jz
+ k
; i
++) // Add q[jz+1] to q[jz+k].
1635 f
[nx
+ i
] = TWO_OVER_PI
[jv
+ i
];
1636 for (j
= 0, fw
= 0; j
<= nx
; j
++)
1637 fw
+= x
[j
] * f
[nx
+ i
- j
];
1647 // Chop off zero terms.
1658 else // Break z into 24-bit if necessary.
1663 fw
= (int) (1 / TWO_24
* z
);
1664 iq
[jz
] = (int) (z
- TWO_24
* fw
);
1673 // Convert integer "bit" chunk to floating-point value.
1675 for (i
= jz
; i
>= 0; i
--)
1681 // Compute PI_OVER_TWO[0,...,jk]*q[jz,...,0].
1682 double[] fq
= new double[20];
1683 for (i
= jz
; i
>= 0; i
--)
1686 for (int k
= 0; k
<= jk
&& k
<= jz
- i
; k
++)
1687 fw
+= PI_OVER_TWO
[k
] * q
[i
+ k
];
1691 // Compress fq[] into y[].
1693 for (i
= jz
; i
>= 0; i
--)
1695 y
[0] = (ih
== 0) ? fw
: -fw
;
1697 for (i
= 1; i
<= jz
; i
++)
1699 y
[1] = (ih
== 0) ? fw
: -fw
;
1704 * Helper method for scaling a double by a power of 2.
1706 * @param x the double
1707 * @param n the scale; |n| < 2048
1710 private static double scale(double x
, int n
)
1712 if (Configuration
.DEBUG
&& abs(n
) >= 2048)
1713 throw new InternalError("Assertion failure");
1714 if (x
== 0 || x
== Double
.NEGATIVE_INFINITY
1715 || ! (x
< Double
.POSITIVE_INFINITY
) || n
== 0)
1717 long bits
= Double
.doubleToLongBits(x
);
1718 int exp
= (int) (bits
>> 52) & 0x7ff;
1719 if (exp
== 0) // Subnormal x.
1722 exp
= ((int) (Double
.doubleToLongBits(x
) >> 52) & 0x7ff) - 54;
1725 if (exp
> 0x7fe) // Overflow.
1726 return Double
.POSITIVE_INFINITY
* x
;
1727 if (exp
> 0) // Normal.
1728 return Double
.longBitsToDouble((bits
& 0x800fffffffffffffL
)
1729 | ((long) exp
<< 52));
1731 return 0 * x
; // Underflow.
1732 exp
+= 54; // Subnormal result.
1733 x
= Double
.longBitsToDouble((bits
& 0x800fffffffffffffL
)
1734 | ((long) exp
<< 52));
1735 return x
* (1 / TWO_54
);
1739 * Helper trig function; computes sin in range [-pi/4, pi/4].
1741 * @param x angle within about pi/4
1742 * @param y tail of x, created by remPiOver2
1745 private static double sin(double x
, double y
)
1747 if (Configuration
.DEBUG
&& abs(x
+ y
) > 0.7854)
1748 throw new InternalError("Assertion failure");
1749 if (abs(x
) < 1 / TWO_27
)
1750 return x
; // If |x| ~< 2**-27, already know answer.
1754 double r
= S2
+ z
* (S3
+ z
* (S4
+ z
* (S5
+ z
* S6
)));
1756 return x
+ v
* (S1
+ z
* r
);
1757 return x
- ((z
* (0.5 * y
- v
* r
) - y
) - v
* S1
);
1761 * Helper trig function; computes cos in range [-pi/4, pi/4].
1763 * @param x angle within about pi/4
1764 * @param y tail of x, created by remPiOver2
1767 private static double cos(double x
, double y
)
1769 if (Configuration
.DEBUG
&& abs(x
+ y
) > 0.7854)
1770 throw new InternalError("Assertion failure");
1773 return 1; // If |x| ~< 2**-27, already know answer.
1776 double r
= z
* (C1
+ z
* (C2
+ z
* (C3
+ z
* (C4
+ z
* (C5
+ z
* C6
)))));
1779 return 1 - (0.5 * z
- (z
* r
- x
* y
));
1781 double qx
= (x
> 0.78125) ?
0.28125 : (x
* 0.25);
1782 return 1 - qx
- ((0.5 * z
- qx
) - (z
* r
- x
* y
));
1786 * Helper trig function; computes tan in range [-pi/4, pi/4].
1788 * @param x angle within about pi/4
1789 * @param y tail of x, created by remPiOver2
1790 * @param invert true iff -1/tan should be returned instead
1793 private static double tan(double x
, double y
, boolean invert
)
1795 // PI/2 is irrational, so no double is a perfect multiple of it.
1796 if (Configuration
.DEBUG
&& (abs(x
+ y
) > 0.7854 || (x
== 0 && invert
)))
1797 throw new InternalError("Assertion failure");
1798 boolean negative
= x
< 0;
1804 if (x
< 1 / TWO_28
) // If |x| ~< 2**-28, already know answer.
1805 return (negative ?
-1 : 1) * (invert ?
-1 / x
: x
);
1809 boolean large
= x
>= 0.6744;
1819 // Break x**5*(T1+x**2*T2+...) into
1820 // x**5(T1+x**4*T3+...+x**20*T11)
1821 // + x**5(x**2*(T2+x**4*T4+...+x**22*T12)).
1822 double r
= T1
+ w
* (T3
+ w
* (T5
+ w
* (T7
+ w
* (T9
+ w
* T11
))));
1823 double v
= z
* (T2
+ w
* (T4
+ w
* (T6
+ w
* (T8
+ w
* (T10
+ w
* T12
)))));
1825 r
= y
+ z
* (s
* (r
+ v
) + y
);
1830 v
= invert ?
-1 : 1;
1831 return (negative ?
-1 : 1) * (v
- 2 * (x
- (w
* w
/ (w
+ v
) - r
)));
1836 // Compute -1.0/(x+r) accurately.
1840 double t
= (float) a
;
1841 return t
+ a
* (1 + t
* z
+ t
* v
);