View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.util;
23  
24  import org.hipparchus.CalculusFieldElement;
25  import org.hipparchus.exception.LocalizedCoreFormats;
26  import org.hipparchus.exception.MathRuntimeException;
27  
28  /**
29   * Faster, more accurate, portable alternative to {@link Math} and
30   * {@link StrictMath} for large scale computation.
31   * <p>
32   * FastMath is a drop-in replacement for both Math and StrictMath. This
33   * means that for any method in Math (say {@code Math.sin(x)} or
34   * {@code Math.cbrt(y)}), user can directly change the class and use the
35   * methods as is (using {@code FastMath.sin(x)} or {@code FastMath.cbrt(y)}
36   * in the previous example).
37   * <p>
38   * FastMath speed is achieved by relying heavily on optimizing compilers
39   * to native code present in many JVMs today and use of large tables.
40   * The larger tables are lazily initialized on first use, so that the setup
41   * time does not penalize methods that don't need them.
42   * <p>
43   * Note that FastMath is
44   * extensively used inside Hipparchus, so by calling some algorithms,
45   * the overhead when the the tables need to be initialized will occur
46   * regardless of the end-user calling FastMath methods directly or not.
47   * Performance figures for a specific JVM and hardware can be evaluated by
48   * running the FastMathTestPerformance tests in the test directory of the source
49   * distribution.
50   * <p>
51   * FastMath accuracy should be mostly independent of the JVM as it relies only
52   * on IEEE-754 basic operations and on embedded tables. Almost all operations
53   * are accurate to about 0.5 ulp throughout the domain range. This statement,
54   * of course is only a rough global observed behavior, it is <em>not</em> a
55   * guarantee for <em>every</em> double numbers input (see William Kahan's <a
56   * href="http://en.wikipedia.org/wiki/Rounding#The_table-maker.27s_dilemma">Table
57   * Maker's Dilemma</a>).
58   * <p>
59   * FastMath additionally implements the following methods not found in Math/StrictMath:
60   * <ul>
61   * <li>{@link #asinh(double)}</li>
62   * <li>{@link #acosh(double)}</li>
63   * <li>{@link #atanh(double)}</li>
64   * </ul>
65   * The following methods are found in Math/StrictMath since 1.6 only, they are provided
66   * by FastMath even in 1.5 Java virtual machines
67   * <ul>
68   * <li>{@link #copySign(double, double)}</li>
69   * <li>{@link #getExponent(double)}</li>
70   * <li>{@link #nextAfter(double,double)}</li>
71   * <li>{@link #nextUp(double)}</li>
72   * <li>{@link #scalb(double, int)}</li>
73   * <li>{@link #copySign(float, float)}</li>
74   * <li>{@link #getExponent(float)}</li>
75   * <li>{@link #nextAfter(float,double)}</li>
76   * <li>{@link #nextUp(float)}</li>
77   * <li>{@link #scalb(float, int)}</li>
78   * </ul>
79   */
80  public class FastMath {
81      /** Archimede's constant PI, ratio of circle circumference to diameter. */
82      public static final double PI = 105414357.0 / 33554432.0 + 1.984187159361080883e-9;
83  
84      /** Napier's constant e, base of the natural logarithm. */
85      public static final double E = 2850325.0 / 1048576.0 + 8.254840070411028747e-8;
86  
87      /** Index of exp(0) in the array of integer exponentials. */
88      static final int EXP_INT_TABLE_MAX_INDEX = 750;
89      /** Length of the array of integer exponentials. */
90      static final int EXP_INT_TABLE_LEN = EXP_INT_TABLE_MAX_INDEX * 2;
91      /** Logarithm table length. */
92      static final int LN_MANT_LEN = 1024;
93      /** Exponential fractions table length. */
94      static final int EXP_FRAC_TABLE_LEN = 1025; // 0, 1/1024, ... 1024/1024
95  
96      /** StrictMath.log(Double.MAX_VALUE): {@value} */
97      private static final double LOG_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);
98  
99      /** 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 
139     /** Sine table (high bits). */
140     private static final double SINE_TABLE_A[] =
141         {
142         +0.0d,
143         +0.1246747374534607d,
144         +0.24740394949913025d,
145         +0.366272509098053d,
146         +0.4794255495071411d,
147         +0.5850973129272461d,
148         +0.6816387176513672d,
149         +0.7675435543060303d,
150         +0.8414709568023682d,
151         +0.902267575263977d,
152         +0.9489846229553223d,
153         +0.9808930158615112d,
154         +0.9974949359893799d,
155         +0.9985313415527344d,
156     };
157 
158     /** Sine table (low bits). */
159     private static final double SINE_TABLE_B[] =
160         {
161         +0.0d,
162         -4.068233003401932E-9d,
163         +9.755392680573412E-9d,
164         +1.9987994582857286E-8d,
165         -1.0902938113007961E-8d,
166         -3.9986783938944604E-8d,
167         +4.23719669792332E-8d,
168         -5.207000323380292E-8d,
169         +2.800552834259E-8d,
170         +1.883511811213715E-8d,
171         -3.5997360512765566E-9d,
172         +4.116164446561962E-8d,
173         +5.0614674548127384E-8d,
174         -1.0129027912496858E-9d,
175     };
176 
177     /** Cosine table (high bits). */
178     private static final double COSINE_TABLE_A[] =
179         {
180         +1.0d,
181         +0.9921976327896118d,
182         +0.9689123630523682d,
183         +0.9305076599121094d,
184         +0.8775825500488281d,
185         +0.8109631538391113d,
186         +0.7316888570785522d,
187         +0.6409968137741089d,
188         +0.5403022766113281d,
189         +0.4311765432357788d,
190         +0.3153223395347595d,
191         +0.19454771280288696d,
192         +0.07073719799518585d,
193         -0.05417713522911072d,
194     };
195 
196     /** Cosine table (low bits). */
197     private static final double COSINE_TABLE_B[] =
198         {
199         +0.0d,
200         +3.4439717236742845E-8d,
201         +5.865827662008209E-8d,
202         -3.7999795083850525E-8d,
203         +1.184154459111628E-8d,
204         -3.43338934259355E-8d,
205         +1.1795268640216787E-8d,
206         +4.438921624363781E-8d,
207         +2.925681159240093E-8d,
208         -2.6437112632041807E-8d,
209         +2.2860509143963117E-8d,
210         -4.813899778443457E-9d,
211         +3.6725170580355583E-9d,
212         +2.0217439756338078E-10d,
213     };
214 
215 
216     /** Tangent table, used by atan() (high bits). */
217     private static final double TANGENT_TABLE_A[] =
218         {
219         +0.0d,
220         +0.1256551444530487d,
221         +0.25534194707870483d,
222         +0.3936265707015991d,
223         +0.5463024377822876d,
224         +0.7214844226837158d,
225         +0.9315965175628662d,
226         +1.1974215507507324d,
227         +1.5574076175689697d,
228         +2.092571258544922d,
229         +3.0095696449279785d,
230         +5.041914939880371d,
231         +14.101419448852539d,
232         -18.430862426757812d,
233     };
234 
235     /** Tangent table, used by atan() (low bits). */
236     private static final double TANGENT_TABLE_B[] =
237         {
238         +0.0d,
239         -7.877917738262007E-9d,
240         -2.5857668567479893E-8d,
241         +5.2240336371356666E-9d,
242         +5.206150291559893E-8d,
243         +1.8307188599677033E-8d,
244         -5.7618793749770706E-8d,
245         +7.848361555046424E-8d,
246         +1.0708593250394448E-7d,
247         +1.7827257129423813E-8d,
248         +2.893485277253286E-8d,
249         +3.1660099222737955E-7d,
250         +4.983191803254889E-7d,
251         -3.356118100840571E-7d,
252     };
253 
254     /** Bits of 1/(2*pi), need for reducePayneHanek(). */
255     private static final long RECIP_2PI[] = {
256         (0x28be60dbL << 32) | 0x9391054aL,
257         (0x7f09d5f4L << 32) | 0x7d4d3770L,
258         (0x36d8a566L << 32) | 0x4f10e410L,
259         (0x7f9458eaL << 32) | 0xf7aef158L,
260         (0x6dc91b8eL << 32) | 0x909374b8L,
261         (0x01924bbaL << 32) | 0x82746487L,
262         (0x3f877ac7L << 32) | 0x2c4a69cfL,
263         (0xba208d7dL << 32) | 0x4baed121L,
264         (0x3a671c09L << 32) | 0xad17df90L,
265         (0x4e64758eL << 32) | 0x60d4ce7dL,
266         (0x272117e2L << 32) | 0xef7e4a0eL,
267         (0xc7fe25ffL << 32) | 0xf7816603L,
268         (0xfbcbc462L << 32) | 0xd6829b47L,
269         (0xdb4d9fb3L << 32) | 0xc9f2c26dL,
270         (0xd3d18fd9L << 32) | 0xa797fa8bL,
271         (0x5d49eeb1L << 32) | 0xfaf97c5eL,
272         (0xcf41ce7dL << 32) | 0xe294a4baL,
273          0x9afed7ecL << 32  };
274 
275     /** Bits of pi/4, need for reducePayneHanek(). */
276     private static final long PI_O_4_BITS[] = {
277         (0xc90fdaa2L << 32) | 0x2168c234L,
278         (0xc4c6628bL << 32) | 0x80dc1cd1L };
279 
280     /** Eighths.
281      * This is used by sinQ, because its faster to do a table lookup than
282      * a multiply in this time-critical routine
283      */
284     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};
285 
286     /** Table of 2^((n+2)/3) */
287     private static final double CBRTTWO[] = { 0.6299605249474366,
288                                             0.7937005259840998,
289                                             1.0,
290                                             1.2599210498948732,
291                                             1.5874010519681994 };
292 
293     /*
294      *  There are 52 bits in the mantissa of a double.
295      *  For additional precision, the code splits double numbers into two parts,
296      *  by clearing the low order 30 bits if possible, and then performs the arithmetic
297      *  on each half separately.
298      */
299 
300     /**
301      * 0x40000000 - used to split a double into two parts, both with the low order bits cleared.
302      * Equivalent to 2^30.
303      */
304     private static final long HEX_40000000 = 0x40000000L; // 1073741824L
305 
306     /** Mask used to clear low order 30 bits */
307     private static final long MASK_30BITS = -1L - (HEX_40000000 -1); // 0xFFFFFFFFC0000000L;
308 
309     /** Mask used to clear the non-sign part of an int. */
310     private static final int MASK_NON_SIGN_INT = 0x7fffffff;
311 
312     /** Mask used to clear the non-sign part of a long. */
313     private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
314 
315     /** Mask used to extract exponent from double bits. */
316     private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
317 
318     /** Mask used to extract mantissa from double bits. */
319     private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
320 
321     /** Mask used to add implicit high order bit for normalized double. */
322     private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
323 
324     /** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
325     private static final double TWO_POWER_52 = 4503599627370496.0;
326 
327     /** Constant: {@value}. */
328     private static final double F_1_3 = 1d / 3d;
329     /** Constant: {@value}. */
330     private static final double F_1_5 = 1d / 5d;
331     /** Constant: {@value}. */
332     private static final double F_1_7 = 1d / 7d;
333     /** Constant: {@value}. */
334     private static final double F_1_9 = 1d / 9d;
335     /** Constant: {@value}. */
336     private static final double F_1_11 = 1d / 11d;
337     /** Constant: {@value}. */
338     private static final double F_1_13 = 1d / 13d;
339     /** Constant: {@value}. */
340     private static final double F_1_15 = 1d / 15d;
341     /** Constant: {@value}. */
342     private static final double F_1_17 = 1d / 17d;
343     /** Constant: {@value}. */
344     private static final double F_3_4 = 3d / 4d;
345     /** Constant: {@value}. */
346     private static final double F_15_16 = 15d / 16d;
347     /** Constant: {@value}. */
348     private static final double F_13_14 = 13d / 14d;
349     /** Constant: {@value}. */
350     private static final double F_11_12 = 11d / 12d;
351     /** Constant: {@value}. */
352     private static final double F_9_10 = 9d / 10d;
353     /** Constant: {@value}. */
354     private static final double F_7_8 = 7d / 8d;
355     /** Constant: {@value}. */
356     private static final double F_5_6 = 5d / 6d;
357     /** Constant: {@value}. */
358     private static final double F_1_2 = 1d / 2d;
359     /** Constant: {@value}. */
360     private static final double F_1_4 = 1d / 4d;
361 
362     /**
363      * Private Constructor
364      */
365     private FastMath() {}
366 
367     // Generic helper methods
368 
369     /**
370      * Get the high order bits from the mantissa.
371      * Equivalent to adding and subtracting HEX_40000 but also works for very large numbers
372      *
373      * @param d the value to split
374      * @return the high order part of the mantissa
375      */
376     private static double doubleHighPart(double d) {
377         if (d > -Precision.SAFE_MIN && d < Precision.SAFE_MIN){
378             return d; // These are un-normalised - don't try to convert
379         }
380         long xl = Double.doubleToRawLongBits(d); // can take raw bits because just gonna convert it back
381         xl &= MASK_30BITS; // Drop low order bits
382         return Double.longBitsToDouble(xl);
383     }
384 
385     /** Compute the square root of a number.
386      * <p><b>Note:</b> this implementation currently delegates to {@link Math#sqrt}
387      * @param a number on which evaluation is done
388      * @return square root of a
389      */
390     public static double sqrt(final double a) {
391         return Math.sqrt(a);
392     }
393 
394     /** Compute the hyperbolic cosine of a number.
395      * @param x number on which evaluation is done
396      * @return hyperbolic cosine of x
397      */
398     public static double cosh(double x) {
399       if (Double.isNaN(x)) {
400           return x;
401       }
402 
403       // cosh[z] = (exp(z) + exp(-z))/2
404 
405       // for numbers with magnitude 20 or so,
406       // exp(-z) can be ignored in comparison with exp(z)
407 
408       if (x > 20) {
409           if (x >= LOG_MAX_VALUE) {
410               // Avoid overflow (MATH-905).
411               final double t = exp(0.5 * x);
412               return (0.5 * t) * t;
413           } else {
414               return 0.5 * exp(x);
415           }
416       } else if (x < -20) {
417           if (x <= -LOG_MAX_VALUE) {
418               // Avoid overflow (MATH-905).
419               final double t = exp(-0.5 * x);
420               return (0.5 * t) * t;
421           } else {
422               return 0.5 * exp(-x);
423           }
424       }
425 
426       final double hiPrec[] = new double[2];
427       if (x < 0.0) {
428           x = -x;
429       }
430       exp(x, 0.0, hiPrec);
431 
432       double ya = hiPrec[0] + hiPrec[1];
433       double yb = -(ya - hiPrec[0] - hiPrec[1]);
434 
435       double temp = ya * HEX_40000000;
436       double yaa = ya + temp - temp;
437       double yab = ya - yaa;
438 
439       // recip = 1/y
440       double recip = 1.0/ya;
441       temp = recip * HEX_40000000;
442       double recipa = recip + temp - temp;
443       double recipb = recip - recipa;
444 
445       // Correct for rounding in division
446       recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
447       // Account for yb
448       recipb += -yb * recip * recip;
449 
450       // y = y + 1/y
451       temp = ya + recipa;
452       yb += -(temp - ya - recipa);
453       ya = temp;
454       temp = ya + recipb;
455       yb += -(temp - ya - recipb);
456       ya = temp;
457 
458       double result = ya + yb;
459       result *= 0.5;
460       return result;
461     }
462 
463     /** Compute the hyperbolic sine of a number.
464      * @param x number on which evaluation is done
465      * @return hyperbolic sine of x
466      */
467     public static double sinh(double x) {
468       boolean negate = false;
469       if (Double.isNaN(x)) {
470           return x;
471       }
472 
473       // sinh[z] = (exp(z) - exp(-z) / 2
474 
475       // for values of z larger than about 20,
476       // exp(-z) can be ignored in comparison with exp(z)
477 
478       if (x > 20) {
479           if (x >= LOG_MAX_VALUE) {
480               // Avoid overflow (MATH-905).
481               final double t = exp(0.5 * x);
482               return (0.5 * t) * t;
483           } else {
484               return 0.5 * exp(x);
485           }
486       } else if (x < -20) {
487           if (x <= -LOG_MAX_VALUE) {
488               // Avoid overflow (MATH-905).
489               final double t = exp(-0.5 * x);
490               return (-0.5 * t) * t;
491           } else {
492               return -0.5 * exp(-x);
493           }
494       }
495 
496       if (x == 0) {
497           return x;
498       }
499 
500       if (x < 0.0) {
501           x = -x;
502           negate = true;
503       }
504 
505       double hiPrec[] = new double[2];
506       double result;
507 
508       if (x > 0.25) {
509           exp(x, 0.0, hiPrec);
510 
511           double ya = hiPrec[0] + hiPrec[1];
512           double yb = -(ya - hiPrec[0] - hiPrec[1]);
513 
514           double temp = ya * HEX_40000000;
515           double yaa = ya + temp - temp;
516           double yab = ya - yaa;
517 
518           // recip = 1/y
519           double recip = 1.0/ya;
520           temp = recip * HEX_40000000;
521           double recipa = recip + temp - temp;
522           double recipb = recip - recipa;
523 
524           // Correct for rounding in division
525           recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
526           // Account for yb
527           recipb += -yb * recip * recip;
528 
529           recipa = -recipa;
530           recipb = -recipb;
531 
532           // y = y - 1/y
533           temp = ya + recipa;
534           yb += -(temp - ya - recipa);
535           ya = temp;
536           temp = ya + recipb;
537           yb += -(temp - ya - recipb);
538           ya = temp;
539 
540           result = ya + yb;
541           result *= 0.5;
542       } else {
543           expm1(x, hiPrec);
544 
545           double ya = hiPrec[0] + hiPrec[1];
546           double yb = -(ya - hiPrec[0] - hiPrec[1]);
547 
548           /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
549           double denom = 1.0 + ya;
550           double denomr = 1.0 / denom;
551           double denomb = -(denom - 1.0 - ya) + yb;
552           double ratio = ya * denomr;
553           double temp = ratio * HEX_40000000;
554           double ra = ratio + temp - temp;
555           double rb = ratio - ra;
556 
557           temp = denom * HEX_40000000;
558           double za = denom + temp - temp;
559           double zb = denom - za;
560 
561           rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
562 
563           // Adjust for yb
564           rb += yb*denomr;                        // numerator
565           rb += -ya * denomb * denomr * denomr;   // denominator
566 
567           // y = y - 1/y
568           temp = ya + ra;
569           yb += -(temp - ya - ra);
570           ya = temp;
571           temp = ya + rb;
572           yb += -(temp - ya - rb);
573           ya = temp;
574 
575           result = ya + yb;
576           result *= 0.5;
577       }
578 
579       if (negate) {
580           result = -result;
581       }
582 
583       return result;
584     }
585 
586     /**
587      * Combined hyperbolic sine and hyperbolic cosine function.
588      *
589      * @param x Argument.
590      * @return [sinh(x), cosh(x)]
591      */
592     public static SinhCosh sinhCosh(double x) {
593       boolean negate = false;
594       if (Double.isNaN(x)) {
595           return new SinhCosh(x, x);
596       }
597 
598       // sinh[z] = (exp(z) - exp(-z) / 2
599       // cosh[z] = (exp(z) + exp(-z))/2
600 
601       // for values of z larger than about 20,
602       // exp(-z) can be ignored in comparison with exp(z)
603 
604       if (x > 20) {
605           final double e;
606           if (x >= LOG_MAX_VALUE) {
607               // Avoid overflow (MATH-905).
608               final double t = exp(0.5 * x);
609               e = (0.5 * t) * t;
610           } else {
611               e = 0.5 * exp(x);
612           }
613           return new SinhCosh(e, e);
614       } else if (x < -20) {
615           final double e;
616           if (x <= -LOG_MAX_VALUE) {
617               // Avoid overflow (MATH-905).
618               final double t = exp(-0.5 * x);
619               e = (-0.5 * t) * t;
620           } else {
621               e = -0.5 * exp(-x);
622           }
623           return new SinhCosh(e, -e);
624       }
625 
626       if (x == 0) {
627           return new SinhCosh(x, 1.0);
628       }
629 
630       if (x < 0.0) {
631           x = -x;
632           negate = true;
633       }
634 
635       double hiPrec[] = new double[2];
636       double resultM;
637       double resultP;
638 
639       if (x > 0.25) {
640           exp(x, 0.0, hiPrec);
641 
642           final double ya = hiPrec[0] + hiPrec[1];
643           final double yb = -(ya - hiPrec[0] - hiPrec[1]);
644 
645           double temp = ya * HEX_40000000;
646           double yaa = ya + temp - temp;
647           double yab = ya - yaa;
648 
649           // recip = 1/y
650           double recip = 1.0/ya;
651           temp = recip * HEX_40000000;
652           double recipa = recip + temp - temp;
653           double recipb = recip - recipa;
654 
655           // Correct for rounding in division
656           recipb += (1.0 - yaa*recipa - yaa*recipb - yab*recipa - yab*recipb) * recip;
657           // Account for yb
658           recipb += -yb * recip * recip;
659 
660           // y = y - 1/y
661           temp = ya - recipa;
662           double ybM = yb - (temp - ya + recipa);
663           double yaM = temp;
664           temp = yaM - recipb;
665           ybM += -(temp - yaM + recipb);
666           yaM = temp;
667           resultM = yaM + ybM;
668           resultM *= 0.5;
669 
670           // y = y + 1/y
671           temp = ya + recipa;
672           double ybP = yb - (temp - ya - recipa);
673           double yaP = temp;
674           temp = yaP + recipb;
675           ybP += -(temp - yaP - recipb);
676           yaP = temp;
677           resultP = yaP + ybP;
678           resultP *= 0.5;
679 
680       } else {
681           expm1(x, hiPrec);
682 
683           final double ya = hiPrec[0] + hiPrec[1];
684           final double yb = -(ya - hiPrec[0] - hiPrec[1]);
685 
686           /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
687           double denom = 1.0 + ya;
688           double denomr = 1.0 / denom;
689           double denomb = -(denom - 1.0 - ya) + yb;
690           double ratio = ya * denomr;
691           double temp = ratio * HEX_40000000;
692           double ra = ratio + temp - temp;
693           double rb = ratio - ra;
694 
695           temp = denom * HEX_40000000;
696           double za = denom + temp - temp;
697           double zb = denom - za;
698 
699           rb += (ya - za*ra - za*rb - zb*ra - zb*rb) * denomr;
700 
701           // Adjust for yb
702           rb += yb*denomr;                        // numerator
703           rb += -ya * denomb * denomr * denomr;   // denominator
704 
705           // y = y - 1/y
706           temp = ya + ra;
707           double ybM = yb - (temp - ya - ra);
708           double yaM = temp;
709           temp = yaM + rb;
710           ybM += -(temp - yaM - rb);
711           yaM = temp;
712           resultM = yaM + ybM;
713           resultM *= 0.5;
714 
715           // y = y + 1/y + 2
716           temp = ya - ra;
717           double ybP = yb - (temp - ya + ra);
718           double yaP = temp;
719           temp = yaP - rb;
720           ybP += -(temp - yaP + rb);
721           yaP = temp;
722           resultP = yaP + ybP + 2;
723           resultP *= 0.5;
724       }
725 
726       if (negate) {
727           resultM = -resultM;
728       }
729 
730       return new SinhCosh(resultM, resultP);
731 
732     }
733 
734     /**
735      * Combined hyperbolic sine and hyperbolic cosine function.
736      *
737      * @param x Argument.
738      * @param <T> the type of the field element
739      * @return [sinh(x), cosh(x)]
740      */
741     public static <T extends CalculusFieldElement<T>> FieldSinhCosh<T> sinhCosh(T x) {
742         return x.sinhCosh();
743     }
744 
745     /** Compute the hyperbolic tangent of a number.
746      * @param x number on which evaluation is done
747      * @return hyperbolic tangent of x
748      */
749     public static double tanh(double x) {
750       boolean negate = false;
751 
752       if (Double.isNaN(x)) {
753           return x;
754       }
755 
756       // tanh[z] = sinh[z] / cosh[z]
757       // = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
758       // = (exp(2x) - 1) / (exp(2x) + 1)
759 
760       // for magnitude > 20, sinh[z] == cosh[z] in double precision
761 
762       if (x > 20.0) {
763           return 1.0;
764       }
765 
766       if (x < -20) {
767           return -1.0;
768       }
769 
770       if (x == 0) {
771           return x;
772       }
773 
774       if (x < 0.0) {
775           x = -x;
776           negate = true;
777       }
778 
779       double result;
780       if (x >= 0.5) {
781           double hiPrec[] = new double[2];
782           // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
783           exp(x*2.0, 0.0, hiPrec);
784 
785           double ya = hiPrec[0] + hiPrec[1];
786           double yb = -(ya - hiPrec[0] - hiPrec[1]);
787 
788           /* Numerator */
789           double na = -1.0 + ya;
790           double nb = -(na + 1.0 - ya);
791           double temp = na + yb;
792           nb += -(temp - na - yb);
793           na = temp;
794 
795           /* Denominator */
796           double da = 1.0 + ya;
797           double db = -(da - 1.0 - ya);
798           temp = da + yb;
799           db += -(temp - da - yb);
800           da = temp;
801 
802           temp = da * HEX_40000000;
803           double daa = da + temp - temp;
804           double dab = da - daa;
805 
806           // ratio = na/da
807           double ratio = na/da;
808           temp = ratio * HEX_40000000;
809           double ratioa = ratio + temp - temp;
810           double ratiob = ratio - ratioa;
811 
812           // Correct for rounding in division
813           ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
814 
815           // Account for nb
816           ratiob += nb / da;
817           // Account for db
818           ratiob += -db * na / da / da;
819 
820           result = ratioa + ratiob;
821       }
822       else {
823           double hiPrec[] = new double[2];
824           // tanh(x) = expm1(2x) / (expm1(2x) + 2)
825           expm1(x*2.0, hiPrec);
826 
827           double ya = hiPrec[0] + hiPrec[1];
828           double yb = -(ya - hiPrec[0] - hiPrec[1]);
829 
830           /* Numerator */
831           double na = ya;
832           double nb = yb;
833 
834           /* Denominator */
835           double da = 2.0 + ya;
836           double db = -(da - 2.0 - ya);
837           double temp = da + yb;
838           db += -(temp - da - yb);
839           da = temp;
840 
841           temp = da * HEX_40000000;
842           double daa = da + temp - temp;
843           double dab = da - daa;
844 
845           // ratio = na/da
846           double ratio = na/da;
847           temp = ratio * HEX_40000000;
848           double ratioa = ratio + temp - temp;
849           double ratiob = ratio - ratioa;
850 
851           // Correct for rounding in division
852           ratiob += (na - daa*ratioa - daa*ratiob - dab*ratioa - dab*ratiob) / da;
853 
854           // Account for nb
855           ratiob += nb / da;
856           // Account for db
857           ratiob += -db * na / da / da;
858 
859           result = ratioa + ratiob;
860       }
861 
862       if (negate) {
863           result = -result;
864       }
865 
866       return result;
867     }
868 
869     /** Compute the inverse hyperbolic cosine of a number.
870      * @param a number on which evaluation is done
871      * @return inverse hyperbolic cosine of a
872      */
873     public static double acosh(final double a) {
874         return FastMath.log(a + FastMath.sqrt(a * a - 1));
875     }
876 
877     /** Compute the inverse hyperbolic sine of a number.
878      * @param a number on which evaluation is done
879      * @return inverse hyperbolic sine of a
880      */
881     public static double asinh(double a) {
882         boolean negative = false;
883         if (a < 0) {
884             negative = true;
885             a = -a;
886         }
887 
888         double absAsinh;
889         if (a > 0.167) {
890             absAsinh = FastMath.log(FastMath.sqrt(a * a + 1) + a);
891         } else {
892             final double a2 = a * a;
893             if (a > 0.097) {
894                 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);
895             } else if (a > 0.036) {
896                 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);
897             } else if (a > 0.0036) {
898                 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);
899             } else {
900                 absAsinh = a * (1 - a2 * (F_1_3 - a2 * F_1_5 * F_3_4) * F_1_2);
901             }
902         }
903 
904         return negative ? -absAsinh : absAsinh;
905     }
906 
907     /** Compute the inverse hyperbolic tangent of a number.
908      * @param a number on which evaluation is done
909      * @return inverse hyperbolic tangent of a
910      */
911     public static double atanh(double a) {
912         boolean negative = false;
913         if (a < 0) {
914             negative = true;
915             a = -a;
916         }
917 
918         double absAtanh;
919         if (a > 0.15) {
920             absAtanh = 0.5 * FastMath.log((1 + a) / (1 - a));
921         } else {
922             final double a2 = a * a;
923             if (a > 0.087) {
924                 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))))))));
925             } else if (a > 0.031) {
926                 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))))));
927             } else if (a > 0.003) {
928                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * (F_1_5 + a2 * (F_1_7 + a2 * F_1_9))));
929             } else {
930                 absAtanh = a * (1 + a2 * (F_1_3 + a2 * F_1_5));
931             }
932         }
933 
934         return negative ? -absAtanh : absAtanh;
935     }
936 
937     /** Compute the signum of a number.
938      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
939      * @param a number on which evaluation is done
940      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
941      */
942     public static double signum(final double a) {
943         return (a < 0.0) ? -1.0 : ((a > 0.0) ? 1.0 : a); // return +0.0/-0.0/NaN depending on a
944     }
945 
946     /** Compute the signum of a number.
947      * The signum is -1 for negative numbers, +1 for positive numbers and 0 otherwise
948      * @param a number on which evaluation is done
949      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
950      */
951     public static float signum(final float a) {
952         return (a < 0.0f) ? -1.0f : ((a > 0.0f) ? 1.0f : a); // return +0.0/-0.0/NaN depending on a
953     }
954 
955     /** Compute next number towards positive infinity.
956      * @param a number to which neighbor should be computed
957      * @return neighbor of a towards positive infinity
958      */
959     public static double nextUp(final double a) {
960         return nextAfter(a, Double.POSITIVE_INFINITY);
961     }
962 
963     /** Compute next number towards positive infinity.
964      * @param a number to which neighbor should be computed
965      * @return neighbor of a towards positive infinity
966      */
967     public static float nextUp(final float a) {
968         return nextAfter(a, Float.POSITIVE_INFINITY);
969     }
970 
971     /** Compute next number towards negative infinity.
972      * @param a number to which neighbor should be computed
973      * @return neighbor of a towards negative infinity
974      */
975     public static double nextDown(final double a) {
976         return nextAfter(a, Double.NEGATIVE_INFINITY);
977     }
978 
979     /** Compute next number towards negative infinity.
980      * @param a number to which neighbor should be computed
981      * @return neighbor of a towards negative infinity
982      */
983     public static float nextDown(final float a) {
984         return nextAfter(a, Float.NEGATIVE_INFINITY);
985     }
986 
987     /** Clamp a value within an interval.
988      * @param value value to clamp
989      * @param inf lower bound of the clamping interval
990      * @param sup upper bound of the clamping interval
991      * @return value clamped within [inf; sup], or value if already within bounds.
992      * @since 3.0
993      */
994     public static int clamp(final int value, final int inf, final int sup) {
995         return max(inf, min(value, sup));
996     }
997 
998     /** Clamp a value within an interval.
999      * @param value value to clamp
1000      * @param inf lower bound of the clamping interval
1001      * @param sup upper bound of the clamping interval
1002      * @return value clamped within [inf; sup], or value if already within bounds.
1003      * @since 3.0
1004      */
1005     public static long clamp(final long value, final long inf, final long sup) {
1006         return max(inf, min(value, sup));
1007     }
1008 
1009     /** Clamp a value within an interval.
1010      * @param value value to clamp
1011      * @param inf lower bound of the clamping interval
1012      * @param sup upper bound of the clamping interval
1013      * @return value clamped within [inf; sup], or value if already within bounds.
1014      * @since 3.0
1015      */
1016     public static int clamp(final long value, final int inf, final int sup) {
1017         return (int) max(inf, min(value, sup));
1018     }
1019 
1020     /** Clamp a value within an interval.
1021      * <p>
1022      * This method assumes -0.0 is below +0.0
1023      * </p>
1024      * @param value value to clamp
1025      * @param inf lower bound of the clamping interval
1026      * @param sup upper bound of the clamping interval
1027      * @return value clamped within [inf; sup], or value if already within bounds.
1028      * @since 3.0
1029      */
1030     public static float clamp(final float value, final float inf, final float sup) {
1031         return max(inf, min(value, sup));
1032     }
1033 
1034     /** Clamp a value within an interval.
1035      * <p>
1036      * This method assumes -0.0 is below +0.0
1037      * </p>
1038      * @param value value to clamp
1039      * @param inf lower bound of the clamping interval
1040      * @param sup upper bound of the clamping interval
1041      * @return value clamped within [inf; sup], or value if already within bounds.
1042      * @since 3.0
1043      */
1044     public static double clamp(final double value, final double inf, final double sup) {
1045         return max(inf, min(value, sup));
1046     }
1047 
1048     /** Returns a pseudo-random number between 0.0 and 1.0.
1049      * <p><b>Note:</b> this implementation currently delegates to {@link Math#random}
1050      * @return a random number between 0.0 and 1.0
1051      */
1052     public static double random() {
1053         return Math.random();
1054     }
1055 
1056     /**
1057      * Exponential function.
1058      *
1059      * Computes exp(x), function result is nearly rounded.   It will be correctly
1060      * rounded to the theoretical value for 99.9% of input values, otherwise it will
1061      * have a 1 ULP error.
1062      *
1063      * Method:
1064      *    Lookup intVal = exp(int(x))
1065      *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
1066      *    Compute z as the exponential of the remaining bits by a polynomial minus one
1067      *    exp(x) = intVal * fracVal * (1 + z)
1068      *
1069      * Accuracy:
1070      *    Calculation is done with 63 bits of precision, so result should be correctly
1071      *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
1072      *
1073      * @param x   a double
1074      * @return double e<sup>x</sup>
1075      */
1076     public static double exp(double x) {
1077         return exp(x, 0.0, null);
1078     }
1079 
1080     /**
1081      * Internal helper method for exponential function.
1082      * @param x original argument of the exponential function
1083      * @param extra extra bits of precision on input (To Be Confirmed)
1084      * @param hiPrec extra bits of precision on output (To Be Confirmed)
1085      * @return exp(x)
1086      */
1087     private static double exp(double x, double extra, double[] hiPrec) {
1088         double intPartA;
1089         double intPartB;
1090         int intVal = (int) x;
1091 
1092         /* Lookup exp(floor(x)).
1093          * intPartA will have the upper 22 bits, intPartB will have the lower
1094          * 52 bits.
1095          */
1096         if (x < 0.0) {
1097 
1098             // We don't check against intVal here as conversion of large negative double values
1099             // may be affected by a JIT bug. Subsequent comparisons can safely use intVal
1100             if (x < -746d) {
1101                 if (hiPrec != null) {
1102                     hiPrec[0] = 0.0;
1103                     hiPrec[1] = 0.0;
1104                 }
1105                 return 0.0;
1106             }
1107 
1108             if (intVal < -709) {
1109                 /* This will produce a subnormal output */
1110                 final double result = exp(x+40.19140625, extra, hiPrec) / 285040095144011776.0;
1111                 if (hiPrec != null) {
1112                     hiPrec[0] /= 285040095144011776.0;
1113                     hiPrec[1] /= 285040095144011776.0;
1114                 }
1115                 return result;
1116             }
1117 
1118             if (intVal == -709) {
1119                 /* exp(1.494140625) is nearly a machine number... */
1120                 final double result = exp(x+1.494140625, extra, hiPrec) / 4.455505956692756620;
1121                 if (hiPrec != null) {
1122                     hiPrec[0] /= 4.455505956692756620;
1123                     hiPrec[1] /= 4.455505956692756620;
1124                 }
1125                 return result;
1126             }
1127 
1128             intVal--;
1129 
1130         } else {
1131             if (intVal > 709) {
1132                 if (hiPrec != null) {
1133                     hiPrec[0] = Double.POSITIVE_INFINITY;
1134                     hiPrec[1] = 0.0;
1135                 }
1136                 return Double.POSITIVE_INFINITY;
1137             }
1138 
1139         }
1140 
1141         intPartA = ExpIntTable.EXP_INT_TABLE_A[EXP_INT_TABLE_MAX_INDEX+intVal];
1142         intPartB = ExpIntTable.EXP_INT_TABLE_B[EXP_INT_TABLE_MAX_INDEX+intVal];
1143 
1144         /* Get the fractional part of x, find the greatest multiple of 2^-10 less than
1145          * x and look up the exp function of it.
1146          * fracPartA will have the upper 22 bits, fracPartB the lower 52 bits.
1147          */
1148         final int intFrac = (int) ((x - intVal) * 1024.0);
1149         final double fracPartA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac];
1150         final double fracPartB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1151 
1152         /* epsilon is the difference in x from the nearest multiple of 2^-10.  It
1153          * has a value in the range 0 <= epsilon < 2^-10.
1154          * Do the subtraction from x as the last step to avoid possible loss of precision.
1155          */
1156         final double epsilon = x - (intVal + intFrac / 1024.0);
1157 
1158         /* Compute z = exp(epsilon) - 1.0 via a minimax polynomial.  z has
1159        full double precision (52 bits).  Since z < 2^-10, we will have
1160        62 bits of precision when combined with the constant 1.  This will be
1161        used in the last addition below to get proper rounding. */
1162 
1163         /* Remez generated polynomial.  Converges on the interval [0, 2^-10], error
1164        is less than 0.5 ULP */
1165         double z = 0.04168701738764507;
1166         z = z * epsilon + 0.1666666505023083;
1167         z = z * epsilon + 0.5000000000042687;
1168         z = z * epsilon + 1.0;
1169         z = z * epsilon + -3.940510424527919E-20;
1170 
1171         /* Compute (intPartA+intPartB) * (fracPartA+fracPartB) by binomial
1172        expansion.
1173        tempA is exact since intPartA and intPartB only have 22 bits each.
1174        tempB will have 52 bits of precision.
1175          */
1176         double tempA = intPartA * fracPartA;
1177         double tempB = intPartA * fracPartB + intPartB * fracPartA + intPartB * fracPartB;
1178 
1179         /* Compute the result.  (1+z)(tempA+tempB).  Order of operations is
1180        important.  For accuracy add by increasing size.  tempA is exact and
1181        much larger than the others.  If there are extra bits specified from the
1182        pow() function, use them. */
1183         final double tempC = tempB + tempA;
1184 
1185         // If tempC is positive infinite, the evaluation below could result in NaN,
1186         // because z could be negative at the same time.
1187         if (tempC == Double.POSITIVE_INFINITY) {
1188             if (hiPrec != null) {
1189                 hiPrec[0] = Double.POSITIVE_INFINITY;
1190                 hiPrec[1] = 0.0;
1191             }
1192             return Double.POSITIVE_INFINITY;
1193         }
1194 
1195         final double result;
1196         if (extra != 0.0) {
1197             result = tempC*extra*z + tempC*extra + tempC*z + tempB + tempA;
1198         } else {
1199             result = tempC*z + tempB + tempA;
1200         }
1201 
1202         if (hiPrec != null) {
1203             // If requesting high precision
1204             hiPrec[0] = tempA;
1205             hiPrec[1] = tempC*extra*z + tempC*extra + tempC*z + tempB;
1206         }
1207 
1208         return result;
1209     }
1210 
1211     /** Compute exp(x) - 1
1212      * @param x number to compute shifted exponential
1213      * @return exp(x) - 1
1214      */
1215     public static double expm1(double x) {
1216       return expm1(x, null);
1217     }
1218 
1219     /** Internal helper method for expm1
1220      * @param x number to compute shifted exponential
1221      * @param hiPrecOut receive high precision result for -1.0 < x < 1.0
1222      * @return exp(x) - 1
1223      */
1224     private static double expm1(double x, double hiPrecOut[]) {
1225         if (Double.isNaN(x) || x == 0.0) { // NaN or zero
1226             return x;
1227         }
1228 
1229         if (x <= -1.0 || x >= 1.0) {
1230             // If not between +/- 1.0
1231             //return exp(x) - 1.0;
1232             double hiPrec[] = new double[2];
1233             exp(x, 0.0, hiPrec);
1234             if (x > 0.0) {
1235                 return -1.0 + hiPrec[0] + hiPrec[1];
1236             } else {
1237                 final double ra = -1.0 + hiPrec[0];
1238                 double rb = -(ra + 1.0 - hiPrec[0]);
1239                 rb += hiPrec[1];
1240                 return ra + rb;
1241             }
1242         }
1243 
1244         double baseA;
1245         double baseB;
1246         double epsilon;
1247         boolean negative = false;
1248 
1249         if (x < 0.0) {
1250             x = -x;
1251             negative = true;
1252         }
1253 
1254         {
1255             int intFrac = (int) (x * 1024.0);
1256             double tempA = ExpFracTable.EXP_FRAC_TABLE_A[intFrac] - 1.0;
1257             double tempB = ExpFracTable.EXP_FRAC_TABLE_B[intFrac];
1258 
1259             double temp = tempA + tempB;
1260             tempB = -(temp - tempA - tempB);
1261             tempA = temp;
1262 
1263             temp = tempA * HEX_40000000;
1264             baseA = tempA + temp - temp;
1265             baseB = tempB + (tempA - baseA);
1266 
1267             epsilon = x - intFrac/1024.0;
1268         }
1269 
1270 
1271         /* Compute expm1(epsilon) */
1272         double zb = 0.008336750013465571;
1273         zb = zb * epsilon + 0.041666663879186654;
1274         zb = zb * epsilon + 0.16666666666745392;
1275         zb = zb * epsilon + 0.49999999999999994;
1276         zb *= epsilon;
1277         zb *= epsilon;
1278 
1279         double za = epsilon;
1280         double temp = za + zb;
1281         zb = -(temp - za - zb);
1282         za = temp;
1283 
1284         temp = za * HEX_40000000;
1285         temp = za + temp - temp;
1286         zb += za - temp;
1287         za = temp;
1288 
1289         /* Combine the parts.   expm1(a+b) = expm1(a) + expm1(b) + expm1(a)*expm1(b) */
1290         double ya = za * baseA;
1291         //double yb = za*baseB + zb*baseA + zb*baseB;
1292         temp = ya + za * baseB;
1293         double yb = -(temp - ya - za * baseB);
1294         ya = temp;
1295 
1296         temp = ya + zb * baseA;
1297         yb += -(temp - ya - zb * baseA);
1298         ya = temp;
1299 
1300         temp = ya + zb * baseB;
1301         yb += -(temp - ya - zb*baseB);
1302         ya = temp;
1303 
1304         //ya = ya + za + baseA;
1305         //yb = yb + zb + baseB;
1306         temp = ya + baseA;
1307         yb += -(temp - baseA - ya);
1308         ya = temp;
1309 
1310         temp = ya + za;
1311         //yb += (ya > za) ? -(temp - ya - za) : -(temp - za - ya);
1312         yb += -(temp - ya - za);
1313         ya = temp;
1314 
1315         temp = ya + baseB;
1316         //yb += (ya > baseB) ? -(temp - ya - baseB) : -(temp - baseB - ya);
1317         yb += -(temp - ya - baseB);
1318         ya = temp;
1319 
1320         temp = ya + zb;
1321         //yb += (ya > zb) ? -(temp - ya - zb) : -(temp - zb - ya);
1322         yb += -(temp - ya - zb);
1323         ya = temp;
1324 
1325         if (negative) {
1326             /* Compute expm1(-x) = -expm1(x) / (expm1(x) + 1) */
1327             double denom = 1.0 + ya;
1328             double denomr = 1.0 / denom;
1329             double denomb = -(denom - 1.0 - ya) + yb;
1330             double ratio = ya * denomr;
1331             temp = ratio * HEX_40000000;
1332             final double ra = ratio + temp - temp;
1333             double rb = ratio - ra;
1334 
1335             temp = denom * HEX_40000000;
1336             za = denom + temp - temp;
1337             zb = denom - za;
1338 
1339             rb += (ya - za * ra - za * rb - zb * ra - zb * rb) * denomr;
1340 
1341             // f(x) = x/1+x
1342             // Compute f'(x)
1343             // Product rule:  d(uv) = du*v + u*dv
1344             // Chain rule:  d(f(g(x)) = f'(g(x))*f(g'(x))
1345             // d(1/x) = -1/(x*x)
1346             // d(1/1+x) = -1/( (1+x)^2) *  1 =  -1/((1+x)*(1+x))
1347             // d(x/1+x) = -x/((1+x)(1+x)) + 1/1+x = 1 / ((1+x)(1+x))
1348 
1349             // Adjust for yb
1350             rb += yb * denomr;                      // numerator
1351             rb += -ya * denomb * denomr * denomr;   // denominator
1352 
1353             // negate
1354             ya = -ra;
1355             yb = -rb;
1356         }
1357 
1358         if (hiPrecOut != null) {
1359             hiPrecOut[0] = ya;
1360             hiPrecOut[1] = yb;
1361         }
1362 
1363         return ya + yb;
1364     }
1365 
1366     /**
1367      * Natural logarithm.
1368      *
1369      * @param x   a double
1370      * @return log(x)
1371      */
1372     public static double log(final double x) {
1373         return log(x, null);
1374     }
1375 
1376     /**
1377      * Internal helper method for natural logarithm function.
1378      * @param x original argument of the natural logarithm function
1379      * @param hiPrec extra bits of precision on output (To Be Confirmed)
1380      * @return log(x)
1381      */
1382     private static double log(final double x, final double[] hiPrec) {
1383         if (x==0) { // Handle special case of +0/-0
1384             return Double.NEGATIVE_INFINITY;
1385         }
1386         long bits = Double.doubleToRawLongBits(x);
1387 
1388         /* Handle special cases of negative input, and NaN */
1389         if (((bits & 0x8000000000000000L) != 0 || Double.isNaN(x)) && x != 0.0) {
1390             if (hiPrec != null) {
1391                 hiPrec[0] = Double.NaN;
1392             }
1393 
1394             return Double.NaN;
1395         }
1396 
1397         /* Handle special cases of Positive infinity. */
1398         if (x == Double.POSITIVE_INFINITY) {
1399             if (hiPrec != null) {
1400                 hiPrec[0] = Double.POSITIVE_INFINITY;
1401             }
1402 
1403             return Double.POSITIVE_INFINITY;
1404         }
1405 
1406         /* Extract the exponent */
1407         int exp = (int)(bits >> 52)-1023;
1408 
1409         if ((bits & 0x7ff0000000000000L) == 0) {
1410             // Subnormal!
1411             if (x == 0) {
1412                 // Zero
1413                 if (hiPrec != null) {
1414                     hiPrec[0] = Double.NEGATIVE_INFINITY;
1415                 }
1416 
1417                 return Double.NEGATIVE_INFINITY;
1418             }
1419 
1420             /* Normalize the subnormal number. */
1421             bits <<= 1;
1422             while ( (bits & 0x0010000000000000L) == 0) {
1423                 --exp;
1424                 bits <<= 1;
1425             }
1426         }
1427 
1428 
1429         if ((exp == -1 || exp == 0) && x < 1.01 && x > 0.99 && hiPrec == null) {
1430             /* The normal method doesn't work well in the range [0.99, 1.01], so call do a straight
1431            polynomial expansion in higer precision. */
1432 
1433             /* Compute x - 1.0 and split it */
1434             double xa = x - 1.0;
1435             double tmp = xa * HEX_40000000;
1436             double aa = xa + tmp - tmp;
1437             double ab = xa - aa;
1438             xa = aa;
1439             double xb = ab;
1440 
1441             final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
1442             double ya = lnCoef_last[0];
1443             double yb = lnCoef_last[1];
1444 
1445             for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
1446                 /* Multiply a = y * x */
1447                 aa = ya * xa;
1448                 ab = ya * xb + yb * xa + yb * xb;
1449                 /* split, so now y = a */
1450                 tmp = aa * HEX_40000000;
1451                 ya = aa + tmp - tmp;
1452                 yb = aa - ya + ab;
1453 
1454                 /* Add  a = y + lnQuickCoef */
1455                 final double[] lnCoef_i = LN_QUICK_COEF[i];
1456                 aa = ya + lnCoef_i[0];
1457                 ab = yb + lnCoef_i[1];
1458                 /* Split y = a */
1459                 tmp = aa * HEX_40000000;
1460                 ya = aa + tmp - tmp;
1461                 yb = aa - ya + ab;
1462             }
1463 
1464             /* Multiply a = y * x */
1465             aa = ya * xa;
1466             ab = ya * xb + yb * xa + yb * xb;
1467             /* split, so now y = a */
1468             tmp = aa * HEX_40000000;
1469             ya = aa + tmp - tmp;
1470             yb = aa - ya + ab;
1471 
1472             return ya + yb;
1473         }
1474 
1475         // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
1476         final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
1477 
1478         /*
1479     double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1480 
1481     epsilon -= 1.0;
1482          */
1483 
1484         // y is the most significant 10 bits of the mantissa
1485         //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
1486         //double epsilon = (x - y) / y;
1487         final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
1488 
1489         double lnza;
1490         double lnzb = 0.0;
1491 
1492         if (hiPrec != null) {
1493             /* split epsilon -> x */
1494             double tmp = epsilon * HEX_40000000;
1495             double aa = epsilon + tmp - tmp;
1496             double ab = epsilon - aa;
1497             double xa = aa;
1498             double xb = ab;
1499 
1500             /* Need a more accurate epsilon, so adjust the division. */
1501             final double numer = bits & 0x3ffffffffffL;
1502             final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
1503             aa = numer - xa*denom - xb * denom;
1504             xb += aa / denom;
1505 
1506             /* Remez polynomial evaluation */
1507             final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
1508             double ya = lnCoef_last[0];
1509             double yb = lnCoef_last[1];
1510 
1511             for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
1512                 /* Multiply a = y * x */
1513                 aa = ya * xa;
1514                 ab = ya * xb + yb * xa + yb * xb;
1515                 /* split, so now y = a */
1516                 tmp = aa * HEX_40000000;
1517                 ya = aa + tmp - tmp;
1518                 yb = aa - ya + ab;
1519 
1520                 /* Add  a = y + lnHiPrecCoef */
1521                 final double[] lnCoef_i = LN_HI_PREC_COEF[i];
1522                 aa = ya + lnCoef_i[0];
1523                 ab = yb + lnCoef_i[1];
1524                 /* Split y = a */
1525                 tmp = aa * HEX_40000000;
1526                 ya = aa + tmp - tmp;
1527                 yb = aa - ya + ab;
1528             }
1529 
1530             /* Multiply a = y * x */
1531             aa = ya * xa;
1532             ab = ya * xb + yb * xa + yb * xb;
1533 
1534             /* split, so now lnz = a */
1535             /*
1536       tmp = aa * 1073741824.0;
1537       lnza = aa + tmp - tmp;
1538       lnzb = aa - lnza + ab;
1539              */
1540             lnza = aa + ab;
1541             lnzb = -(lnza - aa - ab);
1542         } else {
1543             /* High precision not required.  Eval Remez polynomial
1544          using standard double precision */
1545             lnza = -0.16624882440418567;
1546             lnza = lnza * epsilon + 0.19999954120254515;
1547             lnza = lnza * epsilon + -0.2499999997677497;
1548             lnza = lnza * epsilon + 0.3333333333332802;
1549             lnza = lnza * epsilon + -0.5;
1550             lnza = lnza * epsilon + 1.0;
1551             lnza *= epsilon;
1552         }
1553 
1554         /* Relative sizes:
1555          * lnzb     [0, 2.33E-10]
1556          * lnm[1]   [0, 1.17E-7]
1557          * ln2B*exp [0, 1.12E-4]
1558          * lnza      [0, 9.7E-4]
1559          * lnm[0]   [0, 0.692]
1560          * ln2A*exp [0, 709]
1561          */
1562 
1563         /* Compute the following sum:
1564          * lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1565          */
1566 
1567         //return lnzb + lnm[1] + ln2B*exp + lnza + lnm[0] + ln2A*exp;
1568         double a = LN_2_A*exp;
1569         double b = 0.0;
1570         double c = a+lnm[0];
1571         double d = -(c-a-lnm[0]);
1572         a = c;
1573         b += d;
1574 
1575         c = a + lnza;
1576         d = -(c - a - lnza);
1577         a = c;
1578         b += d;
1579 
1580         c = a + LN_2_B*exp;
1581         d = -(c - a - LN_2_B*exp);
1582         a = c;
1583         b += d;
1584 
1585         c = a + lnm[1];
1586         d = -(c - a - lnm[1]);
1587         a = c;
1588         b += d;
1589 
1590         c = a + lnzb;
1591         d = -(c - a - lnzb);
1592         a = c;
1593         b += d;
1594 
1595         if (hiPrec != null) {
1596             hiPrec[0] = a;
1597             hiPrec[1] = b;
1598         }
1599 
1600         return a + b;
1601     }
1602 
1603     /**
1604      * Computes log(1 + x).
1605      *
1606      * @param x Number.
1607      * @return {@code log(1 + x)}.
1608      */
1609     public static double log1p(final double x) {
1610         if (x == -1) {
1611             return Double.NEGATIVE_INFINITY;
1612         }
1613 
1614         if (x == Double.POSITIVE_INFINITY) {
1615             return Double.POSITIVE_INFINITY;
1616         }
1617 
1618         if (x > 1e-6 ||
1619             x < -1e-6) {
1620             final double xpa = 1 + x;
1621             final double xpb = -(xpa - 1 - x);
1622 
1623             final double[] hiPrec = new double[2];
1624             final double lores = log(xpa, hiPrec);
1625             if (Double.isInfinite(lores)) { // Don't allow this to be converted to NaN
1626                 return lores;
1627             }
1628 
1629             // Do a taylor series expansion around xpa:
1630             //   f(x+y) = f(x) + f'(x) y + f''(x)/2 y^2
1631             final double fx1 = xpb / xpa;
1632             final double epsilon = 0.5 * fx1 + 1;
1633             return epsilon * fx1 + hiPrec[1] + hiPrec[0];
1634         } else {
1635             // Value is small |x| < 1e6, do a Taylor series centered on 1.
1636             final double y = (x * F_1_3 - F_1_2) * x + 1;
1637             return y * x;
1638         }
1639     }
1640 
1641     /** Compute the base 10 logarithm.
1642      * @param x a number
1643      * @return log10(x)
1644      */
1645     public static double log10(final double x) {
1646         final double hiPrec[] = new double[2];
1647 
1648         final double lores = log(x, hiPrec);
1649         if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
1650             return lores;
1651         }
1652 
1653         final double tmp = hiPrec[0] * HEX_40000000;
1654         final double lna = hiPrec[0] + tmp - tmp;
1655         final double lnb = hiPrec[0] - lna + hiPrec[1];
1656 
1657         final double rln10a = 0.4342944622039795;
1658         final double rln10b = 1.9699272335463627E-8;
1659 
1660         return rln10b * lnb + rln10b * lna + rln10a * lnb + rln10a * lna;
1661     }
1662 
1663     /**
1664      * Computes the <a href="http://mathworld.wolfram.com/Logarithm.html">
1665      * logarithm</a> in a given base.
1666      *
1667      * Returns {@code NaN} if either argument is negative.
1668      * If {@code base} is 0 and {@code x} is positive, 0 is returned.
1669      * If {@code base} is positive and {@code x} is 0,
1670      * {@code Double.NEGATIVE_INFINITY} is returned.
1671      * If both arguments are 0, the result is {@code NaN}.
1672      *
1673      * @param base Base of the logarithm, must be greater than 0.
1674      * @param x Argument, must be greater than 0.
1675      * @return the value of the logarithm, i.e. the number {@code y} such that
1676      * <code>base<sup>y</sup> = x</code>.
1677      */
1678     public static double log(double base, double x) {
1679         return log(x) / log(base);
1680     }
1681 
1682     /**
1683      * Power function.  Compute x^y.
1684      *
1685      * @param x   a double
1686      * @param y   a double
1687      * @return double
1688      */
1689     public static double pow(final double x, final double y) {
1690 
1691         if (y == 0) {
1692             // y = -0 or y = +0
1693             return 1.0;
1694         } else {
1695 
1696             final long yBits        = Double.doubleToRawLongBits(y);
1697             final int  yRawExp      = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
1698             final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
1699             final long xBits        = Double.doubleToRawLongBits(x);
1700             final int  xRawExp      = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
1701             final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;
1702 
1703             if (yRawExp > 1085) {
1704                 // y is either a very large integral value that does not fit in a long or it is a special number
1705 
1706                 if ((yRawExp == 2047 && yRawMantissa != 0) ||
1707                     (xRawExp == 2047 && xRawMantissa != 0)) {
1708                     // NaN
1709                     return Double.NaN;
1710                 } else if (xRawExp == 1023 && xRawMantissa == 0) {
1711                     // x = -1.0 or x = +1.0
1712                     if (yRawExp == 2047) {
1713                         // y is infinite
1714                         return Double.NaN;
1715                     } else {
1716                         // y is a large even integer
1717                         return 1.0;
1718                     }
1719                 } else {
1720                     // the absolute value of x is either greater or smaller than 1.0
1721 
1722                     // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
1723                     // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
1724                     // accuracy, at this magnitude it behaves just like infinity with regards to x
1725                     if ((y > 0) ^ (xRawExp < 1023)) {
1726                         // either y = +infinity (or large engouh) and abs(x) > 1.0
1727                         // or     y = -infinity (or large engouh) and abs(x) < 1.0
1728                         return Double.POSITIVE_INFINITY;
1729                     } else {
1730                         // either y = +infinity (or large engouh) and abs(x) < 1.0
1731                         // or     y = -infinity (or large engouh) and abs(x) > 1.0
1732                         return +0.0;
1733                     }
1734                 }
1735 
1736             } else {
1737                 // y is a regular non-zero number
1738 
1739                 if (yRawExp >= 1023) {
1740                     // y may be an integral value, which should be handled specifically
1741                     final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
1742                     if (yRawExp < 1075) {
1743                         // normal number with negative shift that may have a fractional part
1744                         final long integralMask = (-1L) << (1075 - yRawExp);
1745                         if ((yFullMantissa & integralMask) == yFullMantissa) {
1746                             // all fractional bits are 0, the number is really integral
1747                             final long l = yFullMantissa >> (1075 - yRawExp);
1748                             return FastMath.pow(x, (y < 0) ? -l : l);
1749                         }
1750                     } else {
1751                         // normal number with positive shift, always an integral value
1752                         // we know it fits in a primitive long because yRawExp > 1085 has been handled above
1753                         final long l =  yFullMantissa << (yRawExp - 1075);
1754                         return FastMath.pow(x, (y < 0) ? -l : l);
1755                     }
1756                 }
1757 
1758                 // y is a non-integral value
1759 
1760                 if (x == 0) {
1761                     // x = -0 or x = +0
1762                     // the integer powers have already been handled above
1763                     return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
1764                 } else if (xRawExp == 2047) {
1765                     if (xRawMantissa == 0) {
1766                         // x = -infinity or x = +infinity
1767                         return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
1768                     } else {
1769                         // NaN
1770                         return Double.NaN;
1771                     }
1772                 } else if (x < 0) {
1773                     // the integer powers have already been handled above
1774                     return Double.NaN;
1775                 } else {
1776 
1777                     // this is the general case, for regular fractional numbers x and y
1778 
1779                     // Split y into ya and yb such that y = ya+yb
1780                     final double tmp = y * HEX_40000000;
1781                     final double ya = (y + tmp) - tmp;
1782                     final double yb = y - ya;
1783 
1784                     /* Compute ln(x) */
1785                     final double lns[] = new double[2];
1786                     final double lores = log(x, lns);
1787                     if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
1788                         return lores;
1789                     }
1790 
1791                     double lna = lns[0];
1792                     double lnb = lns[1];
1793 
1794                     /* resplit lns */
1795                     final double tmp1 = lna * HEX_40000000;
1796                     final double tmp2 = (lna + tmp1) - tmp1;
1797                     lnb += lna - tmp2;
1798                     lna = tmp2;
1799 
1800                     // y*ln(x) = (aa+ab)
1801                     final double aa = lna * ya;
1802                     final double ab = lna * yb + lnb * ya + lnb * yb;
1803 
1804                     lna = aa+ab;
1805                     lnb = -(lna - aa - ab);
1806 
1807                     double z = 1.0 / 120.0;
1808                     z = z * lnb + (1.0 / 24.0);
1809                     z = z * lnb + (1.0 / 6.0);
1810                     z = z * lnb + 0.5;
1811                     z = z * lnb + 1.0;
1812                     z *= lnb;
1813 
1814                     return exp(lna, z, null);
1815 
1816                 }
1817             }
1818 
1819         }
1820 
1821     }
1822 
1823     /**
1824      * Raise a double to an int power.
1825      *
1826      * @param d Number to raise.
1827      * @param e Exponent.
1828      * @return d<sup>e</sup>
1829      */
1830     public static double pow(double d, int e) {
1831         return pow(d, (long) e);
1832     }
1833 
1834     /**
1835      * Raise a double to a long power.
1836      *
1837      * @param d Number to raise.
1838      * @param e Exponent.
1839      * @return d<sup>e</sup>
1840      */
1841     public static double pow(double d, long e) {
1842         if (e == 0) {
1843             return 1.0;
1844         } else if (e > 0) {
1845             return new Split(d).pow(e).full;
1846         } else {
1847             return new Split(d).reciprocal().pow(-e).full;
1848         }
1849     }
1850 
1851     /** Class operator on double numbers split into one 26 bits number and one 27 bits number. */
1852     private static class Split {
1853 
1854         /** Split version of NaN. */
1855         public static final Split NAN = new Split(Double.NaN, 0);
1856 
1857         /** Split version of positive infinity. */
1858         public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
1859 
1860         /** Split version of negative infinity. */
1861         public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
1862 
1863         /** Full number. */
1864         private final double full;
1865 
1866         /** High order bits. */
1867         private final double high;
1868 
1869         /** Low order bits. */
1870         private final double low;
1871 
1872         /** Simple constructor.
1873          * @param x number to split
1874          */
1875         Split(final double x) {
1876             full = x;
1877             high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
1878             low  = x - high;
1879         }
1880 
1881         /** Simple constructor.
1882          * @param high high order bits
1883          * @param low low order bits
1884          */
1885         Split(final double high, final double low) {
1886             this(high == 0.0 ? (low == 0.0 && Double.doubleToRawLongBits(high) == Long.MIN_VALUE /* negative zero */ ? -0.0 : low) : high + low, high, low);
1887         }
1888 
1889         /** Simple constructor.
1890          * @param full full number
1891          * @param high high order bits
1892          * @param low low order bits
1893          */
1894         Split(final double full, final double high, final double low) {
1895             this.full = full;
1896             this.high = high;
1897             this.low  = low;
1898         }
1899 
1900         /** Multiply the instance by another one.
1901          * @param b other instance to multiply by
1902          * @return product
1903          */
1904         public Split multiply(final Split b) {
1905             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1906             final Split  mulBasic  = new Split(full * b.full);
1907             final double mulError  = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
1908             return new Split(mulBasic.high, mulBasic.low + mulError);
1909         }
1910 
1911         /** Compute the reciprocal of the instance.
1912          * @return reciprocal of the instance
1913          */
1914         public Split reciprocal() {
1915 
1916             final double approximateInv = 1.0 / full;
1917             final Split  splitInv       = new Split(approximateInv);
1918 
1919             // if 1.0/d were computed perfectly, remultiplying it by d should give 1.0
1920             // we want to estimate the error so we can fix the low order bits of approximateInvLow
1921             // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
1922             final Split product = multiply(splitInv);
1923             final double error  = (product.high - 1) + product.low;
1924 
1925             // better accuracy estimate of reciprocal
1926             return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
1927 
1928         }
1929 
1930         /** Computes this^e.
1931          * @param e exponent (beware, here it MUST be > 0; the only exclusion is Long.MIN_VALUE)
1932          * @return d^e, split in high and low bits
1933          */
1934         private Split pow(final long e) {
1935 
1936             // prepare result
1937             Split result = new Split(1);
1938 
1939             // d^(2p)
1940             Split d2p = new Split(full, high, low);
1941 
1942             for (long p = e; p != 0; p >>>= 1) {
1943 
1944                 if ((p & 0x1) != 0) {
1945                     // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
1946                     result = result.multiply(d2p);
1947                 }
1948 
1949                 // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
1950                 d2p = d2p.multiply(d2p);
1951 
1952             }
1953 
1954             if (Double.isNaN(result.full)) {
1955                 if (Double.isNaN(full)) {
1956                     return Split.NAN;
1957                 } else {
1958                     // some intermediate numbers exceeded capacity,
1959                     // and the low order bits became NaN (because infinity - infinity = NaN)
1960                     if (FastMath.abs(full) < 1) {
1961                         return new Split(FastMath.copySign(0.0, full), 0.0);
1962                     } else if (full < 0 && (e & 0x1) == 1) {
1963                         return Split.NEGATIVE_INFINITY;
1964                     } else {
1965                         return Split.POSITIVE_INFINITY;
1966                     }
1967                 }
1968             } else {
1969                 return result;
1970             }
1971 
1972         }
1973 
1974     }
1975 
1976     /**
1977      *  Computes sin(x) - x, where |x| < 1/16.
1978      *  Use a Remez polynomial approximation.
1979      *  @param x a number smaller than 1/16
1980      *  @return sin(x) - x
1981      */
1982     private static double polySine(final double x)
1983     {
1984         double x2 = x*x;
1985 
1986         double p = 2.7553817452272217E-6;
1987         p = p * x2 + -1.9841269659586505E-4;
1988         p = p * x2 + 0.008333333333329196;
1989         p = p * x2 + -0.16666666666666666;
1990         //p *= x2;
1991         //p *= x;
1992         p = p * x2 * x;
1993 
1994         return p;
1995     }
1996 
1997     /**
1998      *  Computes cos(x) - 1, where |x| < 1/16.
1999      *  Use a Remez polynomial approximation.
2000      *  @param x a number smaller than 1/16
2001      *  @return cos(x) - 1
2002      */
2003     private static double polyCosine(double x) {
2004         double x2 = x*x;
2005 
2006         double p = 2.479773539153719E-5;
2007         p = p * x2 + -0.0013888888689039883;
2008         p = p * x2 + 0.041666666666621166;
2009         p = p * x2 + -0.49999999999999994;
2010         p *= x2;
2011 
2012         return p;
2013     }
2014 
2015     /**
2016      *  Compute sine over the first quadrant (0 < x < pi/2).
2017      *  Use combination of table lookup and rational polynomial expansion.
2018      *  @param xa number from which sine is requested
2019      *  @param xb extra bits for x (may be 0.0)
2020      *  @return sin(xa + xb)
2021      */
2022     private static double sinQ(double xa, double xb) {
2023         int idx = (int) ((xa * 8.0) + 0.5);
2024         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
2025 
2026         // Table lookups
2027         final double sintA = SINE_TABLE_A[idx];
2028         final double sintB = SINE_TABLE_B[idx];
2029         final double costA = COSINE_TABLE_A[idx];
2030         final double costB = COSINE_TABLE_B[idx];
2031 
2032         // Polynomial eval of sin(epsilon), cos(epsilon)
2033         double sinEpsA = epsilon;
2034         double sinEpsB = polySine(epsilon);
2035         final double cosEpsA = 1.0;
2036         final double cosEpsB = polyCosine(epsilon);
2037 
2038         // Split epsilon   xa + xb = x
2039         final double temp = sinEpsA * HEX_40000000;
2040         double temp2 = (sinEpsA + temp) - temp;
2041         sinEpsB +=  sinEpsA - temp2;
2042         sinEpsA = temp2;
2043 
2044         /* Compute sin(x) by angle addition formula */
2045         double result;
2046 
2047         /* Compute the following sum:
2048          *
2049          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2050          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2051          *
2052          * Ranges of elements
2053          *
2054          * xxxtA   0            PI/2
2055          * xxxtB   -1.5e-9      1.5e-9
2056          * sinEpsA -0.0625      0.0625
2057          * sinEpsB -6e-11       6e-11
2058          * cosEpsA  1.0
2059          * cosEpsB  0           -0.0625
2060          *
2061          */
2062 
2063         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2064         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2065 
2066         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
2067         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
2068         double a = 0;
2069         double b = 0;
2070 
2071         double t = sintA;
2072         double c = a + t;
2073         double d = -(c - a - t);
2074         a = c;
2075         b += d;
2076 
2077         t = costA * sinEpsA;
2078         c = a + t;
2079         d = -(c - a - t);
2080         a = c;
2081         b += d;
2082 
2083         b = b + sintA * cosEpsB + costA * sinEpsB;
2084         /*
2085     t = sintA*cosEpsB;
2086     c = a + t;
2087     d = -(c - a - t);
2088     a = c;
2089     b = b + d;
2090 
2091     t = costA*sinEpsB;
2092     c = a + t;
2093     d = -(c - a - t);
2094     a = c;
2095     b = b + d;
2096          */
2097 
2098         b = b + sintB + costB * sinEpsA + sintB * cosEpsB + costB * sinEpsB;
2099         /*
2100     t = sintB;
2101     c = a + t;
2102     d = -(c - a - t);
2103     a = c;
2104     b = b + d;
2105 
2106     t = costB*sinEpsA;
2107     c = a + t;
2108     d = -(c - a - t);
2109     a = c;
2110     b = b + d;
2111 
2112     t = sintB*cosEpsB;
2113     c = a + t;
2114     d = -(c - a - t);
2115     a = c;
2116     b = b + d;
2117 
2118     t = costB*sinEpsB;
2119     c = a + t;
2120     d = -(c - a - t);
2121     a = c;
2122     b = b + d;
2123          */
2124 
2125         if (xb != 0.0) {
2126             t = ((costA + costB) * (cosEpsA + cosEpsB) -
2127                  (sintA + sintB) * (sinEpsA + sinEpsB)) * xb;  // approximate cosine*xb
2128             c = a + t;
2129             d = -(c - a - t);
2130             a = c;
2131             b += d;
2132         }
2133 
2134         result = a + b;
2135 
2136         return result;
2137     }
2138 
2139     /**
2140      * Compute cosine in the first quadrant by subtracting input from PI/2 and
2141      * then calling sinQ.  This is more accurate as the input approaches PI/2.
2142      *  @param xa number from which cosine is requested
2143      *  @param xb extra bits for x (may be 0.0)
2144      *  @return cos(xa + xb)
2145      */
2146     private static double cosQ(double xa, double xb) {
2147         final double pi2a = 1.5707963267948966;
2148         final double pi2b = 6.123233995736766E-17;
2149 
2150         final double a = pi2a - xa;
2151         double b = -(a - pi2a + xa);
2152         b += pi2b - xb;
2153 
2154         return sinQ(a, b);
2155     }
2156 
2157     /**
2158      *  Compute tangent (or cotangent) over the first quadrant.   0 < x < pi/2
2159      *  Use combination of table lookup and rational polynomial expansion.
2160      *  @param xa number from which sine is requested
2161      *  @param xb extra bits for x (may be 0.0)
2162      *  @param cotanFlag if true, compute the cotangent instead of the tangent
2163      *  @return tan(xa+xb) (or cotangent, depending on cotanFlag)
2164      */
2165     private static double tanQ(double xa, double xb, boolean cotanFlag) {
2166 
2167         int idx = (int) ((xa * 8.0) + 0.5);
2168         final double epsilon = xa - EIGHTHS[idx]; //idx*0.125;
2169 
2170         // Table lookups
2171         final double sintA = SINE_TABLE_A[idx];
2172         final double sintB = SINE_TABLE_B[idx];
2173         final double costA = COSINE_TABLE_A[idx];
2174         final double costB = COSINE_TABLE_B[idx];
2175 
2176         // Polynomial eval of sin(epsilon), cos(epsilon)
2177         double sinEpsA = epsilon;
2178         double sinEpsB = polySine(epsilon);
2179         final double cosEpsA = 1.0;
2180         final double cosEpsB = polyCosine(epsilon);
2181 
2182         // Split epsilon   xa + xb = x
2183         double temp = sinEpsA * HEX_40000000;
2184         double temp2 = (sinEpsA + temp) - temp;
2185         sinEpsB +=  sinEpsA - temp2;
2186         sinEpsA = temp2;
2187 
2188         /* Compute sin(x) by angle addition formula */
2189 
2190         /* Compute the following sum:
2191          *
2192          * result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2193          *          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2194          *
2195          * Ranges of elements
2196          *
2197          * xxxtA   0            PI/2
2198          * xxxtB   -1.5e-9      1.5e-9
2199          * sinEpsA -0.0625      0.0625
2200          * sinEpsB -6e-11       6e-11
2201          * cosEpsA  1.0
2202          * cosEpsB  0           -0.0625
2203          *
2204          */
2205 
2206         //result = sintA + costA*sinEpsA + sintA*cosEpsB + costA*sinEpsB +
2207         //          sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2208 
2209         //result = sintA + sintA*cosEpsB + sintB + sintB * cosEpsB;
2210         //result += costA*sinEpsA + costA*sinEpsB + costB*sinEpsA + costB * sinEpsB;
2211         double a = 0;
2212         double b = 0;
2213 
2214         // Compute sine
2215         double t = sintA;
2216         double c = a + t;
2217         double d = -(c - a - t);
2218         a = c;
2219         b += d;
2220 
2221         t = costA*sinEpsA;
2222         c = a + t;
2223         d = -(c - a - t);
2224         a = c;
2225         b += d;
2226 
2227         b += sintA*cosEpsB + costA*sinEpsB;
2228         b += sintB + costB*sinEpsA + sintB*cosEpsB + costB*sinEpsB;
2229 
2230         double sina = a + b;
2231         double sinb = -(sina - a - b);
2232 
2233         // Compute cosine
2234 
2235         a = 0.0;
2236         b = 0.0;
2237 
2238         t = costA*cosEpsA;
2239         c = a + t;
2240         d = -(c - a - t);
2241         a = c;
2242         b += d;
2243 
2244         t = -sintA*sinEpsA;
2245         c = a + t;
2246         d = -(c - a - t);
2247         a = c;
2248         b += d;
2249 
2250         b += costB*cosEpsA + costA*cosEpsB + costB*cosEpsB;
2251         b -= sintB*sinEpsA + sintA*sinEpsB + sintB*sinEpsB;
2252 
2253         double cosa = a + b;
2254         double cosb = -(cosa - a - b);
2255 
2256         if (cotanFlag) {
2257             double tmp;
2258             tmp = cosa; cosa = sina; sina = tmp;
2259             tmp = cosb; cosb = sinb; sinb = tmp;
2260         }
2261 
2262 
2263         /* estimate and correct, compute 1.0/(cosa+cosb) */
2264         /*
2265     double est = (sina+sinb)/(cosa+cosb);
2266     double err = (sina - cosa*est) + (sinb - cosb*est);
2267     est += err/(cosa+cosb);
2268     err = (sina - cosa*est) + (sinb - cosb*est);
2269          */
2270 
2271         // f(x) = 1/x,   f'(x) = -1/x^2
2272 
2273         double est = sina/cosa;
2274 
2275         /* Split the estimate to get more accurate read on division rounding */
2276         temp = est * HEX_40000000;
2277         double esta = (est + temp) - temp;
2278         double estb =  est - esta;
2279 
2280         temp = cosa * HEX_40000000;
2281         double cosaa = (cosa + temp) - temp;
2282         double cosab =  cosa - cosaa;
2283 
2284         //double err = (sina - est*cosa)/cosa;  // Correction for division rounding
2285         double err = (sina - esta*cosaa - esta*cosab - estb*cosaa - estb*cosab)/cosa;  // Correction for division rounding
2286         err += sinb/cosa;                     // Change in est due to sinb
2287         err += -sina * cosb / cosa / cosa;    // Change in est due to cosb
2288 
2289         if (xb != 0.0) {
2290             // tan' = 1 + tan^2      cot' = -(1 + cot^2)
2291             // Approximate impact of xb
2292             double xbadj = xb + est*est*xb;
2293             if (cotanFlag) {
2294                 xbadj = -xbadj;
2295             }
2296 
2297             err += xbadj;
2298         }
2299 
2300         return est+err;
2301     }
2302 
2303     /** Reduce the input argument using the Payne and Hanek method.
2304      *  This is good for all inputs 0.0 < x < inf
2305      *  Output is remainder after dividing by PI/2
2306      *  The result array should contain 3 numbers.
2307      *  result[0] is the integer portion, so mod 4 this gives the quadrant.
2308      *  result[1] is the upper bits of the remainder
2309      *  result[2] is the lower bits of the remainder
2310      *
2311      * @param x number to reduce
2312      * @param result placeholder where to put the result
2313      */
2314     private static void reducePayneHanek(double x, double result[])
2315     {
2316         /* Convert input double to bits */
2317         long inbits = Double.doubleToRawLongBits(x);
2318         int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
2319 
2320         /* Convert to fixed point representation */
2321         inbits &= 0x000fffffffffffffL;
2322         inbits |= 0x0010000000000000L;
2323 
2324         /* Normalize input to be between 0.5 and 1.0 */
2325         exponent++;
2326         inbits <<= 11;
2327 
2328         /* Based on the exponent, get a shifted copy of recip2pi */
2329         long shpi0;
2330         long shpiA;
2331         long shpiB;
2332         int idx = exponent >> 6;
2333         int shift = exponent - (idx << 6);
2334 
2335         if (shift != 0) {
2336             shpi0 = (idx == 0) ? 0 : (RECIP_2PI[idx-1] << shift);
2337             shpi0 |= RECIP_2PI[idx] >>> (64-shift);
2338             shpiA = (RECIP_2PI[idx] << shift) | (RECIP_2PI[idx+1] >>> (64-shift));
2339             shpiB = (RECIP_2PI[idx+1] << shift) | (RECIP_2PI[idx+2] >>> (64-shift));
2340         } else {
2341             shpi0 = (idx == 0) ? 0 : RECIP_2PI[idx-1];
2342             shpiA = RECIP_2PI[idx];
2343             shpiB = RECIP_2PI[idx+1];
2344         }
2345 
2346         /* Multiply input by shpiA */
2347         long a = inbits >>> 32;
2348         long b = inbits & 0xffffffffL;
2349 
2350         long c = shpiA >>> 32;
2351         long d = shpiA & 0xffffffffL;
2352 
2353         long ac = a * c;
2354         long bd = b * d;
2355         long bc = b * c;
2356         long ad = a * d;
2357 
2358         long prodB = bd + (ad << 32);
2359         long prodA = ac + (ad >>> 32);
2360 
2361         boolean bita = (bd & 0x8000000000000000L) != 0;
2362         boolean bitb = (ad & 0x80000000L ) != 0;
2363         boolean bitsum = (prodB & 0x8000000000000000L) != 0;
2364 
2365         /* Carry */
2366         if ( (bita && bitb) ||
2367                 ((bita || bitb) && !bitsum) ) {
2368             prodA++;
2369         }
2370 
2371         bita = (prodB & 0x8000000000000000L) != 0;
2372         bitb = (bc & 0x80000000L ) != 0;
2373 
2374         prodB += bc << 32;
2375         prodA += bc >>> 32;
2376 
2377         bitsum = (prodB & 0x8000000000000000L) != 0;
2378 
2379         /* Carry */
2380         if ( (bita && bitb) ||
2381                 ((bita || bitb) && !bitsum) ) {
2382             prodA++;
2383         }
2384 
2385         /* Multiply input by shpiB */
2386         c = shpiB >>> 32;
2387         d = shpiB & 0xffffffffL;
2388         ac = a * c;
2389         bc = b * c;
2390         ad = a * d;
2391 
2392         /* Collect terms */
2393         ac += (bc + ad) >>> 32;
2394 
2395         bita = (prodB & 0x8000000000000000L) != 0;
2396         bitb = (ac & 0x8000000000000000L ) != 0;
2397         prodB += ac;
2398         bitsum = (prodB & 0x8000000000000000L) != 0;
2399         /* Carry */
2400         if ( (bita && bitb) ||
2401                 ((bita || bitb) && !bitsum) ) {
2402             prodA++;
2403         }
2404 
2405         /* Multiply by shpi0 */
2406         c = shpi0 >>> 32;
2407         d = shpi0 & 0xffffffffL;
2408 
2409         bd = b * d;
2410         bc = b * c;
2411         ad = a * d;
2412 
2413         prodA += bd + ((bc + ad) << 32);
2414 
2415         /*
2416          * prodA, prodB now contain the remainder as a fraction of PI.  We want this as a fraction of
2417          * PI/2, so use the following steps:
2418          * 1.) multiply by 4.
2419          * 2.) do a fixed point muliply by PI/4.
2420          * 3.) Convert to floating point.
2421          * 4.) Multiply by 2
2422          */
2423 
2424         /* This identifies the quadrant */
2425         int intPart = (int)(prodA >>> 62);
2426 
2427         /* Multiply by 4 */
2428         prodA <<= 2;
2429         prodA |= prodB >>> 62;
2430         prodB <<= 2;
2431 
2432         /* Multiply by PI/4 */
2433         a = prodA >>> 32;
2434         b = prodA & 0xffffffffL;
2435 
2436         c = PI_O_4_BITS[0] >>> 32;
2437         d = PI_O_4_BITS[0] & 0xffffffffL;
2438 
2439         ac = a * c;
2440         bd = b * d;
2441         bc = b * c;
2442         ad = a * d;
2443 
2444         long prod2B = bd + (ad << 32);
2445         long prod2A = ac + (ad >>> 32);
2446 
2447         bita = (bd & 0x8000000000000000L) != 0;
2448         bitb = (ad & 0x80000000L ) != 0;
2449         bitsum = (prod2B & 0x8000000000000000L) != 0;
2450 
2451         /* Carry */
2452         if ( (bita && bitb) ||
2453                 ((bita || bitb) && !bitsum) ) {
2454             prod2A++;
2455         }
2456 
2457         bita = (prod2B & 0x8000000000000000L) != 0;
2458         bitb = (bc & 0x80000000L ) != 0;
2459 
2460         prod2B += bc << 32;
2461         prod2A += bc >>> 32;
2462 
2463         bitsum = (prod2B & 0x8000000000000000L) != 0;
2464 
2465         /* Carry */
2466         if ( (bita && bitb) ||
2467                 ((bita || bitb) && !bitsum) ) {
2468             prod2A++;
2469         }
2470 
2471         /* Multiply input by pio4bits[1] */
2472         c = PI_O_4_BITS[1] >>> 32;
2473         d = PI_O_4_BITS[1] & 0xffffffffL;
2474         ac = a * c;
2475         bc = b * c;
2476         ad = a * d;
2477 
2478         /* Collect terms */
2479         ac += (bc + ad) >>> 32;
2480 
2481         bita = (prod2B & 0x8000000000000000L) != 0;
2482         bitb = (ac & 0x8000000000000000L ) != 0;
2483         prod2B += ac;
2484         bitsum = (prod2B & 0x8000000000000000L) != 0;
2485         /* Carry */
2486         if ( (bita && bitb) ||
2487                 ((bita || bitb) && !bitsum) ) {
2488             prod2A++;
2489         }
2490 
2491         /* Multiply inputB by pio4bits[0] */
2492         a = prodB >>> 32;
2493         b = prodB & 0xffffffffL;
2494         c = PI_O_4_BITS[0] >>> 32;
2495         d = PI_O_4_BITS[0] & 0xffffffffL;
2496         ac = a * c;
2497         bc = b * c;
2498         ad = a * d;
2499 
2500         /* Collect terms */
2501         ac += (bc + ad) >>> 32;
2502 
2503         bita = (prod2B & 0x8000000000000000L) != 0;
2504         bitb = (ac & 0x8000000000000000L ) != 0;
2505         prod2B += ac;
2506         bitsum = (prod2B & 0x8000000000000000L) != 0;
2507         /* Carry */
2508         if ( (bita && bitb) ||
2509                 ((bita || bitb) && !bitsum) ) {
2510             prod2A++;
2511         }
2512 
2513         /* Convert to double */
2514         double tmpA = (prod2A >>> 12) / TWO_POWER_52;  // High order 52 bits
2515         double tmpB = (((prod2A & 0xfffL) << 40) + (prod2B >>> 24)) / TWO_POWER_52 / TWO_POWER_52; // Low bits
2516 
2517         double sumA = tmpA + tmpB;
2518         double sumB = -(sumA - tmpA - tmpB);
2519 
2520         /* Multiply by PI/2 and return */
2521         result[0] = intPart;
2522         result[1] = sumA * 2.0;
2523         result[2] = sumB * 2.0;
2524     }
2525 
2526     /**
2527      * Sine function.
2528      *
2529      * @param x Argument.
2530      * @return sin(x)
2531      */
2532     public static double sin(double x) {
2533         boolean negative = false;
2534         int quadrant = 0;
2535         double xa;
2536         double xb = 0.0;
2537 
2538         /* Take absolute value of the input */
2539         xa = x;
2540         if (x < 0) {
2541             negative = true;
2542             xa = -xa;
2543         }
2544 
2545         /* Check for zero and negative zero */
2546         if (xa == 0.0) {
2547             long bits = Double.doubleToRawLongBits(x);
2548             if (bits < 0) {
2549                 return -0.0;
2550             }
2551             return 0.0;
2552         }
2553 
2554         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2555             return Double.NaN;
2556         }
2557 
2558         /* Perform any argument reduction */
2559         if (xa > 3294198.0) {
2560             // PI * (2**20)
2561             // Argument too big for CodyWaite reduction.  Must use
2562             // PayneHanek.
2563             double reduceResults[] = new double[3];
2564             reducePayneHanek(xa, reduceResults);
2565             quadrant = ((int) reduceResults[0]) & 3;
2566             xa = reduceResults[1];
2567             xb = reduceResults[2];
2568         } else if (xa > 1.5707963267948966) {
2569             final CodyWaite cw = new CodyWaite(xa);
2570             quadrant = cw.getK() & 3;
2571             xa = cw.getRemA();
2572             xb = cw.getRemB();
2573         }
2574 
2575         if (negative) {
2576             quadrant ^= 2;  // Flip bit 1
2577         }
2578 
2579         switch (quadrant) {
2580             case 0:
2581                 return sinQ(xa, xb);
2582             case 1:
2583                 return cosQ(xa, xb);
2584             case 2:
2585                 return -sinQ(xa, xb);
2586             case 3:
2587                 return -cosQ(xa, xb);
2588             default:
2589                 return Double.NaN;
2590         }
2591     }
2592 
2593     /**
2594      * Cosine function.
2595      *
2596      * @param x Argument.
2597      * @return cos(x)
2598      */
2599     public static double cos(double x) {
2600         int quadrant = 0;
2601 
2602         /* Take absolute value of the input */
2603         double xa = x;
2604         if (x < 0) {
2605             xa = -xa;
2606         }
2607 
2608         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2609             return Double.NaN;
2610         }
2611 
2612         /* Perform any argument reduction */
2613         double xb = 0;
2614         if (xa > 3294198.0) {
2615             // PI * (2**20)
2616             // Argument too big for CodyWaite reduction.  Must use
2617             // PayneHanek.
2618             double reduceResults[] = new double[3];
2619             reducePayneHanek(xa, reduceResults);
2620             quadrant = ((int) reduceResults[0]) & 3;
2621             xa = reduceResults[1];
2622             xb = reduceResults[2];
2623         } else if (xa > 1.5707963267948966) {
2624             final CodyWaite cw = new CodyWaite(xa);
2625             quadrant = cw.getK() & 3;
2626             xa = cw.getRemA();
2627             xb = cw.getRemB();
2628         }
2629 
2630         //if (negative)
2631         //  quadrant = (quadrant + 2) % 4;
2632 
2633         switch (quadrant) {
2634             case 0:
2635                 return cosQ(xa, xb);
2636             case 1:
2637                 return -sinQ(xa, xb);
2638             case 2:
2639                 return -cosQ(xa, xb);
2640             case 3:
2641                 return sinQ(xa, xb);
2642             default:
2643                 return Double.NaN;
2644         }
2645     }
2646 
2647     /**
2648      * Combined Sine and Cosine function.
2649      *
2650      * @param x Argument.
2651      * @return [sin(x), cos(x)]
2652      */
2653     public static SinCos sinCos(double x) {
2654         boolean negative = false;
2655         int quadrant = 0;
2656         double xa;
2657         double xb = 0.0;
2658 
2659         /* Take absolute value of the input */
2660         xa = x;
2661         if (x < 0) {
2662             negative = true;
2663             xa = -xa;
2664         }
2665 
2666         /* Check for zero and negative zero */
2667         if (xa == 0.0) {
2668             long bits = Double.doubleToRawLongBits(x);
2669             if (bits < 0) {
2670                 return new SinCos(-0.0, 1.0);
2671             }
2672             return new SinCos(0.0, 1.0);
2673         }
2674 
2675         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2676             return new SinCos(Double.NaN, Double.NaN);
2677         }
2678 
2679         /* Perform any argument reduction */
2680         if (xa > 3294198.0) {
2681             // PI * (2**20)
2682             // Argument too big for CodyWaite reduction.  Must use
2683             // PayneHanek.
2684             double reduceResults[] = new double[3];
2685             reducePayneHanek(xa, reduceResults);
2686             quadrant = ((int) reduceResults[0]) & 3;
2687             xa = reduceResults[1];
2688             xb = reduceResults[2];
2689         } else if (xa > 1.5707963267948966) {
2690             final CodyWaite cw = new CodyWaite(xa);
2691             quadrant = cw.getK() & 3;
2692             xa = cw.getRemA();
2693             xb = cw.getRemB();
2694         }
2695 
2696         switch (quadrant) {
2697             case 0:
2698                 return new SinCos(negative ? -sinQ(xa, xb) :  sinQ(xa, xb),  cosQ(xa, xb));
2699             case 1:
2700                 return new SinCos(negative ? -cosQ(xa, xb) :  cosQ(xa, xb), -sinQ(xa, xb));
2701             case 2:
2702                 return new SinCos(negative ?  sinQ(xa, xb) : -sinQ(xa, xb), -cosQ(xa, xb));
2703             case 3:
2704                 return new SinCos(negative ?  cosQ(xa, xb) : -cosQ(xa, xb),  sinQ(xa, xb));
2705             default:
2706                 return new SinCos(Double.NaN, Double.NaN);
2707         }
2708     }
2709 
2710     /**
2711      * Combined Sine and Cosine function.
2712      *
2713      * @param x Argument.
2714      * @param <T> the type of the field element
2715      * @return [sin(x), cos(x)]
2716      * @since 1.4
2717      */
2718     public static <T extends CalculusFieldElement<T>> FieldSinCos<T> sinCos(T x) {
2719         return x.sinCos();
2720     }
2721 
2722     /**
2723      * Tangent function.
2724      *
2725      * @param x Argument.
2726      * @return tan(x)
2727      */
2728     public static double tan(double x) {
2729         boolean negative = false;
2730         int quadrant = 0;
2731 
2732         /* Take absolute value of the input */
2733         double xa = x;
2734         if (x < 0) {
2735             negative = true;
2736             xa = -xa;
2737         }
2738 
2739         /* Check for zero and negative zero */
2740         if (xa == 0.0) {
2741             long bits = Double.doubleToRawLongBits(x);
2742             if (bits < 0) {
2743                 return -0.0;
2744             }
2745             return 0.0;
2746         }
2747 
2748         if (xa != xa || xa == Double.POSITIVE_INFINITY) {
2749             return Double.NaN;
2750         }
2751 
2752         /* Perform any argument reduction */
2753         double xb = 0;
2754         if (xa > 3294198.0) {
2755             // PI * (2**20)
2756             // Argument too big for CodyWaite reduction.  Must use
2757             // PayneHanek.
2758             double reduceResults[] = new double[3];
2759             reducePayneHanek(xa, reduceResults);
2760             quadrant = ((int) reduceResults[0]) & 3;
2761             xa = reduceResults[1];
2762             xb = reduceResults[2];
2763         } else if (xa > 1.5707963267948966) {
2764             final CodyWaite cw = new CodyWaite(xa);
2765             quadrant = cw.getK() & 3;
2766             xa = cw.getRemA();
2767             xb = cw.getRemB();
2768         }
2769 
2770         if (xa > 1.5) {
2771             // Accuracy suffers between 1.5 and PI/2
2772             final double pi2a = 1.5707963267948966;
2773             final double pi2b = 6.123233995736766E-17;
2774 
2775             final double a = pi2a - xa;
2776             double b = -(a - pi2a + xa);
2777             b += pi2b - xb;
2778 
2779             xa = a + b;
2780             xb = -(xa - a - b);
2781             quadrant ^= 1;
2782             negative ^= true;
2783         }
2784 
2785         double result;
2786         if ((quadrant & 1) == 0) {
2787             result = tanQ(xa, xb, false);
2788         } else {
2789             result = -tanQ(xa, xb, true);
2790         }
2791 
2792         if (negative) {
2793             result = -result;
2794         }
2795 
2796         return result;
2797     }
2798 
2799     /**
2800      * Arctangent function
2801      *  @param x a number
2802      *  @return atan(x)
2803      */
2804     public static double atan(double x) {
2805         return atan(x, 0.0, false);
2806     }
2807 
2808     /** Internal helper function to compute arctangent.
2809      * @param xa number from which arctangent is requested
2810      * @param xb extra bits for x (may be 0.0)
2811      * @param leftPlane if true, result angle must be put in the left half plane
2812      * @return atan(xa + xb) (or angle shifted by {@code PI} if leftPlane is true)
2813      */
2814     private static double atan(double xa, double xb, boolean leftPlane) {
2815         if (xa == 0.0) { // Matches +/- 0.0; return correct sign
2816             return leftPlane ? copySign(Math.PI, xa) : xa;
2817         }
2818 
2819         final boolean negate;
2820         if (xa < 0) {
2821             // negative
2822             xa = -xa;
2823             xb = -xb;
2824             negate = true;
2825         } else {
2826             negate = false;
2827         }
2828 
2829         if (xa > 1.633123935319537E16) { // Very large input
2830             return (negate ^ leftPlane) ? (-Math.PI * F_1_2) : (Math.PI * F_1_2);
2831         }
2832 
2833         /* Estimate the closest tabulated arctan value, compute eps = xa-tangentTable */
2834         final int idx;
2835         if (xa < 1) {
2836             idx = (int) (((-1.7168146928204136 * xa * xa + 8.0) * xa) + 0.5);
2837         } else {
2838             final double oneOverXa = 1 / xa;
2839             idx = (int) (-((-1.7168146928204136 * oneOverXa * oneOverXa + 8.0) * oneOverXa) + 13.07);
2840         }
2841 
2842         final double ttA = TANGENT_TABLE_A[idx];
2843         final double ttB = TANGENT_TABLE_B[idx];
2844 
2845         double epsA = xa - ttA;
2846         double epsB = -(epsA - xa + ttA);
2847         epsB += xb - ttB;
2848 
2849         double temp = epsA + epsB;
2850         epsB = -(temp - epsA - epsB);
2851         epsA = temp;
2852 
2853         /* Compute eps = eps / (1.0 + xa*tangent) */
2854         temp = xa * HEX_40000000;
2855         double ya = xa + temp - temp;
2856         double yb = xb + xa - ya;
2857         xa = ya;
2858         xb += yb;
2859 
2860         //if (idx > 8 || idx == 0)
2861         if (idx == 0) {
2862             /* If the slope of the arctan is gentle enough (< 0.45), this approximation will suffice */
2863             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx] + xb*tangentTableA[idx] + xa*tangentTableB[idx] + xb*tangentTableB[idx]);
2864             final double denom = 1d / (1d + (xa + xb) * (ttA + ttB));
2865             //double denom = 1.0 / (1.0 + xa*tangentTableA[idx]);
2866             ya = epsA * denom;
2867             yb = epsB * denom;
2868         } else {
2869             double temp2 = xa * ttA;
2870             double za = 1d + temp2;
2871             double zb = -(za - 1d - temp2);
2872             temp2 = xb * ttA + xa * ttB;
2873             temp = za + temp2;
2874             zb += -(temp - za - temp2);
2875             za = temp;
2876 
2877             zb += xb * ttB;
2878             ya = epsA / za;
2879 
2880             temp = ya * HEX_40000000;
2881             final double yaa = (ya + temp) - temp;
2882             final double yab = ya - yaa;
2883 
2884             temp = za * HEX_40000000;
2885             final double zaa = (za + temp) - temp;
2886             final double zab = za - zaa;
2887 
2888             /* Correct for rounding in division */
2889             yb = (epsA - yaa * zaa - yaa * zab - yab * zaa - yab * zab) / za;
2890 
2891             yb += -epsA * zb / za / za;
2892             yb += epsB / za;
2893         }
2894 
2895 
2896         epsA = ya;
2897         epsB = yb;
2898 
2899         /* Evaluate polynomial */
2900         final double epsA2 = epsA * epsA;
2901 
2902         /*
2903     yb = -0.09001346640161823;
2904     yb = yb * epsA2 + 0.11110718400605211;
2905     yb = yb * epsA2 + -0.1428571349122913;
2906     yb = yb * epsA2 + 0.19999999999273194;
2907     yb = yb * epsA2 + -0.33333333333333093;
2908     yb = yb * epsA2 * epsA;
2909          */
2910 
2911         yb = 0.07490822288864472;
2912         yb = yb * epsA2 - 0.09088450866185192;
2913         yb = yb * epsA2 + 0.11111095942313305;
2914         yb = yb * epsA2 - 0.1428571423679182;
2915         yb = yb * epsA2 + 0.19999999999923582;
2916         yb = yb * epsA2 - 0.33333333333333287;
2917         yb = yb * epsA2 * epsA;
2918 
2919 
2920         ya = epsA;
2921 
2922         temp = ya + yb;
2923         yb = -(temp - ya - yb);
2924         ya = temp;
2925 
2926         /* Add in effect of epsB.   atan'(x) = 1/(1+x^2) */
2927         yb += epsB / (1d + epsA * epsA);
2928 
2929         final double eighths = EIGHTHS[idx];
2930 
2931         //result = yb + eighths[idx] + ya;
2932         double za = eighths + ya;
2933         double zb = -(za - eighths - ya);
2934         temp = za + yb;
2935         zb += -(temp - za - yb);
2936         za = temp;
2937 
2938         double result = za + zb;
2939 
2940         if (leftPlane) {
2941             // Result is in the left plane
2942             final double resultb = -(result - za - zb);
2943             final double pia = 1.5707963267948966 * 2;
2944             final double pib = 6.123233995736766E-17 * 2;
2945 
2946             za = pia - result;
2947             zb = -(za - pia + result);
2948             zb += pib - resultb;
2949 
2950             result = za + zb;
2951         }
2952 
2953 
2954         if (negate ^ leftPlane) {
2955             result = -result;
2956         }
2957 
2958         return result;
2959     }
2960 
2961     /**
2962      * Two arguments arctangent function
2963      * @param y ordinate
2964      * @param x abscissa
2965      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
2966      */
2967     public static double atan2(double y, double x) {
2968         if (Double.isNaN(x) || Double.isNaN(y)) {
2969             return Double.NaN;
2970         }
2971 
2972         if (y == 0) {
2973             final double result = x * y;
2974             final double invx = 1d / x;
2975 
2976             if (invx == 0) { // X is infinite
2977                 if (x > 0) {
2978                     return y; // return +/- 0.0
2979                 } else {
2980                     return copySign(Math.PI, y);
2981                 }
2982             }
2983 
2984             if (x < 0 || invx < 0) {
2985                 return copySign(Math.PI, y);
2986             } else {
2987                 return result;
2988             }
2989         }
2990 
2991         // y cannot now be zero
2992 
2993         if (y == Double.POSITIVE_INFINITY) {
2994             if (x == Double.POSITIVE_INFINITY) {
2995                 return Math.PI * F_1_4;
2996             }
2997 
2998             if (x == Double.NEGATIVE_INFINITY) {
2999                 return Math.PI * F_3_4;
3000             }
3001 
3002             return Math.PI * F_1_2;
3003         }
3004 
3005         if (y == Double.NEGATIVE_INFINITY) {
3006             if (x == Double.POSITIVE_INFINITY) {
3007                 return -Math.PI * F_1_4;
3008             }
3009 
3010             if (x == Double.NEGATIVE_INFINITY) {
3011                 return -Math.PI * F_3_4;
3012             }
3013 
3014             return -Math.PI * F_1_2;
3015         }
3016 
3017         if (x == Double.POSITIVE_INFINITY) {
3018             return copySign(0d, y);
3019         }
3020 
3021         if (x == Double.NEGATIVE_INFINITY)
3022         {
3023             return copySign(Math.PI, y);
3024         }
3025 
3026         // Neither y nor x can be infinite or NAN here
3027 
3028         if (x == 0) {
3029             return copySign(Math.PI * F_1_2, y);
3030         }
3031 
3032         // Compute ratio r = y/x
3033         final double r = y / x;
3034         if (Double.isInfinite(r)) { // bypass calculations that can create NaN
3035             return atan(r, 0, x < 0);
3036         }
3037 
3038         double ra = doubleHighPart(r);
3039         double rb = r - ra;
3040 
3041         // Split x
3042         final double xa = doubleHighPart(x);
3043         final double xb = x - xa;
3044 
3045         rb += (y - ra * xa - ra * xb - rb * xa - rb * xb) / x;
3046 
3047         final double temp = ra + rb;
3048         rb = -(temp - ra - rb);
3049         ra = temp;
3050 
3051         if (ra == 0) { // Fix up the sign so atan works correctly
3052             ra = copySign(0d, y);
3053         }
3054 
3055         // Call atan
3056         return atan(ra, rb, x < 0);
3057 
3058     }
3059 
3060     /** Compute the arc sine of a number.
3061      * @param x number on which evaluation is done
3062      * @return arc sine of x
3063      */
3064     public static double asin(double x) {
3065       if (Double.isNaN(x)) {
3066           return Double.NaN;
3067       }
3068 
3069       if (x > 1.0 || x < -1.0) {
3070           return Double.NaN;
3071       }
3072 
3073       if (x == 1.0) {
3074           return Math.PI/2.0;
3075       }
3076 
3077       if (x == -1.0) {
3078           return -Math.PI/2.0;
3079       }
3080 
3081       if (x == 0.0) { // Matches +/- 0.0; return correct sign
3082           return x;
3083       }
3084 
3085       /* Compute asin(x) = atan(x/sqrt(1-x*x)) */
3086 
3087       /* Split x */
3088       double temp = x * HEX_40000000;
3089       final double xa = x + temp - temp;
3090       final double xb = x - xa;
3091 
3092       /* Square it */
3093       double ya = xa*xa;
3094       double yb = xa*xb*2.0 + xb*xb;
3095 
3096       /* Subtract from 1 */
3097       ya = -ya;
3098       yb = -yb;
3099 
3100       double za = 1.0 + ya;
3101       double zb = -(za - 1.0 - ya);
3102 
3103       temp = za + yb;
3104       zb += -(temp - za - yb);
3105       za = temp;
3106 
3107       /* Square root */
3108       double y;
3109       y = sqrt(za);
3110       temp = y * HEX_40000000;
3111       ya = y + temp - temp;
3112       yb = y - ya;
3113 
3114       /* Extend precision of sqrt */
3115       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
3116 
3117       /* Contribution of zb to sqrt */
3118       double dx = zb / (2.0*y);
3119 
3120       // Compute ratio r = x/y
3121       double r = x/y;
3122       temp = r * HEX_40000000;
3123       double ra = r + temp - temp;
3124       double rb = r - ra;
3125 
3126       rb += (x - ra*ya - ra*yb - rb*ya - rb*yb) / y;  // Correct for rounding in division
3127       rb += -x * dx / y / y;  // Add in effect additional bits of sqrt.
3128 
3129       temp = ra + rb;
3130       rb = -(temp - ra - rb);
3131       ra = temp;
3132 
3133       return atan(ra, rb, false);
3134     }
3135 
3136     /** Compute the arc cosine of a number.
3137      * @param x number on which evaluation is done
3138      * @return arc cosine of x
3139      */
3140     public static double acos(double x) {
3141       if (Double.isNaN(x)) {
3142           return Double.NaN;
3143       }
3144 
3145       if (x > 1.0 || x < -1.0) {
3146           return Double.NaN;
3147       }
3148 
3149       if (x == -1.0) {
3150           return Math.PI;
3151       }
3152 
3153       if (x == 1.0) {
3154           return 0.0;
3155       }
3156 
3157       if (x == 0) {
3158           return Math.PI/2.0;
3159       }
3160 
3161       /* Compute acos(x) = atan(sqrt(1-x*x)/x) */
3162 
3163       /* Split x */
3164       double temp = x * HEX_40000000;
3165       final double xa = x + temp - temp;
3166       final double xb = x - xa;
3167 
3168       /* Square it */
3169       double ya = xa*xa;
3170       double yb = xa*xb*2.0 + xb*xb;
3171 
3172       /* Subtract from 1 */
3173       ya = -ya;
3174       yb = -yb;
3175 
3176       double za = 1.0 + ya;
3177       double zb = -(za - 1.0 - ya);
3178 
3179       temp = za + yb;
3180       zb += -(temp - za - yb);
3181       za = temp;
3182 
3183       /* Square root */
3184       double y = sqrt(za);
3185       temp = y * HEX_40000000;
3186       ya = y + temp - temp;
3187       yb = y - ya;
3188 
3189       /* Extend precision of sqrt */
3190       yb += (za - ya*ya - 2*ya*yb - yb*yb) / (2.0*y);
3191 
3192       /* Contribution of zb to sqrt */
3193       yb += zb / (2.0*y);
3194       y = ya+yb;
3195       yb = -(y - ya - yb);
3196 
3197       // Compute ratio r = y/x
3198       double r = y/x;
3199 
3200       // Did r overflow?
3201       if (Double.isInfinite(r)) { // x is effectively zero
3202           return Math.PI/2; // so return the appropriate value
3203       }
3204 
3205       double ra = doubleHighPart(r);
3206       double rb = r - ra;
3207 
3208       rb += (y - ra*xa - ra*xb - rb*xa - rb*xb) / x;  // Correct for rounding in division
3209       rb += yb / x;  // Add in effect additional bits of sqrt.
3210 
3211       temp = ra + rb;
3212       rb = -(temp - ra - rb);
3213       ra = temp;
3214 
3215       return atan(ra, rb, x<0);
3216     }
3217 
3218     /** Compute the cubic root of a number.
3219      * @param x number on which evaluation is done
3220      * @return cubic root of x
3221      */
3222     public static double cbrt(double x) {
3223       /* Convert input double to bits */
3224       long inbits = Double.doubleToRawLongBits(x);
3225       int exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
3226       boolean subnormal = false;
3227 
3228       if (exponent == -1023) {
3229           if (x == 0) {
3230               return x;
3231           }
3232 
3233           /* Subnormal, so normalize */
3234           subnormal = true;
3235           x *= 1.8014398509481984E16;  // 2^54
3236           inbits = Double.doubleToRawLongBits(x);
3237           exponent = (int) ((inbits >> 52) & 0x7ff) - 1023;
3238       }
3239 
3240       if (exponent == 1024) {
3241           // Nan or infinity.  Don't care which.
3242           return x;
3243       }
3244 
3245       /* Divide the exponent by 3 */
3246       int exp3 = exponent / 3;
3247 
3248       /* p2 will be the nearest power of 2 to x with its exponent divided by 3 */
3249       double p2 = Double.longBitsToDouble((inbits & 0x8000000000000000L) |
3250                                           (long)(((exp3 + 1023) & 0x7ff)) << 52);
3251 
3252       /* This will be a number between 1 and 2 */
3253       final double mant = Double.longBitsToDouble((inbits & 0x000fffffffffffffL) | 0x3ff0000000000000L);
3254 
3255       /* Estimate the cube root of mant by polynomial */
3256       double est = -0.010714690733195933;
3257       est = est * mant + 0.0875862700108075;
3258       est = est * mant + -0.3058015757857271;
3259       est = est * mant + 0.7249995199969751;
3260       est = est * mant + 0.5039018405998233;
3261 
3262       est *= CBRTTWO[exponent % 3 + 2];
3263 
3264       // est should now be good to about 15 bits of precision.   Do 2 rounds of
3265       // Newton's method to get closer,  this should get us full double precision
3266       // Scale down x for the purpose of doing newtons method.  This avoids over/under flows.
3267       final double xs = x / (p2*p2*p2);
3268       est += (xs - est*est*est) / (3*est*est);
3269       est += (xs - est*est*est) / (3*est*est);
3270 
3271       // Do one round of Newton's method in extended precision to get the last bit right.
3272       double temp = est * HEX_40000000;
3273       double ya = est + temp - temp;
3274       double yb = est - ya;
3275 
3276       double za = ya * ya;
3277       double zb = ya * yb * 2.0 + yb * yb;
3278       temp = za * HEX_40000000;
3279       double temp2 = za + temp - temp;
3280       zb += za - temp2;
3281       za = temp2;
3282 
3283       zb = za * yb + ya * zb + zb * yb;
3284       za *= ya;
3285 
3286       double na = xs - za;
3287       double nb = -(na - xs + za);
3288       nb -= zb;
3289 
3290       est += (na+nb)/(3*est*est);
3291 
3292       /* Scale by a power of two, so this is exact. */
3293       est *= p2;
3294 
3295       if (subnormal) {
3296           est *= 3.814697265625E-6;  // 2^-18
3297       }
3298 
3299       return est;
3300     }
3301 
3302     /**
3303      *  Convert degrees to radians, with error of less than 0.5 ULP
3304      *  @param x angle in degrees
3305      *  @return x converted into radians
3306      */
3307     public static double toRadians(double x)
3308     {
3309         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
3310             return x;
3311         }
3312 
3313         // These are PI/180 split into high and low order bits
3314         final double facta = 0.01745329052209854;
3315         final double factb = 1.997844754509471E-9;
3316 
3317         double xa = doubleHighPart(x);
3318         double xb = x - xa;
3319 
3320         double result = xb * factb + xb * facta + xa * factb + xa * facta;
3321         if (result == 0) {
3322             result *= x; // ensure correct sign if calculation underflows
3323         }
3324         return result;
3325     }
3326 
3327     /**
3328      *  Convert radians to degrees, with error of less than 0.5 ULP
3329      *  @param x angle in radians
3330      *  @return x converted into degrees
3331      */
3332     public static double toDegrees(double x)
3333     {
3334         if (Double.isInfinite(x) || x == 0.0) { // Matches +/- 0.0; return correct sign
3335             return x;
3336         }
3337 
3338         // These are 180/PI split into high and low order bits
3339         final double facta = 57.2957763671875;
3340         final double factb = 3.145894820876798E-6;
3341 
3342         double xa = doubleHighPart(x);
3343         double xb = x - xa;
3344 
3345         return xb * factb + xb * facta + xa * factb + xa * facta;
3346     }
3347 
3348     /**
3349      * Absolute value.
3350      * @param x number from which absolute value is requested
3351      * @return abs(x)
3352      */
3353     public static int abs(final int x) {
3354         final int i = x >>> 31;
3355         return (x ^ (~i + 1)) + i;
3356     }
3357 
3358     /**
3359      * Absolute value.
3360      * @param x number from which absolute value is requested
3361      * @return abs(x)
3362      */
3363     public static long abs(final long x) {
3364         final long l = x >>> 63;
3365         // l is one if x negative zero else
3366         // ~l+1 is zero if x is positive, -1 if x is negative
3367         // x^(~l+1) is x is x is positive, ~x if x is negative
3368         // add around
3369         return (x ^ (~l + 1)) + l;
3370     }
3371 
3372     /**
3373      * Absolute value.
3374      * @param x number from which absolute value is requested
3375      * @return abs(x), or throws an exception for {@code Integer.MIN_VALUE}
3376      * @since 2.0
3377      */
3378     public static int absExact(final int x) {
3379         if (x == Integer.MIN_VALUE) {
3380             throw new ArithmeticException();
3381         }
3382         return abs(x);
3383     }
3384 
3385     /**
3386      * Absolute value.
3387      * @param x number from which absolute value is requested
3388      * @return abs(x), or throws an exception for {@code Long.MIN_VALUE}
3389      * @since 2.0
3390      */
3391     public static long absExact(final long x) {
3392         if (x == Long.MIN_VALUE) {
3393             throw new ArithmeticException();
3394         }
3395         return abs(x);
3396     }
3397 
3398     /**
3399      * Absolute value.
3400      * @param x number from which absolute value is requested
3401      * @return abs(x)
3402      * @since 2.0
3403      */
3404     public static float abs(final float x) {
3405         return Float.intBitsToFloat(MASK_NON_SIGN_INT & Float.floatToRawIntBits(x));
3406     }
3407 
3408     /**
3409      * Absolute value.
3410      * @param x number from which absolute value is requested
3411      * @return abs(x)
3412      */
3413     public static double abs(double x) {
3414         return Double.longBitsToDouble(MASK_NON_SIGN_LONG & Double.doubleToRawLongBits(x));
3415     }
3416 
3417     /**
3418      * Negates the argument.
3419      * @param x number from which opposite value is requested
3420      * @return -x, or throws an exception for {@code Integer.MIN_VALUE}
3421      * @since 2.0
3422      */
3423     public static int negateExact(final int x) {
3424         if (x == Integer.MIN_VALUE) {
3425             throw new ArithmeticException();
3426         }
3427         return -x;
3428     }
3429 
3430     /**
3431      * Negates the argument.
3432      * @param x number from which opposite value is requested
3433      * @return -x, or throws an exception for {@code Long.MIN_VALUE}
3434      * @since 2.0
3435      */
3436     public static long negateExact(final long x) {
3437         if (x == Long.MIN_VALUE) {
3438             throw new ArithmeticException();
3439         }
3440         return -x;
3441     }
3442 
3443     /**
3444      * Compute least significant bit (Unit in Last Position) for a number.
3445      * @param x number from which ulp is requested
3446      * @return ulp(x)
3447      */
3448     public static double ulp(double x) {
3449         if (Double.isInfinite(x)) {
3450             return Double.POSITIVE_INFINITY;
3451         }
3452         return abs(x - Double.longBitsToDouble(Double.doubleToRawLongBits(x) ^ 1));
3453     }
3454 
3455     /**
3456      * Compute least significant bit (Unit in Last Position) for a number.
3457      * @param x number from which ulp is requested
3458      * @return ulp(x)
3459      */
3460     public static float ulp(float x) {
3461         if (Float.isInfinite(x)) {
3462             return Float.POSITIVE_INFINITY;
3463         }
3464         return abs(x - Float.intBitsToFloat(Float.floatToRawIntBits(x) ^ 1));
3465     }
3466 
3467     /**
3468      * Multiply a double number by a power of 2.
3469      * @param d number to multiply
3470      * @param n power of 2
3471      * @return d &times; 2<sup>n</sup>
3472      */
3473     public static double scalb(final double d, final int n) {
3474 
3475         // first simple and fast handling when 2^n can be represented using normal numbers
3476         if ((n > -1023) && (n < 1024)) {
3477             return d * Double.longBitsToDouble(((long) (n + 1023)) << 52);
3478         }
3479 
3480         // handle special cases
3481         if (Double.isNaN(d) || Double.isInfinite(d) || (d == 0)) {
3482             return d;
3483         }
3484         if (n < -2098) {
3485             return (d > 0) ? 0.0 : -0.0;
3486         }
3487         if (n > 2097) {
3488             return (d > 0) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3489         }
3490 
3491         // decompose d
3492         final long bits = Double.doubleToRawLongBits(d);
3493         final long sign = bits & 0x8000000000000000L;
3494         int  exponent   = ((int) (bits >>> 52)) & 0x7ff;
3495         long mantissa   = bits & 0x000fffffffffffffL;
3496 
3497         // compute scaled exponent
3498         int scaledExponent = exponent + n;
3499 
3500         if (n < 0) {
3501             // we are really in the case n <= -1023
3502             if (scaledExponent > 0) {
3503                 // both the input and the result are normal numbers, we only adjust the exponent
3504                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3505             } else if (scaledExponent > -53) {
3506                 // the input is a normal number and the result is a subnormal number
3507 
3508                 // recover the hidden mantissa bit
3509                 mantissa |= 1L << 52;
3510 
3511                 // scales down complete mantissa, hence losing least significant bits
3512                 final long mostSignificantLostBit = mantissa & (1L << (-scaledExponent));
3513                 mantissa >>>= 1 - scaledExponent;
3514                 if (mostSignificantLostBit != 0) {
3515                     // we need to add 1 bit to round up the result
3516                     mantissa++;
3517                 }
3518                 return Double.longBitsToDouble(sign | mantissa);
3519 
3520             } else {
3521                 // no need to compute the mantissa, the number scales down to 0
3522                 return (sign == 0L) ? 0.0 : -0.0;
3523             }
3524         } else {
3525             // we are really in the case n >= 1024
3526             if (exponent == 0) {
3527 
3528                 // the input number is subnormal, normalize it
3529                 while ((mantissa >>> 52) != 1) {
3530                     mantissa <<= 1;
3531                     --scaledExponent;
3532                 }
3533                 ++scaledExponent;
3534                 mantissa &= 0x000fffffffffffffL;
3535 
3536                 if (scaledExponent < 2047) {
3537                     return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3538                 } else {
3539                     return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3540                 }
3541 
3542             } else if (scaledExponent < 2047) {
3543                 return Double.longBitsToDouble(sign | (((long) scaledExponent) << 52) | mantissa);
3544             } else {
3545                 return (sign == 0L) ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
3546             }
3547         }
3548 
3549     }
3550 
3551     /**
3552      * Multiply a float number by a power of 2.
3553      * @param f number to multiply
3554      * @param n power of 2
3555      * @return f &times; 2<sup>n</sup>
3556      */
3557     public static float scalb(final float f, final int n) {
3558 
3559         // first simple and fast handling when 2^n can be represented using normal numbers
3560         if ((n > -127) && (n < 128)) {
3561             return f * Float.intBitsToFloat((n + 127) << 23);
3562         }
3563 
3564         // handle special cases
3565         if (Float.isNaN(f) || Float.isInfinite(f) || (f == 0f)) {
3566             return f;
3567         }
3568         if (n < -277) {
3569             return (f > 0) ? 0.0f : -0.0f;
3570         }
3571         if (n > 276) {
3572             return (f > 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3573         }
3574 
3575         // decompose f
3576         final int bits = Float.floatToIntBits(f);
3577         final int sign = bits & 0x80000000;
3578         int  exponent  = (bits >>> 23) & 0xff;
3579         int mantissa   = bits & 0x007fffff;
3580 
3581         // compute scaled exponent
3582         int scaledExponent = exponent + n;
3583 
3584         if (n < 0) {
3585             // we are really in the case n <= -127
3586             if (scaledExponent > 0) {
3587                 // both the input and the result are normal numbers, we only adjust the exponent
3588                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3589             } else if (scaledExponent > -24) {
3590                 // the input is a normal number and the result is a subnormal number
3591 
3592                 // recover the hidden mantissa bit
3593                 mantissa |= 1 << 23;
3594 
3595                 // scales down complete mantissa, hence losing least significant bits
3596                 final int mostSignificantLostBit = mantissa & (1 << (-scaledExponent));
3597                 mantissa >>>= 1 - scaledExponent;
3598                 if (mostSignificantLostBit != 0) {
3599                     // we need to add 1 bit to round up the result
3600                     mantissa++;
3601                 }
3602                 return Float.intBitsToFloat(sign | mantissa);
3603 
3604             } else {
3605                 // no need to compute the mantissa, the number scales down to 0
3606                 return (sign == 0) ? 0.0f : -0.0f;
3607             }
3608         } else {
3609             // we are really in the case n >= 128
3610             if (exponent == 0) {
3611 
3612                 // the input number is subnormal, normalize it
3613                 while ((mantissa >>> 23) != 1) {
3614                     mantissa <<= 1;
3615                     --scaledExponent;
3616                 }
3617                 ++scaledExponent;
3618                 mantissa &= 0x007fffff;
3619 
3620                 if (scaledExponent < 255) {
3621                     return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3622                 } else {
3623                     return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3624                 }
3625 
3626             } else if (scaledExponent < 255) {
3627                 return Float.intBitsToFloat(sign | (scaledExponent << 23) | mantissa);
3628             } else {
3629                 return (sign == 0) ? Float.POSITIVE_INFINITY : Float.NEGATIVE_INFINITY;
3630             }
3631         }
3632 
3633     }
3634 
3635     /**
3636      * Get the next machine representable number after a number, moving
3637      * in the direction of another number.
3638      * <p>
3639      * The ordering is as follows (increasing):
3640      * </p>
3641      * <ul>
3642      * <li>-INFINITY</li>
3643      * <li>-MAX_VALUE</li>
3644      * <li>-MIN_VALUE</li>
3645      * <li>-0.0</li>
3646      * <li>+0.0</li>
3647      * <li>+MIN_VALUE</li>
3648      * <li>+MAX_VALUE</li>
3649      * <li>+INFINITY</li>
3650      * </ul>
3651      * <p>
3652      * If arguments compare equal, then the second argument is returned.
3653      * </p>
3654      * <p>
3655      * If {@code direction} is greater than {@code d},
3656      * the smallest machine representable number strictly greater than
3657      * {@code d} is returned; if less, then the largest representable number
3658      * strictly less than {@code d} is returned.
3659      * </p>
3660      * <p>
3661      * If {@code d} is infinite and direction does not
3662      * bring it back to finite numbers, it is returned unchanged.
3663      * </p>
3664      *
3665      * @param d base number
3666      * @param direction (the only important thing is whether
3667      * {@code direction} is greater or smaller than {@code d})
3668      * @return the next machine representable number in the specified direction
3669      */
3670     public static double nextAfter(double d, double direction) {
3671 
3672         // handling of some important special cases
3673         if (Double.isNaN(d) || Double.isNaN(direction)) {
3674             return Double.NaN;
3675         } else if (d == direction) {
3676             return direction;
3677         } else if (Double.isInfinite(d)) {
3678             return (d < 0) ? -Double.MAX_VALUE : Double.MAX_VALUE;
3679         } else if (d == 0) {
3680             return (direction < 0) ? -Double.MIN_VALUE : Double.MIN_VALUE;
3681         }
3682         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3683         // are handled just as normal numbers
3684         // can use raw bits since already dealt with infinity and NaN
3685         final long bits = Double.doubleToRawLongBits(d);
3686         final long sign = bits & 0x8000000000000000L;
3687         if ((direction < d) ^ (sign == 0L)) {
3688             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) + 1));
3689         } else {
3690             return Double.longBitsToDouble(sign | ((bits & 0x7fffffffffffffffL) - 1));
3691         }
3692 
3693     }
3694 
3695     /**
3696      * Get the next machine representable number after a number, moving
3697      * in the direction of another number.
3698      * <p>* The ordering is as follows (increasing):</p>
3699      * <ul>
3700      * <li>-INFINITY</li>
3701      * <li>-MAX_VALUE</li>
3702      * <li>-MIN_VALUE</li>
3703      * <li>-0.0</li>
3704      * <li>+0.0</li>
3705      * <li>+MIN_VALUE</li>
3706      * <li>+MAX_VALUE</li>
3707      * <li>+INFINITY</li>
3708      * </ul>
3709      * <p>
3710      * If arguments compare equal, then the second argument is returned.
3711      * </p>
3712      * <p>
3713      * If {@code direction} is greater than {@code f},
3714      * the smallest machine representable number strictly greater than
3715      * {@code f} is returned; if less, then the largest representable number
3716      * strictly less than {@code f} is returned.
3717      * </p>
3718      * <p>
3719      * If {@code f} is infinite and direction does not
3720      * bring it back to finite numbers, it is returned unchanged.
3721      * </p>
3722      *
3723      * @param f base number
3724      * @param direction (the only important thing is whether
3725      * {@code direction} is greater or smaller than {@code f})
3726      * @return the next machine representable number in the specified direction
3727      */
3728     public static float nextAfter(final float f, final double direction) {
3729 
3730         // handling of some important special cases
3731         if (Double.isNaN(f) || Double.isNaN(direction)) {
3732             return Float.NaN;
3733         } else if (f == direction) {
3734             return (float) direction;
3735         } else if (Float.isInfinite(f)) {
3736             return (f < 0f) ? -Float.MAX_VALUE : Float.MAX_VALUE;
3737         } else if (f == 0f) {
3738             return (direction < 0) ? -Float.MIN_VALUE : Float.MIN_VALUE;
3739         }
3740         // special cases MAX_VALUE to infinity and  MIN_VALUE to 0
3741         // are handled just as normal numbers
3742 
3743         final int bits = Float.floatToIntBits(f);
3744         final int sign = bits & 0x80000000;
3745         if ((direction < f) ^ (sign == 0)) {
3746             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) + 1));
3747         } else {
3748             return Float.intBitsToFloat(sign | ((bits & 0x7fffffff) - 1));
3749         }
3750 
3751     }
3752 
3753     /** Get the largest whole number smaller than x.
3754      * @param x number from which floor is requested
3755      * @return a double number f such that f is an integer f &lt;= x &lt; f + 1.0
3756      */
3757     public static double floor(double x) {
3758         long y;
3759 
3760         if (Double.isNaN(x)) {
3761             return x;
3762         }
3763 
3764         if (x >= TWO_POWER_52 || x <= -TWO_POWER_52) {
3765             return x;
3766         }
3767 
3768         y = (long) x;
3769         if (x < 0 && y != x) {
3770             y--;
3771         }
3772 
3773         if (y == 0) {
3774             return x*y;
3775         }
3776 
3777         return y;
3778     }
3779 
3780     /** Get the smallest whole number larger than x.
3781      * @param x number from which ceil is requested
3782      * @return a double number c such that c is an integer c - 1.0 &lt; x &lt;= c
3783      */
3784     public static double ceil(double x) {
3785         double y;
3786 
3787         if (Double.isNaN(x)) {
3788             return x;
3789         }
3790 
3791         y = floor(x);
3792         if (y == x) {
3793             return y;
3794         }
3795 
3796         y += 1.0;
3797 
3798         if (y == 0) {
3799             return x*y;
3800         }
3801 
3802         return y;
3803     }
3804 
3805     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
3806      * @param x number from which nearest whole number is requested
3807      * @return a double number r such that r is an integer r - 0.5 &lt;= x &lt;= r + 0.5
3808      */
3809     public static double rint(double x) {
3810         double y = floor(x);
3811         double d = x - y;
3812 
3813         if (d > 0.5) {
3814             if (y == -1.0) {
3815                 return -0.0; // Preserve sign of operand
3816             }
3817             return y+1.0;
3818         }
3819         if (d < 0.5) {
3820             return y;
3821         }
3822 
3823         /* half way, round to even */
3824         long z = (long) y;
3825         return (z & 1) == 0 ? y : y + 1.0;
3826     }
3827 
3828     /** Get the closest long to x.
3829      * @param x number from which closest long is requested
3830      * @return closest long to x
3831      */
3832     public static long round(double x) {
3833         final long bits = Double.doubleToRawLongBits(x);
3834         final int biasedExp = ((int)(bits>>52)) & 0x7ff;
3835         // Shift to get rid of bits past comma except first one: will need to
3836         // 1-shift to the right to end up with correct magnitude.
3837         final int shift = (52 - 1 + Double.MAX_EXPONENT) - biasedExp;
3838         if ((shift & -64) == 0) {
3839             // shift in [0,63], so unbiased exp in [-12,51].
3840             long extendedMantissa = 0x0010000000000000L | (bits & 0x000fffffffffffffL);
3841             if (bits < 0) {
3842                 extendedMantissa = -extendedMantissa;
3843             }
3844             // If value is positive and first bit past comma is 0, rounding
3845             // to lower integer, else to upper one, which is what "+1" and
3846             // then ">>1" do.
3847             return ((extendedMantissa >> shift) + 1L) >> 1;
3848         } else {
3849             // +-Infinity, NaN, or a mathematical integer.
3850             return (long) x;
3851         }
3852     }
3853 
3854     /** Get the closest int to x.
3855      * @param x number from which closest int is requested
3856      * @return closest int to x
3857      */
3858     public static int round(final float x) {
3859         final int bits = Float.floatToRawIntBits(x);
3860         final int biasedExp = (bits>>23) & 0xff;
3861         // Shift to get rid of bits past comma except first one: will need to
3862         // 1-shift to the right to end up with correct magnitude.
3863         final int shift = (23 - 1 + Float.MAX_EXPONENT) - biasedExp;
3864         if ((shift & -32) == 0) {
3865             // shift in [0,31], so unbiased exp in [-9,22].
3866             int extendedMantissa = 0x00800000 | (bits & 0x007fffff);
3867             if (bits < 0) {
3868                 extendedMantissa = -extendedMantissa;
3869             }
3870             // If value is positive and first bit past comma is 0, rounding
3871             // to lower integer, else to upper one, which is what "+1" and
3872             // then ">>1" do.
3873             return ((extendedMantissa >> shift) + 1) >> 1;
3874         } else {
3875             // +-Infinity, NaN, or a mathematical integer.
3876             return (int) x;
3877         }
3878     }
3879 
3880     /** Compute the minimum of two values
3881      * @param a first value
3882      * @param b second value
3883      * @return a if a is lesser or equal to b, b otherwise
3884      */
3885     public static int min(final int a, final int b) {
3886         return (a <= b) ? a : b;
3887     }
3888 
3889     /** Compute the minimum of two values
3890      * @param a first value
3891      * @param b second value
3892      * @return a if a is lesser or equal to b, b otherwise
3893      */
3894     public static long min(final long a, final long b) {
3895         return (a <= b) ? a : b;
3896     }
3897 
3898     /** Compute the minimum of two values
3899      * @param a first value
3900      * @param b second value
3901      * @return a if a is lesser or equal to b, b otherwise
3902      */
3903     public static float min(final float a, final float b) {
3904         if (a > b) {
3905             return b;
3906         }
3907         if (a < b) {
3908             return a;
3909         }
3910         /* if either arg is NaN, return NaN */
3911         if (a != b) {
3912             return Float.NaN;
3913         }
3914         /* min(+0.0,-0.0) == -0.0 */
3915         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3916         int bits = Float.floatToRawIntBits(a);
3917         if (bits == 0x80000000) {
3918             return a;
3919         }
3920         return b;
3921     }
3922 
3923     /** Compute the minimum of two values
3924      * @param a first value
3925      * @param b second value
3926      * @return a if a is lesser or equal to b, b otherwise
3927      */
3928     public static double min(final double a, final double b) {
3929         if (a > b) {
3930             return b;
3931         }
3932         if (a < b) {
3933             return a;
3934         }
3935         /* if either arg is NaN, return NaN */
3936         if (a != b) {
3937             return Double.NaN;
3938         }
3939         /* min(+0.0,-0.0) == -0.0 */
3940         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
3941         long bits = Double.doubleToRawLongBits(a);
3942         if (bits == 0x8000000000000000L) {
3943             return a;
3944         }
3945         return b;
3946     }
3947 
3948     /** Compute the maximum of two values
3949      * @param a first value
3950      * @param b second value
3951      * @return b if a is lesser or equal to b, a otherwise
3952      */
3953     public static int max(final int a, final int b) {
3954         return (a <= b) ? b : a;
3955     }
3956 
3957     /** Compute the maximum of two values
3958      * @param a first value
3959      * @param b second value
3960      * @return b if a is lesser or equal to b, a otherwise
3961      */
3962     public static long max(final long a, final long b) {
3963         return (a <= b) ? b : a;
3964     }
3965 
3966     /** Compute the maximum of two values
3967      * @param a first value
3968      * @param b second value
3969      * @return b if a is lesser or equal to b, a otherwise
3970      */
3971     public static float max(final float a, final float b) {
3972         if (a > b) {
3973             return a;
3974         }
3975         if (a < b) {
3976             return b;
3977         }
3978         /* if either arg is NaN, return NaN */
3979         if (a != b) {
3980             return Float.NaN;
3981         }
3982         /* min(+0.0,-0.0) == -0.0 */
3983         /* 0x80000000 == Float.floatToRawIntBits(-0.0d) */
3984         int bits = Float.floatToRawIntBits(a);
3985         if (bits == 0x80000000) {
3986             return b;
3987         }
3988         return a;
3989     }
3990 
3991     /** Compute the maximum of two values
3992      * @param a first value
3993      * @param b second value
3994      * @return b if a is lesser or equal to b, a otherwise
3995      */
3996     public static double max(final double a, final double b) {
3997         if (a > b) {
3998             return a;
3999         }
4000         if (a < b) {
4001             return b;
4002         }
4003         /* if either arg is NaN, return NaN */
4004         if (a != b) {
4005             return Double.NaN;
4006         }
4007         /* min(+0.0,-0.0) == -0.0 */
4008         /* 0x8000000000000000L == Double.doubleToRawLongBits(-0.0d) */
4009         long bits = Double.doubleToRawLongBits(a);
4010         if (bits == 0x8000000000000000L) {
4011             return b;
4012         }
4013         return a;
4014     }
4015 
4016     /**
4017      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
4018      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br>
4019      * avoiding intermediate overflow or underflow.
4020      *
4021      * <ul>
4022      * <li> If either argument is infinite, then the result is positive infinity.</li>
4023      * <li> else, if either argument is NaN then the result is NaN.</li>
4024      * </ul>
4025      *
4026      * @param x a value
4027      * @param y a value
4028      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
4029      */
4030     public static double hypot(final double x, final double y) {
4031         if (Double.isInfinite(x) || Double.isInfinite(y)) {
4032             return Double.POSITIVE_INFINITY;
4033         } else if (Double.isNaN(x) || Double.isNaN(y)) {
4034             return Double.NaN;
4035         } else {
4036 
4037             final int expX = getExponent(x);
4038             final int expY = getExponent(y);
4039             if (expX > expY + 27) {
4040                 // y is neglectible with respect to x
4041                 return abs(x);
4042             } else if (expY > expX + 27) {
4043                 // x is neglectible with respect to y
4044                 return abs(y);
4045             } else {
4046 
4047                 // find an intermediate scale to avoid both overflow and underflow
4048                 final int middleExp = (expX + expY) / 2;
4049 
4050                 // scale parameters without losing precision
4051                 final double scaledX = scalb(x, -middleExp);
4052                 final double scaledY = scalb(y, -middleExp);
4053 
4054                 // compute scaled hypotenuse
4055                 final double scaledH = sqrt(scaledX * scaledX + scaledY * scaledY);
4056 
4057                 // remove scaling
4058                 return scalb(scaledH, middleExp);
4059 
4060             }
4061 
4062         }
4063     }
4064 
4065     /**
4066      * Computes the remainder as prescribed by the IEEE 754 standard.
4067      * <p>
4068      * The remainder value is mathematically equal to {@code x - y*n}
4069      * where {@code n} is the mathematical integer closest to the exact mathematical value
4070      * of the quotient {@code x/y}.
4071      * If two mathematical integers are equally close to {@code x/y} then
4072      * {@code n} is the integer that is even.
4073      * </p>
4074      * <ul>
4075      * <li>If either operand is NaN, the result is NaN.</li>
4076      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
4077      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
4078      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
4079      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
4080      * </ul>
4081      * @param dividend the number to be divided
4082      * @param divisor the number by which to divide
4083      * @return the remainder, rounded
4084      */
4085     public static double IEEEremainder(final double dividend, final double divisor) {
4086         if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
4087             // we are in one of the special cases
4088             if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
4089                 return dividend;
4090             } else {
4091                 return Double.NaN;
4092             }
4093         } else {
4094             // we are in the general case
4095             final double n         = FastMath.rint(dividend / divisor);
4096             final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
4097             return (remainder == 0) ? FastMath.copySign(remainder, dividend) : remainder;
4098         }
4099     }
4100 
4101     /** Convert a long to interger, detecting overflows
4102      * @param n number to convert to int
4103      * @return integer with same valie as n if no overflows occur
4104      * @exception MathRuntimeException if n cannot fit into an int
4105      */
4106     public static int toIntExact(final long n) throws MathRuntimeException {
4107         if (n < Integer.MIN_VALUE || n > Integer.MAX_VALUE) {
4108             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW);
4109         }
4110         return (int) n;
4111     }
4112 
4113     /** Increment a number, detecting overflows.
4114      * @param n number to increment
4115      * @return n+1 if no overflows occur
4116      * @exception MathRuntimeException if an overflow occurs
4117      */
4118     public static int incrementExact(final int n) throws MathRuntimeException {
4119 
4120         if (n == Integer.MAX_VALUE) {
4121             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
4122         }
4123 
4124         return n + 1;
4125 
4126     }
4127 
4128     /** Increment a number, detecting overflows.
4129      * @param n number to increment
4130      * @return n+1 if no overflows occur
4131      * @exception MathRuntimeException if an overflow occurs
4132      */
4133     public static long incrementExact(final long n) throws MathRuntimeException {
4134 
4135         if (n == Long.MAX_VALUE) {
4136             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, n, 1);
4137         }
4138 
4139         return n + 1;
4140 
4141     }
4142 
4143     /** Decrement a number, detecting overflows.
4144      * @param n number to decrement
4145      * @return n-1 if no overflows occur
4146      * @exception MathRuntimeException if an overflow occurs
4147      */
4148     public static int decrementExact(final int n) throws MathRuntimeException {
4149 
4150         if (n == Integer.MIN_VALUE) {
4151             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
4152         }
4153 
4154         return n - 1;
4155 
4156     }
4157 
4158     /** Decrement a number, detecting overflows.
4159      * @param n number to decrement
4160      * @return n-1 if no overflows occur
4161      * @exception MathRuntimeException if an overflow occurs
4162      */
4163     public static long decrementExact(final long n) throws MathRuntimeException {
4164 
4165         if (n == Long.MIN_VALUE) {
4166             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, n, 1);
4167         }
4168 
4169         return n - 1;
4170 
4171     }
4172 
4173     /** Add two numbers, detecting overflows.
4174      * @param a first number to add
4175      * @param b second number to add
4176      * @return a+b if no overflows occur
4177      * @exception MathRuntimeException if an overflow occurs
4178      */
4179     public static int addExact(final int a, final int b) throws MathRuntimeException {
4180 
4181         // compute sum
4182         final int sum = a + b;
4183 
4184         // check for overflow
4185         if ((a ^ b) >= 0 && (sum ^ b) < 0) {
4186             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
4187         }
4188 
4189         return sum;
4190 
4191     }
4192 
4193     /** Add two numbers, detecting overflows.
4194      * @param a first number to add
4195      * @param b second number to add
4196      * @return a+b if no overflows occur
4197      * @exception MathRuntimeException if an overflow occurs
4198      */
4199     public static long addExact(final long a, final long b) throws MathRuntimeException {
4200 
4201         // compute sum
4202         final long sum = a + b;
4203 
4204         // check for overflow
4205         if ((a ^ b) >= 0 && (sum ^ b) < 0) {
4206             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_ADDITION, a, b);
4207         }
4208 
4209         return sum;
4210 
4211     }
4212 
4213     /** Subtract two numbers, detecting overflows.
4214      * @param a first number
4215      * @param b second number to subtract from a
4216      * @return a-b if no overflows occur
4217      * @exception MathRuntimeException if an overflow occurs
4218      */
4219     public static int subtractExact(final int a, final int b) {
4220 
4221         // compute subtraction
4222         final int sub = a - b;
4223 
4224         // check for overflow
4225         if ((a ^ b) < 0 && (sub ^ b) >= 0) {
4226             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
4227         }
4228 
4229         return sub;
4230 
4231     }
4232 
4233     /** Subtract two numbers, detecting overflows.
4234      * @param a first number
4235      * @param b second number to subtract from a
4236      * @return a-b if no overflows occur
4237      * @exception MathRuntimeException if an overflow occurs
4238      */
4239     public static long subtractExact(final long a, final long b) {
4240 
4241         // compute subtraction
4242         final long sub = a - b;
4243 
4244         // check for overflow
4245         if ((a ^ b) < 0 && (sub ^ b) >= 0) {
4246             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_SUBTRACTION, a, b);
4247         }
4248 
4249         return sub;
4250 
4251     }
4252 
4253     /** Multiply two numbers, detecting overflows.
4254      * @param a first number to multiply
4255      * @param b second number to multiply
4256      * @return a*b if no overflows occur
4257      * @exception MathRuntimeException if an overflow occurs
4258      */
4259     public static int multiplyExact(final int a, final int b) {
4260         if (((b  >  0)  && (a > Integer.MAX_VALUE / b || a < Integer.MIN_VALUE / b)) ||
4261             ((b  < -1)  && (a > Integer.MIN_VALUE / b || a < Integer.MAX_VALUE / b)) ||
4262             ((b == -1)  && (a == Integer.MIN_VALUE))) {
4263             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
4264         }
4265         return a * b;
4266     }
4267 
4268     /** Multiply two numbers, detecting overflows.
4269      * @param a first number to multiply
4270      * @param b second number to multiply
4271      * @return a*b if no overflows occur
4272      * @exception MathRuntimeException if an overflow occurs
4273      * @since 1.3
4274      */
4275     public static long multiplyExact(final long a, final int b) {
4276         return multiplyExact(a, (long) b);
4277     }
4278 
4279     /** Multiply two numbers, detecting overflows.
4280      * @param a first number to multiply
4281      * @param b second number to multiply
4282      * @return a*b if no overflows occur
4283      * @exception MathRuntimeException if an overflow occurs
4284      */
4285     public static long multiplyExact(final long a, final long b) {
4286         if (((b  >  0l)  && (a > Long.MAX_VALUE / b || a < Long.MIN_VALUE / b)) ||
4287             ((b  < -1l)  && (a > Long.MIN_VALUE / b || a < Long.MAX_VALUE / b)) ||
4288             ((b == -1l)  && (a == Long.MIN_VALUE))) {
4289                 throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_MULTIPLICATION, a, b);
4290             }
4291             return a * b;
4292     }
4293 
4294     /** Multiply two integers and give an exact result without overflow.
4295      * @param a first factor
4296      * @param b second factor
4297      * @return a * b exactly
4298      * @since 1.3
4299      */
4300     public static long multiplyFull(final int a, final int b) {
4301         return ((long) a) * ((long) b);
4302     }
4303 
4304     /** Multiply two long integers and give the 64 most significant bits of the result.
4305      * <p>
4306      * Beware that as Java primitive long are always considered to be signed, there are some
4307      * intermediate values {@code a} and {@code b} for which {@code a * b} exceeds {@code Long.MAX_VALUE}
4308      * but this method will still return 0l. This happens for example for {@code a = 2³¹} and
4309      * {@code b = 2³²} as {@code a * b = 2⁶³ = Long.MAX_VALUE + 1}, so it exceeds the max value
4310      * for a long, but still fits in 64 bits, so this method correctly returns 0l in this case,
4311      * but multiplication result would be considered negative (and in fact equal to {@code Long.MIN_VALUE}
4312      * </p>
4313      * @param a first factor
4314      * @param b second factor
4315      * @return a * b / 2<sup>64</sup>
4316      * @since 1.3
4317      */
4318     public static long multiplyHigh(final long a, final long b) {
4319 
4320         // all computations below are performed on unsigned numbers because we start
4321         // by using logical shifts (and not arithmetic shifts). We will therefore
4322         // need to take care of sign before returning
4323         // a negative long n between -2⁶³ and -1, interpreted as an unsigned long
4324         // corresponds to 2⁶⁴ + n (which is between 2⁶³ and 2⁶⁴-1)
4325         // so if this number is multiplied by p, what we really compute
4326         // is (2⁶⁴ + n) * p = 2⁶⁴ * p + n * p, therefore the part above 64 bits
4327         // will have an extra term p that we will need to remove
4328         final long tobeRemoved = ((a < 0) ? b : 0) + ((b < 0) ? a : 0);
4329 
4330         return unsignedMultiplyHigh(a, b) - tobeRemoved;
4331 
4332     }
4333 
4334     /** Multiply two long unsigned integers and give the 64 most significant bits of the unsigned result.
4335      * <p>
4336      * Beware that as Java primitive long are always considered to be signed, there are some
4337      * intermediate values {@code a} and {@code b} for which {@code a * b} exceeds {@code Long.MAX_VALUE}
4338      * but this method will still return 0l. This happens for example for {@code a = 2³¹} and
4339      * {@code b = 2³²} as {@code a * b = 2⁶³ = Long.MAX_VALUE + 1}, so it exceeds the max value
4340      * for a long, but still fits in 64 bits, so this method correctly returns 0l in this case,
4341      * but multiplication result would be considered negative (and in fact equal to {@code Long.MIN_VALUE}
4342      * </p>
4343      * @param a first factor
4344      * @param b second factor
4345      * @return a * b / 2<sup>64</sup>
4346      * @since 3.0
4347      */
4348     public static long unsignedMultiplyHigh(final long a, final long b) {
4349 
4350         // split numbers in two 32 bits parts
4351         final long aHigh  = a >>> 32;
4352         final long aLow   = a & 0xFFFFFFFFl;
4353         final long bHigh  = b >>> 32;
4354         final long bLow   = b & 0xFFFFFFFFl;
4355 
4356         // ab = aHigh * bHigh * 2⁶⁴ + (aHigh * bLow + aLow * bHigh) * 2³² + aLow * bLow
4357         final long hh     = aHigh * bHigh;
4358         final long hl1    = aHigh * bLow;
4359         final long hl2    = aLow  * bHigh;
4360         final long ll     = aLow  * bLow;
4361 
4362         // adds up everything in the above 64 bit part, taking care to avoid overflow
4363         final long hlHigh = (hl1 >>> 32) + (hl2 >>> 32);
4364         final long hlLow  = (hl1 & 0xFFFFFFFFl) + (hl2 & 0xFFFFFFFFl);
4365         final long carry  = (hlLow + (ll >>> 32)) >>> 32;
4366 
4367         return hh + hlHigh + carry;
4368 
4369     }
4370 
4371     /** Divide two integers, checking for overflow.
4372      * @param x dividend
4373      * @param y divisor
4374      * @return x / y
4375      * @exception MathRuntimeException if an overflow occurs
4376      * @since 3.0
4377      */
4378     public static int divideExact(final int x, final int y) {
4379         if (y == 0) {
4380             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4381         }
4382         if (y == -1 && x == Integer.MIN_VALUE) {
4383             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
4384         }
4385         return x / y;
4386     }
4387 
4388     /** Divide two long integers, checking for overflow.
4389      * @param x dividend
4390      * @param y divisor
4391      * @return x / y
4392      * @exception MathRuntimeException if an overflow occurs
4393      * @since 3.0
4394      */
4395     public static long divideExact(final long x, final long y) {
4396         if (y == 0l) {
4397             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4398         }
4399         if (y == -1l && x == Long.MIN_VALUE) {
4400             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, x, y);
4401         }
4402         return x / y;
4403     }
4404 
4405     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4406      * <p>
4407      * This methods returns the same value as integer division when
4408      * a and b are opposite signs, but returns a different value when
4409      * they are same (i.e. q is positive).
4410      *
4411      * @param a dividend
4412      * @param b divisor
4413      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4414      * @exception MathRuntimeException if b == 0
4415      * @since 3.0
4416      */
4417     public static int ceilDiv(final int a, final int b) throws MathRuntimeException {
4418 
4419         if (b == 0) {
4420             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4421         }
4422 
4423         final int m = a % b;
4424         if ((a ^ b) < 0 || m == 0) {
4425             // a and b have opposite signs, or division is exact
4426             return a / b;
4427         } else {
4428             // a and b have same signs and division is not exact
4429             return (a / b) + 1;
4430         }
4431 
4432     }
4433 
4434     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4435      * <p>
4436      * This methods returns the same value as integer division when
4437      * a and b are opposite signs, but returns a different value when
4438      * they are same (i.e. q is positive).
4439      *
4440      * @param a dividend
4441      * @param b divisor
4442      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4443      * @exception MathRuntimeException if b == 0 or if a == {@code Integer.MIN_VALUE} and b = -1
4444      * @since 3.0
4445      */
4446     public static int ceilDivExact(final int a, final int b) throws MathRuntimeException {
4447 
4448         if (a == Integer.MIN_VALUE && b == -1) {
4449             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4450         }
4451 
4452         return ceilDiv(a, b);
4453 
4454     }
4455 
4456     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4457      * <p>
4458      * This methods returns the same value as integer division when
4459      * a and b are opposite signs, but returns a different value when
4460      * they are same (i.e. q is positive).
4461      *
4462      * @param a dividend
4463      * @param b divisor
4464      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4465      * @exception MathRuntimeException if b == 0
4466      * @since 3.0
4467      */
4468     public static long ceilDiv(final long a, final long b) throws MathRuntimeException {
4469 
4470         if (b == 0l) {
4471             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4472         }
4473 
4474         final long m = a % b;
4475         if ((a ^ b) < 0 || m == 0l) {
4476             // a and b have opposite signs, or division is exact
4477             return a / b;
4478         } else {
4479             // a and b have same signs and division is not exact
4480             return (a / b) + 1l;
4481         }
4482 
4483     }
4484 
4485     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4486      * <p>
4487      * This methods returns the same value as integer division when
4488      * a and b are opposite signs, but returns a different value when
4489      * they are same (i.e. q is positive).
4490      *
4491      * @param a dividend
4492      * @param b divisor
4493      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4494      * @exception MathRuntimeException if b == 0 or if a == {@code Long.MIN_VALUE} and b = -1
4495      * @since 3.0
4496      */
4497     public static long ceilDivExact(final long a, final long b) throws MathRuntimeException {
4498 
4499         if (a == Long.MIN_VALUE && b == -1l) {
4500             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4501         }
4502 
4503         return ceilDiv(a, b);
4504 
4505     }
4506 
4507     /** Finds q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4508      * <p>
4509      * This methods returns the same value as integer division when
4510      * a and b are opposite signs, but returns a different value when
4511      * they are same (i.e. q is positive).
4512      *
4513      * @param a dividend
4514      * @param b divisor
4515      * @return q such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4516      * @exception MathRuntimeException if b == 0
4517      * @since 3.0
4518      */
4519     public static long ceilDiv(final long a, final int b) throws MathRuntimeException {
4520         return ceilDiv(a, (long) b);
4521     }
4522 
4523     /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4524      * <p>
4525      * This methods returns the same value as integer modulo when
4526      * a and b are opposite signs, but returns a different value when
4527      * they are same (i.e. q is positive).
4528      *
4529      * @param a dividend
4530      * @param b divisor
4531      * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4532      * @exception MathRuntimeException if b == 0
4533      * @since 3.0
4534      */
4535     public static int ceilMod(final int a, final int b) throws MathRuntimeException {
4536 
4537         if (b == 0) {
4538             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4539         }
4540 
4541         final int m = a % b;
4542         if ((a ^ b) < 0 || m == 0) {
4543             // a and b have opposite signs, or division is exact
4544             return m;
4545         } else {
4546             // a and b have same signs and division is not exact
4547             return m - b;
4548         }
4549 
4550     }
4551 
4552     /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4553      * <p>
4554      * This methods returns the same value as integer modulo when
4555      * a and b are opposite signs, but returns a different value when
4556      * they are same (i.e. q is positive).
4557      *
4558      * @param a dividend
4559      * @param b divisor
4560      * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4561      * @exception MathRuntimeException if b == 0
4562      * @since 3.0
4563      */
4564     public static int ceilMod(final long a, final int b) throws MathRuntimeException {
4565         return (int) ceilMod(a, (long) b);
4566     }
4567 
4568     /** Finds r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}.
4569      * <p>
4570      * This methods returns the same value as integer modulo when
4571      * a and b are opposite signs, but returns a different value when
4572      * they are same (i.e. q is positive).
4573      *
4574      * @param a dividend
4575      * @param b divisor
4576      * @return r such that {@code a = q b + r} with {@code b < r <= 0} if {@code b > 0} and {@code 0 <= r < b} if {@code  b < 0}
4577      * @exception MathRuntimeException if b == 0
4578      * @since 3.0
4579      */
4580     public static long ceilMod(final long a, final long b) throws MathRuntimeException {
4581 
4582         if (b == 0l) {
4583             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4584         }
4585 
4586         final long m = a % b;
4587         if ((a ^ b) < 0l || m == 0l) {
4588             // a and b have opposite signs, or division is exact
4589             return m;
4590         } else {
4591             // a and b have same signs and division is not exact
4592             return m - b;
4593         }
4594 
4595     }
4596 
4597     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}.
4598      * <p>
4599      * This methods returns the same value as integer division when
4600      * a and b are same signs, but returns a different value when
4601      * they are opposite (i.e. q is negative).
4602      *
4603      * @param a dividend
4604      * @param b divisor
4605      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}
4606      * @exception MathRuntimeException if b == 0
4607      * @see #floorMod(int, int)
4608      */
4609     public static int floorDiv(final int a, final int b) throws MathRuntimeException {
4610 
4611         if (b == 0) {
4612             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4613         }
4614 
4615         final int m = a % b;
4616         if ((a ^ b) >= 0 || m == 0) {
4617             // a and b have same sign, or division is exact
4618             return a / b;
4619         } else {
4620             // a and b have opposite signs and division is not exact
4621             return (a / b) - 1;
4622         }
4623 
4624     }
4625 
4626     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}.
4627      * <p>
4628      * This methods returns the same value as integer division when
4629      * a and b are same signs, but returns a different value when
4630      * they are opposite (i.e. q is negative).
4631      *
4632      * @param a dividend
4633      * @param b divisor
4634      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code  b < 0}
4635      * @exception MathRuntimeException if b == 0 or if a == {@code Integer.MIN_VALUE} and b = -1
4636      * @see #floorMod(int, int)
4637      * @since 3.0
4638      */
4639     public static int floorDivExact(final int a, final int b) throws MathRuntimeException {
4640 
4641         if (a == Integer.MIN_VALUE && b == -1) {
4642             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4643         }
4644 
4645         return floorDiv(a, b);
4646 
4647     }
4648 
4649     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4650      * <p>
4651      * This methods returns the same value as integer division when
4652      * a and b are same signs, but returns a different value when
4653      * they are opposite (i.e. q is negative).
4654      *
4655      * @param a dividend
4656      * @param b divisor
4657      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4658      * @exception MathRuntimeException if b == 0
4659      * @see #floorMod(long, int)
4660      * @since 1.3
4661      */
4662     public static long floorDiv(final long a, final int b) throws MathRuntimeException {
4663         return floorDiv(a, (long) b);
4664     }
4665 
4666     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4667      * <p>
4668      * This methods returns the same value as integer division when
4669      * a and b are same signs, but returns a different value when
4670      * they are opposite (i.e. q is negative).
4671      *
4672      * @param a dividend
4673      * @param b divisor
4674      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4675      * @exception MathRuntimeException if b == 0
4676      * @see #floorMod(long, long)
4677      */
4678     public static long floorDiv(final long a, final long b) throws MathRuntimeException {
4679 
4680         if (b == 0l) {
4681             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4682         }
4683 
4684         final long m = a % b;
4685         if ((a ^ b) >= 0l || m == 0l) {
4686             // a and b have same sign, or division is exact
4687             return a / b;
4688         } else {
4689             // a and b have opposite signs and division is not exact
4690             return (a / b) - 1l;
4691         }
4692 
4693     }
4694 
4695     /** Finds q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4696      * <p>
4697      * This methods returns the same value as integer division when
4698      * a and b are same signs, but returns a different value when
4699      * they are opposite (i.e. q is negative).
4700      *
4701      * @param a dividend
4702      * @param b divisor
4703      * @return q such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4704      * @exception MathRuntimeException if b == 0 or if a == {@code Long.MIN_VALUE} and b = -1
4705      * @see #floorMod(long, long)
4706      * @since 3.0
4707      */
4708     public static long floorDivExact(final long a, final long b) throws MathRuntimeException {
4709 
4710         if (a == Long.MIN_VALUE && b == -1l) {
4711             throw new MathRuntimeException(LocalizedCoreFormats.OVERFLOW_IN_FRACTION, a, b);
4712         }
4713 
4714         return floorDiv(a, b);
4715 
4716     }
4717 
4718     /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4719      * <p>
4720      * This methods returns the same value as integer modulo when
4721      * a and b are same signs, but returns a different value when
4722      * they are opposite (i.e. q is negative).
4723      * </p>
4724      * @param a dividend
4725      * @param b divisor
4726      * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4727      * @exception MathRuntimeException if b == 0
4728      * @see #floorDiv(int, int)
4729      */
4730     public static int floorMod(final int a, final int b) throws MathRuntimeException {
4731 
4732         if (b == 0) {
4733             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4734         }
4735 
4736         final int m = a % b;
4737         if ((a ^ b) >= 0 || m == 0) {
4738             // a and b have same sign, or division is exact
4739             return m;
4740         } else {
4741             // a and b have opposite signs and division is not exact
4742             return b + m;
4743         }
4744 
4745     }
4746 
4747     /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4748      * <p>
4749      * This methods returns the same value as integer modulo when
4750      * a and b are same signs, but returns a different value when
4751      * they are opposite (i.e. q is negative).
4752      * </p>
4753      * @param a dividend
4754      * @param b divisor
4755      * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4756      * @exception MathRuntimeException if b == 0
4757      * @see #floorDiv(long, int)
4758      * @since 1.3
4759      */
4760     public static int floorMod(final long a, final int b) {
4761         return (int) floorMod(a, (long) b);
4762     }
4763 
4764     /** Finds r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}.
4765      * <p>
4766      * This methods returns the same value as integer modulo when
4767      * a and b are same signs, but returns a different value when
4768      * they are opposite (i.e. q is negative).
4769      * </p>
4770      * @param a dividend
4771      * @param b divisor
4772      * @return r such that {@code a = q b + r} with {@code 0 <= r < b} if {@code b > 0} and {@code b < r <= 0} if {@code b < 0}
4773      * @exception MathRuntimeException if b == 0
4774      * @see #floorDiv(long, long)
4775      */
4776     public static long floorMod(final long a, final long b) {
4777 
4778         if (b == 0l) {
4779             throw new MathRuntimeException(LocalizedCoreFormats.ZERO_DENOMINATOR);
4780         }
4781 
4782         final long m = a % b;
4783         if ((a ^ b) >= 0l || m == 0l) {
4784             // a and b have same sign, or division is exact
4785             return m;
4786         } else {
4787             // a and b have opposite signs and division is not exact
4788             return b + m;
4789         }
4790 
4791     }
4792 
4793     /**
4794      * Returns the first argument with the sign of the second argument.
4795      * A NaN {@code sign} argument is treated as positive.
4796      *
4797      * @param magnitude the value to return
4798      * @param sign the sign for the returned value
4799      * @return the magnitude with the same sign as the {@code sign} argument
4800      */
4801     public static double copySign(double magnitude, double sign){
4802         // The highest order bit is going to be zero if the
4803         // highest order bit of m and s is the same and one otherwise.
4804         // So (m^s) will be positive if both m and s have the same sign
4805         // and negative otherwise.
4806         final long m = Double.doubleToRawLongBits(magnitude); // don't care about NaN
4807         final long s = Double.doubleToRawLongBits(sign);
4808         if ((m ^ s) >= 0) {
4809             return magnitude;
4810         }
4811         return -magnitude; // flip sign
4812     }
4813 
4814     /**
4815      * Returns the first argument with the sign of the second argument.
4816      * A NaN {@code sign} argument is treated as positive.
4817      *
4818      * @param magnitude the value to return
4819      * @param sign the sign for the returned value
4820      * @return the magnitude with the same sign as the {@code sign} argument
4821      */
4822     public static float copySign(float magnitude, float sign){
4823         // The highest order bit is going to be zero if the
4824         // highest order bit of m and s is the same and one otherwise.
4825         // So (m^s) will be positive if both m and s have the same sign
4826         // and negative otherwise.
4827         final int m = Float.floatToRawIntBits(magnitude);
4828         final int s = Float.floatToRawIntBits(sign);
4829         if ((m ^ s) >= 0) {
4830             return magnitude;
4831         }
4832         return -magnitude; // flip sign
4833     }
4834 
4835     /**
4836      * Return the exponent of a double number, removing the bias.
4837      * <p>
4838      * For double numbers of the form 2<sup>x</sup>, the unbiased
4839      * exponent is exactly x.
4840      * </p>
4841      * @param d number from which exponent is requested
4842      * @return exponent for d in IEEE754 representation, without bias
4843      */
4844     public static int getExponent(final double d) {
4845         // NaN and Infinite will return 1024 anyhow so can use raw bits
4846         return (int) ((Double.doubleToRawLongBits(d) >>> 52) & 0x7ff) - 1023;
4847     }
4848 
4849     /**
4850      * Return the exponent of a float number, removing the bias.
4851      * <p>
4852      * For float numbers of the form 2<sup>x</sup>, the unbiased
4853      * exponent is exactly x.
4854      * </p>
4855      * @param f number from which exponent is requested
4856      * @return exponent for d in IEEE754 representation, without bias
4857      */
4858     public static int getExponent(final float f) {
4859         // NaN and Infinite will return the same exponent anyhow so can use raw bits
4860         return ((Float.floatToRawIntBits(f) >>> 23) & 0xff) - 127;
4861     }
4862 
4863     /** Compute Fused-multiply-add operation a * b + c.
4864      * <p>
4865      * This method was introduced in the regular {@code Math} and {@code StrictMath}
4866      * methods with Java 9, and then added to Hipparchus for consistency. However,
4867      * a more general method was available in Hipparchus that also allow to repeat
4868      * this computation across several terms: {@link MathArrays#linearCombination(double[], double[])}.
4869      * The linear combination method should probably be preferred in most cases.
4870      * </p>
4871      * @param a first factor
4872      * @param b second factor
4873      * @param c additive term
4874      * @return a * b + c, using extended precision in the multiplication
4875      * @see MathArrays#linearCombination(double[], double[])
4876      * @see MathArrays#linearCombination(double, double, double, double)
4877      * @see MathArrays#linearCombination(double, double, double, double, double, double)
4878      * @see MathArrays#linearCombination(double, double, double, double, double, double, double, double)
4879      * @since 1.3
4880      */
4881     public static double fma(final double a, final double b, final double c) {
4882         return MathArrays.linearCombination(a, b, 1.0, c);
4883     }
4884 
4885     /** Compute Fused-multiply-add operation a * b + c.
4886      * <p>
4887      * This method was introduced in the regular {@code Math} and {@code StrictMath}
4888      * methods with Java 9, and then added to Hipparchus for consistency. However,
4889      * a more general method was available in Hipparchus that also allow to repeat
4890      * this computation across several terms: {@link MathArrays#linearCombination(double[], double[])}.
4891      * The linear combination method should probably be preferred in most cases.
4892      * </p>
4893      * @param a first factor
4894      * @param b second factor
4895      * @param c additive term
4896      * @return a * b + c, using extended precision in the multiplication
4897      * @see MathArrays#linearCombination(double[], double[])
4898      * @see MathArrays#linearCombination(double, double, double, double)
4899      * @see MathArrays#linearCombination(double, double, double, double, double, double)
4900      * @see MathArrays#linearCombination(double, double, double, double, double, double, double, double)
4901      */
4902     public static float fma(final float a, final float b, final float c) {
4903         return (float) MathArrays.linearCombination(a, b, 1.0, c);
4904     }
4905 
4906     /** Compute the square root of a number.
4907      * @param a number on which evaluation is done
4908      * @param <T> the type of the field element
4909      * @return square root of a
4910      * @since 1.3
4911      */
4912     public static <T extends CalculusFieldElement<T>> T sqrt(final T a) {
4913         return a.sqrt();
4914     }
4915 
4916     /** Compute the hyperbolic cosine of a number.
4917      * @param x number on which evaluation is done
4918      * @param <T> the type of the field element
4919      * @return hyperbolic cosine of x
4920      * @since 1.3
4921      */
4922     public static <T extends CalculusFieldElement<T>> T cosh(final T x) {
4923         return x.cosh();
4924     }
4925 
4926     /** Compute the hyperbolic sine of a number.
4927      * @param x number on which evaluation is done
4928      * @param <T> the type of the field element
4929      * @return hyperbolic sine of x
4930      * @since 1.3
4931      */
4932     public static <T extends CalculusFieldElement<T>> T sinh(final T x) {
4933         return x.sinh();
4934     }
4935 
4936     /** Compute the hyperbolic tangent of a number.
4937      * @param x number on which evaluation is done
4938      * @param <T> the type of the field element
4939      * @return hyperbolic tangent of x
4940      * @since 1.3
4941      */
4942     public static <T extends CalculusFieldElement<T>> T tanh(final T x) {
4943         return x.tanh();
4944     }
4945 
4946     /** Compute the inverse hyperbolic cosine of a number.
4947      * @param a number on which evaluation is done
4948      * @param <T> the type of the field element
4949      * @return inverse hyperbolic cosine of a
4950      * @since 1.3
4951      */
4952     public static <T extends CalculusFieldElement<T>> T acosh(final T a) {
4953         return a.acosh();
4954     }
4955 
4956     /** Compute the inverse hyperbolic sine of a number.
4957      * @param a number on which evaluation is done
4958      * @param <T> the type of the field element
4959      * @return inverse hyperbolic sine of a
4960      * @since 1.3
4961      */
4962     public static <T extends CalculusFieldElement<T>> T asinh(final T a) {
4963         return a.asinh();
4964     }
4965 
4966     /** Compute the inverse hyperbolic tangent of a number.
4967      * @param a number on which evaluation is done
4968      * @param <T> the type of the field element
4969      * @return inverse hyperbolic tangent of a
4970      * @since 1.3
4971      */
4972     public static <T extends CalculusFieldElement<T>> T atanh(final T a) {
4973         return a.atanh();
4974     }
4975 
4976     /** Compute the sign of a number.
4977      * The sign is -1 for negative numbers, +1 for positive numbers and 0 otherwise,
4978      * for Complex number, it is extended on the unit circle (equivalent to z/|z|,
4979      * with special handling for 0 and NaN)
4980      * @param a number on which evaluation is done
4981      * @param <T> the type of the field element
4982      * @return -1.0, -0.0, +0.0, +1.0 or NaN depending on sign of a
4983      * @since 2.0
4984      */
4985     public static <T extends CalculusFieldElement<T>> T sign(final T a) {
4986         return a.sign();
4987     }
4988 
4989     /**
4990      * Exponential function.
4991      *
4992      * Computes exp(x), function result is nearly rounded.   It will be correctly
4993      * rounded to the theoretical value for 99.9% of input values, otherwise it will
4994      * have a 1 ULP error.
4995      *
4996      * Method:
4997      *    Lookup intVal = exp(int(x))
4998      *    Lookup fracVal = exp(int(x-int(x) / 1024.0) * 1024.0 );
4999      *    Compute z as the exponential of the remaining bits by a polynomial minus one
5000      *    exp(x) = intVal * fracVal * (1 + z)
5001      *
5002      * Accuracy:
5003      *    Calculation is done with 63 bits of precision, so result should be correctly
5004      *    rounded for 99.9% of input values, with less than 1 ULP error otherwise.
5005      *
5006      * @param x   a double
5007      * @param <T> the type of the field element
5008      * @return double e<sup>x</sup>
5009      * @since 1.3
5010      */
5011     public static <T extends CalculusFieldElement<T>> T exp(final T x) {
5012         return x.exp();
5013     }
5014 
5015     /** Compute exp(x) - 1
5016      * @param x number to compute shifted exponential
5017      * @param <T> the type of the field element
5018      * @return exp(x) - 1
5019      * @since 1.3
5020      */
5021     public static <T extends CalculusFieldElement<T>> T expm1(final T x) {
5022         return x.expm1();
5023     }
5024 
5025     /**
5026      * Natural logarithm.
5027      *
5028      * @param x   a double
5029      * @param <T> the type of the field element
5030      * @return log(x)
5031      * @since 1.3
5032      */
5033     public static <T extends CalculusFieldElement<T>> T log(final T x) {
5034         return x.log();
5035     }
5036 
5037     /**
5038      * Computes log(1 + x).
5039      *
5040      * @param x Number.
5041      * @param <T> the type of the field element
5042      * @return {@code log(1 + x)}.
5043      * @since 1.3
5044      */
5045     public static <T extends CalculusFieldElement<T>> T log1p(final T x) {
5046         return x.log1p();
5047     }
5048 
5049     /** Compute the base 10 logarithm.
5050      * @param x a number
5051      * @param <T> the type of the field element
5052      * @return log10(x)
5053      * @since 1.3
5054      */
5055     public static <T extends CalculusFieldElement<T>> T log10(final T x) {
5056         return x.log10();
5057     }
5058 
5059     /**
5060      * Power function.  Compute x<sup>y</sup>.
5061      *
5062      * @param x   a double
5063      * @param y   a double
5064      * @param <T> the type of the field element
5065      * @return x<sup>y</sup>
5066      * @since 1.3
5067      */
5068     public static <T extends CalculusFieldElement<T>> T pow(final T x, final T y) {
5069         return x.pow(y);
5070     }
5071 
5072     /**
5073      * Power function.  Compute x<sup>y</sup>.
5074      *
5075      * @param x   a double
5076      * @param y   a double
5077      * @param <T> the type of the field element
5078      * @return x<sup>y</sup>
5079      * @since 1.7
5080      */
5081     public static <T extends CalculusFieldElement<T>> T pow(final T x, final double y) {
5082         return x.pow(y);
5083     }
5084 
5085     /**
5086      * Raise a double to an int power.
5087      *
5088      * @param d Number to raise.
5089      * @param e Exponent.
5090      * @param <T> the type of the field element
5091      * @return d<sup>e</sup>
5092      * @since 1.3
5093      */
5094     public static <T extends CalculusFieldElement<T>> T pow(T d, int e) {
5095         return d.pow(e);
5096     }
5097 
5098     /**
5099      * Sine function.
5100      *
5101      * @param x Argument.
5102      * @param <T> the type of the field element
5103      * @return sin(x)
5104      * @since 1.3
5105      */
5106     public static <T extends CalculusFieldElement<T>> T sin(final T x) {
5107         return x.sin();
5108     }
5109 
5110     /**
5111      * Cosine function.
5112      *
5113      * @param x Argument.
5114      * @param <T> the type of the field element
5115      * @return cos(x)
5116      * @since 1.3
5117      */
5118     public static <T extends CalculusFieldElement<T>> T cos(final T x) {
5119         return x.cos();
5120     }
5121 
5122     /**
5123      * Tangent function.
5124      *
5125      * @param x Argument.
5126      * @param <T> the type of the field element
5127      * @return tan(x)
5128      * @since 1.3
5129      */
5130     public static <T extends CalculusFieldElement<T>> T tan(final T x) {
5131         return x.tan();
5132     }
5133 
5134     /**
5135      * Arctangent function
5136      *  @param x a number
5137      * @param <T> the type of the field element
5138      *  @return atan(x)
5139      * @since 1.3
5140      */
5141     public static <T extends CalculusFieldElement<T>> T atan(final T x) {
5142         return x.atan();
5143     }
5144 
5145     /**
5146      * Two arguments arctangent function
5147      * @param y ordinate
5148      * @param x abscissa
5149      * @param <T> the type of the field element
5150      * @return phase angle of point (x,y) between {@code -PI} and {@code PI}
5151      * @since 1.3
5152      */
5153     public static <T extends CalculusFieldElement<T>> T atan2(final T y, final T x) {
5154         return y.atan2(x);
5155     }
5156 
5157     /** Compute the arc sine of a number.
5158      * @param x number on which evaluation is done
5159      * @param <T> the type of the field element
5160      * @return arc sine of x
5161      * @since 1.3
5162      */
5163     public static <T extends CalculusFieldElement<T>> T asin(final T x) {
5164         return x.asin();
5165     }
5166 
5167     /** Compute the arc cosine of a number.
5168      * @param x number on which evaluation is done
5169      * @param <T> the type of the field element
5170      * @return arc cosine of x
5171      * @since 1.3
5172      */
5173     public static <T extends CalculusFieldElement<T>> T acos(final T x) {
5174         return x.acos();
5175     }
5176 
5177     /** Compute the cubic root of a number.
5178      * @param x number on which evaluation is done
5179      * @param <T> the type of the field element
5180      * @return cubic root of x
5181      * @since 1.3
5182      */
5183     public static <T extends CalculusFieldElement<T>> T cbrt(final T x) {
5184         return x.cbrt();
5185     }
5186 
5187     /**
5188      * Norm.
5189      * @param x number from which norm is requested
5190      * @param <T> the type of the field element
5191      * @return norm(x)
5192      * @since 2.0
5193      */
5194     public static <T extends CalculusFieldElement<T>> double norm(final T x) {
5195         return x.norm();
5196     }
5197 
5198     /**
5199      * Absolute value.
5200      * @param x number from which absolute value is requested
5201      * @param <T> the type of the field element
5202      * @return abs(x)
5203      * @since 2.0
5204      */
5205     public static <T extends CalculusFieldElement<T>> T abs(final T x) {
5206         return x.abs();
5207     }
5208 
5209     /**
5210      *  Convert degrees to radians, with error of less than 0.5 ULP
5211      *  @param x angle in degrees
5212      *  @param <T> the type of the field element
5213      *  @return x converted into radians
5214      */
5215     public static <T extends CalculusFieldElement<T>> T toRadians(T x) {
5216         return x.toRadians();
5217     }
5218 
5219     /**
5220      *  Convert radians to degrees, with error of less than 0.5 ULP
5221      *  @param x angle in radians
5222      *  @param <T> the type of the field element
5223      *  @return x converted into degrees
5224      */
5225     public static <T extends CalculusFieldElement<T>> T toDegrees(T x) {
5226         return x.toDegrees();
5227     }
5228 
5229     /**
5230      * Multiply a double number by a power of 2.
5231      * @param d number to multiply
5232      * @param n power of 2
5233      * @param <T> the type of the field element
5234      * @return d &times; 2<sup>n</sup>
5235      * @since 1.3
5236      */
5237     public static <T extends CalculusFieldElement<T>> T scalb(final T d, final int n) {
5238         return d.scalb(n);
5239     }
5240 
5241     /**
5242      * Compute least significant bit (Unit in Last Position) for a number.
5243      * @param x number from which ulp is requested
5244      * @param <T> the type of the field element
5245      * @return ulp(x)
5246      * @since 2.0
5247      */
5248     public static <T extends CalculusFieldElement<T>> T ulp(final T x) {
5249         if (Double.isInfinite(x.getReal())) {
5250             return x.newInstance(Double.POSITIVE_INFINITY);
5251         }
5252         return x.ulp();
5253     }
5254 
5255     /** Get the largest whole number smaller than x.
5256      * @param x number from which floor is requested
5257      * @param <T> the type of the field element
5258      * @return a double number f such that f is an integer f &lt;= x &lt; f + 1.0
5259      * @since 1.3
5260      */
5261     public static <T extends CalculusFieldElement<T>> T floor(final T x) {
5262         return x.floor();
5263     }
5264 
5265     /** Get the smallest whole number larger than x.
5266      * @param x number from which ceil is requested
5267      * @param <T> the type of the field element
5268      * @return a double number c such that c is an integer c - 1.0 &lt; x &lt;= c
5269      * @since 1.3
5270      */
5271     public static <T extends CalculusFieldElement<T>> T ceil(final T x) {
5272         return x.ceil();
5273     }
5274 
5275     /** Get the whole number that is the nearest to x, or the even one if x is exactly half way between two integers.
5276      * @param x number from which nearest whole number is requested
5277      * @param <T> the type of the field element
5278      * @return a double number r such that r is an integer r - 0.5 &lt;= x &lt;= r + 0.5
5279      * @since 1.3
5280      */
5281     public static <T extends CalculusFieldElement<T>> T rint(final T x) {
5282         return x.rint();
5283     }
5284 
5285     /** Get the closest long to x.
5286      * @param x number from which closest long is requested
5287      * @param <T> the type of the field element
5288      * @return closest long to x
5289      * @since 1.3
5290      */
5291     public static <T extends CalculusFieldElement<T>> long round(final T x) {
5292         return x.round();
5293     }
5294 
5295     /** Compute the minimum of two values
5296      * @param a first value
5297      * @param b second value
5298      * @param <T> the type of the field element
5299      * @return a if a is lesser or equal to b, b otherwise
5300      * @since 1.3
5301      */
5302     public static <T extends CalculusFieldElement<T>> T min(final T a, final T b) {
5303         final double aR = a.getReal();
5304         final double bR = b.getReal();
5305         if (aR < bR) {
5306             return a;
5307         } else if (bR < aR) {
5308             return b;
5309         } else {
5310             // either the numbers are equal, or one of them is a NaN
5311             return Double.isNaN(aR) ? a : b;
5312         }
5313     }
5314 
5315     /** Compute the minimum of two values
5316      * @param a first value
5317      * @param b second value
5318      * @param <T> the type of the field element
5319      * @return a if a is lesser or equal to b, b otherwise
5320      * @since 1.3
5321      */
5322     public static <T extends CalculusFieldElement<T>> T min(final T a, final double b) {
5323         final double aR = a.getReal();
5324         if (aR < b) {
5325             return a;
5326         } else if (b < aR) {
5327             return a.getField().getZero().add(b);
5328         } else {
5329             // either the numbers are equal, or one of them is a NaN
5330             return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
5331         }
5332     }
5333 
5334     /** Compute the maximum of two values
5335      * @param a first value
5336      * @param b second value
5337      * @param <T> the type of the field element
5338      * @return b if a is lesser or equal to b, a otherwise
5339      * @since 1.3
5340      */
5341     public static <T extends CalculusFieldElement<T>> T max(final T a, final T b) {
5342         final double aR = a.getReal();
5343         final double bR = b.getReal();
5344         if (aR < bR) {
5345             return b;
5346         } else if (bR < aR) {
5347             return a;
5348         } else {
5349             // either the numbers are equal, or one of them is a NaN
5350             return Double.isNaN(aR) ? a : b;
5351         }
5352     }
5353 
5354     /** Compute the maximum of two values
5355      * @param a first value
5356      * @param b second value
5357      * @param <T> the type of the field element
5358      * @return b if a is lesser or equal to b, a otherwise
5359      * @since 1.3
5360      */
5361     public static <T extends CalculusFieldElement<T>> T max(final T a, final double b) {
5362         final double aR = a.getReal();
5363         if (aR < b) {
5364             return a.getField().getZero().add(b);
5365         } else if (b < aR) {
5366             return a;
5367         } else {
5368             // either the numbers are equal, or one of them is a NaN
5369             return Double.isNaN(aR) ? a : a.getField().getZero().add(b);
5370         }
5371     }
5372 
5373     /**
5374      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
5375      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br>
5376      * avoiding intermediate overflow or underflow.
5377      *
5378      * <ul>
5379      * <li> If either argument is infinite, then the result is positive infinity.</li>
5380      * <li> else, if either argument is NaN then the result is NaN.</li>
5381      * </ul>
5382      *
5383      * @param x a value
5384      * @param y a value
5385      * @param <T> the type of the field element
5386      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
5387      * @since 1.3
5388      */
5389     public static <T extends CalculusFieldElement<T>> T hypot(final T x, final T y) {
5390         return x.hypot(y);
5391     }
5392 
5393     /**
5394      * Computes the remainder as prescribed by the IEEE 754 standard.
5395      * <p>
5396      * The remainder value is mathematically equal to {@code x - y*n}
5397      * where {@code n} is the mathematical integer closest to the exact mathematical value
5398      * of the quotient {@code x/y}.
5399      * If two mathematical integers are equally close to {@code x/y} then
5400      * {@code n} is the integer that is even.
5401      * </p>
5402      * <ul>
5403      * <li>If either operand is NaN, the result is NaN.</li>
5404      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
5405      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
5406      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
5407      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
5408      * </ul>
5409      * @param dividend the number to be divided
5410      * @param divisor the number by which to divide
5411      * @param <T> the type of the field element
5412      * @return the remainder, rounded
5413      * @since 1.3
5414      */
5415     public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final double divisor) {
5416         return dividend.remainder(divisor);
5417     }
5418 
5419     /**
5420      * Computes the remainder as prescribed by the IEEE 754 standard.
5421      * <p>
5422      * The remainder value is mathematically equal to {@code x - y*n}
5423      * where {@code n} is the mathematical integer closest to the exact mathematical value
5424      * of the quotient {@code x/y}.
5425      * If two mathematical integers are equally close to {@code x/y} then
5426      * {@code n} is the integer that is even.
5427      * </p>
5428      * <ul>
5429      * <li>If either operand is NaN, the result is NaN.</li>
5430      * <li>If the result is not NaN, the sign of the result equals the sign of the dividend.</li>
5431      * <li>If the dividend is an infinity, or the divisor is a zero, or both, the result is NaN.</li>
5432      * <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
5433      * <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
5434      * </ul>
5435      * @param dividend the number to be divided
5436      * @param divisor the number by which to divide
5437      * @param <T> the type of the field element
5438      * @return the remainder, rounded
5439      * @since 1.3
5440      */
5441     public static <T extends CalculusFieldElement<T>> T IEEEremainder(final T dividend, final T divisor) {
5442         return dividend.remainder(divisor);
5443     }
5444 
5445     /**
5446      * Returns the first argument with the sign of the second argument.
5447      * A NaN {@code sign} argument is treated as positive.
5448      *
5449      * @param magnitude the value to return
5450      * @param sign the sign for the returned value
5451      * @param <T> the type of the field element
5452      * @return the magnitude with the same sign as the {@code sign} argument
5453      * @since 1.3
5454      */
5455     public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, T sign) {
5456         return magnitude.copySign(sign);
5457     }
5458 
5459     /**
5460      * Returns the first argument with the sign of the second argument.
5461      * A NaN {@code sign} argument is treated as positive.
5462      *
5463      * @param magnitude the value to return
5464      * @param sign the sign for the returned value
5465      * @param <T> the type of the field element
5466      * @return the magnitude with the same sign as the {@code sign} argument
5467      * @since 1.3
5468      */
5469     public static <T extends CalculusFieldElement<T>> T copySign(T magnitude, double sign) {
5470         return magnitude.copySign(sign);
5471     }
5472 
5473 //    /**
5474 //     * Print out contents of arrays, and check the length.
5475 //     * <p>used to generate the preset arrays originally.</p>
5476 //     * @param a unused
5477 //     */
5478 //    public static void main(String[] a) {
5479 //        FastMathCalc.printarray(System.out, "EXP_INT_TABLE_A", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_A);
5480 //        FastMathCalc.printarray(System.out, "EXP_INT_TABLE_B", EXP_INT_TABLE_LEN, ExpIntTable.EXP_INT_TABLE_B);
5481 //        FastMathCalc.printarray(System.out, "EXP_FRAC_TABLE_A", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_A);
5482 //        FastMathCalc.printarray(System.out, "EXP_FRAC_TABLE_B", EXP_FRAC_TABLE_LEN, ExpFracTable.EXP_FRAC_TABLE_B);
5483 //        FastMathCalc.printarray(System.out, "LN_MANT",LN_MANT_LEN, lnMant.LN_MANT);
5484 //        FastMathCalc.printarray(System.out, "SINE_TABLE_A", SINE_TABLE_LEN, SINE_TABLE_A);
5485 //        FastMathCalc.printarray(System.out, "SINE_TABLE_B", SINE_TABLE_LEN, SINE_TABLE_B);
5486 //        FastMathCalc.printarray(System.out, "COSINE_TABLE_A", SINE_TABLE_LEN, COSINE_TABLE_A);
5487 //        FastMathCalc.printarray(System.out, "COSINE_TABLE_B", SINE_TABLE_LEN, COSINE_TABLE_B);
5488 //        FastMathCalc.printarray(System.out, "TANGENT_TABLE_A", SINE_TABLE_LEN, TANGENT_TABLE_A);
5489 //        FastMathCalc.printarray(System.out, "TANGENT_TABLE_B", SINE_TABLE_LEN, TANGENT_TABLE_B);
5490 //    }
5491 
5492     /** Enclose large data table in nested static class so it's only loaded on first access. */
5493     private static class ExpIntTable {
5494         /** Exponential evaluated at integer values,
5495          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX].
5496          */
5497         private static final double[] EXP_INT_TABLE_A;
5498         /** Exponential evaluated at integer values,
5499          * exp(x) =  expIntTableA[x + EXP_INT_TABLE_MAX_INDEX] + expIntTableB[x+EXP_INT_TABLE_MAX_INDEX]
5500          */
5501         private static final double[] EXP_INT_TABLE_B;
5502 
5503         static {
5504             if (RECOMPUTE_TABLES_AT_RUNTIME) {
5505                 EXP_INT_TABLE_A = new double[FastMath.EXP_INT_TABLE_LEN];
5506                 EXP_INT_TABLE_B = new double[FastMath.EXP_INT_TABLE_LEN];
5507 
5508                 final double tmp[] = new double[2];
5509                 final double recip[] = new double[2];
5510 
5511                 // Populate expIntTable
5512                 for (int i = 0; i < FastMath.EXP_INT_TABLE_MAX_INDEX; i++) {
5513                     FastMathCalc.expint(i, tmp);
5514                     EXP_INT_TABLE_A[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[0];
5515                     EXP_INT_TABLE_B[i + FastMath.EXP_INT_TABLE_MAX_INDEX] = tmp[1];
5516 
5517                     if (i != 0) {
5518                         // Negative integer powers
5519                         FastMathCalc.splitReciprocal(tmp, recip);
5520                         EXP_INT_TABLE_A[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[0];
5521                         EXP_INT_TABLE_B[FastMath.EXP_INT_TABLE_MAX_INDEX - i] = recip[1];
5522                     }
5523                 }
5524             } else {
5525                 EXP_INT_TABLE_A = FastMathLiteralArrays.loadExpIntA();
5526                 EXP_INT_TABLE_B = FastMathLiteralArrays.loadExpIntB();
5527             }
5528         }
5529     }
5530 
5531     /** Enclose large data table in nested static class so it's only loaded on first access. */
5532     private static class ExpFracTable {
5533         /** Exponential over the range of 0 - 1 in increments of 2^-10
5534          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
5535          * 1024 = 2^10
5536          */
5537         private static final double[] EXP_FRAC_TABLE_A;
5538         /** Exponential over the range of 0 - 1 in increments of 2^-10
5539          * exp(x/1024) =  expFracTableA[x] + expFracTableB[x].
5540          */
5541         private static final double[] EXP_FRAC_TABLE_B;
5542 
5543         static {
5544             if (RECOMPUTE_TABLES_AT_RUNTIME) {
5545                 EXP_FRAC_TABLE_A = new double[FastMath.EXP_FRAC_TABLE_LEN];
5546                 EXP_FRAC_TABLE_B = new double[FastMath.EXP_FRAC_TABLE_LEN];
5547 
5548                 final double tmp[] = new double[2];
5549 
5550                 // Populate expFracTable
5551                 final double factor = 1d / (EXP_FRAC_TABLE_LEN - 1);
5552                 for (int i = 0; i < EXP_FRAC_TABLE_A.length; i++) {
5553                     FastMathCalc.slowexp(i * factor, tmp);
5554                     EXP_FRAC_TABLE_A[i] = tmp[0];
5555                     EXP_FRAC_TABLE_B[i] = tmp[1];
5556                 }
5557             } else {
5558                 EXP_FRAC_TABLE_A = FastMathLiteralArrays.loadExpFracA();
5559                 EXP_FRAC_TABLE_B = FastMathLiteralArrays.loadExpFracB();
5560             }
5561         }
5562     }
5563 
5564     /** Enclose large data table in nested static class so it's only loaded on first access. */
5565     private static class lnMant {
5566         /** Extended precision logarithm table over the range 1 - 2 in increments of 2^-10. */
5567         private static final double[][] LN_MANT;
5568 
5569         static {
5570             if (RECOMPUTE_TABLES_AT_RUNTIME) {
5571                 LN_MANT = new double[FastMath.LN_MANT_LEN][];
5572 
5573                 // Populate lnMant table
5574                 for (int i = 0; i < LN_MANT.length; i++) {
5575                     final double d = Double.longBitsToDouble( (((long) i) << 42) | 0x3ff0000000000000L );
5576                     LN_MANT[i] = FastMathCalc.slowLog(d);
5577                 }
5578             } else {
5579                 LN_MANT = FastMathLiteralArrays.loadLnMant();
5580             }
5581         }
5582     }
5583 
5584     /** Enclose the Cody/Waite reduction (used in "sin", "cos" and "tan"). */
5585     private static class CodyWaite {
5586         /** k */
5587         private final int finalK;
5588         /** remA */
5589         private final double finalRemA;
5590         /** remB */
5591         private final double finalRemB;
5592 
5593         /**
5594          * @param xa Argument.
5595          */
5596         CodyWaite(double xa) {
5597             // Estimate k.
5598             //k = (int)(xa / 1.5707963267948966);
5599             int k = (int)(xa * 0.6366197723675814);
5600 
5601             // Compute remainder.
5602             double remA;
5603             double remB;
5604             while (true) {
5605                 double a = -k * 1.570796251296997;
5606                 remA = xa + a;
5607                 remB = -(remA - xa - a);
5608 
5609                 a = -k * 7.549789948768648E-8;
5610                 double b = remA;
5611                 remA = a + b;
5612                 remB += -(remA - b - a);
5613 
5614                 a = -k * 6.123233995736766E-17;
5615                 b = remA;
5616                 remA = a + b;
5617                 remB += -(remA - b - a);
5618 
5619                 if (remA > 0) {
5620                     break;
5621                 }
5622 
5623                 // Remainder is negative, so decrement k and try again.
5624                 // This should only happen if the input is very close
5625                 // to an even multiple of pi/2.
5626                 --k;
5627             }
5628 
5629             this.finalK = k;
5630             this.finalRemA = remA;
5631             this.finalRemB = remB;
5632         }
5633 
5634         /**
5635          * @return k
5636          */
5637         int getK() {
5638             return finalK;
5639         }
5640         /**
5641          * @return remA
5642          */
5643         double getRemA() {
5644             return finalRemA;
5645         }
5646         /**
5647          * @return remB
5648          */
5649         double getRemB() {
5650             return finalRemB;
5651         }
5652     }
5653 
5654 }