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