← Back to branch summary

~ubuntu-branches/ubuntu/trusty/monodevelop/trusty-proposed

« back to all changes in this revision

Viewing changes to external/ikvm/classpath/ikvm/internal/JMath.java

  • Committer: Package Import Robot
  • Author(s): Jo Shields
  • Date: 2013-05-12 09:46:03 UTC
  • mto: This revision was merged to the branch mainline in revision 29.
  • Revision ID: package-import@ubuntu.com-20130512094603-mad323bzcxvmcam0
Tags: upstream-4.0.5+dfsg
Import upstream version 4.0.5+dfsg

Show diffs side-by-side

added added

removed removed

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