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