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