View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math3.util;
18  
19  import java.io.PrintStream;
20  
21  /**
22   * Faster, more accurate, portable alternative to {@link Math} and
23   * {@link StrictMath} for large scale computation.
24   * <p>
25   * FastMath is a drop-in replacement for both Math and StrictMath. This
26   * means that for any method in Math (say {@code Math.sin(x)} or
27   * {@code Math.cbrt(y)}), user can directly change the class and use the
28   * methods as is (using {@code FastMath.sin(x)} or {@code FastMath.cbrt(y)}
29   * in the previous example).
30   * </p>
31   * <p>
32   * FastMath speed is achieved by relying heavily on optimizing compilers
33   * to native code present in many JVMs today and use of large tables.
34   * The larger tables are lazily initialised on first use, so that the setup
35   * time does not penalise methods that don't need them.
36   * </p>
37   * <p>
38   * Note that FastMath is
39   * extensively used inside Apache Commons Math, so by calling some algorithms,
40   * the overhead when the the tables need to be intialised will occur
41   * regardless of the end-user calling FastMath methods directly or not.
42   * Performance figures for a specific JVM and hardware can be evaluated by
43   * running the FastMathTestPerformance tests in the test directory of the source
44   * distribution.
45   * </p>
46   * <p>
47   * FastMath accuracy should be mostly independent of the JVM as it relies only
48   * on IEEE-754 basic operations and on embedded tables. Almost all operations
49   * are accurate to about 0.5 ulp throughout the domain range. This statement,
50   * of course is only a rough global observed behavior, it is <em>not</em> a
51   * guarantee for <em>every</em> double numbers input (see William Kahan's <a
52   * href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">Table
53   * Maker's Dilemma</a>).
54   * </p>
55   * <p>
56   * FastMath additionally implements the following methods not found in Math/StrictMath:
57   * <ul>
58   * <li>{@link #asinh(double)}</li>
59   * <li>{@link #acosh(double)}</li>
60   * <li>{@link #atanh(double)}</li>
61   * </ul>
62   * The following methods are found in Math/StrictMath since 1.6 only, they are provided
63   * by FastMath even in 1.5 Java virtual machines
64   * <ul>
65   * <li>{@link #copySign(double, double)}</li>
66   * <li>{@link #getExponent(double)}</li>
67   * <li>{@link #nextAfter(double,double)}</li>
68   * <li>{@link #nextUp(double)}</li>
69   * <li>{@link #scalb(double, int)}</li>
70   * <li>{@link #copySign(float, float)}</li>
71   * <li>{@link #getExponent(float)}</li>
72   * <li>{@link #nextAfter(float,double)}</li>
73   * <li>{@link #nextUp(float)}</li>
74   * <li>{@link #scalb(float, int)}</li>
75   * </ul>
76   * </p>
77   * @version $Id: FastMath.java 1591835 2014-05-02 09:04:01Z tn $
78   * @since 2.2
79   */
80  public class FastMath {
81      /** Archimede's constant PI, ratio of circle circumference to diameter. */
82      public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
83  
84      /** Napier's constant e, base of the natural logarithm. */
85      public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
86  
87      /** Index of exp(0) in the array of integer exponentials. */
88      static final int EXP_INT_TABLE_MAX_INDEX = 750;
89      /** Length of the array of integer exponentials. */
90      static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
91      /** Logarithm table length. */
92      static final int LN_MANT_LEN = 1024;
93      /** Exponential fractions table length. */
94      static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024
95  
96      /** StrictMath.log(Double.MAX_VALUE): {@value} */
97      private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
98  
99      /** Indicator for tables initialization.
100      * <p>
101      * This compile-time constant should be set to true only if one explicitly
102      * wants to compute the tables at class loading time instead of using the
103      * already computed ones provided as literal arrays below.
104      * </p>
105      */
106     private static final boolean RECOMPUTE_TABLES_AT_RUNTIME = false;
107 
108     /** log(2) (high bits). */
109     private static final double LN_2_A = 0.693147063255310059;
110 
111     /** log(2) (low bits). */
112     private static final double LN_2_B = 1.17304635250823482e-7;
113 
114     /** Coefficients for log, when input 0.99 < x < 1.01. */
115     private static final double LN_QUICK_COEF[][] = {
116         {1.0, 5.669184079525E-24},
117         {-0.25, -0.25},
118         {0.3333333134651184, 1.986821492305628E-8},
119         {-0.25, -6.663542893624021E-14},
120         {0.19999998807907104, 1.1921056801463227E-8},
121         {-0.1666666567325592, -7.800414592973399E-9},
122         {0.1428571343421936, 5.650007086920087E-9},
123         {-0.12502530217170715, -7.44321345601866E-11},
124         {0.11113807559013367, 9.219544613762692E-9},
125     };
126 
127     /** Coefficients for log in the range of 1.0 < x < 1.0 + 2^-10. */
128     private static final double LN_HI_PREC_COEF[][] = {
129         {1.0, -6.032174644509064E-23},
130         {-0.25, -0.25},
131         {0.3333333134651184, 1.9868161777724352E-8},
132         {-0.2499999701976776, -2.957007209750105E-8},
133         {0.19999954104423523, 1.5830993332061267E-10},
134         {-0.16624879837036133, -2.6033824355191673E-8}
135     };
136 
137     /** Sine, Cosine, Tangent tables are for 0, 1/8, 2/8, ... 13/8 = PI/2 approx. */
138     private static final int SINE_TABLE_LEN = 14;
139 
140     /** Sine table (high bits). */
141     private static final double SINE_TABLE_A[] =
142         {
143         +0.0d,
144         +0.1246747374534607d,
145         +0.24740394949913025d,
146         +0.366272509098053d,
147         +0.4794255495071411d,
148         +0.5850973129272461d,
149         +0.6816387176513672d,
150         +0.7675435543060303d,
151         +0.8414709568023682d,
152         +0.902267575263977d,
153         +0.9489846229553223d,
154         +0.9808930158615112d,
155         +0.9974949359893799d,
156         +0.9985313415527344d,
157     };
158 
159     /** Sine table (low bits). */
160     private static final double SINE_TABLE_B[] =
161         {
162         +0.0d,
163         -4.068233003401932E-9d,
164         +9.755392680573412E-9d,
165         +1.9987994582857286E-8d,
166         -1.0902938113007961E-8d,
167         -3.9986783938944604E-8d,
168         +4.23719669792332E-8d,
169         -5.207000323380292E-8d,
170         +2.800552834259E-8d,
171         +1.883511811213715E-8d,
172         -3.5997360512765566E-9d,
173         +4.116164446561962E-8d,
174         +5.0614674548127384E-8d,
175         -1.0129027912496858E-9d,
176     };
177 
178     /** Cosine table (high bits). */
179     private static final double COSINE_TABLE_A[] =
180         {
181         +1.0d,
182         +0.9921976327896118d,
183         +0.9689123630523682d,
184         +0.9305076599121094d,
185         +0.8775825500488281d,
186         +0.8109631538391113d,
187         +0.7316888570785522d,
188         +0.6409968137741089d,
189         +0.5403022766113281d,
190         +0.4311765432357788d,
191         +0.3153223395347595d,
192         +0.19454771280288696d,
193         +0.07073719799518585d,
194         -0.05417713522911072d,
195     };
196 
197     /** Cosine table (low bits). */
198     private static final double COSINE_TABLE_B[] =
199         {
200         +0.0d,
201         +3.4439717236742845E-8d,
202         +5.865827662008209E-8d,
203         -3.7999795083850525E-8d,
204         +1.184154459111628E-8d,
205         -3.43338934259355E-8d,
206         +1.1795268640216787E-8d,
207         +4.438921624363781E-8d,
208         +2.925681159240093E-8d,
209         -2.6437112632041807E-8d,
210         +2.2860509143963117E-8d,
211         -4.813899778443457E-9d,
212         +3.6725170580355583E-9d,
213         +2.0217439756338078E-10d,
214     };
215 
216 
217     /** Tangent table, used by atan() (high bits). */
218     private static final double TANGENT_TABLE_A[] =
219         {
220         +0.0d,
221         +0.1256551444530487d,
222         +0.25534194707870483d,
223         +0.3936265707015991d,
224         +0.5463024377822876d,
225         +0.7214844226837158d,
226         +0.9315965175628662d,
227         +1.1974215507507324d,
228         +1.5574076175689697d,
229         +2.092571258544922d,
230         +3.0095696449279785d,
231         +5.041914939880371d,
232         +14.101419448852539d,
233         -18.430862426757812d,
234     };
235 
236     /** Tangent table, used by atan() (low bits). */
237     private static final double TANGENT_TABLE_B[] =
238         {
239         +0.0d,
240         -7.877917738262007E-9d,
241         -2.5857668567479893E-8d,
242         +5.2240336371356666E-9d,
243         +5.206150291559893E-8d,
244         +1.8307188599677033E-8d,
245         -5.7618793749770706E-8d,
246         +7.848361555046424E-8d,
247         +1.0708593250394448E-7d,
248         +1.7827257129423813E-8d,
249         +2.893485277253286E-8d,
250         +3.1660099222737955E-7d,
251         +4.983191803254889E-7d,
252         -3.356118100840571E-7d,
253     };
254 
255     /** Bits of 1/(2*pi), need for reducePayneHanek(). */
256     private static final long RECIP_2PI[] = new long[] {
257         (0x28be60dbL << 32) | 0x9391054aL,
258         (0x7f09d5f4L << 32) | 0x7d4d3770L,
259         (0x36d8a566L << 32) | 0x4f10e410L,
260         (0x7f9458eaL << 32) | 0xf7aef158L,
261         (0x6dc91b8eL << 32) | 0x909374b8L,
262         (0x01924bbaL << 32) | 0x82746487L,
263         (0x3f877ac7L << 32) | 0x2c4a69cfL,
264         (0xba208d7dL << 32) | 0x4baed121L,
265         (0x3a671c09L << 32) | 0xad17df90L,
266         (0x4e64758eL << 32) | 0x60d4ce7dL,
267         (0x272117e2L << 32) | 0xef7e4a0eL,
268         (0xc7fe25ffL << 32) | 0xf7816603L,
269         (0xfbcbc462L << 32) | 0xd6829b47L,
270         (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
271         (0xd3d18fd9L << 32) | 0xa797fa8bL,
272         (0x5d49eeb1L << 32) | 0xfaf97c5eL,
273         (0xcf41ce7dL << 32) | 0xe294a4baL,
274          0x9afed7ecL << 32  };
275 
276     /** Bits of pi/4, need for reducePayneHanek(). */
277     private static final long PI_O_4_BITS[] = new long[] {
278         (0xc90fdaa2L << 32) | 0x2168c234L,
279         (0xc4c6628bL << 32) | 0x80dc1cd1L };
280 
281     /** Eighths.
282      * This is used by sinQ, because its faster to do a table lookup than
283      * a multiply in this time-critical routine
284      */
285     private static final double EIGHTHS[] = {0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625};
286 
287     /** Table of 2^((n+2)/3) */
288     private static final double CBRTTWO[] = { 0.6299605249474366,
289                                             0.7937005259840998,
290                                             1.0,
291                                             1.2599210498948732,
292                                             1.5874010519681994 };
293 
294     /*
295      *  There are 52 bits in the mantissa of a double.
296      *  For additional precision, the code splits double numbers into two parts,
297      *  by clearing the low order 30 bits if possible, and then performs the arithmetic
298      *  on each half separately.
299      */
300 
301     /**
302      * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
303      * Equivalent to 2^30.
304      */
305     private static final long HEX_40000000 = 0x40000000L; // 1073741824L
306 
307     /** Mask used to clear low order 30 bits */
308     private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L;
309 
310     /** Mask used to clear the non-sign part of an int. */
311     private static final int MASK_NON_SIGN_INT = 0x7fffffff;
312 
313     /** Mask used to clear the non-sign part of a long. */
314     private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
315 
316     /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
317     private static final double TWO_POWER_52 = 4503599627370496.0;
318     /** 2^53 - double numbers this large must be even. */
319     private static final double TWO_POWER_53 = 2 * TWO_POWER_52;
320 
321     /** Constant: {@value}. */
322     private static final double F_1_3 = 1d / 3d;
323     /** Constant: {@value}. */
324     private static final double F_1_5 = 1d / 5d;
325     /** Constant: {@value}. */
326     private static final double F_1_7 = 1d / 7d;
327     /** Constant: {@value}. */
328     private static final double F_1_9 = 1d / 9d;
329     /** Constant: {@value}. */
330     private static final double F_1_11 = 1d / 11d;
331     /** Constant: {@value}. */
332     private static final double F_1_13 = 1d / 13d;
333     /** Constant: {@value}. */
334     private static final double F_1_15 = 1d / 15d;
335     /** Constant: {@value}. */
336     private static final double F_1_17 = 1d / 17d;
337     /** Constant: {@value}. */
338     private static final double F_3_4 = 3d / 4d;
339     /** Constant: {@value}. */
340     private static final double F_15_16 = 15d / 16d;
341     /** Constant: {@value}. */
342     private static final double F_13_14 = 13d / 14d;
343     /** Constant: {@value}. */
344     private static final double F_11_12 = 11d / 12d;
345     /** Constant: {@value}. */
346     private static final double F_9_10 = 9d / 10d;
347     /** Constant: {@value}. */
348     private static final double F_7_8 = 7d / 8d;
349     /** Constant: {@value}. */
350     private static final double F_5_6 = 5d / 6d;
351     /** Constant: {@value}. */
352     private static final double F_1_2 = 1d / 2d;
353     /** Constant: {@value}. */
354     private static final double F_1_4 = 1d / 4d;
355 
356     /**
357      * Private Constructor
358      */
359     private FastMath() {}
360 
361     // Generic helper methods
362 
363     /**
364      * Get the high order bits from the mantissa.
365      * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
366      *
367      * @param d the value to split
368      * @return the high order part of the mantissa
369      */
370     private static double doubleHighPart(double d) {
371         if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN){
372             return d; // These are un-normalised - don't try to convert
373         }
374         long xl = Double.doubleToRawLongBits(d); // can take raw bits because just gonna convert it back
375         xl &= MASK_30BITS; // Drop low order bits
376         return Double.longBitsToDouble(xl);
377     }
378 
379     /** Compute the square root of a number.
380      * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
381      * @param a number on which evaluation is done
382      * @return square root of a
383      */
384     public static double sqrt(final double a) {
385         return Math.sqrt(a);
386     }
387 
388     /** Compute the hyperbolic cosine of a number.
389      * @param x number on which evaluation is done
390      * @return hyperbolic cosine of x
391      */
392     public static double cosh(double x) {
393       if (x != x) {
394           return x;
395       }
396 
397       // cosh[z] = (exp(z) + exp(-z))/2
398 
399       // for numbers with magnitude 20 or so,
400       // exp(-z) can be ignored in comparison with exp(z)
401 
402       if (x > 20) {
403           if (x >= LOG_MAX_VALUE) {
404               // Avoid overflow (MATH-905).
405               final double t = exp(0.5 * x);
406               return (0.5 * t) * t;
407           } else {
408               return 0.5 * exp(x);
409           }
410       } else if (x < -20) {
411           if (x <= -LOG_MAX_VALUE) {
412               // Avoid overflow (MATH-905).
413               final double t = exp(-0.5 * x);
414               return (0.5 * t) * t;
415           } else {
416               return 0.5 * exp(-x);
417           }
418       }
419 
420       final double hiPrec[] = new double[2];
421       if (x < 0.0) {
422           x = -x;
423       }
424       exp(x, 0.0, hiPrec);
425 
426       double ya = hiPrec[0] + hiPrec[1];
427       double yb = -(ya - hiPrec[0] - hiPrec[1]);
428 
429       double temp = ya * HEX_40000000;
430       double yaa = ya + temp - temp;
431       double yab = ya - yaa;
432 
433       // recip = 1/y
434       double recip = 1.0/ya;
435       temp = recip * HEX_40000000;
436       double recipa = recip + temp - temp;
437       double recipb = recip - recipa;
438 
439       // Correct for rounding in division
440       recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
441       // Account for yb
442       recipb += -yb * recip * recip;
443 
444       // y = y + 1/y
445       temp = ya + recipa;
446       yb += -(temp - ya - recipa);
447       ya = temp;
448       temp = ya + recipb;
449       yb += -(temp - ya - recipb);
450       ya = temp;
451 
452       double result = ya + yb;
453       result *= 0.5;
454       return result;
455     }
456 
457     /** Compute the hyperbolic sine of a number.
458      * @param x number on which evaluation is done
459      * @return hyperbolic sine of x
460      */
461     public static double sinh(double x) {
462       boolean negate = false;
463       if (x != x) {
464           return x;
465       }
466 
467       // sinh[z] = (exp(z) - exp(-z) / 2
468 
469       // for values of z larger than about 20,
470       // exp(-z) can be ignored in comparison with exp(z)
471 
472       if (x > 20) {
473           if (x >= LOG_MAX_VALUE) {
474               // Avoid overflow (MATH-905).
475               final double t = exp(0.5 * x);
476               return (0.5 * t) * t;
477           } else {
478               return 0.5 * exp(x);
479           }
480       } else if (x < -20) {
481           if (x <= -LOG_MAX_VALUE) {
482               // Avoid overflow (MATH-905).
483               final double t = exp(-0.5 * x);
484               return (-0.5 * t) * t;
485           } else {
486               return -0.5 * exp(-x);
487           }
488       }
489 
490       if (x == 0) {
491           return x;
492       }
493 
494       if (x < 0.0) {
495           x = -x;
496           negate = true;
497       }
498 
499       double result;
500 
501       if (x > 0.25) {
502           double hiPrec[] = new double[2];
503           exp(x, 0.0, hiPrec);
504 
505           double ya = hiPrec[0] + hiPrec[1];
506           double yb = -(ya - hiPrec[0] - hiPrec[1]);
507 
508           double temp = ya * HEX_40000000;
509           double yaa = ya + temp - temp;
510           double yab = ya - yaa;
511 
512           // recip = 1/y
513           double recip = 1.0/ya;
514           temp = recip * HEX_40000000;
515           double recipa = recip + temp - temp;
516           double recipb = recip - recipa;
517 
518           // Correct for rounding in division
519           recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
520           // Account for yb
521           recipb += -yb * recip * recip;
522 
523           recipa = -recipa;
524           recipb = -recipb;
525 
526           // y = y + 1/y
527           temp = ya + recipa;
528           yb += -(temp - ya - recipa);
529           ya = temp;
530           temp = ya + recipb;
531           yb += -(temp - ya - recipb);
532           ya = temp;
533 
534           result = ya + yb;
535           result *= 0.5;
536       }
537       else {
538           double hiPrec[] = new double[2];
539           expm1(x, hiPrec);
540 
541           double ya = hiPrec[0] + hiPrec[1];
542           double yb = -(ya - hiPrec[0] - hiPrec[1]);
543 
544           /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
545           double denom = 1.0 + ya;
546           double denomr = 1.0 / denom;
547           double denomb = -(denom - 1.0 - ya) + yb;
548           double ratio = ya * denomr;
549           double temp = ratio * HEX_40000000;
550           double ra = ratio + temp - temp;
551           double rb = ratio - ra;
552 
553           temp = denom * HEX_40000000;
554           double za = denom + temp - temp;
555           double zb = denom - za;
556 
557           rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
558 
559           // Adjust for yb
560           rb += yb*denomr;                        // numerator
561           rb += -ya * denomb * denomr * denomr;   // denominator
562 
563           // y = y - 1/y
564           temp = ya + ra;
565           yb += -(temp - ya - ra);
566           ya = temp;
567           temp = ya + rb;
568           yb += -(temp - ya - rb);
569           ya = temp;
570 
571           result = ya + yb;
572           result *= 0.5;
573       }
574 
575       if (negate) {
576           result = -result;
577       }
578 
579       return result;
580     }
581 
582     /** Compute the hyperbolic tangent of a number.
583      * @param x number on which evaluation is done
584      * @return hyperbolic tangent of x
585      */
586     public static double tanh(double x) {
587       boolean negate = false;
588 
589       if (x != x) {
590           return x;
591       }
592 
593       // tanh[z] = sinh[z] / cosh[z]
594       // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
595       // = (exp(2x) - 1) / (exp(2x) + 1)
596 
597       // for magnitude > 20, sinh[z] == cosh[z] in double precision
598 
599       if (x > 20.0) {
600           return 1.0;
601       }
602 
603       if (x < -20) {
604           return -1.0;
605       }
606 
607       if (x == 0) {
608           return x;
609       }
610 
611       if (x < 0.0) {
612           x = -x;
613           negate = true;
614       }
615 
616       double result;
617       if (x >= 0.5) {
618           double hiPrec[] = new double[2];
619           // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
620           exp(x*2.0, 0.0, hiPrec);
621 
622           double ya = hiPrec[0] + hiPrec[1];
623           double yb = -(ya - hiPrec[0] - hiPrec[1]);
624 
625           /* Numerator */
626           double na = -1.0 + ya;
627           double nb = -(na + 1.0 - ya);
628           double temp = na + yb;
629           nb += -(temp - na - yb);
630           na = temp;
631 
632           /* Denominator */
633           double da = 1.0 + ya;
634           double db = -(da - 1.0 - ya);
635           temp = da + yb;
636           db += -(temp - da - yb);
637           da = temp;
638 
639           temp = da * HEX_40000000;
640           double daa = da + temp - temp;
641           double dab = da - daa;
642 
643           // ratio = na/da
644           double ratio = na/da;
645           temp = ratio * HEX_40000000;
646           double ratioa = ratio + temp - temp;
647           double ratiob = ratio - ratioa;
648 
649           // Correct for rounding in division
650           ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
651 
652           // Account for nb
653           ratiob += nb / da;
654           // Account for db
655           ratiob += -db * na / da / da;
656 
657           result = ratioa + ratiob;
658       }
659       else {
660           double hiPrec[] = new double[2];
661           // tanh(x) = expm1(2x) / (expm1(2x) + 2)
662           expm1(x*2.0, hiPrec);
663 
664           double ya = hiPrec[0] + hiPrec[1];
665           double yb = -(ya - hiPrec[0] - hiPrec[1]);
666 
667           /* Numerator */
668           double na = ya;
669           double nb = yb;
670 
671           /* Denominator */
672           double da = 2.0 + ya;
673           double db = -(da - 2.0 - ya);
674           double temp = da + yb;
675           db += -(temp - da - yb);
676           da = temp;
677 
678           temp = da * HEX_40000000;
679           double daa = da + temp - temp;
680           double dab = da - daa;
681 
682           // ratio = na/da
683           double ratio = na/da;
684           temp = ratio * HEX_40000000;
685           double ratioa = ratio + temp - temp;
686           double ratiob = ratio - ratioa;
687 
688           // Correct for rounding in division
689           ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
690 
691           // Account for nb
692           ratiob += nb / da;
693           // Account for db
694           ratiob += -db * na / da / da;
695 
696           result = ratioa + ratiob;
697       }
698 
699       if (negate) {
700           result = -result;
701       }
702 
703       return result;
704     }
705 
706     /** Compute the inverse hyperbolic cosine of a number.
707      * @param a number on which evaluation is done
708      * @return inverse hyperbolic cosine of a
709      */
710     public static double acosh(final double a) {
711         return FastMath.log(a + FastMath.sqrt(a * a - 1));
712     }
713 
714     /** Compute the inverse hyperbolic sine of a number.
715      * @param a number on which evaluation is done
716      * @return inverse hyperbolic sine of a
717      */
718     public static double asinh(double a) {
719         boolean negative = false;
720         if (a < 0) {
721             negative = true;
722             a = -a;
723         }
724 
725         double absAsinh;
726         if (a > 0.167) {
727             absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
728         } else {
729             final double a2 = a * a;
730             if (a > 0.097) {
731                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * (F_1_9 - a2 * (F_1_11 - a2 * (F_1_13 - a2 * (F_1_15 - a2 * F_1_17 * F_15_16) * F_13_14) * F_11_12) * F_9_10) * F_7_8) * F_5_6) * F_3_4) * F_1_2);
732             } else if (a > 0.036) {
733                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * (F_1_9 - a2 * (F_1_11 - a2 * F_1_13 * F_11_12) * F_9_10) * F_7_8) * F_5_6) * F_3_4) * F_1_2);
734             } else if (a > 0.0036) {
735                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * (F_1_5 - a2 * (F_1_7 - a2 * F_1_9 * F_7_8) * F_5_6) * F_3_4) * F_1_2);
736             } else {
737                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
738             }
739         }
740 
741         return negative ? -absAsinh : absAsinh;
742     }
743 
744     /** Compute the inverse hyperbolic tangent of a number.
745      * @param a number on which evaluation is done
746      * @return inverse hyperbolic tangent of a
747      */
748     public static double atanh(double a) {
749         boolean negative = false;
750         if (a < 0) {
751             negative = true;
752             a = -a;
753         }
754 
755         double absAtanh;
756         if (a > 0.15) {
757             absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
758         } else {
759             final double a2 = a * a;
760             if (a > 0.087) {
761                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * (F_1_9 + a2 * (F_1_11 + a2 * (F_1_13 + a2 * (F_1_15 + a2 * F_1_17))))))));
762             } else if (a > 0.031) {
763                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * (F_1_9 + a2 * (F_1_11 + a2 * F_1_13))))));
764             } else if (a > 0.003) {
765                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
766             } else {
767                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
768             }
769         }
770 
771         return negative ? -absAtanh : absAtanh;
772     }
773 
774     /** Compute the signum of a number.
775      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
776      * @param a number on which evaluation is done
777      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
778      */
779     public static double signum(final double a) {
780         return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
781     }
782 
783     /** Compute the signum of a number.
784      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
785      * @param a number on which evaluation is done
786      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
787      */
788     public static float signum(final float a) {
789         return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
790     }
791 
792     /** Compute next number towards positive infinity.
793      * @param a number to which neighbor should be computed
794      * @return neighbor of a towards positive infinity
795      */
796     public static double nextUp(final double a) {
797         return nextAfter(a, Double.POSITIVE_INFINITY);
798     }
799 
800     /** Compute next number towards positive infinity.
801      * @param a number to which neighbor should be computed
802      * @return neighbor of a towards positive infinity
803      */
804     public static float nextUp(final float a) {
805         return nextAfter(a, Float.POSITIVE_INFINITY);
806     }
807 
808     /** Returns a pseudo-random number between 0.0 and 1.0.
809      * <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
810      * @return a random number between 0.0 and 1.0
811      */
812     public static double random() {
813         return Math.random();
814     }
815 
816     /**
817      * Exponential function.
818      *
819      * Computes exp(x), function result is nearly rounded.   It will be correctly
820      * rounded to the theoretical value for 99.9% of input values, otherwise it will
821      * have a 1 UPL error.
822      *
823      * Method:
824      *    Lookup intVal = exp(int(x))
825      *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
826      *    Compute z as the exponential of the remaining bits by a polynomial minus one
827      *    exp(x) = intVal * fracVal * (1 + z)
828      *
829      * Accuracy:
830      *    Calculation is done with 63 bits of precision, so result should be correctly
831      *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
832      *
833      * @param x   a double
834      * @return double e<sup>x</sup>
835      */
836     public static double exp(double x) {
837         return exp(x, 0.0, null);
838     }
839 
840     /**
841      * Internal helper method for exponential function.
842      * @param x original argument of the exponential function
843      * @param extra extra bits of precision on input (To Be Confirmed)
844      * @param hiPrec extra bits of precision on output (To Be Confirmed)
845      * @return exp(x)
846      */
847     private static double exp(double x, double extra, double[] hiPrec) {
848         double intPartA;
849         double intPartB;
850         int intVal;
851 
852         /* Lookup exp(floor(x)).
853          * intPartA will have the upper 22 bits, intPartB will have the lower
854          * 52 bits.
855          */
856         if (x < 0.0) {
857             intVal = (int) -x;
858 
859             if (intVal > 746) {
860                 if (hiPrec != null) {
861                     hiPrec[0] = 0.0;
862                     hiPrec[1] = 0.0;
863                 }
864                 return 0.0;
865             }
866 
867             if (intVal > 709) {
868                 /* This will produce a subnormal output */
869                 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
870                 if (hiPrec != null) {
871                     hiPrec[0] /= 285040095144011776.0;
872                     hiPrec[1] /= 285040095144011776.0;
873                 }
874                 return result;
875             }
876 
877             if (intVal == 709) {
878                 /* exp(1.494140625) is nearly a machine number... */
879                 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
880                 if (hiPrec != null) {
881                     hiPrec[0] /= 4.455505956692756620;
882                     hiPrec[1] /= 4.455505956692756620;
883                 }
884                 return result;
885             }
886 
887             intVal++;
888 
889             intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX-intVal];
890             intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX-intVal];
891 
892             intVal = -intVal;
893         } else {
894             intVal = (int) x;
895 
896             if (intVal > 709) {
897                 if (hiPrec != null) {
898                     hiPrec[0] = Double.POSITIVE_INFINITY;
899                     hiPrec[1] = 0.0;
900                 }
901                 return Double.POSITIVE_INFINITY;
902             }
903 
904             intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX+intVal];
905             intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX+intVal];
906         }
907 
908         /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
909          * x and look up the exp function of it.
910          * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
911          */
912         final int intFrac = (int) ((x - intVal) * 1024.0);
913         final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
914         final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
915 
916         /* epsilon is the difference in x from the nearest multiple of 2^-10.  It
917          * has a value in the range 0 <= epsilon < 2^-10.
918          * Do the subtraction from x as the last step to avoid possible loss of percison.
919          */
920         final double epsilon = x - (intVal + intFrac / 1024.0);
921 
922         /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial.  z has
923        full double precision (52 bits).  Since z < 2^-10, we will have
924        62 bits of precision when combined with the contant 1.  This will be
925        used in the last addition below to get proper rounding. */
926 
927         /* Remez generated polynomial.  Converges on the interval [0, 2^-10], error
928        is less than 0.5 ULP */
929         double z = 0.04168701738764507;
930         z = z * epsilon + 0.1666666505023083;
931         z = z * epsilon + 0.5000000000042687;
932         z = z * epsilon + 1.0;
933         z = z * epsilon + -3.940510424527919E-20;
934 
935         /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
936        expansion.
937        tempA is exact since intPartA and intPartB only have 22 bits each.
938        tempB will have 52 bits of precision.
939          */
940         double tempA = intPartA * fracPartA;
941         double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
942 
943         /* Compute the result.  (1+z)(tempA+tempB).  Order of operations is
944        important.  For accuracy add by increasing size.  tempA is exact and
945        much larger than the others.  If there are extra bits specified from the
946        pow() function, use them. */
947         final double tempC = tempB + tempA;
948         final double result;
949         if (extra != 0.0) {
950             result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
951         } else {
952             result = tempC*z + tempB + tempA;
953         }
954 
955         if (hiPrec != null) {
956             // If requesting high precision
957             hiPrec[0] = tempA;
958             hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
959         }
960 
961         return result;
962     }
963 
964     /** Compute exp(x) - 1
965      * @param x number to compute shifted exponential
966      * @return exp(x) - 1
967      */
968     public static double expm1(double x) {
969       return expm1(x, null);
970     }
971 
972     /** Internal helper method for expm1
973      * @param x number to compute shifted exponential
974      * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
975      * @return exp(x) - 1
976      */
977     private static double expm1(double x, double hiPrecOut[]) {
978         if (x != x || x == 0.0) { // NaN or zero
979             return x;
980         }
981 
982         if (x <= -1.0 || x >= 1.0) {
983             // If not between +/- 1.0
984             //return exp(x) - 1.0;
985             double hiPrec[] = new double[2];
986             exp(x, 0.0, hiPrec);
987             if (x > 0.0) {
988                 return -1.0 + hiPrec[0] + hiPrec[1];
989             } else {
990                 final double ra = -1.0 + hiPrec[0];
991                 double rb = -(ra + 1.0 - hiPrec[0]);
992                 rb += hiPrec[1];
993                 return ra + rb;
994             }
995         }
996 
997         double baseA;
998         double baseB;
999         double epsilon;
1000         boolean negative = false;
1001 
1002         if (x < 0.0) {
1003             x = -x;
1004             negative = true;
1005         }
1006 
1007         {
1008             int intFrac = (int) (x * 1024.0);
1009             double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
1010             double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1011 
1012             double temp = tempA + tempB;
1013             tempB = -(temp - tempA - tempB);
1014             tempA = temp;
1015 
1016             temp = tempA * HEX_40000000;
1017             baseA = tempA + temp - temp;
1018             baseB = tempB + (tempA - baseA);
1019 
1020             epsilon = x - intFrac/1024.0;
1021         }
1022 
1023 
1024         /* Compute expm1(epsilon) */
1025         double zb = 0.008336750013465571;
1026         zb = zb * epsilon + 0.041666663879186654;
1027         zb = zb * epsilon + 0.16666666666745392;
1028         zb = zb * epsilon + 0.49999999999999994;
1029         zb *= epsilon;
1030         zb *= epsilon;
1031 
1032         double za = epsilon;
1033         double temp = za + zb;
1034         zb = -(temp - za - zb);
1035         za = temp;
1036 
1037         temp = za * HEX_40000000;
1038         temp = za + temp - temp;
1039         zb += za - temp;
1040         za = temp;
1041 
1042         /* Combine the parts.   expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
1043         double ya = za * baseA;
1044         //double yb = za*baseB + zb*baseA + zb*baseB;
1045         temp = ya + za * baseB;
1046         double yb = -(temp - ya - za * baseB);
1047         ya = temp;
1048 
1049         temp = ya + zb * baseA;
1050         yb += -(temp - ya - zb * baseA);
1051         ya = temp;
1052 
1053         temp = ya + zb * baseB;
1054         yb += -(temp - ya - zb*baseB);
1055         ya = temp;
1056 
1057         //ya = ya + za + baseA;
1058         //yb = yb + zb + baseB;
1059         temp = ya + baseA;
1060         yb += -(temp - baseA - ya);
1061         ya = temp;
1062 
1063         temp = ya + za;
1064         //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
1065         yb += -(temp - ya - za);
1066         ya = temp;
1067 
1068         temp = ya + baseB;
1069         //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
1070         yb += -(temp - ya - baseB);
1071         ya = temp;
1072 
1073         temp = ya + zb;
1074         //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
1075         yb += -(temp - ya - zb);
1076         ya = temp;
1077 
1078         if (negative) {
1079             /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
1080             double denom = 1.0 + ya;
1081             double denomr = 1.0 / denom;
1082             double denomb = -(denom - 1.0 - ya) + yb;
1083             double ratio = ya * denomr;
1084             temp = ratio * HEX_40000000;
1085             final double ra = ratio + temp - temp;
1086             double rb = ratio - ra;
1087 
1088             temp = denom * HEX_40000000;
1089             za = denom + temp - temp;
1090             zb = denom - za;
1091 
1092             rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
1093 
1094             // f(x) = x/1+x
1095             // Compute f'(x)
1096             // Product rule:  d(uv) = du*v + u*dv
1097             // Chain rule:  d(f(g(x)) = f'(g(x))*f(g'(x))
1098             // d(1/x) = -1/(x*x)
1099             // d(1/1+x) = -1/( (1+x)^2) *  1 =  -1/((1+x)*(1+x))
1100             // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))
1101 
1102             // Adjust for yb
1103             rb += yb * denomr;                      // numerator
1104             rb += -ya * denomb * denomr * denomr;   // denominator
1105 
1106             // negate
1107             ya = -ra;
1108             yb = -rb;
1109         }
1110 
1111         if (hiPrecOut != null) {
1112             hiPrecOut[0] = ya;
1113             hiPrecOut[1] = yb;
1114         }
1115 
1116         return ya + yb;
1117     }
1118 
1119     /**
1120      * Natural logarithm.
1121      *
1122      * @param x   a double
1123      * @return log(x)
1124      */
1125     public static double log(final double x) {
1126         return log(x, null);
1127     }
1128 
1129     /**
1130      * Internal helper method for natural logarithm function.
1131      * @param x original argument of the natural logarithm function
1132      * @param hiPrec extra bits of precision on output (To Be Confirmed)
1133      * @return log(x)
1134      */
1135     private static double log(final double x, final double[] hiPrec) {
1136         if (x==0) { // Handle special case of +0/-0
1137             return Double.NEGATIVE_INFINITY;
1138         }
1139         long bits = Double.doubleToRawLongBits(x);
1140 
1141         /* Handle special cases of negative input, and NaN */
1142         if (((bits & 0x8000000000000000L) != 0 || x != x) && x != 0.0) {
1143             if (hiPrec != null) {
1144                 hiPrec[0] = Double.NaN;
1145             }
1146 
1147             return Double.NaN;
1148         }
1149 
1150         /* Handle special cases of Positive infinity. */
1151         if (x == Double.POSITIVE_INFINITY) {
1152             if (hiPrec != null) {
1153                 hiPrec[0] = Double.POSITIVE_INFINITY;
1154             }
1155 
1156             return Double.POSITIVE_INFINITY;
1157         }
1158 
1159         /* Extract the exponent */
1160         int exp = (int)(bits >> 52)-1023;
1161 
1162         if ((bits & 0x7ff0000000000000L) == 0) {
1163             // Subnormal!
1164             if (x == 0) {
1165                 // Zero
1166                 if (hiPrec != null) {
1167                     hiPrec[0] = Double.NEGATIVE_INFINITY;
1168                 }
1169 
1170                 return Double.NEGATIVE_INFINITY;
1171             }
1172 
1173             /* Normalize the subnormal number. */
1174             bits <<= 1;
1175             while ( (bits & 0x0010000000000000L) == 0) {
1176                 --exp;
1177                 bits <<= 1;
1178             }
1179         }
1180 
1181 
1182         if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
1183             /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
1184            polynomial expansion in higer precision. */
1185 
1186             /* Compute x - 1.0 and split it */
1187             double xa = x - 1.0;
1188             double xb = xa - x + 1.0;
1189             double tmp = xa * HEX_40000000;
1190             double aa = xa + tmp - tmp;
1191             double ab = xa - aa;
1192             xa = aa;
1193             xb = ab;
1194 
1195             final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
1196             double ya = lnCoef_last[0];
1197             double yb = lnCoef_last[1];
1198 
1199             for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
1200                 /* Multiply a = y * x */
1201                 aa = ya * xa;
1202                 ab = ya * xb + yb * xa + yb * xb;
1203                 /* split, so now y = a */
1204                 tmp = aa * HEX_40000000;
1205                 ya = aa + tmp - tmp;
1206                 yb = aa - ya + ab;
1207 
1208                 /* Add  a = y + lnQuickCoef */
1209                 final double[] lnCoef_i = LN_QUICK_COEF[i];
1210                 aa = ya + lnCoef_i[0];
1211                 ab = yb + lnCoef_i[1];
1212                 /* Split y = a */
1213                 tmp = aa * HEX_40000000;
1214                 ya = aa + tmp - tmp;
1215                 yb = aa - ya + ab;
1216             }
1217 
1218             /* Multiply a = y * x */
1219             aa = ya * xa;
1220             ab = ya * xb + yb * xa + yb * xb;
1221             /* split, so now y = a */
1222             tmp = aa * HEX_40000000;
1223             ya = aa + tmp - tmp;
1224             yb = aa - ya + ab;
1225 
1226             return ya + yb;
1227         }
1228 
1229         // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
1230         final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1231 
1232         /*
1233     double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1234 
1235     epsilon -= 1.0;
1236          */
1237 
1238         // y is the most significant 10 bits of the mantissa
1239         //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1240         //double epsilon = (x - y) / y;
1241         final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
1242 
1243         double lnza = 0.0;
1244         double lnzb = 0.0;
1245 
1246         if (hiPrec != null) {
1247             /* split epsilon -> x */
1248             double tmp = epsilon * HEX_40000000;
1249             double aa = epsilon + tmp - tmp;
1250             double ab = epsilon - aa;
1251             double xa = aa;
1252             double xb = ab;
1253 
1254             /* Need a more accurate epsilon, so adjust the division. */
1255             final double numer = bits & 0x3ffffffffffL;
1256             final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
1257             aa = numer - xa*denom - xb * denom;
1258             xb += aa / denom;
1259 
1260             /* Remez polynomial evaluation */
1261             final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
1262             double ya = lnCoef_last[0];
1263             double yb = lnCoef_last[1];
1264 
1265             for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
1266                 /* Multiply a = y * x */
1267                 aa = ya * xa;
1268                 ab = ya * xb + yb * xa + yb * xb;
1269                 /* split, so now y = a */
1270                 tmp = aa * HEX_40000000;
1271                 ya = aa + tmp - tmp;
1272                 yb = aa - ya + ab;
1273 
1274                 /* Add  a = y + lnHiPrecCoef */
1275                 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1276                 aa = ya + lnCoef_i[0];
1277                 ab = yb + lnCoef_i[1];
1278                 /* Split y = a */
1279                 tmp = aa * HEX_40000000;
1280                 ya = aa + tmp - tmp;
1281                 yb = aa - ya + ab;
1282             }
1283 
1284             /* Multiply a = y * x */
1285             aa = ya * xa;
1286             ab = ya * xb + yb * xa + yb * xb;
1287 
1288             /* split, so now lnz = a */
1289             /*
1290       tmp = aa * 1073741824.0;
1291       lnza = aa + tmp - tmp;
1292       lnzb = aa - lnza + ab;
1293              */
1294             lnza = aa + ab;
1295             lnzb = -(lnza - aa - ab);
1296         } else {
1297             /* High precision not required.  Eval Remez polynomial
1298          using standard double precision */
1299             lnza = -0.16624882440418567;
1300             lnza = lnza * epsilon + 0.19999954120254515;
1301             lnza = lnza * epsilon + -0.2499999997677497;
1302             lnza = lnza * epsilon + 0.3333333333332802;
1303             lnza = lnza * epsilon + -0.5;
1304             lnza = lnza * epsilon + 1.0;
1305             lnza *= epsilon;
1306         }
1307 
1308         /* Relative sizes:
1309          * lnzb     [0, 2.33E-10]
1310          * lnm[1]   [0, 1.17E-7]
1311          * ln2B*exp [0, 1.12E-4]
1312          * lnza      [0, 9.7E-4]
1313          * lnm[0]   [0, 0.692]
1314          * ln2A*exp [0, 709]
1315          */
1316 
1317         /* Compute the following sum:
1318          * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1319          */
1320 
1321         //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1322         double a = LN_2_A*exp;
1323         double b = 0.0;
1324         double c = a+lnm[0];
1325         double d = -(c-a-lnm[0]);
1326         a = c;
1327         b += d;
1328 
1329         c = a + lnza;
1330         d = -(c - a - lnza);
1331         a = c;
1332         b += d;
1333 
1334         c = a + LN_2_B*exp;
1335         d = -(c - a - LN_2_B*exp);
1336         a = c;
1337         b += d;
1338 
1339         c = a + lnm[1];
1340         d = -(c - a - lnm[1]);
1341         a = c;
1342         b += d;
1343 
1344         c = a + lnzb;
1345         d = -(c - a - lnzb);
1346         a = c;
1347         b += d;
1348 
1349         if (hiPrec != null) {
1350             hiPrec[0] = a;
1351             hiPrec[1] = b;
1352         }
1353 
1354         return a + b;
1355     }
1356 
1357     /**
1358      * Computes log(1 + x).
1359      *
1360      * @param x Number.
1361      * @return {@code log(1 + x)}.
1362      */
1363     public static double log1p(final double x) {
1364         if (x == -1) {
1365             return Double.NEGATIVE_INFINITY;
1366         }
1367 
1368         if (x == Double.POSITIVE_INFINITY) {
1369             return Double.POSITIVE_INFINITY;
1370         }
1371 
1372         if (x > 1e-6 ||
1373             x < -1e-6) {
1374             final double xpa = 1 + x;
1375             final double xpb = -(xpa - 1 - x);
1376 
1377             final double[] hiPrec = new double[2];
1378             final double lores = log(xpa, hiPrec);
1379             if (Double.isInfinite(lores)) { // Don't allow this to be converted to NaN
1380                 return lores;
1381             }
1382 
1383             // Do a taylor series expansion around xpa:
1384             //   f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
1385             final double fx1 = xpb / xpa;
1386             final double epsilon = 0.5 * fx1 + 1;
1387             return epsilon * fx1 + hiPrec[1] + hiPrec[0];
1388         } else {
1389             // Value is small |x| < 1e6, do a Taylor series centered on 1.
1390             final double y = (x * F_1_3 - F_1_2) * x + 1;
1391             return y * x;
1392         }
1393     }
1394 
1395     /** Compute the base 10 logarithm.
1396      * @param x a number
1397      * @return log10(x)
1398      */
1399     public static double log10(final double x) {
1400         final double hiPrec[] = new double[2];
1401 
1402         final double lores = log(x, hiPrec);
1403         if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1404             return lores;
1405         }
1406 
1407         final double tmp = hiPrec[0] * HEX_40000000;
1408         final double lna = hiPrec[0] + tmp - tmp;
1409         final double lnb = hiPrec[0] - lna + hiPrec[1];
1410 
1411         final double rln10a = 0.4342944622039795;
1412         final double rln10b = 1.9699272335463627E-8;
1413 
1414         return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
1415     }
1416 
1417     /**
1418      * Computes the <a href="http://mathworld.wolfram.com/Logarithm.html">
1419      * logarithm</a> in a given base.
1420      *
1421      * Returns {@code NaN} if either argument is negative.
1422      * If {@code base} is 0 and {@code x} is positive, 0 is returned.
1423      * If {@code base} is positive and {@code x} is 0,
1424      * {@code Double.NEGATIVE_INFINITY} is returned.
1425      * If both arguments are 0, the result is {@code NaN}.
1426      *
1427      * @param base Base of the logarithm, must be greater than 0.
1428      * @param x Argument, must be greater than 0.
1429      * @return the value of the logarithm, i.e. the number {@code y} such that
1430      * <code>base<sup>y</sup> = x</code>.
1431      * @since 1.2 (previously in {@code MathUtils}, moved as of version 3.0)
1432      */
1433     public static double log(double base, double x) {
1434         return log(x) / log(base);
1435     }
1436 
1437     /**
1438      * Power function.  Compute x^y.
1439      *
1440      * @param x   a double
1441      * @param y   a double
1442      * @return double
1443      */
1444     public static double pow(double x, double y) {
1445         final double lns[] = new double[2];
1446 
1447         if (y == 0.0) {
1448             return 1.0;
1449         }
1450 
1451         if (x != x) { // X is NaN
1452             return x;
1453         }
1454 
1455 
1456         if (x == 0) {
1457             long bits = Double.doubleToRawLongBits(x);
1458             if ((bits & 0x8000000000000000L) != 0) {
1459                 // -zero
1460                 long yi = (long) y;
1461 
1462                 if (y < 0 && y == yi && (yi & 1) == 1) {
1463                     return Double.NEGATIVE_INFINITY;
1464                 }
1465 
1466                 if (y > 0 && y == yi && (yi & 1) == 1) {
1467                     return -0.0;
1468                 }
1469             }
1470 
1471             if (y < 0) {
1472                 return Double.POSITIVE_INFINITY;
1473             }
1474             if (y > 0) {
1475                 return 0.0;
1476             }
1477 
1478             return Double.NaN;
1479         }
1480 
1481         if (x == Double.POSITIVE_INFINITY) {
1482             if (y != y) { // y is NaN
1483                 return y;
1484             }
1485             if (y < 0.0) {
1486                 return 0.0;
1487             } else {
1488                 return Double.POSITIVE_INFINITY;
1489             }
1490         }
1491 
1492         if (y == Double.POSITIVE_INFINITY) {
1493             if (x * x == 1.0) {
1494                 return Double.NaN;
1495             }
1496 
1497             if (x * x > 1.0) {
1498                 return Double.POSITIVE_INFINITY;
1499             } else {
1500                 return 0.0;
1501             }
1502         }
1503 
1504         if (x == Double.NEGATIVE_INFINITY) {
1505             if (y != y) { // y is NaN
1506                 return y;
1507             }
1508 
1509             if (y < 0) {
1510                 long yi = (long) y;
1511                 if (y == yi && (yi & 1) == 1) {
1512                     return -0.0;
1513                 }
1514 
1515                 return 0.0;
1516             }
1517 
1518             if (y > 0)  {
1519                 long yi = (long) y;
1520                 if (y == yi && (yi & 1) == 1) {
1521                     return Double.NEGATIVE_INFINITY;
1522                 }
1523 
1524                 return Double.POSITIVE_INFINITY;
1525             }
1526         }
1527 
1528         if (y == Double.NEGATIVE_INFINITY) {
1529 
1530             if (x * x == 1.0) {
1531                 return Double.NaN;
1532             }
1533 
1534             if (x * x < 1.0) {
1535                 return Double.POSITIVE_INFINITY;
1536             } else {
1537                 return 0.0;
1538             }
1539         }
1540 
1541         /* Handle special case x<0 */
1542         if (x < 0) {
1543             // y is an even integer in this case
1544             if (y >= TWO_POWER_53 || y <= -TWO_POWER_53) {
1545                 return pow(-x, y);
1546             }
1547 
1548             if (y == (long) y) {
1549                 // If y is an integer
1550                 return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
1551             } else {
1552                 return Double.NaN;
1553             }
1554         }
1555 
1556         /* Split y into ya and yb such that y = ya+yb */
1557         double ya;
1558         double yb;
1559         if (y < 8e298 && y > -8e298) {
1560             double tmp1 = y * HEX_40000000;
1561             ya = y + tmp1 - tmp1;
1562             yb = y - ya;
1563         } else {
1564             double tmp1 = y * 9.31322574615478515625E-10;
1565             double tmp2 = tmp1 * 9.31322574615478515625E-10;
1566             ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
1567             yb = y - ya;
1568         }
1569 
1570         /* Compute ln(x) */
1571         final double lores = log(x, lns);
1572         if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1573             return lores;
1574         }
1575 
1576         double lna = lns[0];
1577         double lnb = lns[1];
1578 
1579         /* resplit lns */
1580         double tmp1 = lna * HEX_40000000;
1581         double tmp2 = lna + tmp1 - tmp1;
1582         lnb += lna - tmp2;
1583         lna = tmp2;
1584 
1585         // y*ln(x) = (aa+ab)
1586         final double aa = lna * ya;
1587         final double ab = lna * yb + lnb * ya + lnb * yb;
1588 
1589         lna = aa+ab;
1590         lnb = -(lna - aa - ab);
1591 
1592         double z = 1.0 / 120.0;
1593         z = z * lnb + (1.0 / 24.0);
1594         z = z * lnb + (1.0 / 6.0);
1595         z = z * lnb + 0.5;
1596         z = z * lnb + 1.0;
1597         z *= lnb;
1598 
1599         final double result = exp(lna, z, null);
1600         //result = result + result * z;
1601         return result;
1602     }
1603 
1604 
1605     /**
1606      * Raise a double to an int power.
1607      *
1608      * @param d Number to raise.
1609      * @param e Exponent.
1610      * @return d<sup>e</sup>
1611      * @since 3.1
1612      */
1613     public static double pow(double d, int e) {
1614 
1615         if (e == 0) {
1616             return 1.0;
1617         } else if (e < 0) {
1618             e = -e;
1619             d = 1.0 / d;
1620         }
1621 
1622         // split d as two 26 bits numbers
1623         // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1624         final int splitFactor = 0x8000001;
1625         final double cd       = splitFactor * d;
1626         final double d1High   = cd - (cd - d);
1627         final double d1Low    = d - d1High;
1628 
1629         // prepare result
1630         double resultHigh = 1;
1631         double resultLow  = 0;
1632 
1633         // d^(2p)
1634         double d2p     = d;
1635         double d2pHigh = d1High;
1636         double d2pLow  = d1Low;
1637 
1638         while (e != 0) {
1639 
1640             if ((e & 0x1) != 0) {
1641                 // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
1642                 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1643                 final double tmpHigh = resultHigh * d2p;
1644                 final double cRH     = splitFactor * resultHigh;
1645                 final double rHH     = cRH - (cRH - resultHigh);
1646                 final double rHL     = resultHigh - rHH;
1647                 final double tmpLow  = rHL * d2pLow - (((tmpHigh - rHH * d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
1648                 resultHigh = tmpHigh;
1649                 resultLow  = resultLow * d2p + tmpLow;
1650             }
1651 
1652             // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
1653             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1654             final double tmpHigh = d2pHigh * d2p;
1655             final double cD2pH   = splitFactor * d2pHigh;
1656             final double d2pHH   = cD2pH - (cD2pH - d2pHigh);
1657             final double d2pHL   = d2pHigh - d2pHH;
1658             final double tmpLow  = d2pHL * d2pLow - (((tmpHigh - d2pHH * d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
1659             final double cTmpH   = splitFactor * tmpHigh;
1660             d2pHigh = cTmpH - (cTmpH - tmpHigh);
1661             d2pLow  = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
1662             d2p     = d2pHigh + d2pLow;
1663 
1664             e >>= 1;
1665 
1666         }
1667 
1668         return resultHigh + resultLow;
1669 
1670     }
1671 
1672     /**
1673      *  Computes sin(x) - x, where |x| < 1/16.
1674      *  Use a Remez polynomial approximation.
1675      *  @param x a number smaller than 1/16
1676      *  @return sin(x) - x
1677      */
1678     private static double polySine(final double x)
1679     {
1680         double x2 = x*x;
1681 
1682         double p = 2.7553817452272217E-6;
1683         p = p * x2 + -1.9841269659586505E-4;
1684         p = p * x2 + 0.008333333333329196;
1685         p = p * x2 + -0.16666666666666666;
1686         //p *= x2;
1687         //p *= x;
1688         p = p * x2 * x;
1689 
1690         return p;
1691     }
1692 
1693     /**
1694      *  Computes cos(x) - 1, where |x| < 1/16.
1695      *  Use a Remez polynomial approximation.
1696      *  @param x a number smaller than 1/16
1697      *  @return cos(x) - 1
1698      */
1699     private static double polyCosine(double x) {
1700         double x2 = x*x;
1701 
1702         double p = 2.479773539153719E-5;
1703         p = p * x2 + -0.0013888888689039883;
1704         p = p * x2 + 0.041666666666621166;
1705         p = p * x2 + -0.49999999999999994;
1706         p *= x2;
1707 
1708         return p;
1709     }
1710 
1711     /**
1712      *  Compute sine over the first quadrant (0 < x < pi/2).
1713      *  Use combination of table lookup and rational polynomial expansion.
1714      *  @param xa number from which sine is requested
1715      *  @param xb extra bits for x (may be 0.0)
1716      *  @return sin(xa + xb)
1717      */
1718     private static double sinQ(double xa, double xb) {
1719         int idx = (int) ((xa * 8.0) + 0.5);
1720         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
1721 
1722         // Table lookups
1723         final double sintA = SINE_TABLE_A[idx];
1724         final double sintB = SINE_TABLE_B[idx];
1725         final double costA = COSINE_TABLE_A[idx];
1726         final double costB = COSINE_TABLE_B[idx];
1727 
1728         // Polynomial eval of sin(epsilon), cos(epsilon)
1729         double sinEpsA = epsilon;
1730         double sinEpsB = polySine(epsilon);
1731         final double cosEpsA = 1.0;
1732         final double cosEpsB = polyCosine(epsilon);
1733 
1734         // Split epsilon   xa + xb = x
1735         final double temp = sinEpsA * HEX_40000000;
1736         double temp2 = (sinEpsA + temp) - temp;
1737         sinEpsB +=  sinEpsA - temp2;
1738         sinEpsA = temp2;
1739 
1740         /* Compute sin(x) by angle addition formula */
1741         double result;
1742 
1743         /* Compute the following sum:
1744          *
1745          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1746          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1747          *
1748          * Ranges of elements
1749          *
1750          * xxxtA   0            PI/2
1751          * xxxtB   -1.5e-9      1.5e-9
1752          * sinEpsA -0.0625      0.0625
1753          * sinEpsB -6e-11       6e-11
1754          * cosEpsA  1.0
1755          * cosEpsB  0           -0.0625
1756          *
1757          */
1758 
1759         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1760         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1761 
1762         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1763         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1764         double a = 0;
1765         double b = 0;
1766 
1767         double t = sintA;
1768         double c = a + t;
1769         double d = -(c - a - t);
1770         a = c;
1771         b += d;
1772 
1773         t = costA * sinEpsA;
1774         c = a + t;
1775         d = -(c - a - t);
1776         a = c;
1777         b += d;
1778 
1779         b = b + sintA * cosEpsB + costA * sinEpsB;
1780         /*
1781     t = sintA*cosEpsB;
1782     c = a + t;
1783     d = -(c - a - t);
1784     a = c;
1785     b = b + d;
1786 
1787     t = costA*sinEpsB;
1788     c = a + t;
1789     d = -(c - a - t);
1790     a = c;
1791     b = b + d;
1792          */
1793 
1794         b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
1795         /*
1796     t = sintB;
1797     c = a + t;
1798     d = -(c - a - t);
1799     a = c;
1800     b = b + d;
1801 
1802     t = costB*sinEpsA;
1803     c = a + t;
1804     d = -(c - a - t);
1805     a = c;
1806     b = b + d;
1807 
1808     t = sintB*cosEpsB;
1809     c = a + t;
1810     d = -(c - a - t);
1811     a = c;
1812     b = b + d;
1813 
1814     t = costB*sinEpsB;
1815     c = a + t;
1816     d = -(c - a - t);
1817     a = c;
1818     b = b + d;
1819          */
1820 
1821         if (xb != 0.0) {
1822             t = ((costA + costB) * (cosEpsA + cosEpsB) -
1823                  (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;  // approximate cosine*xb
1824             c = a + t;
1825             d = -(c - a - t);
1826             a = c;
1827             b += d;
1828         }
1829 
1830         result = a + b;
1831 
1832         return result;
1833     }
1834 
1835     /**
1836      * Compute cosine in the first quadrant by subtracting input from PI/2 and
1837      * then calling sinQ.  This is more accurate as the input approaches PI/2.
1838      *  @param xa number from which cosine is requested
1839      *  @param xb extra bits for x (may be 0.0)
1840      *  @return cos(xa + xb)
1841      */
1842     private static double cosQ(double xa, double xb) {
1843         final double pi2a = 1.5707963267948966;
1844         final double pi2b = 6.123233995736766E-17;
1845 
1846         final double a = pi2a - xa;
1847         double b = -(a - pi2a + xa);
1848         b += pi2b - xb;
1849 
1850         return sinQ(a, b);
1851     }
1852 
1853     /**
1854      *  Compute tangent (or cotangent) over the first quadrant.   0 < x < pi/2
1855      *  Use combination of table lookup and rational polynomial expansion.
1856      *  @param xa number from which sine is requested
1857      *  @param xb extra bits for x (may be 0.0)
1858      *  @param cotanFlag if true, compute the cotangent instead of the tangent
1859      *  @return tan(xa+xb) (or cotangent, depending on cotanFlag)
1860      */
1861     private static double tanQ(double xa, double xb, boolean cotanFlag) {
1862 
1863         int idx = (int) ((xa * 8.0) + 0.5);
1864         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
1865 
1866         // Table lookups
1867         final double sintA = SINE_TABLE_A[idx];
1868         final double sintB = SINE_TABLE_B[idx];
1869         final double costA = COSINE_TABLE_A[idx];
1870         final double costB = COSINE_TABLE_B[idx];
1871 
1872         // Polynomial eval of sin(epsilon), cos(epsilon)
1873         double sinEpsA = epsilon;
1874         double sinEpsB = polySine(epsilon);
1875         final double cosEpsA = 1.0;
1876         final double cosEpsB = polyCosine(epsilon);
1877 
1878         // Split epsilon   xa + xb = x
1879         double temp = sinEpsA * HEX_40000000;
1880         double temp2 = (sinEpsA + temp) - temp;
1881         sinEpsB +=  sinEpsA - temp2;
1882         sinEpsA = temp2;
1883 
1884         /* Compute sin(x) by angle addition formula */
1885 
1886         /* Compute the following sum:
1887          *
1888          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1889          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1890          *
1891          * Ranges of elements
1892          *
1893          * xxxtA   0            PI/2
1894          * xxxtB   -1.5e-9      1.5e-9
1895          * sinEpsA -0.0625      0.0625
1896          * sinEpsB -6e-11       6e-11
1897          * cosEpsA  1.0
1898          * cosEpsB  0           -0.0625
1899          *
1900          */
1901 
1902         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1903         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1904 
1905         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1906         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1907         double a = 0;
1908         double b = 0;
1909 
1910         // Compute sine
1911         double t = sintA;
1912         double c = a + t;
1913         double d = -(c - a - t);
1914         a = c;
1915         b += d;
1916 
1917         t = costA*sinEpsA;
1918         c = a + t;
1919         d = -(c - a - t);
1920         a = c;
1921         b += d;
1922 
1923         b += sintA*cosEpsB + costA*sinEpsB;
1924         b += sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1925 
1926         double sina = a + b;
1927         double sinb = -(sina - a - b);
1928 
1929         // Compute cosine
1930 
1931         a = b = c = d = 0.0;
1932 
1933         t = costA*cosEpsA;
1934         c = a + t;
1935         d = -(c - a - t);
1936         a = c;
1937         b += d;
1938 
1939         t = -sintA*sinEpsA;
1940         c = a + t;
1941         d = -(c - a - t);
1942         a = c;
1943         b += d;
1944 
1945         b += costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
1946         b -= sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB;
1947 
1948         double cosa = a + b;
1949         double cosb = -(cosa - a - b);
1950 
1951         if (cotanFlag) {
1952             double tmp;
1953             tmp = cosa; cosa = sina; sina = tmp;
1954             tmp = cosb; cosb = sinb; sinb = tmp;
1955         }
1956 
1957 
1958         /* estimate and correct, compute 1.0/(cosa+cosb) */
1959         /*
1960     double est = (sina+sinb)/(cosa+cosb);
1961     double err = (sina - cosa*est) + (sinb - cosb*est);
1962     est += err/(cosa+cosb);
1963     err = (sina - cosa*est) + (sinb - cosb*est);
1964          */
1965 
1966         // f(x) = 1/x,   f'(x) = -1/x^2
1967 
1968         double est = sina/cosa;
1969 
1970         /* Split the estimate to get more accurate read on division rounding */
1971         temp = est * HEX_40000000;
1972         double esta = (est + temp) - temp;
1973         double estb =  est - esta;
1974 
1975         temp = cosa * HEX_40000000;
1976         double cosaa = (cosa + temp) - temp;
1977         double cosab =  cosa - cosaa;
1978 
1979         //double err = (sina - est*cosa)/cosa;  // Correction for division rounding
1980         double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa;  // Correction for division rounding
1981         err += sinb/cosa;                     // Change in est due to sinb
1982         err += -sina * cosb / cosa / cosa;    // Change in est due to cosb
1983 
1984         if (xb != 0.0) {
1985             // tan' = 1 + tan^2      cot' = -(1 + cot^2)
1986             // Approximate impact of xb
1987             double xbadj = xb + est*est*xb;
1988             if (cotanFlag) {
1989                 xbadj = -xbadj;
1990             }
1991 
1992             err += xbadj;
1993         }
1994 
1995         return est+err;
1996     }
1997 
1998     /** Reduce the input argument using the Payne and Hanek method.
1999      *  This is good for all inputs 0.0 < x < inf
2000      *  Output is remainder after dividing by PI/2
2001      *  The result array should contain 3 numbers.
2002      *  result[0] is the integer portion, so mod 4 this gives the quadrant.
2003      *  result[1] is the upper bits of the remainder
2004      *  result[2] is the lower bits of the remainder
2005      *
2006      * @param x number to reduce
2007      * @param result placeholder where to put the result
2008      */
2009     private static void reducePayneHanek(double x, double result[])
2010     {
2011         /* Convert input double to bits */
2012         long inbits = Double.doubleToRawLongBits(x);
2013         int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2014 
2015         /* Convert to fixed point representation */
2016         inbits &= 0x000fffffffffffffL;
2017         inbits |= 0x0010000000000000L;
2018 
2019         /* Normalize input to be between 0.5 and 1.0 */
2020         exponent++;
2021         inbits <<= 11;
2022 
2023         /* Based on the exponent, get a shifted copy of recip2pi */
2024         long shpi0;
2025         long shpiA;
2026         long shpiB;
2027         int idx = exponent >> 6;
2028         int shift = exponent - (idx << 6);
2029 
2030         if (shift != 0) {
2031             shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
2032             shpi0 |= RECIP_2PI[idx] >>> (64-shift);
2033             shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
2034             shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
2035         } else {
2036             shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
2037             shpiA = RECIP_2PI[idx];
2038             shpiB = RECIP_2PI[idx+1];
2039         }
2040 
2041         /* Multiply input by shpiA */
2042         long a = inbits >>> 32;
2043         long b = inbits & 0xffffffffL;
2044 
2045         long c = shpiA >>> 32;
2046         long d = shpiA & 0xffffffffL;
2047 
2048         long ac = a * c;
2049         long bd = b * d;
2050         long bc = b * c;
2051         long ad = a * d;
2052 
2053         long prodB = bd + (ad << 32);
2054         long prodA = ac + (ad >>> 32);
2055 
2056         boolean bita = (bd & 0x8000000000000000L) != 0;
2057         boolean bitb = (ad & 0x80000000L ) != 0;
2058         boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2059 
2060         /* Carry */
2061         if ( (bita && bitb) ||
2062                 ((bita || bitb) && !bitsum) ) {
2063             prodA++;
2064         }
2065 
2066         bita = (prodB & 0x8000000000000000L) != 0;
2067         bitb = (bc & 0x80000000L ) != 0;
2068 
2069         prodB += bc << 32;
2070         prodA += bc >>> 32;
2071 
2072         bitsum = (prodB & 0x8000000000000000L) != 0;
2073 
2074         /* Carry */
2075         if ( (bita && bitb) ||
2076                 ((bita || bitb) && !bitsum) ) {
2077             prodA++;
2078         }
2079 
2080         /* Multiply input by shpiB */
2081         c = shpiB >>> 32;
2082         d = shpiB & 0xffffffffL;
2083         ac = a * c;
2084         bc = b * c;
2085         ad = a * d;
2086 
2087         /* Collect terms */
2088         ac += (bc + ad) >>> 32;
2089 
2090         bita = (prodB & 0x8000000000000000L) != 0;
2091         bitb = (ac & 0x8000000000000000L ) != 0;
2092         prodB += ac;
2093         bitsum = (prodB & 0x8000000000000000L) != 0;
2094         /* Carry */
2095         if ( (bita && bitb) ||
2096                 ((bita || bitb) && !bitsum) ) {
2097             prodA++;
2098         }
2099 
2100         /* Multiply by shpi0 */
2101         c = shpi0 >>> 32;
2102         d = shpi0 & 0xffffffffL;
2103 
2104         bd = b * d;
2105         bc = b * c;
2106         ad = a * d;
2107 
2108         prodA += bd + ((bc + ad) << 32);
2109 
2110         /*
2111          * prodA, prodB now contain the remainder as a fraction of PI.  We want this as a fraction of
2112          * PI/2, so use the following steps:
2113          * 1.) multiply by 4.
2114          * 2.) do a fixed point muliply by PI/4.
2115          * 3.) Convert to floating point.
2116          * 4.) Multiply by 2
2117          */
2118 
2119         /* This identifies the quadrant */
2120         int intPart = (int)(prodA >>> 62);
2121 
2122         /* Multiply by 4 */
2123         prodA <<= 2;
2124         prodA |= prodB >>> 62;
2125         prodB <<= 2;
2126 
2127         /* Multiply by PI/4 */
2128         a = prodA >>> 32;
2129         b = prodA & 0xffffffffL;
2130 
2131         c = PI_O_4_BITS[0] >>> 32;
2132         d = PI_O_4_BITS[0] & 0xffffffffL;
2133 
2134         ac = a * c;
2135         bd = b * d;
2136         bc = b * c;
2137         ad = a * d;
2138 
2139         long prod2B = bd + (ad << 32);
2140         long prod2A = ac + (ad >>> 32);
2141 
2142         bita = (bd & 0x8000000000000000L) != 0;
2143         bitb = (ad & 0x80000000L ) != 0;
2144         bitsum = (prod2B & 0x8000000000000000L) != 0;
2145 
2146         /* Carry */
2147         if ( (bita && bitb) ||
2148                 ((bita || bitb) && !bitsum) ) {
2149             prod2A++;
2150         }
2151 
2152         bita = (prod2B & 0x8000000000000000L) != 0;
2153         bitb = (bc & 0x80000000L ) != 0;
2154 
2155         prod2B += bc << 32;
2156         prod2A += bc >>> 32;
2157 
2158         bitsum = (prod2B & 0x8000000000000000L) != 0;
2159 
2160         /* Carry */
2161         if ( (bita && bitb) ||
2162                 ((bita || bitb) && !bitsum) ) {
2163             prod2A++;
2164         }
2165 
2166         /* Multiply input by pio4bits[1] */
2167         c = PI_O_4_BITS[1] >>> 32;
2168         d = PI_O_4_BITS[1] & 0xffffffffL;
2169         ac = a * c;
2170         bc = b * c;
2171         ad = a * d;
2172 
2173         /* Collect terms */
2174         ac += (bc + ad) >>> 32;
2175 
2176         bita = (prod2B & 0x8000000000000000L) != 0;
2177         bitb = (ac & 0x8000000000000000L ) != 0;
2178         prod2B += ac;
2179         bitsum = (prod2B & 0x8000000000000000L) != 0;
2180         /* Carry */
2181         if ( (bita && bitb) ||
2182                 ((bita || bitb) && !bitsum) ) {
2183             prod2A++;
2184         }
2185 
2186         /* Multiply inputB by pio4bits[0] */
2187         a = prodB >>> 32;
2188         b = prodB & 0xffffffffL;
2189         c = PI_O_4_BITS[0] >>> 32;
2190         d = PI_O_4_BITS[0] & 0xffffffffL;
2191         ac = a * c;
2192         bc = b * c;
2193         ad = a * d;
2194 
2195         /* Collect terms */
2196         ac += (bc + ad) >>> 32;
2197 
2198         bita = (prod2B & 0x8000000000000000L) != 0;
2199         bitb = (ac & 0x8000000000000000L ) != 0;
2200         prod2B += ac;
2201         bitsum = (prod2B & 0x8000000000000000L) != 0;
2202         /* Carry */
2203         if ( (bita && bitb) ||
2204                 ((bita || bitb) && !bitsum) ) {
2205             prod2A++;
2206         }
2207 
2208         /* Convert to double */
2209         double tmpA = (prod2A >>> 12) / TWO_POWER_52;  // High order 52 bits
2210         double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
2211 
2212         double sumA = tmpA + tmpB;
2213         double sumB = -(sumA - tmpA - tmpB);
2214 
2215         /* Multiply by PI/2 and return */
2216         result[0] = intPart;
2217         result[1] = sumA * 2.0;
2218         result[2] = sumB * 2.0;
2219     }
2220 
2221     /**
2222      * Sine function.
2223      *
2224      * @param x Argument.
2225      * @return sin(x)
2226      */
2227     public static double sin(double x) {
2228         boolean negative = false;
2229         int quadrant = 0;
2230         double xa;
2231         double xb = 0.0;
2232 
2233         /* Take absolute value of the input */
2234         xa = x;
2235         if (x < 0) {
2236             negative = true;
2237             xa = -xa;
2238         }
2239 
2240         /* Check for zero and negative zero */
2241         if (xa == 0.0) {
2242             long bits = Double.doubleToRawLongBits(x);
2243             if (bits < 0) {
2244                 return -0.0;
2245             }
2246             return 0.0;
2247         }
2248 
2249         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2250             return Double.NaN;
2251         }
2252 
2253         /* Perform any argument reduction */
2254         if (xa > 3294198.0) {
2255             // PI * (2**20)
2256             // Argument too big for CodyWaite reduction.  Must use
2257             // PayneHanek.
2258             double reduceResults[] = new double[3];
2259             reducePayneHanek(xa, reduceResults);
2260             quadrant = ((int) reduceResults[0]) & 3;
2261             xa = reduceResults[1];
2262             xb = reduceResults[2];
2263         } else if (xa > 1.5707963267948966) {
2264             final CodyWaite cw = new CodyWaite(xa);
2265             quadrant = cw.getK() & 3;
2266             xa = cw.getRemA();
2267             xb = cw.getRemB();
2268         }
2269 
2270         if (negative) {
2271             quadrant ^= 2;  // Flip bit 1
2272         }
2273 
2274         switch (quadrant) {
2275             case 0:
2276                 return sinQ(xa, xb);
2277             case 1:
2278                 return cosQ(xa, xb);
2279             case 2:
2280                 return -sinQ(xa, xb);
2281             case 3:
2282                 return -cosQ(xa, xb);
2283             default:
2284                 return Double.NaN;
2285         }
2286     }
2287 
2288     /**
2289      * Cosine function.
2290      *
2291      * @param x Argument.
2292      * @return cos(x)
2293      */
2294     public static double cos(double x) {
2295         int quadrant = 0;
2296 
2297         /* Take absolute value of the input */
2298         double xa = x;
2299         if (x < 0) {
2300             xa = -xa;
2301         }
2302 
2303         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2304             return Double.NaN;
2305         }
2306 
2307         /* Perform any argument reduction */
2308         double xb = 0;
2309         if (xa > 3294198.0) {
2310             // PI * (2**20)
2311             // Argument too big for CodyWaite reduction.  Must use
2312             // PayneHanek.
2313             double reduceResults[] = new double[3];
2314             reducePayneHanek(xa, reduceResults);
2315             quadrant = ((int) reduceResults[0]) & 3;
2316             xa = reduceResults[1];
2317             xb = reduceResults[2];
2318         } else if (xa > 1.5707963267948966) {
2319             final CodyWaite cw = new CodyWaite(xa);
2320             quadrant = cw.getK() & 3;
2321             xa = cw.getRemA();
2322             xb = cw.getRemB();
2323         }
2324 
2325         //if (negative)
2326         //  quadrant = (quadrant + 2) % 4;
2327 
2328         switch (quadrant) {
2329             case 0:
2330                 return cosQ(xa, xb);
2331             case 1:
2332                 return -sinQ(xa, xb);
2333             case 2:
2334                 return -cosQ(xa, xb);
2335             case 3:
2336                 return sinQ(xa, xb);
2337             default:
2338                 return Double.NaN;
2339         }
2340     }
2341 
2342     /**
2343      * Tangent function.
2344      *
2345      * @param x Argument.
2346      * @return tan(x)
2347      */
2348     public static double tan(double x) {
2349         boolean negative = false;
2350         int quadrant = 0;
2351 
2352         /* Take absolute value of the input */
2353         double xa = x;
2354         if (x < 0) {
2355             negative = true;
2356             xa = -xa;
2357         }
2358 
2359         /* Check for zero and negative zero */
2360         if (xa == 0.0) {
2361             long bits = Double.doubleToRawLongBits(x);
2362             if (bits < 0) {
2363                 return -0.0;
2364             }
2365             return 0.0;
2366         }
2367 
2368         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2369             return Double.NaN;
2370         }
2371 
2372         /* Perform any argument reduction */
2373         double xb = 0;
2374         if (xa > 3294198.0) {
2375             // PI * (2**20)
2376             // Argument too big for CodyWaite reduction.  Must use
2377             // PayneHanek.
2378             double reduceResults[] = new double[3];
2379             reducePayneHanek(xa, reduceResults);
2380             quadrant = ((int) reduceResults[0]) & 3;
2381             xa = reduceResults[1];
2382             xb = reduceResults[2];
2383         } else if (xa > 1.5707963267948966) {
2384             final CodyWaite cw = new CodyWaite(xa);
2385             quadrant = cw.getK() & 3;
2386             xa = cw.getRemA();
2387             xb = cw.getRemB();
2388         }
2389 
2390         if (xa > 1.5) {
2391             // Accuracy suffers between 1.5 and PI/2
2392             final double pi2a = 1.5707963267948966;
2393             final double pi2b = 6.123233995736766E-17;
2394 
2395             final double a = pi2a - xa;
2396             double b = -(a - pi2a + xa);
2397             b += pi2b - xb;
2398 
2399             xa = a + b;
2400             xb = -(xa - a - b);
2401             quadrant ^= 1;
2402             negative ^= true;
2403         }
2404 
2405         double result;
2406         if ((quadrant & 1) == 0) {
2407             result = tanQ(xa, xb, false);
2408         } else {
2409             result = -tanQ(xa, xb, true);
2410         }
2411 
2412         if (negative) {
2413             result = -result;
2414         }
2415 
2416         return result;
2417     }
2418 
2419     /**
2420      * Arctangent function
2421      *  @param x a number
2422      *  @return atan(x)
2423      */
2424     public static double atan(double x) {
2425         return atan(x, 0.0, false);
2426     }
2427 
2428     /** Internal helper function to compute arctangent.
2429      * @param xa number from which arctangent is requested
2430      * @param xb extra bits for x (may be 0.0)
2431      * @param leftPlane if true, result angle must be put in the left half plane
2432      * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
2433      */
2434     private static double atan(double xa, double xb, boolean leftPlane) {
2435         if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2436             return leftPlane ? copySign(Math.PI, xa) : xa;
2437         }
2438 
2439         final boolean negate;
2440         if (xa < 0) {
2441             // negative
2442             xa = -xa;
2443             xb = -xb;
2444             negate = true;
2445         } else {
2446             negate = false;
2447         }
2448 
2449         if (xa > 1.633123935319537E16) { // Very large input
2450             return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2451         }
2452 
2453         /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
2454         final int idx;
2455         if (xa < 1) {
2456             idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2457         } else {
2458             final double oneOverXa = 1 / xa;
2459             idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2460         }
2461 
2462         final double ttA = TANGENT_TABLE_A[idx];
2463         final double ttB = TANGENT_TABLE_B[idx];
2464 
2465         double epsA = xa - ttA;
2466         double epsB = -(epsA - xa + ttA);
2467         epsB += xb - ttB;
2468 
2469         double temp = epsA + epsB;
2470         epsB = -(temp - epsA - epsB);
2471         epsA = temp;
2472 
2473         /* Compute eps = eps / (1.0 + xa*tangent) */
2474         temp = xa * HEX_40000000;
2475         double ya = xa + temp - temp;
2476         double yb = xb + xa - ya;
2477         xa = ya;
2478         xb += yb;
2479 
2480         //if (idx > 8 || idx == 0)
2481         if (idx == 0) {
2482             /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2483             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2484             final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2485             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
2486             ya = epsA * denom;
2487             yb = epsB * denom;
2488         } else {
2489             double temp2 = xa * ttA;
2490             double za = 1d + temp2;
2491             double zb = -(za - 1d - temp2);
2492             temp2 = xb * ttA + xa * ttB;
2493             temp = za + temp2;
2494             zb += -(temp - za - temp2);
2495             za = temp;
2496 
2497             zb += xb * ttB;
2498             ya = epsA / za;
2499 
2500             temp = ya * HEX_40000000;
2501             final double yaa = (ya + temp) - temp;
2502             final double yab = ya - yaa;
2503 
2504             temp = za * HEX_40000000;
2505             final double zaa = (za + temp) - temp;
2506             final double zab = za - zaa;
2507 
2508             /* Correct for rounding in division */
2509             yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2510 
2511             yb += -epsA * zb / za / za;
2512             yb += epsB / za;
2513         }
2514 
2515 
2516         epsA = ya;
2517         epsB = yb;
2518 
2519         /* Evaluate polynomial */
2520         final double epsA2 = epsA * epsA;
2521 
2522         /*
2523     yb = -0.09001346640161823;
2524     yb = yb * epsA2 + 0.11110718400605211;
2525     yb = yb * epsA2 + -0.1428571349122913;
2526     yb = yb * epsA2 + 0.19999999999273194;
2527     yb = yb * epsA2 + -0.33333333333333093;
2528     yb = yb * epsA2 * epsA;
2529          */
2530 
2531         yb = 0.07490822288864472;
2532         yb = yb * epsA2 - 0.09088450866185192;
2533         yb = yb * epsA2 + 0.11111095942313305;
2534         yb = yb * epsA2 - 0.1428571423679182;
2535         yb = yb * epsA2 + 0.19999999999923582;
2536         yb = yb * epsA2 - 0.33333333333333287;
2537         yb = yb * epsA2 * epsA;
2538 
2539 
2540         ya = epsA;
2541 
2542         temp = ya + yb;
2543         yb = -(temp - ya - yb);
2544         ya = temp;
2545 
2546         /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
2547         yb += epsB / (1d + epsA * epsA);
2548 
2549         final double eighths = EIGHTHS[idx];
2550 
2551         //result = yb + eighths[idx] + ya;
2552         double za = eighths + ya;
2553         double zb = -(za - eighths - ya);
2554         temp = za + yb;
2555         zb += -(temp - za - yb);
2556         za = temp;
2557 
2558         double result = za + zb;
2559 
2560         if (leftPlane) {
2561             // Result is in the left plane
2562             final double resultb = -(result - za - zb);
2563             final double pia = 1.5707963267948966 * 2;
2564             final double pib = 6.123233995736766E-17 * 2;
2565 
2566             za = pia - result;
2567             zb = -(za - pia + result);
2568             zb += pib - resultb;
2569 
2570             result = za + zb;
2571         }
2572 
2573 
2574         if (negate ^ leftPlane) {
2575             result = -result;
2576         }
2577 
2578         return result;
2579     }
2580 
2581     /**
2582      * Two arguments arctangent function
2583      * @param y ordinate
2584      * @param x abscissa
2585      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
2586      */
2587     public static double atan2(double y, double x) {
2588         if (x != x || y != y) {
2589             return Double.NaN;
2590         }
2591 
2592         if (y == 0) {
2593             final double result = x * y;
2594             final double invx = 1d / x;
2595             final double invy = 1d / y;
2596 
2597             if (invx == 0) { // X is infinite
2598                 if (x > 0) {
2599                     return y; // return +/- 0.0
2600                 } else {
2601                     return copySign(Math.PI, y);
2602                 }
2603             }
2604 
2605             if (x < 0 || invx < 0) {
2606                 if (y < 0 || invy < 0) {
2607                     return -Math.PI;
2608                 } else {
2609                     return Math.PI;
2610                 }
2611             } else {
2612                 return result;
2613             }
2614         }
2615 
2616         // y cannot now be zero
2617 
2618         if (y == Double.POSITIVE_INFINITY) {
2619             if (x == Double.POSITIVE_INFINITY) {
2620                 return Math.PI * F_1_4;
2621             }
2622 
2623             if (x == Double.NEGATIVE_INFINITY) {
2624                 return Math.PI * F_3_4;
2625             }
2626 
2627             return Math.PI * F_1_2;
2628         }
2629 
2630         if (y == Double.NEGATIVE_INFINITY) {
2631             if (x == Double.POSITIVE_INFINITY) {
2632                 return -Math.PI * F_1_4;
2633             }
2634 
2635             if (x == Double.NEGATIVE_INFINITY) {
2636                 return -Math.PI * F_3_4;
2637             }
2638 
2639             return -Math.PI * F_1_2;
2640         }
2641 
2642         if (x == Double.POSITIVE_INFINITY) {
2643             if (y > 0 || 1 / y > 0) {
2644                 return 0d;
2645             }
2646 
2647             if (y < 0 || 1 / y < 0) {
2648                 return -0d;
2649             }
2650         }
2651 
2652         if (x == Double.NEGATIVE_INFINITY)
2653         {
2654             if (y > 0.0 || 1 / y > 0.0) {
2655                 return Math.PI;
2656             }
2657 
2658             if (y < 0 || 1 / y < 0) {
2659                 return -Math.PI;
2660             }
2661         }
2662 
2663         // Neither y nor x can be infinite or NAN here
2664 
2665         if (x == 0) {
2666             if (y > 0 || 1 / y > 0) {
2667                 return Math.PI * F_1_2;
2668             }
2669 
2670             if (y < 0 || 1 / y < 0) {
2671                 return -Math.PI * F_1_2;
2672             }
2673         }
2674 
2675         // Compute ratio r = y/x
2676         final double r = y / x;
2677         if (Double.isInfinite(r)) { // bypass calculations that can create NaN
2678             return atan(r, 0, x < 0);
2679         }
2680 
2681         double ra = doubleHighPart(r);
2682         double rb = r - ra;
2683 
2684         // Split x
2685         final double xa = doubleHighPart(x);
2686         final double xb = x - xa;
2687 
2688         rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
2689 
2690         final double temp = ra + rb;
2691         rb = -(temp - ra - rb);
2692         ra = temp;
2693 
2694         if (ra == 0) { // Fix up the sign so atan works correctly
2695             ra = copySign(0d, y);
2696         }
2697 
2698         // Call atan
2699         final double result = atan(ra, rb, x < 0);
2700 
2701         return result;
2702     }
2703 
2704     /** Compute the arc sine of a number.
2705      * @param x number on which evaluation is done
2706      * @return arc sine of x
2707      */
2708     public static double asin(double x) {
2709       if (x != x) {
2710           return Double.NaN;
2711       }
2712 
2713       if (x > 1.0 || x < -1.0) {
2714           return Double.NaN;
2715       }
2716 
2717       if (x == 1.0) {
2718           return Math.PI/2.0;
2719       }
2720 
2721       if (x == -1.0) {
2722           return -Math.PI/2.0;
2723       }
2724 
2725       if (x == 0.0) { // Matches +/- 0.0; return correct sign
2726           return x;
2727       }
2728 
2729       /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
2730 
2731       /* Split x */
2732       double temp = x * HEX_40000000;
2733       final double xa = x + temp - temp;
2734       final double xb = x - xa;
2735 
2736       /* Square it */
2737       double ya = xa*xa;
2738       double yb = xa*xb*2.0 + xb*xb;
2739 
2740       /* Subtract from 1 */
2741       ya = -ya;
2742       yb = -yb;
2743 
2744       double za = 1.0 + ya;
2745       double zb = -(za - 1.0 - ya);
2746 
2747       temp = za + yb;
2748       zb += -(temp - za - yb);
2749       za = temp;
2750 
2751       /* Square root */
2752       double y;
2753       y = sqrt(za);
2754       temp = y * HEX_40000000;
2755       ya = y + temp - temp;
2756       yb = y - ya;
2757 
2758       /* Extend precision of sqrt */
2759       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2760 
2761       /* Contribution of zb to sqrt */
2762       double dx = zb / (2.0*y);
2763 
2764       // Compute ratio r = x/y
2765       double r = x/y;
2766       temp = r * HEX_40000000;
2767       double ra = r + temp - temp;
2768       double rb = r - ra;
2769 
2770       rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y;  // Correct for rounding in division
2771       rb += -x * dx / y / y;  // Add in effect additional bits of sqrt.
2772 
2773       temp = ra + rb;
2774       rb = -(temp - ra - rb);
2775       ra = temp;
2776 
2777       return atan(ra, rb, false);
2778     }
2779 
2780     /** Compute the arc cosine of a number.
2781      * @param x number on which evaluation is done
2782      * @return arc cosine of x
2783      */
2784     public static double acos(double x) {
2785       if (x != x) {
2786           return Double.NaN;
2787       }
2788 
2789       if (x > 1.0 || x < -1.0) {
2790           return Double.NaN;
2791       }
2792 
2793       if (x == -1.0) {
2794           return Math.PI;
2795       }
2796 
2797       if (x == 1.0) {
2798           return 0.0;
2799       }
2800 
2801       if (x == 0) {
2802           return Math.PI/2.0;
2803       }
2804 
2805       /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
2806 
2807       /* Split x */
2808       double temp = x * HEX_40000000;
2809       final double xa = x + temp - temp;
2810       final double xb = x - xa;
2811 
2812       /* Square it */
2813       double ya = xa*xa;
2814       double yb = xa*xb*2.0 + xb*xb;
2815 
2816       /* Subtract from 1 */
2817       ya = -ya;
2818       yb = -yb;
2819 
2820       double za = 1.0 + ya;
2821       double zb = -(za - 1.0 - ya);
2822 
2823       temp = za + yb;
2824       zb += -(temp - za - yb);
2825       za = temp;
2826 
2827       /* Square root */
2828       double y = sqrt(za);
2829       temp = y * HEX_40000000;
2830       ya = y + temp - temp;
2831       yb = y - ya;
2832 
2833       /* Extend precision of sqrt */
2834       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2835 
2836       /* Contribution of zb to sqrt */
2837       yb += zb / (2.0*y);
2838       y = ya+yb;
2839       yb = -(y - ya - yb);
2840 
2841       // Compute ratio r = y/x
2842       double r = y/x;
2843 
2844       // Did r overflow?
2845       if (Double.isInfinite(r)) { // x is effectively zero
2846           return Math.PI/2; // so return the appropriate value
2847       }
2848 
2849       double ra = doubleHighPart(r);
2850       double rb = r - ra;
2851 
2852       rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x;  // Correct for rounding in division
2853       rb += yb / x;  // Add in effect additional bits of sqrt.
2854 
2855       temp = ra + rb;
2856       rb = -(temp - ra - rb);
2857       ra = temp;
2858 
2859       return atan(ra, rb, x<0);
2860     }
2861 
2862     /** Compute the cubic root of a number.
2863      * @param x number on which evaluation is done
2864      * @return cubic root of x
2865      */
2866     public static double cbrt(double x) {
2867       /* Convert input double to bits */
2868       long inbits = Double.doubleToRawLongBits(x);
2869       int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2870       boolean subnormal = false;
2871 
2872       if (exponent == -1023) {
2873           if (x == 0) {
2874               return x;
2875           }
2876 
2877           /* Subnormal, so normalize */
2878           subnormal = true;
2879           x *= 1.8014398509481984E16;  // 2^54
2880           inbits = Double.doubleToRawLongBits(x);
2881           exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2882       }
2883 
2884       if (exponent == 1024) {
2885           // Nan or infinity.  Don't care which.
2886           return x;
2887       }
2888 
2889       /* Divide the exponent by 3 */
2890       int exp3 = exponent / 3;
2891 
2892       /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
2893       double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
2894                                           (long)(((exp3 + 1023) & 0x7ff)) << 52);
2895 
2896       /* This will be a number between 1 and 2 */
2897       final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
2898 
2899       /* Estimate the cube root of mant by polynomial */
2900       double est = -0.010714690733195933;
2901       est = est * mant + 0.0875862700108075;
2902       est = est * mant + -0.3058015757857271;
2903       est = est * mant + 0.7249995199969751;
2904       est = est * mant + 0.5039018405998233;
2905 
2906       est *= CBRTTWO[exponent % 3 + 2];
2907 
2908       // est should now be good to about 15 bits of precision.   Do 2 rounds of
2909       // Newton's method to get closer,  this should get us full double precision
2910       // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
2911       final double xs = x / (p2*p2*p2);
2912       est += (xs - est*est*est) / (3*est*est);
2913       est += (xs - est*est*est) / (3*est*est);
2914 
2915       // Do one round of Newton's method in extended precision to get the last bit right.
2916       double temp = est * HEX_40000000;
2917       double ya = est + temp - temp;
2918       double yb = est - ya;
2919 
2920       double za = ya * ya;
2921       double zb = ya * yb * 2.0 + yb * yb;
2922       temp = za * HEX_40000000;
2923       double temp2 = za + temp - temp;
2924       zb += za - temp2;
2925       za = temp2;
2926 
2927       zb = za * yb + ya * zb + zb * yb;
2928       za *= ya;
2929 
2930       double na = xs - za;
2931       double nb = -(na - xs + za);
2932       nb -= zb;
2933 
2934       est += (na+nb)/(3*est*est);
2935 
2936       /* Scale by a power of two, so this is exact. */
2937       est *= p2;
2938 
2939       if (subnormal) {
2940           est *= 3.814697265625E-6;  // 2^-18
2941       }
2942 
2943       return est;
2944     }
2945 
2946     /**
2947      *  Convert degrees to radians, with error of less than 0.5 ULP
2948      *  @param x angle in degrees
2949      *  @return x converted into radians
2950      */
2951     public static double toRadians(double x)
2952     {
2953         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2954             return x;
2955         }
2956 
2957         // These are PI/180 split into high and low order bits
2958         final double facta = 0.01745329052209854;
2959         final double factb = 1.997844754509471E-9;
2960 
2961         double xa = doubleHighPart(x);
2962         double xb = x - xa;
2963 
2964         double result = xb * factb + xb * facta + xa * factb + xa * facta;
2965         if (result == 0) {
2966             result *= x; // ensure correct sign if calculation underflows
2967         }
2968         return result;
2969     }
2970 
2971     /**
2972      *  Convert radians to degrees, with error of less than 0.5 ULP
2973      *  @param x angle in radians
2974      *  @return x converted into degrees
2975      */
2976     public static double toDegrees(double x)
2977     {
2978         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2979             return x;
2980         }
2981 
2982         // These are 180/PI split into high and low order bits
2983         final double facta = 57.2957763671875;
2984         final double factb = 3.145894820876798E-6;
2985 
2986         double xa = doubleHighPart(x);
2987         double xb = x - xa;
2988 
2989         return xb * factb + xb * facta + xa * factb + xa * facta;
2990     }
2991 
2992     /**
2993      * Absolute value.
2994      * @param x number from which absolute value is requested
2995      * @return abs(x)
2996      */
2997     public static int abs(final int x) {
2998         final int i = x >>> 31;
2999         return (x ^ (~i + 1)) + i;
3000     }
3001 
3002     /**
3003      * Absolute value.
3004      * @param x number from which absolute value is requested
3005      * @return abs(x)
3006      */
3007     public static long abs(final long x) {
3008         final long l = x >>> 63;
3009         // l is one if x negative zero else
3010         // ~l+1 is zero if x is positive, -1 if x is negative
3011         // x^(~l+1) is x is x is positive, ~x if x is negative
3012         // add around
3013         return (x ^ (~l + 1)) + l;
3014     }
3015 
3016     /**
3017      * Absolute value.
3018      * @param x number from which absolute value is requested
3019      * @return abs(x)
3020      */
3021     public static float abs(final float x) {
3022         return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3023     }
3024 
3025     /**
3026      * Absolute value.
3027      * @param x number from which absolute value is requested
3028      * @return abs(x)
3029      */
3030     public static double abs(double x) {
3031         return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3032     }
3033 
3034     /**
3035      * Compute least significant bit (Unit in Last Position) for a number.
3036      * @param x number from which ulp is requested
3037      * @return ulp(x)
3038      */
3039     public static double ulp(double x) {
3040         if (Double.isInfinite(x)) {
3041             return Double.POSITIVE_INFINITY;
3042         }
3043         return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3044     }
3045 
3046     /**
3047      * Compute least significant bit (Unit in Last Position) for a number.
3048      * @param x number from which ulp is requested
3049      * @return ulp(x)
3050      */
3051     public static float ulp(float x) {
3052         if (Float.isInfinite(x)) {
3053             return Float.POSITIVE_INFINITY;
3054         }
3055         return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
3056     }
3057 
3058     /**
3059      * Multiply a double number by a power of 2.
3060      * @param d number to multiply
3061      * @param n power of 2
3062      * @return d &times; 2<sup>n</sup>
3063      */
3064     public static double scalb(final double d, final int n) {
3065 
3066         // first simple and fast handling when 2^n can be represented using normal numbers
3067         if ((n > -1023) && (n < 1024)) {
3068             return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3069         }
3070 
3071         // handle special cases
3072         if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3073             return d;
3074         }
3075         if (n < -2098) {
3076             return (d > 0) ? 0.0 : -0.0;
3077         }
3078         if (n > 2097) {
3079             return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3080         }
3081 
3082         // decompose d
3083         final long bits = Double.doubleToRawLongBits(d);
3084         final long sign = bits & 0x8000000000000000L;
3085         int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
3086         long mantissa   = bits & 0x000fffffffffffffL;
3087 
3088         // compute scaled exponent
3089         int scaledExponent = exponent + n;
3090 
3091         if (n < 0) {
3092             // we are really in the case n <= -1023
3093             if (scaledExponent > 0) {
3094                 // both the input and the result are normal numbers, we only adjust the exponent
3095                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3096             } else if (scaledExponent > -53) {
3097                 // the input is a normal number and the result is a subnormal number
3098 
3099                 // recover the hidden mantissa bit
3100                 mantissa |= 1L << 52;
3101 
3102                 // scales down complete mantissa, hence losing least significant bits
3103                 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3104                 mantissa >>>= 1 - scaledExponent;
3105                 if (mostSignificantLostBit != 0) {
3106                     // we need to add 1 bit to round up the result
3107                     mantissa++;
3108                 }
3109                 return Double.longBitsToDouble(sign | mantissa);
3110 
3111             } else {
3112                 // no need to compute the mantissa, the number scales down to 0
3113                 return (sign == 0L) ? 0.0 : -0.0;
3114             }
3115         } else {
3116             // we are really in the case n >= 1024
3117             if (exponent == 0) {
3118 
3119                 // the input number is subnormal, normalize it
3120                 while ((mantissa >>> 52) != 1) {
3121                     mantissa <<= 1;
3122                     --scaledExponent;
3123                 }
3124                 ++scaledExponent;
3125                 mantissa &= 0x000fffffffffffffL;
3126 
3127                 if (scaledExponent < 2047) {
3128                     return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3129                 } else {
3130                     return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3131                 }
3132 
3133             } else if (scaledExponent < 2047) {
3134                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3135             } else {
3136                 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3137             }
3138         }
3139 
3140     }
3141 
3142     /**
3143      * Multiply a float number by a power of 2.
3144      * @param f number to multiply
3145      * @param n power of 2
3146      * @return f &times; 2<sup>n</sup>
3147      */
3148     public static float scalb(final float f, final int n) {
3149 
3150         // first simple and fast handling when 2^n can be represented using normal numbers
3151         if ((n > -127) && (n < 128)) {
3152             return f * Float.intBitsToFloat((n + 127) << 23);
3153         }
3154 
3155         // handle special cases
3156         if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3157             return f;
3158         }
3159         if (n < -277) {
3160             return (f > 0) ? 0.0f : -0.0f;
3161         }
3162         if (n > 276) {
3163             return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3164         }
3165 
3166         // decompose f
3167         final int bits = Float.floatToIntBits(f);
3168         final int sign = bits & 0x80000000;
3169         int  exponent  = (bits >>> 23) & 0xff;
3170         int mantissa   = bits & 0x007fffff;
3171 
3172         // compute scaled exponent
3173         int scaledExponent = exponent + n;
3174 
3175         if (n < 0) {
3176             // we are really in the case n <= -127
3177             if (scaledExponent > 0) {
3178                 // both the input and the result are normal numbers, we only adjust the exponent
3179                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3180             } else if (scaledExponent > -24) {
3181                 // the input is a normal number and the result is a subnormal number
3182 
3183                 // recover the hidden mantissa bit
3184                 mantissa |= 1 << 23;
3185 
3186                 // scales down complete mantissa, hence losing least significant bits
3187                 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3188                 mantissa >>>= 1 - scaledExponent;
3189                 if (mostSignificantLostBit != 0) {
3190                     // we need to add 1 bit to round up the result
3191                     mantissa++;
3192                 }
3193                 return Float.intBitsToFloat(sign | mantissa);
3194 
3195             } else {
3196                 // no need to compute the mantissa, the number scales down to 0
3197                 return (sign == 0) ? 0.0f : -0.0f;
3198             }
3199         } else {
3200             // we are really in the case n >= 128
3201             if (exponent == 0) {
3202 
3203                 // the input number is subnormal, normalize it
3204                 while ((mantissa >>> 23) != 1) {
3205                     mantissa <<= 1;
3206                     --scaledExponent;
3207                 }
3208                 ++scaledExponent;
3209                 mantissa &= 0x007fffff;
3210 
3211                 if (scaledExponent < 255) {
3212                     return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3213                 } else {
3214                     return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3215                 }
3216 
3217             } else if (scaledExponent < 255) {
3218                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3219             } else {
3220                 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3221             }
3222         }
3223 
3224     }
3225 
3226     /**
3227      * Get the next machine representable number after a number, moving
3228      * in the direction of another number.
3229      * <p>
3230      * The ordering is as follows (increasing):
3231      * <ul>
3232      * <li>-INFINITY</li>
3233      * <li>-MAX_VALUE</li>
3234      * <li>-MIN_VALUE</li>
3235      * <li>-0.0</li>
3236      * <li>+0.0</li>
3237      * <li>+MIN_VALUE</li>
3238      * <li>+MAX_VALUE</li>
3239      * <li>+INFINITY</li>
3240      * <li></li>
3241      * <p>
3242      * If arguments compare equal, then the second argument is returned.
3243      * <p>
3244      * If {@code direction} is greater than {@code d},
3245      * the smallest machine representable number strictly greater than
3246      * {@code d} is returned; if less, then the largest representable number
3247      * strictly less than {@code d} is returned.</p>
3248      * <p>
3249      * If {@code d} is infinite and direction does not
3250      * bring it back to finite numbers, it is returned unchanged.</p>
3251      *
3252      * @param d base number
3253      * @param direction (the only important thing is whether
3254      * {@code direction} is greater or smaller than {@code d})
3255      * @return the next machine representable number in the specified direction
3256      */
3257     public static double nextAfter(double d, double direction) {
3258 
3259         // handling of some important special cases
3260         if (Double.isNaN(d) || Double.isNaN(direction)) {
3261             return Double.NaN;
3262         } else if (d == direction) {
3263             return direction;
3264         } else if (Double.isInfinite(d)) {
3265             return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3266         } else if (d == 0) {
3267             return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3268         }
3269         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3270         // are handled just as normal numbers
3271         // can use raw bits since already dealt with infinity and NaN
3272         final long bits = Double.doubleToRawLongBits(d);
3273         final long sign = bits & 0x8000000000000000L;
3274         if ((direction < d) ^ (sign == 0L)) {
3275             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3276         } else {
3277             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3278         }
3279 
3280     }
3281 
3282     /**
3283      * Get the next machine representable number after a number, moving
3284      * in the direction of another number.
3285      * <p>
3286      * The ordering is as follows (increasing):
3287      * <ul>
3288      * <li>-INFINITY</li>
3289      * <li>-MAX_VALUE</li>
3290      * <li>-MIN_VALUE</li>
3291      * <li>-0.0</li>
3292      * <li>+0.0</li>
3293      * <li>+MIN_VALUE</li>
3294      * <li>+MAX_VALUE</li>
3295      * <li>+INFINITY</li>
3296      * <li></li>
3297      * <p>
3298      * If arguments compare equal, then the second argument is returned.
3299      * <p>
3300      * If {@code direction} is greater than {@code f},
3301      * the smallest machine representable number strictly greater than
3302      * {@code f} is returned; if less, then the largest representable number
3303      * strictly less than {@code f} is returned.</p>
3304      * <p>
3305      * If {@code f} is infinite and direction does not
3306      * bring it back to finite numbers, it is returned unchanged.</p>
3307      *
3308      * @param f base number
3309      * @param direction (the only important thing is whether
3310      * {@code direction} is greater or smaller than {@code f})
3311      * @return the next machine representable number in the specified direction
3312      */
3313     public static float nextAfter(final float f, final double direction) {
3314 
3315         // handling of some important special cases
3316         if (Double.isNaN(f) || Double.isNaN(direction)) {
3317             return Float.NaN;
3318         } else if (f == direction) {
3319             return (float) direction;
3320         } else if (Float.isInfinite(f)) {
3321             return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3322         } else if (f == 0f) {
3323             return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3324         }
3325         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3326         // are handled just as normal numbers
3327 
3328         final int bits = Float.floatToIntBits(f);
3329         final int sign = bits & 0x80000000;
3330         if ((direction < f) ^ (sign == 0)) {
3331             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3332         } else {
3333             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3334         }
3335 
3336     }
3337 
3338     /** Get the largest whole number smaller than x.
3339      * @param x number from which floor is requested
3340      * @return a double number f such that f is an integer f <= x < f + 1.0
3341      */
3342     public static double floor(double x) {
3343         long y;
3344 
3345         if (x != x) { // NaN
3346             return x;
3347         }
3348 
3349         if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3350             return x;
3351         }
3352 
3353         y = (long) x;
3354         if (x < 0 && y != x) {
3355             y--;
3356         }
3357 
3358         if (y == 0) {
3359             return x*y;
3360         }
3361 
3362         return y;
3363     }
3364 
3365     /** Get the smallest whole number larger than x.
3366      * @param x number from which ceil is requested
3367      * @return a double number c such that c is an integer c - 1.0 < x <= c
3368      */
3369     public static double ceil(double x) {
3370         double y;
3371 
3372         if (x != x) { // NaN
3373             return x;
3374         }
3375 
3376         y = floor(x);
3377         if (y == x) {
3378             return y;
3379         }
3380 
3381         y += 1.0;
3382 
3383         if (y == 0) {
3384             return x*y;
3385         }
3386 
3387         return y;
3388     }
3389 
3390     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3391      * @param x number from which nearest whole number is requested
3392      * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
3393      */
3394     public static double rint(double x) {
3395         double y = floor(x);
3396         double d = x - y;
3397 
3398         if (d > 0.5) {
3399             if (y == -1.0) {
3400                 return -0.0; // Preserve sign of operand
3401             }
3402             return y+1.0;
3403         }
3404         if (d < 0.5) {
3405             return y;
3406         }
3407 
3408         /* half way, round to even */
3409         long z = (long) y;
3410         return (z & 1) == 0 ? y : y + 1.0;
3411     }
3412 
3413     /** Get the closest long to x.
3414      * @param x number from which closest long is requested
3415      * @return closest long to x
3416      */
3417     public static long round(double x) {
3418         return (long) floor(x + 0.5);
3419     }
3420 
3421     /** Get the closest int to x.
3422      * @param x number from which closest int is requested
3423      * @return closest int to x
3424      */
3425     public static int round(final float x) {
3426         return (int) floor(x + 0.5f);
3427     }
3428 
3429     /** Compute the minimum of two values
3430      * @param a first value
3431      * @param b second value
3432      * @return a if a is lesser or equal to b, b otherwise
3433      */
3434     public static int min(final int a, final int b) {
3435         return (a <= b) ? a : b;
3436     }
3437 
3438     /** Compute the minimum of two values
3439      * @param a first value
3440      * @param b second value
3441      * @return a if a is lesser or equal to b, b otherwise
3442      */
3443     public static long min(final long a, final long b) {
3444         return (a <= b) ? a : b;
3445     }
3446 
3447     /** Compute the minimum of two values
3448      * @param a first value
3449      * @param b second value
3450      * @return a if a is lesser or equal to b, b otherwise
3451      */
3452     public static float min(final float a, final float b) {
3453         if (a > b) {
3454             return b;
3455         }
3456         if (a < b) {
3457             return a;
3458         }
3459         /* if either arg is NaN, return NaN */
3460         if (a != b) {
3461             return Float.NaN;
3462         }
3463         /* min(+0.0,-0.0) == -0.0 */
3464         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3465         int bits = Float.floatToRawIntBits(a);
3466         if (bits == 0x80000000) {
3467             return a;
3468         }
3469         return b;
3470     }
3471 
3472     /** Compute the minimum of two values
3473      * @param a first value
3474      * @param b second value
3475      * @return a if a is lesser or equal to b, b otherwise
3476      */
3477     public static double min(final double a, final double b) {
3478         if (a > b) {
3479             return b;
3480         }
3481         if (a < b) {
3482             return a;
3483         }
3484         /* if either arg is NaN, return NaN */
3485         if (a != b) {
3486             return Double.NaN;
3487         }
3488         /* min(+0.0,-0.0) == -0.0 */
3489         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3490         long bits = Double.doubleToRawLongBits(a);
3491         if (bits == 0x8000000000000000L) {
3492             return a;
3493         }
3494         return b;
3495     }
3496 
3497     /** Compute the maximum of two values
3498      * @param a first value
3499      * @param b second value
3500      * @return b if a is lesser or equal to b, a otherwise
3501      */
3502     public static int max(final int a, final int b) {
3503         return (a <= b) ? b : a;
3504     }
3505 
3506     /** Compute the maximum of two values
3507      * @param a first value
3508      * @param b second value
3509      * @return b if a is lesser or equal to b, a otherwise
3510      */
3511     public static long max(final long a, final long b) {
3512         return (a <= b) ? b : a;
3513     }
3514 
3515     /** Compute the maximum of two values
3516      * @param a first value
3517      * @param b second value
3518      * @return b if a is lesser or equal to b, a otherwise
3519      */
3520     public static float max(final float a, final float b) {
3521         if (a > b) {
3522             return a;
3523         }
3524         if (a < b) {
3525             return b;
3526         }
3527         /* if either arg is NaN, return NaN */
3528         if (a != b) {
3529             return Float.NaN;
3530         }
3531         /* min(+0.0,-0.0) == -0.0 */
3532         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3533         int bits = Float.floatToRawIntBits(a);
3534         if (bits == 0x80000000) {
3535             return b;
3536         }
3537         return a;
3538     }
3539 
3540     /** Compute the maximum of two values
3541      * @param a first value
3542      * @param b second value
3543      * @return b if a is lesser or equal to b, a otherwise
3544      */
3545     public static double max(final double a, final double b) {
3546         if (a > b) {
3547             return a;
3548         }
3549         if (a < b) {
3550             return b;
3551         }
3552         /* if either arg is NaN, return NaN */
3553         if (a != b) {
3554             return Double.NaN;
3555         }
3556         /* min(+0.0,-0.0) == -0.0 */
3557         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3558         long bits = Double.doubleToRawLongBits(a);
3559         if (bits == 0x8000000000000000L) {
3560             return b;
3561         }
3562         return a;
3563     }
3564 
3565     /**
3566      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
3567      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br/>
3568      * avoiding intermediate overflow or underflow.
3569      *
3570      * <ul>
3571      * <li> If either argument is infinite, then the result is positive infinity.</li>
3572      * <li> else, if either argument is NaN then the result is NaN.</li>
3573      * </ul>
3574      *
3575      * @param x a value
3576      * @param y a value
3577      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
3578      */
3579     public static double hypot(final double x, final double y) {
3580         if (Double.isInfinite(x) || Double.isInfinite(y)) {
3581             return Double.POSITIVE_INFINITY;
3582         } else if (Double.isNaN(x) || Double.isNaN(y)) {
3583             return Double.NaN;
3584         } else {
3585 
3586             final int expX = getExponent(x);
3587             final int expY = getExponent(y);
3588             if (expX > expY + 27) {
3589                 // y is neglectible with respect to x
3590                 return abs(x);
3591             } else if (expY > expX + 27) {
3592                 // x is neglectible with respect to y
3593                 return abs(y);
3594             } else {
3595 
3596                 // find an intermediate scale to avoid both overflow and underflow
3597                 final int middleExp = (expX + expY) / 2;
3598 
3599                 // scale parameters without losing precision
3600                 final double scaledX = scalb(x, -middleExp);
3601                 final double scaledY = scalb(y, -middleExp);
3602 
3603                 // compute scaled hypotenuse
3604                 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
3605 
3606                 // remove scaling
3607                 return scalb(scaledH, middleExp);
3608 
3609             }
3610 
3611         }
3612     }
3613 
3614     /**
3615      * Computes the remainder as prescribed by the IEEE 754 standard.
3616      * The remainder value is mathematically equal to {@code x - y*n}
3617      * where {@code n} is the mathematical integer closest to the exact mathematical value
3618      * of the quotient {@code x/y}.
3619      * If two mathematical integers are equally close to {@code x/y} then
3620      * {@code n} is the integer that is even.
3621      * <p>
3622      * <ul>
3623      * <li>If either operand is NaN, the result is NaN.</li>
3624      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
3625      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
3626      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
3627      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
3628      * </ul>
3629      * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder}
3630      * @param dividend the number to be divided
3631      * @param divisor the number by which to divide
3632      * @return the remainder, rounded
3633      */
3634     public static double IEEEremainder(double dividend, double divisor) {
3635         return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
3636     }
3637 
3638     /**
3639      * Returns the first argument with the sign of the second argument.
3640      * A NaN {@code sign} argument is treated as positive.
3641      *
3642      * @param magnitude the value to return
3643      * @param sign the sign for the returned value
3644      * @return the magnitude with the same sign as the {@code sign} argument
3645      */
3646     public static double copySign(double magnitude, double sign){
3647         // The highest order bit is going to be zero if the
3648         // highest order bit of m and s is the same and one otherwise.
3649         // So (m^s) will be positive if both m and s have the same sign
3650         // and negative otherwise.
3651         final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
3652         final long s = Double.doubleToRawLongBits(sign);
3653         if ((m^s) >= 0) {
3654             return magnitude;
3655         }
3656         return -magnitude; // flip sign
3657     }
3658 
3659     /**
3660      * Returns the first argument with the sign of the second argument.
3661      * A NaN {@code sign} argument is treated as positive.
3662      *
3663      * @param magnitude the value to return
3664      * @param sign the sign for the returned value
3665      * @return the magnitude with the same sign as the {@code sign} argument
3666      */
3667     public static float copySign(float magnitude, float sign){
3668         // The highest order bit is going to be zero if the
3669         // highest order bit of m and s is the same and one otherwise.
3670         // So (m^s) will be positive if both m and s have the same sign
3671         // and negative otherwise.
3672         final int m = Float.floatToRawIntBits(magnitude);
3673         final int s = Float.floatToRawIntBits(sign);
3674         if ((m^s) >= 0) {
3675             return magnitude;
3676         }
3677         return -magnitude; // flip sign
3678     }
3679 
3680     /**
3681      * Return the exponent of a double number, removing the bias.
3682      * <p>
3683      * For double numbers of the form 2<sup>x</sup>, the unbiased
3684      * exponent is exactly x.
3685      * </p>
3686      * @param d number from which exponent is requested
3687      * @return exponent for d in IEEE754 representation, without bias
3688      */
3689     public static int getExponent(final double d) {
3690         // NaN and Infinite will return 1024 anywho so can use raw bits
3691         return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
3692     }
3693 
3694     /**
3695      * Return the exponent of a float number, removing the bias.
3696      * <p>
3697      * For float numbers of the form 2<sup>x</sup>, the unbiased
3698      * exponent is exactly x.
3699      * </p>
3700      * @param f number from which exponent is requested
3701      * @return exponent for d in IEEE754 representation, without bias
3702      */
3703     public static int getExponent(final float f) {
3704         // NaN and Infinite will return the same exponent anywho so can use raw bits
3705         return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
3706     }
3707 
3708     /**
3709      * Print out contents of arrays, and check the length.
3710      * <p>used to generate the preset arrays originally.</p>
3711      * @param a unused
3712      */
3713     public static void main(String[] a) {
3714         PrintStream out = System.out;
3715         FastMathCalc.printarray(out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
3716         FastMathCalc.printarray(out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
3717         FastMathCalc.printarray(out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
3718         FastMathCalc.printarray(out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
3719         FastMathCalc.printarray(out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
3720         FastMathCalc.printarray(out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
3721         FastMathCalc.printarray(out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
3722         FastMathCalc.printarray(out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
3723         FastMathCalc.printarray(out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
3724         FastMathCalc.printarray(out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
3725         FastMathCalc.printarray(out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
3726     }
3727 
3728     /** Enclose large data table in nested static class so it's only loaded on first access. */
3729     private static class ExpIntTable {
3730         /** Exponential evaluated at integer values,
3731          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
3732          */
3733         private static final double[] EXP_INT_TABLE_A;
3734         /** Exponential evaluated at integer values,
3735          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
3736          */
3737         private static final double[] EXP_INT_TABLE_B;
3738 
3739         static {
3740             if (RECOMPUTE_TABLES_AT_RUNTIME) {
3741                 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
3742                 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
3743 
3744                 final double tmp[] = new double[2];
3745                 final double recip[] = new double[2];
3746 
3747                 // Populate expIntTable
3748                 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
3749                     FastMathCalc.expint(i, tmp);
3750                     EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
3751                     EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
3752 
3753                     if (i != 0) {
3754                         // Negative integer powers
3755                         FastMathCalc.splitReciprocal(tmp, recip);
3756                         EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
3757                         EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
3758                     }
3759                 }
3760             } else {
3761                 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
3762                 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
3763             }
3764         }
3765     }
3766 
3767     /** Enclose large data table in nested static class so it's only loaded on first access. */
3768     private static class ExpFracTable {
3769         /** Exponential over the range of 0 - 1 in increments of 2^-10
3770          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
3771          * 1024 = 2^10
3772          */
3773         private static final double[] EXP_FRAC_TABLE_A;
3774         /** Exponential over the range of 0 - 1 in increments of 2^-10
3775          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
3776          */
3777         private static final double[] EXP_FRAC_TABLE_B;
3778 
3779         static {
3780             if (RECOMPUTE_TABLES_AT_RUNTIME) {
3781                 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
3782                 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
3783 
3784                 final double tmp[] = new double[2];
3785 
3786                 // Populate expFracTable
3787                 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
3788                 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
3789                     FastMathCalc.slowexp(i * factor, tmp);
3790                     EXP_FRAC_TABLE_A[i] = tmp[0];
3791                     EXP_FRAC_TABLE_B[i] = tmp[1];
3792                 }
3793             } else {
3794                 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
3795                 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
3796             }
3797         }
3798     }
3799 
3800     /** Enclose large data table in nested static class so it's only loaded on first access. */
3801     private static class lnMant {
3802         /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
3803         private static final double[][] LN_MANT;
3804 
3805         static {
3806             if (RECOMPUTE_TABLES_AT_RUNTIME) {
3807                 LN_MANT = new double[FastMath.LN_MANT_LEN][];
3808 
3809                 // Populate lnMant table
3810                 for (int i = 0; i < LN_MANT.length; i++) {
3811                     final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
3812                     LN_MANT[i] = FastMathCalc.slowLog(d);
3813                 }
3814             } else {
3815                 LN_MANT = FastMathLiteralArrays.loadLnMant();
3816             }
3817         }
3818     }
3819 
3820     /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
3821     private static class CodyWaite {
3822         /** k */
3823         private final int finalK;
3824         /** remA */
3825         private final double finalRemA;
3826         /** remB */
3827         private final double finalRemB;
3828 
3829         /**
3830          * @param xa Argument.
3831          */
3832         CodyWaite(double xa) {
3833             // Estimate k.
3834             //k = (int)(xa / 1.5707963267948966);
3835             int k = (int)(xa * 0.6366197723675814);
3836 
3837             // Compute remainder.
3838             double remA;
3839             double remB;
3840             while (true) {
3841                 double a = -k * 1.570796251296997;
3842                 remA = xa + a;
3843                 remB = -(remA - xa - a);
3844 
3845                 a = -k * 7.549789948768648E-8;
3846                 double b = remA;
3847                 remA = a + b;
3848                 remB += -(remA - b - a);
3849 
3850                 a = -k * 6.123233995736766E-17;
3851                 b = remA;
3852                 remA = a + b;
3853                 remB += -(remA - b - a);
3854 
3855                 if (remA > 0) {
3856                     break;
3857                 }
3858 
3859                 // Remainder is negative, so decrement k and try again.
3860                 // This should only happen if the input is very close
3861                 // to an even multiple of pi/2.
3862                 --k;
3863             }
3864 
3865             this.finalK = k;
3866             this.finalRemA = remA;
3867             this.finalRemB = remB;
3868         }
3869 
3870         /**
3871          * @return k
3872          */
3873         int getK() {
3874             return finalK;
3875         }
3876         /**
3877          * @return remA
3878          */
3879         double getRemA() {
3880             return finalRemA;
3881         }
3882         /**
3883          * @return remB
3884          */
3885         double getRemB() {
3886             return finalRemB;
3887         }
3888     }
3889 }