View Javadoc

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 }