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         } else if (x != x) { // X is NaN
1466             return x;
1467         } else if (y != y) { // y is NaN
1468             return y;
1469         } else if (x == 0) {
1470             long bits = Double.doubleToRawLongBits(x);
1471             if ((bits & 0x8000000000000000L) != 0) {
1472                 // -zero
1473                 long yi = (long) y;
1474 
1475                 if (y < 0 && y == yi && (yi & 1) == 1) {
1476                     return Double.NEGATIVE_INFINITY;
1477                 }
1478 
1479                 if (y > 0 && y == yi && (yi & 1) == 1) {
1480                     return -0.0;
1481                 }
1482             }
1483 
1484             if (y < 0) {
1485                 return Double.POSITIVE_INFINITY;
1486             }
1487             if (y > 0) {
1488                 return 0.0;
1489             }
1490 
1491             return Double.NaN;
1492         } else if (x == Double.POSITIVE_INFINITY) {
1493             if (y < 0.0) {
1494                 return 0.0;
1495             } else {
1496                 return Double.POSITIVE_INFINITY;
1497             }
1498         } else if (y == Double.POSITIVE_INFINITY) {
1499             if (x * x == 1.0) {
1500                 return Double.NaN;
1501             }
1502 
1503             if (x * x > 1.0) {
1504                 return Double.POSITIVE_INFINITY;
1505             } else {
1506                 return 0.0;
1507             }
1508         } else if (x == Double.NEGATIVE_INFINITY) {
1509             if (y < 0) {
1510                 long yi = (long) y;
1511                 if (y == yi && (yi & 1) == 1) {
1512                     return -0.0;
1513                 }
1514 
1515                 return 0.0;
1516             }
1517 
1518             if (y > 0)  {
1519                 long yi = (long) y;
1520                 if (y == yi && (yi & 1) == 1) {
1521                     return Double.NEGATIVE_INFINITY;
1522                 }
1523 
1524                 return Double.POSITIVE_INFINITY;
1525             }
1526         } else if (y == Double.NEGATIVE_INFINITY) {
1527             if (x * x == 1.0) {
1528                 return Double.NaN;
1529             }
1530 
1531             if (x * x < 1.0) {
1532                 return Double.POSITIVE_INFINITY;
1533             } else {
1534                 return 0.0;
1535             }
1536         } else if (x < 0) { // Handle special case x<0
1537             // y is an even integer in this case
1538             if (y >= TWO_POWER_53 || y <= -TWO_POWER_53) {
1539                 return pow(-x, y);
1540             }
1541 
1542             if (y == (long) y) {
1543                 // If y is an integer
1544                 return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
1545             } else {
1546                 return Double.NaN;
1547             }
1548         }
1549 
1550         /* Split y into ya and yb such that y = ya+yb */
1551         double ya;
1552         double yb;
1553         if (y < 8e298 && y > -8e298) {
1554             double tmp1 = y * HEX_40000000;
1555             ya = y + tmp1 - tmp1;
1556             yb = y - ya;
1557         } else {
1558             double tmp1 = y * 9.31322574615478515625E-10;
1559             double tmp2 = tmp1 * 9.31322574615478515625E-10;
1560             ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
1561             yb = y - ya;
1562         }
1563 
1564         /* Compute ln(x) */
1565         final double lores = log(x, lns);
1566         if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1567             return lores;
1568         }
1569 
1570         double lna = lns[0];
1571         double lnb = lns[1];
1572 
1573         /* resplit lns */
1574         double tmp1 = lna * HEX_40000000;
1575         double tmp2 = lna + tmp1 - tmp1;
1576         lnb += lna - tmp2;
1577         lna = tmp2;
1578 
1579         // y*ln(x) = (aa+ab)
1580         final double aa = lna * ya;
1581         final double ab = lna * yb + lnb * ya + lnb * yb;
1582 
1583         lna = aa+ab;
1584         lnb = -(lna - aa - ab);
1585 
1586         double z = 1.0 / 120.0;
1587         z = z * lnb + (1.0 / 24.0);
1588         z = z * lnb + (1.0 / 6.0);
1589         z = z * lnb + 0.5;
1590         z = z * lnb + 1.0;
1591         z *= lnb;
1592 
1593         final double result = exp(lna, z, null);
1594         //result = result + result * z;
1595         return result;
1596     }
1597 
1598 
1599     /**
1600      * Raise a double to an int power.
1601      *
1602      * @param d Number to raise.
1603      * @param e Exponent.
1604      * @return d<sup>e</sup>
1605      * @since 3.1
1606      */
1607     public static double pow(double d, int e) {
1608 
1609         if (e == 0) {
1610             return 1.0;
1611         } else if (e < 0) {
1612             e = -e;
1613             d = 1.0 / d;
1614         }
1615 
1616         // split d as one 26 bits number and one 27 bits number
1617         // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1618         final double d1High = Double.longBitsToDouble(Double.doubleToRawLongBits(d) & ((-1L) << 27));
1619         final double d1Low  = d - d1High;
1620 
1621         // prepare result
1622         double resultHigh = 1;
1623         double resultLow  = 0;
1624 
1625         // d^(2p)
1626         double d2p     = d;
1627         double d2pHigh = d1High;
1628         double d2pLow  = d1Low;
1629 
1630         while (e != 0) {
1631 
1632             if ((e & 0x1) != 0) {
1633                 // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
1634                 // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1635                 final double tmpHigh = resultHigh * d2p;
1636                 final double rHH     = Double.longBitsToDouble(Double.doubleToRawLongBits(resultHigh) & ((-1L) << 27));
1637                 final double rHL     = resultHigh - rHH;
1638                 final double tmpLow  = rHL * d2pLow - (((tmpHigh - rHH * d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
1639                 resultHigh = tmpHigh;
1640                 resultLow  = resultLow * d2p + tmpLow;
1641             }
1642 
1643             // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
1644             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1645             final double tmpHigh = d2pHigh * d2p;
1646             final double cD2pH   = Double.longBitsToDouble(Double.doubleToRawLongBits(d2pHigh) & ((-1L) << 27));
1647             final double d2pHH   = cD2pH - (cD2pH - d2pHigh);
1648             final double d2pHL   = d2pHigh - d2pHH;
1649             final double tmpLow  = d2pHL * d2pLow - (((tmpHigh - d2pHH * d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
1650             d2pHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(tmpHigh) & ((-1L) << 27));
1651             d2pLow  = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
1652             d2p     = d2pHigh + d2pLow;
1653 
1654             e >>= 1;
1655 
1656         }
1657 
1658         final double result = resultHigh + resultLow;
1659 
1660         if (Double.isNaN(result)) {
1661             if (Double.isNaN(d)) {
1662                 return Double.NaN;
1663             } else {
1664                 // some intermediate numbers exceeded capacity,
1665                 // and the low order bits became NaN (because infinity - infinity = NaN)
1666                 return (d < 0 && (e & 0x1) == 1) ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
1667             }
1668         } else {
1669             return result;
1670         }
1671 
1672     }
1673 
1674     /**
1675      *  Computes sin(x) - x, where |x| < 1/16.
1676      *  Use a Remez polynomial approximation.
1677      *  @param x a number smaller than 1/16
1678      *  @return sin(x) - x
1679      */
1680     private static double polySine(final double x)
1681     {
1682         double x2 = x*x;
1683 
1684         double p = 2.7553817452272217E-6;
1685         p = p * x2 + -1.9841269659586505E-4;
1686         p = p * x2 + 0.008333333333329196;
1687         p = p * x2 + -0.16666666666666666;
1688         //p *= x2;
1689         //p *= x;
1690         p = p * x2 * x;
1691 
1692         return p;
1693     }
1694 
1695     /**
1696      *  Computes cos(x) - 1, where |x| < 1/16.
1697      *  Use a Remez polynomial approximation.
1698      *  @param x a number smaller than 1/16
1699      *  @return cos(x) - 1
1700      */
1701     private static double polyCosine(double x) {
1702         double x2 = x*x;
1703 
1704         double p = 2.479773539153719E-5;
1705         p = p * x2 + -0.0013888888689039883;
1706         p = p * x2 + 0.041666666666621166;
1707         p = p * x2 + -0.49999999999999994;
1708         p *= x2;
1709 
1710         return p;
1711     }
1712 
1713     /**
1714      *  Compute sine over the first quadrant (0 < x < pi/2).
1715      *  Use combination of table lookup and rational polynomial expansion.
1716      *  @param xa number from which sine is requested
1717      *  @param xb extra bits for x (may be 0.0)
1718      *  @return sin(xa + xb)
1719      */
1720     private static double sinQ(double xa, double xb) {
1721         int idx = (int) ((xa * 8.0) + 0.5);
1722         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
1723 
1724         // Table lookups
1725         final double sintA = SINE_TABLE_A[idx];
1726         final double sintB = SINE_TABLE_B[idx];
1727         final double costA = COSINE_TABLE_A[idx];
1728         final double costB = COSINE_TABLE_B[idx];
1729 
1730         // Polynomial eval of sin(epsilon), cos(epsilon)
1731         double sinEpsA = epsilon;
1732         double sinEpsB = polySine(epsilon);
1733         final double cosEpsA = 1.0;
1734         final double cosEpsB = polyCosine(epsilon);
1735 
1736         // Split epsilon   xa + xb = x
1737         final double temp = sinEpsA * HEX_40000000;
1738         double temp2 = (sinEpsA + temp) - temp;
1739         sinEpsB +=  sinEpsA - temp2;
1740         sinEpsA = temp2;
1741 
1742         /* Compute sin(x) by angle addition formula */
1743         double result;
1744 
1745         /* Compute the following sum:
1746          *
1747          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1748          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1749          *
1750          * Ranges of elements
1751          *
1752          * xxxtA   0            PI/2
1753          * xxxtB   -1.5e-9      1.5e-9
1754          * sinEpsA -0.0625      0.0625
1755          * sinEpsB -6e-11       6e-11
1756          * cosEpsA  1.0
1757          * cosEpsB  0           -0.0625
1758          *
1759          */
1760 
1761         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1762         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1763 
1764         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1765         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1766         double a = 0;
1767         double b = 0;
1768 
1769         double t = sintA;
1770         double c = a + t;
1771         double d = -(c - a - t);
1772         a = c;
1773         b += d;
1774 
1775         t = costA * sinEpsA;
1776         c = a + t;
1777         d = -(c - a - t);
1778         a = c;
1779         b += d;
1780 
1781         b = b + sintA * cosEpsB + costA * sinEpsB;
1782         /*
1783     t = sintA*cosEpsB;
1784     c = a + t;
1785     d = -(c - a - t);
1786     a = c;
1787     b = b + d;
1788 
1789     t = costA*sinEpsB;
1790     c = a + t;
1791     d = -(c - a - t);
1792     a = c;
1793     b = b + d;
1794          */
1795 
1796         b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
1797         /*
1798     t = sintB;
1799     c = a + t;
1800     d = -(c - a - t);
1801     a = c;
1802     b = b + d;
1803 
1804     t = costB*sinEpsA;
1805     c = a + t;
1806     d = -(c - a - t);
1807     a = c;
1808     b = b + d;
1809 
1810     t = sintB*cosEpsB;
1811     c = a + t;
1812     d = -(c - a - t);
1813     a = c;
1814     b = b + d;
1815 
1816     t = costB*sinEpsB;
1817     c = a + t;
1818     d = -(c - a - t);
1819     a = c;
1820     b = b + d;
1821          */
1822 
1823         if (xb != 0.0) {
1824             t = ((costA + costB) * (cosEpsA + cosEpsB) -
1825                  (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;  // approximate cosine*xb
1826             c = a + t;
1827             d = -(c - a - t);
1828             a = c;
1829             b += d;
1830         }
1831 
1832         result = a + b;
1833 
1834         return result;
1835     }
1836 
1837     /**
1838      * Compute cosine in the first quadrant by subtracting input from PI/2 and
1839      * then calling sinQ.  This is more accurate as the input approaches PI/2.
1840      *  @param xa number from which cosine is requested
1841      *  @param xb extra bits for x (may be 0.0)
1842      *  @return cos(xa + xb)
1843      */
1844     private static double cosQ(double xa, double xb) {
1845         final double pi2a = 1.5707963267948966;
1846         final double pi2b = 6.123233995736766E-17;
1847 
1848         final double a = pi2a - xa;
1849         double b = -(a - pi2a + xa);
1850         b += pi2b - xb;
1851 
1852         return sinQ(a, b);
1853     }
1854 
1855     /**
1856      *  Compute tangent (or cotangent) over the first quadrant.   0 < x < pi/2
1857      *  Use combination of table lookup and rational polynomial expansion.
1858      *  @param xa number from which sine is requested
1859      *  @param xb extra bits for x (may be 0.0)
1860      *  @param cotanFlag if true, compute the cotangent instead of the tangent
1861      *  @return tan(xa+xb) (or cotangent, depending on cotanFlag)
1862      */
1863     private static double tanQ(double xa, double xb, boolean cotanFlag) {
1864 
1865         int idx = (int) ((xa * 8.0) + 0.5);
1866         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
1867 
1868         // Table lookups
1869         final double sintA = SINE_TABLE_A[idx];
1870         final double sintB = SINE_TABLE_B[idx];
1871         final double costA = COSINE_TABLE_A[idx];
1872         final double costB = COSINE_TABLE_B[idx];
1873 
1874         // Polynomial eval of sin(epsilon), cos(epsilon)
1875         double sinEpsA = epsilon;
1876         double sinEpsB = polySine(epsilon);
1877         final double cosEpsA = 1.0;
1878         final double cosEpsB = polyCosine(epsilon);
1879 
1880         // Split epsilon   xa + xb = x
1881         double temp = sinEpsA * HEX_40000000;
1882         double temp2 = (sinEpsA + temp) - temp;
1883         sinEpsB +=  sinEpsA - temp2;
1884         sinEpsA = temp2;
1885 
1886         /* Compute sin(x) by angle addition formula */
1887 
1888         /* Compute the following sum:
1889          *
1890          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1891          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1892          *
1893          * Ranges of elements
1894          *
1895          * xxxtA   0            PI/2
1896          * xxxtB   -1.5e-9      1.5e-9
1897          * sinEpsA -0.0625      0.0625
1898          * sinEpsB -6e-11       6e-11
1899          * cosEpsA  1.0
1900          * cosEpsB  0           -0.0625
1901          *
1902          */
1903 
1904         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
1905         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1906 
1907         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
1908         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
1909         double a = 0;
1910         double b = 0;
1911 
1912         // Compute sine
1913         double t = sintA;
1914         double c = a + t;
1915         double d = -(c - a - t);
1916         a = c;
1917         b += d;
1918 
1919         t = costA*sinEpsA;
1920         c = a + t;
1921         d = -(c - a - t);
1922         a = c;
1923         b += d;
1924 
1925         b += sintA*cosEpsB + costA*sinEpsB;
1926         b += sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
1927 
1928         double sina = a + b;
1929         double sinb = -(sina - a - b);
1930 
1931         // Compute cosine
1932 
1933         a = b = c = d = 0.0;
1934 
1935         t = costA*cosEpsA;
1936         c = a + t;
1937         d = -(c - a - t);
1938         a = c;
1939         b += d;
1940 
1941         t = -sintA*sinEpsA;
1942         c = a + t;
1943         d = -(c - a - t);
1944         a = c;
1945         b += d;
1946 
1947         b += costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
1948         b -= sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB;
1949 
1950         double cosa = a + b;
1951         double cosb = -(cosa - a - b);
1952 
1953         if (cotanFlag) {
1954             double tmp;
1955             tmp = cosa; cosa = sina; sina = tmp;
1956             tmp = cosb; cosb = sinb; sinb = tmp;
1957         }
1958 
1959 
1960         /* estimate and correct, compute 1.0/(cosa+cosb) */
1961         /*
1962     double est = (sina+sinb)/(cosa+cosb);
1963     double err = (sina - cosa*est) + (sinb - cosb*est);
1964     est += err/(cosa+cosb);
1965     err = (sina - cosa*est) + (sinb - cosb*est);
1966          */
1967 
1968         // f(x) = 1/x,   f'(x) = -1/x^2
1969 
1970         double est = sina/cosa;
1971 
1972         /* Split the estimate to get more accurate read on division rounding */
1973         temp = est * HEX_40000000;
1974         double esta = (est + temp) - temp;
1975         double estb =  est - esta;
1976 
1977         temp = cosa * HEX_40000000;
1978         double cosaa = (cosa + temp) - temp;
1979         double cosab =  cosa - cosaa;
1980 
1981         //double err = (sina - est*cosa)/cosa;  // Correction for division rounding
1982         double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa;  // Correction for division rounding
1983         err += sinb/cosa;                     // Change in est due to sinb
1984         err += -sina * cosb / cosa / cosa;    // Change in est due to cosb
1985 
1986         if (xb != 0.0) {
1987             // tan' = 1 + tan^2      cot' = -(1 + cot^2)
1988             // Approximate impact of xb
1989             double xbadj = xb + est*est*xb;
1990             if (cotanFlag) {
1991                 xbadj = -xbadj;
1992             }
1993 
1994             err += xbadj;
1995         }
1996 
1997         return est+err;
1998     }
1999 
2000     /** Reduce the input argument using the Payne and Hanek method.
2001      *  This is good for all inputs 0.0 < x < inf
2002      *  Output is remainder after dividing by PI/2
2003      *  The result array should contain 3 numbers.
2004      *  result[0] is the integer portion, so mod 4 this gives the quadrant.
2005      *  result[1] is the upper bits of the remainder
2006      *  result[2] is the lower bits of the remainder
2007      *
2008      * @param x number to reduce
2009      * @param result placeholder where to put the result
2010      */
2011     private static void reducePayneHanek(double x, double result[])
2012     {
2013         /* Convert input double to bits */
2014         long inbits = Double.doubleToRawLongBits(x);
2015         int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2016 
2017         /* Convert to fixed point representation */
2018         inbits &= 0x000fffffffffffffL;
2019         inbits |= 0x0010000000000000L;
2020 
2021         /* Normalize input to be between 0.5 and 1.0 */
2022         exponent++;
2023         inbits <<= 11;
2024 
2025         /* Based on the exponent, get a shifted copy of recip2pi */
2026         long shpi0;
2027         long shpiA;
2028         long shpiB;
2029         int idx = exponent >> 6;
2030         int shift = exponent - (idx << 6);
2031 
2032         if (shift != 0) {
2033             shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
2034             shpi0 |= RECIP_2PI[idx] >>> (64-shift);
2035             shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
2036             shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
2037         } else {
2038             shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
2039             shpiA = RECIP_2PI[idx];
2040             shpiB = RECIP_2PI[idx+1];
2041         }
2042 
2043         /* Multiply input by shpiA */
2044         long a = inbits >>> 32;
2045         long b = inbits & 0xffffffffL;
2046 
2047         long c = shpiA >>> 32;
2048         long d = shpiA & 0xffffffffL;
2049 
2050         long ac = a * c;
2051         long bd = b * d;
2052         long bc = b * c;
2053         long ad = a * d;
2054 
2055         long prodB = bd + (ad << 32);
2056         long prodA = ac + (ad >>> 32);
2057 
2058         boolean bita = (bd & 0x8000000000000000L) != 0;
2059         boolean bitb = (ad & 0x80000000L ) != 0;
2060         boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2061 
2062         /* Carry */
2063         if ( (bita && bitb) ||
2064                 ((bita || bitb) && !bitsum) ) {
2065             prodA++;
2066         }
2067 
2068         bita = (prodB & 0x8000000000000000L) != 0;
2069         bitb = (bc & 0x80000000L ) != 0;
2070 
2071         prodB += bc << 32;
2072         prodA += bc >>> 32;
2073 
2074         bitsum = (prodB & 0x8000000000000000L) != 0;
2075 
2076         /* Carry */
2077         if ( (bita && bitb) ||
2078                 ((bita || bitb) && !bitsum) ) {
2079             prodA++;
2080         }
2081 
2082         /* Multiply input by shpiB */
2083         c = shpiB >>> 32;
2084         d = shpiB & 0xffffffffL;
2085         ac = a * c;
2086         bc = b * c;
2087         ad = a * d;
2088 
2089         /* Collect terms */
2090         ac += (bc + ad) >>> 32;
2091 
2092         bita = (prodB & 0x8000000000000000L) != 0;
2093         bitb = (ac & 0x8000000000000000L ) != 0;
2094         prodB += ac;
2095         bitsum = (prodB & 0x8000000000000000L) != 0;
2096         /* Carry */
2097         if ( (bita && bitb) ||
2098                 ((bita || bitb) && !bitsum) ) {
2099             prodA++;
2100         }
2101 
2102         /* Multiply by shpi0 */
2103         c = shpi0 >>> 32;
2104         d = shpi0 & 0xffffffffL;
2105 
2106         bd = b * d;
2107         bc = b * c;
2108         ad = a * d;
2109 
2110         prodA += bd + ((bc + ad) << 32);
2111 
2112         /*
2113          * prodA, prodB now contain the remainder as a fraction of PI.  We want this as a fraction of
2114          * PI/2, so use the following steps:
2115          * 1.) multiply by 4.
2116          * 2.) do a fixed point muliply by PI/4.
2117          * 3.) Convert to floating point.
2118          * 4.) Multiply by 2
2119          */
2120 
2121         /* This identifies the quadrant */
2122         int intPart = (int)(prodA >>> 62);
2123 
2124         /* Multiply by 4 */
2125         prodA <<= 2;
2126         prodA |= prodB >>> 62;
2127         prodB <<= 2;
2128 
2129         /* Multiply by PI/4 */
2130         a = prodA >>> 32;
2131         b = prodA & 0xffffffffL;
2132 
2133         c = PI_O_4_BITS[0] >>> 32;
2134         d = PI_O_4_BITS[0] & 0xffffffffL;
2135 
2136         ac = a * c;
2137         bd = b * d;
2138         bc = b * c;
2139         ad = a * d;
2140 
2141         long prod2B = bd + (ad << 32);
2142         long prod2A = ac + (ad >>> 32);
2143 
2144         bita = (bd & 0x8000000000000000L) != 0;
2145         bitb = (ad & 0x80000000L ) != 0;
2146         bitsum = (prod2B & 0x8000000000000000L) != 0;
2147 
2148         /* Carry */
2149         if ( (bita && bitb) ||
2150                 ((bita || bitb) && !bitsum) ) {
2151             prod2A++;
2152         }
2153 
2154         bita = (prod2B & 0x8000000000000000L) != 0;
2155         bitb = (bc & 0x80000000L ) != 0;
2156 
2157         prod2B += bc << 32;
2158         prod2A += bc >>> 32;
2159 
2160         bitsum = (prod2B & 0x8000000000000000L) != 0;
2161 
2162         /* Carry */
2163         if ( (bita && bitb) ||
2164                 ((bita || bitb) && !bitsum) ) {
2165             prod2A++;
2166         }
2167 
2168         /* Multiply input by pio4bits[1] */
2169         c = PI_O_4_BITS[1] >>> 32;
2170         d = PI_O_4_BITS[1] & 0xffffffffL;
2171         ac = a * c;
2172         bc = b * c;
2173         ad = a * d;
2174 
2175         /* Collect terms */
2176         ac += (bc + ad) >>> 32;
2177 
2178         bita = (prod2B & 0x8000000000000000L) != 0;
2179         bitb = (ac & 0x8000000000000000L ) != 0;
2180         prod2B += ac;
2181         bitsum = (prod2B & 0x8000000000000000L) != 0;
2182         /* Carry */
2183         if ( (bita && bitb) ||
2184                 ((bita || bitb) && !bitsum) ) {
2185             prod2A++;
2186         }
2187 
2188         /* Multiply inputB by pio4bits[0] */
2189         a = prodB >>> 32;
2190         b = prodB & 0xffffffffL;
2191         c = PI_O_4_BITS[0] >>> 32;
2192         d = PI_O_4_BITS[0] & 0xffffffffL;
2193         ac = a * c;
2194         bc = b * c;
2195         ad = a * d;
2196 
2197         /* Collect terms */
2198         ac += (bc + ad) >>> 32;
2199 
2200         bita = (prod2B & 0x8000000000000000L) != 0;
2201         bitb = (ac & 0x8000000000000000L ) != 0;
2202         prod2B += ac;
2203         bitsum = (prod2B & 0x8000000000000000L) != 0;
2204         /* Carry */
2205         if ( (bita && bitb) ||
2206                 ((bita || bitb) && !bitsum) ) {
2207             prod2A++;
2208         }
2209 
2210         /* Convert to double */
2211         double tmpA = (prod2A >>> 12) / TWO_POWER_52;  // High order 52 bits
2212         double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
2213 
2214         double sumA = tmpA + tmpB;
2215         double sumB = -(sumA - tmpA - tmpB);
2216 
2217         /* Multiply by PI/2 and return */
2218         result[0] = intPart;
2219         result[1] = sumA * 2.0;
2220         result[2] = sumB * 2.0;
2221     }
2222 
2223     /**
2224      * Sine function.
2225      *
2226      * @param x Argument.
2227      * @return sin(x)
2228      */
2229     public static double sin(double x) {
2230         boolean negative = false;
2231         int quadrant = 0;
2232         double xa;
2233         double xb = 0.0;
2234 
2235         /* Take absolute value of the input */
2236         xa = x;
2237         if (x < 0) {
2238             negative = true;
2239             xa = -xa;
2240         }
2241 
2242         /* Check for zero and negative zero */
2243         if (xa == 0.0) {
2244             long bits = Double.doubleToRawLongBits(x);
2245             if (bits < 0) {
2246                 return -0.0;
2247             }
2248             return 0.0;
2249         }
2250 
2251         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2252             return Double.NaN;
2253         }
2254 
2255         /* Perform any argument reduction */
2256         if (xa > 3294198.0) {
2257             // PI * (2**20)
2258             // Argument too big for CodyWaite reduction.  Must use
2259             // PayneHanek.
2260             double reduceResults[] = new double[3];
2261             reducePayneHanek(xa, reduceResults);
2262             quadrant = ((int) reduceResults[0]) & 3;
2263             xa = reduceResults[1];
2264             xb = reduceResults[2];
2265         } else if (xa > 1.5707963267948966) {
2266             final CodyWaite cw = new CodyWaite(xa);
2267             quadrant = cw.getK() & 3;
2268             xa = cw.getRemA();
2269             xb = cw.getRemB();
2270         }
2271 
2272         if (negative) {
2273             quadrant ^= 2;  // Flip bit 1
2274         }
2275 
2276         switch (quadrant) {
2277             case 0:
2278                 return sinQ(xa, xb);
2279             case 1:
2280                 return cosQ(xa, xb);
2281             case 2:
2282                 return -sinQ(xa, xb);
2283             case 3:
2284                 return -cosQ(xa, xb);
2285             default:
2286                 return Double.NaN;
2287         }
2288     }
2289 
2290     /**
2291      * Cosine function.
2292      *
2293      * @param x Argument.
2294      * @return cos(x)
2295      */
2296     public static double cos(double x) {
2297         int quadrant = 0;
2298 
2299         /* Take absolute value of the input */
2300         double xa = x;
2301         if (x < 0) {
2302             xa = -xa;
2303         }
2304 
2305         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2306             return Double.NaN;
2307         }
2308 
2309         /* Perform any argument reduction */
2310         double xb = 0;
2311         if (xa > 3294198.0) {
2312             // PI * (2**20)
2313             // Argument too big for CodyWaite reduction.  Must use
2314             // PayneHanek.
2315             double reduceResults[] = new double[3];
2316             reducePayneHanek(xa, reduceResults);
2317             quadrant = ((int) reduceResults[0]) & 3;
2318             xa = reduceResults[1];
2319             xb = reduceResults[2];
2320         } else if (xa > 1.5707963267948966) {
2321             final CodyWaite cw = new CodyWaite(xa);
2322             quadrant = cw.getK() & 3;
2323             xa = cw.getRemA();
2324             xb = cw.getRemB();
2325         }
2326 
2327         //if (negative)
2328         //  quadrant = (quadrant + 2) % 4;
2329 
2330         switch (quadrant) {
2331             case 0:
2332                 return cosQ(xa, xb);
2333             case 1:
2334                 return -sinQ(xa, xb);
2335             case 2:
2336                 return -cosQ(xa, xb);
2337             case 3:
2338                 return sinQ(xa, xb);
2339             default:
2340                 return Double.NaN;
2341         }
2342     }
2343 
2344     /**
2345      * Tangent function.
2346      *
2347      * @param x Argument.
2348      * @return tan(x)
2349      */
2350     public static double tan(double x) {
2351         boolean negative = false;
2352         int quadrant = 0;
2353 
2354         /* Take absolute value of the input */
2355         double xa = x;
2356         if (x < 0) {
2357             negative = true;
2358             xa = -xa;
2359         }
2360 
2361         /* Check for zero and negative zero */
2362         if (xa == 0.0) {
2363             long bits = Double.doubleToRawLongBits(x);
2364             if (bits < 0) {
2365                 return -0.0;
2366             }
2367             return 0.0;
2368         }
2369 
2370         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2371             return Double.NaN;
2372         }
2373 
2374         /* Perform any argument reduction */
2375         double xb = 0;
2376         if (xa > 3294198.0) {
2377             // PI * (2**20)
2378             // Argument too big for CodyWaite reduction.  Must use
2379             // PayneHanek.
2380             double reduceResults[] = new double[3];
2381             reducePayneHanek(xa, reduceResults);
2382             quadrant = ((int) reduceResults[0]) & 3;
2383             xa = reduceResults[1];
2384             xb = reduceResults[2];
2385         } else if (xa > 1.5707963267948966) {
2386             final CodyWaite cw = new CodyWaite(xa);
2387             quadrant = cw.getK() & 3;
2388             xa = cw.getRemA();
2389             xb = cw.getRemB();
2390         }
2391 
2392         if (xa > 1.5) {
2393             // Accuracy suffers between 1.5 and PI/2
2394             final double pi2a = 1.5707963267948966;
2395             final double pi2b = 6.123233995736766E-17;
2396 
2397             final double a = pi2a - xa;
2398             double b = -(a - pi2a + xa);
2399             b += pi2b - xb;
2400 
2401             xa = a + b;
2402             xb = -(xa - a - b);
2403             quadrant ^= 1;
2404             negative ^= true;
2405         }
2406 
2407         double result;
2408         if ((quadrant & 1) == 0) {
2409             result = tanQ(xa, xb, false);
2410         } else {
2411             result = -tanQ(xa, xb, true);
2412         }
2413 
2414         if (negative) {
2415             result = -result;
2416         }
2417 
2418         return result;
2419     }
2420 
2421     /**
2422      * Arctangent function
2423      *  @param x a number
2424      *  @return atan(x)
2425      */
2426     public static double atan(double x) {
2427         return atan(x, 0.0, false);
2428     }
2429 
2430     /** Internal helper function to compute arctangent.
2431      * @param xa number from which arctangent is requested
2432      * @param xb extra bits for x (may be 0.0)
2433      * @param leftPlane if true, result angle must be put in the left half plane
2434      * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
2435      */
2436     private static double atan(double xa, double xb, boolean leftPlane) {
2437         if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2438             return leftPlane ? copySign(Math.PI, xa) : xa;
2439         }
2440 
2441         final boolean negate;
2442         if (xa < 0) {
2443             // negative
2444             xa = -xa;
2445             xb = -xb;
2446             negate = true;
2447         } else {
2448             negate = false;
2449         }
2450 
2451         if (xa > 1.633123935319537E16) { // Very large input
2452             return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2453         }
2454 
2455         /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
2456         final int idx;
2457         if (xa < 1) {
2458             idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2459         } else {
2460             final double oneOverXa = 1 / xa;
2461             idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2462         }
2463 
2464         final double ttA = TANGENT_TABLE_A[idx];
2465         final double ttB = TANGENT_TABLE_B[idx];
2466 
2467         double epsA = xa - ttA;
2468         double epsB = -(epsA - xa + ttA);
2469         epsB += xb - ttB;
2470 
2471         double temp = epsA + epsB;
2472         epsB = -(temp - epsA - epsB);
2473         epsA = temp;
2474 
2475         /* Compute eps = eps / (1.0 + xa*tangent) */
2476         temp = xa * HEX_40000000;
2477         double ya = xa + temp - temp;
2478         double yb = xb + xa - ya;
2479         xa = ya;
2480         xb += yb;
2481 
2482         //if (idx > 8 || idx == 0)
2483         if (idx == 0) {
2484             /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2485             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2486             final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2487             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
2488             ya = epsA * denom;
2489             yb = epsB * denom;
2490         } else {
2491             double temp2 = xa * ttA;
2492             double za = 1d + temp2;
2493             double zb = -(za - 1d - temp2);
2494             temp2 = xb * ttA + xa * ttB;
2495             temp = za + temp2;
2496             zb += -(temp - za - temp2);
2497             za = temp;
2498 
2499             zb += xb * ttB;
2500             ya = epsA / za;
2501 
2502             temp = ya * HEX_40000000;
2503             final double yaa = (ya + temp) - temp;
2504             final double yab = ya - yaa;
2505 
2506             temp = za * HEX_40000000;
2507             final double zaa = (za + temp) - temp;
2508             final double zab = za - zaa;
2509 
2510             /* Correct for rounding in division */
2511             yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2512 
2513             yb += -epsA * zb / za / za;
2514             yb += epsB / za;
2515         }
2516 
2517 
2518         epsA = ya;
2519         epsB = yb;
2520 
2521         /* Evaluate polynomial */
2522         final double epsA2 = epsA * epsA;
2523 
2524         /*
2525     yb = -0.09001346640161823;
2526     yb = yb * epsA2 + 0.11110718400605211;
2527     yb = yb * epsA2 + -0.1428571349122913;
2528     yb = yb * epsA2 + 0.19999999999273194;
2529     yb = yb * epsA2 + -0.33333333333333093;
2530     yb = yb * epsA2 * epsA;
2531          */
2532 
2533         yb = 0.07490822288864472;
2534         yb = yb * epsA2 - 0.09088450866185192;
2535         yb = yb * epsA2 + 0.11111095942313305;
2536         yb = yb * epsA2 - 0.1428571423679182;
2537         yb = yb * epsA2 + 0.19999999999923582;
2538         yb = yb * epsA2 - 0.33333333333333287;
2539         yb = yb * epsA2 * epsA;
2540 
2541 
2542         ya = epsA;
2543 
2544         temp = ya + yb;
2545         yb = -(temp - ya - yb);
2546         ya = temp;
2547 
2548         /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
2549         yb += epsB / (1d + epsA * epsA);
2550 
2551         final double eighths = EIGHTHS[idx];
2552 
2553         //result = yb + eighths[idx] + ya;
2554         double za = eighths + ya;
2555         double zb = -(za - eighths - ya);
2556         temp = za + yb;
2557         zb += -(temp - za - yb);
2558         za = temp;
2559 
2560         double result = za + zb;
2561 
2562         if (leftPlane) {
2563             // Result is in the left plane
2564             final double resultb = -(result - za - zb);
2565             final double pia = 1.5707963267948966 * 2;
2566             final double pib = 6.123233995736766E-17 * 2;
2567 
2568             za = pia - result;
2569             zb = -(za - pia + result);
2570             zb += pib - resultb;
2571 
2572             result = za + zb;
2573         }
2574 
2575 
2576         if (negate ^ leftPlane) {
2577             result = -result;
2578         }
2579 
2580         return result;
2581     }
2582 
2583     /**
2584      * Two arguments arctangent function
2585      * @param y ordinate
2586      * @param x abscissa
2587      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
2588      */
2589     public static double atan2(double y, double x) {
2590         if (x != x || y != y) {
2591             return Double.NaN;
2592         }
2593 
2594         if (y == 0) {
2595             final double result = x * y;
2596             final double invx = 1d / x;
2597             final double invy = 1d / y;
2598 
2599             if (invx == 0) { // X is infinite
2600                 if (x > 0) {
2601                     return y; // return +/- 0.0
2602                 } else {
2603                     return copySign(Math.PI, y);
2604                 }
2605             }
2606 
2607             if (x < 0 || invx < 0) {
2608                 if (y < 0 || invy < 0) {
2609                     return -Math.PI;
2610                 } else {
2611                     return Math.PI;
2612                 }
2613             } else {
2614                 return result;
2615             }
2616         }
2617 
2618         // y cannot now be zero
2619 
2620         if (y == Double.POSITIVE_INFINITY) {
2621             if (x == Double.POSITIVE_INFINITY) {
2622                 return Math.PI * F_1_4;
2623             }
2624 
2625             if (x == Double.NEGATIVE_INFINITY) {
2626                 return Math.PI * F_3_4;
2627             }
2628 
2629             return Math.PI * F_1_2;
2630         }
2631 
2632         if (y == Double.NEGATIVE_INFINITY) {
2633             if (x == Double.POSITIVE_INFINITY) {
2634                 return -Math.PI * F_1_4;
2635             }
2636 
2637             if (x == Double.NEGATIVE_INFINITY) {
2638                 return -Math.PI * F_3_4;
2639             }
2640 
2641             return -Math.PI * F_1_2;
2642         }
2643 
2644         if (x == Double.POSITIVE_INFINITY) {
2645             if (y > 0 || 1 / y > 0) {
2646                 return 0d;
2647             }
2648 
2649             if (y < 0 || 1 / y < 0) {
2650                 return -0d;
2651             }
2652         }
2653 
2654         if (x == Double.NEGATIVE_INFINITY)
2655         {
2656             if (y > 0.0 || 1 / y > 0.0) {
2657                 return Math.PI;
2658             }
2659 
2660             if (y < 0 || 1 / y < 0) {
2661                 return -Math.PI;
2662             }
2663         }
2664 
2665         // Neither y nor x can be infinite or NAN here
2666 
2667         if (x == 0) {
2668             if (y > 0 || 1 / y > 0) {
2669                 return Math.PI * F_1_2;
2670             }
2671 
2672             if (y < 0 || 1 / y < 0) {
2673                 return -Math.PI * F_1_2;
2674             }
2675         }
2676 
2677         // Compute ratio r = y/x
2678         final double r = y / x;
2679         if (Double.isInfinite(r)) { // bypass calculations that can create NaN
2680             return atan(r, 0, x < 0);
2681         }
2682 
2683         double ra = doubleHighPart(r);
2684         double rb = r - ra;
2685 
2686         // Split x
2687         final double xa = doubleHighPart(x);
2688         final double xb = x - xa;
2689 
2690         rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
2691 
2692         final double temp = ra + rb;
2693         rb = -(temp - ra - rb);
2694         ra = temp;
2695 
2696         if (ra == 0) { // Fix up the sign so atan works correctly
2697             ra = copySign(0d, y);
2698         }
2699 
2700         // Call atan
2701         final double result = atan(ra, rb, x < 0);
2702 
2703         return result;
2704     }
2705 
2706     /** Compute the arc sine of a number.
2707      * @param x number on which evaluation is done
2708      * @return arc sine of x
2709      */
2710     public static double asin(double x) {
2711       if (x != x) {
2712           return Double.NaN;
2713       }
2714 
2715       if (x > 1.0 || x < -1.0) {
2716           return Double.NaN;
2717       }
2718 
2719       if (x == 1.0) {
2720           return Math.PI/2.0;
2721       }
2722 
2723       if (x == -1.0) {
2724           return -Math.PI/2.0;
2725       }
2726 
2727       if (x == 0.0) { // Matches +/- 0.0; return correct sign
2728           return x;
2729       }
2730 
2731       /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
2732 
2733       /* Split x */
2734       double temp = x * HEX_40000000;
2735       final double xa = x + temp - temp;
2736       final double xb = x - xa;
2737 
2738       /* Square it */
2739       double ya = xa*xa;
2740       double yb = xa*xb*2.0 + xb*xb;
2741 
2742       /* Subtract from 1 */
2743       ya = -ya;
2744       yb = -yb;
2745 
2746       double za = 1.0 + ya;
2747       double zb = -(za - 1.0 - ya);
2748 
2749       temp = za + yb;
2750       zb += -(temp - za - yb);
2751       za = temp;
2752 
2753       /* Square root */
2754       double y;
2755       y = sqrt(za);
2756       temp = y * HEX_40000000;
2757       ya = y + temp - temp;
2758       yb = y - ya;
2759 
2760       /* Extend precision of sqrt */
2761       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2762 
2763       /* Contribution of zb to sqrt */
2764       double dx = zb / (2.0*y);
2765 
2766       // Compute ratio r = x/y
2767       double r = x/y;
2768       temp = r * HEX_40000000;
2769       double ra = r + temp - temp;
2770       double rb = r - ra;
2771 
2772       rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y;  // Correct for rounding in division
2773       rb += -x * dx / y / y;  // Add in effect additional bits of sqrt.
2774 
2775       temp = ra + rb;
2776       rb = -(temp - ra - rb);
2777       ra = temp;
2778 
2779       return atan(ra, rb, false);
2780     }
2781 
2782     /** Compute the arc cosine of a number.
2783      * @param x number on which evaluation is done
2784      * @return arc cosine of x
2785      */
2786     public static double acos(double x) {
2787       if (x != x) {
2788           return Double.NaN;
2789       }
2790 
2791       if (x > 1.0 || x < -1.0) {
2792           return Double.NaN;
2793       }
2794 
2795       if (x == -1.0) {
2796           return Math.PI;
2797       }
2798 
2799       if (x == 1.0) {
2800           return 0.0;
2801       }
2802 
2803       if (x == 0) {
2804           return Math.PI/2.0;
2805       }
2806 
2807       /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
2808 
2809       /* Split x */
2810       double temp = x * HEX_40000000;
2811       final double xa = x + temp - temp;
2812       final double xb = x - xa;
2813 
2814       /* Square it */
2815       double ya = xa*xa;
2816       double yb = xa*xb*2.0 + xb*xb;
2817 
2818       /* Subtract from 1 */
2819       ya = -ya;
2820       yb = -yb;
2821 
2822       double za = 1.0 + ya;
2823       double zb = -(za - 1.0 - ya);
2824 
2825       temp = za + yb;
2826       zb += -(temp - za - yb);
2827       za = temp;
2828 
2829       /* Square root */
2830       double y = sqrt(za);
2831       temp = y * HEX_40000000;
2832       ya = y + temp - temp;
2833       yb = y - ya;
2834 
2835       /* Extend precision of sqrt */
2836       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
2837 
2838       /* Contribution of zb to sqrt */
2839       yb += zb / (2.0*y);
2840       y = ya+yb;
2841       yb = -(y - ya - yb);
2842 
2843       // Compute ratio r = y/x
2844       double r = y/x;
2845 
2846       // Did r overflow?
2847       if (Double.isInfinite(r)) { // x is effectively zero
2848           return Math.PI/2; // so return the appropriate value
2849       }
2850 
2851       double ra = doubleHighPart(r);
2852       double rb = r - ra;
2853 
2854       rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x;  // Correct for rounding in division
2855       rb += yb / x;  // Add in effect additional bits of sqrt.
2856 
2857       temp = ra + rb;
2858       rb = -(temp - ra - rb);
2859       ra = temp;
2860 
2861       return atan(ra, rb, x<0);
2862     }
2863 
2864     /** Compute the cubic root of a number.
2865      * @param x number on which evaluation is done
2866      * @return cubic root of x
2867      */
2868     public static double cbrt(double x) {
2869       /* Convert input double to bits */
2870       long inbits = Double.doubleToRawLongBits(x);
2871       int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2872       boolean subnormal = false;
2873 
2874       if (exponent == -1023) {
2875           if (x == 0) {
2876               return x;
2877           }
2878 
2879           /* Subnormal, so normalize */
2880           subnormal = true;
2881           x *= 1.8014398509481984E16;  // 2^54
2882           inbits = Double.doubleToRawLongBits(x);
2883           exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2884       }
2885 
2886       if (exponent == 1024) {
2887           // Nan or infinity.  Don't care which.
2888           return x;
2889       }
2890 
2891       /* Divide the exponent by 3 */
2892       int exp3 = exponent / 3;
2893 
2894       /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
2895       double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
2896                                           (long)(((exp3 + 1023) & 0x7ff)) << 52);
2897 
2898       /* This will be a number between 1 and 2 */
2899       final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
2900 
2901       /* Estimate the cube root of mant by polynomial */
2902       double est = -0.010714690733195933;
2903       est = est * mant + 0.0875862700108075;
2904       est = est * mant + -0.3058015757857271;
2905       est = est * mant + 0.7249995199969751;
2906       est = est * mant + 0.5039018405998233;
2907 
2908       est *= CBRTTWO[exponent % 3 + 2];
2909 
2910       // est should now be good to about 15 bits of precision.   Do 2 rounds of
2911       // Newton's method to get closer,  this should get us full double precision
2912       // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
2913       final double xs = x / (p2*p2*p2);
2914       est += (xs - est*est*est) / (3*est*est);
2915       est += (xs - est*est*est) / (3*est*est);
2916 
2917       // Do one round of Newton's method in extended precision to get the last bit right.
2918       double temp = est * HEX_40000000;
2919       double ya = est + temp - temp;
2920       double yb = est - ya;
2921 
2922       double za = ya * ya;
2923       double zb = ya * yb * 2.0 + yb * yb;
2924       temp = za * HEX_40000000;
2925       double temp2 = za + temp - temp;
2926       zb += za - temp2;
2927       za = temp2;
2928 
2929       zb = za * yb + ya * zb + zb * yb;
2930       za *= ya;
2931 
2932       double na = xs - za;
2933       double nb = -(na - xs + za);
2934       nb -= zb;
2935 
2936       est += (na+nb)/(3*est*est);
2937 
2938       /* Scale by a power of two, so this is exact. */
2939       est *= p2;
2940 
2941       if (subnormal) {
2942           est *= 3.814697265625E-6;  // 2^-18
2943       }
2944 
2945       return est;
2946     }
2947 
2948     /**
2949      *  Convert degrees to radians, with error of less than 0.5 ULP
2950      *  @param x angle in degrees
2951      *  @return x converted into radians
2952      */
2953     public static double toRadians(double x)
2954     {
2955         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2956             return x;
2957         }
2958 
2959         // These are PI/180 split into high and low order bits
2960         final double facta = 0.01745329052209854;
2961         final double factb = 1.997844754509471E-9;
2962 
2963         double xa = doubleHighPart(x);
2964         double xb = x - xa;
2965 
2966         double result = xb * factb + xb * facta + xa * factb + xa * facta;
2967         if (result == 0) {
2968             result *= x; // ensure correct sign if calculation underflows
2969         }
2970         return result;
2971     }
2972 
2973     /**
2974      *  Convert radians to degrees, with error of less than 0.5 ULP
2975      *  @param x angle in radians
2976      *  @return x converted into degrees
2977      */
2978     public static double toDegrees(double x)
2979     {
2980         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
2981             return x;
2982         }
2983 
2984         // These are 180/PI split into high and low order bits
2985         final double facta = 57.2957763671875;
2986         final double factb = 3.145894820876798E-6;
2987 
2988         double xa = doubleHighPart(x);
2989         double xb = x - xa;
2990 
2991         return xb * factb + xb * facta + xa * factb + xa * facta;
2992     }
2993 
2994     /**
2995      * Absolute value.
2996      * @param x number from which absolute value is requested
2997      * @return abs(x)
2998      */
2999     public static int abs(final int x) {
3000         final int i = x >>> 31;
3001         return (x ^ (~i + 1)) + i;
3002     }
3003 
3004     /**
3005      * Absolute value.
3006      * @param x number from which absolute value is requested
3007      * @return abs(x)
3008      */
3009     public static long abs(final long x) {
3010         final long l = x >>> 63;
3011         // l is one if x negative zero else
3012         // ~l+1 is zero if x is positive, -1 if x is negative
3013         // x^(~l+1) is x is x is positive, ~x if x is negative
3014         // add around
3015         return (x ^ (~l + 1)) + l;
3016     }
3017 
3018     /**
3019      * Absolute value.
3020      * @param x number from which absolute value is requested
3021      * @return abs(x)
3022      */
3023     public static float abs(final float x) {
3024         return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3025     }
3026 
3027     /**
3028      * Absolute value.
3029      * @param x number from which absolute value is requested
3030      * @return abs(x)
3031      */
3032     public static double abs(double x) {
3033         return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3034     }
3035 
3036     /**
3037      * Compute least significant bit (Unit in Last Position) for a number.
3038      * @param x number from which ulp is requested
3039      * @return ulp(x)
3040      */
3041     public static double ulp(double x) {
3042         if (Double.isInfinite(x)) {
3043             return Double.POSITIVE_INFINITY;
3044         }
3045         return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3046     }
3047 
3048     /**
3049      * Compute least significant bit (Unit in Last Position) for a number.
3050      * @param x number from which ulp is requested
3051      * @return ulp(x)
3052      */
3053     public static float ulp(float x) {
3054         if (Float.isInfinite(x)) {
3055             return Float.POSITIVE_INFINITY;
3056         }
3057         return abs(x - Float.intBitsToFloat(Float.floatToIntBits(x) ^ 1));
3058     }
3059 
3060     /**
3061      * Multiply a double number by a power of 2.
3062      * @param d number to multiply
3063      * @param n power of 2
3064      * @return d &times; 2<sup>n</sup>
3065      */
3066     public static double scalb(final double d, final int n) {
3067 
3068         // first simple and fast handling when 2^n can be represented using normal numbers
3069         if ((n > -1023) && (n < 1024)) {
3070             return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3071         }
3072 
3073         // handle special cases
3074         if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3075             return d;
3076         }
3077         if (n < -2098) {
3078             return (d > 0) ? 0.0 : -0.0;
3079         }
3080         if (n > 2097) {
3081             return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3082         }
3083 
3084         // decompose d
3085         final long bits = Double.doubleToRawLongBits(d);
3086         final long sign = bits & 0x8000000000000000L;
3087         int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
3088         long mantissa   = bits & 0x000fffffffffffffL;
3089 
3090         // compute scaled exponent
3091         int scaledExponent = exponent + n;
3092 
3093         if (n < 0) {
3094             // we are really in the case n <= -1023
3095             if (scaledExponent > 0) {
3096                 // both the input and the result are normal numbers, we only adjust the exponent
3097                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3098             } else if (scaledExponent > -53) {
3099                 // the input is a normal number and the result is a subnormal number
3100 
3101                 // recover the hidden mantissa bit
3102                 mantissa |= 1L << 52;
3103 
3104                 // scales down complete mantissa, hence losing least significant bits
3105                 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3106                 mantissa >>>= 1 - scaledExponent;
3107                 if (mostSignificantLostBit != 0) {
3108                     // we need to add 1 bit to round up the result
3109                     mantissa++;
3110                 }
3111                 return Double.longBitsToDouble(sign | mantissa);
3112 
3113             } else {
3114                 // no need to compute the mantissa, the number scales down to 0
3115                 return (sign == 0L) ? 0.0 : -0.0;
3116             }
3117         } else {
3118             // we are really in the case n >= 1024
3119             if (exponent == 0) {
3120 
3121                 // the input number is subnormal, normalize it
3122                 while ((mantissa >>> 52) != 1) {
3123                     mantissa <<= 1;
3124                     --scaledExponent;
3125                 }
3126                 ++scaledExponent;
3127                 mantissa &= 0x000fffffffffffffL;
3128 
3129                 if (scaledExponent < 2047) {
3130                     return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3131                 } else {
3132                     return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3133                 }
3134 
3135             } else if (scaledExponent < 2047) {
3136                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3137             } else {
3138                 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3139             }
3140         }
3141 
3142     }
3143 
3144     /**
3145      * Multiply a float number by a power of 2.
3146      * @param f number to multiply
3147      * @param n power of 2
3148      * @return f &times; 2<sup>n</sup>
3149      */
3150     public static float scalb(final float f, final int n) {
3151 
3152         // first simple and fast handling when 2^n can be represented using normal numbers
3153         if ((n > -127) && (n < 128)) {
3154             return f * Float.intBitsToFloat((n + 127) << 23);
3155         }
3156 
3157         // handle special cases
3158         if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3159             return f;
3160         }
3161         if (n < -277) {
3162             return (f > 0) ? 0.0f : -0.0f;
3163         }
3164         if (n > 276) {
3165             return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3166         }
3167 
3168         // decompose f
3169         final int bits = Float.floatToIntBits(f);
3170         final int sign = bits & 0x80000000;
3171         int  exponent  = (bits >>> 23) & 0xff;
3172         int mantissa   = bits & 0x007fffff;
3173 
3174         // compute scaled exponent
3175         int scaledExponent = exponent + n;
3176 
3177         if (n < 0) {
3178             // we are really in the case n <= -127
3179             if (scaledExponent > 0) {
3180                 // both the input and the result are normal numbers, we only adjust the exponent
3181                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3182             } else if (scaledExponent > -24) {
3183                 // the input is a normal number and the result is a subnormal number
3184 
3185                 // recover the hidden mantissa bit
3186                 mantissa |= 1 << 23;
3187 
3188                 // scales down complete mantissa, hence losing least significant bits
3189                 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3190                 mantissa >>>= 1 - scaledExponent;
3191                 if (mostSignificantLostBit != 0) {
3192                     // we need to add 1 bit to round up the result
3193                     mantissa++;
3194                 }
3195                 return Float.intBitsToFloat(sign | mantissa);
3196 
3197             } else {
3198                 // no need to compute the mantissa, the number scales down to 0
3199                 return (sign == 0) ? 0.0f : -0.0f;
3200             }
3201         } else {
3202             // we are really in the case n >= 128
3203             if (exponent == 0) {
3204 
3205                 // the input number is subnormal, normalize it
3206                 while ((mantissa >>> 23) != 1) {
3207                     mantissa <<= 1;
3208                     --scaledExponent;
3209                 }
3210                 ++scaledExponent;
3211                 mantissa &= 0x007fffff;
3212 
3213                 if (scaledExponent < 255) {
3214                     return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3215                 } else {
3216                     return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3217                 }
3218 
3219             } else if (scaledExponent < 255) {
3220                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3221             } else {
3222                 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3223             }
3224         }
3225 
3226     }
3227 
3228     /**
3229      * Get the next machine representable number after a number, moving
3230      * in the direction of another number.
3231      * <p>
3232      * The ordering is as follows (increasing):
3233      * <ul>
3234      * <li>-INFINITY</li>
3235      * <li>-MAX_VALUE</li>
3236      * <li>-MIN_VALUE</li>
3237      * <li>-0.0</li>
3238      * <li>+0.0</li>
3239      * <li>+MIN_VALUE</li>
3240      * <li>+MAX_VALUE</li>
3241      * <li>+INFINITY</li>
3242      * <li></li>
3243      * <p>
3244      * If arguments compare equal, then the second argument is returned.
3245      * <p>
3246      * If {@code direction} is greater than {@code d},
3247      * the smallest machine representable number strictly greater than
3248      * {@code d} is returned; if less, then the largest representable number
3249      * strictly less than {@code d} is returned.</p>
3250      * <p>
3251      * If {@code d} is infinite and direction does not
3252      * bring it back to finite numbers, it is returned unchanged.</p>
3253      *
3254      * @param d base number
3255      * @param direction (the only important thing is whether
3256      * {@code direction} is greater or smaller than {@code d})
3257      * @return the next machine representable number in the specified direction
3258      */
3259     public static double nextAfter(double d, double direction) {
3260 
3261         // handling of some important special cases
3262         if (Double.isNaN(d) || Double.isNaN(direction)) {
3263             return Double.NaN;
3264         } else if (d == direction) {
3265             return direction;
3266         } else if (Double.isInfinite(d)) {
3267             return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3268         } else if (d == 0) {
3269             return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3270         }
3271         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3272         // are handled just as normal numbers
3273         // can use raw bits since already dealt with infinity and NaN
3274         final long bits = Double.doubleToRawLongBits(d);
3275         final long sign = bits & 0x8000000000000000L;
3276         if ((direction < d) ^ (sign == 0L)) {
3277             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3278         } else {
3279             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3280         }
3281 
3282     }
3283 
3284     /**
3285      * Get the next machine representable number after a number, moving
3286      * in the direction of another number.
3287      * <p>
3288      * The ordering is as follows (increasing):
3289      * <ul>
3290      * <li>-INFINITY</li>
3291      * <li>-MAX_VALUE</li>
3292      * <li>-MIN_VALUE</li>
3293      * <li>-0.0</li>
3294      * <li>+0.0</li>
3295      * <li>+MIN_VALUE</li>
3296      * <li>+MAX_VALUE</li>
3297      * <li>+INFINITY</li>
3298      * <li></li>
3299      * <p>
3300      * If arguments compare equal, then the second argument is returned.
3301      * <p>
3302      * If {@code direction} is greater than {@code f},
3303      * the smallest machine representable number strictly greater than
3304      * {@code f} is returned; if less, then the largest representable number
3305      * strictly less than {@code f} is returned.</p>
3306      * <p>
3307      * If {@code f} is infinite and direction does not
3308      * bring it back to finite numbers, it is returned unchanged.</p>
3309      *
3310      * @param f base number
3311      * @param direction (the only important thing is whether
3312      * {@code direction} is greater or smaller than {@code f})
3313      * @return the next machine representable number in the specified direction
3314      */
3315     public static float nextAfter(final float f, final double direction) {
3316 
3317         // handling of some important special cases
3318         if (Double.isNaN(f) || Double.isNaN(direction)) {
3319             return Float.NaN;
3320         } else if (f == direction) {
3321             return (float) direction;
3322         } else if (Float.isInfinite(f)) {
3323             return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3324         } else if (f == 0f) {
3325             return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3326         }
3327         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3328         // are handled just as normal numbers
3329 
3330         final int bits = Float.floatToIntBits(f);
3331         final int sign = bits & 0x80000000;
3332         if ((direction < f) ^ (sign == 0)) {
3333             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3334         } else {
3335             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3336         }
3337 
3338     }
3339 
3340     /** Get the largest whole number smaller than x.
3341      * @param x number from which floor is requested
3342      * @return a double number f such that f is an integer f <= x < f + 1.0
3343      */
3344     public static double floor(double x) {
3345         long y;
3346 
3347         if (x != x) { // NaN
3348             return x;
3349         }
3350 
3351         if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3352             return x;
3353         }
3354 
3355         y = (long) x;
3356         if (x < 0 && y != x) {
3357             y--;
3358         }
3359 
3360         if (y == 0) {
3361             return x*y;
3362         }
3363 
3364         return y;
3365     }
3366 
3367     /** Get the smallest whole number larger than x.
3368      * @param x number from which ceil is requested
3369      * @return a double number c such that c is an integer c - 1.0 < x <= c
3370      */
3371     public static double ceil(double x) {
3372         double y;
3373 
3374         if (x != x) { // NaN
3375             return x;
3376         }
3377 
3378         y = floor(x);
3379         if (y == x) {
3380             return y;
3381         }
3382 
3383         y += 1.0;
3384 
3385         if (y == 0) {
3386             return x*y;
3387         }
3388 
3389         return y;
3390     }
3391 
3392     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3393      * @param x number from which nearest whole number is requested
3394      * @return a double number r such that r is an integer r - 0.5 <= x <= r + 0.5
3395      */
3396     public static double rint(double x) {
3397         double y = floor(x);
3398         double d = x - y;
3399 
3400         if (d > 0.5) {
3401             if (y == -1.0) {
3402                 return -0.0; // Preserve sign of operand
3403             }
3404             return y+1.0;
3405         }
3406         if (d < 0.5) {
3407             return y;
3408         }
3409 
3410         /* half way, round to even */
3411         long z = (long) y;
3412         return (z & 1) == 0 ? y : y + 1.0;
3413     }
3414 
3415     /** Get the closest long to x.
3416      * @param x number from which closest long is requested
3417      * @return closest long to x
3418      */
3419     public static long round(double x) {
3420         return (long) floor(x + 0.5);
3421     }
3422 
3423     /** Get the closest int to x.
3424      * @param x number from which closest int is requested
3425      * @return closest int to x
3426      */
3427     public static int round(final float x) {
3428         return (int) floor(x + 0.5f);
3429     }
3430 
3431     /** Compute the minimum of two values
3432      * @param a first value
3433      * @param b second value
3434      * @return a if a is lesser or equal to b, b otherwise
3435      */
3436     public static int min(final int a, final int b) {
3437         return (a <= b) ? a : b;
3438     }
3439 
3440     /** Compute the minimum of two values
3441      * @param a first value
3442      * @param b second value
3443      * @return a if a is lesser or equal to b, b otherwise
3444      */
3445     public static long min(final long a, final long b) {
3446         return (a <= b) ? a : b;
3447     }
3448 
3449     /** Compute the minimum of two values
3450      * @param a first value
3451      * @param b second value
3452      * @return a if a is lesser or equal to b, b otherwise
3453      */
3454     public static float min(final float a, final float b) {
3455         if (a > b) {
3456             return b;
3457         }
3458         if (a < b) {
3459             return a;
3460         }
3461         /* if either arg is NaN, return NaN */
3462         if (a != b) {
3463             return Float.NaN;
3464         }
3465         /* min(+0.0,-0.0) == -0.0 */
3466         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3467         int bits = Float.floatToRawIntBits(a);
3468         if (bits == 0x80000000) {
3469             return a;
3470         }
3471         return b;
3472     }
3473 
3474     /** Compute the minimum of two values
3475      * @param a first value
3476      * @param b second value
3477      * @return a if a is lesser or equal to b, b otherwise
3478      */
3479     public static double min(final double a, final double b) {
3480         if (a > b) {
3481             return b;
3482         }
3483         if (a < b) {
3484             return a;
3485         }
3486         /* if either arg is NaN, return NaN */
3487         if (a != b) {
3488             return Double.NaN;
3489         }
3490         /* min(+0.0,-0.0) == -0.0 */
3491         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3492         long bits = Double.doubleToRawLongBits(a);
3493         if (bits == 0x8000000000000000L) {
3494             return a;
3495         }
3496         return b;
3497     }
3498 
3499     /** Compute the maximum of two values
3500      * @param a first value
3501      * @param b second value
3502      * @return b if a is lesser or equal to b, a otherwise
3503      */
3504     public static int max(final int a, final int b) {
3505         return (a <= b) ? b : a;
3506     }
3507 
3508     /** Compute the maximum of two values
3509      * @param a first value
3510      * @param b second value
3511      * @return b if a is lesser or equal to b, a otherwise
3512      */
3513     public static long max(final long a, final long b) {
3514         return (a <= b) ? b : a;
3515     }
3516 
3517     /** Compute the maximum of two values
3518      * @param a first value
3519      * @param b second value
3520      * @return b if a is lesser or equal to b, a otherwise
3521      */
3522     public static float max(final float a, final float b) {
3523         if (a > b) {
3524             return a;
3525         }
3526         if (a < b) {
3527             return b;
3528         }
3529         /* if either arg is NaN, return NaN */
3530         if (a != b) {
3531             return Float.NaN;
3532         }
3533         /* min(+0.0,-0.0) == -0.0 */
3534         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3535         int bits = Float.floatToRawIntBits(a);
3536         if (bits == 0x80000000) {
3537             return b;
3538         }
3539         return a;
3540     }
3541 
3542     /** Compute the maximum of two values
3543      * @param a first value
3544      * @param b second value
3545      * @return b if a is lesser or equal to b, a otherwise
3546      */
3547     public static double max(final double a, final double b) {
3548         if (a > b) {
3549             return a;
3550         }
3551         if (a < b) {
3552             return b;
3553         }
3554         /* if either arg is NaN, return NaN */
3555         if (a != b) {
3556             return Double.NaN;
3557         }
3558         /* min(+0.0,-0.0) == -0.0 */
3559         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3560         long bits = Double.doubleToRawLongBits(a);
3561         if (bits == 0x8000000000000000L) {
3562             return b;
3563         }
3564         return a;
3565     }
3566 
3567     /**
3568      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
3569      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br/>
3570      * avoiding intermediate overflow or underflow.
3571      *
3572      * <ul>
3573      * <li> If either argument is infinite, then the result is positive infinity.</li>
3574      * <li> else, if either argument is NaN then the result is NaN.</li>
3575      * </ul>
3576      *
3577      * @param x a value
3578      * @param y a value
3579      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
3580      */
3581     public static double hypot(final double x, final double y) {
3582         if (Double.isInfinite(x) || Double.isInfinite(y)) {
3583             return Double.POSITIVE_INFINITY;
3584         } else if (Double.isNaN(x) || Double.isNaN(y)) {
3585             return Double.NaN;
3586         } else {
3587 
3588             final int expX = getExponent(x);
3589             final int expY = getExponent(y);
3590             if (expX > expY + 27) {
3591                 // y is neglectible with respect to x
3592                 return abs(x);
3593             } else if (expY > expX + 27) {
3594                 // x is neglectible with respect to y
3595                 return abs(y);
3596             } else {
3597 
3598                 // find an intermediate scale to avoid both overflow and underflow
3599                 final int middleExp = (expX + expY) / 2;
3600 
3601                 // scale parameters without losing precision
3602                 final double scaledX = scalb(x, -middleExp);
3603                 final double scaledY = scalb(y, -middleExp);
3604 
3605                 // compute scaled hypotenuse
3606                 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
3607 
3608                 // remove scaling
3609                 return scalb(scaledH, middleExp);
3610 
3611             }
3612 
3613         }
3614     }
3615 
3616     /**
3617      * Computes the remainder as prescribed by the IEEE 754 standard.
3618      * The remainder value is mathematically equal to {@code x - y*n}
3619      * where {@code n} is the mathematical integer closest to the exact mathematical value
3620      * of the quotient {@code x/y}.
3621      * If two mathematical integers are equally close to {@code x/y} then
3622      * {@code n} is the integer that is even.
3623      * <p>
3624      * <ul>
3625      * <li>If either operand is NaN, the result is NaN.</li>
3626      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
3627      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
3628      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
3629      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
3630      * </ul>
3631      * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder}
3632      * @param dividend the number to be divided
3633      * @param divisor the number by which to divide
3634      * @return the remainder, rounded
3635      */
3636     public static double IEEEremainder(double dividend, double divisor) {
3637         return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
3638     }
3639 
3640     /** Convert a long to interger, detecting overflows
3641      * @param n number to convert to int
3642      * @return integer with same valie as n if no overflows occur
3643      * @exception MathArithmeticException if n cannot fit into an int
3644      * @since 3.4
3645      */
3646     public static int toIntExact(final long n) throws MathArithmeticException {
3647         if (n < Integer.MIN_VALUE || n > Integer.MAX_VALUE) {
3648             throw new MathArithmeticException(LocalizedFormats.OVERFLOW);
3649         }
3650         return (int) n;
3651     }
3652 
3653     /** Increment a number, detecting overflows.
3654      * @param n number to increment
3655      * @return n+1 if no overflows occur
3656      * @exception MathArithmeticException if an overflow occurs
3657      * @since 3.4
3658      */
3659     public static int incrementExact(final int n) throws MathArithmeticException {
3660 
3661         if (n == Integer.MAX_VALUE) {
3662             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, n, 1);
3663         }
3664 
3665         return n + 1;
3666 
3667     }
3668 
3669     /** Increment a number, detecting overflows.
3670      * @param n number to increment
3671      * @return n+1 if no overflows occur
3672      * @exception MathArithmeticException if an overflow occurs
3673      * @since 3.4
3674      */
3675     public static long incrementExact(final long n) throws MathArithmeticException {
3676 
3677         if (n == Long.MAX_VALUE) {
3678             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, n, 1);
3679         }
3680 
3681         return n + 1;
3682 
3683     }
3684 
3685     /** Decrement a number, detecting overflows.
3686      * @param n number to decrement
3687      * @return n-1 if no overflows occur
3688      * @exception MathArithmeticException if an overflow occurs
3689      * @since 3.4
3690      */
3691     public static int decrementExact(final int n) throws MathArithmeticException {
3692 
3693         if (n == Integer.MIN_VALUE) {
3694             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
3695         }
3696 
3697         return n - 1;
3698 
3699     }
3700 
3701     /** Decrement a number, detecting overflows.
3702      * @param n number to decrement
3703      * @return n-1 if no overflows occur
3704      * @exception MathArithmeticException if an overflow occurs
3705      * @since 3.4
3706      */
3707     public static long decrementExact(final long n) throws MathArithmeticException {
3708 
3709         if (n == Long.MIN_VALUE) {
3710             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
3711         }
3712 
3713         return n - 1;
3714 
3715     }
3716 
3717     /** Add two numbers, detecting overflows.
3718      * @param a first number to add
3719      * @param b second number to add
3720      * @return a+b if no overflows occur
3721      * @exception MathArithmeticException if an overflow occurs
3722      * @since 3.4
3723      */
3724     public static int addExact(final int a, final int b) throws MathArithmeticException {
3725 
3726         // compute sum
3727         final int sum = a + b;
3728 
3729         // check for overflow
3730         if ((a ^ b) >= 0 && (sum ^ b) < 0) {
3731             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, b);
3732         }
3733 
3734         return sum;
3735 
3736     }
3737 
3738     /** Add two numbers, detecting overflows.
3739      * @param a first number to add
3740      * @param b second number to add
3741      * @return a+b if no overflows occur
3742      * @exception MathArithmeticException if an overflow occurs
3743      * @since 3.4
3744      */
3745     public static long addExact(final long a, final long b) throws MathArithmeticException {
3746 
3747         // compute sum
3748         final long sum = a + b;
3749 
3750         // check for overflow
3751         if ((a ^ b) >= 0 && (sum ^ b) < 0) {
3752             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_ADDITION, a, b);
3753         }
3754 
3755         return sum;
3756 
3757     }
3758 
3759     /** Subtract two numbers, detecting overflows.
3760      * @param a first number
3761      * @param b second number to subtract from a
3762      * @return a-b if no overflows occur
3763      * @exception MathArithmeticException if an overflow occurs
3764      * @since 3.4
3765      */
3766     public static int subtractExact(final int a, final int b) {
3767 
3768         // compute subtraction
3769         final int sub = a - b;
3770 
3771         // check for overflow
3772         if ((a ^ b) < 0 && (sub ^ b) >= 0) {
3773             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, a, b);
3774         }
3775 
3776         return sub;
3777 
3778     }
3779 
3780     /** Subtract two numbers, detecting overflows.
3781      * @param a first number
3782      * @param b second number to subtract from a
3783      * @return a-b if no overflows occur
3784      * @exception MathArithmeticException if an overflow occurs
3785      * @since 3.4
3786      */
3787     public static long subtractExact(final long a, final long b) {
3788 
3789         // compute subtraction
3790         final long sub = a - b;
3791 
3792         // check for overflow
3793         if ((a ^ b) < 0 && (sub ^ b) >= 0) {
3794             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_SUBTRACTION, a, b);
3795         }
3796 
3797         return sub;
3798 
3799     }
3800 
3801     /** Multiply two numbers, detecting overflows.
3802      * @param a first number to multiply
3803      * @param b second number to multiply
3804      * @return a*b if no overflows occur
3805      * @exception MathArithmeticException if an overflow occurs
3806      * @since 3.4
3807      */
3808     public static int multiplyExact(final int a, final int b) {
3809         if (((b  >  0)  && (a > Integer.MAX_VALUE / b || a < Integer.MIN_VALUE / b)) ||
3810             ((b  < -1)  && (a > Integer.MIN_VALUE / b || a < Integer.MAX_VALUE / b)) ||
3811             ((b == -1)  && (a == Integer.MIN_VALUE))) {
3812             throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
3813         }
3814         return a * b;
3815     }
3816 
3817     /** Multiply two numbers, detecting overflows.
3818      * @param a first number to multiply
3819      * @param b second number to multiply
3820      * @return a*b if no overflows occur
3821      * @exception MathArithmeticException if an overflow occurs
3822      * @since 3.4
3823      */
3824     public static long multiplyExact(final long a, final long b) {
3825         if (((b  >  0l)  && (a > Long.MAX_VALUE / b || a < Long.MIN_VALUE / b)) ||
3826             ((b  < -1l)  && (a > Long.MIN_VALUE / b || a < Long.MAX_VALUE / b)) ||
3827             ((b == -1l)  && (a == Long.MIN_VALUE))) {
3828                 throw new MathArithmeticException(LocalizedFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
3829             }
3830             return a * b;
3831     }
3832 
3833     /** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
3834      * <p>
3835      * This methods returns the same value as integer division when
3836      * a and b are same signs, but returns a different value when
3837      * they are opposite (i.e. q is negative).
3838      * </p>
3839      * @param a dividend
3840      * @param b divisor
3841      * @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
3842      * @exception MathArithmeticException if b == 0
3843      * @see #floorMod(int, int)
3844      * @since 3.4
3845      */
3846     public static int floorDiv(final int a, final int b) throws MathArithmeticException {
3847 
3848         if (b == 0) {
3849             throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
3850         }
3851 
3852         final int m = a % b;
3853         if ((a ^ b) >= 0 || m == 0) {
3854             // a an b have same sign, or division is exact
3855             return a / b;
3856         } else {
3857             // a and b have opposite signs and division is not exact
3858             return (a / b) - 1;
3859         }
3860 
3861     }
3862 
3863     /** Finds q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
3864      * <p>
3865      * This methods returns the same value as integer division when
3866      * a and b are same signs, but returns a different value when
3867      * they are opposite (i.e. q is negative).
3868      * </p>
3869      * @param a dividend
3870      * @param b divisor
3871      * @return q such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
3872      * @exception MathArithmeticException if b == 0
3873      * @see #floorMod(long, long)
3874      * @since 3.4
3875      */
3876     public static long floorDiv(final long a, final long b) throws MathArithmeticException {
3877 
3878         if (b == 0l) {
3879             throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
3880         }
3881 
3882         final long m = a % b;
3883         if ((a ^ b) >= 0l || m == 0l) {
3884             // a an b have same sign, or division is exact
3885             return a / b;
3886         } else {
3887             // a and b have opposite signs and division is not exact
3888             return (a / b) - 1l;
3889         }
3890 
3891     }
3892 
3893     /** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
3894      * <p>
3895      * This methods returns the same value as integer modulo when
3896      * a and b are same signs, but returns a different value when
3897      * they are opposite (i.e. q is negative).
3898      * </p>
3899      * @param a dividend
3900      * @param b divisor
3901      * @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
3902      * @exception MathArithmeticException if b == 0
3903      * @see #floorDiv(int, int)
3904      * @since 3.4
3905      */
3906     public static int floorMod(final int a, final int b) throws MathArithmeticException {
3907 
3908         if (b == 0) {
3909             throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
3910         }
3911 
3912         final int m = a % b;
3913         if ((a ^ b) >= 0 || m == 0) {
3914             // a an b have same sign, or division is exact
3915             return m;
3916         } else {
3917             // a and b have opposite signs and division is not exact
3918             return b + m;
3919         }
3920 
3921     }
3922 
3923     /** Finds r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0.
3924      * <p>
3925      * This methods returns the same value as integer modulo when
3926      * a and b are same signs, but returns a different value when
3927      * they are opposite (i.e. q is negative).
3928      * </p>
3929      * @param a dividend
3930      * @param b divisor
3931      * @return r such that a = q b + r with 0 <= r < b if b > 0 and b < r <= 0 if b > 0
3932      * @exception MathArithmeticException if b == 0
3933      * @see #floorDiv(long, long)
3934      * @since 3.4
3935      */
3936     public static long floorMod(final long a, final long b) {
3937 
3938         if (b == 0l) {
3939             throw new MathArithmeticException(LocalizedFormats.ZERO_DENOMINATOR);
3940         }
3941 
3942         final long m = a % b;
3943         if ((a ^ b) >= 0l || m == 0l) {
3944             // a an b have same sign, or division is exact
3945             return m;
3946         } else {
3947             // a and b have opposite signs and division is not exact
3948             return b + m;
3949         }
3950 
3951     }
3952 
3953     /**
3954      * Returns the first argument with the sign of the second argument.
3955      * A NaN {@code sign} argument is treated as positive.
3956      *
3957      * @param magnitude the value to return
3958      * @param sign the sign for the returned value
3959      * @return the magnitude with the same sign as the {@code sign} argument
3960      */
3961     public static double copySign(double magnitude, double sign){
3962         // The highest order bit is going to be zero if the
3963         // highest order bit of m and s is the same and one otherwise.
3964         // So (m^s) will be positive if both m and s have the same sign
3965         // and negative otherwise.
3966         final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
3967         final long s = Double.doubleToRawLongBits(sign);
3968         if ((m^s) >= 0) {
3969             return magnitude;
3970         }
3971         return -magnitude; // flip sign
3972     }
3973 
3974     /**
3975      * Returns the first argument with the sign of the second argument.
3976      * A NaN {@code sign} argument is treated as positive.
3977      *
3978      * @param magnitude the value to return
3979      * @param sign the sign for the returned value
3980      * @return the magnitude with the same sign as the {@code sign} argument
3981      */
3982     public static float copySign(float magnitude, float sign){
3983         // The highest order bit is going to be zero if the
3984         // highest order bit of m and s is the same and one otherwise.
3985         // So (m^s) will be positive if both m and s have the same sign
3986         // and negative otherwise.
3987         final int m = Float.floatToRawIntBits(magnitude);
3988         final int s = Float.floatToRawIntBits(sign);
3989         if ((m^s) >= 0) {
3990             return magnitude;
3991         }
3992         return -magnitude; // flip sign
3993     }
3994 
3995     /**
3996      * Return the exponent of a double number, removing the bias.
3997      * <p>
3998      * For double numbers of the form 2<sup>x</sup>, the unbiased
3999      * exponent is exactly x.
4000      * </p>
4001      * @param d number from which exponent is requested
4002      * @return exponent for d in IEEE754 representation, without bias
4003      */
4004     public static int getExponent(final double d) {
4005         // NaN and Infinite will return 1024 anywho so can use raw bits
4006         return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
4007     }
4008 
4009     /**
4010      * Return the exponent of a float number, removing the bias.
4011      * <p>
4012      * For float numbers of the form 2<sup>x</sup>, the unbiased
4013      * exponent is exactly x.
4014      * </p>
4015      * @param f number from which exponent is requested
4016      * @return exponent for d in IEEE754 representation, without bias
4017      */
4018     public static int getExponent(final float f) {
4019         // NaN and Infinite will return the same exponent anywho so can use raw bits
4020         return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
4021     }
4022 
4023     /**
4024      * Print out contents of arrays, and check the length.
4025      * <p>used to generate the preset arrays originally.</p>
4026      * @param a unused
4027      */
4028     public static void main(String[] a) {
4029         PrintStream out = System.out;
4030         FastMathCalc.printarray(out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
4031         FastMathCalc.printarray(out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
4032         FastMathCalc.printarray(out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
4033         FastMathCalc.printarray(out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
4034         FastMathCalc.printarray(out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
4035         FastMathCalc.printarray(out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
4036         FastMathCalc.printarray(out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
4037         FastMathCalc.printarray(out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
4038         FastMathCalc.printarray(out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
4039         FastMathCalc.printarray(out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
4040         FastMathCalc.printarray(out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
4041     }
4042 
4043     /** Enclose large data table in nested static class so it's only loaded on first access. */
4044     private static class ExpIntTable {
4045         /** Exponential evaluated at integer values,
4046          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
4047          */
4048         private static final double[] EXP_INT_TABLE_A;
4049         /** Exponential evaluated at integer values,
4050          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
4051          */
4052         private static final double[] EXP_INT_TABLE_B;
4053 
4054         static {
4055             if (RECOMPUTE_TABLES_AT_RUNTIME) {
4056                 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
4057                 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
4058 
4059                 final double tmp[] = new double[2];
4060                 final double recip[] = new double[2];
4061 
4062                 // Populate expIntTable
4063                 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
4064                     FastMathCalc.expint(i, tmp);
4065                     EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
4066                     EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
4067 
4068                     if (i != 0) {
4069                         // Negative integer powers
4070                         FastMathCalc.splitReciprocal(tmp, recip);
4071                         EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
4072                         EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
4073                     }
4074                 }
4075             } else {
4076                 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
4077                 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
4078             }
4079         }
4080     }
4081 
4082     /** Enclose large data table in nested static class so it's only loaded on first access. */
4083     private static class ExpFracTable {
4084         /** Exponential over the range of 0 - 1 in increments of 2^-10
4085          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
4086          * 1024 = 2^10
4087          */
4088         private static final double[] EXP_FRAC_TABLE_A;
4089         /** Exponential over the range of 0 - 1 in increments of 2^-10
4090          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
4091          */
4092         private static final double[] EXP_FRAC_TABLE_B;
4093 
4094         static {
4095             if (RECOMPUTE_TABLES_AT_RUNTIME) {
4096                 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
4097                 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
4098 
4099                 final double tmp[] = new double[2];
4100 
4101                 // Populate expFracTable
4102                 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
4103                 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
4104                     FastMathCalc.slowexp(i * factor, tmp);
4105                     EXP_FRAC_TABLE_A[i] = tmp[0];
4106                     EXP_FRAC_TABLE_B[i] = tmp[1];
4107                 }
4108             } else {
4109                 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
4110                 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
4111             }
4112         }
4113     }
4114 
4115     /** Enclose large data table in nested static class so it's only loaded on first access. */
4116     private static class lnMant {
4117         /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
4118         private static final double[][] LN_MANT;
4119 
4120         static {
4121             if (RECOMPUTE_TABLES_AT_RUNTIME) {
4122                 LN_MANT = new double[FastMath.LN_MANT_LEN][];
4123 
4124                 // Populate lnMant table
4125                 for (int i = 0; i < LN_MANT.length; i++) {
4126                     final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
4127                     LN_MANT[i] = FastMathCalc.slowLog(d);
4128                 }
4129             } else {
4130                 LN_MANT = FastMathLiteralArrays.loadLnMant();
4131             }
4132         }
4133     }
4134 
4135     /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
4136     private static class CodyWaite {
4137         /** k */
4138         private final int finalK;
4139         /** remA */
4140         private final double finalRemA;
4141         /** remB */
4142         private final double finalRemB;
4143 
4144         /**
4145          * @param xa Argument.
4146          */
4147         CodyWaite(double xa) {
4148             // Estimate k.
4149             //k = (int)(xa / 1.5707963267948966);
4150             int k = (int)(xa * 0.6366197723675814);
4151 
4152             // Compute remainder.
4153             double remA;
4154             double remB;
4155             while (true) {
4156                 double a = -k * 1.570796251296997;
4157                 remA = xa + a;
4158                 remB = -(remA - xa - a);
4159 
4160                 a = -k * 7.549789948768648E-8;
4161                 double b = remA;
4162                 remA = a + b;
4163                 remB += -(remA - b - a);
4164 
4165                 a = -k * 6.123233995736766E-17;
4166                 b = remA;
4167                 remA = a + b;
4168                 remB += -(remA - b - a);
4169 
4170                 if (remA > 0) {
4171                     break;
4172                 }
4173 
4174                 // Remainder is negative, so decrement k and try again.
4175                 // This should only happen if the input is very close
4176                 // to an even multiple of pi/2.
4177                 --k;
4178             }
4179 
4180             this.finalK = k;
4181             this.finalRemA = remA;
4182             this.finalRemB = remB;
4183         }
4184 
4185         /**
4186          * @return k
4187          */
4188         int getK() {
4189             return finalK;
4190         }
4191         /**
4192          * @return remA
4193          */
4194         double getRemA() {
4195             return finalRemA;
4196         }
4197         /**
4198          * @return remB
4199          */
4200         double getRemB() {
4201             return finalRemB;
4202         }
4203     }
4204 }