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