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