1 /*
2 * -------------------------------------------------------------------------
3 * $Id: JMath.java 1712 2004-04-28 17:36:13Z joewhaley $
4 * -------------------------------------------------------------------------
5 * Copyright (c) 1999 Visual Numerics Inc. All Rights Reserved.
6 *
7 * Permission to use, copy, modify, and distribute this software is freely
8 * granted by Visual Numerics, Inc., provided that the copyright notice
9 * above and the following warranty disclaimer are preserved in human
10 * readable form.
11 *
12 * Because this software is licenses free of charge, it is provided
13 * "AS IS", with NO WARRANTY. TO THE EXTENT PERMITTED BY LAW, VNI
14 * DISCLAIMS ALL WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
15 * TO ITS PERFORMANCE, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
16 * VNI WILL NOT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING OUT OF THE USE
17 * OF OR INABILITY TO USE THIS SOFTWARE, INCLUDING BUT NOT LIMITED TO DIRECT,
18 * INDIRECT, SPECIAL, CONSEQUENTIAL, PUNITIVE, AND EXEMPLARY DAMAGES, EVEN
19 * IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
20 *
21 *
22 * This Java code is based on C code in the package fdlibm,
23 * which can be obtained from www.netlib.org.
24 * The original fdlibm C code contains the following notice.
25 *
26 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
27 *
28 * Developed at SunSoft, a Sun Microsystems, Inc. business.
29 * Permission to use, copy, modify, and distribute this
30 * software is freely granted, provided that this notice
31 * is preserved.
32 *
33 *--------------------------------------------------------------------------
34 */
35 package joeq.Support;
36
37 import java.util.Random;
38
39 /*
40 * Pure Java implementation of the standard java.lang.Math class.
41 * This Java code is based on C code in the package fdlibm,
42 * which can be obtained from www.netlib.org.
43 *
44 * @author Sun Microsystems (original C code in fdlibm)
45 * @author John F. Brophy (translated from C to Java)
46 */
47 public final class JMath
48 {
49 public static final double PI = Double.longBitsToDouble(0x400921fb54442d18L); /* 3.14159265358979323846 */
50 public static final double E = 2.7182818284590452354;
51
52 private static Random random;
53
54
55 /***
56 * Returns the absolute value of its argument.
57 * @param x The argument, an integer.
58 * @return Returns |x|.
59 */
60 public static strictfp int abs(int x)
61 {
62 return ((x < 0) ? -x : x);
63 }
64
65
66 /***
67 * Returns the absolute value of its argument.
68 * @param x The argument, a long.
69 * @return Returns |x|.
70 */
71 public static strictfp long abs(long x)
72 {
73 return ((x < 0L) ? -x : x);
74 }
75
76
77 /***
78 * Returns the absolute value of its argument.
79 * @param x The argument, a float.
80 * @return Returns |x|.
81 */
82 public static strictfp float abs(float x)
83 {
84 return ((x <= 0.0f) ? 0.0f-x : x);
85 }
86
87
88 /***
89 * Returns the absolute value of its argument.
90 * @param x The argument, a double.
91 * @return Returns |x|.
92 */
93 public static strictfp double abs(double x)
94 {
95 return ((x <= 0.0) ? 0.0-x : x);
96 }
97
98
99 /***
100 * Returns the smaller of its two arguments.
101 * @param x The first argument, an integer.
102 * @param y The second argument, an integer.
103 * @return Returns the smaller of x and y.
104 */
105 public static strictfp int min(int x, int y)
106 {
107 return ((x < y) ? x : y);
108 }
109
110
111 /***
112 * Returns the smaller of its two arguments.
113 * @param x The first argument, a long.
114 * @param y The second argument, a long.
115 * @return Returns the smaller of x and y.
116 */
117 public static strictfp long min(long x, long y)
118 {
119 return ((x < y) ? x : y);
120 }
121
122
123 /***
124 * Returns the smaller of its two arguments.
125 * @param x The first argument, a float.
126 * @param y The second argument, a float.
127 * @return Returns the smaller of x and y.
128 * This function considers -0.0f to
129 * be less than 0.0f.
130 */
131 public static strictfp float min(float x, float y)
132 {
133 if (Float.isNaN(x)) return x;
134 float ans = ((x <= y) ? x : y);
135 if (ans == 0.0f &&
136 Float.floatToIntBits(y) == 0x80000000)
137 ans = y;
138 return ans;
139 }
140
141
142 /***
143 * Returns the smaller of its two arguments.
144 * @param x The first argument, a double.
145 * @param y The second argument, a double.
146 * @return Returns the smaller of x and y.
147 * This function considers -0.0 to
148 * be less than 0.0.
149 */ public static strictfp double min(double x, double y)
150 {
151 if (Double.isNaN(x)) return x;
152 double ans = ((x <= y) ? x : y);
153 if (x == 0.0 && y == 0.0 &&
154 Double.doubleToLongBits(y) == 0x8000000000000000L)
155 ans = y;
156 return ans;
157 }
158
159
160 /***
161 * Returns the larger of its two arguments.
162 * @param x The first argument, an integer.
163 * @param y The second argument, an integer.
164 * @return Returns the larger of x and y.
165 */
166 public static strictfp int max(int x, int y)
167 {
168 return ((x > y) ? x : y);
169 }
170
171
172 /***
173 * Returns the larger of its two arguments.
174 * @param x The first argument, a long.
175 * @param y The second argument, a long.
176 * @return Returns the larger of x and y.
177 */
178 public static strictfp long max(long x, long y)
179 {
180 return ((x > y) ? x : y);
181 }
182
183
184 /***
185 * Returns the larger of its two arguments.
186 * @param x The first argument, a float.
187 * @param y The second argument, a float.
188 * @return Returns the larger of x and y.
189 * This function considers -0.0f to
190 * be less than 0.0f.
191 */
192 public static strictfp float max(float x, float y)
193 {
194 if (Float.isNaN(x)) return x;
195 float ans = ((x >= y) ? x : y);
196 if (ans == 0.0f && Float.floatToIntBits(x) == 0x80000000)
197 ans = y;
198 return ans;
199 }
200
201
202 /***
203 * Returns the larger of its two arguments.
204 * @param x The first argument, a double.
205 * @param y The second argument, a double.
206 * @return Returns the larger of x and y.
207 * This function considers -0.0 to
208 * be less than 0.0.
209 */
210 public static strictfp double max(double x, double y)
211 {
212 if (Double.isNaN(x)) return x;
213 double ans = ((x >= y) ? x : y);
214 if (x == 0.0 && y == 0.0 && Double.doubleToLongBits(x) == 0x8000000000000000L)
215 ans = y;
216 return ans;
217 }
218
219
220 /***
221 * Returns the integer closest to the arguments.
222 * @param x The argument, a float.
223 * @return Returns the integer closest to x.
224 */
225 public static strictfp int round(float x)
226 {
227 return (int)floor(x+0.5f);
228 }
229
230
231 /***
232 * Returns the long closest to the arguments.
233 * @param x The argument, a double.
234 * @return Returns the long closest to x.
235 */
236 public static strictfp long round(double x)
237 {
238 return (long)floor(x+0.5);
239 }
240
241
242 /***
243 * Returns the random number.
244 * @return Returns a random number from a uniform distribution.
245 */
246 public static synchronized strictfp double random()
247 {
248 if (random == null) random = new Random();
249 return random.nextDouble();
250 }
251
252
253 /*
254 * This following code is derived from fdlibm, which contained
255 * the following notice.
256 * ====================================================
257 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
258 *
259 * Developed at SunSoft, a Sun Microsystems, Inc. business.
260 * Permission to use, copy, modify, and distribute this
261 * software is freely granted, provided that this notice
262 * is preserved.
263 * ====================================================
264 *
265 * Constants:
266 * The hexadecimal values are the intended ones for the following
267 * constants. The decimal values may be used, provided that the
268 * compiler will convert from decimal to binary accurately enough
269 * to produce the hexadecimal values shown.
270 */
271
272 private static final double huge = 1.0e+300;
273 private static final double tiny = 1.0e-300;
274
275
276 /***
277 * Returns the value of its argument rounded toward
278 * positive infinity to an integral value.
279 * @param x The argument, a double.
280 * @return Returns the smallest double, not less than x,
281 * that is an integral value.
282 */
283 public static double ceil(double x)
284 {
285 int exp;
286 long ix;
287
288 if (x == 0) return x;
289
290 ix = Double.doubleToLongBits(x);
291 //sign = (int)((ix>>63)&1);
292 exp = ((int)(ix>>52)&0x7ff) - 0x3ff;
293
294 if (exp < 0) {
295 if (x < 0.0)
296 return NEGATIVE_ZERO;
297 else if (x == 0.0)
298 return x;
299 else
300 return 1.0;
301 } else if (exp < 53) {
302 long mask = (0x000fffffffffffffL>>>exp);
303 if ((mask&ix) == 0) return x; // x is integral
304 if (x > 0.0)
305 ix += 0x0010000000000000L>>exp;
306 ix = ix & (~mask);
307 } else if (exp == 1024) { // infinity
308 return x;
309 }
310
311 return Double.longBitsToDouble(ix);
312 }
313
314 /***
315 * Returns the value of its argument rounded toward
316 * negative infinity to an integral value.
317 * @param x The argument, a double.
318 * @return Returns the smallest double, not greater than x,
319 * that is an integral value.
320 */
321 public static double floor(double x)
322 {
323 int exp;
324 long ix;
325
326 if (x == 0) return x;
327
328 ix = Double.doubleToLongBits(x);
329 //sign = (int)((ix>>63)&1);
330 exp = ((int)(ix>>52)&0x7ff) - 0x3ff;
331
332 if (exp < 0) {
333 if (x < 0.0)
334 return -1.0;
335 else if (x == 0.0)
336 return x;
337 else
338 return 0.0;
339 } else if (exp < 53) {
340 long mask = (0x000fffffffffffffL>>>exp);
341 if ((mask&ix) == 0) return x; // x is integral
342 if (x < 0.0)
343 ix += 0x0010000000000000L>>exp;
344 ix = ix & (~mask);
345 } else if (exp == 1024) { // infinity
346 return x;
347 }
348
349 return Double.longBitsToDouble(ix);
350 }
351
352
353
354 private static final double TWO52[] = {
355 Double.longBitsToDouble(0x4330000000000000L), /* 4.50359962737049600000e+15 */
356 Double.longBitsToDouble(0xc330000000000000L) /* -4.50359962737049600000e+15 */
357 };
358
359 private static final double NEGATIVE_ZERO = Double.longBitsToDouble(0x8000000000000000L);
360
361
362 /***
363 * Returns the value of its argument rounded toward
364 * the closest integral value.
365 * @param x The argument, a double.
366 * @return Returns the double closest to x
367 * that is an integral value.
368 */
369 public static double rint(double x)
370 {
371 int exp, sign;
372 long ix;
373 double w;
374
375 if (x == 0) return x;
376
377 ix = Double.doubleToLongBits(x);
378 sign = (int)((ix>>63)&1);
379 exp = ((int)(ix>>52)&0x7ff) - 0x3ff;
380
381 if (exp < 0) {
382 if (x < -0.5)
383 return -1.0;
384 else if (x > 0.5)
385 return 1.0;
386 else if (sign == 0)
387 return 0.0;
388 else
389 return NEGATIVE_ZERO;
390 } else if (exp < 53) {
391 long mask = (0x000fffffffffffffL>>>exp);
392 if ((mask&ix) == 0) return x; // x is integral
393 } else if (exp == 1024) { // infinity
394 return x;
395 }
396
397 x = Double.longBitsToDouble(ix);
398 w = TWO52[sign] + x;
399 return w - TWO52[sign];
400 }
401
402
403 /***
404 * Returns x REM p = x - [x/p]*p as if in infinite
405 * precise arithmetic, where [x/p] is the (infinite bit)
406 * integer nearest x/p (in half way case choose the even one).
407 * @param x The dividend.
408 * @param p The divisor.
409 * @return The remainder computed according to the IEEE 754 standard.
410 */
411 public static double IEEEremainder(double x, double p)
412 {
413 int hx,hp;
414 int sx,lx,lp;// unsigned
415 double p_half;
416
417 hx = __HI(x); /* high word of x */
418 lx = __LO(x); /* low word of x */
419 hp = __HI(p); /* high word of p */
420 lp = __LO(p); /* low word of p */
421 sx = hx&0x80000000;
422 hp &= 0x7fffffff;
423 hx &= 0x7fffffff;
424
425 /* purge off exception values */
426 if ((hp|lp) == 0) return (x*p)/(x*p); /* p = 0 */
427
428 if ((hx >= 0x7ff00000) || /* x not finite */
429 ((hp >= 0x7ff00000)&& /* p is NaN */
430 (((hp-0x7ff00000)|lp)!=0)))
431 return (x*p)/(x*p);
432
433 if (hp <= 0x7fdfffff) x = x%(p+p); /* now x < 2p */
434
435 if (((hx-hp)|(lx-lp)) == 0) return zero*x;
436
437 x = abs(x);
438 p = abs(p);
439 if (hp < 0x00200000) {
440 if (x+x > p) {
441 x -= p;
442 if (x+x >= p) x -= p;
443 }
444 } else {
445 p_half = 0.5*p;
446 if (x > p_half) {
447 x -= p;
448 if (x >= p_half)
449 x -= p;
450 }
451 }
452 lx = __HI(x);
453 lx ^= sx;
454 return setHI(x, lx);
455 }
456
457
458
459 /* sqrt(x)
460 * Return correctly rounded sqrt.
461 * ------------------------------------------
462 * | Use the hardware sqrt if you have one |
463 * ------------------------------------------
464 * Method:
465 * Bit by bit method using integer arithmetic. (Slow, but portable)
466 * 1. Normalization
467 * Scale x to y in [1,4) with even powers of 2:
468 * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
469 * sqrt(x) = 2^k * sqrt(y)
470 * 2. Bit by bit computation
471 * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
472 * i 0
473 * i+1 2
474 * s = 2*q , and y = 2 * ( y - q ). (1)
475 * i i i i
476 *
477 * To compute q from q , one checks whether
478 * i+1 i
479 *
480 * -(i+1) 2
481 * (q + 2 ) <= y. (2)
482 * i
483 * -(i+1)
484 * If (2) is false, then q = q ; otherwise q = q + 2 .
485 * i+1 i i+1 i
486 *
487 * With some algebric manipulation, it is not difficult to see
488 * that (2) is equivalent to
489 * -(i+1)
490 * s + 2 <= y (3)
491 * i i
492 *
493 * The advantage of (3) is that s and y can be computed by
494 * i i
495 * the following recurrence formula:
496 * if (3) is false
497 *
498 * s = s , y = y ; (4)
499 * i+1 i i+1 i
500 *
501 * otherwise,
502 * -i -(i+1)
503 * s = s + 2 , y = y - s - 2 (5)
504 * i+1 i i+1 i i
505 *
506 * One may easily use induction to prove (4) and (5).
507 * Note. Since the left hand side of (3) contain only i+2 bits,
508 * it does not necessary to do a full (53-bit) comparison
509 * in (3).
510 * 3. Final rounding
511 * After generating the 53 bits result, we compute one more bit.
512 * Together with the remainder, we can decide whether the
513 * result is exact, bigger than 1/2ulp, or less than 1/2ulp
514 * (it will never equal to 1/2ulp).
515 * The rounding mode can be detected by checking whether
516 * huge + tiny is equal to huge, and whether huge - tiny is
517 * equal to huge for some floating point number "huge" and "tiny".
518 *
519 *
520 * Special cases:
521 * sqrt(+-0) = +-0 ... exact
522 * sqrt(inf) = inf
523 * sqrt(-ve) = NaN ... with invalid signal
524 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
525 */
526
527 /***
528 * Returns the square root of its argument.
529 * @param x The argument, a double.
530 * @return Returns the square root of x.
531 */
532 public static double sqrt(double x)
533 {
534 long ix = Double.doubleToLongBits(x);
535
536 /* take care of Inf and NaN */
537 if ((ix&0x7ff0000000000000L)==0x7ff0000000000000L) {
538 /* sqrt(NaN)=NaN, sqrt(+inf)=+inf sqrt(-inf)=sNaN */
539 return x*x + x;
540 }
541
542 /* take care of zero */
543 if (x < 0.0) {
544 return Double.NaN;
545 } else if (x == 0.0) {
546 return x; /* sqrt(+-0) = +-0 */
547 }
548
549 /* normalize x */
550 long m = (ix>>52);
551 ix &= 0x000fffffffffffffL;
552
553 /* add implicit bit, if not sub-normal */
554 if (m != 0) ix |= 0x0010000000000000L;
555
556 m -= 1023L; /* unbias exponent */
557 if ((m&1) != 0) { /* odd m, double x to make it even */
558 ix += ix;
559 }
560 m >>= 1; /* m = [m/2] */
561 m += 1023L;
562
563 /* generate sqrt(x) bit by bit */
564 ix += ix;
565 long q = 0L; /* q = sqrt(x) */
566 long s = 0L;
567 long r = 0x0020000000000000L; /* r = moving bit from right to left */
568
569 while (r != 0) {
570 long t = s + r;
571 if (t <= ix) {
572 s = t + r;
573 ix -= t;
574 q += r;
575 }
576 ix += ix;
577 r >>= 1;
578 }
579
580 /* round */
581 if (ix != 0) q += (q&1L);
582
583 /* assemble result */
584 ix = (m << 52) | (0x000fffffffffffffL&(q>>1));
585 return Double.longBitsToDouble(ix);
586 }
587
588
589
590 private static final double halF[] = {0.5, -0.5};
591 private static final double twom1000 =
592 Double.longBitsToDouble(0x0170000000000000L); /* 2**-1000=9.33263618503218878990e-302 */
593
594 private static final double o_threshold = Double.longBitsToDouble(0x40862e42fefa39efL); /* 7.09782712893383973096e+02 */
595 private static final double u_threshold = Double.longBitsToDouble(0xc0874910d52d3051L); /* -7.45133219101941108420e+02 */
596 private static final double ln2HI[] = {
597 Double.longBitsToDouble(0x3fe62e42fee00000L), /* 6.93147180369123816490e-01 */
598 Double.longBitsToDouble(0xbfe62e42fee00000L)}; /* -6.93147180369123816490e-01 */
599 private static final double ln2LO[] = {
600 Double.longBitsToDouble(0x3dea39ef35793c76L), /* 1.90821492927058770002e-10 */
601 Double.longBitsToDouble(0xbdea39ef35793c76L)}; /* -1.90821492927058770002e-10 */
602 private static final double invln2 = Double.longBitsToDouble(0x3ff71547652b82feL); /* 1.44269504088896338700e+00 */
603 private static final double P1 = Double.longBitsToDouble(0x3fc555555555553eL); /* 1.66666666666666019037e-01 */
604 private static final double P2 = Double.longBitsToDouble(0xbf66c16c16bebd93L); /* -2.77777777770155933842e-03 */
605 private static final double P3 = Double.longBitsToDouble(0x3f11566aaf25de2cL); /* 6.61375632143793436117e-05 */
606 private static final double P4 = Double.longBitsToDouble(0xbebbbd41c5d26bf1L); /* -1.65339022054652515390e-06 */
607 private static final double P5 = Double.longBitsToDouble(0x3e66376972bea4d0L); /* 4.13813679705723846039e-08 */
608
609 /* exp(x)
610 * Returns the exponential of x.
611 *
612 * Method
613 * 1. Argument reduction:
614 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
615 * Given x, find r and integer k such that
616 *
617 * x = k*ln2 + r, |r| <= 0.5*ln2.
618 *
619 * Here r will be represented as r = hi-lo for better
620 * accuracy.
621 *
622 * 2. Approximation of exp(r) by a special rational function on
623 * the interval [0,0.34658]:
624 * Write
625 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
626 * We use a special Reme algorithm on [0,0.34658] to generate
627 * a polynomial of degree 5 to approximate R. The maximum error
628 * of this polynomial approximation is bounded by 2**-59. In
629 * other words,
630 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
631 * (where z=r*r, and the values of P1 to P5 are listed below)
632 * and
633 * | 5 | -59
634 * | 2.0+P1*z+...+P5*z - R(z) | <= 2
635 * | |
636 * The computation of exp(r) thus becomes
637 * 2*r
638 * exp(r) = 1 + -------
639 * R - r
640 * r*R1(r)
641 * = 1 + r + ----------- (for better accuracy)
642 * 2 - R1(r)
643 * where
644 * 2 4 10
645 * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
646 *
647 * 3. Scale back to obtain exp(x):
648 * From step 1, we have
649 * exp(x) = 2^k * exp(r)
650 *
651 * Special cases:
652 * exp(INF) is INF, exp(NaN) is NaN;
653 * exp(-INF) is 0, and
654 * for finite argument, only exp(0)=1 is exact.
655 *
656 * Accuracy:
657 * according to an error analysis, the error is always less than
658 * 1 ulp (unit in the last place).
659 *
660 * Misc. info.
661 * For IEEE double
662 * if x > 7.09782712893383973096e+02 then exp(x) overflow
663 * if x < -7.45133219101941108420e+02 then exp(x) underflow
664 */
665
666 /***
667 * Returns the exponential of its argument.
668 * @param x The argument, a double.
669 * @return Returns e to the power x.
670 */
671 public static double exp(double x)
672 {
673 double y, hi=0, lo=0, c, t;
674 int k=0, xsb;
675 int hx;
676
677 hx = __HI(x); /* high word of x */
678 xsb = (hx>>>31)&1; /* sign bit of x */
679 hx &= 0x7fffffff; /* high word of |x| */
680
681 /* filter out non-finite argument */
682 if (hx >= 0x40862E42) { /* if |x|>=709.78... */
683 if (hx>=0x7ff00000) {
684 if (((hx&0xfffff)|__LO(x)) != 0) {
685 return x+x; /* NaN */
686 } else {
687 return ((xsb==0) ? x : 0.0); /* exp(+-inf)={inf,0} */
688 }
689 }
690 if (x > o_threshold)
691 return huge*huge; /* overflow */
692 if (x < u_threshold)
693 return twom1000*twom1000; /* underflow */
694 }
695
696 /* argument reduction */
697 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
698 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
699 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
700 } else {
701 k = (int)(invln2*x+halF[xsb]);
702 t = k;
703 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
704 lo = t*ln2LO[0];
705 }
706 x = hi - lo;
707 } else if(hx < 0x3e300000) { /* when |x|<2**-28 */
708 if (huge+x > one)
709 return one+x;/* trigger inexact */
710 } else {
711 k = 0;
712 }
713
714 /* x is now in primary range */
715 t = x*x;
716 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
717 if (k==0) {
718 return one-((x*c)/(c-2.0)-x);
719 } else {
720 y = one-((lo-(x*c)/(2.0-c))-hi);
721 }
722
723 long iy = Double.doubleToLongBits(y);
724 if (k >= -1021) {
725 iy += ((long)k<<52);
726 } else {
727 iy += ((k+1000L)<<52);
728 }
729 return Double.longBitsToDouble(iy);
730 }
731
732 private static final double ln2_hi = Double.longBitsToDouble(0x3fe62e42fee00000L); /* 6.93147180369123816490e-01 */
733 private static final double ln2_lo = Double.longBitsToDouble(0x3dea39ef35793c76L); /* 1.90821492927058770002e-10 */
734 private static final double Lg1 = Double.longBitsToDouble(0x3fe5555555555593L); /* 6.666666666666735130e-01 */
735 private static final double Lg2 = Double.longBitsToDouble(0x3fd999999997fa04L); /* 3.999999999940941908e-01 */
736 private static final double Lg3 = Double.longBitsToDouble(0x3fd2492494229359L); /* 2.857142874366239149e-01 */
737 private static final double Lg4 = Double.longBitsToDouble(0x3fcc71c51d8e78afL); /* 2.222219843214978396e-01 */
738 private static final double Lg5 = Double.longBitsToDouble(0x3fc7466496cb03deL); /* 1.818357216161805012e-01 */
739 private static final double Lg6 = Double.longBitsToDouble(0x3fc39a09d078c69fL); /* 1.531383769920937332e-01 */
740 private static final double Lg7 = Double.longBitsToDouble(0x3fc2f112df3e5244L); /* 1.479819860511658591e-01 */
741
742 /*
743 * Return the logrithm of x
744 *
745 * Method :
746 * 1. Argument Reduction: find k and f such that
747 * x = 2^k * (1+f),
748 * where sqrt(2)/2 < 1+f < sqrt(2) .
749 *
750 * 2. Approximation of log(1+f).
751 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
752 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
753 * = 2s + s*R
754 * We use a special Reme algorithm on [0,0.1716] to generate
755 * a polynomial of degree 14 to approximate R The maximum error
756 * of this polynomial approximation is bounded by 2**-58.45. In
757 * other words,
758 * 2 4 6 8 10 12 14
759 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
760 * (the values of Lg1 to Lg7 are listed in the program)
761 * and
762 * | 2 14 | -58.45
763 * | Lg1*s +...+Lg7*s - R(z) | <= 2
764 * | |
765 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
766 * In order to guarantee error in log below 1ulp, we compute log
767 * by
768 * log(1+f) = f - s*(f - R) (if f is not too large)
769 * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
770 *
771 * 3. Finally, log(x) = k*ln2 + log(1+f).
772 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
773 * Here ln2 is split into two floating point number:
774 * ln2_hi + ln2_lo,
775 * where n*ln2_hi is always exact for |n| < 2000.
776 *
777 * Special cases:
778 * log(x) is NaN with signal if x < 0 (including -INF) ;
779 * log(+INF) is +INF; log(0) is -INF with signal;
780 * log(NaN) is that NaN with no signal.
781 *
782 * Accuracy:
783 * according to an error analysis, the error is always less than
784 * 1 ulp (unit in the last place).
785 */
786
787 /***
788 * Returns the natural logarithm of its argument.
789 * @param x The argument, a double.
790 * @return Returns the natural (base e) logarithm of x.
791 */
792 public static double log(double x)
793 {
794 double hfsq,f,s,z,R,w,t1,t2,dk;
795 int k,hx,i,j;
796 int lx;
797
798 hx = __HI(x); /* high word of x */
799 lx = __LO(x); /* low word of x */
800
801 k=0;
802 if (hx < 0x00100000) { /* x < 2**-1022 */
803 if (((hx&0x7fffffff)|lx)==0)
804 return -two54/zero; /* log(+-0)=-inf */
805 if (hx<0)
806 return (x-x)/zero; /* log(-#) = NaN */
807 k -= 54;
808 x *= two54; /* subnormal number, scale up x */
809 hx = __HI(x); /* high word of x */
810 }
811 if (hx >= 0x7ff00000)
812 return x+x;
813 k += (hx>>20)-1023;
814 hx &= 0x000fffff;
815 i = (hx+0x95f64)&0x100000;
816 x = setHI(x, hx|(i^0x3ff00000)); /* normalize x or x/2 */
817 k += (i>>20);
818 f = x-1.0;
819 if ((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
820 if (f == zero) {
821 if (k == 0) {
822 return zero;
823 } else {
824 dk=(double)k;
825 }
826 return dk*ln2_hi+dk*ln2_lo;
827 }
828 R = f*f*(0.5-0.33333333333333333*f);
829 if (k==0) {
830 return f-R;
831 } else {
832 dk = (double)k;
833 return dk*ln2_hi-((R-dk*ln2_lo)-f);
834 }
835 }
836 s = f/(2.0+f);
837 dk = (double)k;
838 z = s*s;
839 i = hx-0x6147a;
840 w = z*z;
841 j = 0x6b851-hx;
842 t1= w*(Lg2+w*(Lg4+w*Lg6));
843 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
844 i |= j;
845 R = t2+t1;
846 if (i > 0) {
847 hfsq = 0.5*f*f;
848 if (k == 0) {
849 return f-(hfsq-s*(hfsq+R));
850 } else {
851 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
852 }
853 } else {
854 if (k==0) {
855 return f-s*(f-R);
856 } else {
857 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
858 }
859 }
860 }
861
862
863 /***
864 * Returns the sine of its argument.
865 * @param x The argument, a double, assumed to be in radians.
866 * @return Returns the sine of x.
867 */
868 public static double sin(double x)
869 {
870 double y[] = new double[2];
871 double z=0.0;
872 int n;
873 int ix = __HI(x);
874
875 ix &= 0x7fffffff; /* |x| ~< pi/4 */
876
877 if (ix <= 0x3fe921fb) {
878 return __kernel_sin(x,z,0);
879 } else if (ix >= 0x7ff00000) {
880 /* sin(Inf or NaN) is NaN */
881 return x-x;
882 } else {
883 /* argument reduction needed */
884 n = __ieee754_rem_pio2(x, y);
885 switch(n&3) {
886 case 0:
887 return __kernel_sin(y[0],y[1], 1);
888 case 1:
889 return __kernel_cos(y[0],y[1]);
890 case 2:
891 return -__kernel_sin(y[0],y[1], 1);
892 default:
893 return -__kernel_cos(y[0],y[1]);
894 }
895 }
896 }
897
898
899 private static double S1 = -1.66666666666666324348e-01; /* 0xBFC55555, 0x55555549 */
900 private static double S2 = 8.33333333332248946124e-03; /* 0x3F811111, 0x1110F8A6 */
901 private static double S3 = -1.98412698298579493134e-04; /* 0xBF2A01A0, 0x19C161D5 */
902 private static double S4 = 2.75573137070700676789e-06; /* 0x3EC71DE3, 0x57B1FE7D */
903 private static double S5 = -2.50507602534068634195e-08; /* 0xBE5AE5E6, 0x8A2B9CEB */
904 private static double S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
905
906 /*
907 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
908 * Input x is assumed to be bounded by ~pi/4 in magnitude.
909 * Input y is the tail of x. * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
910 *
911 * Algorithm
912 * 1. Since sin(-x) = -sin(x), we need only to consider positive x.
913 * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
914 * 3. sin(x) is approximated by a polynomial of degree 13 on
915 * [0,pi/4]
916 * 3 13
917 * sin(x) ~ x + S1*x + ... + S6*x
918 * where
919 *
920 * |sin(x) 2 4 6 8 10 12 | -58
921 * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
922 * | x |
923 *
924 * 4. sin(x+y) = sin(x) + sin'(x')*y
925 * ~ sin(x) + (1-x*x/2)*y
926 * For better accuracy, let
927 * 3 2 2 2 2
928 * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
929 * then 3 2
930 * sin(x) = x + (S1*x + (x *(r-y/2)+y))
931 */
932 static double __kernel_sin(double x, double y, int iy)
933 {
934 double z,r,v;
935 int ix;
936 ix = __HI(x)&0x7fffffff; /* high word of x */
937 if (ix<0x3e400000) { /* |x| < 2**-27 */
938 if((int)x ==0 )
939 return x; /* generate inexact */
940 }
941 z = x*x;
942 v = z*x;
943 r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
944 if (iy == 0)
945 return x+v*(S1+z*r);
946 else
947 return x-((z*(half*y-v*r)-y)-v*S1);
948 }
949
950 /***
951 * Returns the cosine of its argument.
952 * @param x The argument, a double, assumed to be in radians.
953 * @return Returns the cosine of x.
954 */
955 public static double cos(double x)
956 {
957 double z = 0.0;
958 double y[] = new double[2];
959 int n, ix;
960
961 /* High word of x. */
962 ix = __HI(x);
963
964 /* |x| ~< pi/4 */
965 ix &= 0x7fffffff;
966 if (ix <= 0x3fe921fb) {
967 return __kernel_cos(x,z);
968
969 /* cos(Inf or NaN) is NaN */
970 } else if (ix >= 0x7ff00000) {
971 return x-x;
972
973 /* argument reduction needed */
974 } else {
975 n = __ieee754_rem_pio2(x,y);
976 switch(n&3) {
977 case 0: return __kernel_cos(y[0],y[1]);
978 case 1: return -__kernel_sin(y[0],y[1],1);
979 case 2: return -__kernel_cos(y[0],y[1]);
980 default:
981 return __kernel_sin(y[0],y[1],1);
982 }
983 }
984 }
985
986 private static final double one = Double.longBitsToDouble(0x3ff0000000000000L); /* 1.00000000000000000000e+00 */
987 private static final double C1 = Double.longBitsToDouble(0x3fa555555555554cL); /* 4.16666666666666019037e-02 */
988 private static final double C2 = Double.longBitsToDouble(0xbf56c16c16c15177L); /* -1.38888888888741095749e-03 */
989 private static final double C3 = Double.longBitsToDouble(0x3efa01a019cb1590L); /* 2.48015872894767294178e-05 */
990 private static final double C4 = Double.longBitsToDouble(0xbe927e4f809c52adL); /* -2.75573143513906633035e-07 */
991 private static final double C5 = Double.longBitsToDouble(0x3e21ee9ebdb4b1c4L); /* 2.08757232129817482790e-09 */
992 private static final double C6 = Double.longBitsToDouble(0xbda8fae9be8838d4L); /* -1.13596475577881948265e-11 */
993
994 /*
995 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
996 * Input x is assumed to be bounded by ~pi/4 in magnitude.
997 * Input y is the tail of x.
998 *
999 * Algorithm
1000 * 1. Since cos(-x) = cos(x), we need only to consider positive x.
1001 * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
1002 * 3. cos(x) is approximated by a polynomial of degree 14 on
1003 * [0,pi/4]
1004 * 4 14
1005 * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
1006 * where the remez error is
1007 *
1008 * | 2 4 6 8 10 12 14 | -58
1009 * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
1010 * | |
1011 *
1012 * 4 6 8 10 12 14
1013 * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
1014 * cos(x) = 1 - x*x/2 + r
1015 * since cos(x+y) ~ cos(x) - sin(x)*y
1016 * ~ cos(x) - x*y,
1017 * a correction term is necessary in cos(x) and hence
1018 * cos(x+y) = 1 - (x*x/2 - (r - x*y))
1019 * For better accuracy when x > 0.3, let qx = |x|/4 with
1020 * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
1021 * Then
1022 * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
1023 * Note that 1-qx and (x*x/2-qx) is EXACT here, and the
1024 * magnitude of the latter is at least a quarter of x*x/2,
1025 * thus, reducing the rounding error in the subtraction.
1026 */
1027 private static double __kernel_cos(double x, double y)
1028 {
1029 double a, hz, z, r, qx=zero;
1030 int ix;
1031 ix = __HI(x)&0x7fffffff; /* ix = |x|'s high word*/
1032 if (ix < 0x3e400000) {
1033 /* if x < 2**27 */
1034 if (((int)x) == 0) return one; /* generate inexact */
1035 }
1036 z = x*x;
1037 r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
1038 if (ix < 0x3FD33333) {
1039 /* if |x| < 0.3 */
1040 return one - (0.5*z - (z*r - x*y));
1041 } else {
1042 if(ix > 0x3fe90000) { /* x > 0.78125 */
1043 qx = 0.28125;
1044 } else {
1045 qx = set(ix-0x00200000, 0);/* x/4 */
1046 }
1047 hz = 0.5*z-qx;
1048 a = one-qx;
1049 return a - (hz - (z*r-x*y));
1050 }
1051 }
1052
1053 /***
1054 * Returns the tangent of its argument.
1055 * @param x The argument, a double, assumed to be in radians.
1056 * @return Returns the tangent of x.
1057 */
1058 public static double tan(double x)
1059 {
1060 double z = zero;
1061
1062 int n;
1063 int ix = __HI(x);
1064
1065 /* |x| ~< pi/4 */
1066 ix &= 0x7fffffff;
1067 if (ix <= 0x3fe921fb) {
1068 return __kernel_tan(x, z, 1);
1069 } else if (ix >= 0x7ff00000) {
1070 /* tan(Inf or NaN) is NaN */
1071 return x-x; /* NaN */
1072 } else {
1073 /* argument reduction needed */
1074 double y[] = new double[2];
1075 n = __ieee754_rem_pio2(x, y);
1076 /* 1 -- n even -1 -- n odd */
1077 return __kernel_tan(y[0], y[1], 1-((n&1)<<1));
1078 }
1079 }
1080
1081
1082 private static final double pio4 = Double.longBitsToDouble(0x3fe921fb54442d18L); /* 7.85398163397448278999e-01 */
1083 private static final double pio4lo= Double.longBitsToDouble(0x3c81a62633145c07L); /* 3.06161699786838301793e-17 */
1084 private static final double T[] = {
1085 Double.longBitsToDouble(0x3fd5555555555563L), /* 3.33333333333334091986e-01 */
1086 Double.longBitsToDouble(0x3fc111111110fe7aL), /* 1.33333333333201242699e-01 */
1087 Double.longBitsToDouble(0x3faba1ba1bb341feL), /* 5.39682539762260521377e-02 */
1088 Double.longBitsToDouble(0x3f9664f48406d637L), /* 2.18694882948595424599e-02 */
1089 Double.longBitsToDouble(0x3f8226e3e96e8493L), /* 8.86323982359930005737e-03 */
1090 Double.longBitsToDouble(0x3f6d6d22c9560328L), /* 3.59207910759131235356e-03 */
1091 Double.longBitsToDouble(0x3f57dbc8fee08315L), /* 1.45620945432529025516e-03 */
1092 Double.longBitsToDouble(0x3f4344d8f2f26501L), /* 5.88041240820264096874e-04 */
1093 Double.longBitsToDouble(0x3f3026f71a8d1068L), /* 2.46463134818469906812e-04 */
1094 Double.longBitsToDouble(0x3f147e88a03792a6L), /* 7.81794442939557092300e-05 */
1095 Double.longBitsToDouble(0x3f12b80f32f0a7e9L), /* 7.14072491382608190305e-05 */
1096 Double.longBitsToDouble(0xbef375cbdb605373L), /* -1.85586374855275456654e-05 */
1097 Double.longBitsToDouble(0x3efb2a7074bf7ad4L) /* 2.59073051863633712884e-05 */
1098 };
1099
1100 /*
1101 * __kernel_tan( x, y, k )
1102 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
1103 * Input x is assumed to be bounded by ~pi/4 in magnitude.
1104 * Input y is the tail of x. * Input k indicates whether tan (if k=1) or
1105 * -1/tan (if k= -1) is returned. * * Algorithm
1106 * 1. Since tan(-x) = -tan(x), we need only to consider positive x.
1107 * 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
1108 * 3. tan(x) is approximated by a odd polynomial of degree 27 on
1109 * [0,0.67434]
1110 * 3 27
1111 * tan(x) ~ x + T1*x + ... + T13*x
1112 * where
1113 *
1114 * |tan(x) 2 4 26 | -59.2
1115 * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
1116 * | x |
1117 *
1118 * Note: tan(x+y) = tan(x) + tan'(x)*y
1119 * ~ tan(x) + (1+x*x)*y
1120 * Therefore, for better accuracy in computing tan(x+y), let
1121 * 3 2 2 2 2
1122 * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
1123 * then
1124 * 3 2
1125 * tan(x+y) = x + (T1*x + (x *(r+y)+y)) *
1126 * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
1127 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
1128 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
1129 */
1130 private static double __kernel_tan(double x, double y, int iy)
1131 {
1132 double z, r, v, w, s;
1133 int ix,hx;
1134
1135 hx = __HI(x); /* high word of x */
1136 ix = hx&0x7fffffff; /* high word of |x| */
1137 if (ix<0x3e300000) { /* x < 2**-28 */
1138 if ((int)x==0) { /* generate inexact */
1139 if(((ix|__LO(x))|(iy+1))==0)
1140 return one/abs(x);
1141 else return (iy==1)? x: -one/x;
1142 }
1143 }
1144 if (ix>=0x3FE59428) {
1145 /* |x|>=0.6744 */
1146 if (hx<0) {
1147 x = -x;
1148 y = -y;
1149 }
1150 z = pio4-x;
1151 w = pio4lo-y;
1152 x = z+w;
1153 y = 0.0;
1154 }
1155 z = x*x;
1156 w = z*z;
1157
1158 /*
1159 * Break x^5*(T[1]+x^2*T[2]+...) into
1160 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
1161 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
1162 */
1163 r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
1164 v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
1165 s = z*x;
1166 r = y + z*(s*(r+v)+y);
1167 r += T[0]*s;
1168 w = x+r;
1169 if (ix >= 0x3FE59428) {
1170 v = (double)iy;
1171 return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
1172 }
1173 if (iy==1) {
1174 return w;
1175 } else {
1176 /* if allow error up to 2 ulp, simply return -1.0/(x+r) here */
1177 /* compute -1.0/(x+r) accurately */
1178 double a,t;
1179 z = w;
1180 z = setLO(z, 0);
1181 v = r-(z - x);
1182 /* z+v = r+x */
1183 t = a = -1.0/w;
1184 /* a = -1.0/w */
1185 t = setLO(t, 0);
1186 s = 1.0+t*z;
1187 return t+a*(s+t*v);
1188 }
1189 }
1190
1191
1192 private static final double pio2_hi = Double.longBitsToDouble(0x3FF921FB54442D18L); /* 1.57079632679489655800e+00 */
1193 private static final double pio2_lo = Double.longBitsToDouble(0x3C91A62633145C07L); /* 6.12323399573676603587e-17 */
1194 private static final double pio4_hi = Double.longBitsToDouble(0x3FE921FB54442D18L); /* 7.85398163397448278999e-01 */
1195 /* coefficient for R(x^2) */
1196 private static final double pS0 = Double.longBitsToDouble(0x3fc5555555555555L); /* 1.66666666666666657415e-01 */
1197 private static final double pS1 = Double.longBitsToDouble(0xbfd4d61203eb6f7dL); /* -3.25565818622400915405e-01 */
1198 private static final double pS2 = Double.longBitsToDouble(0x3fc9c1550e884455L); /* 2.01212532134862925881e-01 */
1199 private static final double pS3 = Double.longBitsToDouble(0xbfa48228b5688f3bL); /* -4.00555345006794114027e-02 */
1200 private static final double pS4 = Double.longBitsToDouble(0x3f49efe07501b288L); /* 7.91534994289814532176e-04 */
1201 private static final double pS5 = Double.longBitsToDouble(0x3f023de10dfdf709L); /* 3.47933107596021167570e-05 */
1202 private static final double qS1 = Double.longBitsToDouble(0xc0033a271c8a2d4bL); /* -2.40339491173441421878e+00 */
1203 private static final double qS2 = Double.longBitsToDouble(0x40002ae59c598ac8L); /* 2.02094576023350569471e+00 */
1204 private static final double qS3 = Double.longBitsToDouble(0xbfe6066c1b8d0159L); /* -6.88283971605453293030e-01 */
1205 private static final double qS4 = Double.longBitsToDouble(0x3fb3b8c5b12e9282L); /* 7.70381505559019352791e-02 */
1206
1207 /*
1208 * asin(x)
1209 * Method :
1210 * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
1211 * we approximate asin(x) on [0,0.5] by
1212 * asin(x) = x + x*x^2*R(x^2)
1213 * where
1214 * R(x^2) is a rational approximation of (asin(x)-x)/x^3
1215 * and its remez error is bounded by
1216 * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
1217 *
1218 * For x in [0.5,1]
1219 * asin(x) = pi/2-2*asin(sqrt((1-x)/2))
1220 * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
1221 * then for x>0.98
1222 * asin(x) = pi/2 - 2*(s+s*z*R(z))
1223 * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
1224 * For x<=0.98, let pio4_hi = pio2_hi/2, then
1225 * f = hi part of s;
1226 * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
1227 * and
1228 * asin(x) = pi/2 - 2*(s+s*z*R(z))
1229 * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
1230 * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
1231 *
1232 * Special cases:
1233 * if x is NaN, return x itself;
1234 * if |x|>1, return NaN with invalid signal.
1235 *
1236 */
1237
1238 /***
1239 * Returns the inverse (arc) sine of its argument.
1240 * @param x The argument, a double.
1241 * @return Returns the angle, in radians, whose sine is x.
1242 * It is in the range [-pi/2,pi/2].
1243 */
1244 public static double asin(double x)
1245 {
1246 double t=zero,w,p,q,c,r,s;
1247 int hx,ix;
1248 hx = __HI(x);
1249
1250 ix = hx&0x7fffffff;
1251 if (ix>= 0x3ff00000) { /* |x|>= 1 */
1252 if(((ix-0x3ff00000)|__LO(x))==0)
1253 /* asin(1)=+-pi/2 with inexact */
1254 return x*pio2_hi+x*pio2_lo;
1255 return (x-x)/(x-x); /* asin(|x|>1) is NaN */
1256 } else if (ix<0x3fe00000) { /* |x|<0.5 */
1257 if (ix<0x3e400000) { /* if |x| < 2**-27 */
1258 if(huge+x>one) return x;/* return x with inexact if x!=0*/
1259 } else
1260 t = x*x;
1261 p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
1262 q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
1263 w = p/q;
1264 return x+x*w;
1265 }
1266 /* 1> |x|>= 0.5 */
1267 w = one - abs(x);
1268 t = w*0.5;
1269 p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
1270 q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
1271 s = sqrt(t);
1272 if(ix>=0x3FEF3333) { /* if |x| > 0.975 */
1273 w = p/q;
1274 t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
1275 } else {
1276 w = s;
1277 w = setLO(w, 0);
1278 c = (t-w*w)/(s+w);
1279 r = p/q;
1280 p = 2.0*s*r-(pio2_lo-2.0*c);
1281 q = pio4_hi-2.0*w;
1282 t = pio4_hi-(p-q);
1283 }
1284 return ((hx>0) ? t : -t);
1285 }
1286
1287 /*
1288 * Method :
1289 * acos(x) = pi/2 - asin(x)
1290 * acos(-x) = pi/2 + asin(x)
1291 * For |x|<=0.5
1292 * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
1293 * For x>0.5
1294 * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
1295 * = 2asin(sqrt((1-x)/2))
1296 * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
1297 * = 2f + (2c + 2s*z*R(z))
1298 * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
1299 * for f so that f+c ~ sqrt(z).
1300 * For x<-0.5
1301 * acos(x) = pi - 2asin(sqrt((1-|x|)/2))
1302 * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
1303 *
1304 * Special cases:
1305 * if x is NaN, return x itself;
1306 * if |x|>1, return NaN with invalid signal.
1307 *
1308 * Function needed: sqrt
1309 */
1310
1311 /***
1312 * Returns the inverse (arc) cosine of its argument.
1313 * @param x The argument, a double.
1314 * @return Returns the angle, in radians, whose cosine is x.
1315 * It is in the range [0,pi].
1316 */
1317 public static double acos(double x)
1318 {
1319 double z,p,q,r,w,s,c,df;
1320 int hx,ix;
1321 hx = __HI(x);
1322 ix = hx&0x7fffffff;
1323 if (ix >= 0x3ff00000) { /* |x| >= 1 */
1324 if (((ix-0x3ff00000)|__LO(x)) == 0) { /* |x|==1 */
1325 if (hx>0)
1326 return 0.0; /* acos(1) = 0 */
1327 else
1328 return PI+2.0*pio2_lo; /* acos(-1)= pi */
1329 }
1330 return (x-x)/(x-x); /* acos(|x|>1) is NaN */
1331 }
1332 if (ix < 0x3fe00000) { /* |x| < 0.5 */
1333 if (ix <= 0x3c600000)
1334 return pio2_hi+pio2_lo;/*if|x|<2**-57*/
1335 z = x*x;
1336 p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
1337 q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
1338 r = p/q;
1339 return pio2_hi - (x - (pio2_lo-x*r));
1340 } else if (hx<0) { /* x < -0.5 */
1341 z = (one+x)*0.5;
1342 p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
1343 q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
1344 s = sqrt(z);
1345 r = p/q;
1346 w = r*s-pio2_lo;
1347 return PI - 2.0*(s+w);
1348 } else { /* x > 0.5 */
1349 z = (one-x)*0.5;
1350 s = sqrt(z);
1351 df = s;
1352 df = setLO(df, 0);
1353 c = (z-df*df)/(s+df);
1354 p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
1355 q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
1356 r = p/q;
1357 w = r*s+c;
1358 return 2.0*(df+w);
1359 }
1360 }
1361
1362
1363 private static final double atanhi[] = {
1364 Double.longBitsToDouble(0x3fddac670561bb4fL), /* 4.63647609000806093515e-01 atan(0.5)hi */
1365 Double.longBitsToDouble(0x3fe921fb54442d18L), /* 7.85398163397448278999e-01 atan(1.0)hi */
1366 Double.longBitsToDouble(0x3fef730bd281f69bL), /* 9.82793723247329054082e-01 atan(1.5)hi */
1367 Double.longBitsToDouble(0x3ff921fb54442d18L) /* 1.57079632679489655800e+00 atan(inf)hi */
1368 };
1369
1370 private static final double atanlo[] = {
1371 Double.longBitsToDouble(0x3c7a2b7f222f65e2L), /* 2.26987774529616870924e-17 atan(0.5)lo */
1372 Double.longBitsToDouble(0x3c81a62633145c07L), /* 3.06161699786838301793e-17 atan(1.0)lo */
1373 Double.longBitsToDouble(0x3c7007887af0cbbdL), /* 1.39033110312309984516e-17 atan(1.5)lo */
1374 Double.longBitsToDouble(0x3c91a62633145c07L) /* 6.12323399573676603587e-17 atan(inf)lo */
1375 };
1376
1377 private static final double aT[] = {
1378 Double.longBitsToDouble(0x3fd555555555550dL), /* 3.33333333333329318027e-01 */
1379 Double.longBitsToDouble(0xbfc999999998ebc4L), /* -1.99999999998764832476e-01 */
1380 Double.longBitsToDouble(0x3fc24924920083ffL), /* 1.42857142725034663711e-01 */
1381 Double.longBitsToDouble(0xbfbc71c6fe231671L), /* -1.11111104054623557880e-01 */
1382 Double.longBitsToDouble(0x3fb745cdc54c206eL), /* 9.09088713343650656196e-02 */
1383 Double.longBitsToDouble(0xbfb3b0f2af749a6dL), /* -7.69187620504482999495e-02 */
1384 Double.longBitsToDouble(0x3fb10d66a0d03d51L), /* 6.66107313738753120669e-02 */
1385 Double.longBitsToDouble(0xbfadde2d52defd9aL), /* -5.83357013379057348645e-02 */
1386 Double.longBitsToDouble(0x3fa97b4b24760debL), /* 4.97687799461593236017e-02 */
1387 Double.longBitsToDouble(0xbfa2b4442c6a6c2fL), /* -3.65315727442169155270e-02 */
1388 Double.longBitsToDouble(0x3f90ad3ae322da11L) /* 1.62858201153657823623e-02 */
1389 };
1390
1391 /*
1392 * Method
1393 * 1. Reduce x to positive by atan(x) = -atan(-x).
1394 * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
1395 * is further reduced to one of the following intervals and the
1396 * arctangent of t is evaluated by the corresponding formula:
1397 *
1398 * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
1399 * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
1400 * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
1401 * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
1402 * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
1403 *
1404 * Constants:
1405 * The hexadecimal values are the intended ones for the following
1406 * constants. The decimal values may be used, provided that the
1407 * compiler will convert from decimal to binary accurately enough
1408 * to produce the hexadecimal values shown.
1409 */
1410
1411 /***
1412 * Returns the inverse (arc) tangent of its argument.
1413 * @param x The argument, a double.
1414 * @return Returns the angle, in radians, whose tangent is x.
1415 * It is in the range [-pi/2,pi/2].
1416 */
1417 public static double atan(double x)
1418 {
1419 double w,s1,s2,z;
1420 int ix,hx,id;
1421
1422 hx = __HI(x);
1423 ix = hx&0x7fffffff;
1424 if(ix>=0x44100000) { /* if |x| >= 2^66 */
1425 if(ix>0x7ff00000||
1426 (ix==0x7ff00000&&(__LO(x)!=0)))
1427 return x+x; /* NaN */
1428 if(hx>0) return atanhi[3]+atanlo[3];
1429 else return -atanhi[3]-atanlo[3];
1430 } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
1431 if (ix < 0x3e200000) { /* |x| < 2^-29 */
1432 if(huge+x>one) return x; /* raise inexact */
1433 }
1434 id = -1;
1435 } else {
1436 x = abs(x);
1437 if (ix < 0x3ff30000) { /* |x| < 1.1875 */
1438 if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
1439 id = 0; x = (2.0*x-one)/(2.0+x);
1440 } else { /* 11/16<=|x|< 19/16 */
1441 id = 1; x = (x-one)/(x+one);
1442 }
1443 } else {
1444 if (ix < 0x40038000) { /* |x| < 2.4375 */
1445 id = 2; x = (x-1.5)/(one+1.5*x);
1446 } else { /* 2.4375 <= |x| < 2^66 */
1447 id = 3; x = -1.0/x;
1448 }
1449 }
1450 }
1451 /* end of argument reduction */
1452 z = x*x;
1453 w = z*z;
1454 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
1455 s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
1456 s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
1457 if (id<0) return x - x*(s1+s2);
1458 else {
1459 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
1460 return (hx<0)? -z:z;
1461 }
1462 }
1463
1464
1465 private static final double pi_o_4 = Double.longBitsToDouble(0x3fe921fb54442d18L); /* 7.8539816339744827900e-01 */
1466 private static final double pi_o_2 = Double.longBitsToDouble(0x3ff921fb54442d18L); /* 1.5707963267948965580e+00 */
1467 private static final double pi_lo = Double.longBitsToDouble(0x3ca1a62633145c07L); /* 1.2246467991473531772e-16 */
1468
1469 /*
1470 * Method :
1471 * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
1472 * 2. Reduce x to positive by (if x and y are unexceptional):
1473 * ARG (x+iy) = arctan(y/x) ... if x > 0,
1474 * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
1475 *
1476 * Special cases:
1477 *
1478 * ATAN2((anything), NaN ) is NaN;
1479 * ATAN2(NAN , (anything) ) is NaN;
1480 * ATAN2(+-0, +(anything but NaN)) is +-0 ;
1481 * ATAN2(+-0, -(anything but NaN)) is +-pi ;
1482 * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
1483 * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
1484 * ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
1485 * ATAN2(+-INF,+INF ) is +-pi/4 ;
1486 * ATAN2(+-INF,-INF ) is +-3pi/4;
1487 * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
1488 *
1489 * Constants:
1490 * The hexadecimal values are the intended ones for the following
1491 * constants. The decimal values may be used, provided that the
1492 * compiler will convert from decimal to binary accurately enough
1493 * to produce the hexadecimal values shown.
1494 */
1495
1496 /***
1497 * Returns angle corresponding to a Cartesian point.
1498 * @param x The first argument, a double.
1499 * @param y The second argument, a double.
1500 * @return Returns the angle, in radians, the the line
1501 * from (0,0) to (x,y) makes with the x-axis.
1502 * It is in the range [-pi,pi].
1503 */
1504 public static double atan2(double y, double x)
1505 {
1506 double z;
1507 int k,m,hx,hy,ix,iy;
1508 int lx,ly;
1509
1510 hx = __HI(x);
1511 ix = hx&0x7fffffff;
1512 lx = __LO(x);
1513 hy = __HI(y);
1514 iy = hy&0x7fffffff;
1515 ly = __LO(y);
1516 if (((ix|((lx|-lx)>>>31))>0x7ff00000)||
1517 ((iy|((ly|-ly)>>>31))>0x7ff00000)) /* x or y is NaN */
1518 return x+y;
1519 if ((hx-0x3ff00000|lx) == 0)
1520 return atan(y); /* x=1.0 */
1521 m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
1522
1523 /* when y = 0 */
1524 if((iy|ly)==0) {
1525 switch(m) {
1526 case 0:
1527 case 1: return y; /* atan(+-0,+anything)=+-0 */
1528 case 2: return PI+tiny;/* atan(+0,-anything) = pi */
1529 case 3: return -PI-tiny;/* atan(-0,-anything) =-pi */
1530 }
1531 }
1532 /* when x = 0 */
1533 if((ix|lx) == 0)
1534 return ((hy<0) ? -pi_o_2-tiny: pi_o_2+tiny);
1535
1536 /* when x is INF */
1537 if (ix == 0x7ff00000) {
1538 if (iy == 0x7ff00000) {
1539 switch(m) {
1540 case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
1541 case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
1542 case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
1543 case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
1544 }
1545 } else {
1546 switch(m) {
1547 case 0: return zero ; /* atan(+...,+INF) */
1548 case 1: return -zero ; /* atan(-...,+INF) */
1549 case 2: return PI+tiny ; /* atan(+...,-INF) */
1550 case 3: return -PI-tiny ; /* atan(-...,-INF) */
1551 }
1552 }
1553 }
1554 /* when y is INF */
1555 if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
1556
1557 /* compute y/x */
1558 k = (iy-ix)>>20;
1559 if (k > 60) {
1560 /* |y/x| > 2**60 */
1561 z=pi_o_2+0.5*pi_lo;
1562 } else if (hx<0&&k<-60) {
1563 /* |y|/x < -2**60 */
1564 z=0.0;
1565 } else {
1566 /* safe to do y/x */
1567 z=atan(abs(y/x));
1568 }
1569
1570 switch (m) {
1571 case 0:
1572 return z ; /* atan(+,+) */
1573 case 1:
1574 return setHI(z, __HI(z) ^ 0x80000000); /* atan(-,+) */
1575 case 2:
1576 return PI-(z-pi_lo);/* atan(+,-) */
1577 default: /* case 3 */
1578 return (z-pi_lo)-PI;/* atan(-,-) */
1579 }
1580 }
1581
1582
1583
1584 /*
1585 * kernel function:
1586 * __kernel_sin ... sine function on [-pi/4,pi/4]
1587 * __ieee754_rem_pio2 ... argument reduction routine
1588 *
1589 * Method.
1590 * Let S,C and T denote the sin, cos and tan respectively on
1591 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
1592 * in [-pi/4 , +pi/4], and let n = k mod 4.
1593 * We have
1594 *
1595 * n sin(x) cos(x) tan(x)
1596 * ----------------------------------------------------------
1597 * 0 S C T
1598 * 1 C -S -1/T
1599 * 2 -S -C T
1600 * 3 -C S -1/T
1601 * ----------------------------------------------------------
1602 *
1603 * Special cases:
1604 * Let trig be any of sin, cos, or tan.
1605 * trig(+-INF) is NaN, with signals;
1606 * trig(NaN) is that NaN;
1607 *
1608 * Accuracy:
1609 * TRIG(x) returns trig(x) nearly rounded
1610 */
1611
1612
1613
1614 /*
1615 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
1616 */
1617 private static final int two_over_pi[] = {
1618 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
1619 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
1620 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
1621 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
1622 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
1623 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
1624 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
1625 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
1626 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
1627 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
1628 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b};
1629
1630 private static final int npio2_hw[] = {
1631 0x3ff921fb, 0x400921fb, 0x4012d97c, 0x401921fb, 0x401f6a7a, 0x4022d97c,
1632 0x4025fdbb, 0x402921fb, 0x402c463a, 0x402f6a7a, 0x4031475c, 0x4032d97c,
1633 0x40346b9c, 0x4035fdbb, 0x40378fdb, 0x403921fb, 0x403ab41b, 0x403c463a,
1634 0x403dd85a, 0x403f6a7a, 0x40407e4c, 0x4041475c, 0x4042106c, 0x4042d97c,
1635 0x4043a28c, 0x40446b9c, 0x404534ac, 0x4045fdbb, 0x4046c6cb, 0x40478fdb,
1636 0x404858eb, 0x404921fb};
1637
1638 private static final double zero = 0.00000000000000000000e+00; // 0x0000000000000000
1639 private static final double half = Double.longBitsToDouble(0x3fe0000000000000L); /* 5.00000000000000000000e-01 */
1640 private static final double two24 = Double.longBitsToDouble(0x4170000000000000L); /* 1.67772160000000000000e+07 */
1641 private static final double invpio2 = Double.longBitsToDouble(0x3fe45f306dc9c883L); /* 6.36619772367581382433e-01 53 bits of 2/pi */
1642 private static final double pio2_1 = Double.longBitsToDouble(0x3ff921fb54400000L); /* 1.57079632673412561417e+00 first 33 bit of pi/2 */
1643 private static final double pio2_1t = Double.longBitsToDouble(0x3dd0b4611a626331L); /* 6.07710050650619224932e-11 pi/2 - pio2_1 */
1644 private static final double pio2_2 = Double.longBitsToDouble(0x3dd0b4611a600000L); /* 6.07710050630396597660e-11 second 33 bit of pi/2 */
1645 private static final double pio2_2t = Double.longBitsToDouble(0x3ba3198a2e037073L); /* 2.02226624879595063154e-21 pi/2 - (pio2_1+pio2_2) */
1646 private static final double pio2_3 = Double.longBitsToDouble(0x3ba3198a2e000000L); /* 2.02226624871116645580e-21 third 33 bit of pi/2 */
1647 private static final double pio2_3t = Double.longBitsToDouble(0x397b839a252049c1L); /* 8.47842766036889956997e-32 pi/2 - (pio2_1+pio2_2+pio2_3) */
1648
1649
1650 /*
1651 * Return the remainder of x % pi/2 in y[0]+y[1]
1652 */
1653 private static int __ieee754_rem_pio2(double x, double y[])
1654 {
1655 double z = zero, w, t, r, fn;
1656 double tx[] = new double[3];
1657 int i, j, nx, n, ix, hx;
1658
1659 hx = __HI(x); /* high word of x */
1660 ix = hx&0x7fffffff;
1661 if (ix <= 0x3fe921fb) {
1662 /* |x| ~<= pi/4 , no need for reduction */
1663 y[0] = x;
1664 y[1] = 0;
1665 return 0;
1666 }
1667
1668 if (ix < 0x4002d97c) {
1669 /* |x| < 3pi/4, special case with n=+-1 */
1670 if (hx > 0) {
1671 z = x - pio2_1;
1672 if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
1673 y[0] = z - pio2_1t;
1674 y[1] = (z-y[0])-pio2_1t;
1675 } else { /* near pi/2, use 33+33+53 bit pi */
1676 z -= pio2_2;
1677 y[0] = z - pio2_2t;
1678 y[1] = (z-y[0])-pio2_2t;
1679 }
1680 return 1;
1681 } else { /* negative x */
1682 z = x + pio2_1;
1683 if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
1684 y[0] = z + pio2_1t;
1685 y[1] = (z-y[0])+pio2_1t;
1686 } else { /* near pi/2, use 33+33+53 bit pi */
1687 z += pio2_2;
1688 y[0] = z + pio2_2t;
1689 y[1] = (z-y[0])+pio2_2t;
1690 }
1691 return -1;
1692 }
1693 }
1694
1695 if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
1696 t = abs(x);
1697 n = (int) (t*invpio2+half);
1698 fn = (double)n;
1699 r = t-fn*pio2_1;
1700 w = fn*pio2_1t; /* 1st round good to 85 bit */
1701 if (n<32&&ix!=npio2_hw[n-1]) {
1702 y[0] = r-w; /* quick check no cancellation */
1703 } else {
1704 j = ix>>20;
1705 y[0] = r-w;
1706 i = j-(((__HI(y[0]))>>20)&0x7ff);
1707 if (i>16) { /* 2nd iteration needed, good to 118 */
1708 t = r;
1709 w = fn*pio2_2;
1710 r = t-w;
1711 w = fn*pio2_2t-((t-r)-w);
1712 y[0] = r-w;
1713 i = j-(((__HI(y[0]))>>20)&0x7ff);
1714 if (i>49) { /* 3rd iteration need, 151 bits acc */
1715 t = r; /* will cover all possible cases */
1716 w = fn*pio2_3;
1717 r = t-w;
1718 w = fn*pio2_3t-((t-r)-w);
1719 y[0] = r-w;
1720 }
1721 }
1722 }
1723 y[1] = (r-y[0])-w;
1724 if (hx<0) {
1725 y[0] = -y[0];
1726 y[1] = -y[1];
1727 return -n;
1728 } else {
1729 return n;
1730 }
1731 }
1732
1733 /*
1734 * all other (large) arguments
1735 */
1736 if (ix >= 0x7ff00000) {
1737 /* x is inf or NaN */
1738 y[0] = y[1] = x-x;
1739 return 0;
1740 }
1741
1742 /* set z = scalbn(|x|,ilogb(x)-23) */
1743 long lx = Double.doubleToLongBits(x);
1744 long exp = (0x7ff0000000000000L&lx) >> 52;
1745 exp -= 1046;
1746 lx -= (exp<<52);
1747 lx &= 0x7fffffffffffffffL;
1748 z = Double.longBitsToDouble(lx);
1749 for (i = 0; i < 2; i++) {
1750 tx[i] = (double)((int)(z));
1751 z = (z-tx[i])*two24;
1752 }
1753 tx[2] = z;
1754 nx = 3;
1755 while (tx[nx-1] == zero) nx--; /* skip zero term */
1756 n = __kernel_rem_pio2(tx, y, (int)exp, nx);
1757 //System.out.println("KERNEL");
1758 //System.out.println("tx "+tx[0]+" "+tx[1]+" "+tx[2]);
1759 //System.out.println("y "+y[0]+" "+y[1]);
1760 //System.out.println("exp "+(int)exp+" nx "+nx+" n "+n);
1761 if (hx < 0) {
1762 y[0] = -y[0];
1763 y[1] = -y[1];
1764 return -n;
1765 }
1766 return n;
1767 }
1768
1769
1770 /*
1771 * __kernel_rem_pio2(x,y,e0,nx)
1772 * double x[],y[]; int e0,nx; int two_over_pi[];
1773 *
1774 * __kernel_rem_pio2 return the last three digits of N with
1775 * y = x - N*pi/2
1776 * so that |y| < pi/2.
1777 *
1778 * The method is to compute the integer (mod 8) and fraction parts of
1779 * (2/pi)*x without doing the full multiplication. In general we
1780 * skip the part of the product that are known to be a huge integer
1781 * (more accurately, = 0 mod 8 ). Thus the number of operations are
1782 * independent of the exponent of the input.
1783 *
1784 * (2/pi) is represented by an array of 24-bit integers in two_over_pi[].
1785 *
1786 * Input parameters:
1787 * x[] The input value (must be positive) is broken into nx
1788 * pieces of 24-bit integers in double precision format.
1789 * x[i] will be the i-th 24 bit of x. The scaled exponent
1790 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
1791 * match x's up to 24 bits.
1792 *
1793 * Example of breaking a double positive z into x[0]+x[1]+x[2]:
1794 * e0 = ilogb(z)-23
1795 * z = scalbn(z,-e0)
1796 * for i = 0,1,2
1797 * x[i] = floor(z)
1798 * z = (z-x[i])*2**24
1799 *
1800 *
1801 * y[] ouput result in an array of double precision numbers.
1802 * The dimension of y[] is:
1803 * 24-bit precision 1
1804 * 53-bit precision 2
1805 * 64-bit precision 2
1806 * 113-bit precision 3
1807 * The actual value is the sum of them. Thus for 113-bit
1808 * precison, one may have to do something like:
1809 *
1810 * long double t,w,r_head, r_tail;
1811 * t = (long double)y[2] + (long double)y[1];
1812 * w = (long double)y[0];
1813 * r_head = t+w;
1814 * r_tail = w - (r_head - t);
1815 *
1816 * e0 The exponent of x[0]
1817 *
1818 * nx dimension of x[]
1819 *
1820 * prec an integer indicating the precision:
1821 * 0 24 bits (single)
1822 * 1 53 bits (double)
1823 * 2 64 bits (extended)
1824 * 3 113 bits (quad)
1825 *
1826 * two_over_pi[]
1827 * integer array, contains the (24*i)-th to (24*i+23)-th
1828 * bit of 2/pi after binary point. The corresponding
1829 * floating value is
1830 *
1831 * two_over_pi[i] * 2^(-24(i+1)).
1832 *
1833 * External function:
1834 * double scalbn(), floor();
1835 *
1836 *
1837 * Here is the description of some local variables:
1838 *
1839 * jk jk+1 is the initial number of terms of two_over_pi[] needed
1840 * in the computation. The recommended value is 2,3,4,
1841 * 6 for single, double, extended,and quad.
1842 *
1843 * jz local integer variable indicating the number of
1844 * terms of two_over_pi[] used.
1845 *
1846 * jx nx - 1
1847 *
1848 * jv index for pointing to the suitable two_over_pi[] for the
1849 * computation. In general, we want
1850 * ( 2^e0*x[0] * two_over_pi[jv-1]*2^(-24jv) )/8
1851 * is an integer. Thus
1852 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
1853 * Hence jv = max(0,(e0-3)/24).
1854 *
1855 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
1856 *
1857 * q[] double array with integral value, representing the
1858 * 24-bits chunk of the product of x and 2/pi.
1859 *
1860 * q0 the corresponding exponent of q[0]. Note that the
1861 * exponent for q[i] would be q0-24*i.
1862 *
1863 * PIo2[] double precision array, obtained by cutting pi/2
1864 * into 24 bits chunks.
1865 *
1866 * f[] two_over_pi[] in floating point
1867 *
1868 * iq[] integer array by breaking up q[] in 24-bits chunk.
1869 *
1870 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
1871 *
1872 * ih integer. If >0 it indicates q[] is >= 0.5, hence
1873 * it also indicates the *sign* of the result.
1874 *
1875 */
1876
1877 /*
1878 * Constants:
1879 * The hexadecimal values are the intended ones for the following
1880 * constants. The decimal values may be used, provided that the
1881 * compiler will convert from decimal to binary accurately enough
1882 * to produce the hexadecimal values shown.
1883 */
1884
1885 private static final double PIo2[] = {
1886 Double.longBitsToDouble(0x3ff921fb40000000L), /* 1.57079625129699707031e+00 */
1887 Double.longBitsToDouble(0x3e74442d00000000L), /* 7.54978941586159635335e-08 */
1888 Double.longBitsToDouble(0x3cf8469880000000L), /* 5.39030252995776476554e-15 */
1889 Double.longBitsToDouble(0x3b78cc5160000000L), /* 3.28200341580791294123e-22 */
1890 Double.longBitsToDouble(0x39f01b8380000000L), /* 1.27065575308067607349e-29 */
1891 Double.longBitsToDouble(0x387a252040000000L), /* 1.22933308981111328932e-36 */
1892 Double.longBitsToDouble(0x36e3822280000000L), /* 2.73370053816464559624e-44 */
1893 Double.longBitsToDouble(0x3569f31d00000000L) /* 2.16741683877804819444e-51 */
1894 };
1895
1896 private static final double twon24 = Double.longBitsToDouble(0x3E70000000000000L); /* 5.96046447753906250000e-08 */
1897
1898 private static int __kernel_rem_pio2(double x[], double y[], int e0, int nx)
1899 {
1900 int jz,jx,jv,jp,jk,carry,n,i,j,k,m,q0,ih;
1901 double z,fw;
1902 double f[] = new double[20];
1903 double q[] = new double[20];
1904 double fq[] = new double[20];
1905 int iq[] = new int[20];
1906
1907 /* initialize jk*/
1908 jk = 4;
1909 jp = jk;
1910
1911 /* determine jx,jv,q0, note that 3>q0 */
1912 jx = nx - 1;
1913 jv = (e0-3)/24;
1914 if (jv < 0) jv = 0;
1915 q0 = e0-24*(jv+1);
1916
1917 /* set up f[0] to f[jx+jk] where f[jx+jk] = two_over_pi[jv+jk] */
1918 j = jv - jx;
1919 m = jx + jk;
1920 for (i = 0; i <= m; i++ ,j++)
1921 f[i] = ((j<0) ? zero : (double)two_over_pi[j]);
1922
1923 /* compute q[0],q[1],...q[jk] */
1924 for (i = 0; i <= jk; i++) {
1925 for (j = 0, fw = 0.0; j <= jx; j++)
1926 fw += x[j]*f[jx+i-j];
1927 q[i] = fw;
1928 }
1929
1930 jz = jk;
1931
1932 while (true) { // recompute:
1933 /* distill q[] into iq[] reversingly */
1934 for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
1935 fw = (double)((int)(twon24* z));
1936 iq[i] = (int)(z-two24*fw);
1937 z = q[j-1]+fw;
1938 }
1939
1940 /* compute n */
1941 z = scalbn(z, q0); /* actual value of z */
1942 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
1943 n = (int) z;
1944 z -= (double)n;
1945 ih = 0;
1946 if (q0 > 0) { /* need iq[jz-1] to determine n */
1947 i = (iq[jz-1]>>(24-q0));
1948 n += i;
1949 iq[jz-1] -= i<<(24-q0);
1950 ih = iq[jz-1]>>(23-q0);
1951 } else if (q0==0) {
1952 ih = iq[jz-1]>>23;
1953 } else if(z>=0.5) {
1954 ih=2;
1955 }
1956 if (ih>0) { /* q > 0.5 */
1957 n += 1;
1958 carry = 0;
1959 for(i=0;i<jz ;i++) { /* compute 1-q */
1960 j = iq[i];
1961 if (carry == 0) {
1962 if (j != 0) {
1963 carry = 1;
1964 iq[i] = 0x1000000- j;
1965 }
1966 } else iq[i] = 0xffffff - j;
1967 }
1968 if (q0 > 0) { /* rare case: chance is 1 in 12 */
1969 switch (q0) {
1970 case 1:
1971 iq[jz-1] &= 0x7fffff; break;
1972 case 2:
1973 iq[jz-1] &= 0x3fffff; break;
1974 }
1975 }
1976 if (ih == 2) {
1977 z = one - z;
1978 if (carry != 0) z -= scalbn(one,q0);
1979 }
1980 }
1981
1982 /* check if recomputation is needed */
1983 if (z == zero) {
1984 j = 0;
1985 for (i = jz-1; i >= jk; i--)
1986 j |= iq[i];
1987 if (j == 0) { /* need recomputation */
1988 for (k = 1; iq[jk-k] == 0; k++); /* k = no. of terms needed */
1989 for (i = jz+1; i <= jz+k; i++) { /* add q[jz+1] to q[jz+k] */
1990 f[jx+i] = (double) two_over_pi[jv+i];
1991 for (j = 0, fw = 0.0; j <= jx; j++)
1992 fw += x[j]*f[jx+i-j];
1993 q[i] = fw;
1994 }
1995 jz += k;
1996 continue; //goto recompute;
1997 }
1998 }
1999 break;
2000 }
2001
2002 /* chop off zero terms */
2003 if (z == 0.0) {
2004 jz--;
2005 q0 -= 24;
2006 while (iq[jz] == 0) {
2007 jz--;
2008 q0 -= 24;
2009 }
2010 } else { /* break z into 24-bit if necessary */
2011 z = scalbn(z,-q0);
2012 if (z >= two24) {
2013 fw = (double)((int)(twon24*z));
2014 iq[jz] = (int)(z-two24*fw);
2015 jz++;
2016 q0 += 24;
2017 iq[jz] = (int)fw;
2018 } else iq[jz] = (int)z;
2019 }
2020
2021 /* convert integer "bit" chunk to floating-point value */
2022 fw = scalbn(one, q0);
2023 for (i = jz; i >= 0; i--) {
2024 q[i] = fw*(double)iq[i]; fw*=twon24;
2025 }
2026
2027 /* compute PIo2[0,...,jp]*q[jz,...,0] */
2028 for (i = jz; i >= 0; i--) {
2029 for (fw = 0.0, k = 0; k <= jp && k <= jz-i; k++)
2030 fw += PIo2[k]*q[i+k];
2031 fq[jz-i] = fw;
2032 }
2033
2034 /* compress fq[] into y[] */
2035 fw = 0.0;
2036 for (i = jz; i >= 0; i--)
2037 fw += fq[i];
2038 y[0] = (ih==0)? fw: -fw;
2039 fw = fq[0]-fw;
2040 for (i = 1; i <= jz; i++)
2041 fw += fq[i];
2042 y[1] = ((ih==0) ? fw: -fw);
2043 return n&7;
2044 }
2045
2046
2047 private static final double bp[] = {1.0, 1.5,};
2048 private static final double dp_h[] = { 0.0, Double.longBitsToDouble(0x3fe2b80340000000L)}; /* 5.84962487220764160156e-01 */
2049 private static final double dp_l[] = { 0.0, Double.longBitsToDouble(0x3e4cfdeb43cfd006L)}; /* 1.35003920212974897128e-08 */
2050 private static final double two53 = Double.longBitsToDouble(0x4340000000000000L); /* 9007199254740992.0 */
2051 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
2052 private static final double L1 = Double.longBitsToDouble(0x3fe3333333333303L); /* 5.99999999999994648725e-01 */
2053 private static final double L2 = Double.longBitsToDouble(0x3fdb6db6db6fabffL); /* 4.28571428578550184252e-01 */
2054 private static final double L3 = Double.longBitsToDouble(0x3fd55555518f264dL); /* 3.33333329818377432918e-01 */
2055 private static final double L4 = Double.longBitsToDouble(0x3fd17460a91d4101L); /* 2.72728123808534006489e-01 */
2056 private static final double L5 = Double.longBitsToDouble(0x3fcd864a93c9db65L); /* 2.30660745775561754067e-01 */
2057 private static final double L6 = Double.longBitsToDouble(0x3fca7e284a454eefL); /* 2.06975017800338417784e-01 */
2058 private static final double lg2 = Double.longBitsToDouble(0x3fe62e42fefa39efL); /* 6.93147180559945286227e-01 */
2059 private static final double lg2_h = Double.longBitsToDouble(0x3fe62e4300000000L); /* 6.93147182464599609375e-01 */
2060 private static final double lg2_l = -1.90465429995776804525e-09; /* 0xbe205c610ca86c39 */
2061 private static final double ovt = 8.0085662595372944372e-17; /* -(1024-log2(ovfl+.5ulp)) */
2062 private static final double cp = Double.longBitsToDouble(0x3feec709dc3a03fdL); /* 9.61796693925975554329e-01 = 2/(3ln2) */
2063 private static final double cp_h = Double.longBitsToDouble(0x3feec709e0000000L); /* 9.61796700954437255859e-01 = (float)cp */
2064 private static final double cp_l = Double.longBitsToDouble(0xbe3e2fe0145b01f5L); /* -7.02846165095275826516e-09 = tail of cp_h*/
2065 private static final double ivln2 = Double.longBitsToDouble(0x3ff71547652b82feL); /* 1.44269504088896338700e+00 = 1/ln2 */
2066 private static final double ivln2_h = Double.longBitsToDouble(0x3ff7154760000000L); /* 1.44269502162933349609e+00 = 24b 1/ln2*/
2067 private static final double ivln2_l = Double.longBitsToDouble(0x3e54ae0bf85ddf44L); /* 1.92596299112661746887e-08 = 1/ln2 tail*/
2068
2069
2070 /*
2071 * Returns x to the power y
2072 *
2073 * n
2074 * Method: Let x = 2 * (1+f)
2075 * 1. Compute and return log2(x) in two pieces:
2076 * log2(x) = w1 + w2,
2077 * where w1 has 53-24 = 29 bit trailing zeros.
2078 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
2079 * arithmetic, where |y'|<=0.5.
2080 * 3. Return x**y = 2**n*exp(y'*log2)
2081 *
2082 * Special cases:
2083 * 1. (anything) ** 0 is 1
2084 * 2. (anything) ** 1 is itself
2085 * 3. (anything) ** NAN is NAN
2086 * 4. NAN ** (anything except 0) is NAN
2087 * 5. +-(|x| > 1) ** +INF is +INF
2088 * 6. +-(|x| > 1) ** -INF is +0
2089 * 7. +-(|x| < 1) ** +INF is +0
2090 * 8. +-(|x| < 1) ** -INF is +INF
2091 * 9. +-1 ** +-INF is NAN
2092 * 10. +0 ** (+anything except 0, NAN) is +0
2093 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
2094 * 12. +0 ** (-anything except 0, NAN) is +INF
2095 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
2096 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
2097 * 15. +INF ** (+anything except 0,NAN) is +INF
2098 * 16. +INF ** (-anything except 0,NAN) is +0
2099 * 17. -INF ** (anything) = -0 ** (-anything)
2100 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
2101 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
2102 *
2103 * Accuracy:
2104 * pow(x,y) returns x**y nearly rounded. In particular
2105 * pow(integer,integer)
2106 * always returns the correct integer provided it is
2107 * representable.
2108 */
2109
2110 /***
2111 * Returns x to the power y.
2112 * @param x The base.
2113 * @param y The exponent.
2114 * @return x to the power y.
2115 */
2116 public static double pow(double x, double y)
2117 {
2118 double z, ax, z_h, z_l, p_h, p_l;
2119 double y1, t1, t2, r, s, t, u, v, w;
2120 int i, j, k, yisint, n;
2121 int hx, hy, ix, iy;
2122 int lx, ly;
2123
2124 hx = __HI(x);
2125 lx = __LO(x);
2126 hy = __HI(y);
2127 ly = __LO(y);
2128 ix = hx & 0x7fffffff;
2129 iy = hy & 0x7fffffff;
2130
2131 /* y==zero: x**0 = 1 */
2132 if ((iy|ly) == 0) return one;
2133
2134 /* +-NaN return x+y */
2135 if (ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
2136 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
2137 return x+y;
2138
2139 /* determine if y is an odd int when x < 0
2140 * yisint = 0 ... y is not an integer
2141 * yisint = 1 ... y is an odd int
2142 * yisint = 2 ... y is an even int
2143 */
2144 yisint = 0;
2145 if (hx < 0) {
2146 if (iy >= 0x43400000)
2147 yisint = 2; /* even integer y */
2148 else if(iy >= 0x3ff00000) {
2149 k = (iy>>20)-0x3ff; /* exponent */
2150 if (k > 20) {
2151 j = ly>>>(52-k);
2152 if ((j<<(52-k)) == ly)
2153 yisint = 2-(j&1);
2154 } else if (ly == 0) {
2155 j = iy >> (20-k);
2156 if ((j<<(20-k)) == iy)
2157 yisint = 2-(j&1);
2158 }
2159 }
2160 }
2161
2162 /* special value of y */
2163 if (ly==0) {
2164 if (iy == 0x7ff00000) { /* y is +-inf */
2165 if (((ix-0x3ff00000)|lx) == 0)
2166 return y - y; /* inf**+-1 is NaN */
2167 else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
2168 return (hy>=0)? y: zero;
2169 else /* (|x|<1)**-,+inf = inf,0 */
2170 return (hy<0)?-y: zero;
2171 }
2172 if (iy == 0x3ff00000) { /* y is +-1 */
2173 if (hy < 0)
2174 return one/x;
2175 else
2176 return x;
2177 }
2178 if (hy == 0x40000000)
2179 return x*x; /* y is 2 */
2180 if (hy == 0x3fe00000) { /* y is 0.5 */
2181 if(hx>=0) /* x >= +0 */
2182 return sqrt(x);
2183 }
2184 }
2185
2186 ax = abs(x);
2187 /* special value of x */
2188 if (lx == 0) {
2189 if (ix==0x7ff00000 || ix==0 || ix==0x3ff00000) {
2190 z = ax; /*x is +-0,+-inf,+-1*/
2191 if (hy < 0 )
2192 z = one/z; /* z = (1/|x|) */
2193 if (hx < 0) {
2194 if (((ix-0x3ff00000)|yisint) == 0) {
2195 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
2196 } else if (yisint==1)
2197 z = -z; /* (x<0)**odd = -(|x|**odd) */
2198 }
2199 return z;
2200 }
2201 }
2202
2203 /* (x<0)**(non-int) is NaN */
2204 if ((((hx>>31)+1)|yisint) == 0)
2205 return (x-x)/(x-x);
2206
2207 /* |y| is huge */
2208 if (iy > 0x41e00000) { /* if |y| > 2**31 */
2209 if (iy > 0x43f00000){ /* if |y| > 2**64, must o/uflow */
2210 if (ix <= 0x3fefffff)
2211 return ((hy < 0) ? huge*huge : tiny*tiny);
2212 if (ix >= 0x3ff00000)
2213 return ((hy > 0) ? huge*huge : tiny*tiny);
2214 }
2215
2216 /* over/underflow if x is not close to one */
2217 if (ix < 0x3fefffff)
2218 return ((hy < 0) ? huge*huge : tiny*tiny);
2219 if (ix > 0x3ff00000)
2220 return ((hy > 0) ? huge*huge : tiny*tiny);
2221 /* now |1-x| is tiny <= 2**-20, suffice to compute
2222 log(x) by x-x^2/2+x^3/3-x^4/4 */
2223 t = x-1; /* t has 20 trailing zeros */
2224 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
2225 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
2226 v = t*ivln2_l-w*ivln2;
2227 t1 = u+v;
2228 t1 = setLO(t1, 0);
2229 t2 = v-(t1-u);
2230 } else {
2231 double s2,s_h,s_l,t_h,t_l;
2232 n = 0;
2233 /* take care subnormal number */
2234 if (ix<0x00100000) {
2235 ax *= two53;
2236 n -= 53;
2237 ix = __HI(ax);
2238 }
2239 n += ((ix)>>20) - 0x3ff;
2240 j = ix & 0x000fffff;
2241 /* determine interval */
2242 ix = j|0x3ff00000; /* normalize ix */
2243 if (j <= 0x3988E)
2244 k=0; /* |x|<sqrt(3/2) */
2245 else if (j < 0xBB67A)
2246 k=1; /* |x|<sqrt(3) */
2247 else {
2248 k = 0;
2249 n+=1;
2250 ix -= 0x00100000;
2251 }
2252 ax = setHI(ax, ix);
2253
2254 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
2255 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
2256 v = one/(ax+bp[k]);
2257 s = u*v;
2258 s_h = s;
2259 s_h = setLO(s_h, 0);
2260 /* t_h=ax+bp[k] High */
2261 t_h = zero;
2262 t_h = setHI(t_h, ((ix>>1)|0x20000000)+0x00080000+(k<<18));
2263 t_l = ax - (t_h-bp[k]);
2264 s_l = v*((u-s_h*t_h)-s_h*t_l);
2265 /* compute log(ax) */
2266 s2 = s*s;
2267 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
2268 r += s_l*(s_h+s);
2269 s2 = s_h*s_h;
2270 t_h = 3.0+s2+r;
2271 t_h = setLO(t_h, 0);
2272 t_l = r-((t_h-3.0)-s2);
2273 /* u+v = s*(1+...) */
2274 u = s_h*t_h;
2275 v = s_l*t_h+t_l*s;
2276 /* 2/(3log2)*(s+...) */
2277 p_h = u+v;
2278 p_h = setLO(p_h, 0);
2279 p_l = v-(p_h-u);
2280 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
2281 z_l = cp_l*p_h+p_l*cp+dp_l[k];
2282 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
2283 t = (double)n;
2284 t1 = (((z_h+z_l)+dp_h[k])+t);
2285 t1 = setLO(t1, 0);
2286 t2 = z_l-(((t1-t)-dp_h[k])-z_h);
2287 }
2288
2289 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
2290 if((((hx>>31)+1)|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
2291
2292 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
2293 y1 = y;
2294 y1 = setLO(y1, 0);
2295 p_l = (y-y1)*t1+y*t2;
2296 p_h = y1*t1;
2297 z = p_l+p_h;
2298 j = __HI(z);
2299 i = __LO(z);
2300 if (j >= 0x40900000) { /* z >= 1024 */
2301 if (((j-0x40900000)|i) != 0) /* if z > 1024 */
2302 return s*huge*huge; /* overflow */
2303 else {
2304 if (p_l+ovt > z-p_h)
2305 return s*huge*huge; /* overflow */
2306 }
2307 } else if((j&0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */
2308 if (((j-0xc090cc00)|i) != 0) /* z < -1075 */
2309 return s*tiny*tiny; /* underflow */
2310 else {
2311 if (p_l <= z-p_h) return s*tiny*tiny; /* underflow */
2312 }
2313 }
2314 /*
2315 * compute 2**(p_h+p_l)
2316 */
2317 i = j&0x7fffffff;
2318 k = (i>>20)-0x3ff;
2319 n = 0;
2320 if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
2321 n = j + (0x00100000>>(k+1));
2322 k = ((n&0x7fffffff)>>20) - 0x3ff; /* new k for n */
2323 t = zero;
2324 t = setHI(t, (n&~(0x000fffff>>k)));
2325 n = ((n&0x000fffff)|0x00100000)>>(20-k);
2326 if (j < 0) n = -n;
2327 p_h -= t;
2328 }
2329 t = p_l + p_h;
2330 t = setLO(t,0);
2331 u = t*lg2_h;
2332 v = (p_l-(t-p_h))*lg2 + t*lg2_l;
2333 z = u + v;
2334 w = v - (z-u);
2335 t = z*z;
2336 t1= z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
2337 r = (z*t1)/(t1-2.0) - (w+z*w);
2338 z = one - (r-z);
2339 j = __HI(z);
2340 j += (n<<20);
2341 if ((j>>20) <= 0) {
2342 z = scalbn(z,n); /* subnormal output */
2343 } else {
2344 i = __HI(z);
2345 i += (n << 20);
2346 z = setHI(z, i);
2347 }
2348 return s*z;
2349 }
2350
2351
2352 /*
2353 * copysign(double x, double y)
2354 * copysign(x,y) returns a value with the magnitude of x and
2355 * with the sign bit of y.
2356 */
2357
2358 private static double copysign(double x, double y)
2359 {
2360 long ix = Double.doubleToLongBits(x);
2361 long iy = Double.doubleToLongBits(y);
2362 ix = (0x7fffffffffffffffL*ix) | (0x8000000000000000L&iy);
2363 return Double.longBitsToDouble(ix);
2364 }
2365
2366 private static final double two54 = Double.longBitsToDouble(0x4350000000000000L); /* 1.80143985094819840000e+16 */
2367 private static final double twom54 = Double.longBitsToDouble(0x3c90000000000000L); /* 5.55111512312578270212e-17 */
2368
2369
2370 /*
2371 * scalbn (double x, int n)
2372 * scalbn(x,n) returns x* 2**n computed by exponent
2373 * manipulation rather than by actually performing an
2374 * exponentiation or a multiplication.
2375 */
2376 private static double scalbn (double x, int n)
2377 {
2378 int k, hx, lx;
2379 hx = __HI(x);
2380 lx = __LO(x);
2381 k = (hx&0x7ff00000)>>20; /* extract exponent */
2382 if (k == 0) { /* 0 or subnormal x */
2383 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
2384 x *= two54;
2385 hx = __HI(x);
2386 k = ((hx&0x7ff00000)>>20) - 54;
2387 if (n < -50000) return tiny*x; /*underflow*/
2388 }
2389 if (k == 0x7ff) return x+x; /* NaN or Inf */
2390 k = k+n;
2391 if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
2392 if (k > 0) {
2393 /* normal result */
2394 return setHI(x, (hx&0x800fffff)|(k<<20));
2395 }
2396 if (k <= -54) {
2397 if (n > 50000) /* in case integer overflow in n+k */
2398 return huge*copysign(huge,x); /*overflow*/
2399 } else {
2400 return tiny*copysign(tiny,x); /*underflow*/
2401 }
2402 k += 54; /* subnormal result */
2403 return twom54 * setHI(x, (hx&0x800fffff)|(k<<20));
2404 }
2405
2406
2407 private static double set(int newHiPart, int newLowPart)
2408 {
2409 return Double.longBitsToDouble((((long)newHiPart)<<32) | newLowPart);
2410 }
2411
2412
2413 private static double setLO(double x, int newLowPart)
2414 {
2415 long lx = Double.doubleToLongBits(x);
2416 lx &= 0xFFFFFFFF00000000L;
2417 lx |= newLowPart;
2418 return Double.longBitsToDouble(lx);
2419 }
2420
2421
2422 private static double setHI(double x, int newHiPart)
2423 {
2424 long lx = Double.doubleToLongBits(x);
2425 lx &= 0x00000000FFFFFFFFL;
2426 lx |= (((long)newHiPart)<<32);
2427 return Double.longBitsToDouble(lx);
2428 }
2429
2430 private static int __HI(double x)
2431 {
2432 return (int)(0xFFFFFFFF&(Double.doubleToLongBits(x)>>32));
2433 }
2434
2435 private static int __LO(double x)
2436 {
2437 return (int)(0xFFFFFFFF&Double.doubleToLongBits(x));
2438 }
2439
2440 }