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 1517542 2013-08-26 14:24:18Z erans $
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 = 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 = zb * epsilon;
1030         zb = 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 = 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 = b + d;
1328 
1329         c = a + lnza;
1330         d = -(c - a - lnza);
1331         a = c;
1332         b = b + d;
1333 
1334         c = a + LN_2_B*exp;
1335         d = -(c - a - LN_2_B*exp);
1336         a = c;
1337         b = b + d;
1338 
1339         c = a + lnm[1];
1340         d = -(c - a - lnm[1]);
1341         a = c;
1342         b = b + d;
1343 
1344         c = a + lnzb;
1345         d = -(c - a - lnzb);
1346         a = c;
1347         b = 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 = 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 = 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 = b + d;
1772 
1773         t = costA * sinEpsA;
1774         c = a + t;
1775         d = -(c - a - t);
1776         a = c;
1777         b = 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 = 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 = b + d;
1916 
1917         t = costA*sinEpsA;
1918         c = a + t;
1919         d = -(c - a - t);
1920         a = c;
1921         b = b + d;
1922 
1923         b = b + sintA*cosEpsB + costA*sinEpsB;
1924         b = 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 = b + d;
1938 
1939         t = -sintA*sinEpsA;
1940         c = a + t;
1941         d = -(c - a - t);
1942         a = c;
1943         b = b + d;
1944 
1945         b = b + costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
1946         b = 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 = prodB + (bc << 32);
2070         prodA = 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 = 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 = prod2B + (bc << 32);
2156         prod2A = 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 = 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 = 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         boolean negate = false;
2436         int idx;
2437 
2438         if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2439             return leftPlane ? copySign(Math.PI, xa) : xa;
2440         }
2441 
2442         if (xa < 0) {
2443             // negative
2444             xa = -xa;
2445             xb = -xb;
2446             negate = true;
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         if (xa < 1) {
2455             idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2456         } else {
2457             final double oneOverXa = 1 / xa;
2458             idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2459         }
2460         double epsA = xa - TANGENT_TABLE_A[idx];
2461         double epsB = -(epsA - xa + TANGENT_TABLE_A[idx]);
2462         epsB += xb - TANGENT_TABLE_B[idx];
2463 
2464         double temp = epsA + epsB;
2465         epsB = -(temp - epsA - epsB);
2466         epsA = temp;
2467 
2468         /* Compute eps = eps / (1.0 + xa*tangent) */
2469         temp = xa * HEX_40000000;
2470         double ya = xa + temp - temp;
2471         double yb = xb + xa - ya;
2472         xa = ya;
2473         xb += yb;
2474 
2475         //if (idx > 8 || idx == 0)
2476         if (idx == 0) {
2477             /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2478             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2479             final double denom = 1d / (1d + (xa + xb) * (TANGENT_TABLE_A[idx] + TANGENT_TABLE_B[idx]));
2480             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
2481             ya = epsA * denom;
2482             yb = epsB * denom;
2483         } else {
2484             double temp2 = xa * TANGENT_TABLE_A[idx];
2485             double za = 1d + temp2;
2486             double zb = -(za - 1d - temp2);
2487             temp2 = xb * TANGENT_TABLE_A[idx] + xa * TANGENT_TABLE_B[idx];
2488             temp = za + temp2;
2489             zb += -(temp - za - temp2);
2490             za = temp;
2491 
2492             zb += xb * TANGENT_TABLE_B[idx];
2493             ya = epsA / za;
2494 
2495             temp = ya * HEX_40000000;
2496             final double yaa = (ya + temp) - temp;
2497             final double yab = ya - yaa;
2498 
2499             temp = za * HEX_40000000;
2500             final double zaa = (za + temp) - temp;
2501             final double zab = za - zaa;
2502 
2503             /* Correct for rounding in division */
2504             yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2505 
2506             yb += -epsA * zb / za / za;
2507             yb += epsB / za;
2508         }
2509 
2510 
2511         epsA = ya;
2512         epsB = yb;
2513 
2514         /* Evaluate polynomial */
2515         final double epsA2 = epsA * epsA;
2516 
2517         /*
2518     yb = -0.09001346640161823;
2519     yb = yb * epsA2 + 0.11110718400605211;
2520     yb = yb * epsA2 + -0.1428571349122913;
2521     yb = yb * epsA2 + 0.19999999999273194;
2522     yb = yb * epsA2 + -0.33333333333333093;
2523     yb = yb * epsA2 * epsA;
2524          */
2525 
2526         yb = 0.07490822288864472;
2527         yb = yb * epsA2 + -0.09088450866185192;
2528         yb = yb * epsA2 + 0.11111095942313305;
2529         yb = yb * epsA2 + -0.1428571423679182;
2530         yb = yb * epsA2 + 0.19999999999923582;
2531         yb = yb * epsA2 + -0.33333333333333287;
2532         yb = yb * epsA2 * epsA;
2533 
2534 
2535         ya = epsA;
2536 
2537         temp = ya + yb;
2538         yb = -(temp - ya - yb);
2539         ya = temp;
2540 
2541         /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
2542         yb += epsB / (1d + epsA * epsA);
2543 
2544         //result = yb + eighths[idx] + ya;
2545         double za = EIGHTHS[idx] + ya;
2546         double zb = -(za - EIGHTHS[idx] - ya);
2547         temp = za + yb;
2548         zb += -(temp - za - yb);
2549         za = temp;
2550 
2551         double result = za + zb;
2552 
2553         if (leftPlane) {
2554             // Result is in the left plane
2555             final double resultb = -(result - za - zb);
2556             final double pia = 1.5707963267948966 * 2;
2557             final double pib = 6.123233995736766E-17 * 2;
2558 
2559             za = pia - result;
2560             zb = -(za - pia + result);
2561             zb += pib - resultb;
2562 
2563             result = za + zb;
2564         }
2565 
2566 
2567         if (negate ^ leftPlane) {
2568             result = -result;
2569         }
2570 
2571         return result;
2572     }
2573 
2574     /**
2575      * Two arguments arctangent function
2576      * @param y ordinate
2577      * @param x abscissa
2578      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
2579      */
2580     public static double atan2(double y, double x) {
2581         if (x != x || y != y) {
2582             return Double.NaN;
2583         }
2584 
2585         if (y == 0) {
2586             final double result = x * y;
2587             final double invx = 1d / x;
2588             final double invy = 1d / y;
2589 
2590             if (invx == 0) { // X is infinite
2591                 if (x > 0) {
2592                     return y; // return +/- 0.0
2593                 } else {
2594                     return copySign(Math.PI, y);
2595                 }
2596             }
2597 
2598             if (x < 0 || invx < 0) {
2599                 if (y < 0 || invy < 0) {
2600                     return -Math.PI;
2601                 } else {
2602                     return Math.PI;
2603                 }
2604             } else {
2605                 return result;
2606             }
2607         }
2608 
2609         // y cannot now be zero
2610 
2611         if (y == Double.POSITIVE_INFINITY) {
2612             if (x == Double.POSITIVE_INFINITY) {
2613                 return Math.PI * F_1_4;
2614             }
2615 
2616             if (x == Double.NEGATIVE_INFINITY) {
2617                 return Math.PI * F_3_4;
2618             }
2619 
2620             return Math.PI * F_1_2;
2621         }
2622 
2623         if (y == Double.NEGATIVE_INFINITY) {
2624             if (x == Double.POSITIVE_INFINITY) {
2625                 return -Math.PI * F_1_4;
2626             }
2627 
2628             if (x == Double.NEGATIVE_INFINITY) {
2629                 return -Math.PI * F_3_4;
2630             }
2631 
2632             return -Math.PI * F_1_2;
2633         }
2634 
2635         if (x == Double.POSITIVE_INFINITY) {
2636             if (y > 0 || 1 / y > 0) {
2637                 return 0d;
2638             }
2639 
2640             if (y < 0 || 1 / y < 0) {
2641                 return -0d;
2642             }
2643         }
2644 
2645         if (x == Double.NEGATIVE_INFINITY)
2646         {
2647             if (y > 0.0 || 1 / y > 0.0) {
2648                 return Math.PI;
2649             }
2650 
2651             if (y < 0 || 1 / y < 0) {
2652                 return -Math.PI;
2653             }
2654         }
2655 
2656         // Neither y nor x can be infinite or NAN here
2657 
2658         if (x == 0) {
2659             if (y > 0 || 1 / y > 0) {
2660                 return Math.PI * F_1_2;
2661             }
2662 
2663             if (y < 0 || 1 / y < 0) {
2664                 return -Math.PI * F_1_2;
2665             }
2666         }
2667 
2668         // Compute ratio r = y/x
2669         final double r = y / x;
2670         if (Double.isInfinite(r)) { // bypass calculations that can create NaN
2671             return atan(r, 0, x < 0);
2672         }
2673 
2674         double ra = doubleHighPart(r);
2675         double rb = r - ra;
2676 
2677         // Split x
2678         final double xa = doubleHighPart(x);
2679         final double xb = x - xa;
2680 
2681         rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
2682 
2683         final double temp = ra + rb;
2684         rb = -(temp - ra - rb);
2685         ra = temp;
2686 
2687         if (ra == 0) { // Fix up the sign so atan works correctly
2688             ra = copySign(0d, y);
2689         }
2690 
2691         // Call atan
2692         final double result = atan(ra, rb, x < 0);
2693 
2694         return result;
2695     }
2696 
2697     /** Compute the arc sine of a number.
2698      * @param x number on which evaluation is done
2699      * @return arc sine of x
2700      */
2701     public static double asin(double x) {
2702       if (x != x) {
2703           return Double.NaN;
2704       }
2705 
2706       if (x > 1.0 || x < -1.0) {
2707           return Double.NaN;
2708       }
2709 
2710       if (x == 1.0) {
2711           return Math.PI/2.0;
2712       }
2713 
2714       if (x == -1.0) {
2715           return -Math.PI/2.0;
2716       }
2717 
2718       if (x == 0.0) { // Matches +/- 0.0; return correct sign
2719           return x;
2720       }
2721 
2722       /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
2723 
2724       /* Split x */
2725       double temp = x * HEX_40000000;
2726       final double xa = x + temp - temp;
2727       final double xb = x - xa;
2728 
2729       /* Square it */
2730       double ya = xa*xa;
2731       double yb = xa*xb*2.0 + xb*xb;
2732 
2733       /* Subtract from 1 */
2734       ya = -ya;
2735       yb = -yb;
2736 
2737       double za = 1.0 + ya;
2738       double zb = -(za - 1.0 - ya);
2739 
2740       temp = za + yb;
2741       zb += -(temp - za - yb);
2742       za = temp;
2743 
2744       /* Square root */
2745       double y;
2746       y = sqrt(za);
2747       temp = y * HEX_40000000;
2748       ya = y + temp - temp;
2749       yb = y - ya;
2750 
2751       /* Extend precision of sqrt */
2752       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2753 
2754       /* Contribution of zb to sqrt */
2755       double dx = zb / (2.0*y);
2756 
2757       // Compute ratio r = x/y
2758       double r = x/y;
2759       temp = r * HEX_40000000;
2760       double ra = r + temp - temp;
2761       double rb = r - ra;
2762 
2763       rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y;  // Correct for rounding in division
2764       rb += -x * dx / y / y;  // Add in effect additional bits of sqrt.
2765 
2766       temp = ra + rb;
2767       rb = -(temp - ra - rb);
2768       ra = temp;
2769 
2770       return atan(ra, rb, false);
2771     }
2772 
2773     /** Compute the arc cosine of a number.
2774      * @param x number on which evaluation is done
2775      * @return arc cosine of x
2776      */
2777     public static double acos(double x) {
2778       if (x != x) {
2779           return Double.NaN;
2780       }
2781 
2782       if (x > 1.0 || x < -1.0) {
2783           return Double.NaN;
2784       }
2785 
2786       if (x == -1.0) {
2787           return Math.PI;
2788       }
2789 
2790       if (x == 1.0) {
2791           return 0.0;
2792       }
2793 
2794       if (x == 0) {
2795           return Math.PI/2.0;
2796       }
2797 
2798       /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
2799 
2800       /* Split x */
2801       double temp = x * HEX_40000000;
2802       final double xa = x + temp - temp;
2803       final double xb = x - xa;
2804 
2805       /* Square it */
2806       double ya = xa*xa;
2807       double yb = xa*xb*2.0 + xb*xb;
2808 
2809       /* Subtract from 1 */
2810       ya = -ya;
2811       yb = -yb;
2812 
2813       double za = 1.0 + ya;
2814       double zb = -(za - 1.0 - ya);
2815 
2816       temp = za + yb;
2817       zb += -(temp - za - yb);
2818       za = temp;
2819 
2820       /* Square root */
2821       double y = sqrt(za);
2822       temp = y * HEX_40000000;
2823       ya = y + temp - temp;
2824       yb = y - ya;
2825 
2826       /* Extend precision of sqrt */
2827       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2828 
2829       /* Contribution of zb to sqrt */
2830       yb += zb / (2.0*y);
2831       y = ya+yb;
2832       yb = -(y - ya - yb);
2833 
2834       // Compute ratio r = y/x
2835       double r = y/x;
2836 
2837       // Did r overflow?
2838       if (Double.isInfinite(r)) { // x is effectively zero
2839           return Math.PI/2; // so return the appropriate value
2840       }
2841 
2842       double ra = doubleHighPart(r);
2843       double rb = r - ra;
2844 
2845       rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x;  // Correct for rounding in division
2846       rb += yb / x;  // Add in effect additional bits of sqrt.
2847 
2848       temp = ra + rb;
2849       rb = -(temp - ra - rb);
2850       ra = temp;
2851 
2852       return atan(ra, rb, x<0);
2853     }
2854 
2855     /** Compute the cubic root of a number.
2856      * @param x number on which evaluation is done
2857      * @return cubic root of x
2858      */
2859     public static double cbrt(double x) {
2860       /* Convert input double to bits */
2861       long inbits = Double.doubleToRawLongBits(x);
2862       int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2863       boolean subnormal = false;
2864 
2865       if (exponent == -1023) {
2866           if (x == 0) {
2867               return x;
2868           }
2869 
2870           /* Subnormal, so normalize */
2871           subnormal = true;
2872           x *= 1.8014398509481984E16;  // 2^54
2873           inbits = Double.doubleToRawLongBits(x);
2874           exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2875       }
2876 
2877       if (exponent == 1024) {
2878           // Nan or infinity.  Don't care which.
2879           return x;
2880       }
2881 
2882       /* Divide the exponent by 3 */
2883       int exp3 = exponent / 3;
2884 
2885       /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
2886       double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
2887                                           (long)(((exp3 + 1023) & 0x7ff)) << 52);
2888 
2889       /* This will be a number between 1 and 2 */
2890       final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
2891 
2892       /* Estimate the cube root of mant by polynomial */
2893       double est = -0.010714690733195933;
2894       est = est * mant + 0.0875862700108075;
2895       est = est * mant + -0.3058015757857271;
2896       est = est * mant + 0.7249995199969751;
2897       est = est * mant + 0.5039018405998233;
2898 
2899       est *= CBRTTWO[exponent % 3 + 2];
2900 
2901       // est should now be good to about 15 bits of precision.   Do 2 rounds of
2902       // Newton's method to get closer,  this should get us full double precision
2903       // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
2904       final double xs = x / (p2*p2*p2);
2905       est += (xs - est*est*est) / (3*est*est);
2906       est += (xs - est*est*est) / (3*est*est);
2907 
2908       // Do one round of Newton's method in extended precision to get the last bit right.
2909       double temp = est * HEX_40000000;
2910       double ya = est + temp - temp;
2911       double yb = est - ya;
2912 
2913       double za = ya * ya;
2914       double zb = ya * yb * 2.0 + yb * yb;
2915       temp = za * HEX_40000000;
2916       double temp2 = za + temp - temp;
2917       zb += za - temp2;
2918       za = temp2;
2919 
2920       zb = za * yb + ya * zb + zb * yb;
2921       za = za * ya;
2922 
2923       double na = xs - za;
2924       double nb = -(na - xs + za);
2925       nb -= zb;
2926 
2927       est += (na+nb)/(3*est*est);
2928 
2929       /* Scale by a power of two, so this is exact. */
2930       est *= p2;
2931 
2932       if (subnormal) {
2933           est *= 3.814697265625E-6;  // 2^-18
2934       }
2935 
2936       return est;
2937     }
2938 
2939     /**
2940      *  Convert degrees to radians, with error of less than 0.5 ULP
2941      *  @param x angle in degrees
2942      *  @return x converted into radians
2943      */
2944     public static double toRadians(double x)
2945     {
2946         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2947             return x;
2948         }
2949 
2950         // These are PI/180 split into high and low order bits
2951         final double facta = 0.01745329052209854;
2952         final double factb = 1.997844754509471E-9;
2953 
2954         double xa = doubleHighPart(x);
2955         double xb = x - xa;
2956 
2957         double result = xb * factb + xb * facta + xa * factb + xa * facta;
2958         if (result == 0) {
2959             result = result * x; // ensure correct sign if calculation underflows
2960         }
2961         return result;
2962     }
2963 
2964     /**
2965      *  Convert radians to degrees, with error of less than 0.5 ULP
2966      *  @param x angle in radians
2967      *  @return x converted into degrees
2968      */
2969     public static double toDegrees(double x)
2970     {
2971         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2972             return x;
2973         }
2974 
2975         // These are 180/PI split into high and low order bits
2976         final double facta = 57.2957763671875;
2977         final double factb = 3.145894820876798E-6;
2978 
2979         double xa = doubleHighPart(x);
2980         double xb = x - xa;
2981 
2982         return xb * factb + xb * facta + xa * factb + xa * facta;
2983     }
2984 
2985     /**
2986      * Absolute value.
2987      * @param x number from which absolute value is requested
2988      * @return abs(x)
2989      */
2990     public static int abs(final int x) {
2991         final int i = x >>> 31;
2992         return (x ^ (~i + 1)) + i;
2993     }
2994 
2995     /**
2996      * Absolute value.
2997      * @param x number from which absolute value is requested
2998      * @return abs(x)
2999      */
3000     public static long abs(final long x) {
3001         final long l = x >>> 63;
3002         // l is one if x negative zero else
3003         // ~l+1 is zero if x is positive, -1 if x is negative
3004         // x^(~l+1) is x is x is positive, ~x if x is negative
3005         // add around
3006         return (x ^ (~l + 1)) + l;
3007     }
3008 
3009     /**
3010      * Absolute value.
3011      * @param x number from which absolute value is requested
3012      * @return abs(x)
3013      */
3014     public static float abs(final float x) {
3015         return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3016     }
3017 
3018     /**
3019      * Absolute value.
3020      * @param x number from which absolute value is requested
3021      * @return abs(x)
3022      */
3023     public static double abs(double x) {
3024         return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3025     }
3026 
3027     /**
3028      * Compute least significant bit (Unit in Last Position) for a number.
3029      * @param x number from which ulp is requested
3030      * @return ulp(x)
3031      */
3032     public static double ulp(double x) {
3033         if (Double.isInfinite(x)) {
3034             return Double.POSITIVE_INFINITY;
3035         }
3036         return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3037     }
3038 
3039     /**
3040      * Compute least significant bit (Unit in Last Position) for a number.
3041      * @param x number from which ulp is requested
3042      * @return ulp(x)
3043      */
3044     public static float ulp(float x) {
3045         if (Float.isInfinite(x)) {
3046             return Float.POSITIVE_INFINITY;
3047         }
3048         return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
3049     }
3050 
3051     /**
3052      * Multiply a double number by a power of 2.
3053      * @param d number to multiply
3054      * @param n power of 2
3055      * @return d &times; 2<sup>n</sup>
3056      */
3057     public static double scalb(final double d, final int n) {
3058 
3059         // first simple and fast handling when 2^n can be represented using normal numbers
3060         if ((n > -1023) && (n < 1024)) {
3061             return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3062         }
3063 
3064         // handle special cases
3065         if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3066             return d;
3067         }
3068         if (n < -2098) {
3069             return (d > 0) ? 0.0 : -0.0;
3070         }
3071         if (n > 2097) {
3072             return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3073         }
3074 
3075         // decompose d
3076         final long bits = Double.doubleToRawLongBits(d);
3077         final long sign = bits & 0x8000000000000000L;
3078         int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
3079         long mantissa   = bits & 0x000fffffffffffffL;
3080 
3081         // compute scaled exponent
3082         int scaledExponent = exponent + n;
3083 
3084         if (n < 0) {
3085             // we are really in the case n <= -1023
3086             if (scaledExponent > 0) {
3087                 // both the input and the result are normal numbers, we only adjust the exponent
3088                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3089             } else if (scaledExponent > -53) {
3090                 // the input is a normal number and the result is a subnormal number
3091 
3092                 // recover the hidden mantissa bit
3093                 mantissa = mantissa | (1L << 52);
3094 
3095                 // scales down complete mantissa, hence losing least significant bits
3096                 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3097                 mantissa = mantissa >>> (1 - scaledExponent);
3098                 if (mostSignificantLostBit != 0) {
3099                     // we need to add 1 bit to round up the result
3100                     mantissa++;
3101                 }
3102                 return Double.longBitsToDouble(sign | mantissa);
3103 
3104             } else {
3105                 // no need to compute the mantissa, the number scales down to 0
3106                 return (sign == 0L) ? 0.0 : -0.0;
3107             }
3108         } else {
3109             // we are really in the case n >= 1024
3110             if (exponent == 0) {
3111 
3112                 // the input number is subnormal, normalize it
3113                 while ((mantissa >>> 52) != 1) {
3114                     mantissa = mantissa << 1;
3115                     --scaledExponent;
3116                 }
3117                 ++scaledExponent;
3118                 mantissa = mantissa & 0x000fffffffffffffL;
3119 
3120                 if (scaledExponent < 2047) {
3121                     return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3122                 } else {
3123                     return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3124                 }
3125 
3126             } else if (scaledExponent < 2047) {
3127                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3128             } else {
3129                 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3130             }
3131         }
3132 
3133     }
3134 
3135     /**
3136      * Multiply a float number by a power of 2.
3137      * @param f number to multiply
3138      * @param n power of 2
3139      * @return f &times; 2<sup>n</sup>
3140      */
3141     public static float scalb(final float f, final int n) {
3142 
3143         // first simple and fast handling when 2^n can be represented using normal numbers
3144         if ((n > -127) && (n < 128)) {
3145             return f * Float.intBitsToFloat((n + 127) << 23);
3146         }
3147 
3148         // handle special cases
3149         if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3150             return f;
3151         }
3152         if (n < -277) {
3153             return (f > 0) ? 0.0f : -0.0f;
3154         }
3155         if (n > 276) {
3156             return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3157         }
3158 
3159         // decompose f
3160         final int bits = Float.floatToIntBits(f);
3161         final int sign = bits & 0x80000000;
3162         int  exponent  = (bits >>> 23) & 0xff;
3163         int mantissa   = bits & 0x007fffff;
3164 
3165         // compute scaled exponent
3166         int scaledExponent = exponent + n;
3167 
3168         if (n < 0) {
3169             // we are really in the case n <= -127
3170             if (scaledExponent > 0) {
3171                 // both the input and the result are normal numbers, we only adjust the exponent
3172                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3173             } else if (scaledExponent > -24) {
3174                 // the input is a normal number and the result is a subnormal number
3175 
3176                 // recover the hidden mantissa bit
3177                 mantissa = mantissa | (1 << 23);
3178 
3179                 // scales down complete mantissa, hence losing least significant bits
3180                 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3181                 mantissa = mantissa >>> (1 - scaledExponent);
3182                 if (mostSignificantLostBit != 0) {
3183                     // we need to add 1 bit to round up the result
3184                     mantissa++;
3185                 }
3186                 return Float.intBitsToFloat(sign | mantissa);
3187 
3188             } else {
3189                 // no need to compute the mantissa, the number scales down to 0
3190                 return (sign == 0) ? 0.0f : -0.0f;
3191             }
3192         } else {
3193             // we are really in the case n >= 128
3194             if (exponent == 0) {
3195 
3196                 // the input number is subnormal, normalize it
3197                 while ((mantissa >>> 23) != 1) {
3198                     mantissa = mantissa << 1;
3199                     --scaledExponent;
3200                 }
3201                 ++scaledExponent;
3202                 mantissa = mantissa & 0x007fffff;
3203 
3204                 if (scaledExponent < 255) {
3205                     return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3206                 } else {
3207                     return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3208                 }
3209 
3210             } else if (scaledExponent < 255) {
3211                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3212             } else {
3213                 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3214             }
3215         }
3216 
3217     }
3218 
3219     /**
3220      * Get the next machine representable number after a number, moving
3221      * in the direction of another number.
3222      * <p>
3223      * The ordering is as follows (increasing):
3224      * <ul>
3225      * <li>-INFINITY</li>
3226      * <li>-MAX_VALUE</li>
3227      * <li>-MIN_VALUE</li>
3228      * <li>-0.0</li>
3229      * <li>+0.0</li>
3230      * <li>+MIN_VALUE</li>
3231      * <li>+MAX_VALUE</li>
3232      * <li>+INFINITY</li>
3233      * <li></li>
3234      * <p>
3235      * If arguments compare equal, then the second argument is returned.
3236      * <p>
3237      * If {@code direction} is greater than {@code d},
3238      * the smallest machine representable number strictly greater than
3239      * {@code d} is returned; if less, then the largest representable number
3240      * strictly less than {@code d} is returned.</p>
3241      * <p>
3242      * If {@code d} is infinite and direction does not
3243      * bring it back to finite numbers, it is returned unchanged.</p>
3244      *
3245      * @param d base number
3246      * @param direction (the only important thing is whether
3247      * {@code direction} is greater or smaller than {@code d})
3248      * @return the next machine representable number in the specified direction
3249      */
3250     public static double nextAfter(double d, double direction) {
3251 
3252         // handling of some important special cases
3253         if (Double.isNaN(d) || Double.isNaN(direction)) {
3254             return Double.NaN;
3255         } else if (d == direction) {
3256             return direction;
3257         } else if (Double.isInfinite(d)) {
3258             return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3259         } else if (d == 0) {
3260             return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3261         }
3262         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3263         // are handled just as normal numbers
3264         // can use raw bits since already dealt with infinity and NaN
3265         final long bits = Double.doubleToRawLongBits(d);
3266         final long sign = bits & 0x8000000000000000L;
3267         if ((direction < d) ^ (sign == 0L)) {
3268             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3269         } else {
3270             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3271         }
3272 
3273     }
3274 
3275     /**
3276      * Get the next machine representable number after a number, moving
3277      * in the direction of another number.
3278      * <p>
3279      * The ordering is as follows (increasing):
3280      * <ul>
3281      * <li>-INFINITY</li>
3282      * <li>-MAX_VALUE</li>
3283      * <li>-MIN_VALUE</li>
3284      * <li>-0.0</li>
3285      * <li>+0.0</li>
3286      * <li>+MIN_VALUE</li>
3287      * <li>+MAX_VALUE</li>
3288      * <li>+INFINITY</li>
3289      * <li></li>
3290      * <p>
3291      * If arguments compare equal, then the second argument is returned.
3292      * <p>
3293      * If {@code direction} is greater than {@code f},
3294      * the smallest machine representable number strictly greater than
3295      * {@code f} is returned; if less, then the largest representable number
3296      * strictly less than {@code f} is returned.</p>
3297      * <p>
3298      * If {@code f} is infinite and direction does not
3299      * bring it back to finite numbers, it is returned unchanged.</p>
3300      *
3301      * @param f base number
3302      * @param direction (the only important thing is whether
3303      * {@code direction} is greater or smaller than {@code f})
3304      * @return the next machine representable number in the specified direction
3305      */
3306     public static float nextAfter(final float f, final double direction) {
3307 
3308         // handling of some important special cases
3309         if (Double.isNaN(f) || Double.isNaN(direction)) {
3310             return Float.NaN;
3311         } else if (f == direction) {
3312             return (float) direction;
3313         } else if (Float.isInfinite(f)) {
3314             return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3315         } else if (f == 0f) {
3316             return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3317         }
3318         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3319         // are handled just as normal numbers
3320 
3321         final int bits = Float.floatToIntBits(f);
3322         final int sign = bits & 0x80000000;
3323         if ((direction < f) ^ (sign == 0)) {
3324             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3325         } else {
3326             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3327         }
3328 
3329     }
3330 
3331     /** Get the largest whole number smaller than x.
3332      * @param x number from which floor is requested
3333      * @return a double number f such that f is an integer f <= x < f + 1.0
3334      */
3335     public static double floor(double x) {
3336         long y;
3337 
3338         if (x != x) { // NaN
3339             return x;
3340         }
3341 
3342         if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3343             return x;
3344         }
3345 
3346         y = (long) x;
3347         if (x < 0 && y != x) {
3348             y--;
3349         }
3350 
3351         if (y == 0) {
3352             return x*y;
3353         }
3354 
3355         return y;
3356     }
3357 
3358     /** Get the smallest whole number larger than x.
3359      * @param x number from which ceil is requested
3360      * @return a double number c such that c is an integer c - 1.0 < x <= c
3361      */
3362     public static double ceil(double x) {
3363         double y;
3364 
3365         if (x != x) { // NaN
3366             return x;
3367         }
3368 
3369         y = floor(x);
3370         if (y == x) {
3371             return y;
3372         }
3373 
3374         y += 1.0;
3375 
3376         if (y == 0) {
3377             return x*y;
3378         }
3379 
3380         return y;
3381     }
3382 
3383     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3384      * @param x number from which nearest whole number is requested
3385      * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
3386      */
3387     public static double rint(double x) {
3388         double y = floor(x);
3389         double d = x - y;
3390 
3391         if (d > 0.5) {
3392             if (y == -1.0) {
3393                 return -0.0; // Preserve sign of operand
3394             }
3395             return y+1.0;
3396         }
3397         if (d < 0.5) {
3398             return y;
3399         }
3400 
3401         /* half way, round to even */
3402         long z = (long) y;
3403         return (z & 1) == 0 ? y : y + 1.0;
3404     }
3405 
3406     /** Get the closest long to x.
3407      * @param x number from which closest long is requested
3408      * @return closest long to x
3409      */
3410     public static long round(double x) {
3411         return (long) floor(x + 0.5);
3412     }
3413 
3414     /** Get the closest int to x.
3415      * @param x number from which closest int is requested
3416      * @return closest int to x
3417      */
3418     public static int round(final float x) {
3419         return (int) floor(x + 0.5f);
3420     }
3421 
3422     /** Compute the minimum of two values
3423      * @param a first value
3424      * @param b second value
3425      * @return a if a is lesser or equal to b, b otherwise
3426      */
3427     public static int min(final int a, final int b) {
3428         return (a <= b) ? a : b;
3429     }
3430 
3431     /** Compute the minimum of two values
3432      * @param a first value
3433      * @param b second value
3434      * @return a if a is lesser or equal to b, b otherwise
3435      */
3436     public static long min(final long a, final long b) {
3437         return (a <= b) ? a : b;
3438     }
3439 
3440     /** Compute the minimum of two values
3441      * @param a first value
3442      * @param b second value
3443      * @return a if a is lesser or equal to b, b otherwise
3444      */
3445     public static float min(final float a, final float b) {
3446         if (a > b) {
3447             return b;
3448         }
3449         if (a < b) {
3450             return a;
3451         }
3452         /* if either arg is NaN, return NaN */
3453         if (a != b) {
3454             return Float.NaN;
3455         }
3456         /* min(+0.0,-0.0) == -0.0 */
3457         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3458         int bits = Float.floatToRawIntBits(a);
3459         if (bits == 0x80000000) {
3460             return a;
3461         }
3462         return b;
3463     }
3464 
3465     /** Compute the minimum of two values
3466      * @param a first value
3467      * @param b second value
3468      * @return a if a is lesser or equal to b, b otherwise
3469      */
3470     public static double min(final double a, final double b) {
3471         if (a > b) {
3472             return b;
3473         }
3474         if (a < b) {
3475             return a;
3476         }
3477         /* if either arg is NaN, return NaN */
3478         if (a != b) {
3479             return Double.NaN;
3480         }
3481         /* min(+0.0,-0.0) == -0.0 */
3482         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3483         long bits = Double.doubleToRawLongBits(a);
3484         if (bits == 0x8000000000000000L) {
3485             return a;
3486         }
3487         return b;
3488     }
3489 
3490     /** Compute the maximum of two values
3491      * @param a first value
3492      * @param b second value
3493      * @return b if a is lesser or equal to b, a otherwise
3494      */
3495     public static int max(final int a, final int b) {
3496         return (a <= b) ? b : a;
3497     }
3498 
3499     /** Compute the maximum of two values
3500      * @param a first value
3501      * @param b second value
3502      * @return b if a is lesser or equal to b, a otherwise
3503      */
3504     public static long max(final long a, final long b) {
3505         return (a <= b) ? b : a;
3506     }
3507 
3508     /** Compute the maximum of two values
3509      * @param a first value
3510      * @param b second value
3511      * @return b if a is lesser or equal to b, a otherwise
3512      */
3513     public static float max(final float a, final float b) {
3514         if (a > b) {
3515             return a;
3516         }
3517         if (a < b) {
3518             return b;
3519         }
3520         /* if either arg is NaN, return NaN */
3521         if (a != b) {
3522             return Float.NaN;
3523         }
3524         /* min(+0.0,-0.0) == -0.0 */
3525         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3526         int bits = Float.floatToRawIntBits(a);
3527         if (bits == 0x80000000) {
3528             return b;
3529         }
3530         return a;
3531     }
3532 
3533     /** Compute the maximum of two values
3534      * @param a first value
3535      * @param b second value
3536      * @return b if a is lesser or equal to b, a otherwise
3537      */
3538     public static double max(final double a, final double b) {
3539         if (a > b) {
3540             return a;
3541         }
3542         if (a < b) {
3543             return b;
3544         }
3545         /* if either arg is NaN, return NaN */
3546         if (a != b) {
3547             return Double.NaN;
3548         }
3549         /* min(+0.0,-0.0) == -0.0 */
3550         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3551         long bits = Double.doubleToRawLongBits(a);
3552         if (bits == 0x8000000000000000L) {
3553             return b;
3554         }
3555         return a;
3556     }
3557 
3558     /**
3559      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
3560      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br/>
3561      * avoiding intermediate overflow or underflow.
3562      *
3563      * <ul>
3564      * <li> If either argument is infinite, then the result is positive infinity.</li>
3565      * <li> else, if either argument is NaN then the result is NaN.</li>
3566      * </ul>
3567      *
3568      * @param x a value
3569      * @param y a value
3570      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
3571      */
3572     public static double hypot(final double x, final double y) {
3573         if (Double.isInfinite(x) || Double.isInfinite(y)) {
3574             return Double.POSITIVE_INFINITY;
3575         } else if (Double.isNaN(x) || Double.isNaN(y)) {
3576             return Double.NaN;
3577         } else {
3578 
3579             final int expX = getExponent(x);
3580             final int expY = getExponent(y);
3581             if (expX > expY + 27) {
3582                 // y is neglectible with respect to x
3583                 return abs(x);
3584             } else if (expY > expX + 27) {
3585                 // x is neglectible with respect to y
3586                 return abs(y);
3587             } else {
3588 
3589                 // find an intermediate scale to avoid both overflow and underflow
3590                 final int middleExp = (expX + expY) / 2;
3591 
3592                 // scale parameters without losing precision
3593                 final double scaledX = scalb(x, -middleExp);
3594                 final double scaledY = scalb(y, -middleExp);
3595 
3596                 // compute scaled hypotenuse
3597                 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
3598 
3599                 // remove scaling
3600                 return scalb(scaledH, middleExp);
3601 
3602             }
3603 
3604         }
3605     }
3606 
3607     /**
3608      * Computes the remainder as prescribed by the IEEE 754 standard.
3609      * The remainder value is mathematically equal to {@code x - y*n}
3610      * where {@code n} is the mathematical integer closest to the exact mathematical value
3611      * of the quotient {@code x/y}.
3612      * If two mathematical integers are equally close to {@code x/y} then
3613      * {@code n} is the integer that is even.
3614      * <p>
3615      * <ul>
3616      * <li>If either operand is NaN, the result is NaN.</li>
3617      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
3618      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
3619      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
3620      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
3621      * </ul>
3622      * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder}
3623      * @param dividend the number to be divided
3624      * @param divisor the number by which to divide
3625      * @return the remainder, rounded
3626      */
3627     public static double IEEEremainder(double dividend, double divisor) {
3628         return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
3629     }
3630 
3631     /**
3632      * Returns the first argument with the sign of the second argument.
3633      * A NaN {@code sign} argument is treated as positive.
3634      *
3635      * @param magnitude the value to return
3636      * @param sign the sign for the returned value
3637      * @return the magnitude with the same sign as the {@code sign} argument
3638      */
3639     public static double copySign(double magnitude, double sign){
3640         // The highest order bit is going to be zero if the
3641         // highest order bit of m and s is the same and one otherwise.
3642         // So (m^s) will be positive if both m and s have the same sign
3643         // and negative otherwise.
3644         final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
3645         final long s = Double.doubleToRawLongBits(sign);
3646         if ((m^s) >= 0) {
3647             return magnitude;
3648         }
3649         return -magnitude; // flip sign
3650     }
3651 
3652     /**
3653      * Returns the first argument with the sign of the second argument.
3654      * A NaN {@code sign} argument is treated as positive.
3655      *
3656      * @param magnitude the value to return
3657      * @param sign the sign for the returned value
3658      * @return the magnitude with the same sign as the {@code sign} argument
3659      */
3660     public static float copySign(float magnitude, float sign){
3661         // The highest order bit is going to be zero if the
3662         // highest order bit of m and s is the same and one otherwise.
3663         // So (m^s) will be positive if both m and s have the same sign
3664         // and negative otherwise.
3665         final int m = Float.floatToRawIntBits(magnitude);
3666         final int s = Float.floatToRawIntBits(sign);
3667         if ((m^s) >= 0) {
3668             return magnitude;
3669         }
3670         return -magnitude; // flip sign
3671     }
3672 
3673     /**
3674      * Return the exponent of a double number, removing the bias.
3675      * <p>
3676      * For double numbers of the form 2<sup>x</sup>, the unbiased
3677      * exponent is exactly x.
3678      * </p>
3679      * @param d number from which exponent is requested
3680      * @return exponent for d in IEEE754 representation, without bias
3681      */
3682     public static int getExponent(final double d) {
3683         // NaN and Infinite will return 1024 anywho so can use raw bits
3684         return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
3685     }
3686 
3687     /**
3688      * Return the exponent of a float number, removing the bias.
3689      * <p>
3690      * For float numbers of the form 2<sup>x</sup>, the unbiased
3691      * exponent is exactly x.
3692      * </p>
3693      * @param f number from which exponent is requested
3694      * @return exponent for d in IEEE754 representation, without bias
3695      */
3696     public static int getExponent(final float f) {
3697         // NaN and Infinite will return the same exponent anywho so can use raw bits
3698         return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
3699     }
3700 
3701     /**
3702      * Print out contents of arrays, and check the length.
3703      * <p>used to generate the preset arrays originally.</p>
3704      * @param a unused
3705      */
3706     public static void main(String[] a) {
3707         PrintStream out = System.out;
3708         FastMathCalc.printarray(out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
3709         FastMathCalc.printarray(out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
3710         FastMathCalc.printarray(out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
3711         FastMathCalc.printarray(out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
3712         FastMathCalc.printarray(out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
3713         FastMathCalc.printarray(out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
3714         FastMathCalc.printarray(out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
3715         FastMathCalc.printarray(out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
3716         FastMathCalc.printarray(out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
3717         FastMathCalc.printarray(out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
3718         FastMathCalc.printarray(out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
3719     }
3720 
3721     /** Enclose large data table in nested static class so it's only loaded on first access. */
3722     private static class ExpIntTable {
3723         /** Exponential evaluated at integer values,
3724          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
3725          */
3726         private static final double[] EXP_INT_TABLE_A;
3727         /** Exponential evaluated at integer values,
3728          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
3729          */
3730         private static final double[] EXP_INT_TABLE_B;
3731 
3732         static {
3733             if (RECOMPUTE_TABLES_AT_RUNTIME) {
3734                 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
3735                 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
3736 
3737                 final double tmp[] = new double[2];
3738                 final double recip[] = new double[2];
3739 
3740                 // Populate expIntTable
3741                 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
3742                     FastMathCalc.expint(i, tmp);
3743                     EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
3744                     EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
3745 
3746                     if (i != 0) {
3747                         // Negative integer powers
3748                         FastMathCalc.splitReciprocal(tmp, recip);
3749                         EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
3750                         EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
3751                     }
3752                 }
3753             } else {
3754                 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
3755                 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
3756             }
3757         }
3758     }
3759 
3760     /** Enclose large data table in nested static class so it's only loaded on first access. */
3761     private static class ExpFracTable {
3762         /** Exponential over the range of 0 - 1 in increments of 2^-10
3763          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
3764          * 1024 = 2^10
3765          */
3766         private static final double[] EXP_FRAC_TABLE_A;
3767         /** Exponential over the range of 0 - 1 in increments of 2^-10
3768          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
3769          */
3770         private static final double[] EXP_FRAC_TABLE_B;
3771 
3772         static {
3773             if (RECOMPUTE_TABLES_AT_RUNTIME) {
3774                 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
3775                 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
3776 
3777                 final double tmp[] = new double[2];
3778 
3779                 // Populate expFracTable
3780                 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
3781                 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
3782                     FastMathCalc.slowexp(i * factor, tmp);
3783                     EXP_FRAC_TABLE_A[i] = tmp[0];
3784                     EXP_FRAC_TABLE_B[i] = tmp[1];
3785                 }
3786             } else {
3787                 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
3788                 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
3789             }
3790         }
3791     }
3792 
3793     /** Enclose large data table in nested static class so it's only loaded on first access. */
3794     private static class lnMant {
3795         /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
3796         private static final double[][] LN_MANT;
3797 
3798         static {
3799             if (RECOMPUTE_TABLES_AT_RUNTIME) {
3800                 LN_MANT = new double[FastMath.LN_MANT_LEN][];
3801 
3802                 // Populate lnMant table
3803                 for (int i = 0; i < LN_MANT.length; i++) {
3804                     final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
3805                     LN_MANT[i] = FastMathCalc.slowLog(d);
3806                 }
3807             } else {
3808                 LN_MANT = FastMathLiteralArrays.loadLnMant();
3809             }
3810         }
3811     }
3812 
3813     /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
3814     private static class CodyWaite {
3815         /** k */
3816         private final int finalK;
3817         /** remA */
3818         private final double finalRemA;
3819         /** remB */
3820         private final double finalRemB;
3821 
3822         /**
3823          * @param xa Argument.
3824          */
3825         CodyWaite(double xa) {
3826             // Estimate k.
3827             //k = (int)(xa / 1.5707963267948966);
3828             int k = (int)(xa * 0.6366197723675814);
3829 
3830             // Compute remainder.
3831             double remA;
3832             double remB;
3833             while (true) {
3834                 double a = -k * 1.570796251296997;
3835                 remA = xa + a;
3836                 remB = -(remA - xa - a);
3837 
3838                 a = -k * 7.549789948768648E-8;
3839                 double b = remA;
3840                 remA = a + b;
3841                 remB += -(remA - b - a);
3842 
3843                 a = -k * 6.123233995736766E-17;
3844                 b = remA;
3845                 remA = a + b;
3846                 remB += -(remA - b - a);
3847 
3848                 if (remA > 0) {
3849                     break;
3850                 }
3851 
3852                 // Remainder is negative, so decrement k and try again.
3853                 // This should only happen if the input is very close
3854                 // to an even multiple of pi/2.
3855                 --k;
3856             }
3857 
3858             this.finalK = k;
3859             this.finalRemA = remA;
3860             this.finalRemB = remB;
3861         }
3862 
3863         /**
3864          * @return k
3865          */
3866         int getK() {
3867             return finalK;
3868         }
3869         /**
3870          * @return remA
3871          */
3872         double getRemA() {
3873             return finalRemA;
3874         }
3875         /**
3876          * @return remB
3877          */
3878         double getRemB() {
3879             return finalRemB;
3880         }
3881     }
3882 }